THESE/codes/snake/gvfc/gvfc.c

   1 /* NOTE:\r
   2  * \r
   3  * You must have "mex" command working in your matlab.  In matlab, \r
   4  * type "mex gvfc.c" to compile the code. See usage for details of \r
   5  * calling this function.\r
   6  *\r
   7  * Replace GVF(...) with GVFC in the snakedeform.m. The speed is \r
   8  * significantly faster than the Matlab version.\r
   9  */\r
  10 \r
  11 /***************************************************************************\r
  12 Copyright (c) 1996-1999 Chenyang Xu and Jerry Prince.\r
  13 \r
  14 This software is copyrighted by Chenyang Xu and Jerry Prince.  The\r
  15 following terms apply to all files associated with the software unless\r
  16 explicitly disclaimed in individual files.\r
  17 \r
  18 The authors hereby grant permission to use this software and its\r
  19 documentation for any purpose, provided that existing copyright\r
  20 notices are retained in all copies and that this notice is included\r
  21 verbatim in any distributions. Additionally, the authors grant\r
  22 permission to modify this software and its documentation for any\r
  23 purpose, provided that such modifications are not distributed without\r
  24 the explicit consent of the authors and that existing copyright\r
  25 notices are retained in all copies.\r
  26 \r
  27 IN NO EVENT SHALL THE AUTHORS OR DISTRIBUTORS BE LIABLE TO ANY PARTY FOR\r
  28 DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT\r
  29 OF THE USE OF THIS SOFTWARE, ITS DOCUMENTATION, OR ANY DERIVATIVES THEREOF,\r
  30 EVEN IF THE AUTHORS HAVE BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.\r
  31 \r
  32 THE AUTHORS AND DISTRIBUTORS SPECIFICALLY DISCLAIM ANY WARRANTIES, INCLUDING,\r
  33 BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A\r
  34 PARTICULAR PURPOSE, AND NON-INFRINGEMENT.  THIS SOFTWARE IS PROVIDED ON AN\r
  35 "AS IS" BASIS, AND THE AUTHORS AND DISTRIBUTORS HAVE NO OBLIGATION TO PROVIDE\r
  36 MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.\r
  37 ****************************************************************************/\r
  38 \r
  39 #include "mex.h"\r
  40 #include <stdio.h>\r
  41 #include <stdlib.h>\r
  42 #include <math.h>\r
  43 #include <string.h>\r
  44 \r
  45 /* \r
  46  * [u,v] = GVFC(f, mu, ITER)\r
  47  * \r
  48  * assuming dx = dy = 1, and dt = 1, \r
  49  * need to figure out CFL\r
  50  *\r
  51  * Gradient Vector Flow\r
  52  *\r
  53  * Chenyang Xu and Jerry L. Prince \r
  54  * Copyright (c) 1996-99 by Chenyang Xu and Jerry L. Prince\r
  55  * Image Analysis and Communications Lab, Johns Hopkins University\r
  56  *\r
  57  */\r
  58 \r
  59 /* PG Macros */\r
  60 #define PG_MAX(a,b)  (a>b? a: b)\r
  61 #define PG_MIN(a,b)  (a<b? a: b)\r
  62 \r
  63 typedef unsigned char  PGbyte;\r
  64 typedef          char  PGchar;\r
  65 typedef unsigned short PGushort;\r
  66 typedef          short PGshort; \r
  67 typedef unsigned int   PGuint;\r
  68 typedef          int   PGint;\r
  69 typedef float          PGfloat;\r
  70 typedef double         PGdouble;\r
  71 typedef void           PGvoid;\r
  72 \r
  73 /* Prints error message: modified to work in mex environment */\r
  74 PGvoid pgError(char error_text[])\r
  75 {\r
  76    char buf[255];\r
  77    \r
  78    sprintf(buf, "GVFC.MEX: %s", error_text);\r
  79    mexErrMsgTxt(buf);\r
  80    \r
  81    return;\r
  82 }\r
  83 \r
  84 \r
  85 /* Allocates a matrix of doubles with range [nrl..nrh][ncl..nch] */\r
  86 PGdouble **pgDmatrix(nrl,nrh,ncl,nch)\r
  87      int nrl,nrh,ncl,nch;\r
  88 {\r
  89    int j;\r
  90    long bufsize,bufptr;\r
  91    PGdouble **m;\r
  92    \r
  93    bufsize = (nrh-nrl+1)*sizeof(PGdouble*)\r
  94       + (nrh-nrl+1)*(nch-ncl+1)*sizeof(PGdouble);\r
  95    \r
  96    m=(PGdouble **) malloc(bufsize);\r
  97    if (!m) pgError("allocation failure 1 in pgDmatrix()");\r
  98    m -= nrl;\r
  99    \r
 100    bufptr = ((long) (m+nrl)) + (nrh-nrl+1)*sizeof(PGdouble*);\r
 101    for(j=nrl;j<=nrh;j++) {\r
 102       m[j] = ((PGdouble *) (bufptr+(j-nrl)*(nch-ncl+1)*sizeof(PGdouble)));\r
 103       m[j] -= ncl;\r
 104    }\r
 105    \r
 106    return m;\r
 107 }\r
 108 \r
 109 \r
 110 /* Frees a matrix allocated by dmatrix */\r
 111 PGvoid pgFreeDmatrix(m,nrl,nrh,ncl,nch)\r
 112      PGdouble **m;\r
 113      int nrl,nrh,ncl,nch;\r
 114 {\r
 115    free((char*) (m+nrl));\r
 116    return;\r
 117 }\r
 118 \r
 119 \r
 120 void GVFC(int YN, int XN, double *f, double *ou, double *ov, \r
 121           double mu, int ITER)\r
 122 {\r
 123    double mag2, temp, tempx, tempy, fmax, fmin;\r
 124    int count, x, y, XN_1, XN_2, YN_1, YN_2;\r
 125    \r
 126    PGdouble **fx, **fy, **u, **v, **Lu, **Lv, **g, **c1, **c2, **b;\r
 127    \r
 128    /* define constants and create row-major double arrays */\r
 129    XN_1 = XN - 1;\r
 130    XN_2 = XN - 2;\r
 131    YN_1 = YN - 1;\r
 132    YN_2 = YN - 2;\r
 133    fx = pgDmatrix(0,YN_1,0,XN_1);\r
 134    fy = pgDmatrix(0,YN_1,0,XN_1);\r
 135    u  = pgDmatrix(0,YN_1,0,XN_1);\r
 136    v  = pgDmatrix(0,YN_1,0,XN_1);\r
 137    Lu = pgDmatrix(0,YN_1,0,XN_1);\r
 138    Lv = pgDmatrix(0,YN_1,0,XN_1);\r
 139    g  = pgDmatrix(0,YN_1,0,XN_1);\r
 140    c1 = pgDmatrix(0,YN_1,0,XN_1);\r
 141    c2 = pgDmatrix(0,YN_1,0,XN_1);\r
 142    b  = pgDmatrix(0,YN_1,0,XN_1);\r
 143    \r
 144 \r
 145    /************** I: Normalize the edge map to [0,1] **************/\r
 146    fmax = -1e10;\r
 147    fmin = 1e10;\r
 148    for (x=0; x<=YN*XN-1; x++) {\r
 149       fmax = PG_MAX(fmax,f[x]);\r
 150       fmin = PG_MIN(fmin,f[x]);\r
 151    }\r
 152 \r
 153    if (fmax==fmin) \r
 154       mexErrMsgTxt("Edge map is a constant image.");      \r
 155 \r
 156    for (x=0; x<=YN*XN-1; x++) \r
 157       f[x] = (f[x]-fmin)/(fmax-fmin);\r
 158 \r
 159    /**************** II: Compute edge map gradient *****************/\r
 160    /* I.1: Neumann boundary condition: \r
 161     *      zero normal derivative at boundary \r
 162     */\r
 163    /* Deal with corners */\r
 164    fx[0][0] = fy[0][0] = fx[0][XN_1] = fy[0][XN_1] = 0;\r
 165    fx[YN_1][XN_1] = fy[YN_1][XN_1] = fx[YN_1][0] = fy[YN_1][0] = 0;\r
 166 \r
 167    /* Deal with left and right column */\r
 168    for (y=1; y<=YN_2; y++) {\r
 169       fx[y][0] = fx[y][XN_1] = 0;\r
 170       fy[y][0] = 0.5 * (f[y+1] - f[y-1]);\r
 171       fy[y][XN_1] = 0.5 * (f[y+1 + XN_1*YN] - f[y-1 + XN_1*YN]);\r
 172    }\r
 173 \r
 174    /* Deal with top and bottom row */\r
 175    for (x=1; x<=XN_2; x++) {\r
 176       fy[0][x] = fy[YN_1][x] = 0;\r
 177       fx[0][x] = 0.5 * (f[(x+1)*YN] - f[(x-1)*YN]);\r
 178       fx[YN_1][x] = 0.5 * (f[YN_1 + (x+1)*YN] - f[YN_1 + (x-1)*YN]);\r
 179    }\r
 180    \r
 181    /* I.2: Compute interior derivative using central difference */\r
 182    for (y=1; y<=YN_2; y++)\r
 183       for (x=1; x<=XN_2; x++) {\r
 184          /* NOTE: f is stored in column major */\r
 185          fx[y][x] = 0.5 * (f[y + (x+1)*YN] - f[y + (x-1)*YN]);   \r
 186          fy[y][x] = 0.5 * (f[y+1 + x*YN] - f[y-1 + x*YN]);\r
 187       }\r
 188    \r
 189    /******* III: Compute parameters and initializing arrays **********/\r
 190    temp = -1.0/(mu*mu);\r
 191    for (y=0; y<=YN_1; y++)\r
 192       for (x=0; x<=XN_1; x++) {\r
 193          tempx = fx[y][x];\r
 194          tempy = fy[y][x];\r
 195          /* initial GVF vector */\r
 196          u[y][x] = tempx;\r
 197          v[y][x] = tempy;\r
 198          /* gradient magnitude square */\r
 199          mag2 = tempx*tempx + tempy*tempy; \r
 200          \r
 201          g[y][x] = mu;\r
 202          b[y][x] = mag2;\r
 203 \r
 204          c1[y][x] = b[y][x] * tempx;\r
 205          c2[y][x] = b[y][x] * tempy;\r
 206       }\r
 207    \r
 208    /* free memory of fx and fy */\r
 209    pgFreeDmatrix(fx,0,YN_1,0,XN_1);\r
 210    pgFreeDmatrix(fy,0,YN_1,0,XN_1);\r
 211 \r
 212    /************* Solve GVF = (u,v) iteratively ***************/\r
 213    for (count=1; count<=ITER; count++) {\r
 214       /* IV: Compute Laplace operator using Neuman condition */\r
 215       /* IV.1: Deal with corners */\r
 216       Lu[0][0] = (u[0][1] + u[1][0])*0.5 - u[0][0]; \r
 217       Lv[0][0] = (v[0][1] + v[1][0])*0.5 - v[0][0];\r
 218       Lu[0][XN_1] = (u[0][XN_2] + u[1][XN_1])*0.5 - u[0][XN_1];\r
 219       Lv[0][XN_1] = (v[0][XN_2] + v[1][XN_1])*0.5 - v[0][XN_1];\r
 220       Lu[YN_1][0] = (u[YN_1][1] + u[YN_2][0])*0.5 - u[YN_1][0];\r
 221       Lv[YN_1][0] = (v[YN_1][1] + v[YN_2][0])*0.5 - v[YN_1][0];\r
 222       Lu[YN_1][XN_1] = (u[YN_1][XN_2] + u[YN_2][XN_1])*0.5 - u[YN_1][XN_1];\r
 223       Lv[YN_1][XN_1] = (v[YN_1][XN_2] + v[YN_2][XN_1])*0.5 - v[YN_1][XN_1];\r
 224       \r
 225       /* IV.2: Deal with left and right columns */\r
 226       for (y=1; y<=YN_2; y++) { \r
 227          Lu[y][0] = (2*u[y][1] + u[y-1][0] + u[y+1][0])*0.25 - u[y][0];\r
 228          Lv[y][0] = (2*v[y][1] + v[y-1][0] + v[y+1][0])*0.25 - v[y][0];\r
 229          Lu[y][XN_1] = (2*u[y][XN_2] + u[y-1][XN_1] \r
 230                         + u[y+1][XN_1])*0.25 - u[y][XN_1];\r
 231          Lv[y][XN_1] = (2*v[y][XN_2] + v[y-1][XN_1] \r
 232                         + v[y+1][XN_1])*0.25 - v[y][XN_1];\r
 233       }\r
 234       \r
 235       /* IV.3: Deal with top and bottom rows */\r
 236       for (x=1; x<=XN_2; x++) { \r
 237          Lu[0][x] = (2*u[1][x] + u[0][x-1] + u[0][x+1])*0.25 - u[0][x];\r
 238          Lv[0][x] = (2*v[1][x] + v[0][x-1] + v[0][x+1])*0.25 - v[0][x];\r
 239          Lu[YN_1][x] = (2*u[YN_2][x] + u[YN_1][x-1] \r
 240                         + u[YN_1][x+1])*0.25 - u[YN_1][x];\r
 241          Lv[YN_1][x] = (2*v[YN_2][x] + v[YN_1][x-1] \r
 242                         + v[YN_1][x+1])*0.25 - v[YN_1][x];\r
 243       }\r
 244       \r
 245       /* IV.4: Compute interior */\r
 246       for (y=1; y<=YN_2; y++)\r
 247          for (x=1; x<=XN_2; x++) {\r
 248             Lu[y][x] = (u[y][x-1] + u[y][x+1] \r
 249                         + u[y-1][x] + u[y+1][x])*0.25 - u[y][x];\r
 250             Lv[y][x] = (v[y][x-1] + v[y][x+1]\r
 251                         + v[y-1][x] + v[y+1][x])*0.25 - v[y][x];\r
 252          }\r
 253       \r
 254       /******** V: Update GVF ************/\r
 255       for (y=0; y<=YN_1; y++)\r
 256          for (x=0; x<=XN_1; x++) {\r
 257             u[y][x] = (1- b[y][x])*u[y][x] + g[y][x]*Lu[y][x]*4 + c1[y][x];\r
 258             v[y][x] = (1- b[y][x])*v[y][x] + g[y][x]*Lv[y][x]*4 + c2[y][x];\r
 259          }\r
 260       \r
 261       /* print iteration number */\r
 262       printf("%5d",count);\r
 263       if (count%15 == 0)\r
 264          printf("\n");\r
 265    }\r
 266    printf("\n");\r
 267    \r
 268    /* copy u,v to the output in column major order */\r
 269    for (y=0; y<=YN_1; y++)\r
 270       for (x=0; x<=XN_1; x++) {\r
 271          ou[x*YN+y] = u[y][x];\r
 272          ov[x*YN+y] = v[y][x];\r
 273       }\r
 274 \r
 275    /* free all the array memory */\r
 276    pgFreeDmatrix(u,0,YN_1,0,XN_1);\r
 277    pgFreeDmatrix(v,0,YN_1,0,XN_1);\r
 278    pgFreeDmatrix(Lu,0,YN_1,0,XN_1);\r
 279    pgFreeDmatrix(Lv,0,YN_1,0,XN_1);\r
 280    pgFreeDmatrix(g,0,YN_1,0,XN_1);\r
 281    pgFreeDmatrix(c1,0,YN_1,0,XN_1);\r
 282    pgFreeDmatrix(c2,0,YN_1,0,XN_1);\r
 283    pgFreeDmatrix(b,0,YN_1,0,XN_1);\r
 284 \r
 285    return;\r
 286 }\r
 287 \r
 288 void mexFunction( int nlhs, mxArray *plhs[], \r
 289                  int nrhs, const mxArray *prhs[] )\r
 290 {\r
 291    double  *f, *u, *v;\r
 292    int mrows, ncols;\r
 293    int m, n;\r
 294    double mu;\r
 295    int ITER;\r
 296    char buf[80];\r
 297 \r
 298    mexUnlock();\r
 299 \r
 300    /* Check for proper number of arguments */\r
 301    if ( nrhs < 1 || nrhs > 3) \r
 302       mexErrMsgTxt("GVFC.MEX: wrong number of input arguments.");\r
 303    if ( nlhs > 2 )\r
 304       mexErrMsgTxt("GVFC.MEX: wrong number of output arguments.");\r
 305    \r
 306    /* Obtain the dimension information of inputs */\r
 307    mrows = mxGetM(prhs[0]);   ncols = mxGetN(prhs[0]);\r
 308 \r
 309    /* Create matrix for the return argument */\r
 310    plhs[0] = mxCreateDoubleMatrix(mrows, ncols, mxREAL);\r
 311    plhs[1] = mxCreateDoubleMatrix(mrows, ncols, mxREAL);\r
 312    \r
 313    /* Assign pointers to each input and output */\r
 314    f = (double *)mxGetPr(prhs[0]); \r
 315    u = (double *)mxGetPr(plhs[0]);\r
 316    v = (double *)mxGetPr(plhs[1]);\r
 317    \r
 318    if (nrhs >= 2 && nrhs <= 3) {\r
 319       /* The input must be a noncomplex scalar double.*/\r
 320       m = mxGetM(prhs[1]);\r
 321       n = mxGetN(prhs[1]);\r
 322       if( !mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]) ||\r
 323          !(m==1 && n==1) ) {\r
 324          mexErrMsgTxt("Input mu must be a noncomplex scalar double.");\r
 325       }\r
 326       mu = mxGetScalar(prhs[1]);\r
 327    }\r
 328 \r
 329    if (nrhs == 3) {\r
 330       /* The input must be an integer.*/\r
 331       m = mxGetM(prhs[2]);\r
 332       n = mxGetN(prhs[2]);\r
 333       if( !mxIsDouble(prhs[2]) || mxIsComplex(prhs[2]) ||\r
 334          !(m==1 && n==1) ) {\r
 335          mexErrMsgTxt("Input mu must be an integer.");\r
 336       }\r
 337       ITER = mxGetScalar(prhs[2]);\r
 338    }\r
 339 \r
 340 \r
 341    /* set default mu and ITER */\r
 342    if (nrhs <= 1) \r
 343          mu = 0.2;\r
 344    if (nrhs <= 2)\r
 345       ITER = PG_MAX(mrows, ncols);\r
 346       \r
 347    /* Call the distedt subroutine. */\r
 348    GVFC(mrows, ncols, f, u, v, mu, ITER);\r
 349 \r
 350    /* done */\r
 351    return;\r
 352 }\r
 353 \r
 354 \r