fourierdesc.c Source File

00001 #include "sadie.h"
00002 #include "proto.h"
00003 
00004 
00005 #define UP 1
00006 #define DOWN 2
00007 #define LEFT 3
00008 #define RIGHT 4
00009 
00010 
00011 /*----------------------------------------------------------------------------*/
00012 /*-General Information--------------------------------------------------------*/
00013 /*                                                                            */
00014 /* This function computes 8-connected fourier descriptors for an input        */
00015 /*  labeled image.                                                            */
00016 /*                                                                            */
00017 /*----------------------------------------------------------------------------*/
00018 /*-Interface Information------------------------------------------------------*/
00019 void FOURIERDESC (
00020 IMAGE  *in       /*  I   Pointer to the input image.                          */
00021 /*----------------------------------------------------------------------------*/
00022 ) { register short i, j, k, n;
00023     char msg[SLEN];
00024     int regions = 0;
00025     int tempval;
00026     int curlbl;
00027     int foundpix, startx, starty, done, dir;
00028     double fourier[1000], sine[500], cosn[500];
00029     int edgecount;
00030     int option = FORWARD;
00031 
00032     /* if (TIMES) TIMING(T_FOURIERDESC); */
00033     if (NAMES) {
00034         MESSAGE('I',"");
00035         MESSAGE('I',"FOURIERDESC");
00036         MESSAGE('I',"");
00037         sprintf(msg," Input image:                       %s",in->text);
00038         MESSAGE('I',msg);
00039         MESSAGE('I'," ...............");
00040     }
00041 
00042     /* check input */
00043     if (!CHECKIMG(in)) {
00044         MESSAGE('E'," Can't identify image.");
00045         goto the_end;
00046     } else if (in->nbnd > 1) {
00047         MESSAGE('E'," Image must be single-band.");
00048         goto the_end;
00049     }
00050 
00051 
00052 
00053     /* Find out how many objects are in the input labeled image */
00054     /* This assumes that labels are assigned consecutively when 
00055         scanning LTR-TTB across the image.                      */
00056     RANGE(in);
00057     tempval = (int) rnd(in->gmin);
00058     for (j=0; j<in->nlin; j++) {
00059         for (k=0; k<in->npix; k++) {
00060             if ((int) rnd(in->data[0][j][k]) > tempval) {
00061                 tempval = (int) rnd(in->data[0][j][k]);
00062                 regions++;
00063             }
00064         }    
00065     }
00066     
00067     if (regions == 0) {
00068         MESSAGE('E'," Number of regions in image must be greater than zero.");
00069         goto the_end;
00070     }
00071 
00072     for (curlbl = (int)rnd(in->gmin) + 1; curlbl <= (int)rnd(in->gmax); curlbl++) {
00073         /* Find the first pixel in the region */
00074         foundpix = 0;
00075         edgecount = 0;
00076 
00077         for (j=0; j<in->nlin;j++) {
00078             for (k=0; k<in->npix; k++) {
00079                 tempval = (int) rnd(in->data[0][j][k]);
00080                 if (curlbl == tempval) {
00081                     foundpix = 1;
00082                     break;
00083                 }
00084             }
00085             
00086             if (foundpix == 1) {
00087                 break;
00088             }
00089         }
00090         
00091         /* Record the starting location */
00092         startx = k;
00093         starty = j;
00094         
00095         if (foundpix == 1) {
00096         
00097             done = 0;
00098             
00099             /* Find the first neighbor included in the region */
00100             if ((k>0) && ((int) rnd(in->data[0][j][k-1]) == curlbl)) { /*Check n1*/
00101                 dir = UP;
00102                 k = k - 1;
00103             } else if ((j>0) && (k>0) && ((int)rnd(in->data[0][j-1][k-1]) ==curlbl)) { /*Check n2*/
00104                 dir = UP;
00105                 k = k - 1;
00106                 j = j - 1;
00107             } else if ((j>0) && ((int)rnd(in->data[0][j-1][k]) ==curlbl)) { /*Check n3*/
00108                 dir = RIGHT;
00109                 j = j - 1;
00110             } else if ((j>0) && (k<in->npix-1) && ((int)rnd(in->data[0][j-1][k+1]) ==curlbl)) { /*Check n4*/
00111                 dir = RIGHT;
00112                 k = k + 1;
00113                 j = j - 1;
00114             } else if ((k<in->npix-1) && ((int)rnd(in->data[0][j][k+1]) ==curlbl)) { /*Check n5*/
00115                 dir = DOWN;
00116                 k = k + 1;
00117             } else if ((j<in->nlin-1) && (k<in->npix-1) && ((int)rnd(in->data[0][j+1][k+1]) ==curlbl)) { /*Check n6*/
00118                 dir = DOWN;
00119                 k = k + 1;
00120                 j = j + 1;
00121             } else if ((j<in->nlin-1) && ((int)rnd(in->data[0][j+1][k]) ==curlbl)) { /*Check n7*/
00122                 dir = LEFT;
00123                 j = j + 1;
00124             } else if ((j<in->nlin-1) && (k>0) && ((int)rnd(in->data[0][j+1][k-1]) ==curlbl)) { /*Check n8*/
00125                 dir = LEFT;
00126                 k = k - 1;
00127                 j = j + 1;
00128             } else {
00129                 /* This is a one pixel region */
00130                 done = 1;
00131                 sprintf(msg,"Region %d is a one pixel region", curlbl);
00132                 MESSAGE('I',msg);  
00133             }
00134 
00135             if (done != 1) {
00136                 fourier[2*edgecount] = (double) k;
00137                 fourier[(2*edgecount)+1] = (double) j;
00138                 edgecount++;
00139             }
00140 
00141         } else {
00142             /* No foreground objects were found */
00143             done = 1;
00144             sprintf(msg,"No foreground objects were found!");
00145             MESSAGE('I',msg);  
00146         }
00147 
00148         while (done != 1) {
00149             if (edgecount >= 500) {
00150                 sprintf(msg,"Too many edge pixels in region %d!", curlbl);
00151                 MESSAGE('I',msg);  
00152                 goto the_end;
00153             }
00154 
00155             if (dir == RIGHT) {
00156                 if ((k>0) && ((int) rnd(in->data[0][j][k-1]) == curlbl)) { /*Check n1*/
00157                     dir = UP;
00158                     k = k - 1;
00159                 } else if ((j>0) && (k>0) && ((int)rnd(in->data[0][j-1][k-1]) ==curlbl)) { /*Check n2*/
00160                     dir = UP;
00161                     k = k - 1;
00162                     j = j - 1;
00163                 } else if ((j>0) && ((int)rnd(in->data[0][j-1][k]) ==curlbl)) { /*Check n3*/
00164                     dir = RIGHT;
00165                     j = j - 1;
00166                 } else if ((j>0) && (k<in->npix-1) && ((int)rnd(in->data[0][j-1][k+1]) ==curlbl)) { /*Check n4*/
00167                     dir = RIGHT;
00168                     k = k + 1;
00169                     j = j - 1;
00170                 } else if ((k<in->npix-1) && ((int)rnd(in->data[0][j][k+1]) ==curlbl)) { /*Check n5*/
00171                     dir = DOWN;
00172                     k = k + 1;
00173                 } else if ((j<in->nlin-1) && (k<in->npix-1) && ((int)rnd(in->data[0][j+1][k+1]) ==curlbl)) { /*Check n6*/
00174                     dir = DOWN;
00175                     k = k + 1;
00176                     j = j + 1;
00177                 } else if ((j<in->nlin-1) && ((int)rnd(in->data[0][j+1][k]) ==curlbl)) { /*Check n7*/
00178                     dir = LEFT;
00179                     j = j + 1;
00180                 } else if ((j<in->nlin-1) && (k>0) && ((int)rnd(in->data[0][j+1][k-1]) ==curlbl)) { /*Check n8*/
00181                     dir = LEFT;
00182                     k = k - 1;
00183                     j = j + 1;
00184                 } else {
00185                     MESSAGE('E'," Could not find neighbor in region.");
00186                     sprintf(msg,"   curlbl = %d, k = %d, j = %d, dir = %d", curlbl,k,j,dir);
00187                     MESSAGE('E',msg);  
00188                     goto the_end;
00189                 }
00190             } else if (dir == LEFT) {
00191                 if ((k<in->npix-1) && ((int)rnd(in->data[0][j][k+1]) ==curlbl)) { /*Check n1*/
00192                     dir = DOWN;
00193                     k = k + 1;
00194                 } else if ((j<in->nlin-1) && (k<in->npix-1) && ((int)rnd(in->data[0][j+1][k+1]) ==curlbl)) { /*Check n2*/
00195                     dir = DOWN;
00196                     k = k + 1;
00197                     j = j + 1;
00198                 } else if ((j<in->nlin-1) && ((int)rnd(in->data[0][j+1][k]) ==curlbl)) { /*Check n3*/
00199                     dir = LEFT;
00200                     j = j + 1;
00201                 } else if ((j<in->nlin-1) && (k>0) && ((int)rnd(in->data[0][j+1][k-1]) ==curlbl)) { /*Check n4*/
00202                     dir = LEFT;
00203                     k = k - 1;
00204                     j = j + 1;
00205                 } else if ((k>0) && ((int) rnd(in->data[0][j][k-1]) == curlbl)) { /*Check n5*/
00206                     dir = UP;
00207                     k = k - 1;
00208                 } else if ((j>0) && (k>0) && ((int)rnd(in->data[0][j-1][k-1]) ==curlbl)) { /*Check n6*/
00209                     dir = UP;
00210                     k = k - 1;
00211                     j = j - 1;
00212                 } else if ((j>0) && ((int)rnd(in->data[0][j-1][k]) ==curlbl)) { /*Check n7*/
00213                     dir = RIGHT;
00214                     j = j - 1;
00215                 } else if ((j>0) && (k<in->npix-1) && ((int)rnd(in->data[0][j-1][k+1]) ==curlbl)) { /*Check n8*/
00216                     dir = RIGHT;
00217                     k = k + 1;
00218                     j = j - 1;
00219                 } else {
00220                     MESSAGE('E'," Could not find neighbor in region.");
00221                     sprintf(msg,"   curlbl = %d, k = %d, j = %d, dir = %d", curlbl,k,j,dir);
00222                     MESSAGE('E',msg);  
00223                     goto the_end;
00224                 }
00225             } else if (dir == UP) {
00226                 if ((j<in->nlin-1) && ((int)rnd(in->data[0][j+1][k]) ==curlbl)) { /*Check n1*/
00227                     dir = LEFT;
00228                     j = j + 1;
00229                 } else if ((j<in->nlin-1) && (k>0) && ((int)rnd(in->data[0][j+1][k-1]) ==curlbl)) { /*Check n2*/
00230                     dir = LEFT;
00231                     k = k - 1;
00232                     j = j + 1;
00233                 } else if ((k>0) && ((int) rnd(in->data[0][j][k-1]) == curlbl)) { /*Check n3*/
00234                     dir = UP;
00235                     k = k - 1;
00236                 } else if ((j>0) && (k>0) && ((int)rnd(in->data[0][j-1][k-1]) ==curlbl)) { /*Check n4*/
00237                     dir = UP;
00238                     k = k - 1;
00239                     j = j - 1;
00240                 } else if ((j>0) && ((int)rnd(in->data[0][j-1][k]) ==curlbl)) { /*Check n5*/
00241                     dir = RIGHT;
00242                     j = j - 1;
00243                 } else if ((j>0) && (k<in->npix-1) && ((int)rnd(in->data[0][j-1][k+1]) ==curlbl)) { /*Check n6*/
00244                     dir = RIGHT;
00245                     k = k + 1;
00246                     j = j - 1;
00247                 } else if ((k<in->npix-1) && ((int)rnd(in->data[0][j][k+1]) ==curlbl)) { /*Check n7*/
00248                     dir = DOWN;
00249                     k = k + 1;
00250                 } else if ((j<in->nlin-1) && (k<in->npix-1) && ((int)rnd(in->data[0][j+1][k+1]) ==curlbl)) { /*Check n8*/
00251                     dir = DOWN;
00252                     k = k + 1;
00253                     j = j + 1;
00254                 } else {
00255                     MESSAGE('E'," Could not find neighbor in region.");
00256                     sprintf(msg,"   curlbl = %d, k = %d, j = %d, dir = %d", curlbl,k,j,dir);
00257                     MESSAGE('E',msg);  
00258                     goto the_end;
00259                 }
00260             } else {    /* dir == DOWN */
00261                 if ((j>0) && ((int)rnd(in->data[0][j-1][k]) ==curlbl)) { /*Check n1*/
00262                     dir = RIGHT;
00263                     j = j - 1;
00264                 } else if ((j>0) && (k<in->npix-1) && ((int)rnd(in->data[0][j-1][k+1]) ==curlbl)) { /*Check n2*/
00265                     dir = RIGHT;
00266                     k = k + 1;
00267                     j = j - 1;
00268                 } else if ((k<in->npix-1) && ((int)rnd(in->data[0][j][k+1]) ==curlbl)) { /*Check n3*/
00269                     dir = DOWN;
00270                     k = k + 1;
00271                 } else if ((j<in->nlin-1) && (k<in->npix-1) && ((int)rnd(in->data[0][j+1][k+1]) ==curlbl)) { /*Check n4*/
00272                     dir = DOWN;
00273                     k = k + 1;
00274                     j = j + 1;
00275                 } else if ((j<in->nlin-1) && ((int)rnd(in->data[0][j+1][k]) ==curlbl)) { /*Check n5*/
00276                     dir = LEFT;
00277                     j = j + 1;
00278                 } else if ((j<in->nlin-1) && (k>0) && ((int)rnd(in->data[0][j+1][k-1]) ==curlbl)) { /*Check n6*/
00279                     dir = LEFT;
00280                     k = k - 1;
00281                     j = j + 1;
00282                 } else if ((k>0) && ((int) rnd(in->data[0][j][k-1]) == curlbl)) { /*Check n7*/
00283                     dir = UP;
00284                     k = k - 1;
00285                 } else if ((j>0) && (k>0) && ((int)rnd(in->data[0][j-1][k-1]) ==curlbl)) { /*Check n8*/
00286                     dir = UP;
00287                     k = k - 1;
00288                     j = j - 1;
00289                 } else {
00290                     MESSAGE('E'," Could not find neighbor in region.");
00291                     sprintf(msg,"   curlbl = %d, k = %d, j = %d, dir = %d", curlbl,k,j,dir);
00292                     MESSAGE('E',msg);  
00293                     goto the_end;
00294                 }
00295             }
00296 
00297             fourier[2*edgecount] = (double) k;
00298             fourier[(2*edgecount)+1] = (double) j;
00299             edgecount++;
00300 
00301             if ((k==startx) && (j==starty)) {
00302                 done = 1;
00303 
00304                 /* compute the sine and cosine terms */
00305                 for (i=0; i<edgecount; i++) {
00306                     sine[i] = (double)option*sin((2.0*PI*(double)i)/(double)edgecount);
00307                     cosn[i] =                cos((2.0*PI*(double)i)/(double)edgecount);
00308                 }
00309 
00310 
00311                 FFT1D(option,edgecount,&sine,&cosn,&fourier);
00312                 sprintf(msg,"Fourier Descriptors for Region %d:", curlbl);
00313                 MESSAGE('I',msg);
00314 
00315                 strcpy(msg, " ");
00316                 i=strlen(msg);
00317                 for (k=0; k<edgecount; k++,i=strlen(msg)) {
00318                     sprintf(msg+i," a(%d) = %.2f+j%.2f &",k,fourier[2*k],fourier[(2*k)+1]);
00319                 
00320                     if (k>0 && (((k+1)%4) == 0)) {
00321                         MESSAGE('I',msg);
00322                         strcpy(msg, " ");
00323                     }
00324                 }
00325                 MESSAGE('I',msg);
00326             }
00327             
00328         }
00329     }
00330 
00331 
00332     
00333     
00334 
00335     the_end:
00336     /* if (TIMES) TIMING(T_EXIT); */
00337 
00338 }
00339 
00340 
00341 
00342 
00343 
00344