THESE/codes/graphe/Ncut_9/affinityic.cpp

   1 /*================================================================\r
   2 * function w = affinityic(emag,ephase,pi,pj,sigma)\r
   3 * Input:\r
   4 *   emag = edge strength at each pixel\r
   5 *   ephase = edge phase at each pixel\r
   6 *   [pi,pj] = index pair representation for MALTAB sparse matrices\r
   7 *   sigma = sigma for IC energy\r
   8 * Output:\r
   9 *   w = affinity with IC at [pi,pj]\r
  10 *\r
  11 \r
  12 % test sequence\r
  13 f = synimg(10);\r
  14 [i,j] = cimgnbmap(size(f),2);\r
  15 [ex,ey,egx,egy] = quadedgep(f);\r
  16 a = affinityic(ex,ey,egx,egy,i,j)\r
  17 show_dist_w(f,a);\r
  18 \r
  19 * Stella X. Yu, Nov 19, 2001.\r
  20 *=================================================================*/\r
  21 \r
  22 # include "mex.h"\r
  23 # include "math.h"\r
  24 \r
  25 void mexFunction(\r
  26     int nargout,\r
  27     mxArray *out[],\r
  28     int nargin,\r
  29     const mxArray *in[]\r
  30 )\r
  31 {\r
  32     /* declare variables */\r
  33     int nr, nc, np, total;\r
  34     int i, j, k, ix, iy, jx, jy, ii, jj, iip1, jjp1, iip2, jjp2, step;\r
  35     double sigma, di, dj, a, z, maxori, phase1, phase2, slope;\r
  36 //      int *ir, *jc;\r
  37     mwIndex*ir,*jc;\r
  38         unsigned int *pi, *pj;\r
  39         double *emag, *ephase, *w;\r
  40     \r
  41     /* check argument */\r
  42     if (nargin<4) {\r
  43         mexErrMsgTxt("Four input arguments required");\r
  44     }\r
  45     if (nargout>1) {\r
  46         mexErrMsgTxt("Too many output arguments");\r
  47     }\r
  48 \r
  49     /* get edgel information */\r
  50         nr = mxGetM(in[0]);\r
  51         nc = mxGetN(in[0]);\r
  52         if ( nr*nc ==0 || nr != mxGetM(in[1]) || nc != mxGetN(in[1]) ) {\r
  53             mexErrMsgTxt("Edge magnitude and phase shall be of the same image size");\r
  54         }\r
  55     emag = mxGetPr(in[0]);\r
  56     ephase = mxGetPr(in[1]);\r
  57     np = nr * nc;\r
  58     \r
  59     /* get new index pair */\r
  60     if (!mxIsUint32(in[2]) | !mxIsUint32(in[3])) {\r
  61         mexErrMsgTxt("Index pair shall be of type UINT32");\r
  62     }\r
  63     if (mxGetM(in[3]) * mxGetN(in[3]) != np + 1) {\r
  64         mexErrMsgTxt("Wrong index representation");\r
  65     }\r
  66     pi = (unsigned int*)mxGetData(in[2]);\r
  67     pj = (unsigned int*)mxGetData(in[3]);    \r
  68 \r
  69     /* create output */\r
  70     out[0] = mxCreateSparse(np,np,pj[np],mxREAL);\r
  71     if (out[0]==NULL) {\r
  72             mexErrMsgTxt("Not enough memory for the output matrix");\r
  73         }\r
  74         w = mxGetPr(out[0]);\r
  75         ir = mxGetIr(out[0]);\r
  76         jc = mxGetJc(out[0]);\r
  77         \r
  78     /* find my sigma */\r
  79         if (nargin<5) {\r
  80             sigma = 0;\r
  81         for (k=0; k<np; k++) { \r
  82             if (emag[k]>sigma) { sigma = emag[k]; }\r
  83         }\r
  84         sigma = sigma / 6;\r
  85         printf("sigma = %6.5f",sigma);\r
  86         } else {\r
  87             sigma = mxGetScalar(in[4]);\r
  88         }\r
  89         a = 0.5 / (sigma * sigma);\r
  90         \r
  91     /* computation */ \r
  92     total = 0;\r
  93     for (j=0; j<np; j++) {            \r
  94 \r
  95         jc[j] = total;\r
  96         jx = j / nr; /* col */\r
  97         jy = j % nr; /* row */\r
  98         \r
  99         for (k=pj[j]; k<pj[j+1]; k++) {\r
 100         \r
 101             i = pi[k];\r
 102           \r
 103             if (i==j) {\r
 104                 maxori = 1;\r
 105             \r
 106             } else {\r
 107         \r
 108                 ix = i / nr; \r
 109                 iy = i % nr;\r
 110 \r
 111                 /* scan */            \r
 112                 di = (double) (iy - jy);\r
 113                 dj = (double) (ix - jx);\r
 114             \r
 115                 maxori = 0.;\r
 116                     phase1 = ephase[j];\r
 117 \r
 118                        \r
 119                 /* sample in i direction */\r
 120                 if (abs(di) >= abs(dj)) {  \r
 121                     slope = dj / di;\r
 122                     step = (iy>=jy) ? 1 : -1;\r
 123                 \r
 124                     iip1 = jy;\r
 125                     jjp1 = jx;\r
 126         \r
 127                        \r
 128                         for (ii=0;ii<abs(di);ii++){\r
 129                             iip2 = iip1 + step;\r
 130                             jjp2 = (int)(0.5 + slope*(iip2-jy) + jx);\r
 131                   \r
 132                             phase2 = ephase[iip2+jjp2*nr];\r
 133                \r
 134                             if (phase1 != phase2) {\r
 135                                 z = (emag[iip1+jjp1*nr] + emag[iip2+jjp2*nr]);\r
 136                                 if (z > maxori){\r
 137                                     maxori = z;\r
 138                                 }\r
 139                             } \r
 140                      \r
 141                             iip1 = iip2;\r
 142                             jjp1 = jjp2;\r
 143                             phase1 = phase2;\r
 144                         }\r
 145                     \r
 146                     /* sample in j direction */    \r
 147                 } else { \r
 148                         slope = di / dj;\r
 149                         step =  (ix>=jx) ? 1: -1;\r
 150 \r
 151                     jjp1 = jx;\r
 152                         iip1 = jy;                 \r
 153             \r
 154          \r
 155                         for (jj=0;jj<abs(dj);jj++){\r
 156                             jjp2 = jjp1 + step;\r
 157                             iip2 = (int)(0.5+ slope*(jjp2-jx) + jy);\r
 158                   \r
 159                             phase2 = ephase[iip2+jjp2*nr];\r
 160              \r
 161                             if (phase1 != phase2){\r
 162                                 z = (emag[iip1+jjp1*nr] + emag[iip2+jjp2*nr]);\r
 163                                 if (z > maxori){ \r
 164                                     maxori = z; \r
 165                                 }\r
 166                                 \r
 167                             }\r
 168           \r
 169                             iip1 = iip2;\r
 170                             jjp1 = jjp2;\r
 171                             phase1 = phase2;\r
 172                         }\r
 173                 }            \r
 174             \r
 175                 maxori = 0.5 * maxori;\r
 176                 maxori = exp(-maxori * maxori * a);\r
 177             }       \r
 178                     ir[total] = i;\r
 179 \r
 180                     w[total] = maxori;\r
 181                     total = total + 1;\r
 182                         \r
 183                 } /* i */\r
 184     } /* j */\r
 185         \r
 186     jc[np] = total;\r
 187 }  \r