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
00023 REFINE (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   )
00030 {
00031   register short i, j, k, l, m;
00032   char msg[SLEN];
00033   double olddiff, newdiff, mean, pinc, psum;
00034 
00035   if (TIMES)
00036     TIMING (T_REFINE);
00037   if (NAMES)
00038     {
00039       MESSAGE ('I', "");
00040       MESSAGE ('I', "REFINE");
00041       MESSAGE ('I', "");
00042       sprintf (msg, " Reduced label map:   %s", in[0]->text);
00043       MESSAGE ('I', msg);
00044       sprintf (msg, " Reduced image:       %s", in[1]->text);
00045       MESSAGE ('I', msg);
00046       sprintf (msg, " Refined image:       %s", in[2]->text);
00047       MESSAGE ('I', msg);
00048       MESSAGE ('I', " ...............");
00049     }
00050 
00051   /* check input */
00052   for (i = 0; i < 3; i++)
00053     {
00054       if (!CHECKIMG (in[i]))
00055         {
00056           MESSAGE ('E', " Can't identify image.");
00057           goto the_end;
00058         }
00059     }
00060   if (in[0]->npix != in[1]->npix)
00061     {
00062       MESSAGE ('E',
00063                " Number of pixels/line must be identical in reduced image and reduced label map.");
00064       goto the_end;
00065     }
00066   else if (in[0]->nlin != in[1]->nlin)
00067     {
00068       MESSAGE ('E',
00069                " Number of lines must be identical in reduced image and reduced label map.");
00070       goto the_end;
00071     }
00072   else if (in[0]->nbnd != 1)
00073     {
00074       MESSAGE ('E', " Reduced label map must be single-band.");
00075       goto the_end;
00076     }
00077   else if (2 * in[1]->npix != in[2]->npix)
00078     {
00079       MESSAGE ('E',
00080                " Number of pixels/line in reduced image must be half the number of pixels/line in regular image.");
00081       goto the_end;
00082     }
00083   else if (2 * in[1]->nlin != in[2]->nlin)
00084     {
00085       MESSAGE ('E',
00086                " Number of lines in reduced image must be half the number of lines in regular image.");
00087       goto the_end;
00088     }
00089   else if (in[1]->nbnd != in[2]->nbnd)
00090     {
00091       MESSAGE ('E',
00092                " Number of bands must be identical in reduced image and regular image");
00093       goto the_end;
00094     }
00095 
00096   /* create image of appropriate size */
00097   if (!CHECKIMG (*out))
00098     GETMEM (1, in[2]->nlin, in[2]->npix, out);
00099   if (!*out)
00100     goto the_end;
00101 
00102   /* initialize progress indicator */
00103   if (LINES && !PROGRESS (psum = 0.0))
00104     goto the_end;
00105   pinc = 1.0 / (double) (*out)->nlin;
00106 
00107   for (j = 0; j < (*out)->nlin; j++)
00108     {
00109       for (k = 0; k < (*out)->npix; k++)
00110         {
00111           switch (j % 2 << 1 | k % 2)
00112             {
00113 
00114               /* type 1 */
00115             case 0:
00116               if (j > 0
00117                   && (*out)->data[0][j - 1][k] == in[0]->data[0][j / 2][k / 2]
00118                   && (*out)->data[0][j - 1][k] !=
00119                   in[0]->data[0][j / 2 - 1][k / 2] || k > 0
00120                   && (*out)->data[0][j][k - 1] == in[0]->data[0][j / 2][k / 2]
00121                   && (*out)->data[0][j][k - 1] !=
00122                   in[0]->data[0][j / 2][k / 2 - 1] || j == 0 && k == 0)
00123                 {
00124                   (*out)->data[0][j][k] = in[0]->data[0][j / 2][k / 2];
00125                 }
00126               else
00127                 {
00128                   olddiff = -1.0;
00129                   if (j > 0)
00130                     {
00131                       l = (j - 1) / 2;
00132                       m = k / 2;
00133                       newdiff = 0.0;
00134                       for (i = 0; i < in[1]->nbnd; i++)
00135                         {
00136                           if (l < in[1]->nlin - 1 && m > 0)
00137                             {
00138                               mean =
00139                                 (double) (in[1]->data[i][l + 1][m] +
00140                                           in[1]->data[i][l + 1][m - 1] +
00141                                           in[1]->data[i][l][m] +
00142                                           in[1]->data[i][l][m - 1]) / 4.0;
00143                             }
00144                           else if (l < in[1]->nlin - 1)
00145                             {
00146                               mean =
00147                                 (double) (in[1]->data[i][l + 1][m] +
00148                                           in[1]->data[i][l][m]) / 2.0;
00149                             }
00150                           else if (m > 0)
00151                             {
00152                               mean =
00153                                 (double) (in[1]->data[i][l][m] +
00154                                           in[1]->data[i][l][m - 1]) / 2.0;
00155                             }
00156                           else
00157                             {
00158                               mean = (double) (in[1]->data[i][l][m]);
00159                             }
00160                           newdiff +=
00161                             fabs ((double) in[2]->data[i][j][k] -
00162                                   mean) *
00163                             fabs ((double) in[2]->data[i][j][k] - mean);
00164                         }
00165                       olddiff = newdiff;
00166                       (*out)->data[0][j][k] = (*out)->data[0][j - 1][k];
00167                     }
00168                   if (k > 0)
00169                     {
00170                       l = j / 2;
00171                       m = (k - 1) / 2;
00172                       newdiff = 0.0;
00173                       for (i = 0; i < in[1]->nbnd; i++)
00174                         {
00175                           if (l > 0 && m < in[1]->npix - 1)
00176                             {
00177                               mean =
00178                                 (double) (in[1]->data[i][l][m + 1] +
00179                                           in[1]->data[i][l][m] +
00180                                           in[1]->data[i][l - 1][m + 1] +
00181                                           in[1]->data[i][l - 1][m]) / 4.0;
00182                             }
00183                           else if (l > 0)
00184                             {
00185                               mean =
00186                                 (double) (in[1]->data[i][l][m] +
00187                                           in[1]->data[i][l - 1][m]) / 2.0;
00188                             }
00189                           else if (m < in[1]->npix - 1)
00190                             {
00191                               mean =
00192                                 (double) (in[1]->data[i][l][m + 1] +
00193                                           in[1]->data[i][l][m]) / 2.0;
00194                             }
00195                           else
00196                             {
00197                               mean = (double) (in[1]->data[i][l][m]);
00198                             }
00199                           newdiff +=
00200                             fabs ((double) in[2]->data[i][j][k] -
00201                                   mean) *
00202                             fabs ((double) in[2]->data[i][j][k] - mean);
00203                         }
00204                       if (olddiff > newdiff || olddiff < 0.0)
00205                         {
00206                           olddiff = newdiff;
00207                           (*out)->data[0][j][k] = (*out)->data[0][j][k - 1];
00208                         }
00209                     }
00210                   newdiff = 0.0;
00211                   for (i = 0; i < in[1]->nbnd; i++)
00212                     {
00213                       newdiff +=
00214                         fabs ((double)
00215                               (in[2]->data[i][j][k] -
00216                                in[1]->data[i][j / 2][k / 2])) *
00217                         fabs ((double)
00218                               (in[2]->data[i][j][k] -
00219                                in[1]->data[i][j / 2][k / 2]));
00220                     }
00221                   if (olddiff > newdiff || olddiff < 0.0)
00222                     {
00223                       olddiff = newdiff;
00224                       (*out)->data[0][j][k] = in[0]->data[0][j / 2][k / 2];
00225                     }
00226                 }
00227               break;
00228 
00229               /* type 2 */
00230             case 1:
00231               if (j > 0
00232                   && (*out)->data[0][j - 1][k] == in[0]->data[0][j / 2][k / 2]
00233                   && (*out)->data[0][j - 1][k] ==
00234                   in[0]->data[0][j / 2 - 1][k / 2] || j == 0
00235                   && k == (*out)->npix - 1)
00236                 {
00237                   (*out)->data[0][j][k] = in[0]->data[0][j / 2][k / 2];
00238                 }
00239               else
00240                 {
00241                   olddiff = -1.0;
00242                   if (j > 0)
00243                     {
00244                       l = (j - 1) / 2;
00245                       m = k / 2;
00246                       newdiff = 0.0;
00247                       for (i = 0; i < in[1]->nbnd; i++)
00248                         {
00249                           if (l < in[0]->nlin - 1 && m < in[0]->npix - 1)
00250                             {
00251                               mean =
00252                                 (double) (in[1]->data[i][l + 1][m + 1] +
00253                                           in[1]->data[i][l + 1][m] +
00254                                           in[1]->data[i][l][m + 1] +
00255                                           in[1]->data[i][l][m]) / 4.0;
00256                             }
00257                           else if (l < in[0]->nlin - 1)
00258                             {
00259                               mean =
00260                                 (double) (in[1]->data[i][l + 1][m] +
00261                                           in[1]->data[i][l][m]) / 2.0;
00262                             }
00263                           else if (m < in[0]->npix - 1)
00264                             {
00265                               mean =
00266                                 (double) (in[1]->data[i][l][m + 1] +
00267                                           in[1]->data[i][l][m]) / 2.0;
00268                             }
00269                           else
00270                             {
00271                               mean = (double) (in[1]->data[i][l][m]);
00272                             }
00273                           newdiff +=
00274                             fabs ((double) in[2]->data[i][j][k] -
00275                                   mean) *
00276                             fabs ((double) in[2]->data[i][j][k] - mean);
00277                         }
00278                       olddiff = newdiff;
00279                       (*out)->data[0][j][k] = (*out)->data[0][j - 1][k];
00280                     }
00281                   l = j / 2;
00282                   m = (k - 1) / 2;
00283                   newdiff = 0.0;
00284                   for (i = 0; i < in[1]->nbnd; i++)
00285                     {
00286                       if (l > 0 && m > 0)
00287                         {
00288                           mean =
00289                             (double) (in[1]->data[i][l][m] +
00290                                       in[1]->data[i][l][m - 1] +
00291                                       in[1]->data[i][l - 1][m] +
00292                                       in[1]->data[i][l - 1][m - 1]) / 4.0;
00293                         }
00294                       else if (l > 0)
00295                         {
00296                           mean =
00297                             (double) (in[1]->data[i][l][m] +
00298                                       in[1]->data[i][l - 1][m]) / 2.0;
00299                         }
00300                       else if (m > 0)
00301                         {
00302                           mean =
00303                             (double) (in[1]->data[i][l][m] +
00304                                       in[1]->data[i][l][m - 1]) / 2.0;
00305                         }
00306                       else
00307                         {
00308                           mean = (double) (in[1]->data[i][l][m]);
00309                         }
00310                       newdiff +=
00311                         fabs ((double) in[2]->data[i][j][k] -
00312                               mean) * fabs ((double) in[2]->data[i][j][k] -
00313                                             mean);
00314                     }
00315                   if (olddiff > newdiff || olddiff < 0.0)
00316                     {
00317                       olddiff = newdiff;
00318                       (*out)->data[0][j][k] = (*out)->data[0][j][k - 1];
00319                     }
00320                   if (k / 2 < in[1]->npix - 1)
00321                     {
00322                       newdiff = 0.0;
00323                       for (i = 0; i < in[1]->nbnd; i++)
00324                         {
00325                           newdiff +=
00326                             fabs ((double)
00327                                   (in[2]->data[i][j][k] -
00328                                    in[1]->data[i][j / 2][k / 2 +
00329                                                          1])) *
00330                             fabs ((double)
00331                                   (in[2]->data[i][j][k] -
00332                                    in[1]->data[i][j / 2][k / 2 + 1]));
00333                         }
00334                       if (olddiff > newdiff || olddiff < 0.0)
00335                         {
00336                           olddiff = newdiff;
00337                           (*out)->data[0][j][k] =
00338                             in[0]->data[0][j / 2][k / 2 + 1];
00339                         }
00340                     }
00341                   newdiff = 0.0;
00342                   for (i = 0; i < in[1]->nbnd; i++)
00343                     {
00344                       newdiff +=
00345                         fabs ((double)
00346                               (in[2]->data[i][j][k] -
00347                                in[1]->data[i][j / 2][k / 2])) *
00348                         fabs ((double)
00349                               (in[2]->data[i][j][k] -
00350                                in[1]->data[i][j / 2][k / 2]));
00351                     }
00352                   if (olddiff > newdiff || olddiff < 0.0)
00353                     {
00354                       olddiff = newdiff;
00355                       (*out)->data[0][j][k] = in[0]->data[0][j / 2][k / 2];
00356                     }
00357                 }
00358               break;
00359 
00360               /* type 3 */
00361             case 2:
00362               if (k > 0
00363                   && (*out)->data[0][j][k - 1] == in[0]->data[0][j / 2][k / 2]
00364                   && (*out)->data[0][j][k - 1] !=
00365                   in[0]->data[0][j / 2][k / 2 - 1] || j == (*out)->nlin - 1
00366                   && k == 0)
00367                 {
00368                   (*out)->data[0][j][k] = in[0]->data[0][j / 2][k / 2];
00369                 }
00370               else
00371                 {
00372                   olddiff = -1.0;
00373                   l = (j - 1) / 2;
00374                   m = k / 2;
00375                   newdiff = 0.0;
00376                   for (i = 0; i < in[1]->nbnd; i++)
00377                     {
00378                       if (l > 0 && m > 0)
00379                         {
00380                           mean =
00381                             (double) (in[1]->data[i][l][m] +
00382                                       in[1]->data[i][l][m - 1] +
00383                                       in[1]->data[i][l - 1][m] +
00384                                       in[1]->data[i][l - 1][m - 1]) / 4.0;
00385                         }
00386                       else if (l > 0)
00387                         {
00388                           mean =
00389                             (double) (in[1]->data[i][l][m] +
00390                                       in[1]->data[i][l - 1][m]) / 2.0;
00391                         }
00392                       else if (m > 0)
00393                         {
00394                           mean =
00395                             (double) (in[1]->data[i][l][m] +
00396                                       in[1]->data[i][l][m - 1]) / 2.0;
00397                         }
00398                       else
00399                         {
00400                           mean = (double) (in[1]->data[i][l][m]);
00401                         }
00402                       newdiff +=
00403                         fabs ((double) in[2]->data[i][j][k] -
00404                               mean) * fabs ((double) in[2]->data[i][j][k] -
00405                                             mean);
00406                     }
00407                   olddiff = newdiff;
00408                   (*out)->data[0][j][k] = (*out)->data[0][j - 1][k];
00409                   if (k > 0)
00410                     {
00411                       l = j / 2;
00412                       m = (k - 1) / 2;
00413                       newdiff = 0.0;
00414                       for (i = 0; i < in[1]->nbnd; i++)
00415                         {
00416                           if (l < in[0]->nlin - 1 && m < in[0]->npix - 1)
00417                             {
00418                               mean =
00419                                 (double) (in[1]->data[i][l + 1][m + 1] +
00420                                           in[1]->data[i][l + 1][m] +
00421                                           in[1]->data[i][l][m + 1] +
00422                                           in[1]->data[i][l][m]) / 4.0;
00423                             }
00424                           else if (l < in[0]->nlin - 1)
00425                             {
00426                               mean =
00427                                 (double) (in[1]->data[i][l + 1][m] +
00428                                           in[1]->data[i][l][m]) / 2.0;
00429                             }
00430                           else if (m < in[0]->npix - 1)
00431                             {
00432                               mean =
00433                                 (double) (in[1]->data[i][l][m + 1] +
00434                                           in[1]->data[i][l][m]) / 2.0;
00435                             }
00436                           else
00437                             {
00438                               mean = (double) (in[1]->data[i][l][m]);
00439                             }
00440                           newdiff +=
00441                             fabs ((double) in[2]->data[i][j][k] -
00442                                   mean) *
00443                             fabs ((double) in[2]->data[i][j][k] - mean);
00444                         }
00445                       if (olddiff > newdiff || olddiff < 0.0)
00446                         {
00447                           olddiff = newdiff;
00448                           (*out)->data[0][j][k] = (*out)->data[0][j][k - 1];
00449                         }
00450                     }
00451                   newdiff = 0.0;
00452                   for (i = 0; i < in[1]->nbnd; i++)
00453                     {
00454                       newdiff +=
00455                         fabs ((double)
00456                               (in[2]->data[i][j][k] -
00457                                in[1]->data[i][j / 2][k / 2])) *
00458                         fabs ((double)
00459                               (in[2]->data[i][j][k] -
00460                                in[1]->data[i][j / 2][k / 2]));
00461                     }
00462                   if (olddiff > newdiff || olddiff < 0.0)
00463                     {
00464                       olddiff = newdiff;
00465                       (*out)->data[0][j][k] = in[0]->data[0][j / 2][k / 2];
00466                     }
00467                 }
00468               break;
00469 
00470               /* type 4 */
00471             case 3:
00472               if (j == (*out)->nlin - 1 && k == (*out)->npix - 1)
00473                 {
00474                   (*out)->data[0][j][k] = in[0]->data[0][j / 2][k / 2];
00475                 }
00476               else
00477                 {
00478                   olddiff = -1.0;
00479                   l = (j - 1) / 2;
00480                   m = k / 2;
00481                   newdiff = 0.0;
00482                   for (i = 0; i < in[1]->nbnd; i++)
00483                     {
00484                       if (l > 0 && m < in[1]->npix - 1)
00485                         {
00486                           mean =
00487                             (double) (in[1]->data[i][l][m + 1] +
00488                                       in[1]->data[i][l][m] +
00489                                       in[1]->data[i][l - 1][m + 1] +
00490                                       in[1]->data[i][l - 1][m]) / 4.0;
00491                         }
00492                       else if (l > 0)
00493                         {
00494                           mean =
00495                             (double) (in[1]->data[i][l][m] +
00496                                       in[1]->data[i][l - 1][m]) / 2.0;
00497                         }
00498                       else if (m < in[1]->npix - 1)
00499                         {
00500                           mean =
00501                             (double) (in[1]->data[i][l][m + 1] +
00502                                       in[1]->data[i][l][m]) / 2.0;
00503                         }
00504                       else
00505                         {
00506                           mean = (double) (in[1]->data[i][l][m]);
00507                         }
00508                       newdiff +=
00509                         fabs ((double) in[2]->data[i][j][k] -
00510                               mean) * fabs ((double) in[2]->data[i][j][k] -
00511                                             mean);
00512                     }
00513                   olddiff = newdiff;
00514                   (*out)->data[0][j][k] = (*out)->data[0][j - 1][k];
00515                   l = j / 2;
00516                   m = (k - 1) / 2;
00517                   newdiff = 0.0;
00518                   for (i = 0; i < in[1]->nbnd; i++)
00519                     {
00520                       if (l < in[1]->nlin - 1 && m > 0)
00521                         {
00522                           mean =
00523                             (double) (in[1]->data[i][l + 1][m] +
00524                                       in[1]->data[i][l + 1][m - 1] +
00525                                       in[1]->data[i][l][m] +
00526                                       in[1]->data[i][l][m - 1]) / 4.0;
00527                         }
00528                       else if (l < in[1]->nlin - 1)
00529                         {
00530                           mean =
00531                             (double) (in[1]->data[i][l + 1][m] +
00532                                       in[1]->data[i][l][m]) / 2.0;
00533                         }
00534                       else if (m > 0)
00535                         {
00536                           mean =
00537                             (double) (in[1]->data[i][l][m] +
00538                                       in[1]->data[i][l][m - 1]) / 2.0;
00539                         }
00540                       else
00541                         {
00542                           mean = (double) (in[1]->data[i][l][m]);
00543                         }
00544                       newdiff +=
00545                         fabs ((double) in[2]->data[i][j][k] -
00546                               mean) * fabs ((double) in[2]->data[i][j][k] -
00547                                             mean);
00548                     }
00549                   if (olddiff > newdiff || olddiff < 0.0)
00550                     {
00551                       olddiff = newdiff;
00552                       (*out)->data[0][j][k] = (*out)->data[0][j][k - 1];
00553                     }
00554                   if (j / 2 < in[1]->nlin - 1)
00555                     {
00556                       newdiff = 0.0;
00557                       for (i = 0; i < in[1]->nbnd; i++)
00558                         {
00559                           newdiff +=
00560                             fabs ((double)
00561                                   (in[2]->data[i][j][k] -
00562                                    in[1]->data[i][j / 2 +
00563                                                   1][k / 2])) *
00564                             fabs ((double)
00565                                   (in[2]->data[i][j][k] -
00566                                    in[1]->data[i][j / 2 + 1][k / 2]));
00567                         }
00568                       if (olddiff > newdiff || olddiff < 0.0)
00569                         {
00570                           olddiff = newdiff;
00571                           (*out)->data[0][j][k] =
00572                             in[0]->data[0][j / 2 + 1][k / 2];
00573                         }
00574                     }
00575                   if (k / 2 < in[1]->npix - 1)
00576                     {
00577                       newdiff = 0.0;
00578                       for (i = 0; i < in[1]->nbnd; i++)
00579                         {
00580                           newdiff +=
00581                             fabs ((double)
00582                                   (in[2]->data[i][j][k] -
00583                                    in[1]->data[i][j / 2][k / 2 +
00584                                                          1])) *
00585                             fabs ((double)
00586                                   (in[2]->data[i][j][k] -
00587                                    in[1]->data[i][j / 2][k / 2 + 1]));
00588                         }
00589                       if (olddiff > newdiff || olddiff < 0.0)
00590                         {
00591                           olddiff = newdiff;
00592                           (*out)->data[0][j][k] =
00593                             in[0]->data[0][j / 2][k / 2 + 1];
00594                         }
00595                     }
00596                 }
00597               break;
00598 
00599             }
00600         }
00601       if (LINES && !PROGRESS (psum += pinc))
00602         goto the_end;
00603     }
00604 
00605 the_end:
00606   if (TIMES)
00607     TIMING (T_EXIT);
00608 }