edge box

//例如实现一个 add函数在接口

#include “mex.h”
double add(double x, double y)
{
return x+y;
}

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);
}

/*******************************************************************************
* 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 <mex.h>
#include <math.h>
#ifdef USEOMP
#include <omp.h>
#endif

// return I[x,y] via bilinear interpolation
inline 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<omp_get_max_threads() ? nThreads : omp_get_max_threads();
  #pragma omp parallel for num_threads(nThreads)
  #endif
  for( int x=0; x<w; x++ ) for( int y=0; y<h; y++ ) {
    float e=E[x*h+y]=E0[x*h+y]; if(!e) continue; e*=m;
    float coso=cos(O[x*h+y]), sino=sin(O[x*h+y]);
    for( int d=-r; d<=r; d++ ) if( d ) {
    //得到正交方向上的梯度，用插值法得到
      float e0 = interp(E0,h,w,x+d*coso,y+d*sino);
   //如果正交方向上的梯度大于边界值，说明该点在各个方向上都变化很大为独立店，除去此边界值。
      if(e < e0) { E[x*h+y]=0; break; }
    }
  }

  // suppress noisy edge estimates near boundaries
//给边缘点赋值

  s=s>w/2?w/2:s; s=s>h/2? h/2:s;
  for( int x=0; x<s; x++ ) for( int y=0; y<h; y++ ) {
    E[x*h+y]*=x/float(s); E[(w-1-x)*h+y]*=x/float(s); }
  for( int x=0; x<w; x++ ) for( int y=0; y<s; y++ ) {
    E[x*h+y]*=y/float(s); E[x*h+(h-1-y)]*=y/float(s); }
}

/*******************************************************************************
* 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 <algorithm>
#include <vector>
using namespace std;
#define PI 3.14159265f
int clamp( int v, int a, int b ) { return v<a?a:v>b?b:v; }

// trivial array class encapsulating pointer arrays
//定义一个Array的类。
template <class T> 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 typedefs
typedef vector<float> vectorf;
typedef vector<int> vectori;
typedef Array<float> arrayf;
typedef Array<int> arrayi;

// bounding box data structures and routines
typedef struct { int c, r, w, h; float s; } Box;
typedef vector<Box> Boxes;
bool boxesCompare( const Box &a, const Box &b ) { return a.s<b.s; }
float boxesOverlap( Box &a, Box &b );
void boxesNms( Boxes &boxes, float thr, float eta, int maxBoxes );

// main class for generating edge boxes
class EdgeBoxGenerator
{
public:
  // method parameters (must be manually set)
  float _alpha, _beta, _eta, _minScore; int _maxBoxes;
  float _edgeMinMag, _edgeMergeThr, _clusterMinMag;
  float _maxAspectRatio, _minBoxArea, _gamma, _kappa;

  // main external routine (set parameters first)
  void generate( Boxes &boxes, arrayf &E, arrayf &O, arrayf &V );

private:
  // edge segment information (see clusterEdges)
  int h, w;                         // image dimensions
  int _segCnt;                      // total segment count
  arrayi _segIds;                   // segment ids (-1/0 means no segment)
  vectorf _segMag;                  // segment edge magnitude sums
  vectori _segR, _segC;             // segment lower-right pixel
  vector<vectorf> _segAff;          // segment affinities
  vector<vectori> _segAffIdx;       // segment neighbors

  // data structures for efficiency (see prepDataStructs)
  arrayf _segIImg, _magIImg; arrayi _hIdxImg, _vIdxImg;
  vector<vectori> _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 group
void 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<w; c++ ) for( r=0; r<h; r++ ) {
    if( c==0 || r==0 || c==w-1 || r==h-1 || E.val(c,r)<=_edgeMinMag )
      _segIds.val(c,r)=-1; else _segIds.val(c,r)=0;
  }
  for( c=1; c<w-1; c++ ) for( r=1; r<h-1; r++ ) {
    if(_segIds.val(c,r)!=0) continue;
    float sumv=0; int c0=c, r0=r; vectorf vs; vectori cs, rs;
    while( sumv < _edgeMergeThr ) {
      //将种子点分为不同的类，类的标号为_segCnt
      _segIds.val(c0,r0)=_segCnt;
      float o0 = O.val(c0,r0), o1, v; bool found;
       //遍历每个种子点。
      for( cd=-1; cd<=1; cd++ ) for( rd=-1; rd<=1; rd++ ) {
       //遍历种子点周围四邻域的种子点，如果存在种子点加入栈中。
        if( _segIds.val(c0+cd,r0+rd)!=0 ) continue; found=false;
        for( i=0; i<cs.size(); i++ )
        //点已被加入栈中。
          if( cs[i]==c0+cd && rs[i]==r0+rd ) { found=true; break; }
          if( found ) continue; o1=O.val(c0+cd,r0+rd);
          v=fabs(o1-o0)/PI; if(v>.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<vs.size(); i++ ) if( vs[i]<minv ) {
        minv=vs[i]; c0=cs[i]; r0=rs[i]; j=i;
      }
      sumv+=minv; if(minv<1000) vs[j]=1000;
    }
    _segCnt++;
  }

  // merge or remove small segments
  _segMag.resize(_segCnt,0);
  for( c=1; c<w-1; c++ ) for( r=1; r<h-1; r++ )
    if( (j=_segIds.val(c,r))>0 ) _segMag[j]+=E.val(c,r);
  for( c=1; c<w-1; c++ ) for( r=1; r<h-1; r++ )
//如果edge group的梯度模值过小删除。
    if( (j=_segIds.val(c,r))>0 && _segMag[j]<=_clusterMinMag)
      _segIds.val(c,r)=0;
  i=1; while(i>0) {
    i=0; for( c=1; c<w-1; c++ ) for( r=1; r<h-1; r++ ) {
      if( _segIds.val(c,r)!=0 ) continue;
      float o0=O.val(c,r), o1, v, minv=1000; j=0;
      for( cd=-1; cd<=1; cd++ ) for( rd=-1; rd<=1; rd++ ) {
       //合并四邻域内的edge group.
        if( _segIds.val(c+cd,r+rd)<=0 ) continue; o1=O.val(c+cd,r+rd);
        v=fabs(o1-o0)/PI; if(v>.5) v=1-v;
        if( v<minv ) { minv=v; j=_segIds.val(c+cd,r+rd); }
      }
      _segIds.val(c,r)=j; if(j>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<w-1; c++ ) for( r=1; r<h-1; r++ )
    if( (j=_segIds.val(c,r))>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<w-1; c++ ) for( r=1; r<h-1; r++ )
    if( (j=_segIds.val(c,r))>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<w-1; c++ ) for( r=1; r<h-1; r++ ) {
    j=_segIds.val(c,r); if(j<=0) continue;
    float m=E.val(c,r), o=O.val(c,r); _segMag[j]+=m;
    meanOx[j]+=m*cos(2*o); meanOy[j]+=m*sin(2*o);
    meanX[j]+=m*c; meanY[j]+=m*r;
  }
  for( i=0; i<_segCnt; i++ ) if( _segMag[i]>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<w-rad; c++ ) for( r=rad; r<h-rad; r++ ) {
    int s0=_segIds.val(c,r); if( s0<=0 ) continue;
    for( cd=-rad; cd<=rad; cd++ ) for( rd=-rad; rd<=rad; rd++ ) {
      int s1=_segIds.val(c+cd,r+rd); if(s1<=s0) continue;
      bool found = false; for(i=0;i<_segAffIdx[s0].size();i++)
        if(_segAffIdx[s0][i] == s1) { found=true; break; }
      if( found ) continue;
      float o=atan2(meanY[s0]-meanY[s1],meanX[s0]-meanX[s1])+PI/2;
      float a=fabs(cos(meanO[s0]-o)*cos(meanO[s1]-o)); a=pow(a,_gamma);
      _segAff[s0].push_back(a); _segAffIdx[s0].push_back(s1);
      _segAff[s1].push_back(a); _segAffIdx[s1].push_back(s0);
    }
  }

  // compute _segC and _segR
  _segC.resize(_segCnt); _segR.resize(_segCnt);
  for( c=1; c<w-1; c++ ) for( r=1; r<h-1; r++ )
    if( (j=_segIds.val(c,r))>0 ) { _segC[j]=c; _segR[j]=r; }

  // optionally create visualization (assume memory initialized is 3*w*h)
  if( V._x ) for( c=0; c<w; c++ ) for( r=0; r<h; r++ ) {
    i=_segIds.val(c,r);
    V.val(c+w*0,r) = i<=0 ? 1 : ((123*i + 128)%255)/255.0f;
    V.val(c+w*1,r) = i<=0 ? 1 : ((7*i + 3)%255)/255.0f;
    V.val(c+w*2,r) = i<=0 ? 1 : ((174*i + 80)%255)/255.0f;
  }
}

void EdgeBoxGenerator::prepDataStructs( arrayf &E )
{
  int c, r, i;

  // initialize step sizes
  _scStep=sqrt(1/_alpha);
  _arStep=(1+_alpha)/(2*_alpha);
  _rcStepRatio=(1-_alpha)/(1+_alpha);

  // create _scaleNorm
  _scaleNorm.resize(10000);
  for( i=0; i<10000; i++ )
    _scaleNorm[i]=pow(1.f/i,_kappa);

// 生成分割图。
  // create _segIImg
  arrayf E1; E1.init(h,w);
  for( i=0; i<_segCnt; i++ ) if( _segMag[i]>0 ) {
    E1.val(_segC[i],_segR[i]) = _segMag[i];
  }
  _segIImg.init(h+1,w+1);
  for( c=1; c<w; c++ ) for( r=1; r<h; r++ ) {
    _segIImg.val(c+1,r+1) = E1.val(c,r) + _segIImg.val(c,r+1) +
      _segIImg.val(c+1,r) - _segIImg.val(c,r);
  }
//得到边缘图。
  // create _magIImg
  _magIImg.init(h+1,w+1);
  for( c=1; c<w; c++ ) for( r=1; r<h; r++ ) {
    float e = E.val(c,r) > _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<h; r++ ) {
    int s=0, s1; _hIdxs[r].push_back(s);
    for( c=0; c<w; c++ ) {
      s1 = _segIds.val(c,r);
      if( s1!=s ) { s=s1; _hIdxs[r].push_back(s); }
      _hIdxImg.val(c,r) = int(_hIdxs[r].size())-1;
    }
  }
  _vIdxs.resize(w); _vIdxImg.init(h,w);
  for( c=0; c<w; c++ ) {
    int s=0; _vIdxs[c].push_back(s);
    for( r=0; r<h; r++ ) {
      int s1 = _segIds.val(c,r);
      if( s1!=s ) { s=s1; _vIdxs[c].push_back(s); }
      _vIdxImg.val(c,r) = int(_vIdxs[c].size())-1;
    }
  }

  // initialize scoreBox() data structures
  int n=_segCnt+1; _sWts.init(n,1);
  _sDone.init(n,1); _sMap.init(n,1); _sIds.init(n,1);
  for( i=0; i<n; i++ ) _sDone.val(0,i)=-1; _sId=0;
}

void EdgeBoxGenerator::scoreBox( Box &box )
{
  int i, j, k, q, bh, bw, r0, c0, r1, c1, r0m, r1m, c0m, c1m;
  float *sWts=_sWts._x; int sId=_sId++;
  int *sDone=_sDone._x, *sMap=_sMap._x, *sIds=_sIds._x;
  // add edge count inside box
  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);
  bh=box.h=r1-box.r; bh/=2; bw=box.w=c1-box.c; bw/=2;
  float v = _segIImg.val(c0,r0) + _segIImg.val(c1+1,r1+1)
    - _segIImg.val(c1+1,r0) - _segIImg.val(c0,r1+1);
  // subtract middle quarter of edges
  r0m=r0+bh/2; r1m=r0m+bh; c0m=c0+bw/2; c1m=c0m+bw;
  v -= _magIImg.val(c0m, r0m) + _magIImg.val(c1m+1,r1m+1)
    - _magIImg.val(c1m+1,r0m) - _magIImg.val(c0m,r1m+1);
  // short circuit computation if impossible to score highly
  float norm = _scaleNorm[bw+bh]; box.s=v*norm;
  if( box.s<_minScore ) { box.s=0; return; }
  // find interesecting segments along four boundaries
  int cs, ce, rs, re, n=0;
  cs=_hIdxImg.val(c0,r0); ce=_hIdxImg.val(c1,r0); // top
  for( i=cs; i<=ce; i++ ) if( (j=_hIdxs[r0][i])>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<n; i++ ) {
    float w=sWts[i]; j=sIds[i];
    for( k=0; k<int(_segAffIdx[j].size()); k++ ) {
      q=_segAffIdx[j][k]; float wq=w*_segAff[j][k];
      if( wq<.05f ) continue; // short circuit for efficiency
      if( sDone[q]==sId ) {
        if( wq>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<n; i++ ) {
    k = sIds[i];
    if( _segC[k]>=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<w; c++ ) for( r=0; r<h; r++ )
    V.val(c+w*0,r)=V.val(c+w*1,r)=V.val(c+w*2,r)=1;
  for( c=0; c<w; c++ ) for( r=0; r<h; r++ ) {
    i=_segIds.val(c,r); if(i<=0) continue; e = E.val(c,r);
    o = (_sDone._x[i]==sId) ? _sWts._x[_sMap._x[i]] :
      (_segC[i]>=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<scNum; s++ ) {
    int a, r, c, bh, bw, kr, kc, bId=-1; float ar, sc;
    for( a=0; a<2*arRad+1; a++ ) {
      ar=pow(_arStep,float(a-arRad)); sc=minSize*pow(_scStep,float(s));
      bh=int(sc/ar); kr=max(2,int(bh*_rcStepRatio));
      bw=int(sc*ar); kc=max(2,int(bw*_rcStepRatio));
      for( c=0; c<w-bw+kc; c+=kc ) for( r=0; r<h-bh+kr; r+=kr ) {
        Box b; b.r=r; b.c=c; b.h=bh; b.w=bw; boxes.push_back(b);
      }
    }
  }

  // score all boxes, refine top candidates, perform nms
  int i, k=0, m = int(boxes.size());
  for( i=0; i<m; i++ ) {
    scoreBox(boxes[i]);
    if( !boxes[i].s ) continue; k++;
    refineBox(boxes[i]);
  }
  sort(boxes.rbegin(),boxes.rend(),boxesCompare);
  boxes.resize(k); boxesNms(boxes,_beta,_eta,_maxBoxes);
}

float boxesOverlap( Box &a, Box &b ) {
  float areai, areaj, areaij;
  int r0, r1, c0, c1, r1i, c1i, r1j, c1j;
  r1i=a.r+a.h; c1i=a.c+a.w; if( a.r>=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<Boxes> kept; kept.resize(nBin+1);
  int i=0, j, k, n=(int) boxes.size(), m=0, b, d=1;
  while( i<n && m<maxBoxes ) {
    b = boxes[i].w*boxes[i].h; bool keep=1;
    b = clamp(int(ceil(log(float(b))/lstep)),d,nBin-d);
    for( j=b-d; j<=b+d; j++ )
      for( k=0; k<kept[j].size(); k++ ) if( keep )
        keep = boxesOverlap( boxes[i], kept[j][k] ) <= thr;
    if(keep) { kept[b].push_back(boxes[i]); m++; } i++;
    if(keep && eta<1 && thr>.5) { thr*=eta; d=ceil(log(1/thr)/lstep); }
  }
  boxes.resize(m); i=0;
  for( j=0; j<nBin; j++ )
    for( k=0; k<kept[j].size(); k++ )
      boxes[i++]=kept[j][k];
  sort(boxes.rbegin(),boxes.rend(),boxesCompare);
}

////////////////////////////////////////////////////////////////////////////////

// Matlab entry point: bbs = mex( E, O, prm1, prm2, ... )
//matlab 中调用

//edgeBoxesMex(E,O,o.alpha,o.beta,o.eta,o.minScore,o.maxBoxes,...,o.edgeMinMag,o.edgeMergeThr,o.clusterMinMag,...,o.maxAspectRatio,o.minBoxArea,o.gamma,o.kappa);

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<n; i++) {
    bbs[ i + 0*n ] = (float) boxes[i].c+1;
    bbs[ i + 1*n ] = (float) boxes[i].r+1;
    bbs[ i + 2*n ] = (float) boxes[i].w;
    bbs[ i + 3*n ] = (float) boxes[i].h;
    bbs[ i + 4*n ] = boxes[i].s;
  }
}