link_doublethresh.c Source File

00001 #include "sadie.h"
00002 #include "proto.h"
00003 
00004 /* #define LINK_DOUBLETHRESH_DEBUG */
00005 
00006 #define UP 1
00007 #define DOWN 2
00008 #define LEFT 3
00009 #define RIGHT 4
00010 
00011 
00012 /* variables necessary for drawing lines with bresh algorithm */
00013 extern int **x_array_buffer;     /* The data where info is stored      */
00014 extern int point;                /* Counter for points used            */
00015 extern int count_point;          /* Flag decide to count or fill array */
00016 
00017 
00018 typedef struct {                        
00019     short topx;
00020     short topy;
00021     short bottomx;
00022     short bottomy;
00023     short toplinkdone;
00024     short bottomlinkdone;
00025 } SEGMENT_DATA;
00026 
00027 
00028 /*----------------------------------------------------------------------------*/
00029 /*-General Information--------------------------------------------------------*/
00030 /*                                                                            */
00031 /* This function performs double threshold linking.  All edge fragments from  */
00032 /* thresh1 will be kept in the output.  All edge fragments from thresh2 that  */
00033 /* overlap fragments from thresh1 will also be kept in the output.  If        */
00034 /* maxlinkdist is greater than 0, and fragments from thresh2 with an endpoint */
00035 /* that is within maxlinkdist of the endpoint of a fragment from thresh1      */
00036 /* will be kept in the output edge map.                                       */
00037 /*                                                                            */
00038 /*----------------------------------------------------------------------------*/
00039 /*-Interface Information------------------------------------------------------*/
00040 void LINK_DOUBLETHRESH (
00041 IMAGE  *thresh1,    /*  I   Pointer to the input image with T1.               */
00042 IMAGE  *thresh2,    /*  I   Pointer to the input image with T2 < T1.          */
00043 int    maxlinkdist, /*  I   The maximum distance between fragments that       */
00044                     /*         can be linked    (typically 0 or 5)            */
00045 IMAGE  **fillgaps,  /*  O   Address of a pointer to the image representing    */
00046                     /*         ring fragment connections                      */
00047 IMAGE  **out        /*  O   Pointer to the output image address.              */
00048 /*----------------------------------------------------------------------------*/
00049 ) { register short i, j, k, l, n;
00050     char msg[SLEN];
00051     IMAGE *t1=NULL, *t2=NULL;       /* Mutually exclusive versions of thresh1 and thresh2 */
00052     IMAGE *t1lbl=NULL, *t2lbl=NULL; /* Labeled versions of t1 and t2 */
00053     SEGMENT_DATA *t1_data=NULL, *t2_data=NULL;
00054     int* t1include=NULL;
00055     int* t2include=NULL;
00056     
00057     int nbrhoodsize = (maxlinkdist+1)*(maxlinkdist+maxlinkdist+1);
00058     
00059     /* Memory to store information about t2 fragments in upper neighborhood around endpoints */
00060     int t2upnbrchk[nbrhoodsize][2];
00061     int t2upnbrcnt;         /* Number of t2 fragments in upper neighborhood */
00062     
00063     /* Memory to store information about t1 fragments in upper neighborhood around endpoints */
00064     int t1upnbrchk[nbrhoodsize][2];
00065     int t1upnbrcnt;         /* Number of t1 fragments in upper neighborhood */
00066 
00067     /* Memory to store information about t1 fragments in lower neighborhood around endpoints */
00068     int t1downnbrchk[nbrhoodsize][2];
00069     int t1downnbrcnt;         /* Number of t1 fragments in lower neighborhood */
00070 
00071     int min_x, min_y, max_x, max_y, dist_y, dist_x;
00072     int curt1lbl, curt2lbl, closet1lbl, closet2lbl;
00073     int startx, starty, endx, endy; /* Variables for drawing connections between fragments */
00074     int done;
00075     int connect;
00076 
00077 
00078     /* check input */
00079     if (!CHECKIMG(thresh1)) {
00080         MESSAGE('E'," Can't identify image.");
00081         goto the_end;
00082     } else if (!CHECKIMG(thresh2)) {
00083         MESSAGE('E'," Can't identify image.");
00084         goto the_end;
00085     }
00086     
00087     /* if (TIMES) TIMING(T_LINK_DOUBLETHRESH); */
00088     if (NAMES) {
00089         MESSAGE('I',"");
00090         MESSAGE('I',"LINK_DOUBLETHRESH");
00091         MESSAGE('I',"");
00092         sprintf(msg," Threshold 1 image:                     %s",thresh1->text);
00093         MESSAGE('I',msg);
00094         MESSAGE('I',"");
00095         sprintf(msg," Threshold 2 image:                     %s",thresh2->text);
00096         MESSAGE('I',msg);
00097         MESSAGE('I'," ...............");
00098     }
00099 
00100     /* check input */
00101     if (thresh1->nbnd > 1) {
00102         MESSAGE('E'," Only using first band of first threshold image.");
00103     }
00104     if (thresh2->nbnd > 1) {
00105         MESSAGE('E'," Only using first band of second threshold image.");
00106     }
00107     if ((thresh1->nlin != thresh2->nlin) || (thresh1->npix != thresh2->npix)) {
00108         MESSAGE('E'," Input images must be same size.");
00109         goto the_end;
00110     }
00111     
00112     /* Label the two images */
00113     CMPLBL8(thresh1,&t1lbl);
00114     if (!CHECKIMG(t1lbl)) {
00115         MESSAGE('E'," Problem labeling thresh1 image.");
00116         goto the_end;
00117     }
00118     CMPLBL8(thresh2,&t2lbl);
00119     if (!CHECKIMG(t2lbl)) {
00120         MESSAGE('E'," Problem labeling thresh2 image.");
00121         goto the_end;
00122     }
00123     RANGE(t1lbl);
00124     RANGE(t2lbl);
00125     
00126     
00127     
00128     
00129     /* Ensure that the two threshold images are mutually exclusive */
00130     /*   Create t1 and t2 for this purpose */
00131     
00132     /* create t1 and t2 images of appropriate size */
00133     if (!CHECKIMG(t1)) GETMEM(1,thresh1->nlin,thresh1->npix,&t1);
00134     if (!t1) goto the_end;
00135     if (!CHECKIMG(t2)) GETMEM(1,thresh1->nlin,thresh1->npix,&t2);
00136     if (!t2) goto the_end;
00137 
00138     if (!(t1include=(int *)malloc((int)rnd(t2lbl->gmax)*sizeof(int)))) { 
00139         MESSAGE('E'," Could not allocate memory for t1include list.");
00140         goto the_end;
00141     }
00142     
00143     for (curt2lbl = (int)rnd(t2lbl->gmin); curt2lbl < (int)rnd(t2lbl->gmax); curt2lbl++) {
00144         /* Initialize with assumption that nothing from thresh2 is included in thresh1 */
00145         t1include[curt2lbl] = 0;
00146     }
00147     
00148     for (j=0; j<thresh1->nlin; j++) {
00149         for (k=0; k<thresh1->npix; k++) {
00150             if (t1lbl->data[0][j][k] > (PIXEL)0.0) {
00151                 /* curt2lbl from thresh2 is included in thresh1 */
00152                 curt2lbl = (int)rnd(t2lbl->data[0][j][k]) - 1;
00153                 
00154                 if (curt2lbl >= 0) {
00155                     t1include[curt2lbl] = 1;
00156                 }
00157             }
00158         }
00159     }
00160 
00161     for (j=0; j<thresh1->nlin; j++) {
00162         for (k=0; k<thresh1->npix; k++) {
00163             if (t2lbl->data[0][j][k] > (PIXEL)0.0) {
00164                 curt2lbl = (int)rnd(t2lbl->data[0][j][k]) - 1;
00165                 
00166                 if (t1include[curt2lbl] == 1) {
00167                     t1->data[0][j][k] = (PIXEL)1.0;
00168                     t2->data[0][j][k] = (PIXEL)0.0;
00169                 } else {
00170                     t1->data[0][j][k] = (PIXEL)0.0;
00171                     t2->data[0][j][k] = (PIXEL)1.0;
00172                 }
00173             
00174             } else {
00175                 t1->data[0][j][k] = (PIXEL)0.0;
00176                 t2->data[0][j][k] = (PIXEL)0.0;
00177             }
00178         }
00179     }
00180     
00181     if (t1include) free(t1include);
00182     t1include = NULL;
00183     
00184     
00185 
00186 
00187     /* Re-label t1 and t2 */
00188     if (CHECKIMG(t1lbl))  RELIMG(&t1lbl);
00189     if (CHECKIMG(t2lbl))  RELIMG(&t2lbl);
00190     CMPLBL8(t1,&t1lbl);
00191     if (!CHECKIMG(t1lbl)) {
00192         MESSAGE('E'," Problem labeling t1 image.");
00193         goto the_end;
00194     }
00195     CMPLBL8(t2,&t2lbl);
00196     if (!CHECKIMG(t2lbl)) {
00197         MESSAGE('E'," Problem labeling t2 image.");
00198         goto the_end;
00199     }
00200     RANGE(t1lbl);
00201     RANGE(t2lbl);
00202     
00203     
00204     
00205     
00206 
00207 #ifdef LINK_DOUBLETHRESH_DEBUG
00208     printf("Thresh1 has %d fragments, Thresh2 has %d fragments\n",(int)rnd(t1lbl->gmax),(int)rnd(t2lbl->gmax));
00209 #endif    
00210     
00211     /* Allocate memory */
00212     if (!(t1_data=(SEGMENT_DATA *)malloc((int)rnd(t1lbl->gmax)*sizeof(SEGMENT_DATA)))) { 
00213         MESSAGE('E'," Could not allocate memory for t1_data list.");
00214         goto the_end;
00215     }
00216     if (!(t2_data=(SEGMENT_DATA *)malloc((int)rnd(t2lbl->gmax)*sizeof(SEGMENT_DATA)))) { 
00217         MESSAGE('E'," Could not allocate memory for t2_data list.");
00218         goto the_end;
00219     }
00220     if (!(t2include=(int *)malloc((int)rnd(t2lbl->gmax)*sizeof(int)))) { 
00221         MESSAGE('E'," Could not allocate memory for t2include list.");
00222         goto the_end;
00223     }
00224     
00225     
00226 
00227     /* Initialize the details of each ring fragment */
00228     for (curt1lbl = (int)rnd(t1lbl->gmin); curt1lbl < (int)rnd(t1lbl->gmax); curt1lbl++) {
00229         t1_data[curt1lbl].topx = 0;
00230         t1_data[curt1lbl].topy = t1lbl->nlin-1;
00231         t1_data[curt1lbl].bottomx = 0;
00232         t1_data[curt1lbl].bottomy = 0;
00233         t1_data[curt1lbl].toplinkdone = 0;
00234         t1_data[curt1lbl].bottomlinkdone = 0;
00235     }
00236     for (curt2lbl = (int)rnd(t2lbl->gmin); curt2lbl < (int)rnd(t2lbl->gmax); curt2lbl++) {
00237         t2_data[curt2lbl].topx = 0;
00238         t2_data[curt2lbl].topy = t2lbl->nlin-1;
00239         t2_data[curt2lbl].bottomx = 0;
00240         t2_data[curt2lbl].bottomy = 0;
00241         t2_data[curt2lbl].toplinkdone = 0;
00242         t2_data[curt2lbl].bottomlinkdone = 0;
00243         t2include[curt2lbl] = 0;
00244     }
00245 
00246     /* Find the details of each ring fragment */
00247     for (j=0; j<thresh1->nlin; j++) {
00248         for (k=0; k<thresh1->npix; k++) {
00249             if (t1lbl->data[0][j][k] > (PIXEL)0.0) {
00250                 curt1lbl = (int)rnd(t1lbl->data[0][j][k]) - 1;
00251                 if (j<t1_data[curt1lbl].topy) {
00252                     t1_data[curt1lbl].topy = j;
00253                     t1_data[curt1lbl].topx = k;
00254                 }
00255                 if (j>t1_data[curt1lbl].bottomy) {
00256                     t1_data[curt1lbl].bottomy = j;
00257                     t1_data[curt1lbl].bottomx = k;
00258                 }
00259             }
00260 
00261             if (t2lbl->data[0][j][k] > (PIXEL)0.0) {
00262                 curt2lbl = (int)rnd(t2lbl->data[0][j][k]) - 1;
00263                 if (j<t2_data[curt2lbl].topy) {
00264                     t2_data[curt2lbl].topy = j;
00265                     t2_data[curt2lbl].topx = k;
00266                 }
00267                 if (j>t2_data[curt2lbl].bottomy) {
00268                     t2_data[curt2lbl].bottomy = j;
00269                     t2_data[curt2lbl].bottomx = k;
00270                 }
00271             }
00272         }
00273     }
00274     
00275 #ifdef LINK_DOUBLETHRESH_DEBUG
00276     //printf("\n Details for Thresh1:\n");
00277     //for (curt1lbl = (int)rnd(t1lbl->gmin); curt1lbl < (int)rnd(t1lbl->gmax); curt1lbl++) {
00278     //    printf("T1-%d: top(%d,%d), bottom(%d,%d)\n",curt1lbl+1,t1_data[curt1lbl].topx,
00279     //        t1_data[curt1lbl].topy,t1_data[curt1lbl].bottomx,t1_data[curt1lbl].bottomy);
00280     //}
00281     //printf("\n");
00282 #endif    
00283     
00284 
00285     /* create output images of appropriate size */
00286     if (!CHECKIMG(*out)) GETMEM(1,thresh1->nlin,thresh1->npix,out);
00287     if (!*out) goto the_end;
00288     if (!CHECKIMG(*fillgaps)) GETMEM(1,thresh1->nlin,thresh1->npix,fillgaps);
00289     if (!*fillgaps) goto the_end;
00290 
00291     /* Initialize the output images */
00292     for (j=0; j<thresh1->nlin; j++) {
00293         for (k=0; k<thresh1->npix; k++) {
00294             if (t1lbl->data[0][j][k] > (PIXEL)0.0) {
00295                 (*out)->data[0][j][k] = (PIXEL)1.0;
00296             } else {
00297                 (*out)->data[0][j][k] = (PIXEL)0.0;
00298             }
00299             (*fillgaps)->data[0][j][k] = (PIXEL)0.0;
00300         }
00301     }
00302 
00303 
00304 
00305 
00306     
00307     
00308     /* Perform linking */
00309     for (curt1lbl = (int)rnd(t1lbl->gmin); curt1lbl < (int)rnd(t1lbl->gmax); curt1lbl++) {
00310         
00311         /***************/
00312         /* Link upward */
00313         /***************/
00314         
00315         done = 0;
00316         while (!done) {
00317         
00318             /* Initialize memory to track neighboring segments */
00319             for (i=0; i<nbrhoodsize; i++) {
00320                 t1upnbrchk[i][0] = 0;
00321                 t1upnbrchk[i][1] = 0;
00322                 t2upnbrchk[i][0] = 0;
00323                 t2upnbrchk[i][1] = 0;
00324             }
00325             t1upnbrcnt = 0;
00326             t2upnbrcnt = 0;
00327             connect = 0;
00328 
00329             curt2lbl = (int)rnd(t2lbl->data[0][t1_data[curt1lbl].topy][t1_data[curt1lbl].topx])-1;
00330         
00331             if (t1_data[curt1lbl].topy == 0) {
00332             
00333 #ifdef LINK_DOUBLETHRESH_DEBUG
00334                 printf("Exiting T1-%d uplink -- connected to top of ROI\n",curt1lbl+1);
00335 #endif
00336 
00337                 /* This segment is already connected to the top of the ROI, so end here */
00338                 done = 1;
00339                 t1_data[curt1lbl].toplinkdone = 1;
00340                 
00341                 if (curt2lbl >= 0) {
00342                     t2_data[curt2lbl].toplinkdone = 1;
00343                 }
00344 
00345             } else if (t1_data[curt1lbl].toplinkdone == 1) {
00346 
00347 #ifdef LINK_DOUBLETHRESH_DEBUG
00348                 printf("Exiting T1-%d uplink -- already uplinked\n",curt1lbl+1);
00349 #endif
00350 
00351                 /* This segment was previously uplinked, so end here */
00352                 done = 1;
00353 
00354                 if (curt2lbl >= 0) {
00355                     t2_data[curt2lbl].toplinkdone = 1;
00356                 }
00357 
00358             } else if ((curt2lbl >= 0)&&(t2_data[curt2lbl].toplinkdone == 1)) {
00359 
00360 #ifdef LINK_DOUBLETHRESH_DEBUG
00361                 printf("Exiting T1-%d uplink -- already uplinked\n",curt1lbl+1);
00362 #endif
00363 
00364                 /* This segment was previously uplinked, so end here */
00365                 done = 1;
00366                 t1_data[curt1lbl].toplinkdone = 1;
00367 
00368             } else {
00369             
00370                 /* Find limits of neighborhood */
00371                 if (t1_data[curt1lbl].topy >= maxlinkdist) {
00372                     min_y = t1_data[curt1lbl].topy - maxlinkdist;
00373                 } else {
00374                     min_y = 0;
00375                 }
00376 
00377                 max_y = t1_data[curt1lbl].topy;
00378 
00379                 if (t1_data[curt1lbl].topx >= maxlinkdist) {
00380                     min_x = t1_data[curt1lbl].topx - maxlinkdist;
00381                 } else {
00382                     min_x = 0;
00383                 }
00384 
00385                 if (t1_data[curt1lbl].topx <= t1lbl->npix-(maxlinkdist+1)) {
00386                     max_x = t1_data[curt1lbl].topx + maxlinkdist;
00387                 } else {
00388                     max_x = t1lbl->npix-1;
00389                 }
00390 
00391                 /* Find fragments in upper neighborhood */
00392                 for (j=min_y; j<=max_y; j++) {
00393                     for (k=min_x; k<=max_x; k++) {
00394 
00395                         /* Don't include the endpoint of the current segment */
00396                         if ((j != t1_data[curt1lbl].topy)||(k != t1_data[curt1lbl].topx)) {
00397                             if ((t1lbl->data[0][j][k] > (PIXEL)0.0)&&(t1lbl->data[0][j][k] != (PIXEL)(curt1lbl+1))) {
00398                                 t1upnbrchk[t1upnbrcnt][0] = (int)rnd(t1lbl->data[0][j][k]);
00399 
00400                                 if (t1_data[t1upnbrchk[t1upnbrcnt][0]-1].bottomy > t1_data[curt1lbl].topy) {
00401                                     dist_y = t1_data[t1upnbrchk[t1upnbrcnt][0]-1].bottomy - t1_data[curt1lbl].topy;
00402                                 } else {
00403                                     dist_y = t1_data[curt1lbl].topy - t1_data[t1upnbrchk[t1upnbrcnt][0]-1].bottomy;
00404                                 }
00405 
00406                                 if (t1_data[t1upnbrchk[t1upnbrcnt][0]-1].bottomx > t1_data[curt1lbl].topx) {
00407                                     dist_x = t1_data[t1upnbrchk[t1upnbrcnt][0]-1].bottomx - t1_data[curt1lbl].topx;
00408                                 } else {
00409                                     dist_x = t1_data[curt1lbl].topx - t1_data[t1upnbrchk[t1upnbrcnt][0]-1].bottomx;
00410                                 }
00411 
00412                                 t1upnbrchk[t1upnbrcnt][1] = dist_x + dist_y;
00413                                 t1upnbrcnt++;
00414                             }
00415 
00416                             if (t2lbl->data[0][j][k] > (PIXEL)0.0) {
00417                                 t2upnbrchk[t2upnbrcnt][0] = (int)rnd(t2lbl->data[0][j][k]);
00418 
00419                                 if (t2_data[t2upnbrchk[t2upnbrcnt][0]-1].bottomy > t1_data[curt1lbl].topy) {
00420                                     dist_y = t2_data[t2upnbrchk[t2upnbrcnt][0]-1].bottomy - t1_data[curt1lbl].topy;
00421                                 } else {
00422                                     dist_y = t1_data[curt1lbl].topy - t2_data[t2upnbrchk[t2upnbrcnt][0]-1].bottomy;
00423                                 }
00424 
00425                                 if (t2_data[t2upnbrchk[t2upnbrcnt][0]-1].bottomx > t1_data[curt1lbl].topx) {
00426                                     dist_x = t2_data[t2upnbrchk[t2upnbrcnt][0]-1].bottomx - t1_data[curt1lbl].topx;
00427                                 } else {
00428                                     dist_x = t1_data[curt1lbl].topx - t2_data[t2upnbrchk[t2upnbrcnt][0]-1].bottomx;
00429                                 }
00430 
00431                                 t2upnbrchk[t2upnbrcnt][1] = dist_x + dist_y;
00432                                 t2upnbrcnt++;
00433                             }
00434                         }
00435 
00436                     }
00437                 }
00438 
00439                 if (t1upnbrcnt > 1) {
00440                     /* Sort the neighboring fragments */
00441                     quicksort((int *)t1upnbrchk,0,t1upnbrcnt-1);
00442                 }
00443                 if (t2upnbrcnt > 1) {
00444                     /* Sort the neighboring fragments */
00445                     quicksort((int *)t2upnbrchk,0,t2upnbrcnt-1);
00446                 }
00447 
00448 
00449 #ifdef LINK_DOUBLETHRESH_DEBUG
00450                 //printf("\nt2upnbrcnt = %d\n\n",t2upnbrcnt);
00451 #endif
00452 
00453 
00454 
00455                 if (t1upnbrcnt && !t2upnbrcnt) {
00456 
00457                     /* If the T1 neighbor is reasonable, connect to it */
00458                     closet1lbl = t1upnbrchk[t1upnbrcnt-1][0] - 1;
00459 
00460                     /* Check if the closest segment has already been connected downward */
00461                     if (t1_data[closet1lbl].bottomlinkdone) {
00462                         
00463 #ifdef LINK_DOUBLETHRESH_DEBUG
00464                         printf("Exiting T1-%d uplink -- closest segment connected downward\n",curt1lbl+1);
00465 #endif
00466 
00467                         /* the closest segment is already linked downward, so quit here */
00468                         done = 1;
00469                         t1_data[curt1lbl].toplinkdone = 1;
00470 
00471                         if (curt2lbl >= 0) {
00472                             t2_data[curt2lbl].toplinkdone = 1;
00473                         }
00474 
00475                     } else {
00476                     
00477 #ifdef LINK_DOUBLETHRESH_DEBUG
00478                         printf("Up-Linking T1-%d to T1-%d\n",curt1lbl+1,closet1lbl+1);
00479 #endif
00480 
00481                         startx = t1_data[curt1lbl].topx;
00482                         starty = t1_data[curt1lbl].topy;
00483                         endx = t1_data[closet1lbl].bottomx;
00484                         endy = t1_data[closet1lbl].bottomy;
00485                         connect = 1;
00486                         done = 1;
00487 
00488                         t1_data[curt1lbl].topx = t1_data[closet1lbl].topx;
00489                         t1_data[curt1lbl].topy = t1_data[closet1lbl].topy;
00490                         t1_data[closet1lbl].bottomx = t1_data[curt1lbl].bottomx;
00491                         t1_data[closet1lbl].bottomy = t1_data[curt1lbl].bottomy;
00492                         t1_data[closet1lbl].bottomlinkdone = 1;
00493                         t1_data[curt1lbl].toplinkdone = 1;
00494 
00495                         if (curt2lbl >= 0) {
00496                             t2_data[curt2lbl].toplinkdone = 1;
00497                         }
00498 
00499                     }
00500 
00501 
00502                 } else if (t1upnbrcnt && t2upnbrcnt) {
00503                     
00504                     closet1lbl = t1upnbrchk[t1upnbrcnt-1][0] - 1;
00505                     closet2lbl = t2upnbrchk[t2upnbrcnt-1][0] - 1;
00506 
00507 
00508                     if (!t1_data[closet1lbl].bottomlinkdone) {
00509                         
00510                         if (t1upnbrchk[t1upnbrcnt-1][1] <= t2upnbrchk[t2upnbrcnt-1][1]) {
00511                             /* Connect to t1 neighbor, since it is the closest */
00512 
00513 #ifdef LINK_DOUBLETHRESH_DEBUG
00514                             printf("Up-Linking T1-%d to T1-%d\n",curt1lbl+1,closet1lbl+1);
00515 #endif
00516 
00517                             startx = t1_data[curt1lbl].topx;
00518                             starty = t1_data[curt1lbl].topy;
00519                             endx = t1_data[closet1lbl].bottomx;
00520                             endy = t1_data[closet1lbl].bottomy;
00521                             connect = 1;
00522                             done = 1;
00523 
00524                             t1_data[curt1lbl].topx = t1_data[closet1lbl].topx;
00525                             t1_data[curt1lbl].topy = t1_data[closet1lbl].topy;
00526                             t1_data[closet1lbl].bottomx = t1_data[curt1lbl].bottomx;
00527                             t1_data[closet1lbl].bottomy = t1_data[curt1lbl].bottomy;
00528                             t1_data[closet1lbl].bottomlinkdone = 1;
00529                             t1_data[curt1lbl].toplinkdone = 1;
00530 
00531                             if (curt2lbl >= 0) {
00532                                 t2_data[curt2lbl].toplinkdone = 1;
00533                             }
00534                             
00535                         } else {
00536                             /* Check if curt1lbl is the closest label in t1 to closet1lbl */
00537 
00538                             /* Initialize memory to track neighboring segments */
00539                             for (i=0; i<nbrhoodsize; i++) {
00540                                 t1downnbrchk[i][0] = 0;
00541                                 t1downnbrchk[i][1] = 0;
00542                             }
00543                             t1downnbrcnt = 0;
00544 
00545                             /* Find limits of neighborhood */
00546                             min_y = t1_data[closet1lbl].bottomy;
00547 
00548                             if (t1_data[closet1lbl].bottomy <= t1lbl->nlin-(maxlinkdist+1)) {
00549                                 max_y = t1_data[closet1lbl].bottomy + maxlinkdist;
00550                             } else {
00551                                 min_y = t1lbl->nlin-1;
00552                             }
00553 
00554                             if (t1_data[closet1lbl].bottomx >= maxlinkdist) {
00555                                 min_x = t1_data[closet1lbl].bottomx - maxlinkdist;
00556                             } else {
00557                                 min_x = 0;
00558                             }
00559 
00560                             if (t1_data[closet1lbl].bottomx <= t1lbl->npix-(maxlinkdist+1)) {
00561                                 max_x = t1_data[closet1lbl].bottomx + maxlinkdist;
00562                             } else {
00563                                 max_x = t1lbl->npix-1;
00564                             }
00565 
00566                             /* Find fragments in lower neighborhood */
00567                             for (j=min_y; j<=max_y; j++) {
00568                                 for (k=min_x; k<=max_x; k++) {
00569 
00570                                     /* Don't include the endpoint of the current segment */
00571                                     if ((j != t1_data[closet1lbl].bottomy)||(k != t1_data[closet1lbl].bottomx)) {
00572                                         if ((t1lbl->data[0][j][k] > (PIXEL)0.0)&&(t1lbl->data[0][j][k] != (PIXEL)(closet1lbl+1))) {
00573                                             t1downnbrchk[t1downnbrcnt][0] = (int)rnd(t1lbl->data[0][j][k]);
00574 
00575                                             if (t1_data[t1downnbrchk[t1downnbrcnt][0]-1].topy > t1_data[closet1lbl].bottomy) {
00576                                                 dist_y = t1_data[t1downnbrchk[t1downnbrcnt][0]-1].topy - t1_data[closet1lbl].bottomy;
00577                                             } else {
00578                                                 dist_y = t1_data[closet1lbl].bottomy - t1_data[t1downnbrchk[t1downnbrcnt][0]-1].topy;
00579                                             }
00580 
00581                                             if (t1_data[t1downnbrchk[t1downnbrcnt][0]-1].topx > t1_data[closet1lbl].bottomx) {
00582                                                 dist_x = t1_data[t1downnbrchk[t1downnbrcnt][0]-1].topx - t1_data[closet1lbl].bottomx;
00583                                             } else {
00584                                                 dist_x = t1_data[closet1lbl].bottomx - t1_data[t1downnbrchk[t1downnbrcnt][0]-1].topx;
00585                                             }
00586 
00587                                             t1downnbrchk[t1downnbrcnt][1] = dist_x + dist_y;
00588                                             t1downnbrcnt++;
00589                                         }
00590                                     }
00591 
00592                                 }
00593                             }
00594 
00595                             if (t1downnbrcnt > 1) {
00596                                 /* Sort the neighboring fragments */
00597                                 quicksort((int *)t1downnbrchk,0,t1downnbrcnt-1);
00598                             }
00599                             
00600                             
00601                             if ((t1downnbrcnt > 0)&&(t1upnbrchk[t1upnbrcnt-1][1] <= t1downnbrchk[t1downnbrcnt-1][1])) {
00602 #ifdef LINK_DOUBLETHRESH_DEBUG
00603                                 printf("Up-Linking T1-%d to T1-%d\n",curt1lbl+1,closet1lbl+1);
00604 #endif
00605 
00606                                 startx = t1_data[curt1lbl].topx;
00607                                 starty = t1_data[curt1lbl].topy;
00608                                 endx = t1_data[closet1lbl].bottomx;
00609                                 endy = t1_data[closet1lbl].bottomy;
00610                                 connect = 1;
00611                                 done = 1;
00612 
00613                                 t1_data[curt1lbl].topx = t1_data[closet1lbl].topx;
00614                                 t1_data[curt1lbl].topy = t1_data[closet1lbl].topy;
00615                                 t1_data[closet1lbl].bottomx = t1_data[curt1lbl].bottomx;
00616                                 t1_data[closet1lbl].bottomy = t1_data[curt1lbl].bottomy;
00617                                 t1_data[closet1lbl].bottomlinkdone = 1;
00618                                 t1_data[curt1lbl].toplinkdone = 1;
00619 
00620                                 if (curt2lbl >= 0) {
00621                                     t2_data[curt2lbl].toplinkdone = 1;
00622                                 }
00623                             }
00624 
00625                         }
00626                     }
00627 
00628                     
00629                     if (!connect) {
00630                     
00631                         /* Check if the closest t2 segment has already been connected downward */
00632                         if (t2_data[closet2lbl].bottomlinkdone) {
00633 
00634 #ifdef LINK_DOUBLETHRESH_DEBUG
00635                             printf("Exiting T1-%d uplink -- closest segment connected downward\n",curt1lbl+1);
00636 #endif
00637 
00638                             /* the closest segment is already linked downward, so quit here */
00639                             done = 1;
00640                             t1_data[curt1lbl].toplinkdone = 1;
00641 
00642                             if (curt2lbl >= 0) {
00643                                 t2_data[curt2lbl].toplinkdone = 1;
00644                             }
00645 
00646                         } else {
00647 
00648 #ifdef LINK_DOUBLETHRESH_DEBUG
00649                             printf("Up-Linking T1-%d to T2-%d\n",curt1lbl+1,closet2lbl+1);
00650 #endif
00651 
00652                             startx = t1_data[curt1lbl].topx;
00653                             starty = t1_data[curt1lbl].topy;
00654                             endx = t2_data[closet2lbl].bottomx;
00655                             endy = t2_data[closet2lbl].bottomy;
00656                             connect = 1;
00657 
00658                             t1_data[curt1lbl].topx = t2_data[closet2lbl].topx;
00659                             t1_data[curt1lbl].topy = t2_data[closet2lbl].topy;
00660                             t2_data[closet2lbl].bottomx = t1_data[curt1lbl].bottomx;
00661                             t2_data[closet2lbl].bottomy = t1_data[curt1lbl].bottomy;
00662                             t2_data[closet2lbl].bottomlinkdone = 1;
00663                             t2include[closet2lbl] = 1;
00664 
00665                             if (curt2lbl >= 0) {
00666                                 t2_data[curt2lbl].toplinkdone = 1;
00667                             }
00668 
00669                         }
00670 
00671                     }
00672                     
00673                 } else if (!t1upnbrcnt && t2upnbrcnt) {
00674                     
00675                     closet2lbl = t2upnbrchk[t2upnbrcnt-1][0] - 1;
00676 
00677 #ifdef LINK_DOUBLETHRESH_DEBUG
00678                     printf("\n\nUpper Neighbors of T1-%d:\n",curt1lbl+1);
00679                     for (i=0; i<t2upnbrcnt; i++) {
00680                         printf("%d - %d pixels away\n",t2upnbrchk[i][0],t2upnbrchk[i][1]);
00681                     }
00682                     printf("\n");
00683 #endif
00684                     
00685 
00686                     /* Check if the closest segment has already been connected downward */
00687                     if (t2_data[closet2lbl].bottomlinkdone) {
00688                         
00689 #ifdef LINK_DOUBLETHRESH_DEBUG
00690                         printf("Exiting T1-%d uplink -- closest segment connected downward\n",curt1lbl+1);
00691 #endif
00692 
00693                         /* the closest segment is already linked downward, so quit here */
00694                         done = 1;
00695                         t1_data[curt1lbl].toplinkdone = 1;
00696 
00697                         if (curt2lbl >= 0) {
00698                             t2_data[curt2lbl].toplinkdone = 1;
00699                         }
00700 
00701                     } else {
00702                     
00703 #ifdef LINK_DOUBLETHRESH_DEBUG
00704                         printf("Up-Linking T1-%d to T2-%d\n",curt1lbl+1,closet2lbl+1);
00705 #endif
00706 
00707                         startx = t1_data[curt1lbl].topx;
00708                         starty = t1_data[curt1lbl].topy;
00709                         endx = t2_data[closet2lbl].bottomx;
00710                         endy = t2_data[closet2lbl].bottomy;
00711                         connect = 1;
00712 
00713                         t1_data[curt1lbl].topx = t2_data[closet2lbl].topx;
00714                         t1_data[curt1lbl].topy = t2_data[closet2lbl].topy;
00715                         t2_data[closet2lbl].bottomx = t1_data[curt1lbl].bottomx;
00716                         t2_data[closet2lbl].bottomy = t1_data[curt1lbl].bottomy;
00717                         t2_data[closet2lbl].bottomlinkdone = 1;
00718                         t2include[closet2lbl] = 1;
00719 
00720                         if (curt2lbl >= 0) {
00721                             t2_data[curt2lbl].toplinkdone = 1;
00722                         }
00723 
00724                     }
00725 
00726                 } else {
00727 
00728 #ifdef LINK_DOUBLETHRESH_DEBUG
00729                     printf("Exiting T1-%d uplink -- no fragments in upper neighborhood\n",curt1lbl+1);
00730 #endif
00731 
00732                     /* There are no fragments in the upper neighborhood, so end here */
00733                     done = 1;
00734                     t1_data[curt1lbl].toplinkdone = 1;
00735 
00736                     if (curt2lbl >= 0) {
00737                         t2_data[curt2lbl].toplinkdone = 1;
00738                     }
00739 
00740                 }
00741                 
00742                 
00743 
00744                 if (connect) {
00745                     /* Count number of points in line */
00746                     count_point = TRUE;
00747                     point = 0;
00748                     brshnm(startx,starty,endx, endy);
00749 
00750                     /* Allocate memory */
00751                     x_array_buffer = (int **) malloc(point*sizeof(int  * ));
00752                     if (x_array_buffer == (int **) NULL){
00753                         printf("No memory for x_array_buffer\n");
00754                     } else{
00755 
00756                         for  (l = 0; l < point; l++){
00757                             x_array_buffer[l] = (int  *) calloc(2, sizeof(int));
00758                             if (x_array_buffer[l] == (int) NULL){
00759                                 printf ("No memory for element %d\n",x_array_buffer[l]);
00760                             }
00761                         }
00762 
00763                         /* Fill array */ 
00764                         count_point = FALSE;
00765                         point = 0;
00766                         brshnm(startx,starty,endx, endy);
00767 
00768 #ifdef LINK_DOUBLETHRESH_DEBUG
00769                         //printf("startx=%d, starty=%d, endx=%d, endy=%d\n",startx,starty,endx,endy);
00770                         //printf("     x=%d,      y=%d,    x=%d,    y=%d\n",x_array_buffer[0][0],x_array_buffer[0][1],x_array_buffer[point-1][0],x_array_buffer[point-1][1]);
00771 #endif
00772 
00773                             for (l=0; l < point; l++) {
00774 #ifdef LINK_DOUBLETHRESH_DEBUG
00775                                  //printf("filling (%d,%d)\n",x_array_buffer[l][0],x_array_buffer[l][1]);
00776 #endif
00777                                  (*fillgaps)->data[0][x_array_buffer[l][1]][x_array_buffer[l][0]] = (PIXEL)1.0;
00778                             }
00779 
00780                         if(x_array_buffer) {
00781                             /* Free Bresh Array */
00782                             for (l=0; l < point; l++){
00783                                 if(x_array_buffer[l]) {
00784                                     free(x_array_buffer[l]);
00785                                     x_array_buffer[l] = NULL;
00786                                 }
00787                             }
00788 
00789                             free(x_array_buffer);
00790                             x_array_buffer = NULL;
00791                         }
00792 
00793                     }
00794 
00795                 }
00796 
00797                 
00798             }
00799 
00800         }
00801     }
00802     
00803     
00804 
00805     /* Include appropriate segments from thresh2 in output */
00806     for (j=0; j<thresh2->nlin; j++) {
00807         for (k=0; k<thresh2->npix; k++) {
00808             if (t2lbl->data[0][j][k] > (PIXEL)0.0) {
00809                 if (t2include[(int)rnd(t2lbl->data[0][j][k]-1)] > 0) {
00810                     (*out)->data[0][j][k] = (PIXEL)1.0;
00811                 }
00812             }
00813         }
00814     }
00815 
00816     the_end:
00817     /* if (TIMES) TIMING(T_EXIT); */
00818     
00819     if (CHECKIMG(t1))  RELIMG(&t1);
00820     if (CHECKIMG(t2))  RELIMG(&t2);
00821     if (CHECKIMG(t1lbl))  RELIMG(&t1lbl);
00822     if (CHECKIMG(t2lbl))  RELIMG(&t2lbl);
00823     if (t1_data) free(t1_data);
00824     if (t2_data) free(t2_data);
00825     if (t2include) free(t2include);
00826     if (t1include) free(t1include);
00827 }