Hi,

Would someone help in generating, accessing 4-d arrays with NR types? In the book it recommends the following for 3-D arrays.

vector<MatDoub> threedee(17);
for (Int i=0;i<17;i++) threedee[i].resize(19,21);

How to extend it to 4-D arrays?

Thanks

4-d arrays with NR types?What are you going to do with the 4-D arrays?

The reason that I ask is this: If I were going to have arrays of arrays of MatDoubs, then I might do it one way. If I were going to have matrices of MatDoubs, I might do it another way. If I were going to have 4-D arrays of Doubs, I might do it another way.

I mean, the usage might appear similar (you could address the individual Doubs with 4-D index notation), but there would be a big difference in how they would be used, somehow, with NR functions. (Or, are you going to use the 4-D arrays for something else?)


Bottom line:
Give us an example of how you are going to use them, and then we might be able to recommend some typedefs (or other program constructs) to make it "easier" to declare and initialize them with arbitrary sizes.


Regards,

Dave

Hi Dave,

In my application I am taking a 2d transform of an image and what I need is just to store matrices of images. It doesn't seem to matter at all which one I use since access and elementwise addition operations are the ones that will be frequently used. So, I think I should go with the one which is the most efficient speedwise (runtime). Hope it turns out to be the syntactically the easiest one.

Thanks,
Onur

...what I need is just to store matrices of images. ...

So, if each image is a MatDoub, for example, you can try something like the following:

#include "../code/nr3.h"

// A vector of MatDoubs
typedef std::vector<MatDoub> VecMatDoub;

// A vector of vector of MatDoubs is like a "matrix of MatDoubs"
typedef std::vector<VecMatDoub> MatMatDoub;

int main()
{
//
// x3 is a vector of MatDoubs so x3[i] is a MatDoub
//
VecMatDoub x3(5); // Vector with 5 MatDoubs
for (size_t i = 0; i < x3.size(); i++) {
x3[i].resize(19,21);
}
cout << "x3.size() = " << x3.size() << endl
<< "x3[0] (rows,cols) = (" << x3[0].nrows()
<< "," << x3[0].ncols() << ")"
<< endl << endl;

//
// x4 is vector of vector of MatDoubs so x4[i] is a vector
// of MatDoubs and x4[i][j] is a MatDoub
//
MatMatDoub x4(5);
for (size_t i = 0; i < x4.size(); i++) {
x4[i].resize(2); // Two MatDoubs
cout << "x4[" << i << "].size() = " << x4[i].size() << endl;
for (size_t j = 0; j < x4[i].size(); j++) {
x4[i][j].resize(3, 4);
cout << "x4[" << i << "][" << j << "] (rows,cols) = ("
<< x4[i][j].nrows() << ","
<< x4[i][j].ncols()
<< ")" << endl;
}
cout << endl;
}

return 0;
}


Output:

x3.size() = 5
x3[0] (rows,cols) = (19,21)

x4[0].size() = 2
x4[0][0] (rows,cols) = (3,4)
x4[0][1] (rows,cols) = (3,4)

x4[1].size() = 2
x4[1][0] (rows,cols) = (3,4)
x4[1][1] (rows,cols) = (3,4)

x4[2].size() = 2
x4[2][0] (rows,cols) = (3,4)
x4[2][1] (rows,cols) = (3,4)

x4[3].size() = 2
x4[3][0] (rows,cols) = (3,4)
x4[3][1] (rows,cols) = (3,4)

x4[4].size() = 2
x4[4][0] (rows,cols) = (3,4)
x4[4][1] (rows,cols) = (3,4)

[dave@starrx vec4d]$ m z
make TARGET=z
make: `z' is up to date.
[dave@starrx vec4d]$ z
x3.size() = 5
x3[0] (rows,cols) = (19,21)

x4[0].size() = 2
x4[0][0] (rows,cols) = (3,4)
x4[0][1] (rows,cols) = (3,4)

x4[1].size() = 2
x4[1][0] (rows,cols) = (3,4)
x4[1][1] (rows,cols) = (3,4)

x4[2].size() = 2
x4[2][0] (rows,cols) = (3,4)
x4[2][1] (rows,cols) = (3,4)

x4[3].size() = 2
x4[3][0] (rows,cols) = (3,4)
x4[3][1] (rows,cols) = (3,4)

x4[4].size() = 2
x4[4][0] (rows,cols) = (3,4)
x4[4][1] (rows,cols) = (3,4)

Note that in C++ there are more "sophisticated" ways of declaring initial sizes pf objects of vectors and vectors of vectors, but I chose the "naive" method that may be more comforting for some people.

Also note that, since the 3-D and 4-D things are defined in terms of std::vector objects, you could declare them with no size given and use push_back() to add elements instead of resize to create an object of a given size.

Anyhow...

For example, suppose x4[a][b] is one of your MatDoub matrices. To access element on row[c] and column[d] of that Matdoub, you use the notation x4[a][b][c][d]. To copy MatDoub x4[a][b] to x4[x][y], you can simply use assignment:
x4[x][y] = x4[a][b]

See Footnote.

Regards,

Dave
Footnote:
There is kind of an interesting story of why there is no standard matrix class in the current version of Standard C++, but I won't bore you with the details.

Note that the MatDoub class doesn't have any operators except assignment and subscripting.

If you want to do other operations (add, subtract, multiply by constants, etc.) you can either modify the NRmatrix class definition in your copy of nr3.h to add the operators or you can create your own derived class that defines the extra operators that you need. (Or you can simply put the code to do the operations in your application software.)

Note that you can find matrix classes that are compatible with NR3 functions and have some additional stuff here (among other places): http://www.nr.com/contrib

Then, you can simply use the typedefs in my example, using std::vectors and std::vectors of std::vectors of that class instead of MatDoub.

Chacun Ã* son goût!

Thanks for your great support Dave. Things should now work out for me...

Dave or anyone interested in answering,

Is it possible to return such multidimensional arrays from a mex function? For MatDoub type the recommended style is as follows:

MatDoub output(m,n,plhs[0]);

For VecMatDoub type, is it possible to somehow generalize the above approach?

Thanks,
Onur

Forgot to mention that size of VecMatDoub I'm trying to return to MATLAB workspace is variable (bound to a function parameter).

...to return such multidimensional arrays from a mex function?..

Shortest answer: No.

Short answer: Not with the typedefs of my example, because this application requires that the elements of the arrays occupy contiguous memory. This happens for all types of NRVector and NRMatrix, but not for vectors of NRMatrices. (In my example x3[4][0][0] does not occupy the memory location immediately after x3[3][18][20].) See Footnote.

It's-Not-Completely-Hopeless answer:
Perhaps you can study the code for the NRMat3D class in nr3.h and use that as a starting point for developing your very own NRMat4D class that allocates contiguous storage and sets up arrays of pointers that will allow access with indexed notation.


Regards,

Dave

Footnote:
Remember that for multidimensional arrays (no matter whether they are 2-D or higher orders and no matter how you get them) the Matlab and GNU Octave (and Fortran) storage convention is column-major whereas C and C++ store arrays in row-major order, so you will have to do some transposition, somewhere, in order to use them.
For elucidation, see: http://en.wikipedia.org/wiki/Row-major_order