Mex Code for Array Slicing

Hello,

Thank you for your help.

I have a Matlab 3D array (buildSpaceGrid) 300 * 300 * 600. During my algorithm, i need to access subarrays of it several times. It is very inefficient to get such arrays with simple array slicing.

For example, i tried this, where build Space Grid is my 3D array, ixx,iyy, and izz are the starting x,y and z coordinates in the 3 dimensions.

buildSpaceGrid( ixx:ixx+partSize(1)-1,...
            iyy:iyy+partSize(2)-1,...
            izz:izz+partSize(3)-1)   

It seems doing this type of array slicing is like death sentence in Matlab in terms of performance. I have to do this several times during my algorithm for different values of ixx, iyy and izz while my 3D array also updates as the algorithm progresses.

I decided to write a mex routine in c++ to make this more efficient. The routine takes the 3D array, buildSpaceGrid, ixx,iyy,izz and partsize and outputs the subarray.

#include "mex.h"
#include "matrix.h"

void mexFunction(
	int nlhs, mxArray *plhs[],
	int nrhs, const mxArray *prhs[])
{
	mwSize nDimNum;
	const mwSize *pSize;
	bool    *pD;
	mxArray   *Data;
	bool *buildSpaceGrid;
	int *ixx;
	int *iyy;
	int *izz;
	mwSize i, j, k;
	buildSpaceGrid = mxGetLogicals(prhs[0]);
	ixx = (int *)mxGetPr(prhs[1]);
	iyy = (int *)mxGetPr(prhs[1]);
	izz = (int *)mxGetPr(prhs[3]);
	pSize = (mwSize *)mxGetPr(prhs[4]);

	nDimNum = mxGetNumberOfDimensions(prhs[0]);
	plhs[0] = mxCreateLogicalArray(nDimNum, pSize);
	pD = (bool *)mxGetData(plhs[0]);


	for (i = 0; i < pSize[0]; i++)    // NOT i=i !
	{
		for (j = 0; j < pSize[1]; j++)
		{
			for (k = 0; k < pSize[2]; k++)
			{
				pD[i+j*pSize[0]+k*(pSize[0]*pSize[1])] = buildSpaceGrid[ixx-1+i+(iyy-1+j)*pSize[0]+(izz-1+k)*(pSize[0]*pSize[1])];
			}
		}
	}
}

	

I keep getting this wierd error when I try to compile it using mex command in Matlab.

error C2296: '*': illegal, left operand has type 'int *'

Can someone help me fix this problem.

I was also wondering if there is a better (effiient) way of improving the performance of this line of code.

Thank you so much for your help.

Answer 1

[James Tursa]

You are not going to get any performance boost for this in a mex routine. Since the m-code and the mex code both do deep data copying, and since this action represents the bulk of the work, the mex routine would not be expected to do any better than the m-code for performance. Where a mex routine can help you is if you could do the downstream processing inside your mex routine to avoid the deep data copy. You don't appear to have that situation.

What line has the error? And what is the complete error message?

Answer 2

[Maaz Saleem Kapadia]

Hello James,

Thank you so much for your response. I was able to figure it out. It was giving an error related to variable declarations.

Here is the final code that worked (just for reference)

#include "mex.h"
#include "matrix.h"

void mexFunction(
	int nlhs, mxArray *plhs[],
	int nrhs, const mxArray *prhs[])
{
	mwSize nDimNum;
	const mwSize *pSize;
	const mwSize *bSize;
	bool    *pD;
	mxArray   *Data;
	bool *buildSpaceGrid;
	double ixx;
	double iyy;
	double izz;
	mwSize s1, s2, s3;
	buildSpaceGrid = mxGetLogicals(prhs[0]);
	ixx = mxGetScalar(prhs[1]);
	iyy = mxGetScalar(prhs[2]);
	izz = mxGetScalar(prhs[3]);
	pSize = mxGetDimensions(prhs[4]);
	bSize = mxGetDimensions(prhs[0]);
	int idx1;
	int idx2;

	nDimNum = mxGetNumberOfDimensions(prhs[0]);
	plhs[0] = mxCreateLogicalArray(nDimNum, pSize);
	pD = (bool *)mxGetData(plhs[0]); 
	s1 = pSize[0];
	s2 = pSize[1];
	s3 = pSize[2];

	for (int i = 0; i < s1; i++)    // NOT i=i !
	{
		for (int j = 0; j < s2; j++)
		{
			for (int k = 0; k < s3; k++)
			{
				idx1 = i + j * pSize[0] + k * (pSize[0] * pSize[1]);
				idx2 = ixx - 1 + i + (iyy - 1 + j)*bSize[0] + (izz - 1 + k)*(bSize[0] * bSize[1]);
				pD[idx1] = buildSpaceGrid[idx2];
			}
		}
	}
}

	

It worked as i intended but it did not improve the efficiency of the code (as you mentioned above). I guess, i would have to stick with the Matlab version of Array slicing.

Just wondering if do you think, i could somehow make it faster since it is the only bottleneck in my Matlab code.

Thank you so much for help and time.

[James Tursa]

The only advice I can offer is what I had mentioned earlier. If you can move the downstream processing with that slice into the mex routine then you can avoid the up-front data copy which will save you some time. This may or may not be easy or feasible depending on what that downstream code actually does with the slice.

[Maaz Saleem Kapadia]

Thank you so much for your help. Based on your advice, I decided to move the time-consuming downstream processing within the mex routine.

The purpose of this code is just to find the first overlap in the two 3D arrays, buildSpace Grid and VoxelGrid. First the buildSpaceGrid dimensions are made consistent with voxelGrid dimensions and then the overlap index is calculated.

in the mex routine below, instead of first calculating the sliced buildSpaceGrid and then checking the overlap, the overlap is checked in between and the loop breaks whenever the first overlap is found. So this code, in my understanding, is avoiding the full deep copy of buildSpaceGrid and directly checking for the overlap.

This code works but it turns out it is more inefficient than my original version.

In my original version, the deep copy of buildSpaceGrid was done using matlab and the downstream processing to find the overlap was done in a separate mex routine.

Do you see any posibility of improving this mex routine. I dont understand why this is even more inefficient than my original version.

Here is the code of it.

#include "mex.h"
#include "matrix.h"

void mexFunction(
	int nlhs, mxArray *plhs[],
	int nrhs, const mxArray *prhs[])
{
	const mwSize *pSize;
	const mwSize *bSize;
	double temp;
	bool *buildSpaceGrid;
	bool *voxelGrid;
	double ixx;
	double iyy;
	double izz;
	double *out;
	mwSize s1, s2, s3;
	buildSpaceGrid = mxGetLogicals(prhs[0]);
	voxelGrid = mxGetLogicals(prhs[4]);
	ixx = mxGetScalar(prhs[1]);
	iyy = mxGetScalar(prhs[2]);
	izz = mxGetScalar(prhs[3]);
	pSize = mxGetDimensions(prhs[4]);
	bSize = mxGetDimensions(prhs[0]);

	int idx1;
	int idx2;
	int i, j, k;
	plhs[0] = mxCreateDoubleMatrix(1, 4, mxREAL);
	s1 = pSize[0];
	s2 = pSize[1];
	s3 = pSize[2];

	out = mxGetPr(plhs[0]);

	out[0] = 1;
	out[1] = 1;
	out[2] = 1;
	out[3] = 1;
	for (int c = 0; c < s1 * s2 * s3;c++) {
		i = ((c) % (s1 * s2)) % s1;
		j = ((c) % (s1 * s2)) / s1;
		k = (c) / (s1 * s2);
		idx1 = i + j * s1 + k * (s1 * s2);
		idx2 = ixx - 1 + i + (iyy - 1 + j)*bSize[0] + (izz - 1 + k)*(bSize[0] * bSize[1]);
		temp = voxelGrid[idx1] + buildSpaceGrid[idx2];
		if (temp > 1.5) {

			out[0] = 0;
			out[1] = i+1;
			out[2] = j+1;
			out[3] = k + 1;
			break;
		}
		
	}