refine.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 takes as input a label map and two images, where the       */
00009 /*   second image has twice the resolution of the label map and the first     */
00010 /*   image. The output is a "refined" label map with the same resolution      */
00011 /*   as the second image.                                                     */
00012 /*                                                                            */
00013 /*----------------------------------------------------------------------------*/
00014 /*-Background Information-----------------------------------------------------*/
00015 /*                                                                            */
00016 /*   Mehldau, G.:                                                             */
00017 /*   "An Image Segmentation Algorithm."                                       */
00018 /*   University of Arizona, DIAL, Technical Report No. 90-3                   */
00019 /*                                                                            */
00020 /*----------------------------------------------------------------------------*/
00021 /*-Interface Information------------------------------------------------------*/
00022 void REFINE (
00023 IMAGE  *in[3],  /*  I   Array[3] of pointers to the input images.             */
00024                 /*      in[0]:  Reduced label map.                            */
00025                 /*      in[1]:  Reduced image.                                */
00026                 /*      in[2]:  Refined image.                                */
00027 IMAGE  **out    /*  O   Address of a pointer to the output image.             */
00028 /*----------------------------------------------------------------------------*/
00029 ) { register short i, j, k, l, m;
00030     char   msg[SLEN];
00031     double olddiff, newdiff, mean, pinc, psum;
00032 
00033     if (TIMES) TIMING(T_REFINE);
00034     if (NAMES) {
00035         MESSAGE('I',"");
00036         MESSAGE('I',"REFINE");
00037         MESSAGE('I',"");
00038         sprintf(msg," Reduced label map:   %s",in[0]->text);
00039         MESSAGE('I',msg);
00040         sprintf(msg," Reduced image:       %s",in[1]->text);
00041         MESSAGE('I',msg);
00042         sprintf(msg," Refined image:       %s",in[2]->text);
00043         MESSAGE('I',msg);
00044         MESSAGE('I'," ...............");
00045     }
00046 
00047     /* check input */
00048     for (i=0; i<3; i++) {
00049         if (!CHECKIMG(in[i])) {
00050             MESSAGE('E'," Can't identify image.");
00051             goto the_end;
00052         }
00053     }
00054     if (in[0]->npix != in[1]->npix) {
00055         MESSAGE('E'," Number of pixels/line must be identical in reduced image and reduced label map.");
00056         goto the_end;
00057     } else if (in[0]->nlin != in[1]->nlin) {
00058         MESSAGE('E'," Number of lines must be identical in reduced image and reduced label map.");
00059         goto the_end;
00060     } else if (in[0]->nbnd != 1) {
00061         MESSAGE('E'," Reduced label map must be single-band.");
00062         goto the_end;
00063     } else if (2*in[1]->npix != in[2]->npix) {
00064         MESSAGE('E'," Number of pixels/line in reduced image must be half the number of pixels/line in regular image.");
00065         goto the_end;
00066     } else if (2*in[1]->nlin != in[2]->nlin) {
00067         MESSAGE('E'," Number of lines in reduced image must be half the number of lines in regular image.");
00068         goto the_end;
00069     } else if (in[1]->nbnd != in[2]->nbnd) {
00070         MESSAGE('E'," Number of bands must be identical in reduced image and regular image");
00071         goto the_end;
00072     }
00073 
00074     /* create image of appropriate size */
00075     if (!CHECKIMG(*out)) GETMEM(1,in[2]->nlin,in[2]->npix,out);
00076     if (!*out) goto the_end;
00077 
00078     /* initialize progress indicator */
00079     if (LINES  &&  !PROGRESS(psum=0.0)) goto the_end;
00080     pinc = 1.0/(double)(*out)->nlin;
00081   
00082     for (j=0; j<(*out)->nlin; j++) {
00083         for (k=0; k<(*out)->npix; k++) {
00084             switch (j%2 << 1  |  k%2) {
00085 
00086             /* type 1 */
00087             case 0:
00088                 if (j > 0  &&  (*out)->data[0][j-1][k] == in[0]->data[0][j/2][k/2]  &&  (*out)->data[0][j-1][k] != in[0]->data[0][j/2-1][k/2]   ||
00089                     k > 0  &&  (*out)->data[0][j][k-1] == in[0]->data[0][j/2][k/2]  &&  (*out)->data[0][j][k-1] != in[0]->data[0][j/2][k/2-1]   ||   j == 0  &&  k == 0) {
00090                     (*out)->data[0][j][k] = in[0]->data[0][j/2][k/2];
00091                 } else {
00092                     olddiff = -1.0;
00093                     if (j > 0) {
00094                         l = (j-1)/2;
00095                         m = k/2;
00096                         newdiff = 0.0;
00097                         for (i=0; i<in[1]->nbnd; i++) {
00098                             if (l < in[1]->nlin-1  &&  m > 0) {
00099                                 mean = (double)(in[1]->data[i][l+1][m]+in[1]->data[i][l+1][m-1]+in[1]->data[i][l][m]+in[1]->data[i][l][m-1]) / 4.0;
00100                             } else if (l < in[1]->nlin-1) {
00101                                 mean = (double)(in[1]->data[i][l+1][m]+in[1]->data[i][l][m]) / 2.0;
00102                             } else if (m > 0) {
00103                                 mean = (double)(in[1]->data[i][l][m]+in[1]->data[i][l][m-1]) / 2.0;
00104                             } else {
00105                                 mean = (double)(in[1]->data[i][l][m]);
00106                             }
00107                             newdiff += fabs((double)in[2]->data[i][j][k]-mean)*fabs((double)in[2]->data[i][j][k]-mean);
00108                         }
00109                         olddiff = newdiff;
00110                         (*out)->data[0][j][k] = (*out)->data[0][j-1][k];
00111                     }
00112                     if (k > 0) {
00113                         l = j/2;
00114                         m = (k-1)/2;
00115                         newdiff = 0.0;
00116                         for (i=0; i<in[1]->nbnd; i++) {
00117                             if (l > 0  &&  m < in[1]->npix-1) {
00118                                 mean = (double)(in[1]->data[i][l][m+1]+in[1]->data[i][l][m]+in[1]->data[i][l-1][m+1]+in[1]->data[i][l-1][m]) / 4.0;
00119                             } else if (l > 0) {
00120                                 mean = (double)(in[1]->data[i][l][m]+in[1]->data[i][l-1][m]) / 2.0;
00121                             } else if (m < in[1]->npix-1) {
00122                                 mean = (double)(in[1]->data[i][l][m+1]+in[1]->data[i][l][m]) / 2.0;
00123                             } else {
00124                                 mean = (double)(in[1]->data[i][l][m]);
00125                             }
00126                             newdiff += fabs((double)in[2]->data[i][j][k]-mean)*fabs((double)in[2]->data[i][j][k]-mean);
00127                         }
00128                         if (olddiff > newdiff  ||  olddiff < 0.0) {
00129                             olddiff = newdiff;
00130                             (*out)->data[0][j][k] = (*out)->data[0][j][k-1];
00131                         }
00132                     }
00133                     newdiff = 0.0;
00134                     for (i=0; i<in[1]->nbnd; i++) {
00135                         newdiff += fabs((double)(in[2]->data[i][j][k]-in[1]->data[i][j/2][k/2]))*fabs((double)(in[2]->data[i][j][k]-in[1]->data[i][j/2][k/2]));
00136                     }
00137                     if (olddiff > newdiff  ||  olddiff < 0.0) {
00138                         olddiff = newdiff;
00139                         (*out)->data[0][j][k] = in[0]->data[0][j/2][k/2];
00140                     }
00141                 }
00142                 break;
00143 
00144             /* type 2 */
00145             case 1:
00146                 if (j > 0  &&  (*out)->data[0][j-1][k] == in[0]->data[0][j/2][k/2]  &&  (*out)->data[0][j-1][k] == in[0]->data[0][j/2-1][k/2]   ||  j == 0  &&  k == (*out)->npix-1) {
00147                     (*out)->data[0][j][k] = in[0]->data[0][j/2][k/2];
00148                 } else {
00149                     olddiff = -1.0;
00150                     if (j > 0) {
00151                         l = (j-1)/2;
00152                         m = k/2;
00153                         newdiff = 0.0;
00154                         for (i=0; i<in[1]->nbnd; i++) {
00155                             if (l < in[0]->nlin-1  &&  m < in[0]->npix-1) {
00156                                 mean = (double)(in[1]->data[i][l+1][m+1]+in[1]->data[i][l+1][m]+in[1]->data[i][l][m+1]+in[1]->data[i][l][m]) / 4.0;
00157                             } else if (l < in[0]->nlin-1) {
00158                                 mean = (double)(in[1]->data[i][l+1][m]+in[1]->data[i][l][m]) / 2.0;
00159                             } else if (m < in[0]->npix-1) {
00160                                 mean = (double)(in[1]->data[i][l][m+1]+in[1]->data[i][l][m]) / 2.0;
00161                             } else {
00162                                 mean = (double)(in[1]->data[i][l][m]);
00163                             }
00164                             newdiff += fabs((double)in[2]->data[i][j][k]-mean)*fabs((double)in[2]->data[i][j][k]-mean);
00165                         }
00166                         olddiff = newdiff;
00167                         (*out)->data[0][j][k] = (*out)->data[0][j-1][k];
00168                     }
00169                     l = j/2;
00170                     m = (k-1)/2;
00171                     newdiff = 0.0;
00172                     for (i=0; i<in[1]->nbnd; i++) {
00173                         if (l > 0  &&  m > 0) {
00174                             mean = (double)(in[1]->data[i][l][m]+in[1]->data[i][l][m-1]+in[1]->data[i][l-1][m]+in[1]->data[i][l-1][m-1]) / 4.0;
00175                         } else if (l > 0) {
00176                             mean = (double)(in[1]->data[i][l][m]+in[1]->data[i][l-1][m]) / 2.0;
00177                         } else if (m > 0) {
00178                             mean = (double)(in[1]->data[i][l][m]+in[1]->data[i][l][m-1]) / 2.0;
00179                         } else {
00180                             mean = (double)(in[1]->data[i][l][m]);
00181                         }
00182                         newdiff += fabs((double)in[2]->data[i][j][k]-mean)*fabs((double)in[2]->data[i][j][k]-mean);
00183                     }
00184                     if (olddiff > newdiff  ||  olddiff < 0.0) {
00185                         olddiff = newdiff;
00186                         (*out)->data[0][j][k] = (*out)->data[0][j][k-1];
00187                     }
00188                     if (k/2 < in[1]->npix-1) {
00189                         newdiff = 0.0;
00190                         for (i=0; i<in[1]->nbnd; i++) {
00191                             newdiff += fabs((double)(in[2]->data[i][j][k]-in[1]->data[i][j/2][k/2+1]))*fabs((double)(in[2]->data[i][j][k]-in[1]->data[i][j/2][k/2+1]));
00192                         }
00193                         if (olddiff > newdiff  ||  olddiff < 0.0) {
00194                             olddiff = newdiff;
00195                             (*out)->data[0][j][k] = in[0]->data[0][j/2][k/2+1];
00196                         }
00197                     }
00198                     newdiff = 0.0;
00199                     for (i=0; i<in[1]->nbnd; i++) {
00200                         newdiff += fabs((double)(in[2]->data[i][j][k]-in[1]->data[i][j/2][k/2]))*fabs((double)(in[2]->data[i][j][k]-in[1]->data[i][j/2][k/2]));
00201                     }
00202                     if (olddiff > newdiff  ||  olddiff < 0.0) {
00203                         olddiff = newdiff;
00204                         (*out)->data[0][j][k] = in[0]->data[0][j/2][k/2];
00205                     }
00206                 }
00207                 break;
00208 
00209             /* type 3 */
00210             case 2:
00211                 if (k > 0  &&  (*out)->data[0][j][k-1] == in[0]->data[0][j/2][k/2]  &&  (*out)->data[0][j][k-1] != in[0]->data[0][j/2][k/2-1]  ||  j == (*out)->nlin-1  &&  k == 0) {
00212                     (*out)->data[0][j][k] = in[0]->data[0][j/2][k/2];
00213                 } else {
00214                     olddiff = -1.0;
00215                     l = (j-1)/2;
00216                     m = k/2;
00217                     newdiff = 0.0;
00218                     for (i=0; i<in[1]->nbnd; i++) {
00219                         if (l > 0  &&  m > 0) {
00220                             mean = (double)(in[1]->data[i][l][m]+in[1]->data[i][l][m-1]+in[1]->data[i][l-1][m]+in[1]->data[i][l-1][m-1]) / 4.0;
00221                         } else if (l > 0) {
00222                             mean = (double)(in[1]->data[i][l][m]+in[1]->data[i][l-1][m]) / 2.0;
00223                         } else if (m > 0) {
00224                             mean = (double)(in[1]->data[i][l][m]+in[1]->data[i][l][m-1]) / 2.0;
00225                         } else {
00226                             mean = (double)(in[1]->data[i][l][m]);
00227                         }
00228                         newdiff += fabs((double)in[2]->data[i][j][k]-mean)*fabs((double)in[2]->data[i][j][k]-mean);
00229                     }
00230                     olddiff = newdiff;
00231                     (*out)->data[0][j][k] = (*out)->data[0][j-1][k];
00232                     if (k > 0) {
00233                         l = j/2;
00234                         m = (k-1)/2;
00235                         newdiff = 0.0;
00236                         for (i=0; i<in[1]->nbnd; i++) {
00237                             if (l < in[0]->nlin-1  &&  m < in[0]->npix-1) {
00238                                 mean = (double)(in[1]->data[i][l+1][m+1]+in[1]->data[i][l+1][m]+in[1]->data[i][l][m+1]+in[1]->data[i][l][m]) / 4.0;
00239                             } else if (l < in[0]->nlin-1) {
00240                                 mean = (double)(in[1]->data[i][l+1][m]+in[1]->data[i][l][m]) / 2.0;
00241                             } else if (m < in[0]->npix-1) {
00242                                 mean = (double)(in[1]->data[i][l][m+1]+in[1]->data[i][l][m]) / 2.0;
00243                             } else {
00244                                 mean = (double)(in[1]->data[i][l][m]);
00245                             }
00246                             newdiff += fabs((double)in[2]->data[i][j][k]-mean)*fabs((double)in[2]->data[i][j][k]-mean);
00247                         }
00248                         if (olddiff > newdiff  ||  olddiff < 0.0) {
00249                             olddiff = newdiff;
00250                             (*out)->data[0][j][k] = (*out)->data[0][j][k-1];
00251                         }
00252                     }
00253                     newdiff = 0.0;
00254                     for (i=0; i<in[1]->nbnd; i++) {
00255                         newdiff += fabs((double)(in[2]->data[i][j][k]-in[1]->data[i][j/2][k/2]))*fabs((double)(in[2]->data[i][j][k]-in[1]->data[i][j/2][k/2]));
00256                     }
00257                     if (olddiff > newdiff  ||  olddiff < 0.0) {
00258                         olddiff = newdiff;
00259                         (*out)->data[0][j][k] = in[0]->data[0][j/2][k/2];
00260                     }
00261                 }
00262                 break;
00263 
00264             /* type 4 */
00265             case 3:
00266                 if (j == (*out)->nlin-1  &&  k == (*out)->npix-1) {
00267                     (*out)->data[0][j][k] = in[0]->data[0][j/2][k/2];
00268                 } else {
00269                     olddiff = -1.0;
00270                     l = (j-1)/2;
00271                     m = k/2;
00272                     newdiff = 0.0;
00273                     for (i=0; i<in[1]->nbnd; i++) {
00274                         if (l > 0  &&  m < in[1]->npix-1) {
00275                             mean = (double)(in[1]->data[i][l][m+1]+in[1]->data[i][l][m]+in[1]->data[i][l-1][m+1]+in[1]->data[i][l-1][m]) / 4.0;
00276                         } else if (l > 0) {
00277                             mean = (double)(in[1]->data[i][l][m]+in[1]->data[i][l-1][m]) / 2.0;
00278                         } else if (m < in[1]->npix-1) {
00279                             mean = (double)(in[1]->data[i][l][m+1]+in[1]->data[i][l][m]) / 2.0;
00280                         } else {
00281                             mean = (double)(in[1]->data[i][l][m]);
00282                         }
00283                         newdiff += fabs((double)in[2]->data[i][j][k]-mean)*fabs((double)in[2]->data[i][j][k]-mean);
00284                     }
00285                     olddiff = newdiff;
00286                     (*out)->data[0][j][k] = (*out)->data[0][j-1][k];
00287                     l = j/2;
00288                     m = (k-1)/2;
00289                     newdiff = 0.0;
00290                     for (i=0; i<in[1]->nbnd; i++) {
00291                         if (l < in[1]->nlin-1  &&  m > 0) {
00292                             mean = (double)(in[1]->data[i][l+1][m]+in[1]->data[i][l+1][m-1]+in[1]->data[i][l][m]+in[1]->data[i][l][m-1]) / 4.0;
00293                         } else if (l < in[1]->nlin-1) {
00294                             mean = (double)(in[1]->data[i][l+1][m]+in[1]->data[i][l][m]) / 2.0;
00295                         } else if (m > 0) {
00296                             mean = (double)(in[1]->data[i][l][m]+in[1]->data[i][l][m-1]) / 2.0;
00297                         } else {
00298                             mean = (double)(in[1]->data[i][l][m]);
00299                         }
00300                         newdiff += fabs((double)in[2]->data[i][j][k]-mean)*fabs((double)in[2]->data[i][j][k]-mean);
00301                     }
00302                     if (olddiff > newdiff  ||  olddiff < 0.0) {
00303                         olddiff = newdiff;
00304                         (*out)->data[0][j][k] = (*out)->data[0][j][k-1];
00305                     }
00306                     if (j/2 < in[1]->nlin-1) {
00307                         newdiff = 0.0;
00308                         for (i=0; i<in[1]->nbnd; i++) {
00309                             newdiff += fabs((double)(in[2]->data[i][j][k]-in[1]->data[i][j/2+1][k/2]))*fabs((double)(in[2]->data[i][j][k]-in[1]->data[i][j/2+1][k/2]));
00310                         }
00311                         if (olddiff > newdiff  ||  olddiff < 0.0) {
00312                             olddiff = newdiff;
00313                             (*out)->data[0][j][k] = in[0]->data[0][j/2+1][k/2];
00314                         }
00315                     }
00316                     if (k/2 < in[1]->npix-1) {
00317                         newdiff = 0.0;
00318                         for (i=0; i<in[1]->nbnd; i++) {
00319                             newdiff += fabs((double)(in[2]->data[i][j][k]-in[1]->data[i][j/2][k/2+1]))*fabs((double)(in[2]->data[i][j][k]-in[1]->data[i][j/2][k/2+1]));
00320                         }
00321                         if (olddiff > newdiff  ||  olddiff < 0.0) {
00322                             olddiff = newdiff;
00323                             (*out)->data[0][j][k] = in[0]->data[0][j/2][k/2+1];
00324                         }
00325                     }
00326                 }
00327                 break;
00328 
00329             }
00330         }
00331         if (LINES  &&  !PROGRESS(psum+=pinc)) goto the_end;
00332     }
00333 
00334     the_end:
00335     if (TIMES)   TIMING(T_EXIT);
00336 }