MATLAB: How to convert a 4D matrix from Matlab to C++

Question

I am trying to convert a 4D double matrix from Matlab to C++, using the below code however, it looks to me that mxGetPr() does not support linear mapping! are there any mapping indexing formula that can be used with this function (i.e ptr1[linear indexing formula] in the code below)?

MATFile *pmat;
int dt1,dt2,dt3,dt4;
  dt1=3; dt2=4; dt3=5;  dt4=6; % diemntions of the 4D matrix
    vector<vector<vector<vector<double>>>> array4D1(dt1, vector<vector<vector<double>>>
                                           (dt2, vector<vector<double>>
                                           (dt3, vector<double>
 pa2 = matGetVariable(pmat, "G_wp");    % G_wp is a 3x4x5x6 matlab matrix                                         (dt4))));
ptr1 = mxGetPr( pa2 ); 
int idx1=0;
for (int i=0; i<dt1; i++ )
  {for (int j=0; j<dt2; j++)
    {for (int k=0; k<dt3; k++)
      {for (int l=0; l<dt4; l++)
                           {array4D1[i][j][k][l]=ptr1[idx1];idx1++;}
                      }
              }
      }

Answer 1

C array storage order is reversed from MATLAB storage order. This holds regardless of the number of dimensions. And you certainly can use "linear indexing" on the MATLAB data memory. E.g., a 2D example:

#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
/* Assume 2x3 double input, no argument checking */
    double x[3][2]; /* C array dimensions are reverse order from MATLAB */
    double *pr;
    mwSize i, j, m, n;
      pr = mxGetPr(prhs[0]);
      m = mxGetM(prhs[0]); /* 2 */
      n = mxGetN(prhs[0]); /* 3 */
  /* for loops & limits in reverse index order */
      for( j=0; j<n; j++ ) {
          for( i=0; i<m; i++ ) {
              x[j][i] = *pr++;
          }
      }
      for( i=0; i<m; i++ ) {
          for( j=0; j<n; j++ ) {
              mexPrintf("x[%d][%d] = %f\n",i,j,x[j][i]);
          }
      }
      mexCallMATLAB(0,NULL,1,prhs,"display");
  }

And at the command line:

>> mex matrix2d.c
>> x = [1 2 3;4 5 6]
x =
     1     2     3
     4     5     6
>> matrix2d(x)
x[0][0] = 1.000000
x[0][1] = 2.000000
x[0][2] = 3.000000
x[1][0] = 4.000000
x[1][1] = 5.000000
x[1][2] = 6.000000
x =
     1     2     3
     4     5     6

So, everything is as expected as long as you get the reverse order indexing coded correctly.

In your code above, you don't post enough information for me to be able to tell what the dimensions of the mxArray pa2 are, you simply have hard coded the following:

   dt1=3; dt2=4; dt3=5;  dt4=6;

So I can't say for sure how to correct your code. But if you just extend my 2D example to 4D it should work for you. Just do the dimensions on the C++ side in reverse order from MATLAB.

EDIT 7/26/2016:

E.g., here is some C++ code to read a variable from a mat file and use your 4D vector approach:

#include <vector>
using namespace std;
#include "mex.h"
#include "mat.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    MATFile *pmat;
    const mwSize *dims;
    mwSize ndim;
    mwSize dt1, dt2, dt3, dt4;
    mxArray *pa2;
    double *ptr1;
    int a,b,c,d;
      pmat = matOpen("matGwp.mat","r");
      if( !pmat ) {
          mexErrMsgTxt("Cannot open matGwp.mat");
      }
      pa2 = matGetVariable(pmat, "G_wp");
      matClose(pmat);
      if( pa2 == NULL ) {
          mexErrMsgTxt("Cannot get G_wp");
      }
      dims = mxGetDimensions(pa2);
      ndim = mxGetNumberOfDimensions(pa2);
      if( ndim != 4 ) {
          mexErrMsgTxt("G_wp needs to have 4 explicit dimensions");
      }
  // Set up array dimensions in reverse order
      dt4 = dims[0];
      dt3 = dims[1];
      dt2 = dims[2];
      dt1 = dims[3];
      vector<vector<vector<vector<double>>>> array4D1(dt1, vector<vector<vector<double>>>
                                                     (dt2, vector<vector<double>>
                                                     (dt3, vector<double>
                                                     (dt4)               )       )       );
  // Copy elements in reverse order
      ptr1 = mxGetPr( pa2 ); 
      int idx1=0;
      for (int i=0; i<dt1; i++ ) {
          for (int j=0; j<dt2; j++) {
              for (int k=0; k<dt3; k++) {
                  for (int l=0; l<dt4; l++) {
                      array4D1[i][j][k][l] = ptr1[idx1];
                      idx1++;
                  }
              }
          }
      }
  // Some test output
      a = 1; b = 1; c = 1; d = 1;
      mexPrintf("G_wp(%d,%d,%d,%d) = array4D1[%d][%d][%d][%d] = %f\n",a,b,c,d,d-1,c-1,b-1,a-1,array4D1[d-1][c-1][b-1][a-1]);
      a = 2; b = 4; c = 1; d = 6;
      mexPrintf("G_wp(%d,%d,%d,%d) = array4D1[%d][%d][%d][%d] = %f\n",a,b,c,d,d-1,c-1,b-1,a-1,array4D1[d-1][c-1][b-1][a-1]);
      a = 3; b = 2; c = 3; d = 2;
      mexPrintf("G_wp(%d,%d,%d,%d) = array4D1[%d][%d][%d][%d] = %f\n",a,b,c,d,d-1,c-1,b-1,a-1,array4D1[d-1][c-1][b-1][a-1]);
      a = 3; b = 4; c = 5; d = 6;
      mexPrintf("G_wp(%d,%d,%d,%d) = array4D1[%d][%d][%d][%d] = %f\n",a,b,c,d,d-1,c-1,b-1,a-1,array4D1[d-1][c-1][b-1][a-1]);
  }

And here is the example using that code:

>> G_wp = reshape(1:3*4*5*6,3,4,5,6);
>> save matGwp
>> mex mat2vector.cpp
>> mat2vector
G_wp(1,1,1,1) = array4D1[0][0][0][0] = 1.000000
G_wp(2,4,1,6) = array4D1[5][0][3][1] = 311.000000
G_wp(3,2,3,2) = array4D1[1][2][1][2] = 90.000000
G_wp(3,4,5,6) = array4D1[5][4][3][2] = 360.000000
>> G_wp(1,1,1,1)
ans =
     1
>> G_wp(2,4,1,6)
ans =
   311
>> G_wp(3,2,3,2)
ans =
    90
>> G_wp(3,4,5,6)
ans =
   360

So, again, as long as you get the dimensions in reverse order you can make it look like the C example.

---------------------------------

If you want the dimensions to be in the same order in C++ as they are in MATLAB, and are willing to accept that they will be in a different order in memory, you can do this:

#include <vector>
using namespace std;
#include "mex.h"
#include "mat.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    MATFile *pmat;
    const mwSize *dims;
    mwSize ndim;
    mwSize dt1, dt2, dt3, dt4;
    mxArray *pa2;
    double *ptr1;
    int a,b,c,d;
      pmat = matOpen("matGwp.mat","r");
      if( !pmat ) {
          mexErrMsgTxt("Cannot open matGwp.mat");
      }
      pa2 = matGetVariable(pmat, "G_wp");
      matClose(pmat);
      if( pa2 == NULL ) {
          mexErrMsgTxt("Cannot get G_wp");
      }
      dims = mxGetDimensions(pa2);
      ndim = mxGetNumberOfDimensions(pa2);
      if( ndim != 4 ) {
          mexErrMsgTxt("G_wp needs to have 4 explicit dimensions");
      }
  // Set up array dimensions in same order
      dt4 = dims[3];
      dt3 = dims[2];
      dt2 = dims[1];
      dt1 = dims[0];
      vector<vector<vector<vector<double>>>> array4D1(dt1, vector<vector<vector<double>>>
                                                     (dt2, vector<vector<double>>
                                                     (dt3, vector<double>
                                                     (dt4)               )       )       );
  // Copy elements in same order (but will be reversed in memory compared to MATLAB)
      ptr1 = mxGetPr( pa2 ); 
      int idx1=0;
      for (int l=0; l<dt4; l++) {
          for (int k=0; k<dt3; k++) {
              for (int j=0; j<dt2; j++) {
                  for (int i=0; i<dt1; i++ ) {
                      array4D1[i][j][k][l] = ptr1[idx1];
                      idx1++;
                  }
              }
          }
      }
  // Some test output (dimensions in same order)
      a = 1; b = 1; c = 1; d = 1;
      mexPrintf("G_wp(%d,%d,%d,%d) = array4D1[%d][%d][%d][%d] = %f\n",a,b,c,d,a-1,b-1,c-1,d-1,array4D1[a-1][b-1][c-1][d-1]);
      a = 2; b = 4; c = 1; d = 6;
      mexPrintf("G_wp(%d,%d,%d,%d) = array4D1[%d][%d][%d][%d] = %f\n",a,b,c,d,a-1,b-1,c-1,d-1,array4D1[a-1][b-1][c-1][d-1]);
      a = 3; b = 2; c = 3; d = 2;
      mexPrintf("G_wp(%d,%d,%d,%d) = array4D1[%d][%d][%d][%d] = %f\n",a,b,c,d,a-1,b-1,c-1,d-1,array4D1[a-1][b-1][c-1][d-1]);
      a = 3; b = 4; c = 5; d = 6;
      mexPrintf("G_wp(%d,%d,%d,%d) = array4D1[%d][%d][%d][%d] = %f\n",a,b,c,d,a-1,b-1,c-1,d-1,array4D1[a-1][b-1][c-1][d-1]);
  }

And the command line results are again as expected:

>> mex mat2vector2.cpp
>> mat2vector2
G_wp(1,1,1,1) = array4D1[0][0][0][0] = 1.000000
G_wp(2,4,1,6) = array4D1[1][3][0][5] = 311.000000
G_wp(3,2,3,2) = array4D1[2][1][2][1] = 90.000000
G_wp(3,4,5,6) = array4D1[2][3][4][5] = 360.000000

------------------------------------------------

CAUTION: All of this assumes that you want to go through all of this vector stuff simply in order to preserve the 4D coding [][][][] syntax. You do gain the syntax, but you lose other capabilities because of this. E.g., you can no longer simply pass the pointer to a 2D slice of a double matrix to another routine for processing (e.g., BLAS or LAPACK).