label.c Source File

00001 #include        "sadie.h"
00002 
00003 /*-Copyright Information------------------------------------------------------*/
00004 /* Copyright (c) 1990 by the University of Arizona Digital Image Analysis Lab */
00005 /*----------------------------------------------------------------------------*/
00006 /*-General Information--------------------------------------------------------*/
00007 /*                                                                            */
00008 /*   This function segments an image into regions of similar graylevels       */
00009 /*   and creates a label map. Two regions with mean values mean1 and mean2,   */
00010 /*   and a combined standard deviation dev, are merged if                     */
00011 /*                                                                            */
00012 /*           abs(mean1 - mean2)  <=  thresh * (1 - dev/sigma)                 */
00013 /*                                                                            */
00014 /*   where thresh and sigma are input parameters. Assuming the image is       */
00015 /*   segmented into N regions, the resulting label map consists of N          */
00016 /*   segments, numbered consecutively from 0 to N-1.                          */
00017 /*                                                                            */
00018 /*----------------------------------------------------------------------------*/
00019 /*-Background Information-----------------------------------------------------*/
00020 /*                                                                            */
00021 /*   Mehldau, G.:                                                             */
00022 /*   "An Image Segmentation Algorithm."                                       */
00023 /*   University of Arizona, DIAL, Technical Report No. 90-3                   */
00024 /*                                                                            */
00025 /*----------------------------------------------------------------------------*/
00026 /*-Interface Information------------------------------------------------------*/
00027 void LABEL (
00028 IMAGE  *in,     /*  I   Pointer to the input image.                           */
00029 double sigma,   /*  I   Greylevel variance threshold.                         */
00030 PIXEL  thresh,  /*  I   Graylevel difference threshold.                       */
00031 IMAGE  **out    /*  O   Address of a pointer to the output image.             */
00032 /*----------------------------------------------------------------------------*/
00033 ) { register short t, l, i, j, k;
00034     register double p, q, r, s, topn, lftn, topsum, lftsum, pinc, psum;
00035     char   msg[SLEN];
00036     long   total, *cnt=0, *eqv=0, top, lft, newLabel, label=0L;
00037     float  *mean=0, *dev=0;
00038     double *topmean=0, *topdev=0, *lftmean=0, *lftdev=0;
00039 
00040     if (TIMES) TIMING(T_LABEL);
00041     if (NAMES) {
00042         MESSAGE('I',"");
00043         MESSAGE('I',"LABEL");
00044         MESSAGE('I',"");
00045         sprintf(msg," Input image:                      %s",in->text);
00046         MESSAGE('I',msg);
00047         sprintf(msg," Graylevel difference threshold: %12.4e",thresh);
00048         MESSAGE('I',msg);
00049         sprintf(msg," Graylevel variance threshold:   %12.4e",sigma);
00050         MESSAGE('I',msg);
00051         MESSAGE('I'," ...............");
00052     }
00053 
00054     /* check input */
00055     if (!CHECKIMG(in)) {
00056         MESSAGE('E'," Can't identify image.");
00057         goto the_end;
00058     } else if (sigma < 0.0) {
00059         MESSAGE('E'," Threshold  must be greater than or equal to zero.");
00060         goto the_end;
00061     } else if (thresh < 0.0) {
00062         MESSAGE('E'," Threshold  must be greater than or equal to zero.");
00063         goto the_end;
00064     }
00065 
00066     /* allocate space */
00067     total = (long)in->nlin*(long)in->npix;
00068     if (!(cnt=(long *)malloc(total*sizeof(long)))) {
00069        MESSAGE('E'," Memory request failed.");
00070        goto the_end;
00071     }
00072     if (!(eqv=(long *)malloc(total*sizeof(long)))) {
00073        MESSAGE('E'," Memory request failed.");
00074        goto the_end;
00075     }
00076     total *= (long)in->nbnd;
00077     if (!(mean=(float *)malloc(total*sizeof(float)))) {
00078        MESSAGE('E'," Memory request failed.");
00079        goto the_end;
00080     }
00081     if (!(dev=(float *)malloc(total*sizeof(float)))) {
00082        MESSAGE('E'," Memory request failed.");
00083        goto the_end;
00084     }
00085     if (!(topmean=(double *)malloc(in->nbnd*sizeof(double)))) {
00086        MESSAGE('E'," Memory request failed.");
00087        goto the_end;
00088     }
00089     if (!(topdev=(double *)malloc(in->nbnd*sizeof(double)))) {
00090        MESSAGE('E'," Memory request failed.");
00091        goto the_end;
00092     }
00093     if (!(lftmean=(double *)malloc(in->nbnd*sizeof(double)))) {
00094        MESSAGE('E'," Memory request failed.");
00095        goto the_end;
00096     }
00097     if (!(lftdev=(double *)malloc(in->nbnd*sizeof(double)))) {
00098        MESSAGE('E'," Memory request failed.");
00099        goto the_end;
00100     }
00101 
00102     /* create image of appropriate size */
00103     if (!CHECKIMG(*out)) GETMEM(1,in->nlin,in->npix,out);
00104     if (!*out) goto the_end;
00105 
00106     /* initialize progress indicator */
00107     if (LINES  &&  !PROGRESS(psum=0.0)) goto the_end;
00108     pinc = 1.0/(double)in->nlin;
00109 
00110     for (j=0; j<in->nlin; j++) {
00111         for (k=0; k<in->npix; k++) {
00112 
00113             /* look at adjacent segments */
00114             if (t = j) {
00115                 for (top=eqv[(long)(*out)->data[0][j-1][k]]; top!=eqv[top]; top=eqv[top]);
00116                 eqv[(long)(*out)->data[0][j-1][k]] = top;
00117                 topn = (double)cnt[top];
00118                 top *= in->nbnd;
00119                 topsum = 0.0;
00120             }
00121             if (l = k) {
00122                 for (lft=eqv[(long)(*out)->data[0][j][k-1]]; lft!=eqv[lft]; lft=eqv[lft]);
00123                 eqv[(long)(*out)->data[0][j][k-1]] = lft;
00124                 lftn = (double)cnt[lft];
00125                 lft *= in->nbnd;
00126                 lftsum = 0.0;
00127             }
00128 
00129             /* merge pixel with adjacent segment? */
00130             for (i=0; (t || l)  &&  i<in->nbnd; i++) {
00131                 p = in->data[i][j][k];
00132                 if (t) {
00133                     q = mean[top+i];
00134                     r = dev[top+i];
00135                     topmean[i] = (p + topn*q)/(topn+1.0);
00136                     topdev[i]  = sqrt(topn*(r*r + (p-q)*(p-q)/(topn+1.0))/(topn+1.0));
00137                     t = fabs(p-q) <= thresh*(1.0 - topdev[i]/sigma);
00138                     topsum += (p-q)*(p-q);
00139                 }
00140                 if (l) {
00141                     q = mean[lft+i];
00142                     r = dev[lft+i];
00143                     lftmean[i] = (p + lftn*q)/(lftn+1.0);
00144                     lftdev[i]  = sqrt(lftn*(r*r + (p-q)*(p-q)/(lftn+1.0))/(lftn+1.0));
00145                     l = fabs(p-q) <= thresh*(1.0 - lftdev[i]/sigma);
00146                     lftsum += (p-q)*(p-q);
00147                 }
00148             }
00149 
00150             /* assign pixel label */
00151             if (t && l) {
00152                 l = !(t = topsum < lftsum);
00153             }
00154             if (t) {
00155                 for (i=0; i<in->nbnd; i++) {
00156                     mean[top+i] = topmean[i];
00157                     dev[top+i]  = topdev[i];
00158                 }
00159                 cnt[top/in->nbnd] += 1L;
00160                 (*out)->data[0][j][k] = (*out)->data[0][j-1][k];
00161             } else if (l) {
00162                 for (i=0; i<in->nbnd; i++) {
00163                     mean[lft+i] = lftmean[i];
00164                     dev[lft+i]  = lftdev[i];
00165                 }
00166                 cnt[lft/in->nbnd] += 1L;
00167                 (*out)->data[0][j][k] = (*out)->data[0][j][k-1];
00168             } else {
00169                 newLabel = label*in->nbnd;
00170                 for (i=0; i<in->nbnd; i++) {
00171                     mean[newLabel+i] = in->data[i][j][k];
00172                     dev[newLabel+i]  = 0.0;
00173                 }
00174                 cnt[label] = 1L;
00175                 eqv[label] = label;
00176                 (*out)->data[0][j][k] = (PIXEL)label++;
00177             }
00178 
00179             /* merge adjacent segments? */
00180             if ((t  ||  l)   &&  top != lft  &&  j  &&  k) {
00181                 for (t=TRUE,i=0; t && i<in->nbnd; i++) {
00182                     p = mean[top+i];
00183                     q = mean[lft+i];
00184                     r = dev[top+i];
00185                     s = dev[lft+i];
00186                     topmean[i] = (topn*p + lftn*q)/(topn+lftn);
00187                     topdev[i]  = sqrt((topn*(r*r + p*p) + lftn*(q*q + s*s))/(topn+lftn) - topmean[i]*topmean[i]);
00188                     t = fabs(p-q) <= thresh*(1.0-topdev[i]/sigma);
00189                 }
00190 
00191                 /* merge adjacent segments */
00192                 if (t) {
00193                     if (top < lft) {
00194                         for (i=0; i<in->nbnd; i++) {
00195                             mean[top+i] = topmean[i];
00196                             dev[top+i]  = topdev[i];
00197                         }
00198                         top /= in->nbnd;
00199                         lft /= in->nbnd;
00200                         cnt[top] += cnt[lft];
00201                         eqv[lft] = top;
00202                     } else {
00203                         for (i=0; i<in->nbnd; i++) {
00204                             mean[lft+i] = topmean[i];
00205                             dev[lft+i]  = topdev[i];
00206                         }
00207                         top /= in->nbnd;
00208                         lft /= in->nbnd;
00209                         cnt[lft] += cnt[top];
00210                         eqv[top] = lft;
00211                     }
00212                 }
00213             }
00214 
00215         }
00216         if (LINES  &&  !PROGRESS(psum+=pinc)) goto the_end;
00217     }
00218 
00219     /* assign consecutive labels */
00220     for (top=0L; top<label; top++) {
00221         if (eqv[top] != top) {
00222             for (newLabel=eqv[top]; newLabel!=eqv[newLabel]; newLabel=eqv[newLabel]);
00223             eqv[top] = newLabel;
00224         }
00225     }
00226     for (newLabel=top=0L; top<label; top++) {
00227         eqv[top] = (eqv[top] == top ? newLabel++ : eqv[eqv[top]]);
00228     }
00229     for (j=0; j<in->nlin; j++) {
00230         for (k=0; k<in->npix; k++) {
00231             (*out)->data[0][j][k] = (PIXEL)eqv[(long)(*out)->data[0][j][k]];
00232         }
00233     }
00234     (*out)->gmin = (PIXEL)0L;
00235     (*out)->gmax = (PIXEL)(newLabel-1L);
00236     MESSAGE('I',"");
00237     sprintf(msg," Number of segments = %ld",newLabel);
00238     MESSAGE('I',msg);
00239 
00240     the_end:
00241     if (cnt)     free(cnt);
00242     if (eqv)     free(eqv);
00243     if (mean)    free(mean);
00244     if (dev)     free(dev);
00245     if (topmean) free(topmean);
00246     if (topdev)  free(topdev);
00247     if (lftmean) free(lftmean);
00248     if (lftdev)  free(lftdev);
00249     if (TIMES)   TIMING(T_EXIT);
00250 }