先介绍一下matlab与c混合编程

主要步骤：

• 使用c语言编写函数
• 利用mexFunction()函数创建C与matlab接口
• 从Matlab中编译函数

# include <mex.h>: 为了实现matlab与c混合编程的接口

//例如实现一个 add函数在接口#include “mex.h”double add(double x, double y){return x+y;}

这个代码是算法真正实现的地方。

然后创建接口。mex函数库中的mexFunction()函数，相当于c语言中的main()函数。mex源文件没有main，而是以一个mexFunction()代替。

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray * prhs[]){double  a,b;a=mxGetScalar(prhs[0]);b=mxGetScalar(prhs[1]);plhs[0]=mxCreateDoubleMatrix(1,1,mxREAL);double *c;c= mxGetPr(plhs[0]);*c=add(a,b);}

其中，nlhs为输出函数的个数，nrhs为输入参数的个数。plhs[] 和prhs[] 都是指针数组，每个元素为一个指针，指针指向的类型为mxArray. mxArray 为matlab中的矩阵。

prhs[],为输入参数， prhs[0] 指向a, prhs[1]指向b。不能用*(prhs[0])得到a的值。因为指向mxArray类型。把值以一个常用的数值形式赋给一个标量（Scalar）.

plhs[] 为输出的参数，c应该是以mxArray的类型出现的，为了得到mxArray类型的输出量，要使用 mxCreateDoubleMatrix()函数，它创建一个指向mxArray类型的指针。参数(1, 1, mxREAL)定义了对应c的尺寸类型，MATLAB中c是以1×1的实数矩阵形式报保存的。而使用mxGetPr()函数可以得到plhs[0]指向的mxArray的第一个double类型的指针。

编译add.c 为mex 文件: mex add.c

 

正式结束 edgeBox toolbox：

首先第一个函数 edgesNmsMex.cpp

/******************************************************************************** Structured Edge Detection Toolbox      Version 3.01* Code written by Piotr Dollar, 2014.* Licensed under the MSR-LA Full Rights License [see license.txt]*******************************************************************************/#include 
     
      #include 
      
       #ifdef USEOMP#include 
       
        #endif// return I[x,y] via bilinear interpolationinline float interp( float *I, int h, int w, float x, float y ) {  x = x<0 ? 0 : (x>w-1.001 ? w-1.001 : x);  y = y<0 ? 0 : (y>h-1.001 ? h-1.001 : y);  int x0=int(x), y0=int(y), x1=x0+1, y1=y0+1;  float dx0=x-x0, dy0=y-y0, dx1=1-dx0, dy1=1-dy0;  return I[x0*h+y0]*dx1*dy1 + I[x1*h+y0]*dx0*dy1 +    I[x0*h+y1]*dx1*dy0 + I[x1*h+y1]*dx0*dy0;}// E = mexFunction(E,O,r,s,m,nThreads)void mexFunction( int nl, mxArray *pl[], int nr, const mxArray *pr[] ){  float *E0 = (float*) mxGetData(pr[0]);    // original edge map  float *O = (float*) mxGetData(pr[1]);     // orientation map  int r = (int) mxGetScalar(pr[2]);         // radius for nms supr  int s = (int) mxGetScalar(pr[3]);         // radius for supr boundaries  float m = (float) mxGetScalar(pr[4]);     // multiplier for conservative supr  int nThreads = (int) mxGetScalar(pr[5]);  // number of threads for evaluation  int h=(int) mxGetM(pr[0]), w=(int) mxGetN(pr[0]);  pl[0] = mxCreateNumericMatrix(h,w,mxSINGLE_CLASS,mxREAL);  float *E = (float*) mxGetData(pl[0]);  // suppress edges where edge is stronger in orthogonal direction  #ifdef USEOMP  nThreads = nThreads
        
         w/2?w/2:s; s=s>h/2? h/2:s;  for( int x=0; x
        
       
      
     

edgeBoxesMex.cpp

/******************************************************************************** Structured Edge Detection Toolbox      Version 3.01* Code written by Piotr Dollar and Larry Zitnick, 2014.* Licensed under the MSR-LA Full Rights License [see license.txt]*******************************************************************************/#include "mex.h"#include "math.h"#include 
     
      #include 
      
       using namespace std;#define PI 3.14159265fint clamp( int v, int a, int b ) { return v
       
        b?b:v; }// trivial array class encapsulating pointer arrays //定义一个Array的类。 template 
        
          class Array{public:  Array() { _h=_w=0; _x=0; _free=0; }  virtual ~Array() { clear(); }  void clear() { if(_free) delete [] _x; _h=_w=0; _x=0; _free=0; }  void init(int h, int w) { clear(); _h=h; _w=w; _x=new T[h*w](); _free=1; }//初始化矩阵  T& val(size_t c, size_t r) { return _x[c*_h+r]; }//取矩阵中的值val(i,j)  int _h, _w; T *_x; bool _free;};// convenient typedefstypedef vector
         
           vectorf;typedef vector
          
            vectori;typedef Array
           
             arrayf;typedef Array
            
              arrayi;// bounding box data structures and routinestypedef struct { int c, r, w, h; float s; } Box;typedef vector
             
               Boxes;bool boxesCompare( const Box &a, const Box &b ) { return a.s
              
                _segAff; // segment affinities vector
               
                 _segAffIdx; // segment neighbors // data structures for efficiency (see prepDataStructs) arrayf _segIImg, _magIImg; arrayi _hIdxImg, _vIdxImg; vector
                
                  _hIdxs, _vIdxs; vectorf _scaleNorm; float _scStep, _arStep, _rcStepRatio; // data structures for efficiency (see scoreBox) arrayf _sWts; arrayi _sDone, _sMap, _sIds; int _sId; // helper routines void clusterEdges( arrayf &E, arrayf &O, arrayf &V ); void prepDataStructs( arrayf &E ); void scoreAllBoxes( Boxes &boxes ); void scoreBox( Box &box ); void refineBox( Box &box ); void drawBox( Box &box, arrayf &E, arrayf &V );};////
                
               
              
             
            
           
          
         
        
       
      
     

//在mexFunction中定义一个EdgeBoxGenerator对象edgeBoxGen和Boxes对象EdgeBoxGenerator edgeBoxGen; Boxes boxes;

//调用edgeBoxGen对象的generat函数：  edgeBoxGen.generate( boxes, E, O, V );

void EdgeBoxGenerator::generate( Boxes &boxes, arrayf &E, arrayf &O, arrayf &V ){  clusterEdges(E,O,V); prepDataStructs(E); scoreAllBoxes(boxes);}// 先将edge聚类为edge groupvoid EdgeBoxGenerator::clusterEdges( arrayf &E, arrayf &O, arrayf &V ){
    int c, r, cd, rd, i, j; h=E._h; w=E._w;  // greedily merge connected edge pixels into clusters (create _segIds)  _segIds.init(h,w); _segCnt=1; //初始化一个种子矩阵用定义的Array class，将edge梯度模值大于阈值的点作为seed。  for( c=0; c
     
      .5) v=1-v;          vs.push_back(v); cs.push_back(c0+cd); rs.push_back(r0+rd);      }      float minv=1000; j=0;      for( i=0; i
      
       0 ) _segMag[j]+=E.val(c,r);  for( c=1; c
       
        0 && _segMag[j]<=_clusterMinMag)      _segIds.val(c,r)=0;  i=1; while(i>0) {    i=0; for( c=1; c
        
         .5) v=1-v;        if( v
         
          0) i++;    }  }  // compactify representation //上一步合并之后 _segCnt 的值不是0,1,2,3，。。。， 中间有些点被合并，重新从0,1,2,3,...排起。  _segMag.assign(_segCnt,0); vectori map(_segCnt,0); _segCnt=1;  for( c=1; c
          
           0 ) _segMag[j]+=E.val(c,r); for( i=0; i<_segMag.size(); i++ ) if( _segMag[i]>0 ) map[i]=_segCnt++; for( c=1; c
           
            0 ) _segIds.val(c,r)=map[j];//重新计算梯度方向与位置。 // compute positional means and recompute _segMag _segMag.assign(_segCnt,0); vectorf meanX(_segCnt,0), meanY(_segCnt,0); vectorf meanOx(_segCnt,0), meanOy(_segCnt,0), meanO(_segCnt,0); for( c=1; c
            
             0 ) { float m=_segMag[i]; meanX[i]/=m; meanY[i]/=m; meanO[i]=atan2(meanOy[i]/m,meanOx[i]/m)/2; } // compute segment affinities //计算分割的edge group 之间的相似度在-2 到2 的范围内， 每个group有一个相似groupindex 的stack和 一个记录对应的相似group相似值得stack。 _segAff.resize(_segCnt); _segAffIdx.resize(_segCnt); for(i=0; i<_segCnt; i++) _segAff[i].resize(0); for(i=0; i<_segCnt; i++) _segAffIdx[i].resize(0); const int rad = 2; for( c=rad; c
             
              0 ) { _segC[j]=c; _segR[j]=r; } // optionally create visualization (assume memory initialized is 3*w*h) if( V._x ) for( c=0; c
              
               0 ) { E1.val(_segC[i],_segR[i]) = _segMag[i]; } _segIImg.init(h+1,w+1); for( c=1; c
               
                 _edgeMinMag ? E.val(c,r) : 0; _magIImg.val(c+1,r+1) = e + _magIImg.val(c,r+1) + _magIImg.val(c+1,r) - _magIImg.val(c,r); } // create remaining data structures _hIdxs.resize(h); _hIdxImg.init(h,w); for( r=0; r
                
                 0 && sDone[j]!=sId ) { sIds[n]=j; sWts[n]=1; sDone[j]=sId; sMap[j]=n++; } cs=_hIdxImg.val(c0,r1); ce=_hIdxImg.val(c1,r1); // bottom for( i=cs; i<=ce; i++ ) if( (j=_hIdxs[r1][i])>0 && sDone[j]!=sId ) { sIds[n]=j; sWts[n]=1; sDone[j]=sId; sMap[j]=n++; } rs=_vIdxImg.val(c0,r0); re=_vIdxImg.val(c0,r1); // left for( i=rs; i<=re; i++ ) if( (j=_vIdxs[c0][i])>0 && sDone[j]!=sId ) { sIds[n]=j; sWts[n]=1; sDone[j]=sId; sMap[j]=n++; } rs=_vIdxImg.val(c1,r0); re=_vIdxImg.val(c1,r1); // right for( i=rs; i<=re; i++ ) if( (j=_vIdxs[c1][i])>0 && sDone[j]!=sId ) { sIds[n]=j; sWts[n]=1; sDone[j]=sId; sMap[j]=n++; } // follow connected paths and set weights accordingly (w=1 means remove) for( i=0; i
                 
                  sWts[sMap[q]] ) { sWts[sMap[q]]=wq; i=min(i,sMap[q]-1); } } else if(_segC[q]>=c0 && _segC[q]<=c1 && _segR[q]>=r0 && _segR[q]<=r1) { sIds[n]=q; sWts[n]=wq; sDone[q]=sId; sMap[q]=n++; } } } // finally remove segments connected to boundaries for( i=0; i
                  
                   =c0 && _segC[k]<=c1 && _segR[k]>=r0 && _segR[k]<=r1 ) v -= sWts[i]*_segMag[k]; } v*=norm; if(v<_minScore) v=0; box.s=v;}void EdgeBoxGenerator::refineBox( Box &box ){ int rStep = int(box.h*_rcStepRatio); int cStep = int(box.w*_rcStepRatio); while( 1 ) { // prepare for iteration rStep/=2; cStep/=2; if( rStep<=2 && cStep<=2 ) break; rStep=max(1,rStep); cStep=max(1,cStep); Box B; // search over r start B=box; B.r=box.r-rStep; B.h=B.h+rStep; scoreBox(B); if(B.s<=box.s) { B=box; B.r=box.r+rStep; B.h=B.h-rStep; scoreBox(B); } if(B.s>box.s) box=B; // search over r end B=box; B.h=B.h+rStep; scoreBox(B); if(B.s<=box.s) { B=box; B.h=B.h-rStep; scoreBox(B); } if(B.s>box.s) box=B; // search over c start B=box; B.c=box.c-cStep; B.w=B.w+cStep; scoreBox(B); if(B.s<=box.s) { B=box; B.c=box.c+cStep; B.w=B.w-cStep; scoreBox(B); } if(B.s>box.s) box=B; // search over c end B=box; B.w=B.w+cStep; scoreBox(B); if(B.s<=box.s) { B=box; B.w=B.w-cStep; scoreBox(B); } if(B.s>box.s) box=B; }}void EdgeBoxGenerator::drawBox( Box &box, arrayf &E, arrayf &V ){ // score box and draw color coded edges (red=out, green=in) int i, c, r; float e, o; if( !V._x ) return; int sId=_sId; scoreBox(box); int c0, r0, c1, r1; r1=clamp(box.r+box.h,0,h-1); r0=box.r=clamp(box.r,0,h-1); c1=clamp(box.c+box.w,0,w-1); c0=box.c=clamp(box.c,0,w-1); for( c=0; c
                   
                    =c0 && _segC[i]<=c1 && _segR[i]>=r0 && _segR[i]<=r1 ) ? 0 : 1; V.val(c+w*0,r)=1-e+e*o; V.val(c+w*1,r)=1-e*o; V.val(c+w*2,r)=1-e; } // finally draw bounding box r=r0; for(c=c0; c<=c1; c++) V.val(c+w*0,r)=V.val(c+w*1,r)=V.val(c+w*2,r)=0; r=r1; for(c=c0; c<=c1; c++) V.val(c+w*0,r)=V.val(c+w*1,r)=V.val(c+w*2,r)=0; c=c0; for(r=r0; r<=r1; r++) V.val(c+w*0,r)=V.val(c+w*1,r)=V.val(c+w*2,r)=0; c=c1; for(r=r0; r<=r1; r++) V.val(c+w*0,r)=V.val(c+w*1,r)=V.val(c+w*2,r)=0;}void EdgeBoxGenerator::scoreAllBoxes( Boxes &boxes ){ // get list of all boxes roughly distributed in grid boxes.resize(0); int arRad, scNum; float minSize=sqrt(_minBoxArea); arRad = int(log(_maxAspectRatio)/log(_arStep*_arStep)); scNum = int(ceil(log(max(w,h)/minSize)/log(_scStep))); for( int s=0; s
                    
                     =r1i || a.c>=c1i ) return 0; r1j=b.r+b.h; c1j=b.c+b.w; if( a.r>=r1j || a.c>=c1j ) return 0; areai = (float) a.w*a.h; r0=max(a.r,b.r); r1=min(r1i,r1j); areaj = (float) b.w*b.h; c0=max(a.c,b.c); c1=min(c1i,c1j); areaij = (float) max(0,r1-r0)*max(0,c1-c0); return areaij / (areai + areaj - areaij);}void boxesNms( Boxes &boxes, float thr, float eta, int maxBoxes ){ sort(boxes.rbegin(),boxes.rend(),boxesCompare); if( thr>.99 ) return; const int nBin=10000; const float step=1/thr, lstep=log(step); vector
                     
                       kept; kept.resize(nBin+1); int i=0, j, k, n=(int) boxes.size(), m=0, b, d=1; while( i
                      
                       .5) { thr*=eta; d=ceil(log(1/thr)/lstep); } } boxes.resize(m); i=0; for( j=0; j
                      
                     
                    
                   
                  
                 
                
               
              
             
            
           
          
         
        
       
      
     

void mexFunction( int nl, mxArray *pl[], int nr, const mxArray *pr[] ){  // check and get inputs  if(nr != 14) mexErrMsgTxt("Fourteen inputs required.");  if(nl > 2) mexErrMsgTxt("At most two outputs expected.");  if(mxGetClassID(pr[0])!=mxSINGLE_CLASS) mexErrMsgTxt("E must be a float*");  if(mxGetClassID(pr[1])!=mxSINGLE_CLASS) mexErrMsgTxt("O must be a float*");  arrayf E; E._x = (float*) mxGetData(pr[0]);  arrayf O; O._x = (float*) mxGetData(pr[1]);  int h = (int) mxGetM(pr[0]); O._h=E._h=h;  int w = (int) mxGetN(pr[0]); O._w=E._w=w;  // optionally create memory for visualization  arrayf V; if( nl>1 ) {    const int ds[3] = {h,w,3};    pl[1] = mxCreateNumericArray(3,ds,mxSINGLE_CLASS,mxREAL);    V._x = (float*) mxGetData(pl[1]); V._h=h; V._w=w;  }  // setup and run EdgeBoxGenerator  EdgeBoxGenerator edgeBoxGen; Boxes boxes;  edgeBoxGen._alpha = float(mxGetScalar(pr[2]));  edgeBoxGen._beta = float(mxGetScalar(pr[3]));  edgeBoxGen._eta = float(mxGetScalar(pr[4]));  edgeBoxGen._minScore = float(mxGetScalar(pr[5]));  edgeBoxGen._maxBoxes = int(mxGetScalar(pr[6]));  edgeBoxGen._edgeMinMag = float(mxGetScalar(pr[7]));  edgeBoxGen._edgeMergeThr = float(mxGetScalar(pr[8]));  edgeBoxGen._clusterMinMag = float(mxGetScalar(pr[9]));  edgeBoxGen._maxAspectRatio = float(mxGetScalar(pr[10]));  edgeBoxGen._minBoxArea = float(mxGetScalar(pr[11]));  edgeBoxGen._gamma = float(mxGetScalar(pr[12]));  edgeBoxGen._kappa = float(mxGetScalar(pr[13]));  edgeBoxGen.generate( boxes, E, O, V );  // create output bbs and output to Matlab  int n = (int) boxes.size();  pl[0] = mxCreateNumericMatrix(n,5,mxSINGLE_CLASS,mxREAL);  float *bbs = (float*) mxGetData(pl[0]);  for(int i=0; i