00001 #include "mrilib.h"
00002 #include "thd_brainormalize.h"
00003 
00004 #undef  IJK
00005 #define IJK(i,j,k) ((i)+(j)*nx+(k)*nxy)
00006 
00007 static int verb = 0 ;
00008 void mri_brainormalize_verbose( int v ){ verb = v ; }
00009 static float thd_bn_dxyz = 1.0;
00010 static int thd_bn_nx     = 167;
00011 static int thd_bn_ny     = 212;
00012 static int thd_bn_nz     = 175;
00013 void mri_brainormalize_initialize(float dx, float dy, float dz)
00014 {
00015    /* Set the reinterpolation resolution to the smallest delta */
00016    thd_bn_dxyz = MIN(fabs(dx), fabs(dy)); thd_bn_dxyz = MIN(thd_bn_dxyz, fabs(dz));
00017    thd_bn_nx = (int)( (float)thd_bn_nx / thd_bn_dxyz );
00018    thd_bn_ny = (int)( (float)thd_bn_ny / thd_bn_dxyz );
00019    thd_bn_nz = (int)( (float)thd_bn_nz / thd_bn_dxyz );
00020    return;
00021 }
00022 float THD_BN_dxyz()
00023 {
00024    return thd_bn_dxyz;
00025 }
00026 int THD_BN_nx()
00027 {
00028    return thd_bn_nx;
00029 }
00030 int THD_BN_ny()
00031 {
00032    return thd_bn_ny;
00033 }
00034 int THD_BN_nz()
00035 {
00036    return thd_bn_nz;
00037 }
00038 
00039 
00040 /*---------------------------------------------------------------------*/
00041 /*! Count number of nonzeros in mask array */
00042 
00043 static int mask_count( int nvox , byte *mmm )
00044 {
00045    int ii , nn ;
00046    for( nn=ii=0 ; ii < nvox ; ii++ ) nn += (mmm[ii] != 0) ;
00047    return nn ;
00048 }
00049 
00050 /*------------------------------------------------------------------*/
00051 
00052 #undef  DALL
00053 #define DALL 4096  /* Allocation size for cluster arrays */
00054 
00055 /*--------------------------------------------------------------------------*/
00056 /*! Put (i,j,k) into the cluster, if it is nonzero. */
00057 
00058 #undef  DPUT
00059 #define DPUT(i,j,k,d)                                               \
00060   do{ ijk = (i)+(j)*nx+(k)*nxy ;                                    \
00061       if( mmm[ijk] == 0 ) break ;                                   \
00062       if( nnow == nall ){ /* increase array lengths */              \
00063         nall += DALL ;                                              \
00064         inow = (short *) realloc((void *)inow,sizeof(short)*nall) ; \
00065         jnow = (short *) realloc((void *)jnow,sizeof(short)*nall) ; \
00066         know = (short *) realloc((void *)know,sizeof(short)*nall) ; \
00067       }                                                             \
00068       inow[nnow] = (i); jnow[nnow] = (j); know[nnow] = (k);         \
00069       nnow++ ; mmm[ijk] = 0 ; ddd[ijk] = (d) ;                      \
00070     } while(0)
00071 
00072 /*--------------------------------------------------------------------------*/
00073 
00074 short * THD_mask_distize( int nx, int ny, int nz, byte *mmm, byte *ccc )
00075 {
00076    short *ddd , dnow ;
00077    int ii,jj,kk , nxy=nx*ny , nxyz=nx*ny*nz , ijk ;
00078    int ip,jp,kp , im,jm,km , icl ;
00079    int nccc,nmmm , nall,nnow ;
00080    short *inow , *jnow , *know ;
00081    float drat ;
00082 
00083    if( mmm == NULL || ccc == NULL ) return NULL ;
00084 
00085    ddd = (short *)malloc( sizeof(short)*nxyz ) ;
00086    nccc = nmmm = 0 ;
00087    for( ii=0 ; ii < nxyz ; ii++ ){
00088           if( ccc[ii] ){ ddd[ii] =  1; nccc++; nmmm++; }
00089      else if( mmm[ii] ){ ddd[ii] = -1; nmmm++; }
00090      else              { ddd[ii] =  0; }
00091    }
00092    if( nccc == 0 ){ free((void *)ddd); return NULL; }
00093 
00094    nall  = nccc+DALL ;                            /* # allocated pts */
00095    inow  = (short *) malloc(sizeof(short)*nall) ; /* coords of pts */
00096    jnow  = (short *) malloc(sizeof(short)*nall) ;
00097    know  = (short *) malloc(sizeof(short)*nall) ;
00098    nnow  = 0 ;
00099 
00100    for( ii=0 ; ii < nxyz ; ii++ ){
00101      if( ccc[ii] ){
00102        inow[nnow] = ii % nx ;
00103        jnow[nnow] = (ii%nxy)/nx ;
00104        know[nnow] = ii / nxy ;
00105        mmm[ii]    = 0 ;
00106        nnow++ ;
00107      }
00108    }
00109 
00110    for( icl=0 ; icl < nnow ; icl++ ){
00111      ii = inow[icl] ; jj = jnow[icl] ; kk = know[icl] ; ijk = ii+jj*nx+kk*nxy ;
00112      im = ii-1      ; jm = jj-1      ; km = kk-1      ;
00113      ip = ii+1      ; jp = jj+1      ; kp = kk+1      ; dnow = ddd[ijk]+1 ;
00114 
00115      if( im >= 0 ) DPUT(im,jj,kk,dnow) ;
00116      if( ip < nx ) DPUT(ip,jj,kk,dnow) ;
00117      if( jm >= 0 ) DPUT(ii,jm,kk,dnow) ;
00118      if( jp < ny ) DPUT(ii,jp,kk,dnow) ;
00119      if( km >= 0 ) DPUT(ii,jj,km,dnow) ;
00120      if( kp < nz ) DPUT(ii,jj,kp,dnow) ;
00121    }
00122 
00123    for( ii=0 ; ii < nxyz ; ii++ ) mmm[ii] = (ddd[ii] > 0) ;
00124 
00125    free((void *)inow); free((void *)jnow); free((void *)know);
00126 
00127    return ddd ;
00128 }
00129 
00130 /*--------------------------------------------------------------------------*/
00131 /*! Put (i,j,k) into the current cluster, if it is nonzero. */
00132 
00133 #undef  CPUT
00134 #define CPUT(i,j,k)                                                   \
00135   do{ ijk = (i)+(j)*nx+(k)*nxy ;                                      \
00136       if( mmm[ijk] ){                                                 \
00137         if( nnow == nall ){ /* increase array lengths */              \
00138           nall += DALL ;                                              \
00139           inow = (short *) realloc((void *)inow,sizeof(short)*nall) ; \
00140           jnow = (short *) realloc((void *)jnow,sizeof(short)*nall) ; \
00141           know = (short *) realloc((void *)know,sizeof(short)*nall) ; \
00142         }                                                             \
00143         inow[nnow] = (i); jnow[nnow] = (j); know[nnow] = (k);         \
00144         nnow++ ; mmm[ijk] = 0 ;                                       \
00145       } } while(0)
00146 
00147 /*--------------------------------------------------------------------------*/
00148 /*! Remove clusters below csize from the binary mask mmm. */
00149 
00150 static void clustedit3D( int nx, int ny, int nz, byte *mmm, int csize )
00151 {
00152    int ii,jj,kk, icl , nxy,nxyz , ijk , ijk_last ;
00153    int ip,jp,kp , im,jm,km ;
00154    int nnow , nall , nsav , nkill ;
00155    short *inow , *jnow , *know ;
00156    short *isav , *jsav , *ksav ;
00157 
00158    if( nx < 1 || ny < 1 || nz < 1 || mmm == NULL || csize < 2 ) return ;
00159 
00160    nxy = nx*ny ; nxyz = nxy*nz ;
00161 
00162    nall  = 8 ;                                    /* # allocated pts */
00163    inow  = (short *) malloc(sizeof(short)*nall) ; /* coords of pts */
00164    jnow  = (short *) malloc(sizeof(short)*nall) ;
00165    know  = (short *) malloc(sizeof(short)*nall) ;
00166 
00167    nsav  = nkill = 0 ;
00168    isav  = (short *) malloc(sizeof(short)) ;
00169    jsav  = (short *) malloc(sizeof(short)) ;
00170    ksav  = (short *) malloc(sizeof(short)) ;
00171 
00172    /*--- scan through array, find nonzero point, build a cluster, ... ---*/
00173 
00174 #if 0
00175    if( verb ) fprintf(stderr," + clustedit3D: threshold size=%d voxels\n",csize) ;
00176 #endif
00177 
00178    ijk_last = 0 ;
00179    while(1) {
00180      /* find next nonzero point */
00181 
00182      for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( mmm[ijk] ) break ;
00183      if( ijk == nxyz ) break ;  /* didn't find any! */
00184 
00185      ijk_last = ijk+1 ;         /* start here next time */
00186 
00187      /* init current cluster list with this point */
00188 
00189      mmm[ijk] = 0 ;                                /* clear found point */
00190      nnow     = 1 ;                                /* # pts in cluster */
00191      inow[0]  = ijk % nx ;
00192      jnow[0]  = (ijk%nxy)/nx ;
00193      know[0]  = ijk / nxy ;
00194 
00195      /*--
00196         for each point in cluster:
00197            check 6 neighboring points for nonzero entries in mmm
00198            enter those into cluster (and clear them in mmm)
00199            continue until end of cluster is reached
00200            (note that cluster size nnow is expanding as we progress)
00201      --*/
00202 
00203      for( icl=0 ; icl < nnow ; icl++ ){
00204        ii = inow[icl] ; jj = jnow[icl] ; kk = know[icl] ;
00205        im = ii-1      ; jm = jj-1      ; km = kk-1      ;
00206        ip = ii+1      ; jp = jj+1      ; kp = kk+1      ;
00207 
00208        if( im >= 0 ) CPUT(im,jj,kk) ;
00209        if( ip < nx ) CPUT(ip,jj,kk) ;
00210        if( jm >= 0 ) CPUT(ii,jm,kk) ;
00211        if( jp < ny ) CPUT(ii,jp,kk) ;
00212        if( km >= 0 ) CPUT(ii,jj,km) ;
00213        if( kp < nz ) CPUT(ii,jj,kp) ;
00214      }
00215 
00216      /* see if now cluster is large enough to live */
00217 
00218      if( nnow >= csize ){
00219        kk = nsav+nnow ;
00220        isav = (short *)realloc((void *)isav,sizeof(short)*kk) ;
00221        jsav = (short *)realloc((void *)jsav,sizeof(short)*kk) ;
00222        ksav = (short *)realloc((void *)ksav,sizeof(short)*kk) ;
00223        memcpy(isav+nsav,inow,sizeof(short)*nnow) ;
00224        memcpy(jsav+nsav,jnow,sizeof(short)*nnow) ;
00225        memcpy(ksav+nsav,know,sizeof(short)*nnow) ;
00226        nsav = kk ;
00227 #if 0
00228        if( verb )
00229          fprintf(stderr," + clustedit3D: saved cluster with %d voxels\n",nnow);
00230 #endif
00231      } else {
00232        nkill += nnow ;
00233      }
00234 
00235    } /* loop ends when all nonzero points are clustered */
00236 
00237    free((void *)inow); free((void *)jnow); free((void *)know);
00238 
00239    /* copy saved points back into mmm */
00240 
00241    for( ii=0 ; ii < nsav ; ii++ )
00242      mmm[ IJK(isav[ii],jsav[ii],ksav[ii]) ] = 1 ;
00243 
00244    free((void *)isav); free((void *)jsav); free((void *)ksav) ;
00245 
00246 #if 0
00247    if( verb )
00248      fprintf(stderr," + clustedit3D totals:"
00249                     " saved=%d killed=%d nxyz=%d\n",nsav,nkill,nxyz) ;
00250 #endif
00251    return ;
00252 }
00253 
00254 /*--------------------------------------------------------------------------*/
00255 /*! Find the cliplevel in just part of the image,
00256     with i=ibot..itop (inclusive), et cetera.     */
00257 
00258 static float partial_cliplevel( MRI_IMAGE *im , float mfrac ,
00259                                 int ibot,int itop ,
00260                                 int jbot,int jtop , int kbot,int ktop )
00261 {
00262    int ncut,ib , qq,nold ;
00263    int ii,jj,kk , nx,ny,nz,nxy ;
00264    int *hist ;
00265    short *sar ;
00266    byte *bar ;
00267    int nhist , npos , nhalf , val ;
00268    double dsum ;
00269 
00270 ENTRY("partial_cliplevel") ;
00271    if( im == NULL || im->kind != MRI_short ) RETURN(1.0) ;
00272 
00273    if( mfrac <= 0.0 || mfrac >= 0.99 ) mfrac = 0.50 ;
00274 
00275    nx = im->nx ; ny = im->ny ; nz = im->nz ; nxy = nx*ny ;
00276 
00277    if( ibot <  0  ) ibot = 0 ;
00278    if( jbot <  0  ) jbot = 0 ;
00279    if( kbot <  0  ) kbot = 0 ;
00280    if( itop >= nx ) itop = nx-1 ;
00281    if( jtop >= ny ) jtop = ny-1 ;
00282    if( ktop >= nz ) ktop = nz-1 ;
00283 
00284    if( itop < ibot || jtop < jbot || ktop < kbot ) RETURN(1.0) ;
00285 
00286    /*-- allocate histogram --*/
00287 
00288    nhist = 32767 ;
00289    hist  = (int *) calloc(sizeof(int),nhist+1) ;
00290 
00291    /*-- make histogram of positive entries --*/
00292 
00293    dsum = 0.0 ;  /* will be sum of squares */
00294    npos = 0 ;    /* will be number of positive values */
00295    sar  =  MRI_SHORT_PTR(im) ;
00296    for( kk=kbot ; kk <= ktop ; kk++ ){
00297     for( jj=jbot ; jj <= jtop ; jj++ ){
00298      for( ii=ibot ; ii <= itop ; ii++ ){
00299        val = sar[IJK(ii,jj,kk)] ;
00300        if( val > 0 && val <= nhist ){
00301          hist[val]++ ;
00302          dsum += (double)(val)*(double)(val); npos++;
00303        }
00304    }}}
00305 
00306    if( npos <= 999 ){ free((void *)hist); RETURN(1.0); }
00307 
00308    /*-- initialize cut position to include upper 65% of positive voxels --*/
00309 
00310    qq = 0.65 * npos ; ib = rint(0.5*sqrt(dsum/npos)) ;
00311    for( kk=0,ii=nhist-1 ; ii >= ib && kk < qq ; ii-- ) kk += hist[ii] ;
00312 
00313    /*-- algorithm --*/
00314 
00315    ncut = ii ;   /* initial cut position */
00316    qq   = 0 ;    /* iteration count */
00317    do{
00318      for( npos=0,ii=ncut ; ii < nhist ; ii++ ) npos += hist[ii]; /* number >= cut */
00319      nhalf = npos/2 ;
00320      for( kk=0,ii=ncut ; ii < nhist && kk < nhalf ; ii++ )  /* find median of */
00321        kk += hist[ii] ;                                     /* valuess >= cut */
00322      nold = ncut ;                                          /* last cut */
00323      ncut = mfrac * ii ;                                    /* new cut */
00324      qq++ ;
00325   } while( qq < 20 && ncut != nold ) ; /* iterate until done, or at most 20 times */
00326 
00327    free((void *)hist) ;
00328    RETURN( (float)(ncut) ) ;
00329 }
00330 
00331 /*--------------------------------------------------------------------------*/
00332 
00333 typedef struct {
00334    float clip_000, clip_100, clip_010, clip_110,
00335          clip_001, clip_101, clip_011, clip_111 ;
00336    float x0,x1,dxi , y0,y1,dyi , z0,z1,dzi ;
00337    float clip_min , clip_max ;
00338 } clipvec ;
00339 
00340 /*! Get the cliplevel for each octant about the center-of-mass. */
00341 
00342 static clipvec get_octant_clips( MRI_IMAGE *im , float mfrac )
00343 {
00344    float xcm,ycm,zcm , sum,val , clip_min , clip_max ;
00345    int ii,jj,kk , nx,ny,nz,nxy , ic,jc,kc , it,jt,kt , ijk ;
00346    short *sar ;
00347    clipvec cv ;
00348 
00349 ENTRY("get_octant_clips") ;
00350 
00351    cv.clip_000 = -1 ;  /* flags error return */
00352 
00353    if( im == NULL || im->kind != MRI_short ) RETURN(cv) ;
00354 
00355    nx = im->nx ; ny = im->ny ; nz = im->nz ; nxy = nx*ny ;
00356    it = nx-1   ; jt = ny-1   ; kt = nz-1   ;
00357 
00358    /* compute CM of image */
00359 
00360    sar = MRI_SHORT_PTR(im) ; if( sar == NULL ) RETURN(cv) ;
00361 
00362    xcm = ycm = zcm = sum = 0.0 ;
00363    for( ijk=kk=0 ; kk < nz ; kk++ ){
00364     for( jj=0 ; jj < ny ; jj++ ){
00365      for( ii=0 ; ii < nx ; ii++,ijk++ ){
00366        val = (float)sar[ijk]; if( val <= 0.0 ) continue;
00367        sum += val ;
00368        xcm += val * ii ;
00369        ycm += val * jj ;
00370        zcm += val * kk ;
00371    }}}
00372    if( sum == 0.0 ) RETURN(cv) ;
00373    ic = (int)rint(xcm/sum); jc = (int)rint(ycm/sum); kc = (int)rint(zcm/sum);
00374 
00375    /* compute cliplevel in each octant about the CM */
00376 
00377    val = 0.5 * partial_cliplevel( im,mfrac , 0,it , 0,jt , 0,kt ) ;
00378 
00379    cv.clip_000 = partial_cliplevel( im,mfrac,  0  ,ic+2,  0  ,jc+2,  0  ,kc+2 );
00380    cv.clip_100 = partial_cliplevel( im,mfrac, ic-2,it  ,  0  ,jc+2,  0  ,kc+2 );
00381    cv.clip_010 = partial_cliplevel( im,mfrac,  0  ,ic+2, jc-2,jt  ,  0  ,kc+2 );
00382    cv.clip_110 = partial_cliplevel( im,mfrac, ic-2,it  , jc-2,jt  ,  0  ,kc+2 );
00383    cv.clip_001 = partial_cliplevel( im,mfrac,  0  ,ic+2,  0  ,jc+2, kc-2,kt   );
00384    cv.clip_101 = partial_cliplevel( im,mfrac, ic-2,it  ,  0  ,jc+2, kc-2,kt   );
00385    cv.clip_011 = partial_cliplevel( im,mfrac,  0  ,ic+2, jc-2,jt  , kc-2,kt   );
00386    cv.clip_111 = partial_cliplevel( im,mfrac, ic-2,it  , jc-2,jt  , kc-2,kt   );
00387 
00388    if( cv.clip_000 < val ) cv.clip_000 = val ;  /* don't let   */
00389    if( cv.clip_100 < val ) cv.clip_100 = val ;  /* clip levels */
00390    if( cv.clip_010 < val ) cv.clip_010 = val ;  /* get too     */
00391    if( cv.clip_110 < val ) cv.clip_110 = val ;  /* small!      */
00392    if( cv.clip_001 < val ) cv.clip_001 = val ;
00393    if( cv.clip_101 < val ) cv.clip_101 = val ;
00394    if( cv.clip_011 < val ) cv.clip_011 = val ;
00395    if( cv.clip_111 < val ) cv.clip_111 = val ;
00396 
00397    clip_min =              cv.clip_000  ;
00398    clip_min = MIN(clip_min,cv.clip_100) ;
00399    clip_min = MIN(clip_min,cv.clip_010) ;
00400    clip_min = MIN(clip_min,cv.clip_110) ;
00401    clip_min = MIN(clip_min,cv.clip_001) ;
00402    clip_min = MIN(clip_min,cv.clip_101) ;
00403    clip_min = MIN(clip_min,cv.clip_011) ;
00404    clip_min = MIN(clip_min,cv.clip_111) ;  cv.clip_min = clip_min ;
00405 
00406    clip_max =              cv.clip_000  ;
00407    clip_max = MAX(clip_max,cv.clip_100) ;
00408    clip_max = MAX(clip_max,cv.clip_010) ;
00409    clip_max = MAX(clip_max,cv.clip_110) ;
00410    clip_max = MAX(clip_max,cv.clip_001) ;
00411    clip_max = MAX(clip_max,cv.clip_101) ;
00412    clip_max = MAX(clip_max,cv.clip_011) ;
00413    clip_max = MAX(clip_max,cv.clip_111) ;  cv.clip_max = clip_max ;
00414 
00415    /* (x0,y0,z0) = center of lowest octant
00416       (x1,y1,z1) = center of highest octant */
00417 
00418    cv.x0  = 0.5*ic ; cv.x1 = 0.5*(ic+it) ;
00419    cv.y0  = 0.5*jc ; cv.y1 = 0.5*(jc+jt) ;
00420    cv.z0  = 0.5*kc ; cv.z1 = 0.5*(kc+kt) ;
00421    cv.dxi = (cv.x1 > cv.x0) ? 1.0/(cv.x1-cv.x0) : 0.0 ;
00422    cv.dyi = (cv.y1 > cv.y0) ? 1.0/(cv.y1-cv.y0) : 0.0 ;
00423    cv.dzi = (cv.z1 > cv.z0) ? 1.0/(cv.z1-cv.z0) : 0.0 ;
00424 
00425 #if 0
00426    if( verb )
00427     fprintf(stderr," + get_octant_clips:  min clip=%.1f\n"
00428                    "   clip_000=%.1f  clip_100=%.1f  clip_010=%.1f  clip_110=%.1f\n"
00429                    "   clip_001=%.1f  clip_101=%.1f  clip_011=%.1f  clip_111=%.1f\n"
00430                    "   (x0,y0,z0)=(%.1f,%.1f,%.1f) (x1,y1,z1)=(%.1f,%.1f,%.1f)\n" ,
00431             val ,
00432             cv.clip_000 , cv.clip_100 , cv.clip_010 , cv.clip_110 ,
00433             cv.clip_001 , cv.clip_101 , cv.clip_011 , cv.clip_111 ,
00434             cv.x0 , cv.y0 , cv.z0 , cv.x1 , cv.y1 , cv.z1  ) ;
00435 #endif
00436 
00437    RETURN(cv) ;
00438 }
00439 
00440 /*--------------------------------------------------------------------------*/
00441 /*! Return the cliplevel at any point,
00442     tri-linearly interpolated between octant centers. */
00443 
00444 static INLINE float pointclip( int ii, int jj, int kk , clipvec *cv )
00445 {
00446    float x1,y1,z1 , x0,y0,z0 , val ;
00447 
00448    /* get relative position in box defined by octant centers */
00449 
00450    x1 = (ii-cv->x0)*cv->dxi; if(x1 < 0.0) x1=0.0; else if(x1 > 1.0) x1=1.0;
00451    y1 = (jj-cv->y0)*cv->dyi; if(y1 < 0.0) y1=0.0; else if(y1 > 1.0) y1=1.0;
00452    z1 = (kk-cv->z0)*cv->dzi; if(z1 < 0.0) z1=0.0; else if(z1 > 1.0) z1=1.0;
00453 
00454    x0 = 1.0-x1 ; y0 = 1.0-y1 ; z0 = 1.0-z1 ;
00455 
00456    val =  cv->clip_000 * x0*y0*z0 + cv->clip_100 * x1*y0*z0
00457         + cv->clip_010 * x0*y1*z0 + cv->clip_110 * x1*y1*z0
00458         + cv->clip_001 * x0*y0*z1 + cv->clip_101 * x1*y0*z1
00459         + cv->clip_011 * x0*y1*z1 + cv->clip_111 * x1*y1*z1 ;
00460    return val ;
00461 }
00462 
00463 /*--------------------------------------------------------------------------*/
00464 /*! Make a mask from the 3D image, based on gradually changing cliplevels.
00465     This feature is to allow for gradual shading in the MRI volume.        */
00466 
00467 static byte * mri_short2mask( MRI_IMAGE *im )
00468 {
00469    int ii,jj,kk,ijk , nx,ny,nz,nxy,nxyz ;
00470    clipvec bvec , tvec ;
00471    short *sar ;
00472    byte *mask ;
00473    float bval , tval ;
00474 
00475 ENTRY("mri_short2mask") ;
00476    if( im == NULL || im->kind != MRI_short ) RETURN(NULL) ;
00477    sar = MRI_SHORT_PTR(im) ; if( sar == NULL ) RETURN(NULL) ;
00478 
00479    nx = im->nx ; ny = im->ny ; nz = im->nz ; nxy = nx*ny ; nxyz = nxy*nz ;
00480 
00481    bvec = get_octant_clips( im , 0.40f ) ;
00482    if( bvec.clip_000 < 0.0 ) RETURN(NULL) ;
00483 
00484    tvec = bvec ;
00485    tvec.clip_000 *= 9.91 ;
00486    tvec.clip_100 *= 9.91 ;
00487    tvec.clip_010 *= 9.91 ;
00488    tvec.clip_110 *= 9.91 ;
00489    tvec.clip_001 *= 9.91 ;
00490    tvec.clip_101 *= 9.91 ;
00491    tvec.clip_011 *= 9.91 ;
00492    tvec.clip_111 *= 9.91 ;
00493 
00494    /* create mask, clipping at a level that varies spatially */
00495 
00496 #if 0
00497    if( verb ) fprintf(stderr," + mri_short2mask: clipping\n") ;
00498 #endif
00499 
00500    mask = (byte *) malloc( sizeof(byte)*nxyz ) ;
00501    for( ijk=kk=0 ; kk < nz ; kk++ ){
00502     for( jj=0 ; jj < ny ; jj++ ){
00503      for( ii=0 ; ii < nx ; ii++,ijk++ ){
00504        bval = pointclip( ii,jj,kk , &bvec ) ; /* cliplevel here */
00505 #if 0
00506        tval = pointclip( ii,jj,kk , &tvec ) ; /* cliplevel here */
00507        mask[ijk] = (sar[ijk] >= bval && sar[ijk] <= tval) ; /* binarize */
00508 #else
00509        mask[ijk] = (sar[ijk] >= bval) ; /* binarize */
00510 #endif
00511    }}}
00512 
00513    /* remove small clusters */
00514 
00515    clustedit3D( nx,ny,nz , mask , (int)rint(0.02*nxyz) ) ;
00516 
00517    if( verb ) fprintf(stderr," + mri_short2mask: %d voxels survive\n",
00518                              mask_count(nxyz,mask) ) ;
00519 
00520    RETURN(mask) ;
00521 }
00522 
00523 /*--------------------------------------------------------------------------*/
00524 
00525 typedef struct { short i,j,k,val; int basin; } shortvox ;
00526 
00527 /*--------------------------------------------------------------------------*/
00528 /*! Sort array of shortvox into increasing order (decreasing if dec != 0). */
00529 
00530 static sort_shortvox( int n , shortvox *ar , int dec , float botperc, float topperc )
00531 {
00532    int ii,jj , sbot,stop,nsv , sval , pbot,ptop ;
00533    int *hsv , *csv ;
00534    shortvox *tar ;
00535 
00536    if( n < 2 || ar == NULL ) return ;
00537 
00538    /* decreasing order desired?  flip values */
00539 
00540    if( dec ){
00541      float tmp ;
00542      for( ii=0 ; ii < n ; ii++ ) ar[ii].val = -ar[ii].val ;
00543      tmp = botperc ; botperc = topperc ; topperc = tmp ;
00544    }
00545 
00546    /* find range of values */
00547 
00548    sbot = stop = ar[0].val ;
00549    for( ii=1 ; ii < n ; ii++ ){
00550      sval = ar[ii].val ;
00551           if( sval < sbot ) sbot = sval ;
00552      else if( sval > stop ) stop = sval ;
00553    }
00554    nsv = stop-sbot+1 ;   /* number of distinct values */
00555    if( nsv <= 1 ) return ;
00556 
00557    /* build hsv[i] = how many have value = sbot+i
00558             csv[i] = how many have value < sbot+i, i=0..nsv-1 */
00559 
00560    hsv = (int *)calloc(sizeof(int),nsv) ;
00561    csv = (int *)calloc(sizeof(int),nsv+1) ;
00562    for( ii=0 ; ii <  n   ; ii++ ) hsv[ar[ii].val-sbot]++ ;
00563    for( ii=1 ; ii <= nsv ; ii++ ) csv[ii] = csv[ii-1]+hsv[ii-1] ;
00564    free((void *)hsv) ;
00565 
00566    if( botperc > 0.0 && botperc < 50.0 ){
00567      jj = (int)rint(0.01*botperc*n) ;
00568      for( ii=0 ; ii < nsv && csv[ii] <= jj ; ii++ ) ;
00569      pbot = ii+sbot ;
00570      if( verb ) fprintf(stderr," + sort_shortvox: sbot=%d pbot=%d\n",sbot,pbot) ;
00571    } else {
00572      pbot = sbot ;
00573    }
00574 
00575    if( topperc > 0.0 && topperc < 50.0 ){
00576      jj = (int)rint(0.01*(100.0-topperc)*n) ;
00577      for( ii=0 ; ii < nsv && csv[ii] <= jj ; ii++ ) ;
00578      ptop = ii+sbot ;
00579      if( verb ) fprintf(stderr," + sort_shortvox: stop=%d ptop=%d\n",stop,ptop) ;
00580    } else {
00581      ptop = stop ;
00582    }
00583 
00584    /* copy from ar into temp array tar,
00585       putting each one into its place as given by csv */
00586 
00587    tar = (shortvox *)calloc(sizeof(shortvox),n) ;
00588    for( ii=0 ; ii < n ; ii++ ){
00589      sval = ar[ii].val - sbot ;   /* sval is in 0..nsv-1 now */
00590      tar[ csv[sval] ] = ar[ii] ;
00591      csv[sval]++ ;
00592    }
00593 
00594    if( pbot > sbot ){
00595      for( ii=0 ; ii < n ; ii++ )
00596        if( tar[ii].val < pbot ) tar[ii].val = pbot ;
00597    }
00598    if( ptop < stop ){
00599      for( ii=0 ; ii < n ; ii++ )
00600        if( tar[ii].val > ptop ) tar[ii].val = ptop ;
00601    }
00602 
00603    /* copy back into ar */
00604 
00605    memcpy( ar , tar , sizeof(shortvox)*n ) ;
00606    free((void *)tar) ; free((void *)csv) ;
00607 
00608    /* unflip? */
00609 
00610    if( dec )
00611      for( ii=0 ; ii < n ; ii++ ) ar[ii].val = -ar[ii].val ;
00612 
00613    return ;
00614 }
00615 
00616 /*--------------------------------------------------------------------------*/
00617 
00618 typedef struct { int num, nall, depth, *ivox; } basin ;
00619 
00620 #define DBALL 32768
00621 
00622 #define BDEP(i) (baslist[i]->depth)
00623 
00624 #define INIT_BASIN(iv)                                       \
00625  { register int qb=nbtop;                                    \
00626    if( qb >= nball ){                                        \
00627      register int qqb=nball+DBALL,zb ;                       \
00628      baslist = (basin **)realloc((void *)baslist,            \
00629                                  sizeof(basin *)*qqb) ;      \
00630      for( zb=nball ; zb < qqb ; zb++ ) baslist[zb] = NULL ;  \
00631      nball = qqb ;                                           \
00632    }                                                         \
00633    baslist[qb] = (basin *) malloc(sizeof(basin)) ;           \
00634    baslist[qb]->num     = 1 ;                                \
00635    baslist[qb]->nall    = 1 ;                                \
00636    baslist[qb]->depth   = svox[iv].val ;                     \
00637    baslist[qb]->ivox    = (int *)malloc(sizeof(int)) ;       \
00638    baslist[qb]->ivox[0] = (iv) ;                             \
00639    svox[iv].basin       = qb ; nbtop++ ;                     \
00640  }
00641 
00642 #define KILL_BASIN(ib)                                       \
00643  { if( baslist[ib] != NULL ){                                \
00644      free((void *)baslist[ib]->ivox) ;                       \
00645      free((void *)baslist[ib]) ;                             \
00646      baslist[ib] = NULL ; }                                  \
00647  }
00648 
00649 #define ADDTO_BASIN(ib,iv)                                   \
00650  { register basin *bb = baslist[ib] ;                        \
00651    if( bb->num == bb->nall ){                                \
00652      bb->nall = (int)(1.2*bb->nall)+32 ;                     \
00653      bb->ivox = (int *)realloc( (void *)bb->ivox ,           \
00654                                  sizeof(int)*bb->nall ) ;    \
00655    }                                                         \
00656    bb->ivox[bb->num] = (iv) ; bb->num++ ;                    \
00657    svox[iv].basin = (ib) ; }
00658 
00659 #define MERGE_BASIN(ib,ic)                                   \
00660  { register basin *bb = baslist[ib], *cc = baslist[ic] ;     \
00661    int zz = bb->num + cc->num ;                              \
00662    if( bb->nall < zz ){                                      \
00663      bb->nall = zz+1 ;                                       \
00664      bb->ivox = (int *)realloc( (void *)bb->ivox ,           \
00665                                   sizeof(int)*bb->nall ) ;   \
00666    }                                                         \
00667    memcpy( bb->ivox + bb->num ,                              \
00668            cc->ivox , sizeof(int) * cc->num ) ;              \
00669    bb->num = zz ;                                            \
00670    for( zz=0 ; zz < cc->num ; zz++ )                         \
00671      svox[cc->ivox[zz]].basin = (ib) ;                       \
00672    KILL_BASIN(ic) ; }
00673 
00674 /*--------------------------------------------------------------------------*/
00675 
00676 MRI_IMAGE * mri_watershedize( MRI_IMAGE *sim , float prefac )
00677 {
00678    MRI_IMAGE *tim ;
00679    int ii,jj,kk , pp,qq , nx,ny,nz,nxy,nxyz , nvox ;
00680    int ip,jp,kp , im,jm,km ;
00681    short *sar , *tar ;
00682    shortvox *svox ;
00683    int *isvox , *bcount,*bname ;
00684    int nb,vb,mb,m,mu,mq,mz , bp[6] , hpf ;
00685 
00686    basin **baslist ;
00687    int nball , nbtop ;
00688 
00689 ENTRY("watershedize") ;
00690 
00691    if( sim == NULL || sim->kind != MRI_short ) RETURN(NULL) ;
00692    sar = MRI_SHORT_PTR(sim) ; if( sar == NULL ) RETURN(NULL) ;
00693 
00694    nx = sim->nx; ny = sim->ny; nz = sim->nz; nxy = nx*ny; nxyz = nxy*nz;
00695 
00696    /* count number of voxels > 0 */
00697 
00698    for( nvox=0,pp=0 ; pp < nxyz ; pp++ ) if( sar[pp] > 0 ) nvox++ ;
00699    if( nvox <= 999 ) RETURN(NULL) ;
00700 
00701    if( verb ) fprintf(stderr," + mri_watershedize: %d voxels input\n",nvox) ;
00702 
00703    /* create voxel lists */
00704 
00705    svox  = (shortvox *) malloc( sizeof(shortvox)* nvox ) ;
00706    isvox = (int *)      malloc( sizeof(int)     * nxyz ) ;
00707    for( qq=pp=0 ; pp < nxyz ; pp++ ){
00708      if( sar[pp] > 0 ){                  /* save this one: */
00709        ii             = pp % nx ;        /* spatial indexes */
00710        jj             = (pp%nxy) / nx ;
00711        kk             = pp / nxy ;
00712        svox[qq].i     = ii ;
00713        svox[qq].j     = jj ;
00714        svox[qq].k     = kk ;
00715        svox[qq].val   = sar[pp] ;        /* value */
00716        svox[qq].basin = -1 ;             /* which basin */
00717        qq++ ;
00718        isvox[pp]      = qq ;             /* where in list */
00719      } else {
00720        isvox[pp] = -1 ;                  /* voxel not in list */
00721      }
00722    }
00723 
00724    /* sort voxel list into descending order */
00725 
00726    if( verb ) fprintf(stderr," + mri_watershedize: sorting voxels\n") ;
00727 
00728    sort_shortvox( nvox , svox , 1 , 0.00 , 0.02 ) ;
00729 
00730    /* create basin for first (deepest) voxel */
00731 
00732    nball    = DBALL ;
00733    nbtop    = 0 ;
00734    baslist  = (basin **) calloc(sizeof(basin *),nball) ;
00735 
00736    INIT_BASIN(0) ;
00737 
00738    hpf      = (int)rint(prefac*svox[0].val) ;      /* preflood */
00739 
00740    /* scan voxels as they get shallower, and basinate them */
00741 
00742    if( verb ){
00743      fprintf(stderr," + mri_watershedize: basinating voxels\n") ;
00744      fprintf(stderr,"  data range: %d..%d preflood_height=%d\n",
00745              svox[nvox-1].val , svox[0].val , hpf ) ;
00746    }
00747 
00748    for( pp=1 ; pp < nvox ; pp++ ){
00749 
00750      ii = svox[pp].i; jj = svox[pp].j; kk = svox[pp].k;  /* where */
00751      ip = ii+1 ; jp = jj+1 ; kp = kk+1 ;                 /* nbhrs */
00752      im = ii-1 ; jm = jj-1 ; km = kk-1 ;
00753 
00754      if( verb && pp%100000 == 0 ) fprintf(stderr, (pp%1000000)?".":"!") ;
00755 
00756      /* macro checks if (a,b,c) voxel is in the list;
00757         if so and it is already in a basin, then
00758         make a list of basins encountered:
00759           nb = number of unique basins encountered (0..6)
00760           mb = index of deepest basin encountered (0..nb-1)
00761           vb = value (depth) of deepest basin encountered
00762           bp[m] = index of m-th basin encountered (m=0..nb-1) */
00763 
00764 #undef  BASECHECK
00765 #define BASECHECK(a,b,c)                                   \
00766     { qq = isvox[IJK(a,b,c)] ;                             \
00767       if( qq >= 0 && svox[qq].basin >= 0 ){                \
00768         qq = svox[qq].basin ;                              \
00769         for( m=0 ; m < nb && bp[m] != qq ; m++ ) ;         \
00770         if( m == nb ){                                     \
00771           bp[nb] = qq ;                                    \
00772           if( BDEP(qq) > vb ){ mb = nb; vb = BDEP(qq); }   \
00773           nb++ ;                                           \
00774         }                                                  \
00775       }                                                    \
00776     }
00777 
00778      nb = 0 ; vb = -1 ; mb = -1 ;         /* initialize counters */
00779      if( ip < nx ) BASECHECK(ip,jj,kk) ;  /* check each neighbor */
00780      if( im >= 0 ) BASECHECK(im,jj,kk) ;  /* for basin-ositiness */
00781      if( jp < ny ) BASECHECK(ii,jp,kk) ;
00782      if( jm >= 0 ) BASECHECK(ii,jm,kk) ;
00783      if( kp < nz ) BASECHECK(ii,jj,kp) ;
00784      if( km >= 0 ) BASECHECK(ii,jj,km) ;
00785 
00786      if( nb == 0 ){  /*** this voxel is isolated ==> create new basin ****/
00787 
00788        INIT_BASIN(pp) ;
00789 
00790      } else {        /*** this voxel has deeper neighbors ***/
00791 
00792        mq = bp[mb] ;                      /* assign voxel to best basin */
00793        ADDTO_BASIN( mq , pp ) ;
00794 
00795                        /* if have more than one neighbor, other */
00796        if( nb > 1 ){   /* basins could be merged with the best  */
00797          mz = svox[pp].val ;          /* depth of this voxel */
00798          for( m=0 ; m < nb ; m++ ){
00799            if( m == mb ) continue ;        /* can't merge with itself */
00800            mu = bp[m] ;
00801            if( BDEP(mu)-mz <= hpf ){       /* basin not TOO much deeper */
00802              MERGE_BASIN(mq,mu) ;
00803            }
00804          }
00805        }
00806      }
00807    } /* end of loop over voxels */
00808 
00809    /* at this point, all voxels in svox are assigned to a basin */
00810 
00811    free((void *)isvox) ;
00812 
00813    /* count number of basines left */
00814 
00815    for( mu=m=0 ; m < nbtop ; m++ )
00816      if( baslist[m] != NULL ) mu++ ;
00817 
00818    if( verb ) fprintf(stderr,"\n++ %d active basins left, out of %d\n",mu,nbtop) ;
00819 
00820    bcount = (int *) calloc(sizeof(int),mu) ;     /* number in each basin */
00821    bname  = (int *) calloc(sizeof(int),mu) ;
00822    isvox  = (int *) calloc(sizeof(int),nbtop) ;  /* new index */
00823 
00824    for( m=ii=0 ; m < nbtop ; m++ )
00825      if( baslist[m] != NULL ){ isvox[m] = ii; bname[ii] = ii; ii++; KILL_BASIN(m); }
00826    free((void *)baslist) ;
00827 
00828    for( pp=0 ; pp < nvox ; pp++ ){
00829      m  = svox[pp].basin ;           /* old basin name for this voxel */
00830      ii = isvox[m] ;                 /* new basin name for this voxel */
00831      svox[pp].basin = ii ;           /* reassign name in this voxel */
00832      bcount[ii]++ ;                  /* count number in this basin */
00833    }
00834 
00835    tim = mri_new_conforming( sim , MRI_short ) ;  /* output image */
00836    MRI_COPY_AUX(tim,sim) ;
00837    tar = MRI_SHORT_PTR(tim) ;
00838 
00839    for( ii=0 ; ii < mu ; ii++ ) bcount[ii] = -bcount[ii] ;
00840    qsort_intint( mu , bcount , bname ) ;  /* sort into decreasing order */
00841    for( ii=0 ; ii < mu ; ii++ ) bcount[ii] = -bcount[ii] ;
00842 
00843    if( verb )
00844      fprintf(stderr," + top 9 basin counts: %d %d %d %d %d %d %d %d %d\n",
00845              bcount[0] , bcount[1] , bcount[2] , bcount[3] ,
00846              bcount[4] , bcount[5] , bcount[6] , bcount[7] , bcount[8] ) ;
00847 
00848    for( ii=0 ; ii < mu ; ii++ ) isvox[ii] = ii ;
00849    qsort_intint( mu , bname , isvox ) ;
00850 
00851    for( pp=0 ; pp < nvox ; pp++ ){
00852      m = svox[pp].basin ; jj = isvox[m]+1 ; if( jj > 32767 ) jj = 32767 ;
00853      tar[IJK(svox[pp].i,svox[pp].j,svox[pp].k)] = jj ;
00854    }
00855 
00856    free((void *)isvox) ; free((void *)svox );
00857    free((void *)bcount); free((void *)bname);
00858 
00859    return tim ;
00860 }
00861 
00862 /*------------------------------------------------------------------------*/
00863 
00864 static byte * make_peel_mask( int nx, int ny, int nz , byte *mmm, int pdepth )
00865 {
00866    int nxy=nx*ny,nxyz=nxy*nz , ii,jj,kk , ijk , bot,top , pd=pdepth ;
00867    byte *ppp ;
00868 
00869    if( mmm == NULL || pdepth <= 1 ) return NULL ;
00870 
00871    ppp = (byte *)calloc(sizeof(byte),nxyz) ;
00872 
00873    for( kk=0 ; kk < nz ; kk++ ){
00874     for( jj=0 ; jj < ny ; jj++ ){
00875       ijk = jj*nx + kk*nxy ;
00876       for( bot=0 ; bot < nx && !mmm[bot+ijk]; bot++ ) ;
00877       top = bot+pd ; if( top >= nx ) continue ;
00878       for( ii=bot+1 ; ii <= top && mmm[ii+ijk] ; ii++ ) ;
00879       if( ii <= top ){
00880         top = ii; for( ii=bot ; ii <= top ; ii++ ) ppp[ii+ijk] = 1;
00881       }
00882    }}
00883 
00884    for( kk=0 ; kk < nz ; kk++ ){
00885     for( jj=0 ; jj < ny ; jj++ ){
00886       ijk = jj*nx + kk*nxy ;
00887       for( top=nx-1 ; top >= 0 && !mmm[top+ijk]; top-- ) ;
00888       bot = top-pd ; if( bot < 0 ) continue ;
00889       for( ii=top-1 ; ii >= bot && mmm[ii+ijk] ; ii-- ) ;
00890       if( ii >= bot ){
00891         bot = ii; for( ii=top ; ii >= bot ; ii-- ) ppp[ii+ijk] = 1;
00892       }
00893    }}
00894 
00895    for( kk=0 ; kk < nz ; kk++ ){
00896     for( ii=0 ; ii < nx ; ii++ ){
00897       ijk = ii + kk*nxy ;
00898       for( bot=0 ; bot < ny && !mmm[bot*nx+ijk] ; bot++ ) ;
00899       top = bot+pd ;
00900       if( top >= ny ) continue ;
00901       for( jj=bot+1 ; jj <= top && mmm[jj*nx+ijk] ; jj++ ) ;
00902       if( jj <= top ){
00903         top = jj; for( jj=bot ; jj <= top ; jj++ ) ppp[jj*nx+ijk] = 1;
00904       }
00905    }}
00906 
00907    for( kk=0 ; kk < nz ; kk++ ){
00908     for( ii=0 ; ii < nx ; ii++ ){
00909       ijk = ii + kk*nxy ;
00910       for( top=ny-1 ; top >= 0 && !mmm[top*nx+ijk] ; top-- ) ;
00911       bot = top-pd ; if( bot < 0 ) continue ;
00912       for( jj=top-1 ; jj >= bot && mmm[jj*nx+ijk] ; jj-- ) ;
00913       if( jj >= bot ){
00914         bot = jj; for( jj=top ; jj >= bot ; jj-- ) ppp[jj*nx+ijk] = 1;
00915       }
00916    }}
00917 
00918 #if 1
00919     for( jj=0 ; jj < ny ; jj++ ){
00920      for( ii=0 ; ii < nx ; ii++ ){
00921        ijk = ii + jj*nx ;
00922        for( top=nz-1 ; top >= 0 && !mmm[top*nxy+ijk] ; top-- ) ;
00923        bot = top-pd ; if( bot < 0 ) continue ;
00924        for( kk=top-1 ; kk >= bot && mmm[kk*nxy+ijk] ; kk-- ) ;
00925        if( kk >= bot ){
00926          bot = kk; for( kk=top ; kk >= bot ; kk-- ) ppp[kk*nxy+ijk] = 1;
00927        }
00928     }}
00929 #endif
00930 
00931    kk = mask_count(nxyz,ppp) ;
00932    if( kk == 0 ){ free((void *)ppp) ; return NULL ; }
00933    if( verb ) fprintf(stderr," + Initial peel mask has %d voxels\n",kk ) ;
00934    THD_mask_erode( nx,ny,nz, ppp ) ;
00935    THD_mask_clust( nx,ny,nz, ppp ) ;
00936    kk = mask_count(nxyz,ppp) ;
00937    if( kk == 0 ){ free((void *)ppp) ; return NULL ; }
00938    if( verb ) fprintf(stderr," + Final   peel mask has %d voxels\n",kk ) ;
00939 
00940    return ppp ;
00941 }
00942 
00943 /*------------------------------------------------------------------------*/
00944 
00945 static void zedit_mask( int nx, int ny, int nz, byte *mmm, int zdepth, int zbot )
00946 {
00947    int nxy=nx*ny,nxyz=nxy*nz , ii,jj,kk , ijk , bot,top ;
00948    int zd=zdepth , zb=zbot , zt , zslab , zz ;
00949    byte *ppp , *zzz ;
00950 
00951    if( mmm == NULL ) return ;
00952 
00953    if( zd < 1  ) zd = 1  ;
00954    if( zb < zd ) zb = zd ;
00955    zslab = 2*zd+1 ;
00956 
00957    for( kk=nz-1 ; kk >= zb ; kk-- ){
00958      jj = mask_count( nxy , mmm+kk*nxy ) ;
00959      if( jj > 0.005*nxy ) break ;
00960    }
00961    zt = kk-zd ; if( zt < zb ) return ;
00962 
00963    ppp = (byte *)calloc(sizeof(byte),nxyz) ;
00964    zzz = (byte *)calloc(sizeof(byte),nxy*zslab) ;
00965 
00966    for( zz=zb ; zz <= zt ; zz++ ){
00967      memcpy( zzz , mmm+(zz-zd)*nxy , nxy*zslab ) ;
00968      THD_mask_erode( nx,ny,zslab, zzz ) ;
00969      THD_mask_clust( nx,ny,zslab, zzz ) ;
00970      memcpy( ppp+zz*nxy , zzz+zd*nxy , nxy ) ;
00971    }
00972    free((void *)zzz) ;
00973    memcpy( mmm+zb*nxy , ppp+zb*nxy , (zt-zb+1)*nxy ) ;
00974    free((void *)ppp) ;
00975    THD_mask_erode( nx,ny,nz, mmm ) ;
00976    THD_mask_clust( nx,ny,nz, mmm ) ;
00977 }
00978 
00979 /*======================================================================*/
00980    /** ZSS: Definitions moved to thd_brainormalize.h **/
00981 /*----------------------------------------------------------------------*/
00982 /*! Index warping function for mri_warp3D() call. */
00983 
00984 static float ai,bi , aj,bj , ak,bk ;  /* set below */
00985 
00986 static void ijkwarp( float  i, float  j, float  k ,
00987                      float *x, float *y, float *z  )
00988 {
00989   *x = ai*i + bi ;
00990   *y = aj*j + bj ;
00991   *z = ak*k + bk ;
00992 }
00993 
00994 static void ijk_invwarp(   float x, float y, float z ,  
00995                            float  *i, float  *j, float  *k)
00996 {
00997   *i = ( x - bi ) / ai;
00998   *j = ( y - bj ) / aj;
00999   *k = ( z - bk ) / ak;
01000 }
01001 
01002 
01003 /*! 
01004    \brief takes in voxel indices into the Spat Normed volume (RAI) and
01005    returns voxel indices and coordinates in the original volume. Used
01006    to figure out shift to apply to surface model to align it with original volume.
01007    
01008    \param ispat (float) 3D i index into spat norm volume
01009    \param jspat (float) 3D j index into spat norm volume
01010    \param kspat (float) 3D k index into spat norm volume
01011    \param iorig (float*) 3D i index into original volume
01012    \param jorig (float*) 3D j index into original volume
01013    \param korig (float*) 3D k index into original volume
01014    \param origset (THD_3dim_dataset *) Le original dataset
01015    \param *xrai_orig (float *) X coordinate in original volume (dicomm)
01016    \param *yrai_orig (float *) Y coordinate in original volume (dicomm)
01017    \param *zrai_orig (float *) Z coordinate in original volume (dicomm)
01018    
01019    ZSS Sometime in April 05
01020    
01021    \sa brainnormalize_inv_coord
01022 */
01023 void brainnormalize_coord( float  ispat, float  jspat, float  kspat ,
01024                            float *iorig, float *jorig, float *korig ,
01025                            THD_3dim_dataset *origset,
01026                            float *xrai_orig, float *yrai_orig, float *zrai_orig)
01027 {
01028    THD_dataxes * daxes ;
01029    THD_fvec3     fv, fvdic ;
01030    float irai, jrai, krai;
01031    
01032    /* find out corresponding indices in original dset */
01033    ijkwarp(ispat, jspat, kspat , &irai, &jrai, &krai);
01034    
01035 
01036     
01037    /* These indices assume an RAI dset orientation */
01038       /* Find out what these indices should be in origset's orientation */
01039       switch( origset->daxes->xxorient ){
01040         case ORI_R2L_TYPE: *iorig =  irai ; break ;
01041         case ORI_L2R_TYPE: *iorig =  origset->daxes->nxx - irai ; break ;
01042         case ORI_P2A_TYPE: *iorig =  origset->daxes->nxx - jrai  ; break ;
01043         case ORI_A2P_TYPE: *iorig =  jrai ; break ;
01044         case ORI_I2S_TYPE: *iorig =  krai ; break ;
01045         case ORI_S2I_TYPE: *iorig =  origset->daxes->nxx - krai ; break ;
01046       }
01047       switch( origset->daxes->yyorient ){
01048         case ORI_R2L_TYPE: *jorig =  irai ; break ;
01049         case ORI_L2R_TYPE: *jorig =  origset->daxes->nyy - irai  ; break ;
01050         case ORI_P2A_TYPE: *jorig =  origset->daxes->nyy - jrai ; break ;
01051         case ORI_A2P_TYPE: *jorig =  jrai ; break ;
01052         case ORI_I2S_TYPE: *jorig =  krai ; break ;
01053         case ORI_S2I_TYPE: *jorig =  origset->daxes->nyy - krai ; break ;
01054       }
01055       switch( origset->daxes->zzorient ){
01056         case ORI_R2L_TYPE: *korig =  irai ; break ;
01057         case ORI_L2R_TYPE: *korig =  origset->daxes->nzz - irai ; break ;
01058         case ORI_P2A_TYPE: *korig =  origset->daxes->nzz - jrai ; break ;
01059         case ORI_A2P_TYPE: *korig =  jrai ; break ;
01060         case ORI_I2S_TYPE: *korig =  krai ; break ;
01061         case ORI_S2I_TYPE: *korig =  origset->daxes->nzz - krai ; break ;
01062       }
01063       
01064             
01065       /* change indices into mm coords in orig dset*/
01066       daxes = CURRENT_DAXES(origset) ;
01067 
01068       fv.xyz[0] = daxes->xxorg + *iorig * daxes->xxdel ;  /* 3dfind_to_3dmm */
01069       fv.xyz[1] = daxes->yyorg + *jorig * daxes->yydel ;
01070       fv.xyz[2] = daxes->zzorg + *korig * daxes->zzdel ;
01071 
01072       fvdic = THD_3dmm_to_dicomm(origset,   fv );                /* 3dmm_to_dicomm  */
01073       *xrai_orig = fvdic.xyz[0];
01074       *yrai_orig = fvdic.xyz[1]; 
01075       *zrai_orig = fvdic.xyz[2];
01076        
01077    /* report for sanity */
01078    if (0) {
01079       fprintf(stderr,   "brainnormalize_coord:\n"
01080                      " ijk_spat_rai = [%f %f %f]\n"
01081                      " ijk_orig_rai = [%f %f %f] (in rai order, not native to iset!)\n"
01082                      " ijk_orig     = [%f %f %f] (in native order)\n"
01083                      " XYZ_orig     = [%f %f %f]\n"
01084                      " Origin spat = [%f %f %f]\n", 
01085                      ispat, jspat, kspat,
01086                      irai, jrai, krai ,
01087                      *iorig, *jorig, *korig ,
01088                      *xrai_orig, *yrai_orig, *zrai_orig,
01089                      THD_BN_XORG, THD_BN_YORG, THD_BN_ZORG);   
01090    }      
01091    return;
01092 } 
01093 
01094 /*----------------------------------------------------------------------
01095    (a) shortize input and flip brick so that orientation is RAI
01096    (b) find clip levels and create a binary mask
01097    (c) find S-most slice that has at least 10% above clip level;
01098        zero out mask above that slice and also more than 160 mm below
01099    (d) apply mask to image volume
01100    (e) resample to master dataset grid, with CM at (0,20,0)
01101    \param imout_origp (MRI_IMAGE **) a contraption by ZSS to return a SpatNormed
01102    version that will eventually be rewritten in im's orientation. Pass NULL if
01103    you want nothing to do with it.
01104 ------------------------------------------------------------------------*/
01105 
01106 MRI_IMAGE * mri_brainormalize( MRI_IMAGE *im, int xxor, int yyor, int zzor, MRI_IMAGE **imout_origp , MRI_IMAGE **imout_edge)
01107 {
01108    MRI_IMAGE *sim , *tim , *bim ;
01109    short *sar , sval ;
01110    int ii,jj,kk,ijk,ktop,kbot , nx,ny,nz,nxy,nxyz ;
01111    float val , icm,jcm,kcm,sum , dx,dy,dz ;
01112    byte *mask , *bar ;
01113    float sim_dx, sim_dy, sim_dz, sim_xo, sim_yo, sim_zo;
01114    int sim_nx, sim_ny, sim_nz;
01115    int *zcount , *hist,*gist , z1,z2,z3 ;
01116    MRI_IMAGE *imout_orig = NULL;
01117    
01118 ENTRY("mri_brainormalize") ;
01119 
01120    if( im == NULL || xxor < 0 || xxor > LAST_ORIENT_TYPE ||
01121                      yyor < 0 || yyor > LAST_ORIENT_TYPE ||
01122                      zzor < 0 || zzor > LAST_ORIENT_TYPE   ) RETURN(NULL) ;
01123 
01124    if( im->nx < 16 || im->ny < 16 || im->nz < 16 ) RETURN(NULL) ;
01125 
01126    val = mri_maxabs(im) ; if( val <= 0.0 ) RETURN(NULL) ;
01127 
01128    /* make a short copy */
01129 
01130    if( verb ) fprintf(stderr,"++mri_brainormalize: copying input\n") ;
01131 
01132    if( im->kind == MRI_short || im->kind == MRI_byte || im->kind == MRI_rgb )
01133      tim = mri_to_short( 1.0 , im ) ;
01134    else
01135      tim = mri_to_short( 32767.0/val , im ) ;
01136 
01137    /* flip to RAI orientation */
01138 
01139    ii = jj = kk = 0 ;
01140    switch( xxor ){
01141      case ORI_R2L_TYPE: ii =  1 ; break ;
01142      case ORI_L2R_TYPE: ii = -1 ; break ;
01143      case ORI_P2A_TYPE: jj = -1 ; break ;
01144      case ORI_A2P_TYPE: jj =  1 ; break ;
01145      case ORI_I2S_TYPE: kk =  1 ; break ;
01146      case ORI_S2I_TYPE: kk = -1 ; break ;
01147    }
01148    switch( yyor ){
01149      case ORI_R2L_TYPE: ii =  2 ; break ;
01150      case ORI_L2R_TYPE: ii = -2 ; break ;
01151      case ORI_P2A_TYPE: jj = -2 ; break ;
01152      case ORI_A2P_TYPE: jj =  2 ; break ;
01153      case ORI_I2S_TYPE: kk =  2 ; break ;
01154      case ORI_S2I_TYPE: kk = -2 ; break ;
01155    }
01156    switch( zzor ){
01157      case ORI_R2L_TYPE: ii =  3 ; break ;
01158      case ORI_L2R_TYPE: ii = -3 ; break ;
01159      case ORI_P2A_TYPE: jj = -3 ; break ;
01160      case ORI_A2P_TYPE: jj =  3 ; break ;
01161      case ORI_I2S_TYPE: kk =  3 ; break ;
01162      case ORI_S2I_TYPE: kk = -3 ; break ;
01163    }
01164 
01165    if( ii==1 && jj==2 && kk==3 ){      /* no flip needed */
01166      sim = tim ;
01167    } else {                            /* flipization */
01168      if( verb )
01169        fprintf(stderr,"++mri_brainormalize: flipping to RAI orientation\n") ;
01170      sim = mri_flip3D( ii,jj,kk , tim ) ;
01171      mri_free(tim) ;
01172      if( sim == NULL ) RETURN(NULL) ;  /* bad orientation codes? */
01173    }
01174 
01175    sar = MRI_SHORT_PTR(sim) ;
01176    if( sar == NULL ){ mri_free(sim); RETURN(NULL); }  /* bad image? */
01177 
01178    /* make a binary mask */
01179 
01180    nx = sim->nx ; ny = sim->ny ; nz = sim->nz ; nxy = nx*ny ; nxyz = nxy*nz ;
01181    dx = fabs(sim->dx) ; if( dx == 0.0 ) dx = 1.0 ;
01182    dy = fabs(sim->dy) ; if( dy == 0.0 ) dy = 1.0 ;
01183    dz = fabs(sim->dz) ; if( dz == 0.0 ) dz = 1.0 ;
01184       
01185    /* save some info to create an output image with the same number of slices as original image*/
01186    sim_dx = sim->dx; sim_dy = sim->dy; sim_dz = sim->dz;
01187    sim_xo = 0.0; sim_yo = 0.0; sim_zo = 0.0;  /* origins are added after this function returns.*/
01188    sim_nx = sim->nx; sim_ny = sim->ny; sim_nz = sim->nz; 
01189    
01190    if( verb ) fprintf(stderr,"++mri_brainormalize: making mask\n") ;
01191    mask = mri_short2mask( sim ) ;
01192 
01193    if( mask == NULL ){ mri_free(sim); RETURN(NULL); }
01194 
01195    /* fill in any isolated holes in mask */
01196 
01197    (void) THD_mask_fillin_once( nx,ny,nz , mask , 2 ) ;  /* thd_automask.c */
01198           THD_mask_dilate     ( nx,ny,nz , mask , 5 ) ;
01199           THD_mask_dilate     ( nx,ny,nz , mask , 5 ) ;
01200 
01201    kk = mask_count(nxyz,mask) ;
01202    if( verb )
01203      fprintf(stderr,"++mri_brainormalize: filled in mask now has %d voxels\n",kk) ;
01204 
01205    if( kk <= 999 ){ free((void *)mask); mri_free(sim); RETURN(NULL); }
01206 
01207    /* descending from Superior:
01208         count biggest blob in each slice
01209         find Superiormost location that has 3
01210           slices in a row with "a lot" of stuff
01211         zero out all stuff out above that slice */
01212 
01213    zcount = (int *)  malloc( sizeof(int) *nz  ) ;  /* slice counts */
01214    for( kk=nz-1 ; kk >= 0 ; kk-- ){
01215      zcount[kk] = mask_count( nxy , mask+kk*nxy ) ;
01216    }
01217 
01218 #if 0
01219    if( verb ){
01220      fprintf(stderr,"++mri_brainormalize: zcount from top slice #%d\n",nz-1) ;
01221      for( kk=nz-1 ; kk >= 0 ; kk-- ){
01222        fprintf(stderr," %.3f",((double)(zcount[kk]))/((double)nxy) ) ;
01223        if( (nz-kk)%10 == 0 && kk > 0 ) fprintf(stderr,"\n") ;
01224      }
01225      fprintf(stderr,"\n") ;
01226    }
01227 #endif
01228 
01229    /* search down for topmost slice that meets the criterion */
01230 
01231    z1 = (int)(0.010*nxy) ;
01232    z2 = (int)(0.015*nxy) ;
01233    z3 = (int)(0.020*nxy) ;
01234    for( kk=nz-1 ; kk > 2 ; kk-- )
01235      if( zcount[kk] >= z1 && zcount[kk-1] >= z2 && zcount[kk-2] >= z3 ) break ;
01236 
01237    free((void *)zcount) ;
01238    if( kk <= 2 ){ free((void *)mask); mri_free(sim); RETURN(NULL); }
01239 
01240    /* zero out all above the slice we just found */
01241 
01242    ktop = kk ;
01243    if( ktop < nz-1 ){
01244      if( verb )
01245        fprintf(stderr,"++mri_brainormalize: top clip above slice %d\n",ktop) ;
01246      memset( mask+(ktop+1)*nxy , 0 , nxy*(nz-1-ktop)*sizeof(byte) ) ;
01247    }
01248 
01249    /* find slice index THD_BN_ZHEIGHT mm below that top slice */
01250 
01251    jj = (int)( ktop-THD_BN_ZHEIGHT/dz ) ;
01252    if( jj >= 0 ){
01253      if( verb )
01254        fprintf(stderr,"++mri_brainormalize: bot clip below slice %d\n",jj) ;
01255      memset( mask , 0 , nxy*(jj+1)*sizeof(byte) ) ;
01256    }
01257 
01258    kk = mask_count(nxyz,mask) ;
01259    if( kk <= 999 ){ free((void *)mask); mri_free(sim); RETURN(NULL); }
01260 
01261    /* apply mask to image (will also remove any negative values) */
01262 
01263    if( verb )
01264      fprintf(stderr,"++mri_brainormalize: applying mask to image; %d voxels\n",kk) ;
01265    for( ii=0 ; ii < nxyz ; ii++ ) if( !mask[ii] ) sar[ii] = 0 ;
01266 
01267    free((void *)mask) ;  /* done with this mask */
01268 
01269    /* compute CM of masked image (indexes, not mm) */
01270 
01271    icm = jcm = kcm = sum = 0.0 ;
01272 #ifndef THD_BN_CMTOP
01273    kbot = 0 ;
01274    ktop = nz-1 ;
01275 #else
01276    kbot = (int)rint( ktop-110.0/dz ); if( kbot < 0 ) kbot = 0;
01277 #endif
01278    for( ijk=kbot*nxy,kk=kbot ; kk <= ktop ; kk++ ){
01279     for( jj=0 ; jj < ny ; jj++ ){
01280      for( ii=0 ; ii < nx ; ii++,ijk++ ){
01281        val = (float)sar[ijk] ;
01282        sum += val ;
01283        icm += val * ii ;
01284        jcm += val * jj ;
01285        kcm += val * kk ;
01286    }}}
01287    if( sum == 0.0 ){ mri_free(sim); RETURN(NULL); }  /* huh? */
01288 
01289    ai = thd_bn_dxyz/dx ; bi = icm/sum - ai*(THD_BN_XCM-THD_BN_XORG)/thd_bn_dxyz ;
01290    aj = thd_bn_dxyz/dy ; bj = jcm/sum - aj*(THD_BN_YCM-THD_BN_YORG)/thd_bn_dxyz ;
01291    ak = thd_bn_dxyz/dz ; bk = kcm/sum - ak*(THD_BN_ZCM-THD_BN_ZORG)/thd_bn_dxyz ;
01292 
01293    if( verb ) fprintf(stderr,"++mri_brainormalize: warping to standard grid\n a = [%f %f %f], b = [%f %f %f]\n", ai, aj, ak, bi, bj, bk) ;
01294 
01295    mri_warp3D_method( MRI_CUBIC ) ;
01296    tim = mri_warp3D( sim , thd_bn_nx,thd_bn_ny,thd_bn_nz , ijkwarp ) ;
01297    mri_free(sim) ;
01298 
01299    tim->dx = tim->dy = tim->dz = thd_bn_dxyz ;
01300    tim->xo = THD_BN_XORG ;
01301    tim->yo = THD_BN_YORG ;
01302    tim->zo = THD_BN_ZORG ;
01303    
01304    nx = tim->nx ; ny = tim->ny ; nz = tim->nz ; nxy = nx*ny ; nxyz = nxy*nz ;
01305    sar = MRI_SHORT_PTR(tim) ;
01306 
01307    /*-- rescale to partially uniformize --*/
01308 
01309    { clipvec bvec ; float bval , sv ;
01310      bvec = get_octant_clips( tim , 0.40f ) ;
01311      if( bvec.clip_000 > 0.0f ){
01312        for( ijk=kk=0 ; kk < nz ; kk++ ){
01313         for( jj=0 ; jj < ny ; jj++ ){
01314          for( ii=0 ; ii < nx ; ii++,ijk++ ){
01315            bval = pointclip( ii,jj,kk , &bvec ) ; /* cliplevel here */
01316            sv   = 1000.0f * sar[ijk] / bval ;
01317            sar[ijk] = SHORTIZE(sv) ;
01318          }
01319      }}}
01320    }
01321 
01322    /*-- build another mask now --*/
01323    if( !AFNI_noenv("REMASK") ){
01324      int sbot,stop , nwid , cbot,ctop , ibot,itop ;
01325      float dsum , ws , *wt ;
01326      int pval[128] , wval[128] , npk , tval , nmask,nhalf ;
01327      float pdif ;
01328      short mbot,mtop ;
01329 
01330      /* build histogram */
01331 
01332      hist = (int *) calloc(sizeof(int),32768) ;
01333      gist = (int *) calloc(sizeof(int),32768) ;
01334 
01335      memset( hist , 0 , sizeof(int)*32768 ) ;
01336      for( ii=0 ; ii < nxyz ; ii++ ) hist[sar[ii]]++ ;
01337      for( sbot=1 ; sbot < 32768 && hist[sbot]==0 ; sbot++ ) ; /* nada */
01338      if( sbot == 32768 ) goto Remask_Done ;
01339      for( stop=32768-1 ; stop > sbot && hist[stop]==0 ; stop-- ) ; /* nada */
01340      if( stop == sbot ) goto Remask_Done ;
01341 
01342      /* find median */
01343      
01344      nmask = 0 ;
01345      for( ii=sbot ; ii <= stop ; ii++ ) nmask += hist[ii] ;
01346      nhalf = nmask / 2 ; nmask = 0 ;
01347      for( ii=sbot ; ii <= stop && nmask < nhalf ; ii++ ) nmask += hist[ii] ;
01348      cbot = 0.40 * ii ;
01349      ctop = 1.60 * ii ;
01350 
01351 #if 0
01352      /* smooth histogram */
01353 
01354      nwid = rint(0.10*cbot) ;
01355 
01356      if( nwid <= 0 ){
01357        memcpy( gist , hist , sizeof(int)*32768 ) ;
01358      } else {
01359        ws = 0.0f ;
01360        wt = (float *)malloc(sizeof(float)*(2*nwid+1)) ;
01361        for( ii=0 ; ii <= 2*nwid ; ii++ ){
01362          wt[ii] = nwid-abs(nwid-ii) + 0.5f ;
01363          ws += wt[ii] ;
01364        }
01365        for( ii=0 ; ii <= 2*nwid ; ii++ ) wt[ii] /= ws ;
01366        for( jj=cbot ; jj <= ctop ; jj++ ){
01367          ibot = jj-nwid ; if( ibot < sbot ) ibot = sbot ;
01368          itop = jj+nwid ; if( itop > stop ) itop = stop ;
01369          ws = 0.0 ;
01370          for( ii=ibot ; ii <= itop ; ii++ )
01371            ws += wt[nwid-jj+ii] * hist[ii] ;
01372          gist[jj] = rint(ws) ;
01373        }
01374        free(wt) ;
01375      }
01376 
01377      /* scan for peaks */
01378 
01379      npk = 0 ;
01380      for( ii=cbot+2 ; ii <= ctop-2 ; ii++ ){
01381        if( gist[ii] > gist[ii-1] &&
01382            gist[ii] > gist[ii-2] &&
01383            gist[ii] > gist[ii+1] &&
01384            gist[ii] > gist[ii+2]   ){
01385              pval[npk]=ii; wval[npk++] = gist[ii];
01386            }
01387 
01388        else if( gist[ii] == gist[ii+1] &&   /* very special case */
01389                 gist[ii] >  gist[ii-1] &&
01390                 gist[ii] >  gist[ii-2] &&
01391                 gist[ii] >  gist[ii+2]   ){
01392                   pval[npk]=ii+0.5; wval[npk++] = gist[ii];
01393                 }
01394 
01395        else if( gist[ii] == gist[ii+1] &&   /* super special case */
01396                 gist[ii] == gist[ii-1] &&
01397                 gist[ii] >  gist[ii-2] &&
01398                 gist[ii] >  gist[ii+2]   ){
01399                   pval[npk]=ii; wval[npk++] = gist[ii];
01400                 }
01401      }
01402 
01403      if( npk > 2 ){  /* find largest two peaks and keep only them */
01404        float pval_top, pval_2nd, wval_top, wval_2nd , www; int iii,itop ;
01405        www = wval[0] ; iii = 0 ;
01406        for( ii=1 ; ii < npk ; ii++ ){
01407          if( wval[ii] > www ){ www = wval[ii] ; iii = ii ; }
01408        }
01409        pval_top = pval[iii] ; wval_top = www ; itop = iii ;
01410        www = -1 ; iii = -1 ;
01411        for( ii=0 ; ii < npk ; ii++ ){
01412          if( ii != itop && wval[ii] > www ){ www = wval[ii] ; iii = ii ; }
01413        }
01414        pval_2nd = pval[iii] ; wval_2nd = www ;
01415 
01416        /* make sure peaks are increasing in pval */
01417 
01418        if( pval_top < pval_2nd ){
01419          pval[0] = pval_top ; wval[0] = wval_top ;
01420          pval[1] = pval_2nd ; wval[1] = wval_2nd ;
01421        } else {
01422          pval[0] = pval_2nd ; wval[0] = wval_2nd ;
01423          pval[1] = pval_top ; wval[1] = wval_top ;
01424        }
01425        npk = 2 ;
01426      }
01427 
01428      if( npk == 2 ){
01429        jj = gist[pval[0]] ; tval = pval[0] ;
01430        for( ii=pval[0]+1 ; ii < pval[1] ; ii++ ){
01431          if( gist[ii] < jj ){ tval = ii ; jj = gist[ii] ; }
01432        }
01433 
01434        pdif = 1.5f * (tval-pval[0]) ;
01435        if( pdif < pval[1]-pval[0] ) pdif = pval[1]-pval[0] ;
01436        mbot = (short)(pval[0]-pdif) ;
01437        if( mbot < cbot ) mbot = cbot ;
01438 
01439        pdif = 1.5f * (pval[1]-tval) ;
01440        if( pdif < pval[1]-pval[0] ) pdif = pval[1]-pval[0] ;
01441        mtop = (short)(pval[1]+pdif) ;
01442        if( mtop > ctop ) mtop = ctop ;
01443      } else {
01444        mbot = cbot ; mtop = ctop ;
01445      }
01446      mtop = stop+1 ;  /* effectively, no threshold here */
01447 #endif
01448 
01449      mbot = cbot ; mtop = 32767 ;
01450 
01451      if( verb )
01452       fprintf(stderr,"++mri_brainormalize: masking standard image %d..%d\n",mbot,mtop) ;
01453 
01454      mask = (byte *) malloc( sizeof(byte)*nxyz ) ;
01455      for( ii=0 ; ii < nxyz ; ii++ )
01456        mask[ii] = (sar[ii] > mbot) && (sar[ii] < mtop) ;
01457 
01458      THD_mask_erode( nx,ny,nz, mask ) ;
01459      THD_mask_clust( nx,ny,nz, mask ) ;
01460      for( ii=0 ; ii < nxyz ; ii++ ) mask[ii] = !mask[ii] ;
01461      THD_mask_clust( nx,ny,nz, mask ) ;
01462      for( ii=0 ; ii < nxyz ; ii++ ) mask[ii] = !mask[ii] ;
01463 
01464      for( ii=0 ; ii < nxyz ; ii++ ) if( !mask[ii] ) sar[ii] = 0 ;
01465      free((void *)mask) ;
01466 
01467    Remask_Done:
01468      free((void *)hist) ; free((void *)gist) ;
01469 
01470    }
01471 #if 0
01472    else
01473 #endif
01474    {
01475      /*-- clip top 1% of values that have survived --*/
01476 
01477      hist = (int *) calloc(sizeof(int),32768) ;
01478      for( ii=0 ; ii < nxyz ; ii++ ) hist[sar[ii]]++ ;
01479      for( ii=kk=0 ; ii < 32767 ; ii++ ) kk += hist[ii] ;
01480      kk = (int)(0.01*kk) ; ktop = 0 ;
01481      for( jj=0,ii=32767 ; ii > 0 && jj < kk ; ii-- ){
01482        jj += hist[ii] ; if( hist[ii] > 0 && ktop == 0 ) ktop = ii ;
01483      }
01484      jj = ii ;
01485      if( verb ) fprintf(stderr," + 99%% clipping at %d (from %d)\n",jj,ktop) ;
01486      for( ii=0 ; ii < nxyz ; ii++ ) if( sar[ii] > jj ) sar[ii] = jj ;
01487 
01488      free((void *)hist) ;
01489    }
01490 
01491    /* distize? */
01492 
01493    if( AFNI_yesenv("DISTIZE") ){
01494      byte *ccc = (byte *)calloc(sizeof(byte),nxyz);
01495      short *ddd ;
01496      int kbot=(int)rint(0.45*nz) , ktop=(int)rint(0.65*nz) ,
01497          jbot=(int)rint(0.30*ny) , jtop=(int)rint(0.70*ny) ,
01498          ibot=(int)rint(0.30*nx) , itop=(int)rint(0.70*nx)  ;
01499 
01500      mask = (byte *)malloc( sizeof(byte)*nxyz ) ;
01501      for( ii=0 ; ii < nxyz ; ii++ ) mask[ii] = (sar[ii] > 0) ;
01502      for( kk=kbot ; kk <= ktop ; kk++ ){
01503       for( jj=jbot ; jj <= jtop ; jj++ ){
01504        for( ii=ibot ; ii <= itop ; ii++ ){
01505          ijk = ii + jj*nx + kk*nxy ;
01506          ccc[ijk] = mask[ijk] ;
01507      }}}
01508      if( verb ) fprintf(stderr," + distizing\n") ;
01509      ddd = THD_mask_distize( nx,ny,nz , mask , ccc ) ;
01510      if( ddd != NULL ){
01511        int id,jd,kd , ijk , dijk ; float ff ;
01512        for( ijk=0 ; ijk < nxyz ; ijk++ ){
01513          if( ddd[ijk] > 0 ){
01514            ii = ijk % nx ; jj = (ijk%nxy)/nx ; kk = ijk / nxy ;
01515                 if( ii < ibot ) id = ibot-ii ;
01516            else if( ii > itop ) id = ii-itop ; else id = 0 ;
01517                 if( jj < jbot ) jd = jbot-jj ;
01518            else if( jj > jtop ) jd = jj-jtop ; else jd = 0 ;
01519                 if( kk < kbot ) kd = kbot-kk ;
01520            else if( kk > ktop ) kd = kk-ktop ; else kd = 0 ;
01521            dijk = id+jd+kd+1 ;
01522            ff = (100.0f * ddd[ijk]) / (float)dijk - 98.9f ;
01523            if( ff > 255.0f ) ff = 255.0f ;
01524            sar[ijk] = (short)ff ;
01525          } else {
01526            sar[ijk] = 0 ;
01527          }
01528        }
01529        free((void *)ddd) ;
01530      }
01531      free((void *)mask); free((void *)ccc);
01532    }
01533 
01534    /* create a spat norm of the original volume ZSS */
01535    if (imout_origp) {
01536       mri_warp3D_method( MRI_NN ) ;
01537       if (1 && verb) fprintf(stderr,"thd_brainormalize (ZSS):\n n: %d %d %d\n d: %f %f %f\n o: %f %f %f\n ", sim_nx, sim_ny, sim_nz, sim_dx, sim_dy, sim_dz, sim_xo, sim_yo, sim_zo);
01538       imout_orig = mri_warp3D( tim, sim_nx, sim_ny, sim_nz, ijk_invwarp );
01539       imout_orig->dx = sim_dx; imout_orig->dy = sim_dy; imout_orig->dz = sim_dz; 
01540       imout_orig->xo = sim_xo; imout_orig->yo = sim_yo; imout_orig->zo = sim_zo; 
01541       *imout_origp = imout_orig;
01542    }
01543    
01544    if (imout_edge) { /* create an edge version NOT EXISTING YET */
01545       *imout_edge = NULL;
01546    }
01547       
01548    /*-- convert output to bytes --*/
01549 
01550    bim = mri_new_conforming( tim , MRI_byte ) ;
01551    MRI_COPY_AUX(bim,tim) ;
01552    bar = MRI_BYTE_PTR(bim) ;
01553 
01554    jj = 0 ;
01555    for( ii=0 ; ii < nxyz ; ii++ ) if( sar[ii] > jj ) jj = sar[ii] ;
01556 
01557    if( jj > 255 ){
01558      float fac = 255.0 / jj ;
01559      if( verb ) fprintf(stderr," + scaling by fac=%g\n",fac) ;
01560      for( ii=0 ; ii < nxyz ; ii++ ) bar[ii] = (byte)(fac*sar[ii]+0.49) ;
01561    } else {
01562      for( ii=0 ; ii < nxyz ; ii++ ) bar[ii] = (byte)sar[ii] ;
01563    }
01564    mri_free(tim) ;
01565 
01566    /*-- done!!! --*/
01567 
01568    RETURN(bim) ;
01569 }