By following the "Tutorial on using NR3 code from within MATLAB", I was able to work with vector and 2D matrix between matlab and C++. But how about 3D matrix?
Suppose I have a 3D matrix in Matlab: A(4x5x6) , how to access it by using 'Mat3DDoub'?

#include "nr3matlab.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
// Int zi,zj,zk, zm,zn,zo;
// matrix rhs arguments and lhs return values
// Mat3DDoub b(prhs[0]); // bind a rhs argument
Mat3DDoub b(?); // how to initialize b, so b=prhs

}

Thanks for your help!

...
Suppose I have a 3D matrix in Matlab:... using 'Mat3DDoub'?

I can give an example that works with GNU Octave. As far as I know, it "should" work with Matlab, but I have no way of testing.

Notes:
1. Arrays in Matlab and Octave are Fortran-style arrays. That is, elements are stored in contiguous memory locations in column-major order.

2. Mat3DDoub objects have data that is stored C++-style (row-major order) with pointers to pointers set up so the data can be accessed with C++-style indexed notation. (So, if you try something like memcpy() to copy the arrays, you will have to "transpose" from one storage scheme to the other, even though the Mat3DDoub data values are stored in contiguous memory.)

3. The Mat3DDoub class, as supplied in the Numerical Recipes code, does not have a copy constructor or an overloaded assignment operator, so you can't just set one equal to another, and functions with Mat3DDoub arguments should always call by reference. I mean, the compiler won't complain if you ignore my warnings (it's not illegal in C++ to do such things), but I'm betting that the default copy constructor and assignment operator really, really (really) won't do what you have in mind.

Here's the cpp function source of a simple test program:

/*
* test3d.cpp illustrates use of GNU Octave 3D arrays used
* in calling c++ functions
*
* As far as I know, it should work for Matlab, but...
* (Your Mileage May Vary)
*
* davekw7x
*/

/* This includes just about all of the usual C++ headers, among other things */
#include "nr3matlab.h"


/* Function prototypes */

/* Copy colum-major array contents to C++ 3D array */
void Fortran_to_C(Mat3DDoub & dest, Doub * src);

/* Copy C++ 3D array to column-major array contents */
void C_to_Fortran(Doub * dest, Mat3DDoub & src);

/* Here's the action function */
void Foo(Mat3DDoub_I & x, Mat3DDoub & y);


/* This function always has this signature */
void mexFunction(Int nlhs, mxArray * plhs[], Int nrhs, const mxArray * prhs[])
{
Int number_of_dimensions;
const Int *dimension_array;

Doub *input_pointer;
Doub *output_pointer;

/* Check for proper number of input and output arguments. */
if (nrhs != 1) {
mexErrMsgTxt("One input argument is required.");
}
if (nlhs > 1) {
mexErrMsgTxt("Only one output argument is allowed.");
}

if (!(mxIsDouble(prhs[0]))) {
mexErrMsgTxt("Input array must be of type double.");
}

number_of_dimensions = mxGetNumberOfDimensions(prhs[0]);
if (number_of_dimensions != 3) {
mexErrMsgTxt("Input must be 3-D\n");
}
/*
* The dimension_array pointer gives the size of each dimension
*/
dimension_array = mxGetDimensions(prhs[0]);

/*
* Declare and allocate storage for C++ 3D input and output arrays.
* The C++ arrays are the same shape/size as the Fortran arrays.
*/
Mat3DDoub nrin(dimension_array[0], dimension_array[1], dimension_array[2]);
Mat3DDoub nrout(dimension_array[0], dimension_array[1], dimension_array[2]);

/*
* Allocate the space for the return argument.
* This will be same shape/size as the input argument.
*/
plhs[0] = mxCreateNumericArray(
number_of_dimensions,
dimension_array,
mxDOUBLE_CLASS,
mxREAL
);
/*
* Get addresses of start of data regions of octave arguments
*/
input_pointer = (Doub *)mxGetData(prhs[0]);
output_pointer = (Doub *)mxGetData(plhs[0]);

/* Read input data into nrin */
Fortran_to_C(nrin, input_pointer);

/* Work on nrin with some C++ function */
Foo(nrin, nrout);

/*
* Now Convert the resulting array back to a column-major
* equivalent and copy to output array region
*/
C_to_Fortran(output_pointer, nrout);

#define _DEBUG__
#ifdef _DEBUG__
/* For debugging purposes: Print out the input array */
mexPrintf("In test3d.mex:\n");
for (Int i = 0; i < nrin.dim1(); i++) {
for (Int j = 0; j < nrin.dim2(); j++) {
for (Int k = 0; k < nrin.dim3(); k++) {
mexPrintf("nrin[%d][%d][%d] = %2.0f ", i, j, k, nrin[i][j][k]);
}
mexPrintf("\n");
}
mexPrintf("\n");
}
for (Int i = 0; i < nrout.dim1(); i++) {
for (Int j = 0; j < nrout.dim2(); j++) {
for (Int k = 0; k < nrout.dim3(); k++) {
mexPrintf("nrout[%d][%d][%d] = %2.0f ", i, j, k, nrout[i][j][k]);
}
mexPrintf("\n");
}
mexPrintf("\n");
}
mexPrintf("Leaving test3d.mex\n\n");
#endif
}

/*
* Convert a Fortran array(W,X,Y) to a C++ array[W][X][Y]
*/
void Fortran_to_C(Mat3DDoub & dest, Doub * src)
{
for (Int k = 0; k < dest.dim3(); k++) {
for (Int j = 0; j < dest.dim2(); j++) {
for (Int i = 0; i < dest.dim1(); i++) {
dest[i][j][k] = *(src++);
}
}
}
}

/*
* Convert a C++ array[W][X][Y] to a Fortran array(W,X,Y)
*/
void C_to_Fortran(Doub * dest, Mat3DDoub & src)
{
for (Int k = 0; k < src.dim3(); k++) {
for (Int j = 0; j < src.dim2(); j++) {
for (Int i = 0; i < src.dim1(); i++) {
*(dest++) = src[i][j][k];
}
}
}
}

/*------------------------------*/

/* Here's where you would call a Numerical Recipes or other
* C++ function that has a Mat3DDoub argument
*
* I'll just do something simple, like changing each element
* by a different value.
*/
void Foo(Mat3DDoub_I & x, Mat3DDoub & y)
{
Doub adder = 1;
for (Int i = 0; i < x.dim1(); i++) {
for (Int j = 0; j < x.dim2(); j++){
for (Int k = 0; k < x.dim3(); k++) {
y[i][j][k] = x[i][j][k] + adder;
adder += 2;
}
}
}
}

/*------------------------------*/

Here's a little Octave test script:

#
# xtest3d.m
#
# Create a 3-D array, call a mex function to process it
#
Dim1 = 2;
Dim2 = 3;
Dim3 = 4;
value = 1.0;
for i = 1:Dim1
for j = 1:Dim2
for k = 1:Dim3
x(i,j,k) = value;
value = value+1.0;
end
end
end
printf('Originally\n');
print3d(x);

y=test3d(x);

printf('After test3d\n');
print3d(y);

function print3d(m)
for i = 1:size(m,1)
for j = 1:size(m,2)
for k = 1:size(m,3)
printf('m(%d,%d,%d) = %3.0f ', i, j, k, m(i,j,k));
end
printf('\n');
end
printf('\n');
end


And here's the output:

octave:1> xtest3d
Originally
m(1,1,1) = 1 m(1,1,2) = 2 m(1,1,3) = 3 m(1,1,4) = 4
m(1,2,1) = 5 m(1,2,2) = 6 m(1,2,3) = 7 m(1,2,4) = 8
m(1,3,1) = 9 m(1,3,2) = 10 m(1,3,3) = 11 m(1,3,4) = 12

m(2,1,1) = 13 m(2,1,2) = 14 m(2,1,3) = 15 m(2,1,4) = 16
m(2,2,1) = 17 m(2,2,2) = 18 m(2,2,3) = 19 m(2,2,4) = 20
m(2,3,1) = 21 m(2,3,2) = 22 m(2,3,3) = 23 m(2,3,4) = 24

In test3d.mex:
nrin[0][0][0] = 1 nrin[0][0][1] = 2 nrin[0][0][2] = 3 nrin[0][0][3] = 4
nrin[0][1][0] = 5 nrin[0][1][1] = 6 nrin[0][1][2] = 7 nrin[0][1][3] = 8
nrin[0][2][0] = 9 nrin[0][2][1] = 10 nrin[0][2][2] = 11 nrin[0][2][3] = 12

nrin[1][0][0] = 13 nrin[1][0][1] = 14 nrin[1][0][2] = 15 nrin[1][0][3] = 16
nrin[1][1][0] = 17 nrin[1][1][1] = 18 nrin[1][1][2] = 19 nrin[1][1][3] = 20
nrin[1][2][0] = 21 nrin[1][2][1] = 22 nrin[1][2][2] = 23 nrin[1][2][3] = 24

nrout[0][0][0] = 2 nrout[0][0][1] = 5 nrout[0][0][2] = 8 nrout[0][0][3] = 11
nrout[0][1][0] = 14 nrout[0][1][1] = 17 nrout[0][1][2] = 20 nrout[0][1][3] = 23
nrout[0][2][0] = 26 nrout[0][2][1] = 29 nrout[0][2][2] = 32 nrout[0][2][3] = 35

nrout[1][0][0] = 38 nrout[1][0][1] = 41 nrout[1][0][2] = 44 nrout[1][0][3] = 47
nrout[1][1][0] = 50 nrout[1][1][1] = 53 nrout[1][1][2] = 56 nrout[1][1][3] = 59
nrout[1][2][0] = 62 nrout[1][2][1] = 65 nrout[1][2][2] = 68 nrout[1][2][3] = 71

Leaving test3d.mex

After test3d
m(1,1,1) = 2 m(1,1,2) = 5 m(1,1,3) = 8 m(1,1,4) = 11
m(1,2,1) = 14 m(1,2,2) = 17 m(1,2,3) = 20 m(1,2,4) = 23
m(1,3,1) = 26 m(1,3,2) = 29 m(1,3,3) = 32 m(1,3,4) = 35

m(2,1,1) = 38 m(2,1,2) = 41 m(2,1,3) = 44 m(2,1,4) = 47
m(2,2,1) = 50 m(2,2,2) = 53 m(2,2,3) = 56 m(2,2,4) = 59
m(2,3,1) = 62 m(2,3,2) = 65 m(2,3,3) = 68 m(2,3,4) = 71


Regards,

Dave

Footnote: I have shown some care in testing input argument characteristics. I have found that if mexFunction() is fed stuff that it doesn't like, Octave is likely to crash. I don't know whether Matlab actually crashes or not, but letting mexErrMsgTxt() print some useful, intelligible information and gracefully abort the function and return to the Octave prompt is much less unsettling. I'm funny that way.

[QUOTE=davekw7x;3503]I can give an example that works with GNU Octave. As far as I know, it "should" work with Matlab, but I have no way of testing.

Dave,

I just tested with matlab, and it works very well.

Another stupid question: due to memory and speed limit, i have to use 3d matrix of single precision, instead of double precision. Do we have "Mat3DSingle" or the similar?

Thank you very much for your patient explanations and helpful codes!

Drew

[...3d matrix of single precision...To the best of my knowledge and experience GNU Octave doesn't have a single precision data type. I have no way of testing Matlab, so maybe you can look at Matlab documention. I know there is a "single" operator in Matlab that makes numerical quantities to be single precision, but I don't know about passing single precision arguments to mexFunction(). I'm guessing it would work, but...

Do we have "Mat3DSingle" or the similar? If by "we" you mean people who use unmodified Numerical Recipes distributed code, the answer is no.

However...

The Mat3DDoub class is defined with a typedef using the NRMat3D templated class, as you can see in nr3.h and nr3matlab.h.

So: If you really (really) need single precision 3D arrays (in C++) that are similar to Numerical Recipes Mat3DDouble objects, you could try your own typedefs to create a Mat3DSingle class.

I have not tested it, but maybe you could try:


#include "nr3matlab.h"
typedef float Float;
typedef const NRMat3d<Float> Mat3DFloat_I;
typedef NRMat3d<Float> Mat3DFloat, Mat3DFloat_O, Mat3DFloat_IO;


Or some such thing.


Regards,

Dave