1 files changed, 158 insertions, 0 deletions

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