00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #include "mrilib.h"
00008 
00009 /*----------------------------------------------------------------
00010    find clusters of points !=0; return an array of clusters
00011    [the input array fim is destroyed in this
00012     computation; it may be recreated by MCW_cluster_to_vol]
00013 
00014    Nov 1995: fim array may now be short, byte, or float,
00015              signified by new ftype argument.
00016 
00017    Coordinates of voxels in clusters are now stored as 3 separate
00018    short integers, to correct error due to abiguity in
00019    identification of clusters.
00020    BDW  06 March 1997
00021 ------------------------------------------------------------------*/
00022 
00023 MCW_cluster_array * MCW_find_clusters(
00024                        int nx, int ny, int nz,
00025                        float dx, float dy, float dz,
00026                        int ftype , void * fim ,
00027                        float max_dist )
00028 {
00029    MCW_cluster_array * clust_arr ;
00030    MCW_cluster       * clust , * mask ;
00031    int ii,jj,kk ,  nxy,nxyz , ijk , ijk_last , mnum ;
00032    int icl , jma , ijkcl , ijkma , did_one ;
00033    float fimv ;
00034    short * sfar ;
00035    float * ffar ;
00036    byte  * bfar ;
00037    short ic, jc, kc;
00038    short im, jm, km;
00039 
00040 ENTRY("MCW_find_clusters") ;
00041 
00042    if( fim == NULL ) RETURN(NULL) ;
00043 
00044    switch( ftype ){
00045       default: RETURN(NULL) ;
00046       case MRI_short:  sfar = (short *) fim ; break ;
00047       case MRI_byte :  bfar = (byte  *) fim ; break ;
00048       case MRI_float:  ffar = (float *) fim ; break ;
00049    }
00050 
00051    /*--- make a cluster that is a mask of points closer than max_dist ---*/
00052 
00053    mask = MCW_build_mask (nx, ny, nz, dx, dy, dz, max_dist);
00054    if (mask == NULL)
00055    {
00056       fprintf (stderr, "Unable to build mask in MCW_find_clusters");
00057       RETURN(NULL) ;
00058    }
00059 
00060    nxy = nx*ny ; nxyz = nxy * nz ;
00061 
00062    mnum = mask->num_pt ;
00063 
00064 
00065    /*--- scan through array, find nonzero point, build a cluster, ... ---*/
00066 
00067    INIT_CLARR(clust_arr) ;
00068 
00069    ijk_last = 0 ;
00070    do {
00071       switch( ftype ){
00072          case MRI_short:
00073             for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( sfar[ijk] != 0 ) break ;
00074             if( ijk < nxyz ){
00075                fimv = sfar[ijk] ; sfar[ijk] = 0 ;  /* save found point */
00076             }
00077          break ;
00078 
00079          case MRI_byte:
00080             for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( bfar[ijk] != 0 ) break ;
00081             if( ijk < nxyz ){
00082                fimv = bfar[ijk] ; bfar[ijk] = 0 ;  /* save found point */
00083             }
00084          break ;
00085 
00086          case MRI_float:
00087             for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( ffar[ijk] != 0.0 ) break ;
00088             if( ijk < nxyz ){
00089                fimv = ffar[ijk] ; ffar[ijk] = 0.0 ;  /* save found point */
00090             }
00091          break ;
00092       }
00093       if( ijk == nxyz ) break ;  /* didn't find any! */
00094 
00095 #ifdef CLUST_DEBUG
00096 printf("  starting cluster at ijk=%d\n",ijk) ;
00097 #endif
00098 
00099       ijk_last = ijk+1 ;         /* start here next time */
00100 
00101       INIT_CLUSTER(clust) ;                  /* make a new cluster */
00102       IJK_TO_THREE(ijk,ic,jc,kc,nx,nxy) ;
00103       ADDTO_CLUSTER( clust , ic, jc, kc, fimv ) ;  /* start it off */
00104 
00105       /*--
00106         for each point in cluster:
00107            check points offset by the mask for nonzero entries in fim
00108            enter those into cluster
00109            continue until end of cluster is reached
00110              (note that cluster is expanding as we progress)
00111       --*/
00112 
00113       switch( ftype ){
00114          case MRI_short:
00115             for( icl=0 ; icl < clust->num_pt ; icl++ ){
00116                ic = clust->i[icl];
00117                jc = clust->j[icl];
00118                kc = clust->k[icl];
00119 
00120                for( jma=0 ; jma < mnum ; jma++ ){
00121                   im = ic + mask->i[jma];
00122                   jm = jc + mask->j[jma];
00123                   km = kc + mask->k[jma];
00124                   if( im < 0 || im >= nx ||
00125                       jm < 0 || jm >= ny || km < 0 || km >= nz ) continue ;
00126 
00127                   ijkma = THREE_TO_IJK (im, jm, km, nx, nxy);
00128                   if( ijkma < ijk_last || ijkma >= nxyz || sfar[ijkma] == 0 ) continue ;
00129 
00130                   ADDTO_CLUSTER( clust , im, jm, km, sfar[ijkma] ) ;
00131                   sfar[ijkma] = 0 ;
00132                }
00133             }
00134          break ;
00135 
00136          case MRI_byte:
00137             for( icl=0 ; icl < clust->num_pt ; icl++ ){
00138                ic = clust->i[icl];
00139                jc = clust->j[icl];
00140                kc = clust->k[icl];
00141 
00142                for( jma=0 ; jma < mnum ; jma++ ){
00143                   im = ic + mask->i[jma];
00144                   jm = jc + mask->j[jma];
00145                   km = kc + mask->k[jma];
00146                   if( im < 0 || im >= nx ||
00147                       jm < 0 || jm >= ny || km < 0 || km >= nz ) continue ;
00148 
00149                   ijkma = THREE_TO_IJK (im, jm, km, nx, nxy);
00150                   if( ijkma < ijk_last || ijkma >= nxyz || bfar[ijkma] == 0 ) continue ;
00151 
00152                   ADDTO_CLUSTER( clust , im, jm, km, bfar[ijkma] ) ;
00153                   bfar[ijkma] = 0 ;
00154                }
00155             }
00156          break ;
00157 
00158          case MRI_float:
00159             for( icl=0 ; icl < clust->num_pt ; icl++ ){
00160                ic = clust->i[icl];
00161                jc = clust->j[icl];
00162                kc = clust->k[icl];
00163 
00164                for( jma=0 ; jma < mnum ; jma++ ){
00165                   im = ic + mask->i[jma];
00166                   jm = jc + mask->j[jma];
00167                   km = kc + mask->k[jma];
00168                   if( im < 0 || im >= nx ||
00169                       jm < 0 || jm >= ny || km < 0 || km >= nz ) continue ;
00170 
00171                   ijkma = THREE_TO_IJK (im, jm, km, nx, nxy);
00172                   if( ijkma < ijk_last || ijkma >= nxyz || ffar[ijkma] == 0.0 ) continue ;
00173 
00174                   ADDTO_CLUSTER( clust , im, jm, km, ffar[ijkma] ) ;
00175                   ffar[ijkma] = 0.0 ;
00176                }
00177             }
00178          break ;
00179       }
00180 
00181       ADDTO_CLARR(clust_arr,clust) ;
00182    } while( 1 ) ;
00183 
00184    KILL_CLUSTER(mask) ;
00185 
00186    if( clust_arr->num_clu <= 0 ){ DESTROY_CLARR(clust_arr) ; }
00187 
00188    RETURN(clust_arr) ;
00189 }
00190 
00191 /*---------------------------------------------------------------
00192   Write the points stored in a cluster back into a volume
00193 
00194    Coordinates of voxels in clusters are now stored as 3 separate
00195    short integers, to correct error due to abiguity in
00196    identification of clusters.
00197    BDW  06 March 1997
00198 -----------------------------------------------------------------*/
00199 
00200 void MCW_cluster_to_vol( int nx , int ny , int nz ,
00201                          int ftype , void * fim , MCW_cluster * clust )
00202 {
00203    int icl, ijk ;
00204    int nxy ;
00205    short * sfar ;
00206    float * ffar ;
00207    byte  * bfar ;
00208 
00209 ENTRY("MCW_cluster_to_vol") ;
00210 
00211    if( fim == NULL || clust == NULL ) EXRETURN ;
00212 
00213    nxy = nx * ny;
00214 
00215    switch( ftype ){
00216       case MRI_short:
00217          sfar = (short *) fim ;
00218          for( icl=0 ; icl < clust->num_pt ; icl++ )
00219            {
00220              ijk = THREE_TO_IJK (clust->i[icl], clust->j[icl], clust->k[icl],
00221                                  nx, nxy);
00222              sfar[ijk] = clust->mag[icl] ;
00223            }
00224       EXRETURN ;
00225 
00226       case MRI_byte:
00227          bfar = (byte *) fim ;
00228          for( icl=0 ; icl < clust->num_pt ; icl++ )
00229            {
00230              ijk = THREE_TO_IJK (clust->i[icl], clust->j[icl], clust->k[icl],
00231                                  nx, nxy);
00232              bfar[ijk] = clust->mag[icl] ;
00233            }
00234       EXRETURN ;
00235 
00236       case MRI_float:
00237          ffar = (float *) fim ;
00238          for( icl=0 ; icl < clust->num_pt ; icl++ )
00239            {
00240              ijk = THREE_TO_IJK (clust->i[icl], clust->j[icl], clust->k[icl],
00241                                  nx, nxy);
00242              ffar[ijk] = clust->mag[icl] ;
00243            }
00244       EXRETURN ;
00245    }
00246 
00247    EXRETURN ;  /* should not be reached */
00248 }
00249 
00250 
00251 /*----------------------------------------------------------------------------
00252   Erosion and dilation of 3d clusters.
00253 
00254   Purpose:  Sometimes a cluster consists of several main bodies connected by
00255             a narrow path.  The objective is to sever the connecting path,
00256             while keeping the main bodies intact, thereby producing separate
00257             clusters.
00258 
00259   Method:   Erosion is applied to eliminate the outer layer of voxels.  This
00260             should cut the connecting path.  The original main bodies are
00261             then restored by dilating the result.
00262 
00263   Author:   B. Douglas Ward
00264 
00265   Date:     16 June 1998
00266 
00267 
00268 -----------------------------------------------------------------------------*/
00269 
00270 void MCW_erode_clusters
00271 (
00272   int nx, int ny, int nz,           /* dimensions of volume fim */
00273   float dx, float dy, float dz,     /* voxel dimensions */
00274   int ftype,                        /* data type */
00275   void * fim,                       /* volume data */
00276   float max_dist,                   /* voxel connectivity radius */
00277   float pv,                         /* fraction pv of voxels within max_dist must
00278                                        be in the cluster */
00279   int dilate                        /* boolean for perform dilation phase */
00280 )
00281 
00282 {
00283   MCW_cluster * mask = NULL;        /* mask determines nbhd membership */
00284   int minimum;                      /* minimum number of voxels in nbhd */
00285   int count;                        /* count of voxels in neighborhood */
00286   int nxy, nxyz;                    /* numbers of voxels */
00287   int ijk, iv, jv, kv;              /* voxel indices */
00288   int ijkm, im, jm, km;             /* voxel indices */
00289   int imask, nmask;                 /* mask indices */
00290   short * sfar;                     /* pointer to short data */
00291   byte  * bfar;                     /* pointer to byte data */
00292   float * ffar;                     /* pointer to float data */
00293   float * efim = NULL;              /* copy of eroded voxels */
00294 
00295 ENTRY("MCW_erode_clusters") ;
00296 
00297 
00298   /*----- Just in case -----*/
00299   if ( fim == NULL )  EXRETURN;
00300 
00301 
00302   /*----- Initialize local variables -----*/
00303   nxy = nx * ny;   nxyz = nxy * nz;
00304 
00305 
00306   /*----- Set pointer to input data -----*/
00307   switch (ftype)
00308     {
00309     default:  EXRETURN;
00310     case MRI_short:  sfar = (short *) fim;  break;
00311     case MRI_byte :  bfar = (byte  *) fim;  break;
00312     case MRI_float:  ffar = (float *) fim;  break;
00313     }
00314 
00315 
00316   /*----- Initialization for copy of eroded voxels -----*/
00317   efim = (float *) malloc (sizeof(float) * nxyz);
00318   if (efim == NULL)
00319     {
00320       fprintf (stderr, "Unable to allocate memory in MCW_erode_clusters");
00321       EXRETURN;
00322     }
00323   for (ijk = 0;  ijk < nxyz;  ijk++)
00324     efim[ijk] = 0.0;
00325 
00326 
00327   /*--- Make a cluster that is a mask of points closer than max_dist ---*/
00328   mask = MCW_build_mask (nx, ny, nz, dx, dy, dz, max_dist);
00329   if (mask == NULL)
00330     {
00331       fprintf (stderr, "Unable to build mask in MCW_erode_clusters");
00332       EXRETURN;
00333     }
00334 
00335 
00336   /*----- Calculate minimum number of voxels in nbhd. for non-erosion -----*/
00337   nmask = mask->num_pt ;
00338   minimum = floor(pv*nmask + 0.99);
00339   if (minimum <= 0)  EXRETURN;     /*----- Nothing will be eroded -----*/
00340 
00341 
00342   /*----- Step 1:  Identify voxels to be eroded -----*/
00343   switch (ftype)
00344     {
00345     case MRI_short:
00346       for (ijk = 0;  ijk < nxyz;  ijk++)
00347         {       
00348           if (sfar[ijk] == 0)  continue;
00349           IJK_TO_THREE (ijk, iv, jv, kv, nx, nxy);
00350 
00351           /*----- Count number of active voxels in the neighborhood -----*/
00352           count = 0;
00353           for (imask = 0;  imask < nmask;  imask++)
00354             {
00355               im = iv + mask->i[imask];
00356               jm = jv + mask->j[imask];
00357               km = kv + mask->k[imask];
00358               if ( im < 0 || jm < 0 || km < 0 ||
00359                    im >= nx || jm >= ny || km >= nz )  continue;
00360               ijkm = THREE_TO_IJK (im, jm, km, nx, nxy);
00361               if (sfar[ijkm] != 0)   count++;
00362             }
00363 
00364           /*----- Record voxel to be eroded -----*/
00365           if (count < minimum)  efim[ijk] = (float) sfar[ijk];
00366         }
00367       break;
00368 
00369     case MRI_byte:
00370       for (ijk = 0;  ijk < nxyz;  ijk++)
00371         {       
00372           if (bfar[ijk] == 0)  continue;
00373           IJK_TO_THREE (ijk, iv, jv, kv, nx, nxy);
00374 
00375           /*----- Count number of active voxels in the neighborhood -----*/
00376           count = 0;
00377           for (imask = 0;  imask < nmask;  imask++)
00378             {
00379               im = iv + mask->i[imask];
00380               jm = jv + mask->j[imask];
00381               km = kv + mask->k[imask];
00382               if ( im < 0 || jm < 0 || km < 0 ||
00383                    im >= nx || jm >= ny || km >= nz )  continue;
00384               ijkm = THREE_TO_IJK (im, jm, km, nx, nxy);
00385               if (bfar[ijkm] != 0)   count++;
00386             }
00387 
00388           /*----- Record voxel to be eroded -----*/
00389           if (count < minimum)  efim[ijk] = (float) bfar[ijk];
00390         }
00391       break;
00392 
00393     case MRI_float:
00394       for (ijk = 0;  ijk < nxyz;  ijk++)
00395         {       
00396           if (ffar[ijk] == 0.0)  continue;
00397           IJK_TO_THREE (ijk, iv, jv, kv, nx, nxy);
00398 
00399           /*----- Count number of active voxels in the neighborhood -----*/
00400           count = 0;
00401           for (imask = 0;  imask < nmask;  imask++)
00402             {
00403               im = iv + mask->i[imask];
00404               jm = jv + mask->j[imask];
00405               km = kv + mask->k[imask];
00406               if ( im < 0 || jm < 0 || km < 0 ||
00407                    im >= nx || jm >= ny || km >= nz )  continue;
00408               ijkm = THREE_TO_IJK (im, jm, km, nx, nxy);
00409               if (ffar[ijkm] != 0.0)   count++;
00410             }
00411 
00412           /*----- Record voxel to be eroded -----*/
00413           if (count < minimum)  efim[ijk] = ffar[ijk];
00414         }
00415       break;
00416 
00417     }
00418 
00419 
00420   /*----- Step 2:  Erode voxels -----*/
00421   switch (ftype)
00422     {
00423     case MRI_short:
00424       for (ijk = 0;  ijk < nxyz;  ijk++)
00425         if (efim[ijk] != 0.0)  sfar[ijk] = 0;
00426       break;
00427 
00428     case MRI_byte:
00429       for (ijk = 0;  ijk < nxyz;  ijk++)
00430         if (efim[ijk] != 0.0)  bfar[ijk] = 0;
00431       break;
00432 
00433     case MRI_float:
00434       for (ijk = 0;  ijk < nxyz;  ijk++)
00435         if (efim[ijk] != 0.0)  ffar[ijk] = 0.0;
00436       break;
00437     }
00438 
00439 
00440   /*----- Proceed with dilation phase? -----*/
00441   if (dilate)
00442     {
00443 
00444 
00445       /*----- Step 3:  Identify voxels to be restored -----*/
00446       switch (ftype)
00447         {
00448         case MRI_short:
00449           for (ijk = 0;  ijk < nxyz;  ijk++)
00450             {   
00451               if (efim[ijk] == 0.0)  continue;
00452               IJK_TO_THREE (ijk, iv, jv, kv, nx, nxy);
00453         
00454               /*---- Determine if any active voxels in the neighborhood ----*/
00455               for (imask = 0;  imask < nmask;  imask++)
00456                 {
00457                   im = iv + mask->i[imask];
00458                   jm = jv + mask->j[imask];
00459                   km = kv + mask->k[imask];
00460                   if ( im <  0  || jm <  0  || km <  0 ||
00461                        im >= nx || jm >= ny || km >= nz  )  continue;
00462                   ijkm = THREE_TO_IJK (im, jm, km, nx, nxy);
00463                   if (sfar[ijkm] != 0)  break;
00464                 }
00465         
00466               /*----- Reset voxel not to be restored -----*/
00467               if (imask == nmask)  efim[ijk] = 0.0;
00468             }
00469           break;
00470         
00471         case MRI_byte:
00472           for (ijk = 0;  ijk < nxyz;  ijk++)
00473             {   
00474               if (efim[ijk] == 0.0)  continue;
00475               IJK_TO_THREE (ijk, iv, jv, kv, nx, nxy);
00476         
00477               /*---- Determine if any active voxels in the neighborhood ----*/
00478               for (imask = 0;  imask < nmask;  imask++)
00479                 {
00480                   im = iv + mask->i[imask];
00481                   jm = jv + mask->j[imask];
00482                   km = kv + mask->k[imask];
00483                   if ( im < 0 || jm < 0 || km < 0 ||
00484                        im >= nx || jm >= ny || km >= nz )  continue;
00485                   ijkm = THREE_TO_IJK (im, jm, km, nx, nxy);
00486                   if (bfar[ijkm] != 0)  break;
00487                 }
00488         
00489               /*----- Reset voxel not to be restored -----*/
00490               if (imask == nmask)  efim[ijk] = 0.0;
00491             }
00492           break;
00493         
00494         case MRI_float:
00495           for (ijk = 0;  ijk < nxyz;  ijk++)
00496             {   
00497               if (efim[ijk] == 0.0)  continue;
00498               IJK_TO_THREE (ijk, iv, jv, kv, nx, nxy);
00499         
00500               /*---- Determine if any active voxels in the neighborhood ----*/
00501               for (imask = 0;  imask < nmask;  imask++)
00502                 {
00503                   im = iv + mask->i[imask];
00504                   jm = jv + mask->j[imask];
00505                   km = kv + mask->k[imask];
00506                   if ( im < 0 || jm < 0 || km < 0 ||
00507                        im >= nx || jm >= ny || km >= nz )  continue;
00508                   ijkm = THREE_TO_IJK (im, jm, km, nx, nxy);
00509                   if (ffar[ijkm] != 0.0)  break;
00510                 }
00511         
00512               /*----- Reset voxel not to be restored -----*/
00513               if (imask == nmask)  efim[ijk] = 0.0;
00514             }
00515           break;
00516         }
00517 
00518 
00519       /*----- Step 4:  Restore voxels -----*/
00520       switch (ftype)
00521         {
00522         case MRI_short:
00523           for (ijk = 0;  ijk < nxyz;  ijk++)
00524             if (efim[ijk] != 0.0)  sfar[ijk] = (short) efim[ijk];
00525             break;
00526         
00527         case MRI_byte:
00528           for (ijk = 0;  ijk < nxyz;  ijk++)
00529             if (efim[ijk] != 0.0)  bfar[ijk] = (byte) efim[ijk];
00530             break;
00531 
00532         case MRI_float:
00533           for (ijk = 0;  ijk < nxyz;  ijk++)
00534             if (efim[ijk] != 0.0)  ffar[ijk] = efim[ijk];
00535           break;
00536         }
00537 
00538     }   /*  if (dilate)  */
00539 
00540 
00541   /*----- Release memory -----*/
00542   KILL_CLUSTER(mask) ;
00543   free (efim);   efim = NULL;
00544   EXRETURN ;
00545 }