THESE/codes/graphe/Ncut_9/cimgnbmap.cpp

   1 /*================================================================\r
   2 * function [i,j] = cimgnbmap([nr,nc], nb_r, sample_rate)\r
   3 *   computes the neighbourhood index matrix of an image,\r
   4 *   with each neighbourhood sampled.\r
   5 * Input:\r
   6 *   [nr,nc] = image size\r
   7 *   nb_r = neighbourhood radius, could be [r_i,r_j] for i,j\r
   8 *   sample_rate = sampling rate, default = 1\r
   9 * Output:\r
  10 *   [i,j] = each is a column vector, give indices of neighbour pairs\r
  11 *     UINT32 type\r
  12 *       i is of total length of valid elements, 0 for first row\r
  13 *       j is of length nr * nc + 1\r
  14 *\r
  15 * See also: imgnbmap.c, id2cind.m\r
  16 *\r
  17 * Examples: \r
  18 *   [i,j] = imgnbmap(10, 20); % [10,10] are assumed\r
  19 *\r
  20 * Stella X. Yu, Nov 12, 2001.\r
  21 \r
  22 % test sequence:\r
  23 nr = 15;\r
  24 nc = 15;\r
  25 nbr = 1;\r
  26 [i,j] = cimgnbmap([nr,nc], nbr);\r
  27 mask = csparse(i,j,ones(length(i),1),nr*nc);\r
  28 show_dist_w(rand(nr,nc),mask)\r
  29 \r
  30 *=================================================================*/\r
  31 \r
  32 # include "mex.h"\r
  33 # include "math.h"\r
  34 # include <time.h>\r
  35 \r
  36 void mexFunction(\r
  37     int nargout,\r
  38     mxArray *out[],\r
  39     int nargin,\r
  40     const mxArray *in[]\r
  41 )\r
  42 {\r
  43     /* declare variables */\r
  44     int nr, nc, np, nb, total;\r
  45         double *dim, sample_rate;\r
  46     int r_i, r_j, a1, a2, b1, b2, self, neighbor;\r
  47     int i, j, k, s, t, nsamp, th_rand, no_sample;\r
  48     unsigned long *p;\r
  49     \r
  50     \r
  51     /* check argument */\r
  52     if (nargin < 2) {\r
  53         mexErrMsgTxt("Two input arguments required");\r
  54     }\r
  55     if (nargout> 2) {\r
  56         mexErrMsgTxt("Too many output arguments.");\r
  57     }\r
  58 \r
  59     /* get image size */\r
  60     i = mxGetM(in[0]);\r
  61     j = mxGetN(in[0]);\r
  62     dim = (double *)mxGetData(in[0]);\r
  63     nr = (int)dim[0];\r
  64     if (j>1 || i>1) {\r
  65         nc = (int)dim[1];\r
  66     } else {\r
  67         nc = nr;\r
  68     }\r
  69     np = nr * nc;\r
  70     \r
  71     /* get neighbourhood size */\r
  72     i = mxGetM(in[1]);\r
  73     j = mxGetN(in[1]);\r
  74     dim = (double*)mxGetData(in[1]);\r
  75     r_i = (int)dim[0];\r
  76     if (j>1 || i>1) {\r
  77         r_j = (int)dim[1];              \r
  78     } else {\r
  79         r_j = r_i;\r
  80     }\r
  81     if (r_i<0) { r_i = 0; }\r
  82         if (r_j<0) { r_j = 0; }\r
  83 \r
  84         /* get sample rate */\r
  85         if (nargin==3) {                \r
  86                 sample_rate = (mxGetM(in[2])==0) ? 1: mxGetScalar(in[2]);\r
  87     } else {\r
  88                 sample_rate = 1;\r
  89     }\r
  90         /* prepare for random number generator */\r
  91         if (sample_rate<1) {\r
  92         srand( (unsigned)time( NULL ) );\r
  93         th_rand = (int)ceil((double)RAND_MAX * sample_rate);\r
  94         no_sample = 0;\r
  95     } else {\r
  96                 sample_rate = 1;\r
  97         th_rand = RAND_MAX;\r
  98         no_sample = 1;\r
  99     }\r
 100     \r
 101         /* figure out neighbourhood size */\r
 102 \r
 103     nb = (r_i + r_i + 1) * (r_j + r_j + 1); \r
 104     if (nb>np) {\r
 105         nb = np;\r
 106     }\r
 107     nb = (int)ceil((double)nb * sample_rate);    \r
 108 \r
 109         /* intermediate data structure */\r
 110         p = (unsigned long *)mxCalloc(np * (nb+1), sizeof(unsigned long));\r
 111         if (p==NULL) {\r
 112         mexErrMsgTxt("Not enough space for my computation.");\r
 113         }\r
 114         \r
 115     /* computation */    \r
 116         total = 0;\r
 117     for (j=0; j<nc; j++) {\r
 118     for (i=0; i<nr; i++) {\r
 119 \r
 120                 self = i + j * nr;\r
 121 \r
 122                 /* put self in, otherwise the index is not ordered */\r
 123                 p[self] = p[self] + 1;\r
 124                 p[self+p[self]*np] = self;\r
 125 \r
 126         /* j range */\r
 127                 b1 = j;\r
 128         b2 = j + r_j;\r
 129         if (b2>=nc) { b2 = nc-1; }                \r
 130     \r
 131                 /* i range */\r
 132         a1 = i - r_i;\r
 133                 if (a1<0) { a1 = 0; }\r
 134         a2 = i + r_i;\r
 135         if (a2>=nr) { a2 = nr-1; }\r
 136        \r
 137                 /* number of more samples needed */\r
 138                 nsamp = nb - p[self];\r
 139 \r
 140                 k = 0;          \r
 141                 t = b1;\r
 142                 s = i + 1;\r
 143                 if (s>a2) {\r
 144                         s = a1;\r
 145                         t = t + 1;\r
 146                 }\r
 147                 while (k<nsamp && t<=b2) {\r
 148                         if (no_sample || (rand()<th_rand)) {\r
 149                                 k = k + 1;\r
 150                                 neighbor = s + t * nr;\r
 151                                 \r
 152                                 p[self] = p[self] + 1;                                  \r
 153                                 p[self+p[self]*np] = neighbor;\r
 154                         \r
 155                                 p[neighbor] = p[neighbor] + 1;\r
 156                                 p[neighbor+p[neighbor]*np] = self;\r
 157                         }\r
 158                         s = s + 1;\r
 159                         if (s>a2) {\r
 160                 s = a1;\r
 161                                 t = t + 1;\r
 162                         }\r
 163                 } /* k */\r
 164 \r
 165                 total = total + p[self];\r
 166         } /* i */\r
 167     } /* j */\r
 168     \r
 169     /* i, j */\r
 170     out[0] = mxCreateNumericMatrix(total, 1, mxUINT32_CLASS, mxREAL);\r
 171         out[1] = mxCreateNumericMatrix(np+1,  1, mxUINT32_CLASS, mxREAL);\r
 172     unsigned int *qi = (unsigned int *)mxGetData(out[0]);\r
 173         unsigned int *qj = (unsigned int *)mxGetData(out[1]);\r
 174         if (out[0]==NULL || out[1]==NULL) {\r
 175             mexErrMsgTxt("Not enough space for the output matrix.");\r
 176         }\r
 177 \r
 178         total = 0;\r
 179     for (j=0; j<np; j++) {\r
 180                 qj[j] = total;\r
 181                 s = j + np;\r
 182                 for (t=0; t<p[j]; t++) {\r
 183                     qi[total] = p[s];\r
 184                         total = total + 1;\r
 185                         s = s + np;\r
 186                 }\r
 187     }\r
 188         qj[np] = total;\r
 189 \r
 190         mxFree(p);\r
 191 }  \r