75 files changed, 4657 insertions, 0 deletions
diff --git a/SD-VBS/common/toolbox/MultiNcut/MNcut.m b/SD-VBS/common/toolbox/MultiNcut/MNcut.m
new file mode 100755
index 0000000..5486080
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/MNcut.m
@@ -0,0 +1,93 @@
+function [NcutDiscretes,eigenVectors,eigenValues] = MNcut(I,nsegs);
+%
+%  [NcutDiscrete,eigenVectors,eigenValues] = MNcut(I,nsegs);
+%  
+%
+[nr,nc,nb] = size(I);
+max_image_size = max(nr,nc);
+% modified by song, 06/13/2005
+%   test parameters
+if (1) % original settings
+    if (max_image_size>120) & (max_image_size<=500),
+        % use 3 levels,
+        data.layers.number=3;
+        data.layers.dist=3;
+        data.layers.weight=[3000,4000,10000];
+        data.W.scales=[1,2,3];%[1,2,3];
+        data.W.radius=[2,3,7];%[2,3,7];
+    elseif (max_image_size >500),
+        % use 4 levels,
+        data.layers.number=4;
+        data.layers.dist=3;
+        data.layers.weight=[3000,4000,10000,20000];
+        data.W.scales=[1,2,3,3];
+        data.W.radius=[2,3,4,6];
+    elseif (max_image_size <=120)
+        data.layers.number=2;
+        data.layers.dist=3;
+        data.layers.weight=[3000,10000];
+        data.W.scales=[1,2];
+        data.W.radius=[2,6];
+    end
+else % test setting
+    if (max_image_size>200) & (max_image_size<=500),
+        % use 3 levels,
+        data.layers.number=3;
+        data.layers.dist=3;
+        data.layers.weight=[3000,4000,10000];
+        data.W.scales=[1,2,3];%[1,2,3];
+        data.W.radius=[2,3,7];%[2,3,7];
+    elseif (max_image_size >500),
+        % use 4 levels,
+        data.layers.number=4;
+        data.layers.dist=3;
+        data.layers.weight=[3000,4000,10000,20000];
+        data.W.scales=[1,2,3,3];
+        data.W.radius=[2,3,4,6];
+    elseif (max_image_size <=200)
+        data.layers.number=2;
+        data.layers.dist=3;
+        data.layers.weight=[3000,10000];
+        data.W.scales=[1,2];
+        data.W.radius=[2,4];
+    end
+end;
+data.W.edgeVariance=0.1; %0.1
+data.W.gridtype='square';
+data.W.sigmaI=0.12;%0.12
+data.W.sigmaX=1000;
+data.W.mode='mixed';    
+data.W.p=0;
+data.W.q=0;
+%eigensolver
+data.dataGraphCut.offset = 100;% 10; %valeur sur diagonale de W (mieux vaut 10 pour valeurs negatives de W)
+data.dataGraphCut.maxiterations=50;% voir
+data.dataGraphCut.eigsErrorTolerance=1e-2;%1e-6;
+data.dataGraphCut.valeurMin=1e-6;%1e-5;% utilise pour tronquer des valeurs et sparsifier des matrices
+data.dataGraphCut.verbose = 0;
+data.dataGraphCut.nbEigenValues=max(nsegs);
+disp('computeEdge');
+[multiWpp,ConstraintMat, Wind,data,emag,ephase]= computeMultiW (I,data);
+disp('Ncut');
+[eigenVectors,eigenValues]= eigSolve (multiWpp,ConstraintMat,data);
+%NcutDiscretes = zeros(nr,nc,length(nsegs));
+NcutDiscretes = zeros(nr,nc,(nsegs));
+for j=1:length(nsegs),
+        nseg = nsegs(j);
+        [nr,nc,nb] = size(eigenVectors(:,:,1:nseg));
+        [NcutDiscrete,evrotated] =discretisation(reshape(eigenVectors(:,:,1:nb),nr*nc,nb),nr,nc);
+         NcutDiscretes(:,:,j) = NcutDiscrete;
+end
diff --git a/SD-VBS/common/toolbox/MultiNcut/MNcutDemo.m b/SD-VBS/common/toolbox/MultiNcut/MNcutDemo.m
new file mode 100755
index 0000000..972a4eb
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/MNcutDemo.m
@@ -0,0 +1,34 @@
+% MNcutDemo.m
+% created by song, 06/13/2005
+%  an exmaple of how to use and display MNcut
+num_segs = [20];
+imageSize = 800;
+img_filename = '/u/ikkjin/Benchmark/stitch/data/test/capitol/img1.jpg';
+I=readimage(img_filename,imageSize);
+[SegLabel,eigenVectors,eigenValues]= MNcut(I,num_segs);
+for j=1:size(SegLabel,3),
+    [gx,gy] = gradient(SegLabel(:,:,j));
+    bw = (abs(gx)>0.1) + (abs(gy) > 0.1);
+    figure(1);clf; J1=showmask(double(I),bw); imagesc(J1);axis image; axis off;
+    set(gca, 'Position', [0 0 1 1]);
+    %      cm = sprintf('print -djpeg %s/file%.4d-%.2d.jpg',OutputDir,id,num_segs(j)); disp(cm);eval(cm);
+    %      figure(10);imagesc(SegLabel(:,:,j));axis image; axis off;
+    %      set(gca, 'Position', [0 0 1 1]);
+    %      cm = sprintf('print -djpeg %s/Seg%.4d-%.2d.jpg',OutputDir,id,num_segs(j));disp(cm);eval(cm);
+    % pause;
+end
+% fname = files(id).name;
+%cm = sprintf('save %s/SegLabl%.4d.mat I SegLabel fname',OutputDir,id); disp(cm);    eval(cm);
+%cm = sprintf('save %s/SegEig%.4d.mat eigenVectors eigenValues',OutputDir,id);disp(cm); eval(cm);
diff --git a/SD-VBS/common/toolbox/MultiNcut/README.tex b/SD-VBS/common/toolbox/MultiNcut/README.tex
new file mode 100755
index 0000000..5970fb2
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/README.tex
@@ -0,0 +1,9 @@
+1) You need to first compile the .c files,type
+>> compileAll('.');
+2) the top level function is called MNcut.m
diff --git a/SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.c b/SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.c
new file mode 100755
index 0000000..25def92
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.c
@@ -0,0 +1,405 @@
+/*================================================================
+a_times_b_cmplx.c = used by a couple of mex functions 
+provide Matrix vector multiplications,
+and solve triangular systems
+(sparse matrix and full vector)
+CSC_CmplxVecMult_CAB_double, CSR_CmplxVecMult_CAB_double, 
+CSCsymm_CmplxVecMult_CAB_double added by Mirko Visontai (10/24/2003)
+*=================================================================*/
+# include "math.h"
+/*C<-a*A*B+C*/
+void CSC_VecMult_CaABC_double(
+                 const int m,  const int k, const double alpha,
+                 const double *val, const int *indx, 
+                 const int *pntrb,
+                 const double *b,  
+                 double *c)
+{
+  int i,j,jb,je;
+    for (i=0;i!=k;i++){
+      jb = pntrb[i];
+      je = pntrb[i+1];
+      for (j=jb;j!=je;j++)
+        c[indx[j]] += alpha * b[i] * val[j];
+    }
+}
+/*C<-a*A'*B+C*/
+void CSR_VecMult_CaABC_double(
+                 const int k,  const int m, const double alpha,
+                 const double *val, const int *indx, 
+                 const int *pntrb,
+                 const double *b,  
+                 double *c)
+{
+  double t;
+  const double *pval;
+  int i,j,jb,je;
+    pval = val;
+    for (i=0;i!=m;i++) {
+      t = 0;
+      jb = pntrb[i];
+      je = pntrb[i+1];
+      for (j=jb;j!=je;j++)
+        t += alpha * b[indx[j]] * (*pval++);
+      c[i] += t;
+    }
+}
+/*C<-A*b */
+void CSC_VecMult_CAB_double(
+                 const int m,  const int k, /*nb_rows, nb_columns*/
+                 const double *val, const int *indx, 
+                 const int *pntrb,                  
+                 const double *b,  
+                 double *c                 
+                 )
+{
+  int i,j,jb,je;
+  double *pc=c;                                       
+    for (i=0;i!=m;i++) *pc++ = 0;                 
+    for (i=0;i!=k;i++){
+      jb = pntrb[i];
+      je = pntrb[i+1];
+      for (j=jb;j!=je;j++)
+        c[indx[j]] +=  b[i] * val[j];
+    }
+}
+/*C<-A*b (complex)*/
+void CSC_CmplxVecMult_CAB_double(
+                 const int m,  const int k, 
+                 const double *valr, const double *vali,
+                                 const int *indx, 
+                 const int *pntrb,                  
+                 const double *br, const double *bi,  
+                 double *cr, double *ci                 
+                 )
+{
+  int i,j,jb,je;
+  double *pcr=cr;                                       
+  double *pci=ci;                                       
+    for (i=0;i!=m;i++){
+                *pcr++ = 0.0;
+                *pci++ = 0.0;
+        }
+    for (i=0;i!=k;i++){
+      jb = pntrb[i];
+      je = pntrb[i+1];
+      for (j=jb;j!=je;j++){
+        cr[indx[j]] +=  (br[i] * valr[j]) - (bi[i] * vali[j]);
+        ci[indx[j]] +=  (br[i] * vali[j]) + (bi[i] * valr[j]);
+          }
+    }
+}
+/*C<-A'*b 
+ plus rapide que CSC_VecMult_CAB_double */
+void CSR_VecMult_CAB_double(
+                 const int k,  const int m,
+                 const double *val, const int *indx, 
+                 const int *pntrb,
+                 const double *b,  
+                 double *c
+                 )
+{
+  double t;
+  const double *pval;
+  double *pc=c;                                       
+  int i,j,jb,je;
+ 
+    for (i=0;i!=m;i++) *pc++ = 0;                 
+    pval = val;
+    for (i=0;i!=m;i++) {
+      t = 0;
+      jb = pntrb[i];
+      je = pntrb[i+1];
+      for (j=jb;j!=je;j++)
+        t +=  b[indx[j]] * (*pval++);
+      c[i] += t;
+    }
+}
+/*C<-A'*b (complex) 
+ plus rapide que CSC_VecMult_CAB_double */
+void CSR_CmplxVecMult_CAB_double(
+                 const int k,  const int m, 
+                 const double *valr, const double *vali,
+                                 const int *indx, 
+                 const int *pntrb,
+                 const double *br, const double *bi,  
+                 double *cr, double *ci
+                 )
+{
+  double tr, ti;
+  const double *pvalr;
+  const double *pvali;
+  double *pcr=cr;                                       
+  double *pci=ci;                                       
+  int i,j,jb,je;
+ 
+    for (i=0;i!=m;i++){
+                *pcr++ = 0.0;                 
+                *pci++ = 0.0;
+        }               
+    pvalr = valr;
+    pvali = vali;
+    for (i=0;i!=m;i++) {
+      tr = 0.0;
+          ti = 0.0;
+      jb = pntrb[i];
+      je = pntrb[i+1];
+      for (j=jb;j!=je;j++){
+        tr +=  (br[indx[j]] * (*pvalr)) - (bi[indx[j]] * (*pvali));
+        ti +=  (br[indx[j]] * (*pvali++)) + (bi[indx[j]] * (*pvalr++));
+          }
+      cr[i] += tr;
+          ci[i] += ti;
+    }
+}
+/* C<-A*b (A is symmetric) */
+void CSRsymm_VecMult_CAB_double(
+                 const int k,  const int m, 
+                 const double *val, const int *indx, 
+                 const int *pntrb,  
+                 const double *b,  
+                 double *c 
+                 )
+{
+  const double *pval;
+  double *pc=c;                                       
+  int i,j;
+  int jj;
+  int rpntrb, rpntre;
+  int index, nvals;
+  
+      
+    for (i=0;i!=m;i++) *pc++ = 0;                 
+    pval = val;
+    for (j=0;j!=k;j++){
+      rpntrb = pntrb[j];
+      rpntre = pntrb[j+1];
+      for (jj=rpntrb;jj!=rpntre;jj++) {
+        index = indx[jj];
+        if ( index == j ) {
+          c[j] +=  b[j] * (*pval++);
+          continue;
+        }
+                if ( index > j ) {
+                c[index] +=  b[j] * (*pval);
+                c[j] +=  b[index] * (*pval++);
+                }
+                else {
+                        pval++;
+                }
+      }
+    }
+}
+/* C<-A*b (A is symmetric and complex) */
+void CSRsymm_CmplxVecMult_CAB_double(
+                 const int k,  const int m, 
+                 const double *valr, const double *vali,
+                                 const int *indx, 
+                 const int *pntrb,  
+                 const double *br, const double *bi,  
+                 double *cr, double *ci
+                 )
+{
+  const double *pvalr, *pvali;
+  double *pcr=cr;                                       
+  double *pci=ci;                                       
+  int i,j;
+  int jj;
+  int rpntrb, rpntre;
+  int index, nvals;
+  
+      
+    for (i=0;i!=m;i++){
+                *pcr++ = 0.0;                 
+                *pci++ = 0.0;                 
+        }
+    pvalr = valr;
+    pvali = vali;
+    for (j=0;j!=k;j++){
+      rpntrb = pntrb[j];
+      rpntre = pntrb[j+1];
+      for (jj=rpntrb;jj!=rpntre;jj++) {
+        index = indx[jj];
+        if ( index == j ) {
+          cr[j] +=  (br[j] * (*pvalr)) - (bi[j] * (*pvali));
+          ci[j] +=  (br[j] * (*pvali++)) + (bi[j] * (*pvalr++));
+          continue;
+        }
+        if ( index > j ) {
+                cr[index] +=  (br[j] * (*pvalr)) - (bi[j] * (*pvali));
+                ci[index] +=  (br[j] * (*pvali)) + (bi[j] * (*pvalr));
+        
+                cr[j] +=  (br[index] * (*pvalr)) - (bi[index] * (*pvali));
+                ci[j] +=  (br[index] * (*pvali++)) + (bi[index] * (*pvalr++));
+                }
+                else {
+                        pvalr++;
+                        pvali++;
+                }
+        
+      }
+    }
+}
+/*C<-A\B; with Lower triangular A*/
+void CSC_VecTriangSlvLD_CAB_double( 
+                 const int m,  
+                 const double *val,
+                 const int *indx, const int *pntrb, 
+                 const double *b,  
+                 double *c)
+{
+  int i, j, jb, je;
+  double *pc=c;
+  double z; 
+    for (i=0;i!=m;i++){
+      *pc = b[i];
+      pc++;
+    }                                     
+    pc=c;
+    for (i=0;i!=m;i++) {
+      jb = pntrb[i];
+      je = pntrb[i+1];
+      z =  pc[i] / val[jb];
+      pc[i] = z;
+      for (j=jb+1;j<je;j++) {
+        c[indx[j]] -= z*val[j];
+      }
+    }
+}
+/*C<-A\B; with Upper triangular A*/
+void CSC_VecTriangSlvUD_CAB_double(
+                 const int m,  
+                   const double *val,
+                 const int *indx, const int *pntrb,
+                 const double *b,  
+                 double *c)
+{
+  int i, j, jb, je, index;
+  double *pc=c;
+  double z; 
+    for (i=0;i!=m;i++){
+      *pc = b[i];
+      pc++;
+    }                                     
+    pc=c;
+    for (i=m-1;i!=-1;i--) {
+      jb = pntrb[i];
+      je = pntrb[i+1]-1;
+      z = pc[i] /val[je];
+      pc[i] = z;
+      for (j=jb;j<je;j++) {
+          c[indx[j]] -= z * val[j];
+      }
+    }
+}
+/*C<-A'\B; where A is upper (little slower than CSC)*/
+void CSR_VecTriangSlvLD_CAB_double(
+                 const int m,   
+                 const double *val, 
+                 const int *indx, const int *pntrb,
+                 const double *b,   
+                 double *c)
+{
+  int i, j, jb, je, index;
+  double *pc=c;
+  double z; 
+  double valtmp;
+    pc=c;
+    for (i=0;i!=m;i++) {
+      z = 0;
+      jb = pntrb[i];
+      je = pntrb[i+1];
+      for (j=jb;j<je;j++) {
+        index = indx[j];
+        if ( index == i ) {
+          valtmp = val[j];
+        } else {
+          z += c[index] * val[j];
+        }
+      }
+      pc[i] = (b[i] - z) / valtmp;
+    }
+}
+/*C<-A'\B; where A is lower  (little slower than CSC)*/
+void CSR_VecTriangSlvUD_CAB_double(
+                 const int m,   
+                   const double *val, 
+                 const int *indx, const int *pntrb, 
+                 const double *b,   
+                 double *c)
+{
+  int i, j, jb, je, index;
+  double *pc=c;
+  double valtmp;
+  double z; 
+    pc=c;
+    for (i=m-1;i!=-1; i--) {
+      z = 0;
+      jb = pntrb[i];
+      je =  pntrb[i+1];
+      for (j=jb+1; j<je; j++) {
+        z += c[indx[j]] * val[j];
+      }
+      pc[i] = (b[i] - z) / val[jb];
+    }
+}
+/*C<-A*B, where A is (m,k) and B is (k,n)*/
+void CSC_MatMult_CAB_double(
+                 const int m, const int n, const int k, 
+                 const double *val, const int *indx, 
+                 const int *pntrb,
+                 const double *b, const int ldb, 
+                 double *c, const int ldc)
+{
+  int i,j,jb,je;
+  double *pc=c;                                       
+  int l;                                               
+  for (l=0;l!=n;l++)                           
+    for (i=0;i!=m;i++) *pc++ = 0;                 
+  for (l=0;l!=n;l++) {                                 
+    for (i=0;i!=k;i++){
+      jb = pntrb[i];
+      je = pntrb[i+1];
+      for (j=jb;j!=je;j++)
+        c[indx[j]] +=  b[i] * val[j];
+    }
+    /*c += ldc; b += ldb;  */                              
+    c += m; b += m;                                
+  }                                                    
+}
diff --git a/SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.dll b/SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.dll
new file mode 100755
index 0000000..5656f92
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.dll
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.mexa64 b/SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.mexa64
new file mode 100755
index 0000000..1ab8e8c
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.mexa64
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.mexglx b/SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.mexglx
new file mode 100755
index 0000000..4138128
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.mexglx
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.mexmac b/SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.mexmac
new file mode 100755
index 0000000..65232e8
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.mexmac
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/affinityic.c b/SD-VBS/common/toolbox/MultiNcut/affinityic.c
new file mode 100755
index 0000000..d78761b
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/affinityic.c
@@ -0,0 +1,187 @@
+/*================================================================
+* function w = affinityic(emag,ephase,pi,pj,sigma)
+* Input:
+*   emag = edge strength at each pixel
+*   ephase = edge phase at each pixel
+*   [pi,pj] = index pair representation for MALTAB sparse matrices
+*   sigma = sigma for IC energy
+* Output:
+*   w = affinity with IC at [pi,pj]
+*
+% test sequence
+f = synimg(10);
+[i,j] = cimgnbmap(size(f),2);
+[ex,ey,egx,egy] = quadedgep(f);
+a = affinityic(ex,ey,egx,egy,i,j)
+show_dist_w(f,a);
+* Stella X. Yu, Nov 19, 2001.
+*=================================================================*/
+# include "mex.h"
+# include "math.h"
+void mexFunction(
+    int nargout,
+    mxArray *out[],
+    int nargin,
+    const mxArray *in[]
+)
+{
+    /* declare variables */
+    int nr, nc, np, total;
+    int i, j, k, ix, iy, jx, jy, ii, jj, iip1, jjp1, iip2, jjp2, step;
+    double sigma, di, dj, a, z, maxori, phase1, phase2, slope;
+        int *ir, *jc;
+    unsigned int *pi,*pj;
+    /* unsigned long *pi, *pj;
+ */
+        double *emag, *ephase, *w;
+    
+    /* check argument */
+    if (nargin<4) {
+        mexErrMsgTxt("Four input arguments required");
+    }
+    if (nargout>1) {
+        mexErrMsgTxt("Too many output arguments");
+    }
+    /* get edgel information */
+        nr = mxGetM(in[0]);
+        nc = mxGetN(in[0]);
+        if ( nr*nc ==0 || nr != mxGetM(in[1]) || nc != mxGetN(in[1]) ) {
+            mexErrMsgTxt("Edge magnitude and phase shall be of the same image size");
+        }
+    emag = mxGetPr(in[0]);
+    ephase = mxGetPr(in[1]);
+    np = nr * nc;
+    
+    /* get new index pair */
+    if (!mxIsUint32(in[2]) | !mxIsUint32(in[3])) {
+        mexErrMsgTxt("Index pair shall be of type UINT32");
+    }
+    if (mxGetM(in[3]) * mxGetN(in[3]) != np + 1) {
+        mexErrMsgTxt("Wrong index representation");
+    }
+    pi = mxGetData(in[2]);
+    pj = mxGetData(in[3]);    
+    /* create output */
+    out[0] = mxCreateSparse(np,np,pj[np],mxREAL);
+    if (out[0]==NULL) {
+            mexErrMsgTxt("Not enough memory for the output matrix");
+        }
+        w = mxGetPr(out[0]);
+        ir = mxGetIr(out[0]);
+        jc = mxGetJc(out[0]);
+        
+    /* find my sigma */
+        if (nargin<5) {
+            sigma = 0;
+        for (k=0; k<np; k++) { 
+            if (emag[k]>sigma) { sigma = emag[k]; }
+        }
+        sigma = sigma / 6;
+        /* printf("sigma = %6.5f",sigma); */
+        } else {
+            sigma = mxGetScalar(in[4]);
+        }
+        a = 0.5 / (sigma * sigma);
+        
+    /* computation */ 
+    total = 0;
+    for (j=0; j<np; j++) {            
+        jc[j] = total;
+        jx = j / nr; /* col */
+        jy = j % nr; /* row */
+        
+        for (k=pj[j]; k<pj[j+1]; k++) {
+        
+            i = pi[k];
+          
+            if (i==j) {
+                maxori = 1;
+            
+            } else {
+        
+                ix = i / nr; 
+                iy = i % nr;
+                /* scan */            
+                di = (double) (iy - jy);
+                dj = (double) (ix - jx);
+            
+                maxori = 0.;
+                    phase1 = ephase[j];
+                       
+                /* sample in i direction */
+                if (abs(di) >= abs(dj)) {  
+                    slope = dj / di;
+                    step = (iy>=jy) ? 1 : -1;
+                
+                    iip1 = jy;
+                    jjp1 = jx;
+        
+                       
+                        for (ii=0;ii<abs(di);ii++){
+                            iip2 = iip1 + step;
+                            jjp2 = (int)(0.5 + slope*(iip2-jy) + jx);
+                  
+                            phase2 = ephase[iip2+jjp2*nr];
+               
+                            if (phase1 != phase2) {
+                                z = (emag[iip1+jjp1*nr] + emag[iip2+jjp2*nr]);
+                                if (z > maxori){
+                                    maxori = z;
+                                }
+                            } 
+                     
+                            iip1 = iip2;
+                            jjp1 = jjp2;
+                            phase1 = phase2;
+                        }
+                    
+                    /* sample in j direction */    
+                } else { 
+                        slope = di / dj;
+                        step =  (ix>=jx) ? 1: -1;
+                    jjp1 = jx;
+                        iip1 = jy;                 
+            
+         
+                        for (jj=0;jj<abs(dj);jj++){
+                            jjp2 = jjp1 + step;
+                            iip2 = (int)(0.5+ slope*(jjp2-jx) + jy);
+                  
+                            phase2 = ephase[iip2+jjp2*nr];
+             
+                            if (phase1 != phase2){
+                                z = (emag[iip1+jjp1*nr] + emag[iip2+jjp2*nr]);
+                                if (z > maxori){ 
+                                    maxori = z; 
+                                }
+                                
+                            }
+          
+                            iip1 = iip2;
+                            jjp1 = jjp2;
+                            phase1 = phase2;
+                        }
+                }            
+            
+                maxori = 0.5 * maxori;
+                maxori = exp(-maxori * maxori * a);
+            }       
+                    ir[total] = i;
+                    w[total] = maxori;
+                    total = total + 1;
+                        
+                } /* i */
+    } /* j */
+        
+    jc[np] = total;
+}  
diff --git a/SD-VBS/common/toolbox/MultiNcut/affinityic.dll b/SD-VBS/common/toolbox/MultiNcut/affinityic.dll
new file mode 100755
index 0000000..67e4d64
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/affinityic.dll
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/affinityic.mexa64 b/SD-VBS/common/toolbox/MultiNcut/affinityic.mexa64
new file mode 100755
index 0000000..2648387
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/affinityic.mexa64
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/affinityic.mexglx b/SD-VBS/common/toolbox/MultiNcut/affinityic.mexglx
new file mode 100755
index 0000000..c296845
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/affinityic.mexglx
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/batch_MNcut.m b/SD-VBS/common/toolbox/MultiNcut/batch_MNcut.m
new file mode 100755
index 0000000..d4dcb3c
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/batch_MNcut.m
@@ -0,0 +1,48 @@
+data_dir = '/data/insecure/qihuizhu/baseball/Gray/train/';
+save_dir = '/home/songgang/project/MultiNcut/batch_result_MNcut';
+num_segs = [20];
+imageSize = 200;
+filelist = dir(fullfile(data_dir, '*.tif'));
+nb_file = max(size(filelist));
+tic;
+for ii = 1:nb_file
+    fprintf(2, 'Segmenting image: %s ...\n', filelist(ii).name);
+    img_filename = fullfile(data_dir, filelist(ii).name);
+    I=readimage(img_filename,imageSize);
+    
+    
+    [SegLabel,eigenVectors,eigenValues]= MNcut(I,num_segs);
+    for j=1:size(SegLabel,3),
+        [gx,gy] = gradient(SegLabel(:,:,j));
+        bw = (abs(gx)>0.1) + (abs(gy) > 0.1);
+        figure(1);clf; J1=showmask(double(I),bw); imagesc(J1);axis image; axis off;
+        set(gca, 'Position', [0 0 1 1]);
+        set(gca, 'Position', [0 0 1 1]);
+        [PATHSTR,NAME,EXT,VERSN] = fileparts(filelist(ii).name);
+        print('-f1', '-djpeg90', fullfile(save_dir, sprintf('%s%s-%d.jpg', NAME,'-out', num_segs(j))));
+        %      cm = sprintf('print -djpeg %s/file%.4d-%.2d.jpg',OutputDir,id,num_segs(j)); disp(cm);eval(cm);
+        %      figure(10);imagesc(SegLabel(:,:,j));axis image; axis off;
+        %      set(gca, 'Position', [0 0 1 1]);
+        %      cm = sprintf('print -djpeg %s/Seg%.4d-%.2d.jpg',OutputDir,id,num_segs(j));disp(cm);eval(cm);
+%        keyboard;
+    end
+end;
+toc;
+fprintf(2, ' %d files done\n', nb_file);
diff --git a/SD-VBS/common/toolbox/MultiNcut/cimgnbmap.c b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap.c
new file mode 100755
index 0000000..44af715
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap.c
@@ -0,0 +1,198 @@
+/*================================================================
+* function [i,j] = cimgnbmap([nr,nc], nb_r, sample_rate)
+*   computes the neighbourhood index matrix of an image,
+*   with each neighbourhood sampled.
+* Input:
+*   [nr,nc] = image size
+*   nb_r = neighbourhood radius, could be [r_i,r_j] for i,j
+*   sample_rate = sampling rate, default = 1
+* Output:
+*   [i,j] = each is a column vector, give indices of neighbour pairs
+*     UINT32 type
+*       i is of total length of valid elements, 0 for first row
+*       j is of length nr * nc + 1
+*
+* See also: imgnbmap.c, id2cind.m
+*
+* Examples: 
+*   [i,j] = imgnbmap(10, 20); % [10,10] are assumed
+*
+* Stella X. Yu, Nov 12, 2001.
+% test sequence:
+nr = 15;
+nc = 15;
+nbr = 1;
+[i,j] = cimgnbmap([nr,nc], nbr);
+mask = csparse(i,j,ones(length(i),1),nr*nc);
+show_dist_w(rand(nr,nc),mask)
+*=================================================================*/
+# include "mex.h"
+# include "math.h"
+void mexFunction(
+    int nargout,
+    mxArray *out[],
+    int nargin,
+    const mxArray *in[]
+)
+{
+    /* declare variables */
+    int nr, nc, np, nb, total;
+        double *dim, sample_rate;
+    int r_i, r_j, a1, a2, b1, b2, self, neighbor;
+    int i, j, k, s, t, nsamp, th_rand, no_sample;
+    /* unsigned long *p, *qi, *qj; */
+    unsigned int *p, *qi, *qj;
+    
+    /* check argument */
+    if (nargin < 2) {
+        mexErrMsgTxt("Two input arguments required");
+    }
+    if (nargout> 2) {
+        mexErrMsgTxt("Too many output arguments.");
+    }
+    /* get image size */
+    i = mxGetM(in[0]);
+    j = mxGetN(in[0]);
+    dim = mxGetData(in[0]);
+    nr = (int)dim[0];
+    if (j>1 || i>1) {
+        nc = (int)dim[1];
+    } else {
+        nc = nr;
+    }
+    np = nr * nc;
+    
+    /* get neighbourhood size */
+    i = mxGetM(in[1]);
+    j = mxGetN(in[1]);
+    dim = mxGetData(in[1]);
+    r_i = (int)dim[0];
+    if (j>1 || i>1) {
+        r_j = (int)dim[1];              
+    } else {
+        r_j = r_i;
+    }
+    if (r_i<0) { r_i = 0; }
+        if (r_j<0) { r_j = 0; }
+        /* get sample rate */
+        if (nargin==3) {                
+                sample_rate = (mxGetM(in[2])==0) ? 1: mxGetScalar(in[2]);
+    } else {
+                sample_rate = 1;
+    }
+        /* prepare for random number generator */
+        if (sample_rate<1) {
+        srand( (unsigned)time( NULL ) );
+        th_rand = (int)ceil((double)RAND_MAX * sample_rate);
+        no_sample = 0;
+    } else {
+                sample_rate = 1;
+        th_rand = RAND_MAX;
+        no_sample = 1;
+    }
+    
+        /* figure out neighbourhood size */
+    nb = (r_i + r_i + 1) * (r_j + r_j + 1); 
+    if (nb>np) {
+        nb = np;
+    }
+    nb = (int)ceil((double)nb * sample_rate);    
+/*printf("nb=%d\n",nb);*/
+        /* intermediate data structure */
+    /* p = mxCalloc(np * (nb+1), sizeof(unsigned long));
+ */
+        p = mxCalloc(np * (nb+1), sizeof(unsigned int));
+        if (p==NULL) {
+        mexErrMsgTxt("Not enough space for my computation.");
+        }
+        
+    /* computation */    
+        total = 0;
+        for (j=0; j<nc; j++) {
+    /*printf("j=%d\n",j);*/
+    for (i=0; i<nr; i++) {
+                self = i + j * nr;
+                /* put self in, otherwise the index is not ordered */
+                p[self] = p[self] + 1;
+                p[self+p[self]*np] = self;
+   
+        /* j range */
+                b1 = j;
+        b2 = j + r_j;
+        if (b2>=nc) { b2 = nc-1; }                
+ 
+                /* i range */
+        a1 = i - r_i;
+                if (a1<0) { a1 = 0; }
+        a2 = i + r_i;
+        if (a2>=nr) { a2 = nr-1; }
+     
+                /* number of more samples needed */
+                nsamp = nb - p[self];
+        /*if (nsamp<0)
+        {printf("nsamp=%d\n",nsamp);}*/
+                k = 0;          
+                t = b1;
+                s = i + 1;
+                if (s>a2) {
+                        s = a1;
+                        t = t + 1;
+                }
+                   
+                
+                while (k<nsamp && t<=b2) {
+                
+                        if (no_sample || (rand()<th_rand)) {
+                                k = k + 1;
+                                neighbor = s + t * nr;
+                                p[self] = p[self] + 1;                                  
+                                p[self+p[self]*np] = neighbor;
+                                p[neighbor] = p[neighbor] + 1;
+                                p[neighbor+p[neighbor]*np] = self;
+                                        }
+                                            
+                        s = s + 1;
+                        if (s>a2) {
+                s = a1;
+                                t = t + 1;
+                        }
+                } /* k */
+                total = total + p[self];
+        } /* i */
+           
+    } /* j */
+   
+    /* i, j */
+    
+    out[0] = mxCreateNumericMatrix(total, 1, mxUINT32_CLASS, mxREAL);
+        out[1] = mxCreateNumericMatrix(np+1,  1, mxUINT32_CLASS, mxREAL);
+        qi = mxGetData(out[0]);
+        qj = mxGetData(out[1]);
+          
+        if (out[0]==NULL || out[1]==NULL) {
+            mexErrMsgTxt("Not enough space for the output matrix.");
+        }
+        total = 0;
+    for (j=0; j<np; j++) {
+                qj[j] = total;
+                s = j + np;
+                for (t=0; t<p[j]; t++) {
+                    qi[total] = p[s];
+                        total = total + 1;
+                        s = s + np;
+                }
+    }
+        qj[np] = total;
+        mxFree(p);
+        
+}  
diff --git a/SD-VBS/common/toolbox/MultiNcut/cimgnbmap.dll b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap.dll
new file mode 100755
index 0000000..368d76c
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap.dll
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/cimgnbmap.mexa64 b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap.mexa64
new file mode 100755
index 0000000..6564652
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap.mexa64
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/cimgnbmap.mexglx b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap.mexglx
new file mode 100755
index 0000000..655849a
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap.mexglx
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/cimgnbmap_cross.c b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap_cross.c
new file mode 100755
index 0000000..7b1434c
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap_cross.c
@@ -0,0 +1,197 @@
+/*================================================================
+* function [i,j] = cimgnbmap_cross([nr,nc], nb_r, sample_rate)
+*   computes the neighbourhood index matrix of an image,
+*   with each neighbourhood sampled.
+* Input:
+*   [nr,nc] = image size
+*   nb_r = neighbourhood radius, could be [r_i,r_j] for i,j
+*   sample_rate = sampling rate, default = 1
+* Output:
+*   [i,j] = each is a column vector, give indices of neighbour pairs
+*     UINT32 type
+*       i is of total length of valid elements, 0 for first row
+*       j is of length nr * nc + 1
+*
+* See also: imgnbmap.c, id2cind.m
+*
+* Examples: 
+*   [i,j] = imgnbmap(10, 20); % [10,10] are assumed
+*
+* Stella X. Yu, Nov 12, 2001.
+% test sequence:
+nr = 15;
+nc = 15;
+nbr = 1;
+[i,j] = cimgnbmap([nr,nc], nbr);
+mask = csparse(i,j,ones(length(i),1),nr*nc);
+show_dist_w(rand(nr,nc),mask)
+*=================================================================*/
+# include "mex.h"
+# include "math.h"
+void mexFunction(
+    int nargout,
+    mxArray *out[],
+    int nargin,
+    const mxArray *in[]
+)
+{
+    /* declare variables */
+    int nr, nc, np, nb, total;
+        double *dim, sample_rate;
+    int r_i, r_j, a1, a2, b1, b2, self, neighbor;
+    int i, j, k, s, t, nsamp, th_rand, no_sample;
+    /*     unsigned long *p, *qi, *qj;
+ */
+    unsigned int *p, *qi, *qj;
+    
+    /* check argument */
+    if (nargin < 2) {
+        mexErrMsgTxt("Two input arguments required");
+    }
+    if (nargout> 2) {
+        mexErrMsgTxt("Too many output arguments.");
+    }
+    /* get image size */
+    i = mxGetM(in[0]);
+    j = mxGetN(in[0]);
+    dim = mxGetData(in[0]);
+    nr = (int)dim[0];
+    if (j>1 || i>1) {
+        nc = (int)dim[1];
+    } else {
+        nc = nr;
+    }
+    np = nr * nc;
+    
+    /* get neighbourhood size */
+    i = mxGetM(in[1]);
+    j = mxGetN(in[1]);
+    dim = mxGetData(in[1]);
+    r_i = (int)dim[0];
+    if (j>1 || i>1) {
+        r_j = (int)dim[1];              
+    } else {
+        r_j = r_i;
+    }
+    if (r_i<0) { r_i = 0; }
+        if (r_j<0) { r_j = 0; }
+        /* get sample rate */
+        if (nargin==3) {                
+                sample_rate = (mxGetM(in[2])==0) ? 1: mxGetScalar(in[2]);
+    } else {
+                sample_rate = 1;
+    }
+        /* prepare for random number generator */
+        if (sample_rate<1) {
+        srand( (unsigned)time( NULL ) );
+        th_rand = (int)ceil((double)RAND_MAX * sample_rate);
+        no_sample = 0;
+    } else {
+                sample_rate = 1;
+        th_rand = RAND_MAX;
+        no_sample = 1;
+    }
+    
+        /* figure out neighbourhood size */
+    nb = (r_i + r_i) * (r_j + r_j)+1; 
+    if (nb>np) {
+        nb = np;
+    }
+    nb = (int)ceil((double)nb * sample_rate);    
+/*printf("nb=%d\n",nb);*/
+        /* intermediate data structure */
+    /*  p = mxCalloc(np * (nb+1), sizeof(unsigned long));
+ */
+        p = mxCalloc(np * (nb+1), sizeof(unsigned int));
+        if (p==NULL) {
+        mexErrMsgTxt("Not enough space for my computation.");
+        }
+        
+    /* computation */    
+        total = 0;
+        for (j=0; j<nc; j++) {
+    /*printf("j=%d\n",j);*/
+    for (i=0; i<nr; i++) {
+                self = i + j * nr;
+                /* put self in, otherwise the index is not ordered */
+                p[self] = p[self] + 1;
+                p[self+p[self]*np] = self;
+   
+        /* j range */
+                b1 = j;
+        b2 = j + r_j;
+        if (b2>=nc) { b2 = nc-1; }                
+ 
+                /* i range */
+        /*a1 = i - r_i;
+                if (a1<0) { a1 = 0; }*/
+        a2 = i + r_i;
+        if (a2>=nr) { a2 = nr-1; }
+     
+                /* number of more samples needed */
+                nsamp = nb - p[self];
+        /*if (nsamp<0)
+        {printf("nsamp=%d\n",nsamp);}*/
+                k = 0;          
+                t = b1;
+                s = i + 1;
+                if (s>a2) {
+                        s = i;
+                        t = t + 1;
+                }
+                   
+                
+                while (k<nsamp && t<=b2) {
+                
+                        if (no_sample || (rand()<th_rand)) {
+                                k = k + 1;
+                                neighbor = s + t * nr;
+                                p[self] = p[self] + 1;                                  
+                                p[self+p[self]*np] = neighbor;
+                                p[neighbor] = p[neighbor] + 1;
+                                p[neighbor+p[neighbor]*np] = self;
+                                        }
+                                            
+            if (s!=i){
+                            s = s + 1;
+                         if (s>a2) {
+                    s = i;
+                                 t = t + 1;
+                                      }
+                                  }
+                        else {t=t+1;}
+                } /* k */
+                total = total + p[self];
+        } /* i */
+           
+    } /* j */
+   
+    /* i, j */
+    
+    out[0] = mxCreateNumericMatrix(total, 1, mxUINT32_CLASS, mxREAL);
+        out[1] = mxCreateNumericMatrix(np+1,  1, mxUINT32_CLASS, mxREAL);
+        qi = mxGetData(out[0]);
+        qj = mxGetData(out[1]);
+          
+        if (out[0]==NULL || out[1]==NULL) {
+            mexErrMsgTxt("Not enough space for the output matrix.");
+        }
+        total = 0;
+    for (j=0; j<np; j++) {
+                qj[j] = total;
+                s = j + np;
+                for (t=0; t<p[j]; t++) {
+                    qi[total] = p[s];
+                        total = total + 1;
+                        s = s + np;
+                }
+    }
+        qj[np] = total;
+        mxFree(p);
+        
+}  
diff --git a/SD-VBS/common/toolbox/MultiNcut/cimgnbmap_cross.dll b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap_cross.dll
new file mode 100755
index 0000000..d66e578
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap_cross.dll
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/cimgnbmap_cross.mexa64 b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap_cross.mexa64
new file mode 100755
index 0000000..905aae5
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap_cross.mexa64
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/cimgnbmap_cross.mexglx b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap_cross.mexglx
new file mode 100755
index 0000000..a7fc50d
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap_cross.mexglx
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/cimgnbmap_star.c b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap_star.c
new file mode 100755
index 0000000..e8a37d6
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap_star.c
@@ -0,0 +1,294 @@
+/*================================================================
+* function [i,j] = cimgnbmap_star([nr,nc], nb_r, sample_rate)
+*   computes the neighbourhood index matrix of an image,
+*   with each neighbourhood sampled.
+* Input:
+*   [nr,nc] = image size
+*   nb_r = neighbourhood radius, could be [r_i,r_j] for i,j
+*   sample_rate = sampling rate, default = 1
+* Output:
+*   [i,j] = each is a column vector, give indices of neighbour pairs
+*     UINT32 type
+*       i is of total length of valid elements, 0 for first row
+*       j is of length nr * nc + 1
+*
+* See also: imgnbmap.c, id2cind.m
+*
+* Examples: 
+*   [i,j] = imgnbmap(10, 20); % [10,10] are assumed
+*
+* Stella X. Yu, Nov 12, 2001.
+% test sequence:
+nr = 15;
+nc = 15;
+nbr = 1;
+[i,j] = cimgnbmap([nr,nc], nbr);
+mask = csparse(i,j,ones(length(i),1),nr*nc);
+show_dist_w(rand(nr,nc),mask)
+*=================================================================*/
+# include "mex.h"
+# include "math.h"
+void mexFunction(
+    int nargout,
+    mxArray *out[],
+    int nargin,
+    const mxArray *in[]
+)
+{
+    /* declare variables */
+    int nr, nc, np, nb, total;
+        double *dim, sample_rate;
+    int r_i, r_j, a1, a2, b1, b2, self, neighbor;
+    int i, j, k, s, t, nsamp, th_rand, no_sample;
+    /*     unsigned long *p, *qi, *qj; */
+    unsigned int *p, *qi, *qj;
+    
+    /* check argument */
+    if (nargin < 2) {
+        mexErrMsgTxt("Two input arguments required");
+    }
+    if (nargout> 2) {
+        mexErrMsgTxt("Too many output arguments.");
+    }
+    /* get image size */
+    i = mxGetM(in[0]);
+    j = mxGetN(in[0]);
+    dim = mxGetData(in[0]);
+    nr = (int)dim[0];
+    if (j>1 || i>1) {
+        nc = (int)dim[1];
+    } else {
+        nc = nr;
+    }
+    np = nr * nc;
+    
+    /* get neighbourhood size */
+    i = mxGetM(in[1]);
+    j = mxGetN(in[1]);
+    dim = mxGetData(in[1]);
+    r_i = (int)dim[0];
+    if (j>1 || i>1) {
+        r_j = (int)dim[1];              
+    } else {
+        r_j = r_i;
+    }
+    if (r_i<0) { r_i = 0; }
+        if (r_j<0) { r_j = 0; }
+        /* get sample rate */
+        if (nargin==3) {                
+                sample_rate = (mxGetM(in[2])==0) ? 1: mxGetScalar(in[2]);
+    } else {
+                sample_rate = 1;
+    }
+        /* prepare for random number generator */
+        if (sample_rate<1) {
+        srand( (unsigned)time( NULL ) );
+        th_rand = (int)ceil((double)RAND_MAX * sample_rate);
+        no_sample = 0;
+    } else {
+                sample_rate = 1;
+        th_rand = RAND_MAX;
+        no_sample = 1;
+    }
+    
+        /* figure out neighbourhood size */
+    nb = (4*r_i) + (4*r_j)+1; 
+    if (nb>np) {
+        nb = np;
+    }
+    nb = (int)ceil((double)nb * sample_rate);    
+/*printf("nb=%d\n",nb);*/
+        /* intermediate data structure */
+    /*  p = mxCalloc(np * (nb+1), sizeof(unsigned long));*/
+        p = mxCalloc(np * (nb+1), sizeof(unsigned int));
+        if (p==NULL) {
+        mexErrMsgTxt("Not enough space for my computation.");
+        }
+        
+    /* computation */    
+        total = 0;
+        for (j=0; j<nc; j++) {
+    /*printf("j=%d\n",j);*/
+    for (i=0; i<nr; i++) {
+                self = i + j * nr;
+                /* put self in, otherwise the index is not ordered */
+                p[self] = p[self] + 1;
+                p[self+p[self]*np] = self;
+   
+        /* j range */
+                b1 = j;
+        b2 = j + r_j;
+        if (b2>=nc) { b2 = nc-1; }                
+ 
+                /* i range */
+        /*a1 = i - r_i;
+                if (a1<0) { a1 = 0; }*/
+        a2 = i + r_i;
+        if (a2>=nr) { a2 = nr-1; }
+     
+                /* number of more samples needed */
+                nsamp = nb - p[self];
+        /*if (nsamp<0)
+        {printf("nsamp=%d\n",nsamp);}*/
+                k = 0;          
+                t = b1;
+                s = i + 1;
+                if (s>a2) {
+                        s = i;
+                        t = t + 1;
+                }
+                   
+                
+                while (k<nsamp && t<=b2) {
+                
+                        if (no_sample || (rand()<th_rand)) {
+                                k = k + 1;
+                                neighbor = s + t * nr;
+                                p[self] = p[self] + 1;                                  
+                                p[self+p[self]*np] = neighbor;
+                                p[neighbor] = p[neighbor] + 1;
+                                p[neighbor+p[neighbor]*np] = self;
+                                        }
+                                            
+            if (t==b1){
+                                    s = s + 1;
+                                if (s>a2) {
+                                    t = t + 1;
+                                    if (i+j-t>=0)
+                                                {s = i+j-t;}
+                                                else
+                                                    {s=i;}
+                                                }
+                                  }
+                        else { 
+                                if (s==i+j-t) {s=i;}
+                                else{ if (s==i && s+t-j<nr) {s=s+t-j;}
+                                       else {
+                                                t = t + 1;
+                                    if (i+j-t>=0)
+                                                {s = i+j-t;}
+                                                else
+                                                    {s=i;}
+                                             }
+                                 
+                                 }
+                            }
+                } /* k */
+                total = total + p[self];
+        } /* i */
+           
+    } /* j */
+   
+    /* i, j */
+    
+    out[0] = mxCreateNumericMatrix(total, 1, mxUINT32_CLASS, mxREAL);
+        out[1] = mxCreateNumericMatrix(np+1,  1, mxUINT32_CLASS, mxREAL);
+        qi = mxGetData(out[0]);
+        qj = mxGetData(out[1]);
+          
+        if (out[0]==NULL || out[1]==NULL) {
+            mexErrMsgTxt("Not enough space for the output matrix.");
+        }
+        total = 0;
+    for (j=0; j<np; j++) {
+                qj[j] = total;
+                s = j + np;
+                for (t=0; t<p[j]; t++) {
+                    qi[total] = p[s];
+                        total = total + 1;
+                        s = s + np;
+                }
+    }
+        qj[np] = total;
+        mxFree(p);
+        
+}  
diff --git a/SD-VBS/common/toolbox/MultiNcut/cimgnbmap_star.dll b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap_star.dll
new file mode 100755
index 0000000..a7cddfa
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap_star.dll
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/cimgnbmap_star.mexa64 b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap_star.mexa64
new file mode 100755
index 0000000..d687d91
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap_star.mexa64
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/cimgnbmap_star.mexglx b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap_star.mexglx
new file mode 100755
index 0000000..25da17f
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap_star.mexglx
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/compileAll.m b/SD-VBS/common/toolbox/MultiNcut/compileAll.m
new file mode 100755
index 0000000..16fe77b
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/compileAll.m
@@ -0,0 +1,10 @@
+function compileAll(Cdir);
+files = dir([Cdir,'/*.c']);
+for j=1:length(files),
+        
+        cm = sprintf('mex %s',files(j).name);
+disp(cm);
+eval(cm);
+end
diff --git a/SD-VBS/common/toolbox/MultiNcut/computeMultiW.m b/SD-VBS/common/toolbox/MultiNcut/computeMultiW.m
new file mode 100755
index 0000000..b2cb7ea
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/computeMultiW.m
@@ -0,0 +1,245 @@
+%inputs : image, number of layers, distance defining the subgrid, the edge filter scales for each layer, radius for each layer,
+%edge variance for filter, shape of the neighbourhood layout ('square', 'star', 'cross'), sigma for intensity affinity,
+% sigma for distance influence in affinity, weight coefficients for Wpps in the multiscale matrix.
+%output : multiscale affinity matrix , extern constraint matrix, affinity matrices of each layer seperatly.
+function [multiWpp,constraintMat, Wind,data,emag,ephase]= computeMultiW (image,data);
+%variables
+if  isempty(data.layers.number)
+    n=2;
+else
+    n=data.layers.number;
+end
+if isempty(data.layers.dist)
+    dist_grid=3;
+else
+    dist_grid=data.layers.dist;
+end
+if isempty(data.W.scales)
+    s=1:n;
+elseif (length(data.W.scales)==n)
+    s=data.W.scales;
+else
+    s=1:n;
+end
+if isempty(data.W.radius)
+    r(1)=2;
+    for i=2:n
+        r(i)=10;
+    end
+else
+    r=data.W.radius;
+end
+if isempty(data.W.edgeVariance)
+    data.W.edgeVariance=0.1;
+end
+if isempty(data.W.gridtype)
+    data.W.gridtype='square';
+end
+if isempty(data.W.sigmaI)
+    data.W.sigmaI=0.12;
+end
+if isempty(data.W.sigmaX)
+    data.W.sigmaX=10;
+end
+if isempty(data.layers.weight)
+    coef(1)=5;
+    coef(2:n)=200;
+elseif (length(data.layers.weight)==n)
+    coef=data.layers.weight;
+else
+    coef(1)=5;
+    coef(2:n)=100;  %200  
+end
+if isempty(data.W.mode)
+    data.W.mode=mixed;
+end
+[p1,q1,ignore]=size(image);
+image=image(:,:,1);
+filter_par = [4,30,4];  %[9,30,4]
+[x,y,gx,gy,par,threshold,emag,ephase,g,FIe,FIo,mago] = quadedgep2(image,filter_par,data,0.001);
+minW=10^(-2); %-3
+% function [multiWpp,constraintMat,p,q,Wppp,subgrid] = computemultiWpp(n,imageX,r,dist_grid,s,dataWpp,emag,ephase,minW,mode,facteurMul,contrainte,tt,gridtype,colormode,imageOriginale,subgridImageReduite,pG,qG)
+p= p1*ones(n,1);
+q= q1*ones(n,1);
+d= dist_grid*ones(n,1);
+d(1)=1;
+for (i=2:n)
+    d(i)=d(i)*3^(i-2);
+end
+p=ceil(p1./d);
+q=ceil(q1./d);
+%computation of the subgrids (the first pixel is coded by one). S{i,j}(k) is the index of
+%the kth pixel of the jth grid in the ith grid.
+for i=1:n-1
+    for j=i+1:n
+        a=[0:p(j)*q(j)-1];
+        subgrid{i,j}=p(i)*(floor(a/p(j)))*d(j)/d(i)+(1+mod(a,p(j))*d(j)/d(i));
+    end
+end
+%computation of the independent W matrix for each layer Wind{i} 1=<i=<n.
+[w1i,w1j]=cimgnbmap([p1,q1], r(1), 1);
+if strcmp(data.W.mode,'mixed')
+    rMin = 0;
+    imageXX=double(image(:));
+    sigmaI= (std(imageXX(:)) + 1e-10 )* data.W.sigmaI;
+    Wpp{1} = multiIntensityFirstLayer(double(image),w1i,w1j,rMin,data.W.sigmaX,sigmaI,minW);
+    Wpp2= affinityic(emag(:,:,s(1)),ephase(:,:,s(1)),w1i,w1j,max(max(emag(:,:,s(1)))) * data.W.edgeVariance);
+    Wpp{1} = sqrt(Wpp{1} .* Wpp2)+0.1*Wpp2;
+    
+elseif strcmp(data.W.mode,'notmixed')
+    Wpp{1}= affinityic(emag(:,:,s(1)),ephase(:,:,s(1)),w1i,w1j,max(max(emag(:,:,s(1)))) * data.W.edgeVariance);
+    
+elseif strcmp(data.W.mode,'intensity')
+    rMin = 0;
+    imageXX=double(image(:));
+    sigmaI= (std(imageXX(:)) + 1e-10 )* data.W.sigmaI;
+    Wpp{1} = multiIntensityFirstLayer(double(image),w1i,w1j,rMin,data.W.sigmaX,sigmaI,minW);  
+    
+end
+Wpp{1}=coef(1)*(Wpp{1}+Wpp{1}')/2;
+%Wpp{1}= coef(1)*Wpp{1};
+Wind{1}=Wpp{1};
+for i=2:n 
+    if strcmp(data.W.gridtype,'square')
+        [wwi,wwj]=cimgnbmap([p(i),q(i)], r(i), 1); 
+    elseif strcmp(data.W.gridtype,'star')
+        [wwi,wwj]=cimgnbmap_star([p(i),q(i)], r(i), 1);
+    elseif strcmp(data.W.gridtype,'cross')
+        [wwi,wwj]=cimgnbmap_cross([p(i),q(i)], r(i), 1);
+    end
+    wwi=double(wwi);
+    wiInOriginalImage=(p1*(floor(wwi/p(i)))*d(i))+(1+mod(wwi,p(i))*d(i));
+        wiInOriginalImage=(p1*(floor(wwi/p(i)))*d(i))+(1+mod(wwi,p(i))*d(i));
+    wiInOriginalImage= uint32(wiInOriginalImage);
+if strcmp(data.W.mode,'mixed')
+    Wpp2= multiAffinityic(emag(:,:,i),ephase(:,:,i),wiInOriginalImage,wwj,subgrid{1,i},p(i),q(i),uint32(wwi),max(max(emag(:,:,i))) * data.W.edgeVariance);
+    a=floor(d(i)/d(i-1));
+    if (mod(a,2)==0)
+        a=a+1;
+    end
+%     Wpp{i} = multiIntensityWppc(double(imageX),wiInOriginalImage,wwj,rMin,dataWpp.sigmaX,sigmaI,minW,subgrid{1,i},p(i),q(i),wi{i});
+    
+    Wpp{i} = multiIntensityWppc(double(image),wiInOriginalImage,wwj,rMin,data.W.sigmaX,sigmaI,minW,subgrid{1,i},p(i),q(i),uint32(wwi));
+    Wpp{i} = sqrt(Wpp{i} .* Wpp2)+0.1*Wpp2;
+elseif strcmp(data.W.mode,'notmixed')
+    Wpp{i}= multiAffinityic(emag(:,:,i),ephase(:,:,i),wiInOriginalImage,wwj,subgrid{1,i},p(i),q(i),uint32(wwi),max(max(emag(:,:,i))) * data.W.edgeVariance);
+elseif strcmp(data.W.mode,'intensity')
+    Wpp{i} = multiIntensityWppc(double(image),wiInOriginalImage,wwj,rMin,data.W.sigmaX,sigmaI,minW,subgrid{1,i},p(i),q(i),uint32(wwi));
+end
+Wpp{i}= coef(i)*(Wpp{i}+Wpp{i}')/2;
+%Wpp{i}= coef(i)*Wpp{i};
+Wind{i}=Wpp{i};
+end
+%computation of the intern contraint matrices C{i,j}.
+for i=1:n-1
+    r=floor(d(i+1)/(d(i)*2));
+    [wwi,wwj]=cimgnbmap([p(i),q(i)], r, 1);
+    wi{i}=wwi;
+    wj{i}=wwj;
+end
+for i=1:n-1
+    for j=i+1:n
+        C{i,j}=sparse(p(i)*q(i),p(j)*q(j));
+        firstPointer=double(wj{i}(subgrid{i,j}))+1;   
+        lastPointer=double(wj{i}(subgrid{i,j}+1));
+        invNbNeighbours=1./(lastPointer-firstPointer+1);
+        for (k=1:p(j)*q(j))
+            for (m=firstPointer(k):lastPointer(k))
+                C{i,j}(double(wi{i}(m))+1,k)=invNbNeighbours(k);
+            end
+        end
+    end
+end
+            
+%Assembling the built matrices to make up multiWpp.
+for i=1:n
+    if (i>1)
+        for j=i-1:-1:1
+            Wpp{i}=[C{j,i}',Wpp{i}];
+        end
+    end
+    if (i<n)
+        for j=i+1:n
+            Wpp{i}=[Wpp{i},C{i,j}];
+        end
+    end
+end
+% %Assembling the built matrices to make up Wpp without intern constrains.
+% for i=1:n
+%     if (i>1)
+%         for j=i-1:-1:1
+%             Wpp{i}=[sparse(p(i)*q(i),p(j)*q(j)),Wpp{i}];
+%         end
+%     end
+%     if (i<n)
+%         for j=i+1:n
+%             Wpp{i}=[Wpp{i},sparse(p(i)*q(i),p(j)*q(j))];
+%         end
+%     end
+% end
+clear Wind;Wind = 1;
+multiWpp=Wpp{1}; clear Wpp{1}
+for i=2:n
+    multiWpp=[multiWpp;Wpp{i}];clear Wpp{i}
+end
+ 
+% Computing the average extern constraint
+ pq=sum(p(2:n).*q(2:n));
+ p2q2=p(2)*q(2);
+ constraintMat=[-C{1,2};speye(p2q2);sparse(pq-p2q2,p2q2)];
+ if n>2
+     for i=3:n
+         piqi=p(i)*q(i);
+         if i~=n
+             constraintMat=[constraintMat,[sparse(sum(p(1:i-2).*q(1:i-2)),piqi);-C{i-1,i};speye(piqi);sparse(pq-sum(p(2:i).*q(2:i)),piqi)]];
+         else
+             constraintMat=[constraintMat,[sparse(sum(p(1:i-2).*q(1:i-2)),piqi);-C{i-1,i};speye(piqi)]];
+         end
+     end
+ end
+ % saving useful information
+ %subgrids, p and q
+ data.subgrid=subgrid;
+ data.p=p;
+ data.q=q;
diff --git a/SD-VBS/common/toolbox/MultiNcut/discretisation.m b/SD-VBS/common/toolbox/MultiNcut/discretisation.m
new file mode 100755
index 0000000..70b5650
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/discretisation.m
@@ -0,0 +1,49 @@
+function [SegLabel,EigenVectors]=discretisation(EigenVectors,nr,nc)
+% 
+% EigenvectorsDiscrete=discretisation(EigenVectors)
+% 
+% Input: EigenVectors = continuous Ncut vector, size = ndata x nbEigenvectors 
+% Output EigenvectorsDiscrete = discrete Ncut vector, size = ndata x nbEigenvectors
+%
+% Timothee Cour, Stella Yu, Jianbo Shi, 2004
+[n,k]=size(EigenVectors);
+vm = sqrt(sum(EigenVectors.*EigenVectors,2));
+EigenVectors = EigenVectors./repmat(vm,1,k);
+R=zeros(k);
+R(:,1)=EigenVectors(1+round(rand(1)*(n-1)),:)';
+c=zeros(n,1);
+for j=2:k
+    c=c+abs(EigenVectors*R(:,j-1));
+    [minimum,i]=min(c);
+    R(:,j)=EigenVectors(i,:)';
+end
+lastObjectiveValue=0;
+exitLoop=0;
+nbIterationsDiscretisation = 0;
+nbIterationsDiscretisationMax = 20;%voir
+while exitLoop== 0 
+    nbIterationsDiscretisation = nbIterationsDiscretisation + 1 ;   
+    EigenvectorsDiscrete = discretisationEigenVectorData(EigenVectors*R);
+    [U,S,V] = svd(EigenvectorsDiscrete'*EigenVectors,0);    
+    NcutValue=2*(n-trace(S));
+    
+    if abs(NcutValue-lastObjectiveValue) < eps | nbIterationsDiscretisation > nbIterationsDiscretisationMax
+        exitLoop=1;
+    else
+        lastObjectiveValue = NcutValue;
+        R=V*U';
+    end
+end
+%%%%
+SegLabel = zeros(nr,nc);
+for j=1:size(EigenvectorsDiscrete,2),
+        SegLabel = SegLabel + j*reshape(EigenvectorsDiscrete(:,j),nr,nc);
+end
+EigenVectors = reshape(EigenVectors,nr,nc,size(EigenVectors,2));
diff --git a/SD-VBS/common/toolbox/MultiNcut/discretisationEigenVectorData.m b/SD-VBS/common/toolbox/MultiNcut/discretisationEigenVectorData.m
new file mode 100755
index 0000000..4626e3d
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/discretisationEigenVectorData.m
@@ -0,0 +1,12 @@
+function Y = discretisationEigenVectorData(EigenVector)
+% Y = discretisationEigenVectorData(EigenVector)
+%
+% discretizes previously rotated eigenvectors in discretisation
+% Timothee Cour, Stella Yu, Jianbo Shi, 2004
+[n,k]=size(EigenVector);
+[Maximum,J]=max(EigenVector');
+ 
+Y=sparse(1:n,J',1,n,k);    
diff --git a/SD-VBS/common/toolbox/MultiNcut/doog1.m b/SD-VBS/common/toolbox/MultiNcut/doog1.m
new file mode 100755
index 0000000..dd8e87b
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/doog1.m
@@ -0,0 +1,32 @@
+function H=doog1(sig,r,th,N);
+% H=doog1(sig,r,th,N);
+% by Serge Belongie
+no_pts=N;  % no. of points in x,y grid
+[x,y]=meshgrid(-(N/2)+1/2:(N/2)-1/2,-(N/2)+1/2:(N/2)-1/2);
+phi=pi*th/180;
+sigy=sig;
+sigx=r*sig;
+R=[cos(phi) -sin(phi); sin(phi) cos(phi)];
+C=R*diag([sigx,sigy])*R';
+X=[x(:) y(:)];
+Gb=gaussian(X,[0 0]',C);
+Gb=reshape(Gb,N,N);
+m=R*[0 sig]';
+a=1;
+b=-1;
+% make odd-symmetric filter
+Ga=gaussian(X,m/2,C);
+Ga=reshape(Ga,N,N);
+Gb=rot90(Ga,2);
+H=a*Ga+b*Gb; 
diff --git a/SD-VBS/common/toolbox/MultiNcut/doog2.m b/SD-VBS/common/toolbox/MultiNcut/doog2.m
new file mode 100755
index 0000000..a0511cb
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/doog2.m
@@ -0,0 +1,38 @@
+function G=doog2(sig,r,th,N);
+% G=doog2(sig,r,th,N);
+% Make difference of offset gaussians kernel
+% theta is in degrees
+% (see Malik & Perona, J. Opt. Soc. Amer., 1990)
+%
+% Example:
+% >> imagesc(doog2(1,12,0,64,1))
+% >> colormap(gray)
+% by Serge Belongie
+no_pts=N;  % no. of points in x,y grid
+[x,y]=meshgrid(-(N/2)+1/2:(N/2)-1/2,-(N/2)+1/2:(N/2)-1/2);
+phi=pi*th/180;
+sigy=sig;
+sigx=r*sig;
+R=[cos(phi) -sin(phi); sin(phi) cos(phi)];
+C=R*diag([sigx,sigy])*R';
+X=[x(:) y(:)];
+Gb=gaussian(X,[0 0]',C);
+Gb=reshape(Gb,N,N);
+m=R*[0 sig]';
+Ga=gaussian(X,m,C);
+Ga=reshape(Ga,N,N);
+Gc=rot90(Ga,2);
+a=-1;
+b=2;
+c=-1;
+G = a*Ga + b*Gb + c*Gc;
diff --git a/SD-VBS/common/toolbox/MultiNcut/eigSolve.m b/SD-VBS/common/toolbox/MultiNcut/eigSolve.m
new file mode 100755
index 0000000..cfb64fc
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/eigSolve.m
@@ -0,0 +1,5 @@
+function [eigenVectors,s]= eigSolve (multiWpp,ConstraintMat,data)
+[v,s] = quickNcutHardBiais2(multiWpp,ConstraintMat,data.dataGraphCut.nbEigenValues,data.dataGraphCut);
+v=v(1:data.p(1)*data.q(1),:);
+eigenVectors=reshape(v,data.p(1),data.q(1),data.dataGraphCut.nbEigenValues);
diff --git a/SD-VBS/common/toolbox/MultiNcut/fft_filt_2.m b/SD-VBS/common/toolbox/MultiNcut/fft_filt_2.m
new file mode 100755
index 0000000..9c84e96
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/fft_filt_2.m
@@ -0,0 +1,29 @@
+function FI=fft_filt_2(V,FB,sf);
+% FI=fft_filt_2(V,FB,sf);
+% fft-based filtering
+% requires image to be called "V"
+% and filters to be in FB
+% sf is the subsampling factor
+%
+% FI is the result
+[M1,M2,N3]=size(FB);
+% prepare FFT of image for filtering
+[N1,N2]=size(V);
+I=zeros(size(V)+[M1-1 M2-1]);
+I(1:N1,1:N2)=V;
+N1s=length(1:sf:N1);
+N2s=length(1:sf:N2);
+IF=fft2(I);
+FI=zeros(N1s,N2s,N3);
+% apply filters
+for n=1:N3;
+   f=rot90(FB(:,:,n),2);
+   fF=fft2(f,N1+M1-1,N2+M2-1);
+   IfF=IF.*fF;
+   If=real(ifft2(IfF));
+   If=If(ceil((M1+1)/2):ceil((M1+1)/2)+N1-1,ceil((M2+1)/2):ceil((M2+1)/2)+N2-1);
+   FI(:,:,n)=If(1:sf:N1,1:sf:N2);
+end
diff --git a/SD-VBS/common/toolbox/MultiNcut/gaussian.m b/SD-VBS/common/toolbox/MultiNcut/gaussian.m
new file mode 100755
index 0000000..509b129
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/gaussian.m
@@ -0,0 +1,31 @@
+function p=gaussian(x,m,C);
+% p=gaussian(x,m,C);
+%
+% Evaluate the multi-variate density with mean vector m and covariance
+% matrix C for the input vector x.
+% 
+% p=gaussian(X,m,C);
+% 
+% Vectorized version: Here X is a matrix of column vectors, and p is 
+% a vector of probabilities for each vector.
+d=length(m);
+if size(x,1)~=d
+   x=x';
+end
+N=size(x,2);
+detC = det(C);
+if rcond(C)<eps
+%   fprintf(1,'Covariance matrix close to singular. (gaussian.m)\n');
+   p = zeros(N,1);
+else
+   m=m(:);
+   M=m*ones(1,N);
+   denom=(2*pi)^(d/2)*sqrt(abs(detC));
+   mahal=sum(((x-M)'*inv(C)).*(x-M)',2);   % Chris Bregler's trick
+   numer=exp(-0.5*mahal);
+   p=numer/denom;
+end
diff --git a/SD-VBS/common/toolbox/MultiNcut/make_filterbank_even2.m b/SD-VBS/common/toolbox/MultiNcut/make_filterbank_even2.m
new file mode 100755
index 0000000..f7f4527
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/make_filterbank_even2.m
@@ -0,0 +1,45 @@
+function FB = make_filterbank(num_ori,filter_scales,wsz,enlong)
+%
+%  F = make_filterbank(num_ori,num_scale,wsz)
+%
+if nargin<4,
+    enlong = 3;
+end
+enlong = 2*enlong;
+% definine filterbank
+%num_ori=6;
+%num_scale=3;
+num_scale = length(filter_scales);
+M1=wsz; % size in pixels
+M2=M1;
+ori_incr=180/num_ori;
+ori_offset=ori_incr/2; % helps with equalizing quantiz. error across filter set
+FB=zeros(M1,M2,num_ori,num_scale);
+% elongated filter set
+counter = 1;
+for m=1:num_scale  
+   for n=1:num_ori
+      % r=12 here is equivalent to Malik's r=3;
+      f=doog2(filter_scales(m),enlong,ori_offset+(n-1)*ori_incr,M1);
+      FB(:,:,n,m)=f;
+   end
+end
+FB=reshape(FB,M1,M2,num_scale*num_ori);
+total_num_filt=size(FB,3);
+for j=1:total_num_filt,
+  F = FB(:,:,j);
+  a = sum(sum(abs(F)));
+  FB(:,:,j) = FB(:,:,j)/a;
+end
diff --git a/SD-VBS/common/toolbox/MultiNcut/make_filterbank_odd2.m b/SD-VBS/common/toolbox/MultiNcut/make_filterbank_odd2.m
new file mode 100755
index 0000000..0059dca
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/make_filterbank_odd2.m
@@ -0,0 +1,46 @@
+function FB = make_filterbank(num_ori,filter_scales,wsz,enlong)
+%
+%  F = make_filterbank(num_ori,num_scale,wsz)
+%
+if nargin<4,
+    enlong = 3;
+end
+enlong = enlong*2;
+% definine filterbank
+%num_ori=6;
+%num_scale=3;
+num_scale = length(filter_scales);
+M1=wsz; % size in pixels
+M2=M1;
+ori_incr=180/num_ori;
+ori_offset=ori_incr/2; % helps with equalizing quantiz. error across filter set
+FB=zeros(M1,M2,num_ori,num_scale);
+% elongated filter set
+counter = 1;
+for m=1:num_scale  
+   for n=1:num_ori
+      % r=12 here is equivalent to Malik's r=3;
+      f=doog1(filter_scales(m),enlong,ori_offset+(n-1)*ori_incr,M1);
+      FB(:,:,n,m)=f;
+   end
+end
+FB=reshape(FB,M1,M2,num_scale*num_ori);
+total_num_filt=size(FB,3);
+for j=1:total_num_filt,
+  F = FB(:,:,j);
+  a = sum(sum(abs(F)));
+  FB(:,:,j) = FB(:,:,j)/a;
+end
diff --git a/SD-VBS/common/toolbox/MultiNcut/mex_projection_QR.c b/SD-VBS/common/toolbox/MultiNcut/mex_projection_QR.c
new file mode 100755
index 0000000..5ac4820
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/mex_projection_QR.c
@@ -0,0 +1,82 @@
+/*================================================================
+* mex_projection_QR.c = used by quickNcutHardBiais.m in eigensolver.
+* mex_projection_QR(X,P,Ubar,R) == (eye(length(P))-Ubar*(R'*R)^{-1}*Ubar') * P * (eye(length(P))-Ubar*(R'*R)^{-1}*Ubar') * X ;
+* (R'*R)^{-1} is solved by solving a triangular system
+* P, Ubar, R are sparse, but X is full
+* R is upper triangular
+% test sequence:
+*=================================================================*/
+# include "math.h"
+# include "mex.h"
+# include "a_times_b_cmplx.c"
+/*# include "a_times_b.c"*/
+void mexFunction(
+    int nargout,
+    mxArray *out[], 
+    int nargin,
+    const mxArray *in[]
+)
+{
+        int *ir[4], *jc[4], m[4], n[4];
+        double *y, *y1,*y2,*y3,*y4,*y5,*y6, *pr[4];
+        double *y2bis, *y5bis;    
+        int k;  
+        for (k=0; k<4; k++) {
+                m[k] = mxGetM(in[k]);
+                n[k] = mxGetN(in[k]);
+                pr[k] = mxGetPr(in[k]);
+                if (k>0) {
+                        ir[k] = mxGetIr(in[k]);
+                        jc[k] = mxGetJc(in[k]);
+                }
+        }
+        out[0] = mxCreateDoubleMatrix(m[1],1,mxREAL);
+        y = mxGetPr(out[0]);
+    
+        y1 = mxCalloc(n[2], sizeof(double));
+        y2 = mxCalloc(m[3], sizeof(double));    
+        y2bis = mxCalloc(m[3], sizeof(double));
+        y3 = mxCalloc(m[1], sizeof(double));
+        y4 = mxCalloc(m[1], sizeof(double));
+        y5 = mxCalloc(n[2], sizeof(double));
+        y5bis = mxCalloc(n[2], sizeof(double));
+        y6 = mxCalloc(n[2], sizeof(double));
+        CSR_VecMult_CAB_double(m[2],n[2],pr[2],ir[2],jc[2],pr[0],y1);
+        CSR_VecTriangSlvLD_CAB_double(m[3],pr[3],ir[3],jc[3],y1,y2bis);
+        CSC_VecTriangSlvUD_CAB_double(m[3],pr[3],ir[3],jc[3],y2bis,y2);
+        for(k=0;k<m[1];k++) {
+                y3[k]=pr[0][k];
+        }
+        CSC_VecMult_CaABC_double(m[2],n[2],-1,pr[2],ir[2],jc[2],y2,y3);
+        
+        CSR_VecMult_CAB_double(m[1],n[1],pr[1],ir[1],jc[1],y3,y4);
+        CSR_VecMult_CAB_double(m[2],n[2],pr[2],ir[2],jc[2],y4,y5);
+        
+        CSR_VecTriangSlvLD_CAB_double(m[3],pr[3],ir[3],jc[3],y5,y5bis);
+        CSC_VecTriangSlvUD_CAB_double(m[3],pr[3],ir[3],jc[3],y5bis,y6);
+        for(k=0;k<m[1];k++) {
+                y[k]=y4[k];
+        }
+        CSC_VecMult_CaABC_double(m[2],n[2],-1,pr[2],ir[2],jc[2],y6,y);
+        mxFree(y1);
+        mxFree(y2);
+        mxFree(y2bis);
+        mxFree(y3);
+        mxFree(y4);
+        mxFree(y5);
+        mxFree(y5bis);
+        mxFree(y6);
+} 
diff --git a/SD-VBS/common/toolbox/MultiNcut/mex_projection_QR.dll b/SD-VBS/common/toolbox/MultiNcut/mex_projection_QR.dll
new file mode 100755
index 0000000..d77f509
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/mex_projection_QR.dll
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/mex_projection_QR.mexa64 b/SD-VBS/common/toolbox/MultiNcut/mex_projection_QR.mexa64
new file mode 100755
index 0000000..bb7b7bd
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/mex_projection_QR.mexa64
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/mex_projection_QR.mexglx b/SD-VBS/common/toolbox/MultiNcut/mex_projection_QR.mexglx
new file mode 100755
index 0000000..dd63169
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/mex_projection_QR.mexglx
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/mex_projection_QR_symmetric.c b/SD-VBS/common/toolbox/MultiNcut/mex_projection_QR_symmetric.c
new file mode 100755
index 0000000..5915959
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/mex_projection_QR_symmetric.c
@@ -0,0 +1,83 @@
+/*================================================================
+* mex_projection_QR_symmetric.c = used by quickNcutHardBiais.m in eigensolver.
+* mex_projection_QR_symmetric(X,P,Ubar,R) == (eye(length(P))-Ubar*(R'*R)^{-1}*Ubar') * P * (eye(length(P))-Ubar*(R'*R)^{-1}*Ubar') * X ;
+* (R'*R)^{-1} is solved by solving a triangular system
+* P, Ubar, R are sparse, but X is full
+* R is upper triangular
+* P is symmetric
+% test sequence:
+*=================================================================*/
+# include "math.h"
+# include "mex.h"
+# include "a_times_b_cmplx.c"
+/*# include "a_times_b.c"*/
+void mexFunction(
+    int nargout,
+    mxArray *out[], 
+    int nargin,
+    const mxArray *in[]
+)
+{
+        int *ir[4], *jc[4], m[4], n[4];
+        double *y, *y1,*y2,*y3,*y4,*y5,*y6, *pr[4];
+        double *y2bis, *y5bis;    
+        int k;  
+        for (k=0; k<4; k++) {
+                m[k] = mxGetM(in[k]);
+                n[k] = mxGetN(in[k]);
+                pr[k] = mxGetPr(in[k]);
+                if (k>0) {
+                        ir[k] = mxGetIr(in[k]);
+                        jc[k] = mxGetJc(in[k]);
+                }
+        }
+        out[0] = mxCreateDoubleMatrix(m[1],1,mxREAL);
+        y = mxGetPr(out[0]);
+    
+        y1 = mxCalloc(n[2], sizeof(double));
+        y2 = mxCalloc(m[3], sizeof(double));    
+        y2bis = mxCalloc(m[3], sizeof(double));
+        y3 = mxCalloc(m[1], sizeof(double));
+        y4 = mxCalloc(m[1], sizeof(double));
+        y5 = mxCalloc(n[2], sizeof(double));
+        y5bis = mxCalloc(n[2], sizeof(double));
+        y6 = mxCalloc(n[2], sizeof(double));
+        CSR_VecMult_CAB_double(m[2],n[2],pr[2],ir[2],jc[2],pr[0],y1);
+        CSR_VecTriangSlvLD_CAB_double(m[3],pr[3],ir[3],jc[3],y1,y2bis);
+        CSC_VecTriangSlvUD_CAB_double(m[3],pr[3],ir[3],jc[3],y2bis,y2);
+        for(k=0;k<m[1];k++) {
+                y3[k]=pr[0][k];
+        }
+        CSC_VecMult_CaABC_double(m[2],n[2],-1,pr[2],ir[2],jc[2],y2,y3);
+        
+        CSRsymm_VecMult_CAB_double(m[1],n[1],pr[1],ir[1],jc[1],y3,y4);
+        CSR_VecMult_CAB_double(m[2],n[2],pr[2],ir[2],jc[2],y4,y5);
+        
+        CSR_VecTriangSlvLD_CAB_double(m[3],pr[3],ir[3],jc[3],y5,y5bis);
+        CSC_VecTriangSlvUD_CAB_double(m[3],pr[3],ir[3],jc[3],y5bis,y6);
+        for(k=0;k<m[1];k++) {
+                y[k]=y4[k];
+        }
+        CSC_VecMult_CaABC_double(m[2],n[2],-1,pr[2],ir[2],jc[2],y6,y);
+        mxFree(y1);
+        mxFree(y2);
+        mxFree(y2bis);
+        mxFree(y3);
+        mxFree(y4);
+        mxFree(y5);
+        mxFree(y5bis);
+        mxFree(y6);
+} 
diff --git a/SD-VBS/common/toolbox/MultiNcut/mex_projection_QR_symmetric.dll b/SD-VBS/common/toolbox/MultiNcut/mex_projection_QR_symmetric.dll
new file mode 100755
index 0000000..0270d17
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/mex_projection_QR_symmetric.dll
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/mex_projection_QR_symmetric.mexa64 b/SD-VBS/common/toolbox/MultiNcut/mex_projection_QR_symmetric.mexa64
new file mode 100755
index 0000000..201f6e3
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/mex_projection_QR_symmetric.mexa64
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/mex_projection_QR_symmetric.mexglx b/SD-VBS/common/toolbox/MultiNcut/mex_projection_QR_symmetric.mexglx
new file mode 100755
index 0000000..5b80e13
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/mex_projection_QR_symmetric.mexglx
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/mex_w_times_x_symmetric.mexglx b/SD-VBS/common/toolbox/MultiNcut/mex_w_times_x_symmetric.mexglx
new file mode 100755
index 0000000..bfad956
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/mex_w_times_x_symmetric.mexglx
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/mex_w_times_x_symmetric.mexmac b/SD-VBS/common/toolbox/MultiNcut/mex_w_times_x_symmetric.mexmac
new file mode 100755
index 0000000..54be540
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/mex_w_times_x_symmetric.mexmac
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/multiAffinityic.c b/SD-VBS/common/toolbox/MultiNcut/multiAffinityic.c
new file mode 100755
index 0000000..6919db9
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/multiAffinityic.c
@@ -0,0 +1,216 @@
+/*================================================================
+* function w = affinityic(emag,ephase,pi,pj,subgrid,nrSubgrid,ncSubgrid,subpi,sigma)
+* Input:
+*   emag = edge strength at each pixel
+*   ephase = edge phase at each pixel
+*   [pi,pj] = index pair representation for MALTAB sparse matrices, pi index in original grid, size(pj)=nrW*ncW+1
+*   subgrid= index of the subgrid in the original image (first pixel's index is one!)
+*   [nrSubgrid,ncSubgrid]= number of rows et colums of the subgrid
+*   subpi = pi but with index in subgrid
+*   sigma = sigma for IC energy
+* Output:
+*   w = affinity with IC at [pi,pj]
+*
+% test sequence
+f = synimg(10);
+[i,j] = cimgnbmap(size(f),2);
+[ex,ey,egx,egy] = quadedgep(f);
+a = affinityic(ex,ey,egx,egy,i,j)
+show_dist_w(f,a);
+* Stella X. Yu, Nov 19, 2001.
+*=================================================================*/
+# include "mex.h"
+# include "math.h"
+void mexFunction(
+    int nargout,
+    mxArray *out[],
+    int nargin,
+    const mxArray *in[]
+)
+{
+    /* declare variables */
+    int nr, nc, np,  nW, total;
+    int i, j, k, t, ix, iy, jx, jy, ii, jj, iip1, jjp1, iip2, jjp2, step,nrSubgrid, ncSubgrid;
+    double sigma, di, dj, a, z, maxori, phase1, phase2, slope;
+        int *ir, *jc;
+    /*  unsigned long *pi, *pj, *subpi;
+ */
+    unsigned int *pi, *pj, *subpi;
+        double *emag, *ephase, *w,*tmp,*subgrid;
+        
+    
+    /* check argument */
+    if (nargin<8) {
+        mexErrMsgTxt("Eight input arguments required");
+    }
+    if (nargout>1) {
+        mexErrMsgTxt("Too many output arguments");
+    }
+    /* get edgel information */
+        nr = mxGetM(in[0]);
+        nc = mxGetN(in[0]);
+        if ( nr*nc ==0 || nr != mxGetM(in[1]) || nc != mxGetN(in[1]) ) {
+            mexErrMsgTxt("Edge magnitude and phase shall be of the same image size");
+        }
+    emag = mxGetPr(in[0]);
+    ephase = mxGetPr(in[1]);
+    np = nr * nc;
+    
+    /*get subgrid information*/
+    
+    tmp = mxGetData(in[5]);
+    nrSubgrid = (int)tmp[0];
+    tmp = mxGetData(in[6]);
+    ncSubgrid = (int)tmp[0];
+    
+    /* printf("nrSubgrid=%d\n",nrSubgrid);
+ 
+       printf("ncSubgrid=%d\n",ncSubgrid);
+ */
+    if (nrSubgrid* ncSubgrid != mxGetM(in[4])*mxGetN(in[4])) {
+          mexErrMsgTxt("Error in the size of the subgrid");
+          }
+    subgrid = mxGetData(in[4]);
+    nW = nrSubgrid * ncSubgrid;
+    
+    /* get new index pair */
+    if (!mxIsUint32(in[2]) | !mxIsUint32(in[3])) {
+        mexErrMsgTxt("Index pair shall be of type UINT32");
+    }
+    if (mxGetM(in[3]) * mxGetN(in[3]) != nW + 1) {
+        mexErrMsgTxt("Wrong index representation");
+    }
+    pi = mxGetData(in[2]);
+    pj = mxGetData(in[3]);
+    subpi = mxGetData(in[7]);
+    /*{printf("pi[50] = %d\n",pi[50]);}
+    {printf("subpi[5] = %d\n",subpi[5]);}
+    {printf("subpi[6] = %d\n",subpi[6]);}
+    {printf("subpi[4] = %d\n",subpi[4]);}*/
+    /* create output */          /*!!!!!!!!!!!!!!!!!!!!!!!changer taille output!!!!!!!!!!*/
+    out[0] = mxCreateSparse(nW,nW,pj[nW],mxREAL);
+    if (out[0]==NULL) {
+            mexErrMsgTxt("Not enough memory for the output matrix");
+        }
+        w = mxGetPr(out[0]);
+        ir = mxGetIr(out[0]);
+        jc = mxGetJc(out[0]);
+        
+    /* find my sigma */
+        if (nargin<9) {
+            sigma = 0;
+        for (k=0; k<np; k++) { 
+            if (emag[k]>sigma) { sigma = emag[k]; }
+        }
+        sigma = sigma / 6;
+        printf("sigma = %6.5f",sigma);
+        } else {
+            sigma = mxGetScalar(in[8]);
+        }
+        a = 0.5 / (sigma * sigma);
+        
+    /* computation */ 
+    total = 0;
+    
+    for (j=0; j<nW; j++){
+        t= (int)subgrid[j]-1;  /*on parcourt tous les pixels de la sous-grille dans la grille d'origine*/        
+        jc[j] = total;    /* total represente le nombre voisins  du pixel j*/
+        jx = t / nr; /* col */  
+        jy = t % nr; /* row */
+        
+        for (k=pj[j]; k<pj[j+1]; k++) {   /*k represente les indices correspondant au pixel j dans pi*/
+        
+            i = pi[k]-1;                     /*i est un voisin de j a considerer*/
+          
+            if (i==j) {
+                maxori = 1;
+            
+            } else {
+        
+                ix = i / nr; 
+                iy = i % nr;
+                /* scan */            
+                di = (double) (iy - jy);
+                dj = (double) (ix - jx);
+            
+                maxori = 0.;
+                    phase1 = ephase[j];
+                       
+                /* sample in i direction */
+                if (abs(di) >= abs(dj)) {  
+                    slope = dj / di;
+                    step = (iy>=jy) ? 1 : -1;
+                
+                    iip1 = jy;
+                    jjp1 = jx;
+        
+                       
+                        for (ii=0;ii<abs(di);ii++){
+                            iip2 = iip1 + step;
+                            jjp2 = (int)(0.5 + slope*(iip2-jy) + jx);
+                  
+                            phase2 = ephase[iip2+jjp2*nr];
+               
+                            if (phase1 != phase2) {
+                                z = (emag[iip1+jjp1*nr] + emag[iip2+jjp2*nr]);
+                                if (z > maxori){
+                                    maxori = z;
+                                }
+                            } 
+                     
+                            iip1 = iip2;
+                            jjp1 = jjp2;
+                            phase1 = phase2;
+                        }
+                    
+                    /* sample in j direction */    
+                } else { 
+                        slope = di / dj;
+                        step =  (ix>=jx) ? 1: -1;
+                    jjp1 = jx;
+                        iip1 = jy;                 
+            
+         
+                        for (jj=0;jj<abs(dj);jj++){
+                            jjp2 = jjp1 + step;
+                            iip2 = (int)(0.5+ slope*(jjp2-jx) + jy);
+                  
+                            phase2 = ephase[iip2+jjp2*nr];
+             
+                            if (phase1 != phase2){
+                                z = (emag[iip1+jjp1*nr] + emag[iip2+jjp2*nr]);
+                                if (z > maxori){ 
+                                    maxori = z; 
+                                }
+                                
+                            }
+          
+                            iip1 = iip2;
+                            jjp1 = jjp2;
+                            phase1 = phase2;
+                        }
+                }            
+            
+                maxori = 0.5 * maxori;
+                maxori = exp(-maxori * maxori * a);
+            }  
+            /*if (total<20) {printf("subpi[k] = %d\n",subpi[k]);}*/
+                    ir[total] = (int)subpi[k];
+/*if (total<20) {printf("ir[total] = %d\n",ir[total]);}*/
+                    w[total] = maxori;
+                    total = total + 1;
+                        
+                } /* i */
+       }  /*j*/
+     /*printf("total = %d\n",total);*/   
+      /*printf("ir[100] = %d\n",ir[100]);*/ 
+    jc[nW] = total;
+}  
diff --git a/SD-VBS/common/toolbox/MultiNcut/multiAffinityic.dll b/SD-VBS/common/toolbox/MultiNcut/multiAffinityic.dll
new file mode 100755
index 0000000..f9b6bc9
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/multiAffinityic.dll
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/multiAffinityic.mexa64 b/SD-VBS/common/toolbox/MultiNcut/multiAffinityic.mexa64
new file mode 100755
index 0000000..9a2a3cc
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/multiAffinityic.mexa64
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/multiAffinityic.mexglx b/SD-VBS/common/toolbox/MultiNcut/multiAffinityic.mexglx
new file mode 100755
index 0000000..ec11a0e
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/multiAffinityic.mexglx
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/multiIntensityFirstLayer.c b/SD-VBS/common/toolbox/MultiNcut/multiIntensityFirstLayer.c
new file mode 100755
index 0000000..0cae523
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/multiIntensityFirstLayer.c
@@ -0,0 +1,126 @@
+/*================================================================
+* function w = intensityWppc(image,pi,pj,rMin,sigmaX,sigmaIntensite,valeurMinW )
+* Input:
+*   [pi,pj] = index pair representation for MALTAB sparse matrices
+* Output:
+*   w = affinity with IC at [pi,pj]
+*
+pixels i and j (corresponding to the sampling in pi,pj) are fully connected when d(i,j) <= rmin;
+% test sequence
+f = synimg(10);
+[i,j] = cimgnbmap(size(f),2);
+[ex,ey,egx,egy] = quadedgep(f);
+a = affinityic(ex,ey,egx,egy,i,j)
+show_dist_w(f,a);
+* Stella X. Yu, Nov 19, 2001.
+*=================================================================*/
+# include "mex.h"
+# include "math.h"
+void mexFunction(
+    int nargout,
+    mxArray *out[],
+    int nargin,
+    const mxArray *in[]
+)
+{
+    /* declare variables */
+    int nr, nc, np, total;
+    int i, j, k, ix, iy, jx, jy;
+        int *ir, *jc;
+        int squareDistance;
+    /* unsigned long *pi, *pj; */
+unsigned int *pi, *pj;
+        double *w;
+        double temp,a1,a2,wij;
+        double rMin;
+        double sigmaX, sigmaIntensite,valeurMinW;
+        double *image;
+    
+    /* check argument */
+    if (nargin<7) {
+        mexErrMsgTxt("Four input arguments required");
+    }
+    if (nargout>1) {
+        mexErrMsgTxt("Too many output arguments");
+    }
+    /* get edgel information */
+        nr = mxGetM(in[0]);
+        nc = mxGetN(in[0]);
+        np = nr * nc;
+    image = mxGetPr(in[0]);
+    
+    /* get new index pair */
+    if (!mxIsUint32(in[1]) | !mxIsUint32(in[2])) {
+        mexErrMsgTxt("Index pair shall be of type UINT32");
+    }
+    if (mxGetM(in[2]) * mxGetN(in[2]) != np + 1) {
+        mexErrMsgTxt("Wrong index representation");
+    }
+    pi = mxGetData(in[1]);
+    pj = mxGetData(in[2]);    
+    /* create output */
+    out[0] = mxCreateSparse(np,np,pj[np],mxREAL);
+    if (out[0]==NULL) {
+            mexErrMsgTxt("Not enough memory for the output matrix");
+        }
+        w = mxGetPr(out[0]);
+        ir = mxGetIr(out[0]);
+        jc = mxGetJc(out[0]);
+        
+        rMin = mxGetScalar(in[3]);
+        sigmaX = mxGetScalar(in[4]);
+        sigmaIntensite= mxGetScalar(in[5]);
+        valeurMinW = mxGetScalar(in[6]);
+        a1 = 1.0/ (sigmaX*sigmaX);
+        a2 = 1.0 / (sigmaIntensite*sigmaIntensite );
+    /* computation */ 
+    total = 0;
+    for (j=0; j<np; j++) {            
+        jc[j] = total;
+        jx = j / nr; /* col */
+        jy = j % nr; /* row */
+        
+        for (k=pj[j]; k<pj[j+1]; k++) {
+        
+            i = pi[k];
+          
+            if (i==j) {
+                wij= 1; /*voir*/
+            
+            } else {        
+                ix = i / nr; 
+                iy = i % nr;
+                squareDistance = (ix-jx)*(ix-jx)+(iy-jy)*(iy-jy);
+                
+                   temp = image[i]-image[j];
+                   wij = exp(- squareDistance * a1 - temp*temp * a2 );
+                   /*if(wij < valeurMinW)
+                      wij = -0.1;*/
+            }
+                ir[total] = i;
+                w[total] = wij;
+                total = total + 1;
+                        
+          } /* i */
+    } /* j */
+        
+    jc[np] = total;
+}
+  
diff --git a/SD-VBS/common/toolbox/MultiNcut/multiIntensityFirstLayer.dll b/SD-VBS/common/toolbox/MultiNcut/multiIntensityFirstLayer.dll
new file mode 100755
index 0000000..a47e22d
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/multiIntensityFirstLayer.dll
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/multiIntensityFirstLayer.mexa64 b/SD-VBS/common/toolbox/MultiNcut/multiIntensityFirstLayer.mexa64
new file mode 100755
index 0000000..06b2f57
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/multiIntensityFirstLayer.mexa64
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/multiIntensityFirstLayer.mexglx b/SD-VBS/common/toolbox/MultiNcut/multiIntensityFirstLayer.mexglx
new file mode 100755
index 0000000..c6eed07
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/multiIntensityFirstLayer.mexglx
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/multiIntensityWppc.c b/SD-VBS/common/toolbox/MultiNcut/multiIntensityWppc.c
new file mode 100755
index 0000000..e01e516
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/multiIntensityWppc.c
@@ -0,0 +1,158 @@
+/*================================================================
+* function w = multiIntensityWppc(image,pi,pj,rMin,sigmaX,sigmaIntensite,valeurMinW,
+*                                 subgrid,nrSubgrid,ncSubgrid,subpi)
+* Input:
+*   [pi,pj] = index pair representation for MALTAB sparse matrices
+* Output:
+*   w = affinity with IC at [pi,pj]
+*
+imageX,wiInOriginalImage,wwj,rMin,dataWpp.sigmaX,sigmaI,minW,subgrid{1,i},p(i),q(i),wi{i}
+pixels i and j (corresponding to the sampling in pi,pj) are fully connected when d(i,j) <= rmin;
+% test sequence
+f = synimg(10);
+[i,j] = cimgnbmap(size(f),2);
+[ex,ey,egx,egy] = quadedgep(f);
+a = affinityic(ex,ey,egx,egy,i,j)
+show_dist_w(f,a);
+* Stella X. Yu, Nov 19, 2001.
+*=================================================================*/
+# include "mex.h"
+# include "math.h"
+void mexFunction(
+    int nargout,
+    mxArray *out[],
+    int nargin,
+    const mxArray *in[]
+)
+{
+    /* declare variables */
+    int nr, nc, np, nW, total;
+    int i, j, k, t, ix, iy, jx, jy, nrSubgrid, ncSubgrid;
+        int *ir, *jc;
+        int squareDistance;
+    /*  unsigned long *pi, *pj, *subpi; */
+        unsigned int *pi, *pj, *subpi;
+        double *w, *subgrid, *tmp;
+        double temp,a1,a2,wij;
+        double rMin;
+        double sigmaX, sigmaIntensite,valeurMinW;
+        double *image;
+    
+    /* check argument */
+    if (nargin<11) {
+        mexErrMsgTxt("Eleven input arguments required");
+    }
+    if (nargout>1) {
+        mexErrMsgTxt("Too many output arguments");
+    }
+    /* get edgel information */
+        nr = mxGetM(in[0]);
+        nc = mxGetN(in[0]);
+        np = nr * nc;
+    /*printf("size: %d, %d, %d\n", nc, nr, np); */
+    image = mxGetPr(in[0]);
+ 
+    /*get subgrid information*/
+    
+    tmp = mxGetData(in[8]);
+    nrSubgrid = (int)tmp[0];
+    /* printf("image end = %f ", image[np-1]); */
+    
+    tmp = mxGetData(in[9]);
+    ncSubgrid = (int)tmp[0];
+    
+    if (nrSubgrid* ncSubgrid != mxGetM(in[7])*mxGetN(in[7])) {
+          mexErrMsgTxt("Error in the size of the subgrid");
+          }
+    subgrid = mxGetData(in[7]);
+    nW = nrSubgrid * ncSubgrid;
+    
+    
+    
+    
+    /* get new index pair */
+    if (!mxIsUint32(in[1]) | !mxIsUint32(in[2])) {
+        mexErrMsgTxt("Index pair shall be of type UINT32");
+    }
+    if (mxGetM(in[2]) * mxGetN(in[2]) != nW + 1) {
+        mexErrMsgTxt("Wrong index representation");
+    }
+    pi = mxGetData(in[1]);
+    pj = mxGetData(in[2]);    
+    subpi = mxGetData(in[10]);
+    
+    /* create output */
+    out[0] = mxCreateSparse(nW,nW,pj[nW],mxREAL);
+    if (out[0]==NULL) {
+            mexErrMsgTxt("Not enough memory for the output matrix");
+        }
+        
+        w = mxGetPr(out[0]);
+        ir = mxGetIr(out[0]);
+        jc = mxGetJc(out[0]);
+        
+        
+        rMin = mxGetScalar(in[3]);
+        sigmaX = mxGetScalar(in[4]);
+        sigmaIntensite= mxGetScalar(in[5]);
+        valeurMinW = mxGetScalar(in[6]);
+        a1 = 1.0/ (sigmaX*sigmaX);
+        a2 = 1.0 / (sigmaIntensite*sigmaIntensite );
+ 
+ 
+    /* computation */ 
+    total = 0;
+    for (j=0; j<nW; j++) { 
+     
+        t= (int)subgrid[j]-1;  /*on parcourt tous les pixels de la sous-grille dans la grille d'origine*/ 
+        if ( (t<0) || (t>np-1)) {printf("badddddd!");}
+        /*     printf("t = %d\n",t);
+               printf("j = %d\n",j); */
+        jc[j] = total;
+        jx = t / nr; /* col */
+        jy = t % nr; /* row */
+     /*  printf("pj[j+1] = %d\n",pj[j+1]); */
+       
+        for (k=pj[j]; k<pj[j+1]; k++) {
+         /* printf("k = %d\n",k); */
+            i = pi[k]-1;
+            ix = i / nr; 
+            iy = i % nr;
+            squareDistance = (ix-jx)*(ix-jx)+(iy-jy)*(iy-jy);/*abs(ix-jx)+abs(iy-jy);*/
+            if (squareDistance <= rMin) { wij = 1;}
+            else {
+              temp = image[i]-image[t];
+           wij = exp(- squareDistance * a1 - temp*temp * a2 );
+                   /*if(wij < valeurMinW)
+                      wij = 0;*/
+              /*wij = exp( - temp*temp * a2 );*/
+                    }
+          
+            ir[total] = (int)subpi[k];
+            
+           /* if (ir[total] >5000*5000 ) {printf("trouble! [%d,%d]\n",k,(int)subpi[k]);} */
+                
+                     w[total] = wij;
+                   
+                     total = total + 1;
+                                
+          } /* i */
+ 
+    } /* j */
+        
+    jc[nW] = total;
+}  
diff --git a/SD-VBS/common/toolbox/MultiNcut/multiIntensityWppc.dll b/SD-VBS/common/toolbox/MultiNcut/multiIntensityWppc.dll
new file mode 100755
index 0000000..16146c1
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/multiIntensityWppc.dll
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/multiIntensityWppc.mexa64 b/SD-VBS/common/toolbox/MultiNcut/multiIntensityWppc.mexa64
new file mode 100755
index 0000000..525ca62
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/multiIntensityWppc.mexa64
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/multiIntensityWppc.mexglx b/SD-VBS/common/toolbox/MultiNcut/multiIntensityWppc.mexglx
new file mode 100755
index 0000000..a17fa03
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/multiIntensityWppc.mexglx
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/quadedgep2.m b/SD-VBS/common/toolbox/MultiNcut/quadedgep2.m
new file mode 100755
index 0000000..5041377
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/quadedgep2.m
@@ -0,0 +1,188 @@
+% function [xs,ys,gx,gy,par,threshold,mag,mage,g,FIe,FIo,mago] = quadedgep(I,par,threshold);
+% Input:
+%    I = image
+%    par = vector for 4 parameters
+%      [number of filter orientations, number of scale, filter size, elongation]
+%      To use default values, put 0.
+%    threshold = threshold on edge strength
+% Output:
+%    [x,y,gx,gy] = locations and gradients of an ordered list of edgels
+%       x,y could be horizontal or vertical or 45 between pixel sites
+%       but it is guaranteed that there [floor(y) + (floor(x)-1)*nr] 
+%       is ordered and unique.  In other words, each edgel has a unique pixel id.
+%    par = actual par used
+%    threshold = actual threshold used
+%    mag = edge magnitude
+%    mage = phase map
+%    g = gradient map at each pixel
+%    [FIe,FIo] = odd and even filter outputs
+%    mago = odd filter output of optimum orientation
+% Stella X. Yu, 2001
+% This is the multi scale version of the filtering
+% For the moment the parameters are defined by default. See line 34
+function [x,y,gx,gy,par,threshold,mag_s,mage,g,FIe,FIo,mago] = quadedgep2(I,par,data,threshold);
+if nargin<4 | isempty(threshold),
+    threshold = 0.1;
+end
+[r,c] = size(I);
+def_par = [4,30,3];
+display_on = 1;
+% take care of parameters, any missing value is substituted by a default value
+if nargin<2 | isempty(par),
+   par = def_par;
+end
+% par(end+1:4)=0;
+% par = par(:);
+% j = (par>0);
+% have_value = [ j, 1-j ];
+% j = 1; n_filter = have_value(j,:) * [par(j); def_par(j)];
+% j = 2; n_scale  = have_value(j,:) * [par(j); def_par(j)];
+% j = 3; winsz    = have_value(j,:) * [par(j); def_par(j)];
+% j = 4; enlong   = have_value(j,:) * [par(j); def_par(j)];
+n_ori = par(1);      %if it doesn't work, par<-def_par
+winsz = par(2);
+enlong = par(3);
+% always make filter size an odd number so that the results will not be skewed
+j = winsz/2;
+if not(j > fix(j) + 0.1),
+    winsz = winsz + 1;
+end
+% filter the image with quadrature filters 
+if (isempty(data.W.scales))
+        error ('no scales entered');
+end
+n_scale=length(data.W.scales);
+filter_scales=data.W.scales;
+%     
+% if strcmp(data.dataWpp.mode,'multiscale')
+%     n_scale=size((data.dataWpp.scales),2);
+%     filter_scales=data.dataWpp.scales;
+% else
+%     filter_scales=data.dataWpp.scales(1);
+%     n_scale=1;
+% end
+% if n_scale>0&&strcmp(data.dataWpp.mode,'multiscale')
+%     if (~isempty(data.dataWpp.scales))
+%         filter_scales=data.dataWpp.scales;
+%     else
+%         filter_scales=zeros(1,n_scale);
+%         
+%         for i=1:n_scale,
+%         filter_scales(i)=(sqrt(2))^(i-1);
+%         end
+%         data.dataWpp.scales=filter_scales;
+%     end
+% else filter_scale=1;
+%     data.dataWpp.scales=filter_scales;
+% end
+% 
+% %%%%%%% juste pour que ca tourne
+% if ~strcmp(data.dataWpp.mode,'multiscale')
+%     filter_scales=data.dataWpp.scales(1);
+%     n_scale=4;
+% end
+% %%%%%%%%%%%%
+FBo = make_filterbank_odd2(n_ori,filter_scales,winsz,enlong);
+FBe = make_filterbank_even2(n_ori,filter_scales,winsz,enlong);
+n = ceil(winsz/2);
+f = [fliplr(I(:,2:n+1)), I, fliplr(I(:,c-n:c-1))];
+f = [flipud(f(2:n+1,:)); f; flipud(f(r-n:r-1,:))];
+FIo = fft_filt_2(f,FBo,1); 
+FIo = FIo(n+[1:r],n+[1:c],:);
+FIe = fft_filt_2(f,FBe,1);
+FIe = FIe(n+[1:r],n+[1:c],:);
+% compute the orientation energy and recover a smooth edge map
+% pick up the maximum energy across scale and orientation
+% even filter's output: as it is the second derivative, zero cross localize the edge
+% odd filter's output: orientation
+[nr,nc,nb] = size(FIe);
+FIe = reshape(FIe, nr,nc,n_ori,length(filter_scales));
+FIo = reshape(FIo, nr,nc,n_ori,length(filter_scales));
+mag_s = zeros(nr,nc,n_scale);
+mag_a = zeros(nr,nc,n_ori);
+mage = zeros(nr,nc,n_scale);
+mago = zeros(nr,nc,n_scale);
+mage = zeros(nr,nc,n_scale);
+mago = zeros(nr,nc,n_scale);
+for i = 1:n_scale,
+    mag_s(:,:,i) = sqrt(sum(FIo(:,:,:,i).^2,3)+sum(FIe(:,:,:,i).^2,3));
+    mag_a = sqrt(FIo(:,:,:,i).^2+FIe(:,:,:,i).^2);
+    [tmp,max_id] = max(mag_a,[],3);
+ 
+    base_size = nr * nc;
+    id = [1:base_size]';
+    mage(:,:,i) = reshape(FIe(id+(max_id(:)-1)*base_size+(i-1)*base_size*n_ori),[nr,nc]);
+    mago(:,:,i) = reshape(FIo(id+(max_id(:)-1)*base_size+(i-1)*base_size*n_ori),[nr,nc]);
+    
+    mage(:,:,i) = (mage(:,:,i)>0) - (mage(:,:,i)<0);
+    if display_on, 
+    ori_incr=pi/n_ori; % to convert jshi's coords to conventional image xy
+    ori_offset=ori_incr/2;
+    theta = ori_offset+([1:n_ori]-1)*ori_incr; % orientation detectors
+    % [gx,gy] are image gradient in image xy coords, winner take all
+    
+    ori = theta(max_id);
+    ori = ori .* (mago(:,:,i)>0) + (ori + pi).*(mago(:,:,i)<0);
+    gy{i} = mag_s(:,:,i) .* cos(ori);
+    gx{i} = -mag_s(:,:,i) .* sin(ori);
+    g = cat(3,gx{i},gy{i});
+    % phase map: edges are where the phase changes
+    mag_th = max(max(mag_s(:,:,i))) * threshold;
+    eg = (mag_s(:,:,i)>mag_th);
+    h = eg & [(mage(2:r,:,i) ~= mage(1:r-1,:,i)); zeros(1,nc)];
+    v = eg & [(mage(:,2:c,i) ~= mage(:,1:c-1,i)), zeros(nr,1)];
+    [y{i},x{i}] = find(h | v);
+    k = y{i} + (x{i}-1) * nr;
+    h = h(k);
+    v = v(k);
+    y{i} = y{i} + h * 0.5; % i
+    x{i} = x{i} + v * 0.5; % j
+    t = h + v * nr;
+    gx{i} = g(k) + g(k+t);
+    k = k + (nr * nc);
+    gy{i} = g(k) + g(k+t);
+    
+%     figure(1);
+%     clf;
+%     imagesc(I);colormap(gray);
+%     hold on;
+%     quiver(x,y,gx,gy); hold off;
+%     title(sprintf('scale = %d, press return',i));
+    
+   % pause;
+    0;
+else
+    x = [];
+    y = [];
+    gx = [];
+    gy =[];
+     g= [];
+ end
+end
diff --git a/SD-VBS/common/toolbox/MultiNcut/quickNcutHardBiais2.m b/SD-VBS/common/toolbox/MultiNcut/quickNcutHardBiais2.m
new file mode 100755
index 0000000..3ca1046
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/quickNcutHardBiais2.m
@@ -0,0 +1,187 @@
+% function [v,s] = quickNcutHardBiais(W,U,nbEigenValues,dataGraphCut)
+%ligne 35 : changement tim
+%[v,s] = ncut(W,nbEigenValues,[],offset); 
+%devient :
+%[v,s] = tim_ncut_rapide(W,nbEigenValues,[],offset);
+%et eigs devient tim_eigs
+% Input:
+%   W = affinity matrix
+%   U = hard constraint matrix, could be a cell of partial grouping
+%   nbEigenValues = number of eigenvectors
+%   offset = regularization factor, default = 0
+% Output:
+%   v = eigenvector
+%   s = eigenvalue of (W,d), s.t. U' * x = 0.
+% call eigs using my own * operation
+% Stella X. Yu, Jan 2002.
+function [v,s] = quickNcutHardBiais2(W,U,nbEigenValues,dataGraphCut)%voir : rajouter sigma
+n = size(W,1);
+nbEigenValues = min(nbEigenValues,n);
+offset = dataGraphCut.offset;
+%offset = 2;
+% degrees and regularization
+d = sum(abs(W),2);
+dr = 0.5 * (d - sum(W,2));
+d = d + offset * 2;
+dr = dr + offset;
+W = W + spdiags(dr,0,n,n);
+clear dr
+% normalize
+Dinvsqrt = 1./sqrt(d+eps);
+P = spmtimesd(W,Dinvsqrt,Dinvsqrt);
+clear W;
+% if max(max(P-P')) < 1e-10 %ou eps
+%      %S = sparse(1:n,1:n,0.5);
+%     P =max(P,P');
+%      % P=S*(P+P');
+%     %P=0.5*(P+P');
+%     options.issym = 1;
+% end
+P = sparsifyc(P,dataGraphCut.valeurMin);
+options.issym = 1;
+Ubar = spmtimesd(U,Dinvsqrt,[]);
+%Ubar = sparsifyc(Ubar,dataGraphCut.valeurMin); %voir
+options.disp = dataGraphCut.verbose;
+options.maxit = dataGraphCut.maxiterations;
+options.tol = dataGraphCut.eigsErrorTolerance;
+options.v0 = ones(size(P,1),1);%voir
+options.p = max(35,2*nbEigenValues); %voir
+options.p = min(options.p , n);
+% nouvelle idee : factorisation de Cholesky
+C=Ubar'*Ubar;
+%permutation = symamd(C);
+%R = cholinc(C(permutation,permutation));
+t_chol_Ubar = cputime;
+[R,ooo] = cholinc(C,'0');
+t_chol_Ubar = cputime - t_chol_Ubar;
+%if error occurs, check if Ubar = sparsifyc(Ubar,dataGraphCut.valeurMin);
+%sparsifies too much
+% compute  H = (Ubar'*Ubar)^(-1)
+% t_inv_H = cputime;
+% H = inv(sparsifyc(Ubar' * Ubar,dataGraphCut.valeurMin)); %changer
+% t_inv_H = cputime - t_inv_H;
+% H = sparsifyc(H,dataGraphCut.valeurMin);
+% tEigs = cputime;
+% if options.issym & max(max(H-H')) < 1e-10
+%     [vbar,s,convergence] = tim_eigs(@mex_projection_inv_symmetric,n,nbEigenValues,'lm',options,triu(P),Ubar,triu(H)); 
+% else
+%     [vbar,s,convergence] = tim_eigs(@mex_projection_inv,n,nbEigenValues,'lm',options,P,Ubar,H); 
+% end
+% tEigs = cputime - tEigs;
+% 
+R = sparsifyc(R,dataGraphCut.valeurMin);
+tEigs = cputime;
+if options.issym
+    [vbar,s,convergence] = tim_eigs(@mex_projection_QR_symmetric,n,nbEigenValues,'lm',options,tril(P),Ubar,R); 
+else
+    [vbar,s,convergence] = tim_eigs(@mex_projection_QR,n,nbEigenValues,'lm',options,P,Ubar,R); 
+end
+tEigs = cputime - tEigs;
+%afficheTexte(sprintf('\n\nTemps ecoule pendant eigs : %g',tEigs),dataGraphCut.verbose,2);
+%afficheTexte(sprintf('\nTemps ecoule pendant chol(Ubar''*Ubar) : %g',t_chol_Ubar),dataGraphCut.verbose);
+if convergence~=0
+    afficheTexte(sprintf('  (Non-convergence)'),dataGraphCut.verbose);
+end
+%disp(sprintf('nnz(P) : %f\n',nnz(P)));
+%disp(sprintf('nnz(Ubar) : %f\n',nnz(Ubar)));
+%disp(sprintf('nnz(R) : %f\n',nnz(R)));
+%disp(sprintf('nnz(global) : %f\n',nnz(P) + 4 * nnz(Ubar) + 4*nnz(R)));
+s = real(diag(s));
+[x,y] = sort(-s); 
+s = -x;
+vbar = vbar(:,y);
+v = spdiags(Dinvsqrt,0,n,n) * vbar;
+for  i=1:size(v,2)
+    %v(:,i) = v(:,i) / max(abs(v(:,i)));
+    v(:,i) = (v(:,i) / norm(v(:,i))  )*norm(ones(n,1));
+    if v(1,i)~=0
+        v(:,i) = - v(:,i) * sign(v(1,i));
+    end
+end
+% % nouvelle idee : factorisation de Cholesky
+% t_chol_Ubar = cputime;
+% R = chol(Ubar' * Ubar);
+% t_chol_Ubar = cputime - t_chol_Ubar;
+% R = sparsifyc(R,dataGraphCut.valeurMin);
+% tEigs = cputime;
+% if options.issym
+%     [vbar,s,convergence] = tim_eigs(@mex_projection_QR_symmetric,n,nbEigenValues,'lm',options,triu(P),Ubar,R); 
+% else
+%     [vbar,s,convergence] = tim_eigs(@mex_projection_QR,n,nbEigenValues,'lm',options,P,Ubar,R); 
+% end
+% tEigs = cputime - tEigs;
+% % compute H = (Ubar'*Ubar)^(-1)
+% t_inv_H = cputime;
+% H = inv(sparsifyc(Ubar' * Ubar,dataGraphCut.valeurMin)); %changer
+% t_inv_H = cputime - t_inv_H;
+% H = sparsifyc(H,dataGraphCut.valeurMin);
+% tEigs = cputime;
+% if options.issym & max(max(H-H')) < 1e-10
+%     [vbar,s,convergence] = tim_eigs(@mex_projection_inv_symmetric,n,nbEigenValues,'lm',options,triu(P),Ubar,triu(H)); 
+% else
+%     [vbar,s,convergence] = tim_eigs(@mex_projection_inv,n,nbEigenValues,'lm',options,P,Ubar,H); 
+% end
+% tEigs = cputime - tEigs;
+% % idee de mon rapport... semble pas marcher car R = qr(Ubar,0) est plus
+% % lent que H = inv(sparsifyc(Ubar' * Ubar,dataGraphCut.valeurMin));
+% t_qr_Ubar = cputime;
+% R = qr(Ubar,0);
+% t_qr_Ubar = cputime - t_qr_Ubar;
+% R = sparsifyc(R,dataGraphCut.valeurMin);
+% tEigs2 = cputime;
+% if options.issym
+%     [vbar2,s2,convergence] = tim_eigs(@mex_projection_QR_symmetric,n,nbEigenValues,'lm',options,triu(P),Ubar,R); 
+% else
+%     [vbar2,s2,convergence] = tim_eigs(@mex_projection_QR,n,nbEigenValues,'lm',options,P,Ubar,R); 
+% end
+% tEigs2 = cputime - tEigs2;
+% idee de Jianbo... semble pas marcher car on a besoin de prendre k maximal
+% dans [A,S,B] = svds(Ubar,k);
+% 
+% [A,S,B] = svds(Ubar,300);
+% A = sparsifyc(A,dataGraphCut.valeurMin); 
+% tEigs = cputime;
+% [vbar,s,convergence] = tim_eigs(@mex_projection_svd,n,nbEigenValues,'lm',options,P,A); 
+% afficheTexte(sprintf('\ninv(H) : %g',t_inv_H),dataGraphCut.verbose);
+% afficheTexte(sprintf('\n\nTemps ecoule pendant eigs : %g',tEigs2),dataGraphCut.verbose,2);
+% afficheTexte(sprintf('\nqr(Ubar) : %g',t_qr_Ubar),dataGraphCut.verbose);
+% disp(sprintf('nnz(H) : %f\n',nnz(H)));
+% disp(sprintf('nnz(global) : %f\n',nnz(P) + 4 * nnz(Ubar) + 2*nnz(H)));
diff --git a/SD-VBS/common/toolbox/MultiNcut/read_data.m b/SD-VBS/common/toolbox/MultiNcut/read_data.m
new file mode 100755
index 0000000..c0929c0
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/read_data.m
@@ -0,0 +1,13 @@
+%fnames = dir('/home/jshi/Results_DLIB/SegLabl*.mat');
+fnames = dir('/data/jshi/DLIB/Results/Results_DLIB/SegLabl*.mat');
+for j=1:length(fnames),
+        cm = sprintf('load /data/jshi/DLIB/Results/Results_DLIB/%s',fnames(j).name);
+        disp(cm);eval(cm);
+figure(1);imagesc(I); colormap(gray); axis image;
+figure(2); imagesc(SegLabel); axis image;
+pause;
+end
diff --git a/SD-VBS/common/toolbox/MultiNcut/readimage.m b/SD-VBS/common/toolbox/MultiNcut/readimage.m
new file mode 100755
index 0000000..e45fa19
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/readimage.m
@@ -0,0 +1,15 @@
+function I = readimage(fn,maxSize);
+Io = imread(fn);
+[nr,nc,nb] = size(Io);
+if nb>1,
+    I = rgb2gray(Io);
+else
+    I= Io;
+end
+%maxSize = 400;
+if max(nr,nc) > maxSize,
+    I = imresize(I,maxSize/max(nr,nc));
+end
diff --git a/SD-VBS/common/toolbox/MultiNcut/run_script.m b/SD-VBS/common/toolbox/MultiNcut/run_script.m
new file mode 100755
index 0000000..3d2e8c6
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/run_script.m
@@ -0,0 +1,60 @@
+%% set path for the MNcut code
+if 1,
+%   MNcutDir = '/home/jshi/Matlab/Toolbox/MultiNcut';
+  MNcutDir = 'C:\qihuizhu\Checkout\Human\Source\MultiNcut_new\MultiNcut';
+  path(path,MNcutDir);
+end
+%% set the image input and output dir.
+% imagedir = '/data/jshi/DLIB/image.cd';
+% imagedir = 'C:\qihuizhu\Checkout\Human\Source\Data\test';
+imagedir = 'C:\qihuizhu\Checkout\Human\Data\Current\baby_case5';
+% imageformat = 'ppm';
+imageformat = 'tif';
+% OutputDir = '/home/jshi/Results_DLIB';
+% OutputDir = 'C:\qihuizhu\Checkout\Human\Source\Data\test';
+OutputDir = 'C:\qihuizhu\Checkout\Human\Result\Segmentation\MultiNcut_new_03.07';
+a = dir(OutputDir);
+if (length(a) == 0), 
+  cm = sprintf('mkdir %s',OutputDir); 
+  disp(cm); eval(cm);
+end
+files = dir(sprintf('%s/*.%s',imagedir,imageformat));
+%% image size definition
+imageSize = 400;
+% for id =11:200,       
+for id = 1:length(files)
+   %for id = 19:19,
+    I=readimage(sprintf('%s/%s',imagedir,files(id).name),imageSize);
+     
+    num_segs = [10, 20];
+    
+    tic
+    [SegLabel,eigenVectors,eigenValues]= MNcut(I,num_segs);
+    toc
+    for j=1:size(SegLabel,3),
+      [gx,gy] = gradient(SegLabel(:,:,j));
+      bw = (abs(gx)>0.1) + (abs(gy) > 0.1);
+    
+      figure(1);clf; J1=showmask(double(I),bw); imagesc(J1);axis image;
+      cm = sprintf('print -djpeg %s/file%.4d-%.2d.jpg',OutputDir,id,num_segs(j)); disp(cm);eval(cm);
+    
+      figure(10);imagesc(SegLabel(:,:,j));axis image;
+      cm = sprintf('print -djpeg %s/Seg%.4d-%.2d.jpg',OutputDir,id,num_segs(j));disp(cm);eval(cm);
+% pause;
+    end
+    fname = files(id).name;
+    %cm = sprintf('save %s/SegLabl%.4d.mat I SegLabel fname',OutputDir,id); disp(cm);    eval(cm);
+    %cm = sprintf('save %s/SegEig%.4d.mat eigenVectors eigenValues',OutputDir,id);disp(cm); eval(cm);
+    
+end
diff --git a/SD-VBS/common/toolbox/MultiNcut/showmask.m b/SD-VBS/common/toolbox/MultiNcut/showmask.m
new file mode 100755
index 0000000..6fd1142
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/showmask.m
@@ -0,0 +1,65 @@
+% function RGB=showmask(V,M,hue);
+% Input:
+%    V = image
+%    M = nonnegative mask
+%    hue = a number in [0,1], red,yellow,green,cyan,...,red
+%          a char, 'r','g','b','y','c','m'
+%          or a matrix of the same size of image
+%          eg. hue = mask1 * 0.7 + mask2 * 1;
+%   
+% Output:
+%    RGB = an RGB image with V as shades and M as saturated hues
+%    If no output is required, this image is displayed.
+% Stella X. YU, 2000.  Based on Jianbo Shi's version.
+function RGB=showmask(V,M,hue);
+if nargin<3 | isempty(hue),
+   hue = 0;
+end
+if ischar(hue),
+    switch hue,
+        case 'r', hue = 1.0;
+        case 'g', hue = 0.3;
+        case 'b', hue = 0.7;
+        case 'y', hue = 0.15;
+        case 'c', hue = 0.55;
+        case 'm', hue = 0.85;
+    end
+end
+        
+V=V-min(V(:));
+V=V/max(V(:));
+V=.25+0.75*V; %brighten things up a bit
+M = double(M);
+M = M-min(M(:));
+M = M/max(M(:));
+H = hue+zeros(size(V));
+S = M;
+RGB = hsv2rgb(H,S,V);
+if nargout>0,   
+   return;
+end
+hold off; image(RGB); axis('image');
+s = cell(1,2);
+if isempty(inputname(1)),
+   s{1} = 'image';
+else
+   s{1} = inputname(1);
+end
+if isempty(inputname(2)),
+   s{2} = 'mask';
+else
+   s{2} = inputname(2);
+end
+title(sprintf('%s and colored %s',s{1},s{2}));
+if nargout==0,
+    clear RGB; 
+end
diff --git a/SD-VBS/common/toolbox/MultiNcut/sparsifyc.c b/SD-VBS/common/toolbox/MultiNcut/sparsifyc.c
new file mode 100755
index 0000000..82fca98
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/sparsifyc.c
@@ -0,0 +1,232 @@
+/*=================================================================
+* syntax: SPMX = SPARSIFY(MX, THRES)
+*
+* SPARSIFY  - sparsify the input matrix, i.e. ignore the values 
+*                       of the matrix which     are below a threshold
+*                       
+*                       Input:  - MX: m-by-n matrix (sparse or full)
+*                                       - THRES: threshold value (double)
+*
+*                       Output: - SPMX: m-by-n sparse matrix only with values 
+*                                       whose absolut value is above the given threshold
+*
+*                       Written by Mirko Visontai (10/24/2003)
+*=================================================================*/
+#include <math.h>
+#include "mex.h"
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+    /* Declare variable */
+    int i,m,n,nzmax,newnnz,col,processed,passed;
+        int starting_row_index, current_row_index, stopping_row_index;
+    double *in_pr,*in_pi,*out_pr,*out_pi;
+    int *in_ir,*in_jc,*out_ir,*out_jc;
+        double thres;
+    /* Check for proper number of input and output arguments */    
+    if ((nlhs != 1) || (nrhs != 2)){
+                mexErrMsgTxt("usage: SPMX = SPARSIFY(MX, THRES).");
+    } 
+        /* if matrix is complex threshold the norm of the numbers */
+        if (mxIsComplex(prhs[0])){
+                /* Check data type of input argument  */
+                if (mxIsSparse(prhs[0])){
+                        /* read input */
+                        in_pr   = mxGetPr(prhs[0]);
+                        in_pi   = mxGetPi(prhs[0]);
+                        in_ir   = mxGetIr(prhs[0]);
+                        in_jc   = mxGetJc(prhs[0]);
+                        nzmax   = mxGetNzmax(prhs[0]);
+                        m               = mxGetM(prhs[0]);
+                        n               = mxGetN(prhs[0]);
+                        thres   = mxGetScalar(prhs[1]);
+                        /* Count new nonzeros */
+                        newnnz=0;
+                        for(i=0; i<nzmax; i++){
+                                if (sqrt(in_pr[i]*in_pr[i] + in_pi[i]*in_pi[i])>thres) {newnnz++;}
+                        }
+                        if (newnnz>0){
+                                /* create output */
+                                plhs[0]         = mxCreateSparse(m,n,newnnz,mxCOMPLEX);
+                                if (plhs[0]==NULL)
+                                        mexErrMsgTxt("Could not allocate enough memory!\n");
+                                out_pr          = mxGetPr(plhs[0]);
+                                out_pi          = mxGetPr(plhs[0]);
+                                out_ir          = mxGetIr(plhs[0]);
+                                out_jc          = mxGetJc(plhs[0]);
+                                passed          = 0;
+                                out_jc[0]       = 0;
+                                for (col=0; col<n; col++){
+                                        starting_row_index = in_jc[col];
+                                        stopping_row_index = in_jc[col+1];
+                                        out_jc[col+1] = out_jc[col];
+                                        if (starting_row_index == stopping_row_index)
+                                                continue;
+                                        else {
+                                                for (current_row_index = starting_row_index;
+                                                        current_row_index < stopping_row_index;
+                                                        current_row_index++)  {
+                                                                if (sqrt(in_pr[current_row_index]*in_pr[current_row_index] + 
+                                                                             in_pi[current_row_index]*in_pi[current_row_index] ) > thres){
+                                                                        out_pr[passed]=in_pr[current_row_index];        
+                                                                        out_pi[passed]=in_pi[current_row_index];        
+                                                                        out_ir[passed]=in_ir[current_row_index];        
+                                                                        out_jc[col+1] = out_jc[col+1]+1;
+                                                                        passed++;
+                                                                }
+                                                }
+                                        }
+                                }
+                        }
+                        else{
+                                plhs[0] = mxCreateSparse(m,n,0,mxCOMPLEX);
+                        }
+                }
+                else{ /* for full matrices */
+                        /* read input */
+                        in_pr   = mxGetPr(prhs[0]);
+                        in_pi   = mxGetPr(prhs[0]);
+                        m               = mxGetM(prhs[0]);
+                        n               = mxGetN(prhs[0]);
+                        thres   = mxGetScalar(prhs[1]);
+                        /* Count new nonzeros */
+                        newnnz=0;
+                        for(i=0; i<m*n; i++){
+                                if (sqrt(in_pr[i]*in_pr[i] + in_pi[i]*in_pi[i])>thres) {newnnz++;}
+                        }
+                        if (newnnz>0){
+                                /* create output */
+                                plhs[0]         = mxCreateSparse(m,n,newnnz,mxCOMPLEX);
+                                if (plhs[0]==NULL)
+                                        mexErrMsgTxt("Could not allocate enough memory!\n");
+                                out_pr          = mxGetPr(plhs[0]);
+                                out_pi          = mxGetPi(plhs[0]);
+                                out_ir          = mxGetIr(plhs[0]);
+                                out_jc          = mxGetJc(plhs[0]);
+                                passed          = 0;
+                                out_jc[0]       = 0;
+                                for (col=0; col<n; col++){
+                                        out_jc[col+1] = out_jc[col];
+                                        for (current_row_index=0; current_row_index<m; current_row_index++){
+                                                        if (sqrt(in_pr[current_row_index+m*col]*in_pr[current_row_index+m*col] +
+                                                                         in_pi[current_row_index+m*col]*in_pi[current_row_index+m*col]) > thres){
+                                                                
+                                                                out_pr[passed]=in_pr[current_row_index+m*col];  
+                                                                out_ir[passed]=current_row_index;       
+                                                                out_jc[col+1] = out_jc[col+1]+1;
+                                                                passed++;
+                                                        }
+                                        }
+                                }
+                        }
+                        else{
+                                plhs[0] = mxCreateSparse(m,n,0,mxCOMPLEX);
+                        }
+                }
+        }
+        else { 
+        /* Check data type of input argument  */
+        if (mxIsSparse(prhs[0])){
+                        /* read input */
+                        in_pr   = mxGetPr(prhs[0]);
+                        in_ir   = mxGetIr(prhs[0]);
+                in_jc   = mxGetJc(prhs[0]);
+                        nzmax   = mxGetNzmax(prhs[0]);
+                        n               = mxGetN(prhs[0]);
+                        m               = mxGetM(prhs[0]);
+                        thres   = mxGetScalar(prhs[1]);
+                        /* Count new nonzeros */
+                        newnnz=0;
+                        for(i=0; i<nzmax; i++){
+                                if ((fabs(in_pr[i]))>thres) {newnnz++;}
+                        }
+                        if (newnnz>0){
+                                /* create output */
+                                plhs[0]         = mxCreateSparse(m,n,newnnz,mxREAL);
+                                if (plhs[0]==NULL)
+                                        mexErrMsgTxt("Could not allocate enough memory!\n");
+                                out_pr          = mxGetPr(plhs[0]);
+                        out_ir          = mxGetIr(plhs[0]);
+                        out_jc          = mxGetJc(plhs[0]);
+                                passed          = 0;
+                                out_jc[0]       = 0;
+                                for (col=0; col<n; col++){
+                                        starting_row_index = in_jc[col];
+                                        stopping_row_index = in_jc[col+1];
+                                        out_jc[col+1] = out_jc[col];
+                                        if (starting_row_index == stopping_row_index)
+                                                continue;
+                                        else {
+                                                for (current_row_index = starting_row_index;
+                                                        current_row_index < stopping_row_index;
+                                                        current_row_index++)  {
+                                                                if (fabs(in_pr[current_row_index])>thres){
+                                                                        out_pr[passed]=in_pr[current_row_index];        
+                                                                        out_ir[passed]=in_ir[current_row_index];        
+                                                                        out_jc[col+1] = out_jc[col+1]+1;
+                                                                        passed++;
+                                                                }
+                                                }
+                                        }
+                                }
+                        }
+                        else{
+                        plhs[0] = mxCreateSparse(m,n,0,mxREAL);
+                        }
+                }
+                else{ /* for full matrices */
+                        /* read input */
+                        in_pr   = mxGetPr(prhs[0]);
+                        n               = mxGetN(prhs[0]);
+                        m               = mxGetM(prhs[0]);
+                        thres   = mxGetScalar(prhs[1]);
+                        /* Count new nonzeros */
+                        newnnz=0;
+                        for(i=0; i<m*n; i++){
+                                if ((fabs(in_pr[i]))>thres) {newnnz++;}
+                        }
+                        if (newnnz>0){
+                                /* create output */
+                                plhs[0]         = mxCreateSparse(m,n,newnnz,mxREAL);
+                                if (plhs[0]==NULL)
+                                        mexErrMsgTxt("Could not allocate enough memory!\n");
+                                out_pr          = mxGetPr(plhs[0]);
+                        out_ir          = mxGetIr(plhs[0]);
+                        out_jc          = mxGetJc(plhs[0]);
+                                passed          = 0;
+                        out_jc[0]       = 0;
+                                for (col=0; col<n; col++){
+                                        out_jc[col+1] = out_jc[col];
+                                        for (current_row_index=0; current_row_index<m; current_row_index++){
+                                                        if (fabs(in_pr[current_row_index+m*col])>thres){
+                                                                out_pr[passed]=in_pr[current_row_index+m*col];  
+                                                                out_ir[passed]=current_row_index;       
+                                                                out_jc[col+1] = out_jc[col+1]+1;
+                                                                passed++;
+                                                        }
+                                        }
+                                }
+                        }
+                        else{
+                        plhs[0] = mxCreateSparse(m,n,0,mxREAL);
+                        }
+                }
+        }
+}
+ 
diff --git a/SD-VBS/common/toolbox/MultiNcut/sparsifyc.dll b/SD-VBS/common/toolbox/MultiNcut/sparsifyc.dll
new file mode 100755
index 0000000..cf832a6
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/sparsifyc.dll
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/sparsifyc.mexa64 b/SD-VBS/common/toolbox/MultiNcut/sparsifyc.mexa64
new file mode 100755
index 0000000..2f5ed26
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/sparsifyc.mexa64
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/sparsifyc.mexglx b/SD-VBS/common/toolbox/MultiNcut/sparsifyc.mexglx
new file mode 100755
index 0000000..0ba41b6
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/sparsifyc.mexglx
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/sparsifyc.mexmac b/SD-VBS/common/toolbox/MultiNcut/sparsifyc.mexmac
new file mode 100755
index 0000000..caa2306
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/sparsifyc.mexmac
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/spmtimesd.c b/SD-VBS/common/toolbox/MultiNcut/spmtimesd.c
new file mode 100755
index 0000000..a98dc0a
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/spmtimesd.c
@@ -0,0 +1,141 @@
+/*================================================================
+* spmtimesd.c
+* This routine computes a sparse matrix times a diagonal matrix
+* whose diagonal entries are stored in a full vector.
+*
+* Examples: 
+*     spmtimesd(m,d,[]) = diag(d) * m,
+*     spmtimesd(m,[],d) = m * diag(d)
+*     spmtimesd(m,d1,d2) = diag(d1) * m * diag(d2)
+*     m could be complex, but d is assumed real
+*
+* Stella X. Yu's first MEX function, Nov 9, 2001.
+% test sequence:
+m = 1000;
+n = 2000;
+a=sparse(rand(m,n));
+d1 = rand(m,1);
+d2 = rand(n,1);
+tic; b=spmtimesd(a,d1,d2); toc
+tic; bb = spdiags(d1,0,m,m) * a * spdiags(d2,0,n,n); toc
+e = (bb-b);
+max(abs(e(:)))
+*=================================================================*/
+# include "mex.h"
+void mexFunction(
+    int nargout,
+    mxArray *out[],
+    int nargin,
+    const mxArray *in[]
+)
+{
+    /* declare variables */
+    int i, j, k, m, n, nzmax, cmplx, xm, yn;
+    int *pir, *pjc, *qir, *qjc;
+    double *x, *y, *pr, *pi, *qr, *qi;
+    
+    /* check argument */
+    if (nargin != 3) {
+        mexErrMsgTxt("Three input arguments required");
+    }
+    if (nargout>1) {
+        mexErrMsgTxt("Too many output arguments.");
+    }
+    if (!(mxIsSparse(in[0]))) {
+        mexErrMsgTxt("Input argument #1 must be of type sparse");
+    }
+    if ( mxIsSparse(in[1]) || mxIsSparse(in[2]) ) {
+        mexErrMsgTxt("Input argument #2 & #3 must be of type full");
+    }
+  
+    /* computation starts */
+    m = mxGetM(in[0]);
+    n = mxGetN(in[0]);
+    pr = mxGetPr(in[0]);
+    pi = mxGetPi(in[0]);
+    pir = mxGetIr(in[0]);
+    pjc = mxGetJc(in[0]);
+ 
+    i = mxGetM(in[1]); 
+    j = mxGetN(in[1]);
+    xm = ((i>j)? i: j);
+    i = mxGetM(in[2]); 
+    j = mxGetN(in[2]);
+    yn = ((i>j)? i: j);
+   
+    if ( xm>0 && xm != m) {
+        mexErrMsgTxt("Row multiplication dimension mismatch.");
+    }
+    if ( yn>0 && yn != n) {
+        mexErrMsgTxt("Column multiplication dimension mismatch.");
+    }
+    
+ 
+    nzmax = mxGetNzmax(in[0]);
+    cmplx = (pi==NULL ? 0 : 1);    
+    out[0] = mxCreateSparse(m,n,nzmax,cmplx);
+    if (out[0]==NULL) {
+        mexErrMsgTxt("Not enough space for the output matrix.");
+    }    
+   
+    qr = mxGetPr(out[0]);
+    qi = mxGetPi(out[0]);
+    qir = mxGetIr(out[0]);
+    qjc = mxGetJc(out[0]);
+    
+    /* left multiplication */
+    x = mxGetPr(in[1]);
+    if (yn==0) {
+        for (j=0; j<n; j++) {
+            qjc[j] = pjc[j];
+            for (k=pjc[j]; k<pjc[j+1]; k++) {
+                i = pir[k];   
+                qir[k] = i;
+                 qr[k] = x[i] * pr[k];
+                 if (cmplx) {
+                     qi[k] = x[i] * pi[k];
+                 }
+            }
+        }
+        qjc[n] = k;
+        return;
+    }
+    
+    /* right multiplication */
+    y = mxGetPr(in[2]);
+    if (xm==0) {
+        for (j=0; j<n; j++) {
+            qjc[j] = pjc[j];
+            for (k=pjc[j]; k<pjc[j+1]; k++) {
+                qir[k] = pir[k];
+                qr[k]  = pr[k] * y[j];
+                if (cmplx) {
+                    qi[k] = qi[k] * y[j];
+                }
+            }
+       }
+       qjc[n] = k;
+       return;
+   }
+    
+   /* both sides */
+   for (j=0; j<n; j++) {
+       qjc[j] = pjc[j];
+       for (k=pjc[j]; k<pjc[j+1]; k++) {
+           i = pir[k];
+           qir[k]= i;
+           qr[k] = x[i] * pr[k] * y[j];
+           if (cmplx) {
+               qi[k] = x[i] * qi[k] * y[j];      
+           }
+       }
+       qjc[n] = k;
+   }
+   
+}   
diff --git a/SD-VBS/common/toolbox/MultiNcut/spmtimesd.dll b/SD-VBS/common/toolbox/MultiNcut/spmtimesd.dll
new file mode 100755
index 0000000..f78a650
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/spmtimesd.dll
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/spmtimesd.mexa64 b/SD-VBS/common/toolbox/MultiNcut/spmtimesd.mexa64
new file mode 100755
index 0000000..43c59e3
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/spmtimesd.mexa64
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/spmtimesd.mexglx b/SD-VBS/common/toolbox/MultiNcut/spmtimesd.mexglx
new file mode 100755
index 0000000..5bc9ec4
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/spmtimesd.mexglx
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/spmtimesd.mexmac b/SD-VBS/common/toolbox/MultiNcut/spmtimesd.mexmac
new file mode 100755
index 0000000..014113f
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/spmtimesd.mexmac
Binary files differ
diff --git a/SD-VBS/common/toolbox/MultiNcut/tim_eigs.m b/SD-VBS/common/toolbox/MultiNcut/tim_eigs.m
new file mode 100755
index 0000000..4b80262
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/tim_eigs.m
@@ -0,0 +1,1084 @@
+function  varargout = tim_eigs(varargin)
+nombre_A_times_X = 0; %tim
+nombreIterations = 0; %tim
+%seule difference avec eigs : 
+%       arguments_Afun = varargin{7-Amatrix-Bnotthere:end};
+%(pour l'instant : n'accepte que 2 arguments dans le cas de Afun : Afun(W,X))
+%permet d'aller plus vite en minimisant les acces a varargin 
+%(Timothee)
+%EIGS  Find a few eigenvalues and eigenvectors of a matrix using ARPACK.
+%   D = EIGS(A) returns a vector of A's 6 largest magnitude eigenvalues.
+%   A must be square and should be large and sparse.
+%
+%   [V,D] = EIGS(A) returns a diagonal matrix D of A's 6 largest magnitude
+%   eigenvalues and a matrix V whose columns are the corresponding eigenvectors.
+%
+%   [V,D,FLAG] = EIGS(A) also returns a convergence flag.  If FLAG is 0
+%   then all the eigenvalues converged; otherwise not all converged.
+%
+%   EIGS(A,B) solves the generalized eigenvalue problem A*V == B*V*D.  B must
+%   be symmetric (or Hermitian) positive definite and the same size as A.
+%   EIGS(A,[],...) indicates the standard eigenvalue problem A*V == V*D.
+%
+%   EIGS(A,K) and EIGS(A,B,K) return the K largest magnitude eigenvalues.
+%
+%   EIGS(A,K,SIGMA) and EIGS(A,B,K,SIGMA) return K eigenvalues based on SIGMA:
+%      'LM' or 'SM' - Largest or Smallest Magnitude
+%   For real symmetric problems, SIGMA may also be:
+%      'LA' or 'SA' - Largest or Smallest Algebraic
+%      'BE' - Both Ends, one more from high end if K is odd
+%   For nonsymmetric and complex problems, SIGMA may also be:
+%      'LR' or 'SR' - Largest or Smallest Real part
+%      'LI' or 'SI' - Largest or Smallest Imaginary part
+%   If SIGMA is a real or complex scalar including 0, EIGS finds the eigenvalues
+%   closest to SIGMA.  For scalar SIGMA, and also when SIGMA = 'SM' which uses
+%   the same algorithm as SIGMA = 0, B need only be symmetric (or Hermitian)
+%   positive semi-definite since it is not Cholesky factored as in the other cases.
+%
+%   EIGS(A,K,SIGMA,OPTS) and EIGS(A,B,K,SIGMA,OPTS) specify options:
+%   OPTS.issym: symmetry of A or A-SIGMA*B represented by AFUN [{0} | 1]
+%   OPTS.isreal: complexity of A or A-SIGMA*B represented by AFUN [0 | {1}]
+%   OPTS.tol: convergence: Ritz estimate residual <= tol*NORM(A) [scalar | {eps}]
+%   OPTS.maxit: maximum number of iterations [integer | {300}]
+%   OPTS.p: number of Lanczos vectors: K+1<p<=N [integer | {2K}]
+%   OPTS.v0: starting vector [N-by-1 vector | {randomly generated by ARPACK}]
+%   OPTS.disp: diagnostic information display level [0 | {1} | 2]
+%   OPTS.cholB: B is actually its Cholesky factor CHOL(B) [{0} | 1]
+%   OPTS.permB: sparse B is actually CHOL(B(permB,permB)) [permB | {1:N}]
+%
+%   EIGS(AFUN,N) accepts the function AFUN instead of the matrix A.
+%   Y = AFUN(X) should return
+%      A*X            if SIGMA is not specified, or is a string other than 'SM'
+%      A\X            if SIGMA is 0 or 'SM'
+%      (A-SIGMA*I)\X  if SIGMA is a nonzero scalar (standard eigenvalue problem)
+%      (A-SIGMA*B)\X  if SIGMA is a nonzero scalar (generalized eigenvalue problem)
+%   N is the size of A. The matrix A, A-SIGMA*I or A-SIGMA*B represented by AFUN is
+%   assumed to be real and nonsymmetric unless specified otherwise by OPTS.isreal
+%   and OPTS.issym. In all these EIGS syntaxes, EIGS(A,...) may be replaced by
+%   EIGS(AFUN,N,...).
+%
+%   EIGS(AFUN,N,K,SIGMA,OPTS,P1,P2,...) and EIGS(AFUN,N,B,K,SIGMA,OPTS,P1,P2,...)
+%   provide for additional arguments which are passed to AFUN(X,P1,P2,...).
+%
+%   Examples:
+%      A = delsq(numgrid('C',15));  d1 = eigs(A,5,'SM');
+%   Equivalently, if dnRk is the following one-line function:
+%      function y = dnRk(x,R,k)
+%      y = (delsq(numgrid(R,k))) \ x;
+%   then pass dnRk's additional arguments, 'C' and 15, to EIGS:
+%      n = size(A,1);  opts.issym = 1;  d2 = eigs(@dnRk,n,5,'SM',opts,'C',15);
+%
+%   See also EIG, SVDS, ARPACKC.
+%   Copyright 1984-2002 The MathWorks, Inc.
+%   $Revision: 1.45 $  $Date: 2002/05/14 18:50:58 $
+cputms = zeros(5,1);
+t0 = cputime; % start timing pre-processing
+if (nargout > 3)
+    error('Too many output arguments.')
+end
+% Process inputs and do error-checking
+if isa(varargin{1},'double')
+    A = varargin{1};
+    Amatrix = 1;
+else
+    A = fcnchk(varargin{1});
+    Amatrix = 0;
+end
+isrealprob = 1; % isrealprob = isreal(A) & isreal(B) & isreal(sigma)
+if Amatrix
+    isrealprob = isreal(A);
+end
+issymA = 0;
+if Amatrix
+    issymA = isequal(A,A');
+end
+if Amatrix
+    [m,n] = size(A);
+    if (m ~= n)
+        error('A must be a square matrix or a function which computes A*x.')
+    end
+else
+    n = varargin{2};
+    nstr = 'Size of problem, ''n'', must be a positive integer.';
+    if ~isequal(size(n),[1,1]) | ~isreal(n)
+        error(nstr)
+    end
+    if (round(n) ~= n)
+        warning('MATLAB:eigs:NonIntegerSize',['%s\n         ' ...
+            'Rounding input size.'],nstr)
+        n = round(n);
+    end
+    if issparse(n)
+        n = full(n);
+    end
+end
+Bnotthere = 0;
+Bstr = sprintf(['Generalized matrix B must be the same size as A and' ...
+        ' either a symmetric positive (semi-)definite matrix or' ...
+        ' its Cholesky factor.']);
+if (nargin < (3-Amatrix-Bnotthere))
+    B = [];
+    Bnotthere = 1;
+else
+    Bk = varargin{3-Amatrix-Bnotthere};
+    if isempty(Bk) % allow eigs(A,[],k,sigma,opts);
+        B = Bk;
+    else
+        if isequal(size(Bk),[1,1]) & (n ~= 1)
+            B = [];
+            k = Bk;
+            Bnotthere = 1;
+        else % eigs(9,8,...) assumes A=9, B=8, ... NOT A=9, k=8, ...
+            B = Bk;
+            if ~isa(B,'double') | ~isequal(size(B),[n,n])
+                error(Bstr)
+            end
+            isrealprob = isrealprob & isreal(B);
+        end
+    end
+end
+if Amatrix & ((nargin - ~Bnotthere)>4)
+    error('Too many inputs.')
+end
+if (nargin < (4-Amatrix-Bnotthere))
+    k = min(n,6);
+else
+    k = varargin{4-Amatrix-Bnotthere};
+end
+kstr = ['Number of eigenvalues requested, k, must be a' ...
+        ' positive integer <= n.'];
+if ~isa(k,'double') | ~isequal(size(k),[1,1]) | ~isreal(k) | (k>n)
+    error(kstr)
+end
+if issparse(k)
+    k = full(k);
+end
+if (round(k) ~= k)
+    warning('MATLAB:eigs:NonIntegerEigQty',['%s\n         ' ...
+            'Rounding number of eigenvalues.'],kstr)
+    k = round(k);
+end
+whchstr = sprintf(['Eigenvalue range sigma must be a valid 2-element string.']);
+if (nargin < (5-Amatrix-Bnotthere))
+    % default: eigs('LM') => ARPACK which='LM', sigma=0
+    eigs_sigma = 'LM';
+    whch = 'LM';
+    sigma = 0;
+else
+    eigs_sigma = varargin{5-Amatrix-Bnotthere};
+    if isstr(eigs_sigma)
+        % eigs(string) => ARPACK which=string, sigma=0
+        if ~isequal(size(eigs_sigma),[1,2])
+            whchstr = [whchstr sprintf(['\nFor real symmetric A, the choices are ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.'], ...
+                    'LM','SM','LA','SA','BE')];
+            whchstr = [whchstr sprintf(['\nFor non-symmetric or complex A, the choices are ''%s'', ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.\n'], ...
+                    'LM','SM','LR','SR','LI','SI')];
+            error(whchstr)
+        end
+        eigs_sigma = upper(eigs_sigma);
+        if isequal(eigs_sigma,'SM')
+            % eigs('SM') => ARPACK which='LM', sigma=0
+            whch = 'LM';
+        else
+            % eigs(string), where string~='SM' => ARPACK which=string, sigma=0
+            whch = eigs_sigma;
+        end
+        sigma = 0;
+    else
+        % eigs(scalar) => ARPACK which='LM', sigma=scalar
+        if ~isa(eigs_sigma,'double') | ~isequal(size(eigs_sigma),[1,1])
+            error('Eigenvalue shift sigma must be a scalar.')
+        end
+        sigma = eigs_sigma;
+        if issparse(sigma)
+            sigma = full(sigma);
+        end
+        isrealprob = isrealprob & isreal(sigma);
+        whch = 'LM';
+    end
+end
+tol = eps; % ARPACK's minimum tolerance is eps/2 (DLAMCH's EPS)
+maxit = [];
+p = [];
+info = int32(0); % use a random starting vector
+display = 1;
+cholB = 0;
+if (nargin >= (6-Amatrix-Bnotthere))
+    opts = varargin{6-Amatrix-Bnotthere};
+    if ~isa(opts,'struct')
+        error('Options argument must be a structure.')
+    end
+    
+    if isfield(opts,'issym') & ~Amatrix
+        issymA = opts.issym;
+        if (issymA ~= 0) & (issymA ~= 1)
+            error('opts.issym must be 0 or 1.')
+        end
+    end
+    
+    if isfield(opts,'isreal') & ~Amatrix
+        if (opts.isreal ~= 0) & (opts.isreal ~= 1)
+            error('opts.isreal must be 0 or 1.')
+        end
+        isrealprob = isrealprob & opts.isreal;
+    end
+    
+    if ~isempty(B) & (isfield(opts,'cholB') | isfield(opts,'permB'))
+        if isfield(opts,'cholB')
+            cholB = opts.cholB;
+            if (cholB ~= 0) & (cholB ~= 1)
+                error('opts.cholB must be 0 or 1.')
+            end
+            if isfield(opts,'permB')
+                if issparse(B) & cholB
+                    permB = opts.permB;
+                    if ~isequal(sort(permB),(1:n)) & ...
+                            ~isequal(sort(permB),(1:n)')
+                        error('opts.permB must be a permutation of 1:n.')
+                    end
+                else
+                    warning('MATLAB:eigs:IgnoredOptionPermB', ...
+                            ['Ignoring opts.permB since B is not its sparse' ...
+                            ' Cholesky factor.'])
+                end
+            else
+                permB = 1:n;
+            end
+        end
+    end
+    
+    if isfield(opts,'tol')
+        if ~isequal(size(opts.tol),[1,1]) | ~isreal(opts.tol) | (opts.tol<=0)
+            error(['Convergence tolerance opts.tol must be a strictly' ...
+                    ' positive real scalar.'])
+        else
+            tol = full(opts.tol);
+        end
+    end
+    
+    if isfield(opts,'p')
+        p = opts.p;
+        pstr = ['Number of basis vectors opts.p must be a positive' ...
+                ' integer <= n.'];
+        if ~isequal(size(p),[1,1]) | ~isreal(p) | (p<=0) | (p>n)
+            error(pstr)
+        end
+        if issparse(p)
+            p = full(p);
+        end
+        if (round(p) ~= p)
+            warning('MATLAB:eigs:NonIntegerVecQty',['%s\n         ' ...
+                    'Rounding number of basis vectors.'],pstr)
+            p = round(p);
+        end
+    end
+    
+    if isfield(opts,'maxit')
+        maxit = opts.maxit;
+        str = ['Maximum number of iterations opts.maxit must be' ...
+                ' a positive integer.'];
+        if ~isequal(size(maxit),[1,1]) | ~isreal(maxit) | (maxit<=0)
+            error(str)
+        end
+        if issparse(maxit)
+            maxit = full(maxit);
+        end
+        if (round(maxit) ~= maxit)
+            warning('MATLAB:eigs:NonIntegerIterationQty',['%s\n         ' ...
+                    'Rounding number of iterations.'],str)
+            maxit = round(maxit);
+        end
+    end
+    
+    if isfield(opts,'v0')
+        if ~isequal(size(opts.v0),[n,1])
+            error('Start vector opts.v0 must be n-by-1.')
+        end
+        if isrealprob
+            if ~isreal(opts.v0)
+                error(['Start vector opts.v0 must be real for real problems.'])
+            end
+            resid = full(opts.v0);
+        else
+            resid(1:2:(2*n-1),1) = full(real(opts.v0));
+            resid(2:2:2*n,1) = full(imag(opts.v0));
+        end
+        info = int32(1); % use resid as starting vector
+    end
+    
+    if isfield(opts,'disp')
+        display = opts.disp;
+        dispstr = 'Diagnostic level opts.disp must be an integer.';
+        if (~isequal(size(display),[1,1])) | (~isreal(display)) | (display<0)
+            error(dispstr)
+        end
+        if (round(display) ~= display)
+            warning('MATLAB:eigs:NonIntegerDiagnosticLevel', ...
+                    '%s\n         Rounding diagnostic level.',dispstr)
+            display = round(display);
+        end
+    end
+    
+    if isfield(opts,'cheb')
+        warning('MATLAB:eigs:ObsoleteOptionCheb', ...
+                ['Ignoring polynomial acceleration opts.cheb' ...
+                ' (no longer an option).']);
+    end
+    if isfield(opts,'stagtol')
+        warning('MATLAB:eigs:ObsoleteOptionStagtol', ...
+                ['Ignoring stagnation tolerance opts.stagtol' ...
+                ' (no longer an option).']);
+    end
+    
+end
+% Now we know issymA, isrealprob, cholB, and permB
+if isempty(p)
+    if isrealprob & ~issymA
+        p = min(max(2*k+1,20),n);
+    else
+        p = min(max(2*k,20),n);
+    end
+end
+if isempty(maxit)
+    maxit = max(300,ceil(2*n/max(p,1)));
+end
+if (info == int32(0))
+    if isrealprob
+        resid = zeros(n,1);
+    else
+        resid = zeros(2*n,1);
+    end
+end
+if ~isempty(B) % B must be symmetric (Hermitian) positive (semi-)definite
+    if cholB
+        if ~isequal(triu(B),B)
+            error(Bstr)
+        end
+    else
+        if ~isequal(B,B')
+            error(Bstr)
+        end
+    end
+end
+useeig = 0;
+if isrealprob & issymA
+    knstr = sprintf(['For real symmetric problems, must have number' ...
+            ' of eigenvalues k < n.\n']);
+else
+    knstr = sprintf(['For nonsymmetric and complex problems, must have' ...
+            ' number of eigenvalues k < n-1.\n']);
+end
+if isempty(B)
+    knstr = [knstr sprintf(['Using eig(full(A)) instead.'])];
+else
+    knstr = [knstr sprintf(['Using eig(full(A),full(B)) instead.'])];
+end
+if (k == 0)
+    useeig = 1;
+end
+if isrealprob & issymA
+    if (k > n-1)
+        if (n >= 6)
+            warning('MATLAB:eigs:TooManyRequestedEigsForRealSym', ...
+                '%s',knstr)
+        end
+        useeig = 1;
+    end
+else
+    if (k > n-2)
+        if (n >= 7)
+            warning('MATLAB:eigs:TooManyRequestedEigsForComplexNonsym', ...
+                '%s',knstr)
+        end
+        useeig = 1;
+    end
+end
+if isrealprob & issymA
+    if ~isreal(sigma)
+        error(['For real symmetric problems, eigenvalue shift sigma must' ...
+                ' be real.'])
+    end
+else
+    if ~isrealprob & issymA & ~isreal(sigma)
+        warning('MATLAB:eigs:ComplexShiftForRealProblem', ...
+            ['Complex eigenvalue shift sigma on a Hermitian problem' ...
+            ' (all real eigenvalues).'])
+    end
+end
+if isrealprob & issymA
+    if strcmp(whch,'LR')
+        whch = 'LA';
+        warning('MATLAB:eigs:SigmaChangedToLA', ... 
+            ['For real symmetric problems, sigma value ''LR''' ...
+            ' (Largest Real) is now ''LA'' (Largest Algebraic).'])
+    end
+    if strcmp(whch,'SR')
+        whch = 'SA';
+        warning('MATLAB:eigs:SigmaChangedToSA', ...
+            ['For real symmetric problems, sigma value ''SR''' ...
+            ' (Smallest Real) is now ''SA'' (Smallest Algebraic).'])
+    end
+    if ~ismember(whch,{'LM', 'SM', 'LA', 'SA', 'BE'})
+        whchstr = [whchstr sprintf(['\nFor real symmetric A, the choices are ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.'], ...
+                'LM','SM','LA','SA','BE')];
+        error(whchstr)
+    end
+else
+    if strcmp(whch,'BE')
+        warning('MATLAB:eigs:SigmaChangedToLM', ...
+            ['Sigma value ''BE'' is now only available for real' ...
+            ' symmetric problems.  Computing ''LM'' eigenvalues instead.'])
+        whch = 'LM';
+    end
+    if ~ismember(whch,{'LM', 'SM', 'LR', 'SR', 'LI', 'SI'})
+         whchstr = [whchstr sprintf(['\nFor non-symmetric or complex A, the choices are ''%s'', ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.\n'], ...
+                 'LM','SM','LR','SR','LI','SI')];
+        error(whchstr)
+    end
+end
+% Now have enough information to do early return on cases eigs does not handle
+if useeig
+    if (nargout <= 1)
+        varargout{1} = eigs2(A,n,B,k,whch,sigma,cholB, ...
+            varargin{7-Amatrix-Bnotthere:end});
+    else
+        [varargout{1},varargout{2}] = eigs2(A,n,B,k,whch,sigma,cholB, ...
+            varargin{7-Amatrix-Bnotthere:end});
+    end
+    if (nargout == 3)
+        varargout{3} = 0; % flag indicates "convergence"
+    end
+    return
+end
+if isrealprob & ~issymA
+    sigmar = real(sigma);
+    sigmai = imag(sigma);
+end
+if isrealprob & issymA
+    if (p <= k)
+        error(['For real symmetric problems, must have number of' ...
+                ' basis vectors opts.p > k.'])
+    end
+else
+    if (p <= k+1)
+        error(['For nonsymmetric and complex problems, must have number of' ...
+                ' basis vectors opts.p > k+1.'])
+    end
+end
+if isequal(whch,'LM') & ~isequal(eigs_sigma,'LM')
+    % A*x = lambda*M*x, M symmetric (positive) semi-definite
+    % => OP = inv(A - sigma*M)*M and B = M
+    % => shift-and-invert mode
+    mode = 3;
+elseif isempty(B)
+    % A*x = lambda*x
+    % => OP = A and B = I
+    mode = 1;
+else % B is not empty
+    % Do not use mode=2.
+    % Use mode = 1 with OP = R'\(A*(R\x)) and B = I
+    % where R is B's upper triangular Cholesky factor: B = R'*R.
+    % Finally, V = R\V returns the actual generalized eigenvectors of A and B.
+    mode = 1;
+end
+if cholB
+    pB = 0;
+    RB = B;
+    RBT = B';
+end
+if (mode == 3) & Amatrix % need lu(A-sigma*B)
+    if issparse(A) & (isempty(B) | issparse(B))
+        if isempty(B)
+            AsB = A - sigma * speye(n);
+        else
+            if cholB
+                AsB = A - sigma * RBT * RB;
+            else
+                AsB = A - sigma * B;
+            end
+        end
+        [L,U,P,Q] = lu(AsB);
+        [perm,dummy] = find(Q);
+    else
+        if isempty(B)
+            AsB = A - sigma * eye(n);
+        else
+            if cholB
+                AsB = A - sigma * RBT * RB;
+            else
+                AsB = A - sigma * B;
+            end
+        end
+        [L,U,P] = lu(AsB);
+    end
+    dU = diag(U);
+    rcondestU = full(min(abs(dU)) / max(abs(dU)));
+    if (rcondestU < eps)
+        if isempty(B)
+            ds = sprintf(['(A-sigma*I) has small reciprocal condition' ...
+                    ' estimate: %f\n'],rcondestU);
+        else
+            ds = sprintf(['(A-sigma*B) has small reciprocal condition' ...
+                    ' estimate: %f\n'],rcondestU);
+        end
+        ds = [ds sprintf(['indicating that sigma is near an exact' ...
+                    ' eigenvalue. The\nalgorithm may not converge unless' ...
+                    ' you try a new value for sigma.\n'])];
+        warning('MATLAB:eigs:SigmaNearExactEig',ds)
+    end
+end
+if (mode == 1) & ~isempty(B) & ~cholB % need chol(B)
+    if issparse(B)
+        permB = symamd(B);
+        [RB,pB] = chol(B(permB,permB));
+    else
+        [RB,pB] = chol(B);
+    end
+    if (pB == 0)
+        RBT = RB';
+    else
+        error(Bstr)
+    end
+end
+% Allocate outputs and ARPACK work variables
+if isrealprob
+    if issymA % real and symmetric
+        prefix = 'ds';
+        v = zeros(n,p);
+        ldv = int32(size(v,1));
+        ipntr = int32(zeros(15,1));
+        workd = zeros(n,3);
+        lworkl = p*(p+8);
+        workl = zeros(lworkl,1);
+        lworkl = int32(lworkl);
+        d = zeros(k,1);
+    else % real but not symmetric
+        prefix = 'dn';
+        v = zeros(n,p);
+        ldv = int32(size(v,1));
+        ipntr = int32(zeros(15,1));
+        workd = zeros(n,3);
+        lworkl = 3*p*(p+2);
+        workl = zeros(lworkl,1);
+        lworkl = int32(lworkl);
+        workev = zeros(3*p,1);
+        d = zeros(k+1,1);
+        di = zeros(k+1,1);
+    end
+else % complex
+    prefix = 'zn';
+    zv = zeros(2*n*p,1);
+    ldv = int32(n);
+    ipntr = int32(zeros(15,1));
+    workd = complex(zeros(n,3));
+    zworkd = zeros(2*prod(size(workd)),1);
+    lworkl = 3*p^2+5*p;
+    workl = zeros(2*lworkl,1);
+    lworkl = int32(lworkl);
+    workev = zeros(2*2*p,1);
+    zd = zeros(2*(k+1),1);
+    rwork = zeros(p,1);
+end
+ido = int32(0); % reverse communication parameter
+if isempty(B) | (~isempty(B) & (mode == 1))
+    bmat = 'I'; % standard eigenvalue problem
+else
+    bmat = 'G'; % generalized eigenvalue problem
+end
+nev = int32(k); % number of eigenvalues requested
+ncv = int32(p); % number of Lanczos vectors
+iparam = int32(zeros(11,1));
+iparam([1 3 7]) = int32([1 maxit mode]);
+select = int32(zeros(p,1));
+cputms(1) = cputime - t0; % end timing pre-processing
+iter = 0;
+ariter = 0;
+%tim
+indexArgumentsAfun = 7-Amatrix-Bnotthere:length(varargin);
+nbArgumentsAfun = length(indexArgumentsAfun);
+if nbArgumentsAfun >=1
+    arguments_Afun = varargin{7-Amatrix-Bnotthere};
+end
+if nbArgumentsAfun >=2
+    arguments_Afun2 = varargin{7-Amatrix-Bnotthere+1};
+end
+if nbArgumentsAfun >=3
+    arguments_Afun3 = varargin{7-Amatrix-Bnotthere+2};
+end
+%fin tim
+while (ido ~= 99)
+        
+    t0 = cputime; % start timing ARPACK calls **aupd
+        
+    if isrealprob
+        arpackc( [prefix 'aupd'], ido, ...
+            bmat, int32(n), whch, nev, tol, resid, ncv, ...
+            v, ldv, iparam, ipntr, workd, workl, lworkl, info);
+    else
+        zworkd(1:2:end-1) = real(workd);
+        zworkd(2:2:end) = imag(workd);
+        arpackc( 'znaupd', ido, ...
+            bmat, int32(n), whch, nev, tol, resid, ncv, zv, ...
+            ldv, iparam, ipntr, zworkd, workl, lworkl, ...
+            rwork, info );
+        workd = reshape(complex(zworkd(1:2:end-1),zworkd(2:2:end)),[n,3]);
+    end
+    if (info < 0)
+        es = sprintf('Error with ARPACK routine %saupd: info = %d',...
+           prefix,full(info));
+        error(es)
+    end
+     
+    cputms(2) = cputms(2) + (cputime-t0); % end timing ARPACK calls **aupd
+    t0 = cputime; % start timing MATLAB OP(X)
+    
+    % Compute which columns of workd ipntr references
+    
+    
+    
+    
+    
+    %[row,col1] = ind2sub([n,3],double(ipntr(1)));
+    %tim
+    row = mod(double(ipntr(1))-1,n) + 1 ;
+    col1 = floor((double(ipntr(1))-1)/n) + 1;
+    
+    
+    if (row ~= 1)
+        str = sprintf(['ipntr(1)=%d does not refer to the start of a' ...
+                ' column of the %d-by-3 array workd.'],full(ipntr(1)),n);
+        error(str)
+    end
+    
+    
+    
+    %[row,col2] = ind2sub([n,3],double(ipntr(2)));
+    %tim
+    row = mod(double(ipntr(2))-1,n) + 1 ;
+    col2 = floor((double(ipntr(2))-1)/n) + 1;
+    
+    
+    if (row ~= 1)
+        str = sprintf(['ipntr(2)=%d does not refer to the start of a' ...
+                ' column of the %d-by-3 array workd.'],full(ipntr(2)),n);
+        error(str)
+    end
+    if ~isempty(B) & (mode == 3) & (ido == 1)
+        [row,col3] = ind2sub([n,3],double(ipntr(3)));
+        if (row ~= 1)
+            str = sprintf(['ipntr(3)=%d does not refer to the start of a' ...
+                    ' column of the %d-by-3 array workd.'],full(ipntr(3)),n);
+            error(str)
+        end
+    end
+    switch (ido)
+        case {-1,1}
+        if Amatrix
+            if (mode == 1)
+                if isempty(B)
+                    % OP = A*x
+                    workd(:,col2) = A * workd(:,col1);
+                else
+                    % OP = R'\(A*(R\x))
+                    if issparse(B)
+                        workd(permB,col2) = RB \ workd(:,col1);
+                        workd(:,col2) = A * workd(:,col2);
+                        workd(:,col2) = RBT \ workd(permB,col2);
+                    else
+                        workd(:,col2) = RBT \ (A * (RB \ workd(:,col1)));
+                    end
+                end
+            elseif (mode == 3)
+                if isempty(B)
+                    if issparse(A)
+                        workd(perm,col2) = U \ (L \ (P * workd(:,col1)));
+                    else
+                        workd(:,col2) = U \ (L \ (P * workd(:,col1)));
+                    end
+                else % B is not empty
+                    if (ido == -1)
+                        if cholB
+                            workd(:,col2) = RBT * (RB * workd(:,col1));
+                        else
+                            workd(:,col2) = B * workd(:,col1);
+                        end
+                        if issparse(A) & issparse(B)
+                            workd(perm,col2) = U \ (L \ (P * workd(:,col1)));
+                        else
+                            workd(:,col2) = U \ (L \ (P * workd(:,col1)));
+                        end
+                    elseif (ido == 1)
+                        if issparse(A) & issparse(B)
+                            workd(perm,col2) = U \ (L \ (P * workd(:,col3)));
+                        else
+                            workd(:,col2) = U \ (L \ (P * workd(:,col3)));
+                        end
+                    end
+                end
+            else % mode is not 1 or 3
+                error(['Unknown mode returned from ' prefix 'aupd.'])
+            end
+        else % A is not a matrix
+            if (mode == 1)
+                if isempty(B)
+                    % OP = A*x
+                    %workd(:,col2) = feval(A,workd(:,col1),varargin{7-Amatrix-Bnotthere:end});                    
+                    
+                    
+                    
+                    
+                    
+                    nombre_A_times_X = nombre_A_times_X + 1;
+                    
+                    
+                    
+                    pause(0); %voir
+                    
+                    if nbArgumentsAfun == 1
+                        workd(:,col2) = feval(A,workd(:,col1),arguments_Afun);
+                        %workd(:,col2) = max(workd(:,col2),0);
+                    elseif nbArgumentsAfun == 2                        
+                        workd(:,col2) = feval(A,workd(:,col1),arguments_Afun,arguments_Afun2);
+                    elseif nbArgumentsAfun == 3                        
+                        workd(:,col2) = feval(A,workd(:,col1),arguments_Afun,arguments_Afun2,arguments_Afun3);
+                    else
+                        workd(:,col2) = feval(A,workd(:,col1),varargin{indexArgumentsAfun});                    
+                    end                      
+                    %workd(:,col2) = tim_w_times_x_c(workd(:,col1),arguments_Afun); %slower
+                   
+                else
+                    % OP = R'\(A*(R\x))
+                    if issparse(B)
+                        workd(permB,col2) = RB \ workd(:,col1);
+                        workd(:,col2) = feval(A,workd(:,col2),arguments_Afun);
+                        workd(:,col2) = RBT \ workd(permB,col2);
+                    else
+                        workd(:,col2) = RBT \ feval(A,(RB\workd(:,col1)),arguments_Afun);
+                    end
+                end
+            elseif (mode == 3)
+                if isempty(B)
+                    workd(:,col2) = feval(A,workd(:,col1), arguments_Afun);
+                else
+                    if (ido == -1)
+                        if cholB
+                            workd(:,col2) = RBT * (RB * workd(:,col1));
+                        else
+                            workd(:,col2) = B * workd(:,col1);
+                        end
+                        workd(:,col2) = feval(A,workd(:,col2), arguments_Afun);
+                    elseif (ido == 1)
+                        workd(:,col2) = feval(A,workd(:,col3), arguments_Afun);
+                    end
+                end
+            else % mode is not 1 or 3
+                error(['Unknown mode returned from ' prefix 'aupd.'])
+            end
+        end % if Amatrix
+    case 2
+        if (mode == 3)
+            if cholB
+                workd(:,col2) = RBT * (RB * workd(:,col1));
+            else
+                workd(:,col2) = B * workd(:,col1);
+            end
+        else
+            error(['Unknown mode returned from ' prefix 'aupd.'])
+        end
+    case 3
+        % setting iparam(1) = ishift = 1 ensures this never happens
+        warning('MATLAB:eigs:WorklShiftsUnsupported', ...
+            ['EIGS does not yet support computing the shifts in workl.' ...
+            ' Setting reverse communication parameter to 99 and returning.'])
+        ido = int32(99);
+    case 99
+    otherwise
+        error(['Unknown value of reverse communication parameter' ...
+                ' returned from ' prefix 'aupd.'])      
+    end
+    
+    cputms(3) = cputms(3) + (cputime-t0); % end timing MATLAB OP(X)
+    
+    %tim
+    if nombreIterations ~= double(ipntr(15))
+        nombreIterations = double(ipntr(15));
+    end
+    if display >= 1 && display <=2 
+        iter = double(ipntr(15));
+        if (iter > ariter) & (ido ~= 99)
+            ariter = iter;
+            ds = sprintf(['Iteration %d: a few Ritz values of the' ...
+                    ' %d-by-%d matrix:'],iter,p,p);
+            disp(ds)
+            if isrealprob
+                if issymA
+                    dispvec = [workl(double(ipntr(6))+(0:p-1))];
+                    if isequal(whch,'BE')
+                        % roughly k Large eigenvalues and k Small eigenvalues
+                        disp(dispvec(max(end-2*k+1,1):end))
+                    else
+                        % k eigenvalues
+                        disp(dispvec(max(end-k+1,1):end))
+                    end
+                else
+                    dispvec = [complex(workl(double(ipntr(6))+(0:p-1)), ...
+                            workl(double(ipntr(7))+(0:p-1)))];
+                    % k+1 eigenvalues (keep complex conjugate pairs together)
+                    disp(dispvec(max(end-k,1):end))
+                end
+            else
+                dispvec = [complex(workl(2*double(ipntr(6))-1+(0:2:2*(p-1))), ...
+                        workl(2*double(ipntr(6))+(0:2:2*(p-1))))];
+                disp(dispvec(max(end-k+1,1):end))
+            end
+        end
+    end
+    
+end % while (ido ~= 99)
+t0 = cputime; % start timing post-processing
+flag = 0;
+if (info < 0)
+    es = sprintf('Error with ARPACK routine %saupd: info = %d',prefix,full(info));
+    error(es)
+else
+    if (nargout >= 2)
+        rvec = int32(1); % compute eigenvectors
+    else
+        rvec = int32(0); % do not compute eigenvectors
+    end
+    
+    if isrealprob
+        if issymA
+            arpackc( 'dseupd', rvec, 'A', select, ...
+                d, v, ldv, sigma, ...
+                bmat, int32(n), whch, nev, tol, resid, ncv, ...
+                v, ldv, iparam, ipntr, workd, workl, lworkl, info );
+            if isequal(whch,'LM') | isequal(whch,'LA')
+                d = flipud(d);
+                if (rvec == 1)
+                    v(:,1:k) = v(:,k:-1:1);
+                end
+            end
+            if ((isequal(whch,'SM') | isequal(whch,'SA')) & (rvec == 0))
+                d = flipud(d);
+            end
+        else
+            arpackc( 'dneupd', rvec, 'A', select, ...
+                d, di, v, ldv, sigmar, sigmai, workev, ...
+                bmat, int32(n), whch, nev, tol, resid, ncv, ...
+                v, ldv, iparam, ipntr, workd, workl, lworkl, info );
+            d = complex(d,di);
+            if rvec
+                d(k+1) = [];
+            else
+                zind = find(d == 0);
+                if isempty(zind)
+                    d = d(k+1:-1:2);
+                else
+                    d(max(zind)) = [];
+                    d = flipud(d);
+                end
+            end
+        end
+    else
+        zsigma = [real(sigma); imag(sigma)];
+        arpackc( 'zneupd', rvec, 'A', select, ...
+            zd, zv, ldv, zsigma, workev, ...
+            bmat, int32(n), whch, nev, tol, resid, ncv, zv, ...
+            ldv, iparam, ipntr, zworkd, workl, lworkl, ...
+            rwork, info );
+        if issymA
+            d = zd(1:2:end-1);
+        else
+            d = complex(zd(1:2:end-1),zd(2:2:end));
+        end
+        v = reshape(complex(zv(1:2:end-1),zv(2:2:end)),[n p]);
+    end
+    
+    if (info ~= 0)
+        es = ['Error with ARPACK routine ' prefix 'eupd: '];
+        switch double(info)
+            case 2
+            ss = sum(select);
+            if (ss < k)
+            es = [es ...
+                    '  The logical variable select was only set with ' int2str(ss) ...
+                    ' 1''s instead of nconv=' int2str(double(iparam(5))) ...
+                    ' (k=' int2str(k) ').' ...
+                    ' Please contact the ARPACK authors at arpack@caam.rice.edu.'];
+            else
+            es = [es ...
+                    'The LAPACK reordering routine ' prefix(1) ...
+                    'trsen did not return all ' int2str(k) ' eigenvalues.']
+            end
+            case 1
+            es = [es ...
+                    'The Schur form could not be reordered by the LAPACK routine ' ...
+                    prefix(1) 'trsen.' ...
+                    ' Please contact the ARPACK authors at arpack@caam.rice.edu.'];
+            case -14
+            es = [es prefix ...
+                    'aupd did not find any eigenvalues to sufficient accuracy.'];
+            otherwise
+            es = [es sprintf('info = %d',full(info))];
+        end
+        error(es)
+    else
+        nconv = double(iparam(5));
+        if (nconv == 0)
+            if (nargout < 3)
+                warning('MATLAB:eigs:NoEigsConverged', ...
+                    'None of the %d requested eigenvalues converged.',k)
+            else
+                flag = 1;
+            end
+        elseif (nconv < k)
+            if (nargout < 3)
+                warning('MATLAB:eigs:NotAllEigsConverged', ...
+                    'Only %d of the %d requested eigenvalues converged.',nconv,k)
+            else
+                flag = 1;
+            end
+        end
+    end % if (info ~= 0)
+end % if (info < 0)
+if (issymA) | (~isrealprob)
+    if (nargout <= 1)
+        if isrealprob
+            varargout{1} = d;
+        else
+            varargout{1} = d(k:-1:1,1);
+        end
+    else
+        varargout{1} = v(:,1:k);
+        varargout{2} = diag(d(1:k,1));
+        if (nargout >= 3)
+            varargout{3} = flag;
+        end
+    end
+else
+    if (nargout <= 1)
+        varargout{1} = d;
+    else
+        cplxd = find(di ~= 0);
+        % complex conjugate pairs of eigenvalues occur together
+        cplxd = cplxd(1:2:end);
+        v(:,[cplxd cplxd+1]) = [complex(v(:,cplxd),v(:,cplxd+1)) ...
+                complex(v(:,cplxd),-v(:,cplxd+1))];
+        varargout{1} = v(:,1:k);
+        varargout{2} = diag(d);
+        if (nargout >= 3)
+            varargout{3} = flag;
+        end
+    end
+end
+if (nargout >= 2) & (mode == 1) & ~isempty(B)
+    if issparse(B)
+        varargout{1}(permB,:) = RB \ varargout{1};
+    else
+        varargout{1} = RB \ varargout{1};
+    end
+end
+cputms(4) = cputime-t0; % end timing post-processing
+cputms(5) = sum(cputms(1:4)); % total time
+if (display >= 2) %tim
+    if (mode == 1)
+        innerstr = sprintf(['Compute A*X:' ...
+                '                               %f\n'],cputms(3));
+    elseif (mode == 2)
+        innerstr = sprintf(['Compute A*X and solve B*X=Y for X:' ...
+                '         %f\n'],cputms(3));
+    elseif (mode == 3)
+        if isempty(B)
+            innerstr = sprintf(['Solve (A-SIGMA*I)*X=Y for X:' ...
+                    '               %f\n'],cputms(3));
+        else
+            innerstr = sprintf(['Solve (A-SIGMA*B)*X=B*Y for X:' ...
+                    '             %f\n'],cputms(3));
+        end
+    end
+    if ((mode == 3) & (Amatrix))
+        if isempty(B)
+            prepstr = sprintf(['Pre-processing, including lu(A-sigma*I):' ...
+                    '   %f\n'],cputms(1));
+        else
+            prepstr = sprintf(['Pre-processing, including lu(A-sigma*B):' ...
+                    '   %f\n'],cputms(1));
+        end
+    elseif ((mode == 2) & (~cholB))
+        prepstr = sprintf(['Pre-processing, including chol(B):' ...
+                '         %f\n'],cputms(1));
+    else
+        prepstr = sprintf(['Pre-processing:' ...
+                '                            %f\n'],cputms(1));
+    end
+    sstr = sprintf(['***********CPU Timing Results in seconds***********']);
+    ds = sprintf(['\n' sstr '\n' ...
+            prepstr ...
+            'ARPACK''s %saupd:                           %f\n' ...
+            innerstr ...
+            'Post-processing with ARPACK''s %seupd:      %f\n' ...
+            '***************************************************\n' ...
+            'Total:                                     %f\n' ...
+            sstr '\n'], ...
+        prefix,cputms(2),prefix,cputms(4),cputms(5));
+    disp(ds)
+    if nombre_A_times_X > 0      %tim
+        disp(sprintf('Number of iterations : %f\n',nombreIterations));
+        disp(sprintf('Number of times A*X was computed : %f\n',nombre_A_times_X));
+        disp(sprintf('Average time for A*X : %f\n',cputms(3)/nombre_A_times_X));          
+    end
+end