diff --git a/src/MatlabReader.cpp b/src/MatlabReader.cpp
index 8218ff5..2654f2a 100644
--- a/src/MatlabReader.cpp
+++ b/src/MatlabReader.cpp
@@ -67,7 +67,32 @@ Aurora::Matrix MatlabReader::read(const std::string& aFieldName)
     Aurora::Matrix result;
     if(!isComplex)
     {
-        result = Aurora::Matrix::copyFromRawData(reinterpret_cast<double*>(matvar->data),rows,columns,slices);
+        size_t matrixSize =  matvar->nbytes / static_cast<unsigned long long>(matvar->data_size);
+        double* matrixData = new double[matrixSize];
+        switch(matvar->data_type )
+        {
+            case MAT_T_INT16:
+            {
+                std::copy((short*)matvar->data, (short*)matvar->data + matrixSize, matrixData);
+                break;
+            }
+            case MAT_T_UINT16:
+            {
+                std::copy((unsigned short*)matvar->data, (unsigned short*)matvar->data + matrixSize, matrixData);
+                break;
+            }
+            case MAT_T_INT32:
+            {
+                std::copy((int*)matvar->data, (int*)matvar->data + matrixSize, matrixData);
+                break;
+            }
+            case MAT_T_DOUBLE:
+            {
+                std::copy((double*)matvar->data, (double*)matvar->data + matrixSize, matrixData);
+                break;
+            }
+        }
+        result = Aurora::Matrix::fromRawData(matrixData,rows,columns,slices);
     }
     else
     {
@@ -87,3 +112,42 @@ Aurora::Matrix MatlabReader::read(const std::string& aFieldName)
     Mat_VarFree(matvar);
     return result;
 }
+
+std::vector<Aurora::Matrix> MatlabReader::read4d(const std::string& aFieldName)
+{
+    if(mMatFile == nullptr)
+    {
+        return std::vector<Aurora::Matrix>();
+    }
+
+    matvar_t* matvar = Mat_VarRead(mMatFile, aFieldName.c_str());
+    if (matvar == nullptr)
+    {
+        return std::vector<Aurora::Matrix>();
+    }
+    int rank = matvar->rank;
+    if (rank != 4)
+    {
+        return std::vector<Aurora::Matrix>();
+    }
+    bool isComplex = matvar->isComplex;
+    int rows = matvar->dims[0];
+    int columns = matvar->dims[1];
+    int slices = matvar->dims[2];
+    int dims4 = matvar->dims[3];
+    std::vector<Aurora::Matrix> result;
+    if(!isComplex)
+    {
+        double* matrix3d = reinterpret_cast<double*>(matvar->data);
+        for(int i=0; i<dims4; ++i)
+        {
+            result.push_back(Aurora::Matrix::copyFromRawData(matrix3d,rows,columns,slices));
+            matrix3d+=rows*columns*slices;
+        }
+    }
+    else
+    {
+        return std::vector<Aurora::Matrix>();
+    }
+    return result;
+}
diff --git a/src/MatlabReader.h b/src/MatlabReader.h
index 0504ac8..0b77892 100644
--- a/src/MatlabReader.h
+++ b/src/MatlabReader.h
@@ -12,6 +12,7 @@ public:
     
     explicit MatlabReader(const std::string& aFilePath);
     ~MatlabReader();
+    std::vector<Aurora::Matrix> read4d(const std::string& aFieldName);
     Aurora::Matrix read(const std::string& aFieldName);
     std::shared_ptr<short> readint16(const std::string& aFieldName,size_t& size);