Matlab 和 C++ （1）

recbio

用了Matlab 为什么还要C++？
主要的问题是：
1）速度。
2）为了以后可以独立于Matlab
最近一直在用Matlab处理实验的图形。碰到一个问题就是速度。我们的图形是16bit灰度，一般为512（长）*512（宽）*300 - 1000（幅）的一个三维矩阵；300-1000是因为有300-1000个时间点和深度。为了内存的关系，我们不得不用short int来保存图形。图形处理中的一步是得到这些图的平均值图。一般用Matlab 的mean(pic_Array, 3)来完成。可是非常慢。所以，就用C写了一个函数。顺便贴在这里，大如果有同样的要求，可以用用。
函数叫mymean3意思是求第三维的平均。对于16bit 矩阵，速度可以提高10倍以上。
执行结果：
>>A=uint16(rand(200,200,1000)*1000);       (这一步是产生一个随机三维16bit矩阵)
>>tic;B=mean(A,3);toc                      (这一步是Matlab 求平均)
Elapsed time is 4.087521 seconds.
>> tic;C=mymean3(A);toc                    (这一步我们自己的函数求平均，怎么样快了约一个半数量级！)
Elapsed time is 0.129360 seconds.
>>sum(sum(B-C,2),1)        (这一步是一个简单的检验，证明结果和matlab一致)

ans =

0

------------------------------------------------------------
// file name: myMean3.cpp
#include "mex.h"
#include "math.h"
#define MY_START 0
//#include <windows.h>
#define MY_MAX_DIM  10

//#define MATLAB_CLASS_NAME "com.mathworks.mde.desk.MLMainFrame"
//#define MATLAB_WND_NAME "MATLAB"


/*----------------------------------------------------------------------------------------------*/
int mySqueeze(const mxArray *pInArray, int *pOutDim, int nMax){
int i,j,k,l=0;
const int *pAddr = mxGetDimensions(pInArray);
char lpString[100];

j=mxGetNumberOfDimensions(pInArray);
for(i=0;i<j;i++){
    if((k=pAddr)==1){
        //sprintf(lpString,"Dimention %d is 1, it will be queezed !",i+1);   //GetTopWindow(NULL)
        //MessageBox( FindWindow(MATLAB_CLASS_NAME, MATLAB_WND_NAME),lpString, "Squeeze One Dimention !", MB_OK|MB_ICONINFORMATION);
    }else{
        if(l<nMax)l++;
        else break;
        pOutDim[l]=k;
    }
}
return(pOutDim[0]=l);
}
/*----------------------------------------------------------------------------------------------*/
void myCalMean3(unsigned short *pInArray, int nM, int nN, int nL, double *dOutArray){
int nPlaneL= nM * nN, nSum, *pnSum, *pnSumA, i,j;
unsigned short *pB=pInArray;
double *pdOut=dOutArray;

pnSumA = pnSum = (int *)malloc(nPlaneL * sizeof (int));
if(pnSum==NULL)mexErrMsgTxt("Out of Memory");
for(i=0;i<nPlaneL;i++) *pnSumA++ = *pB++;

for(j=1;j<nL;j++){
   pnSumA = pnSum;
   for(i=0;i<nPlaneL;i++) *pnSumA++ += *pB++;
}
pnSumA = pnSum;
for(i=0;i<nPlaneL;i++)
  *pdOut++ = (double)(*pnSumA++)/nL;
free(pnSum);
}
/*----------------------------------------------------------------------------------------------*/
void myCalMean3(unsigned int *pInArray, int nM, int nN, int nL, double *dOutArray){
int nPlaneL=nM*nN;
int nSum,i,j;
unsigned int *pB=pInArray;
double *pdSum, *pdSumA, *pdOut=dOutArray;

pdSumA = pdSum = (double *) malloc(nPlaneL * sizeof (double));
if(pdSum==NULL)mexErrMsgTxt("Out of Memory");
for(i=0;i<nPlaneL;i++) *pdSumA++ = *pB++;

for(j=1;j<nL;j++){
   pdSumA = pdSum;
   for(i=0;i<nPlaneL;i++) *pdSumA++ += *pB++;
}
pdSumA = pdSum;
for(i=0;i<nPlaneL;i++)
  *pdOut++ = (double)(*pdSumA++)/nL;
free(pdSum);
}
/*----------------------------------------------------------------------------------------------*/
void myCalMean3(double *pInArray, int nM, int nN, int nL, double *dOutArray){
int nPlaneL=nM*nN;
int nSum,i,j;
double *pdSum, *pdSumA, *pdOut=dOutArray, *pB=pInArray;

pdSumA = pdSum = (double *)malloc(nPlaneL * sizeof (double));
if(pdSum==NULL)mexErrMsgTxt("Out of Memory");
for(i=0;i<nPlaneL;i++) *pdSumA++ = *pB++;

for(j=1;j<nL;j++){
   pdSumA = pdSum;
   for(i=0;i<nPlaneL;i++) *pdSumA++ += *pB++;
}
pdSumA = pdSum;
for(i=0;i<nPlaneL;i++)
  *pdOut++ = (double)(*pdSumA++)/nL;
free(pdSum);
}
/*----------------------------------------------------------------------------------------------*/
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
int nArrayID, nArrayDim, nFun, nSize;
int i=0,j;
int pnDimCount[MY_MAX_DIM+1];
int pnDimOut[MY_MAX_DIM];

  /*  Check for proper number of arguments. */
  /* NOTE: You do not need an else statement when using
     mexErrMsgTxt within an if statement. It will never
     get to the else statement if mexErrMsgTxt is executed.
     (mexErrMsgTxt breaks you out of the MEX-file.)
  */
  if (nrhs!=1)
    mexErrMsgTxt("Using: aOut = myMean3( aIn ); aIn is 3 DIM matrix ");
  nArrayID = mxGetClassID(prhs[0]);
  nArrayDim = mySqueeze(prhs[0], pnDimCount, MY_MAX_DIM);
  for(i=0,j=1;j<=nArrayDim;j++){
    if(j==3)continue;
    pnDimOut[i++]= pnDimCount[j];
  }
  //printf("i=%d, %d, %d, %d, %d\n",i,pnDimOut[0], pnDimOut[1], pnDimOut[2], pnDimOut[3]);  // for error check only
  if(nArrayDim <3){
    mexErrMsgTxt("aIn must be an array larger than 2 dimension" );
  }
  //plhs[0] = mxCreateNumericArray( i , pnDimOut, mxDOUBLE_CLASS, mxREAL); // for array > 3
  plhs[0] = mxCreateNumericArray( 2 , pnDimOut, mxDOUBLE_CLASS, mxREAL);  // for array ==3
  switch(nArrayID){
    case mxINT16_CLASS:
    case mxUINT16_CLASS:
        myCalMean3((unsigned short *)mxGetData (prhs[0]), pnDimCount[1], pnDimCount[2], pnDimCount[3], mxGetPr(plhs[0]));
        break;
    case mxINT32_CLASS:
    case mxUINT32_CLASS:
        myCalMean3((unsigned int *)mxGetData (prhs[0]), pnDimCount[1], pnDimCount[2], pnDimCount[3], mxGetPr(plhs[0]));
        break;
    case mxDOUBLE_CLASS:
        myCalMean3((double *)mxGetData (prhs[0]), pnDimCount[1], pnDimCount[2], pnDimCount[3], mxGetPr(plhs[0]));
        break;
    default:
        printf("ClassID= %d, Name= %s, is not supported now\n", nArrayID, mxGetClassName(prhs[0]));
        mexErrMsgTxt("aIn must be INT or Double");
        break;

}

}
---------------------------------

注，以上默认的是LCC，您可以用 -o2来优化代码。
如果是用VC 或者BC编译，打"//"的语句就可以用了。这样，可以给您的c函数阿增加一个squeeze的窗口。详细可以看函数中打"//"的各行。这里的matlab class name 对应的是matlab 2006ab 2007a 的class name 主窗口是一个 com.mathworks.mde.desk.MLMainFrame class

eisenstange

曾经看过Mex和C的接口，不过没有仔细研究。多谢楼主！