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 #include "parser.h"
00009 
00010 /*-----------------------------------------------------------------------------
00011    Routine to filter volume data.
00012 
00013       nx, ny, nz   = number of voxels along each axis
00014       dx, dy, dz   = voxel dimensions
00015       fim_type     = volume data type
00016       vfim         = volume data itself
00017       filter_opt   = indicates filtering method to be applied
00018       filter_rmm   = radius of "influence" of neighboring voxels
00019 
00020       fmask        = if non-NULL, is a mask of allowable voxels
00021       fexpr        = character string for FCFLAG_EXPR
00022 
00023    The filtered data is returned in vfim.
00024 
00025    Author :  B. D. Ward
00026    Date   :  11 September 1996
00027    Mod    :  2 October 1996        Changed memory deallocation.
00028    Mod    :  9 October 1996        Made changes to improve speed.
00029 
00030    To correct error due to abiguity in identification of clusters,
00031    voxel coordinates are now stored as 3 separate short integers.
00032    BDW  06 March 1997
00033 
00034    Added fmask and fexpr: RWCox - 09 Aug 2000
00035    Added Winsor filter:   RWCox - 11 Sep 2000
00036 -----------------------------------------------------------------------------*/
00037 
00038 void EDIT_filter_volume (int nx, int ny, int nz, float dx, float dy, float dz,
00039                    int fim_type, void * vfim, int filter_opt, float filter_rmm,
00040                    byte * fmask , char * fexpr )
00041 {
00042    MCW_cluster * mask;                   /* mask for filtering */
00043    int nxy, nxyz;                        /* dimensions of volume data */
00044    int mnum;                             /* number of points in mask */
00045    int i, j, k, ii, jj, kk;              /* voxel indices */
00046    int ijkvox, ijkma, jma;               /* more voxel indices */
00047    float * ffim, * ffim_out;             /* floating point fim's */
00048 
00049    float
00050       mag,                 /* voxel intensity */
00051       sum,                 /* sum of voxel intensities */
00052       sumnz,               /* sum of non-zero voxel intensities */
00053       mean,                /* mean of voxel intensities */
00054       max,                 /* maximum of voxel intensities */
00055       amax,                /* maximum of absolute voxel intensities */
00056       smax;                /* signed maximum of absolute voxel intensities */
00057    int
00058       npts, nznpts;        /* number of points in average */
00059 
00060    float wtsum ;           /* 09 Aug 2000: stuff for FCFLAG_EXPR */
00061    float * wt=NULL ;
00062    PARSER_code * pcode ;
00063 
00064    int nw , nnw , iw ;     /* 11 Sep 2000: Winsor stuff */
00065    float * sw=NULL , vw ;
00066 
00067    nxy = nx*ny;  nxyz = nxy*nz;
00068 
00069 #define GOODVOX(ijk) (fmask==NULL || fmask[ijk]!=0)
00070 #define BADVOX(ijk)  (fmask!=NULL && fmask[ijk]==0)
00071 
00072 ENTRY("EDIT_filter_volume") ;
00073 
00074    /* 09 Aug 2000: can't use AVER code if mask is in place */
00075 
00076    if( fmask != NULL && filter_opt == FCFLAG_AVER ) filter_opt = FCFLAG_MEAN ;
00077 
00078    /***--- 07 Jan 1998 ---***/
00079 
00080    if( filter_opt == FCFLAG_AVER ){
00081       if( fim_type != MRI_float ){
00082          ffim = (float *) malloc (sizeof(float) * nxyz);
00083          if( ffim == NULL ){
00084             fprintf(stderr,"EDIT_filter_volume: no workspace!\n") ;
00085             EXIT(1) ;
00086          }
00087          EDIT_coerce_type (nxyz, fim_type, vfim, MRI_float, ffim);
00088       } else {
00089          ffim = (float *) vfim ;
00090       }
00091       EDIT_aver_fvol(nx,ny,nz,dx,dy,dz,ffim,filter_rmm) ;
00092       if( ffim != vfim ){
00093          EDIT_coerce_autoscale(nxyz, MRI_float, ffim, fim_type, vfim);
00094          free(ffim) ;
00095       }
00096       EXRETURN ;
00097    }
00098 
00099    /***--- end 07 Jan 1998 ---***/
00100 
00101    /*--- Make a cluster that is a mask of points closer than filter_rmm ---*/
00102 
00103    mask = MCW_build_mask (nx, ny, nz, dx, dy, dz, filter_rmm);
00104    if (mask == NULL)
00105    {
00106       fprintf (stderr, "Warning: Filter option has no effect. \n");
00107       EXRETURN;
00108    }
00109    mnum = mask->num_pt;
00110 
00111    /* 09 Aug 2000: evaluate expression weights into wt */
00112 
00113    if( filter_opt == FCFLAG_EXPR ){
00114       double atoz[26] ;
00115 
00116       if( fexpr == NULL ){
00117          fprintf(stderr,"*** EDIT_filter_volume: no fexpr for FCFLAG_EXPR!\n");
00118          EXIT(1) ;
00119       }
00120 
00121       pcode = PARSER_generate_code( fexpr ) ;
00122       if( pcode == NULL ){
00123          fprintf(stderr,"*** EDIT_filter_volume: illegal fexpr!\n"); EXIT(1);
00124       }
00125 
00126       wt = (float *) malloc(sizeof(float)*(mnum+1)) ;
00127 
00128 #define II  8  /* a=0 b=1 ... i=8 ... z=25 */
00129 #define JJ  9
00130 #define KK 10
00131 #define RR 17
00132 #define XX 23
00133 #define YY 24
00134 #define ZZ 25
00135 
00136       for( ii=0 ; ii < 26 ; ii++ ) atoz[ii] = 0.0 ;
00137 
00138       wt[0] = PARSER_evaluate_one( pcode , atoz ) ;  /* weight at center */
00139 
00140       for (jma = 0; jma < mnum; jma++){              /* rest of weights */
00141          atoz[II] = mask->i[jma] ;
00142          atoz[JJ] = mask->j[jma] ;
00143          atoz[KK] = mask->k[jma] ;
00144          atoz[XX] = atoz[II] * dx ;
00145          atoz[YY] = atoz[JJ] * dy ;
00146          atoz[ZZ] = atoz[KK] * dz ;
00147          atoz[RR] = sqrt(atoz[XX]*atoz[XX] + atoz[YY]*atoz[YY] + atoz[ZZ]*atoz[ZZ]) ;
00148          wt[jma+1] = PARSER_evaluate_one( pcode , atoz ) ;
00149       }
00150       free(pcode) ;
00151    }
00152 
00153    /* 11 Sep 2000: setup for Winsorizing */
00154 
00155    if( filter_opt > FCFLAG_WINSOR                 &&
00156        filter_opt < FCFLAG_WINSOR+FCFLAG_ONE_STEP   ){
00157 
00158       static int first=1 ;
00159 
00160       nw  = filter_opt - FCFLAG_WINSOR ;
00161       nnw = mnum - nw ;
00162       filter_opt = FCFLAG_WINSOR ;
00163 
00164       fprintf(stderr,"++ Winsor filter: N=%d nw=%d\n",mnum+1,nw) ;
00165       if( first || nnw < nw ){
00166         first = 0 ;
00167         if( nnw < nw ){
00168           fprintf(stderr,"** Illegal Winsor parameters - skipping!\n") ;
00169           EXRETURN ;
00170         }
00171       }
00172 
00173       sw = (float *) malloc(sizeof(float)*(mnum+1)) ;
00174 
00175       fmask = NULL ;  /* must disable mask for Winsor */
00176    }
00177 
00178    /*--- Allocate space for floating point data ---*/
00179 
00180    ffim = (float *) malloc (sizeof(float) * nxyz);
00181    if (ffim == NULL)
00182    {
00183       fprintf (stderr, "\n Error: cannot allocate filter workspace! \n");
00184       EXIT(1);
00185    }
00186    ffim_out = (float *) malloc (sizeof(float) * nxyz);
00187    if (ffim_out == NULL)
00188    {
00189       fprintf (stderr, "\n Error: cannot allocate filter workspace! \n");
00190       EXIT(1);
00191    }
00192 
00193    /*--- Convert vfim to floating point data ---*/
00194 
00195    EDIT_coerce_type (nxyz, fim_type, vfim, MRI_float, ffim);
00196 
00197    /*--- Iteration over all voxels in vfim ---*/
00198 
00199    for (k = 0; k < nz; k++)
00200    {
00201       for (j = 0; j < ny; j++)
00202       {
00203          for (i = 0; i < nx; i++)
00204          {
00205             /*--- Initialization for filtering of voxel #(i,j,k) ---*/
00206 
00207             ijkvox = THREE_TO_IJK (i, j, k, nx, nxy);
00208             npts = nznpts = 0 ;
00209             sum = sumnz = max = amax = smax = wtsum = 0.0 ;
00210             if( GOODVOX(ijkvox) ){
00211               mag = ffim[ijkvox];
00212               switch (filter_opt)
00213               {
00214                  case FCFLAG_MEAN:
00215                     sum = mag;  npts = 1;  break;
00216                  case FCFLAG_NZMEAN:
00217                     if (mag != 0.0)
00218                        {sumnz = mag;   nznpts = 1;}
00219                     break;
00220                  case FCFLAG_MAX:
00221                     max = mag; npts = 1 ;  break;
00222                  case FCFLAG_AMAX:
00223                     amax = fabs(mag); npts = 1 ;  break;
00224                  case FCFLAG_SMAX:
00225                     smax = mag; npts = 1 ;  break;
00226                  case FCFLAG_EXPR:
00227                     sum = wt[0]*mag ; wtsum = wt[0] ; npts = 1 ; break ;
00228 
00229                  case FCFLAG_WINSOR:
00230                     for( iw=0 ; iw <= mnum ; iw++ ) sw[iw] = mag ;
00231                     vw = mag ; break ;
00232               }
00233             }
00234 
00235             /*--- Now iterate over the positions in the mask ---*/
00236 
00237             switch (filter_opt)
00238               {
00239               case FCFLAG_MEAN:
00240                 for (jma = 0; jma < mnum; jma++)
00241                   {
00242                     ii = i + mask->i[jma];
00243                     jj = j + mask->j[jma];
00244                     kk = k + mask->k[jma];
00245                     if (ii < 0   || jj < 0   || kk < 0 ||
00246                         ii >= nx || jj >= ny || kk >= nz)  continue;
00247                     ijkma = THREE_TO_IJK (ii, jj, kk, nx, nxy);
00248                     if( BADVOX(ijkma) ) continue ;
00249                     mag = ffim[ijkma];
00250                     sum += mag;  npts++;
00251                   }
00252                   break;
00253               case FCFLAG_NZMEAN:
00254                 for (jma = 0; jma < mnum; jma++)
00255                   {
00256                     ii = i + mask->i[jma];
00257                     jj = j + mask->j[jma];
00258                     kk = k + mask->k[jma];
00259                     if (ii < 0   || jj < 0   || kk < 0 ||
00260                         ii >= nx || jj >= ny || kk >= nz)  continue;
00261                     ijkma = THREE_TO_IJK (ii, jj, kk, nx, nxy);
00262                     if( BADVOX(ijkma) ) continue ;
00263                     mag = ffim[ijkma];
00264                     if (mag != 0.0)  {sumnz += mag;  nznpts++;}
00265                   }
00266                 break;
00267               case FCFLAG_EXPR:
00268                 for (jma = 0; jma < mnum; jma++)
00269                   {
00270                     ii = i + mask->i[jma];
00271                     jj = j + mask->j[jma];
00272                     kk = k + mask->k[jma];
00273                     if (ii < 0   || jj < 0   || kk < 0 ||
00274                         ii >= nx || jj >= ny || kk >= nz)  continue;
00275                     ijkma = THREE_TO_IJK (ii, jj, kk, nx, nxy);
00276                     if( BADVOX(ijkma) ) continue ;
00277                     mag = ffim[ijkma];
00278                     sum += wt[jma+1]*mag;  npts++; wtsum += wt[jma+1] ;
00279                   }
00280               break ;
00281               case FCFLAG_MAX:
00282                 for (jma = 0; jma < mnum; jma++)
00283                   {
00284                     ii = i + mask->i[jma];
00285                     jj = j + mask->j[jma];
00286                     kk = k + mask->k[jma];
00287                     if (ii < 0   || jj < 0   || kk < 0 ||
00288                         ii >= nx || jj >= ny || kk >= nz)  continue;
00289                     ijkma = THREE_TO_IJK (ii, jj, kk, nx, nxy);
00290                     if( BADVOX(ijkma) ) continue ;
00291                     mag = ffim[ijkma];
00292                     if (npts == 0 || mag > max)  max = mag;
00293                     npts++ ;
00294                   }
00295                   break;
00296               case FCFLAG_AMAX:
00297                 for (jma = 0; jma < mnum; jma++)
00298                   {
00299                     ii = i + mask->i[jma];
00300                     jj = j + mask->j[jma];
00301                     kk = k + mask->k[jma];
00302                     if (ii < 0   || jj < 0   || kk < 0 ||
00303                         ii >= nx || jj >= ny || kk >= nz)  continue;
00304                     ijkma = THREE_TO_IJK (ii, jj, kk, nx, nxy);
00305                     if( BADVOX(ijkma) ) continue ;
00306                     mag = ffim[ijkma];
00307                     if (npts == 0 || fabs(mag) > amax)  amax = fabs(mag);
00308                     npts++ ;
00309                   }
00310                 break;
00311               case FCFLAG_SMAX:
00312                 for (jma = 0; jma < mnum; jma++)
00313                   {
00314                     ii = i + mask->i[jma];
00315                     jj = j + mask->j[jma];
00316                     kk = k + mask->k[jma];
00317                     if (ii < 0   || jj < 0   || kk < 0 ||
00318                         ii >= nx || jj >= ny || kk >= nz)  continue;
00319                     ijkma = THREE_TO_IJK (ii, jj, kk, nx, nxy);
00320                     if( BADVOX(ijkma) ) continue ;
00321                     mag = ffim[ijkma];
00322                     if (npts == 0 || fabs(mag) > fabs(smax))  smax = mag;
00323                     npts++ ;
00324                   }
00325                 break;
00326               case FCFLAG_WINSOR:
00327                 for (jma = 0; jma < mnum; jma++)
00328                   {
00329                     ii = i + mask->i[jma];
00330                     jj = j + mask->j[jma];
00331                     kk = k + mask->k[jma];
00332                     if (ii < 0   || jj < 0   || kk < 0 ||
00333                         ii >= nx || jj >= ny || kk >= nz)  continue;
00334                     ijkma = THREE_TO_IJK (ii, jj, kk, nx, nxy);
00335                     sw[jma+1] = ffim[ijkma];
00336                   }
00337                   break;
00338 
00339               default:  break;
00340               }
00341 
00342             /*--- Save statistic for this voxel ---*/
00343 
00344             switch (filter_opt)
00345               {
00346               case FCFLAG_MEAN:
00347                 ffim_out[ijkvox] = (npts > 0)     ? sum/npts     : 0.0 ; break;
00348 
00349               case FCFLAG_NZMEAN:
00350                 ffim_out[ijkvox] = (nznpts > 0)   ? sumnz/nznpts : 0.0 ; break ;
00351 
00352               case FCFLAG_EXPR:
00353                 ffim_out[ijkvox] = (wtsum != 0.0) ? sum/wtsum    : 0.0 ; break ;
00354 
00355               case FCFLAG_MAX:  ffim_out[ijkvox] = max ;  break;
00356               case FCFLAG_AMAX: ffim_out[ijkvox] = amax;  break;
00357               case FCFLAG_SMAX: ffim_out[ijkvox] = smax;  break;
00358 
00359               case FCFLAG_WINSOR:
00360                  qsort_float( mnum+1 , sw ) ;
00361                  ffim_out[ijkvox] = (vw < sw[nw]) ? sw[nw]
00362                                                   : (vw > sw[nnw]) ? sw[nnw]
00363                                                                    : vw      ;
00364                  break ;
00365 
00366               default:  break;
00367               }
00368          }  /* i */
00369       }  /* j */
00370    }  /* k */
00371 
00372    /*--- Return the filtered data in the original data type. ---*/
00373    EDIT_coerce_autoscale(nxyz, MRI_float, ffim_out, fim_type, vfim);
00374 
00375    /*--- clean up ---*/
00376    KILL_CLUSTER (mask);
00377    free (ffim);
00378    free (ffim_out);
00379    if( wt != NULL ) free(wt) ;
00380    if( sw != NULL ) free(sw) ;
00381 
00382    EXRETURN;
00383 }
00384 
00385 /*------------------------------------------------------------------------
00386    07 Jan 1998:  A more efficient way to take a local average.
00387                  Creates a larger volume, copies data in, and then
00388                  can average without worrying about edge overrun.
00389 --------------------------------------------------------------------------*/
00390 
00391 void EDIT_aver_fvol( int   nx, int   ny, int   nz,
00392                      float dx, float dy, float dz, float * fim , float rmm )
00393 {
00394    MCW_cluster * mask ;
00395    int i, j, k , ij , ii ;
00396    int jk,jkadd , nxadd,nyadd,nzadd , nxyz_add , mnum ;
00397    float * ffim ;
00398    int * madd ;
00399    float fac , sum ;
00400 
00401 ENTRY("EDIT_aver_fvol") ;
00402 
00403    /*--- Make a cluster that is a mask of points closer than rmm ---*/
00404 
00405    mask = MCW_build_mask(nx,ny,nz, dx,dy,dz, rmm) ;
00406    if( mask == NULL || mask->num_pt < 2 ){
00407       fprintf(stderr,"Warning: EDIT_aver_volume has no effect.\n") ;
00408       EXRETURN ;
00409    }
00410    mnum = mask->num_pt ;
00411 
00412    /*--- Allocate workspaces ---*/
00413 
00414 #if 1
00415    nxadd = nyadd = nzadd = 1 ;
00416    for( ii=0 ; ii < mnum ; ii++ ){
00417       i = abs((int)mask->i[ii]) ; nxadd = MAX(i,nxadd) ;
00418       j = abs((int)mask->j[ii]) ; nyadd = MAX(j,nyadd) ;
00419       k = abs((int)mask->k[ii]) ; nzadd = MAX(k,nzadd) ;
00420    }
00421 #else
00422    nxadd    = (int)(rmm/dx) ;
00423    nyadd    = (int)(rmm/dy) ;
00424    nzadd    = (int)(rmm/dz) ;
00425 #endif
00426 
00427    nxyz_add = (nx+2*nxadd) * (ny+2*nyadd) * (nz+2*nzadd) ;
00428 
00429    ffim = (float *) malloc(sizeof(float) * nxyz_add) ;
00430    if(ffim == NULL){
00431       fprintf(stderr,"*** EDIT_aver_volume can't malloc workspace!\n") ;
00432       fprintf(stderr,"nx=%d ny=%d nz=%d nxadd=%d nyadd=%d nzadd=%d\n",
00433               nx,ny,nz , nxadd,nyadd,nzadd ) ;
00434       EXIT(1) ;
00435    }
00436    for( i=0 ; i < nxyz_add ; i++ ) ffim[i] = 0.0 ;
00437 
00438    madd = (int *) malloc( sizeof(int) * (mnum+1) ) ;
00439    if( madd == NULL ){
00440       fprintf(stderr,"*** EDIT_aver_volume can't malloc workspace!\n") ;
00441       EXIT(1) ;
00442    }
00443    madd[0] = 0 ;
00444    for( ii=0 ; ii < mnum ; ii++ ){
00445       madd[ii+1] = mask->i[ii] +
00446                    mask->j[ii] * (nx+2*nxadd) +
00447                    mask->k[ii] * ((nx+2*nxadd)*(ny+2*nyadd)) ;
00448    }
00449    mnum++ ; fac = 1.0 / mnum ;
00450 
00451    KILL_CLUSTER(mask) ;
00452 
00453    /*-- copy data into workspace --*/
00454 
00455    for( k=0 ; k < nz ; k++ ){
00456       for( j=0 ; j < ny ; j++ ){
00457          jkadd = j * (nx+2*nxadd) + k * ((nx+2*nxadd)*(ny+2*nyadd)) ;
00458          jk    = j * nx + k * (nx * ny) ;
00459          for( i=0 ; i < nx ; i++ ) ffim[i+jkadd] = fim[i+jk] ;
00460       }
00461    }
00462 
00463    /*-- average data from workspace back into original array --*/
00464 
00465    for( k=0 ; k < nz ; k++ ){
00466       for( j=0 ; j < ny ; j++ ){
00467          jkadd = j * (nx+2*nxadd) + k * ((nx+2*nxadd)*(ny+2*nyadd)) ;
00468          jk    = j * nx + k * (nx * ny) ;
00469          for( i=0 ; i < nx ; i++ ){
00470             sum = 0.0 ; ij = i+jkadd ;
00471             for( ii=0 ; ii < mnum ; ii++ ) sum += ffim[ij+madd[ii]] ;
00472             fim[i+jk] = fac * sum ;
00473          }
00474       }
00475    }
00476 
00477    free(ffim); free(madd);
00478    EXRETURN;
00479 }