Sadie_TREES_Ext.c Source File

00001 #include <tcl.h>
00002 #include <tk.h>
00003 #include "sadie.h"
00004 #include "proto.h"
00005 #include "sadie_byte.h"
00006 #include "sadie_short.h"
00007 #include "trees.h"
00008 #include "mosaic.h"
00009 #include "chain.h"
00010 #include "histogram.h"
00011 #include <sys/time.h>
00012 
00013 
00014 /*----------------------------------------------------------------------------*/
00015 /*-General Information--------------------------------------------------------*/
00016 /*                                                                            */
00017 /*   Global Sadie variables that must be set from within Tcl/Tk.              */
00018 /*                                                                            */
00019 /*                                                                            */
00020 /*----------------------------------------------------------------------------*/
00021 extern  short   nlev;
00022 extern  short   csize;
00023 extern  double  weight;
00024 extern  double  *count;
00025 extern  PIXEL   gain;
00026 extern  PIXEL   bias;
00027 extern  PIXEL   gmin;
00028 extern  PIXEL   gmax;
00029 extern  PIXEL   thresh;
00030 extern  PIXEL   gbrk[2][4];
00031 extern  PIXEL   *table;
00032 
00033 
00034 /*----------------------------------------------------------------------------*/
00035 /*-General Information--------------------------------------------------------*/
00036 /*                                                                            */
00037 /*   Variables required for analyst interaction.                              */
00038 /*                                                                            */
00039 /*                                                                            */
00040 /*----------------------------------------------------------------------------*/
00041 int  locx;      /* Column coordinate of the point of user click.              */
00042 int  locy;      /* Row coordinate of the point of user click.                 */
00043 int  zoomx;     /* Column coordinate of top-left corner of zoom window.       */
00044 int  zoomy;     /* Row coordinate of top-left corner of zoom window.          */
00045 int  retain;    /* Variable that tells whether or not to retain the ring.     */
00046 
00047 
00048 /*----------------------------------------------------------------------------*/
00049 /*-General Information--------------------------------------------------------*/
00050 /*                                                                            */
00051 /*   Variables necessary for computing vstart and vend with bresh algorithm.  */
00052 /*                                                                            */
00053 /*                                                                            */
00054 /*----------------------------------------------------------------------------*/
00055 extern int **x_array_buffer; /* The data where info is stored.                */
00056 extern int point;            /* Counter for points used.                      */
00057 extern int count_point;      /* Flag to decide whether to count or fill array.*/
00058 
00059 
00060 /*----------------------------------------------------------------------------*/
00061 /*-General Information--------------------------------------------------------*/
00062 /*                                                                            */
00063 /*   This procedure finds the rings in the mosaic image.                      */
00064 /*                                                                            */
00065 /*                                                                            */
00066 /*----------------------------------------------------------------------------*/
00067 int Sadie_TREES_Ext_TrackRingsCmd(ClientData client_data, Tcl_Interp* interp, int argc, char *argv[])
00068 {
00069     /* General declarations */
00070     int   addrstr;
00071     char  msg[SLEN];
00072     char *mode = NULL;
00073     struct  timeval start, end, t1, t2;
00074     HISTOGRAM  globalhist;
00075     HISTOGRAM  **localhist;
00076     IMAGE_BYTE  *byteimg=NULL;
00077     IMAGE_BYTE  *iMag=NULL;
00078     IMAGE_BYTE  *iDir=NULL;
00079     MOSAIC_INDEX  *index=NULL;
00080     
00081     /* Declarations for edge tracking go in here */
00082     register int i, j, k, l, m, n;
00083     int  max, found_max, done, nmax;
00084     int  startx, starty;
00085     IMAGE_BYTE  *out = NULL;
00086 
00087     /* Declarations for chain-coding go in here */
00088     int  discard, x, y, count;
00089     int  numlinks, pixcount, maxx, maxy;
00090     int  prevx, prevy, currx, curry, reset;
00091     unsigned char  code;
00092     CHAIN  *chainmap;
00093     LIST_NODE  *tempLnode;
00094     LIST_NODE  **prevLnode;
00095     LIST_NODE  *traverseLnode;
00096     CHAIN_NODE  **prevCnode;
00097     CHAIN_NODE  *tempCnode;
00098 
00099 
00100     /* Get the arguments */
00101     if( argc <= 3 ) {
00102         sprintf(msg, "Usage: Sadie_TREES_TrackRings byteimg index mode");
00103         Tcl_SetResult(interp, msg, TCL_VOLATILE);
00104         return TCL_ERROR;
00105     }
00106     
00107     sscanf(argv[1], "%x", &addrstr);
00108     byteimg = (IMAGE_BYTE *) addrstr;
00109     
00110     sscanf(argv[2], "%x", &addrstr);
00111     index = (MOSAIC_INDEX *) addrstr;
00112 
00113     mode = argv[3];
00114 
00115     /* start timing */
00116     gettimeofday(&start, NULL);
00117     
00118     /* This will return the gradient magnitude and direction mosaics */
00119     gettimeofday(&t1, NULL);
00120     MOSAIC_BYTE_TRACKRINGS(byteimg, index, &iMag, &iDir);
00121     gettimeofday(&t2, NULL);
00122     printf("Total time for gradient = %ld ms.\n\n", delay(t1, t2));
00123     
00124     /* Now we'll suppress noise edge pixels */
00125     printf("Post-gradient processing\n");
00126     gettimeofday(&t1, NULL);
00127     MOSAIC_BYTE_SUPPRESS(byteimg, index, iMag, iDir, &globalhist, &localhist);
00128     gettimeofday(&t2, NULL);
00129     printf("Total post-gradient time = %ld ms.\n\n", delay(t1, t2));
00130 
00131 
00132     /* check input image */
00133     if (!CHECKIMG_BYTE(byteimg)) {
00134         sprintf(msg," Can't identify mosaic.");
00135         Tcl_SetResult(interp, msg, TCL_VOLATILE);
00136         return TCL_ERROR;
00137     }
00138     if (!CHECKIMG_BYTE(iMag)) {
00139         sprintf(msg," Can't identify magnitude mosaic.");
00140         Tcl_SetResult(interp, msg, TCL_VOLATILE);
00141         return TCL_ERROR;
00142     }
00143     if (!CHECKIMG_BYTE(iDir)) {
00144         sprintf(msg," Can't identify direction mosaic.");
00145         Tcl_SetResult(interp, msg, TCL_VOLATILE);
00146         return TCL_ERROR;
00147     }
00148     if (!CHECKIMG_MOSAIC_INDEX(index)) {
00149         sprintf(msg," Can't identify mosaic index");
00150         Tcl_SetResult(interp, msg, TCL_VOLATILE);
00151         return TCL_ERROR;
00152     }
00153     if ((iMag->nlin != iDir->nlin) || (iMag->npix != iDir->npix)) {
00154         sprintf(msg," Dimensions of Ring Mosaic and Link Mosaic don't match.");
00155         Tcl_SetResult(interp, msg, TCL_VOLATILE);
00156         return TCL_ERROR;
00157     }
00158 
00159 
00160     gettimeofday(&t1, NULL);
00161     
00162     /* Allocate memory for storing the ring image */
00163     GETMEM_BYTE(iMag->nbnd, iMag->nlin, iMag->npix, &out);
00164     
00165     /* Allocate memory for chainmap */
00166     chainmap = (CHAIN *)malloc(sizeof(CHAIN));
00167     
00168     /* Initialize all the pixels to zero */
00169     for(i = 0; i < out->nbnd; i++)
00170         for(j = 0; j < out->nlin; j++)
00171             for(k = 0; k < out->npix; k++)
00172                 out->data[i][j][k] = (BYTE)0;
00173     
00174     
00175     /* Write chain header */
00176     chainmap->nlin = index->nlin;
00177     chainmap->npix = index->npix;
00178     memcpy((char *)(chainmap->header), (char *)byteimg->text, TLEN - 1);
00179     chainmap->list = NULL;
00180     prevLnode = (LIST_NODE **)&(chainmap->list);
00181 
00182 
00183     /* Edge tracking */
00184     max = 5;
00185     found_max = 0;
00186     count = -1;
00187     for(k = index->hstart + 1; k <= index->hend; k++) {
00188         j = index->vstart[k] - index->voffset[k] + 5;
00189 
00190         /* Identify starting point of a ring */
00191         if (iMag->data[0][j][k] > max) {
00192             max = iMag->data[0][j][k];
00193             startx = index->vstart[k] + 0;
00194             starty = k;
00195             found_max = 1;
00196         }
00197         
00198         /* Start tracking the ring */
00199         if (found_max && iMag->data[0][j][k] <= 5) {
00200             found_max = 0;
00201             max = 5;
00202             done = 0;
00203             discard = numlinks = 0;
00204             pixcount = 1;
00205 
00206             tempLnode = (LIST_NODE *)malloc(sizeof(LIST_NODE));
00207             tempLnode->startj = startx;
00208             tempLnode->startk = starty;
00209             tempLnode->ring_type = NORMAL;
00210             tempLnode->chain = NULL;
00211             tempLnode->next = NULL;
00212             prevCnode = (CHAIN_NODE **)&(tempLnode->chain);
00213 
00214             /* For as long as the ring is not fully coded */
00215             while (!done && !discard) {
00216                     
00217                 /* Check for the end of ROI */
00218                 if (startx == index->vend[starty] - 5 ||
00219                     starty == index->hstart + 5 || starty == index->hend - 5) {
00220                     done = 1;
00221                     break;
00222                 }
00223                 
00224                 /* Find the neighbour with the highest gradient magnitude */
00225                 if (HIST_FIND_MAX_NEIGH(iMag, byteimg, index, startx, starty, &x, &y)) {
00226                     
00227                     /* Check for an intersection with an already coded ring */
00228                     if (out->data[0][x-index->voffset[y]][y] == 255) {
00229                         discard = 1;
00230                         break;
00231                         
00232                     } else { /* No intersection, then we code the pixel */
00233 
00234                         out->data[0][x - index->voffset[y]][y] = 255;
00235                         CHAIN_COMPUTE_CODE(startx, starty, x, y, &code);
00236                         
00237                         startx = x;
00238                         starty = y;
00239                         tempCnode = (CHAIN_NODE *)malloc(sizeof(CHAIN_NODE));
00240                         tempCnode->code = code;
00241                         tempCnode->valid = 1;
00242                         tempCnode->next = NULL;
00243                         (*prevCnode) = tempCnode;
00244                         prevCnode = (CHAIN_NODE **)&(tempCnode->next);
00245                         pixcount++;
00246                         
00247                     }
00248                     
00249                     
00250                 } else { /* No neighbours */
00251 
00252                     if ( (float)numlinks/pixcount > 0.25 && pixcount > 15) { 
00253                         /* Too many links....discard ring */
00254                         discard = 1;
00255                         break;
00256                         
00257                     } else { /* Try bridging a long gap */
00258 
00259                         nmax = 0;
00260                         for(m = starty - 5; m <= starty + 5; m++) {
00261                             if (m < index->hstart || m > index->hend)
00262                                 continue;
00263                             for(n = startx + 1; n <= startx + 10; n++) {
00264                                 if (n < index->vstart[m] || n > index->vend[m])
00265                                     continue;
00266                                 if (iMag->data[0][n - index->voffset[m]][m] > nmax) {
00267                                     nmax = iMag->data[0][n - index->voffset[m]][m];
00268                                     maxx = n;
00269                                     maxy = m;
00270                                 }
00271                             }
00272                         }
00273                         
00274                         if (nmax == 0) { /* If unable to brigde */
00275                             
00276                             /* End if the current pixel is close to the frame boundary */
00277                             if (startx >= index->vend[starty] - 5 ||
00278                                 starty == index->hstart + 5 || starty == index->hend - 5) {
00279                                 done = 1;
00280                                 break;
00281                             } else if (mode[0] == 'I' && pixcount > 50) { /* User Interaction */
00282 
00283                                 /* Call proc to display region in zoom panel */
00284                                 sprintf(msg, "interaction_zoom %d %d", starty, startx);
00285                                 if (Tcl_Eval(interp, msg) != TCL_OK)
00286                                     return TCL_ERROR;
00287 
00288                                 /* Decode the chain so far and display on zoom window */
00289                                 traverseLnode = tempLnode;
00290                                 prevx = traverseLnode->startk;
00291                                 prevy = traverseLnode->startj;
00292                                 tempCnode = (CHAIN_NODE *)traverseLnode->chain;
00293                                 while(tempCnode != NULL) {
00294 
00295                                     CHAIN_COMPUTE_COORDINATES(prevx, prevy, tempCnode->code, &currx, &curry);
00296                                     sprintf(msg, "%s create line %d %d %d %d -width 1 -fill %s -tag %s", ".trees_main.zoompanel.zoomframe.canvas", prevx-zoomx, prevy-zoomy, currx-zoomx, curry-zoomy, "red", "interact");
00297 
00298                                     if (Tcl_Eval(interp, msg) != TCL_OK)
00299                                         return TCL_ERROR;
00300 
00301                                     prevx = currx;
00302                                     prevy = curry;
00303                                     tempCnode = (CHAIN_NODE *)tempCnode->next;
00304                                 }
00305 
00306                                 /* Make sure the ring is visible */
00307                                 sprintf(msg, "update idletasks");
00308                                 if (Tcl_Eval(interp, msg) != TCL_OK)
00309                                     return TCL_ERROR;
00310 
00311                                 /* Prompt and wait for user input */
00312                                 sprintf(msg, "trees_interaction_init");
00313                                 if (Tcl_Eval(interp, msg) != TCL_OK)
00314                                     return TCL_ERROR;
00315 
00316                                 /* Clear up the zoom window */
00317                                 sprintf(msg, "catch {.trees_main.zoompanel.zoomframe.canvas delete all}");
00318                                 
00319                                 maxy = locx + zoomx;
00320                                 maxx = locy + zoomy;
00321                                 
00322                                 /* Discard on analyst instruction */
00323                                 if (!retain) {
00324                                     discard = 1;
00325                                     break;
00326                                 }
00327                                 
00328                             } else {
00329                                 discard = 1;
00330                                 break;
00331                             }
00332                         }
00333                         
00334                         /* Count number of points in line */
00335                         count_point = TRUE;
00336                         point = 0;
00337                         brshnm(starty, startx, maxy, maxx);
00338                         
00339                         /* Allocate memory */
00340                         x_array_buffer = (int **) malloc((long)point*sizeof(int *));
00341                         for  (l = 0; l < point; l++)
00342                             x_array_buffer[l] = (int  *) calloc(2, sizeof(int));
00343                         
00344                         /* Fill array */ 
00345                         count_point = FALSE;
00346                         point = 0;
00347                         brshnm(starty, startx, maxy, maxx);
00348                         
00349                         /* Link */
00350                         for (l = 0; l < point; l++) {
00351                             
00352                             x = x_array_buffer[l][1];
00353                             y = x_array_buffer[l][0];
00354                             out->data[0][x - index->voffset[y]][y] = 255;
00355                             CHAIN_COMPUTE_CODE(startx, starty, x, y, &code);
00356                             
00357                             startx = x;
00358                             starty = y;
00359                             tempCnode = (CHAIN_NODE *)malloc(sizeof(CHAIN_NODE));
00360                             tempCnode->code = code;
00361                             tempCnode->valid = 0;
00362                             tempCnode->next = NULL;
00363                             (*prevCnode) = tempCnode;
00364                             prevCnode = (CHAIN_NODE **)&(tempCnode->next);
00365                             pixcount++;
00366                             numlinks++;
00367                             
00368                         }
00369                         
00370                     }
00371                     
00372                 }
00373                 
00374                 
00375             } /* While loop */
00376             
00377             if (!discard && pixcount > 15) { /* Retain ring */
00378                 (*prevLnode) = tempLnode;
00379                 prevLnode = (LIST_NODE **)&(tempLnode->next);
00380                 tempLnode->ring_id = count++;
00381             } else {
00382                 /* Reset the current chain to zero */
00383                 do {
00384                     reset = 0;
00385                     out->data[0][startx - index->voffset[starty]][starty] = 0;
00386                     for(l = starty + 1; l > starty - 2; l--) {
00387                         if (l < index->hstart || l > index->hend)
00388                             continue;
00389                         for(m = startx - 1; m < startx + 1; m++) {
00390                             if (m < index->vstart[l] || m > index->vend[l])
00391                                 continue;
00392                             if (out->data[0][m - index->voffset[l]][l] == 255) {
00393                                 out->data[0][m - index->voffset[l]][l] = 0;
00394                                 startx = m;
00395                                 starty = l;
00396                                 reset = 1;
00397                                 break;
00398                             }
00399                         }
00400                         if (reset) break;
00401                     }
00402                     if (!reset) break;
00403                 } while (startx != index->vstart[starty] + 1);
00404                 
00405             }
00406             
00407         } /* if condition */
00408         
00409     } /* Outer for loop */
00410     
00411     gettimeofday(&t2, NULL);
00412     printf("Tracking and coding rings = %ld ms.\n\n", delay(t1, t2));
00413 
00414     gettimeofday(&end, NULL);
00415     printf("Total time for finding rings = %ld ms.\n\n", delay(start, end));
00416 
00417     printf("Encoded %d rings\n", count+1);
00418 
00419     /* Return only the chainmap */
00420     sprintf(msg, "%x", (int)chainmap);
00421     Tcl_SetResult(interp, msg, TCL_VOLATILE);
00422 
00423     
00424     return TCL_OK;
00425 }
00426 
00427 
00428 /*----------------------------------------------------------------------------*/
00429 /*-General Information--------------------------------------------------------*/
00430 /*                                                                            */
00431 /*   This procedure conveys the user inputs to the algorithm.                 */
00432 /*                                                                            */
00433 /*                                                                            */
00434 /*----------------------------------------------------------------------------*/
00435 int Sadie_TREES_Ext_SetUserInputCmd(ClientData client_data, Tcl_Interp* interp, int argc, char *argv[])
00436 {
00437     double  row, col;
00438     char msg[SLEN];
00439     
00440     if( argc <= 5 ) {
00441         sprintf(msg, "Usage: Sadie_TREES_SetUserInput retain locx locy zoomx zoomy");
00442         Tcl_SetResult(interp, msg, TCL_VOLATILE);
00443         return TCL_ERROR;
00444     }
00445     
00446     if( Tcl_GetInt (interp, argv[1], &retain) != TCL_OK ) {
00447         return TCL_ERROR;
00448     }
00449     if( Tcl_GetDouble (interp, argv[2], &col) != TCL_OK ) {
00450         return TCL_ERROR;
00451     }
00452     if( Tcl_GetDouble (interp, argv[3], &row) != TCL_OK ) {
00453         return TCL_ERROR;
00454     }
00455     if( Tcl_GetInt (interp, argv[4], &zoomx) != TCL_OK ) {
00456         return TCL_ERROR;
00457     }
00458     if( Tcl_GetInt (interp, argv[5], &zoomy) != TCL_OK ) {
00459         return TCL_ERROR;
00460     }
00461       
00462     /* User Input convey protocol */
00463     /* locx, locy, discard and read are globals */
00464     locx = (int)col;
00465     locy = (int)row;
00466     
00467     return TCL_OK;
00468 }
00469 
00470 
00471 /*----------------------------------------------------------------------------*/
00472 /*-General Information--------------------------------------------------------*/
00473 /*                                                                            */
00474 /*   This procedure initializes the the TREES Extension  procedures when      */
00475 /*   this package is loaded.                                                  */
00476 /*                                                                            */
00477 /*----------------------------------------------------------------------------*/
00478 int Sadie_TREES_Ext_Init(Tcl_Interp *interp) {
00479     
00480   Tcl_CreateCommand(interp, "Sadie_TREES_Ext_TrackRings", Sadie_TREES_Ext_TrackRingsCmd, 
00481         (ClientData) NULL, NULL);
00482   Tcl_CreateCommand(interp, "Sadie_TREES_Ext_SetUserInput", Sadie_TREES_Ext_SetUserInputCmd, 
00483         (ClientData) NULL, NULL);
00484 
00485   return TCL_OK;
00486 }
00487