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 #define CFRAC      0.20
00010 #define LIM(x,b,t) ((x<b) ? b : (x>t) ? t : x)
00011 
00012 void qsort_sh( int n , short * a ) ;  /* at end of file */
00013 
00014 /*-------- 06 Jul 2000 - RWCox ----------------------------------------*/
00015 
00016 THD_3dim_dataset * WINsorize( THD_3dim_dataset *inset ,
00017                               int nrep , int cbot , int ctop ,
00018                               float irad , char *prefix ,
00019                               int keep_zero , int clipval , byte *mask )
00020 {
00021    THD_3dim_dataset *outset ;
00022    short *shin , *shout , *di,*dj,*dk , *tmp , val,nval ;
00023    MCW_cluster *cl ;
00024    int jj,kk , krep,kdiff, nx,ny,nz,nxy,nxyz , nd,dd ;
00025    int ip,jp,kp , nx1,ny1,nz1 , verb=1 ;
00026    int nrep_until ;
00027    register int ii,ijk ;
00028 
00029    /*- check inputs -*/
00030 
00031    if( inset == NULL || DSET_BRICK_TYPE(inset,0) != MRI_short ) return NULL ;
00032    DSET_load(inset) ;
00033    if( DSET_BRICK_ARRAY(inset,0) == NULL ) return NULL ;
00034 
00035    if( nrep == 0 ) return NULL ;
00036 
00037    if( nrep < 0 ){ nrep_until = abs(nrep) ; nrep = 999 ; }
00038    else          { nrep_until = 2 ; }
00039 
00040    if( irad < 0.0 ){ verb=0 ; irad = -irad ; }
00041 
00042    if( irad < 1.01 ) irad = 1.01 ;
00043    if( !THD_filename_ok(prefix) ) prefix = "Winsor" ;
00044 
00045    /*- build list of points to use -*/
00046 
00047    cl = MCW_build_mask( 0,0,0 , 1.0,1.0,1.0 , irad ) ;
00048 
00049    if( cl == NULL || cl->num_pt < 6 ){ KILL_CLUSTER(cl); return NULL; }
00050 
00051    ADDTO_CLUSTER(cl,0,0,0,0) ;
00052 
00053    di = cl->i ; dj = cl->j ; dk = cl->k ; nd = cl->num_pt ;
00054 
00055    if( verb ) fprintf(stderr,"+++ WINsorize irad=%f nbhd=%d\n",irad,nd) ;
00056 
00057    /*- make output array -*/
00058 
00059    nx = DSET_NX(inset) ; nx1 = nx-1 ;
00060    ny = DSET_NY(inset) ; ny1 = ny-1 ;
00061    nz = DSET_NZ(inset) ; nz1 = nz-1 ; nxy = nx*ny ; nxyz = nxy*nz ;
00062 
00063    shout = (short *) malloc(sizeof(short)*nxyz) ;
00064    tmp   = (short *) malloc(sizeof(short)*nd) ;
00065 
00066    if( nrep > 1 ){
00067       shin  = (short *) malloc(sizeof(short)*nxyz) ;
00068       memcpy( shin , DSET_BRICK_ARRAY(inset,0) , sizeof(short)*nxyz ) ;
00069    } else {
00070       shin = DSET_BRICK_ARRAY(inset,0) ;
00071    }
00072 
00073    if( cbot <= 0 || cbot >= nd-1 ){
00074       cbot = rint( CFRAC*nd ) ;
00075       if( cbot <= 0 ) cbot = 1 ;
00076       if( verb ) fprintf(stderr,"+++ WINsorize cbot=%d\n",cbot) ;
00077    }
00078    if( ctop <= cbot || cbot >= nd-1 ){
00079       ctop = nd-1-cbot ;
00080       if( verb ) fprintf(stderr,"+++ WINsorize ctop=%d\n",ctop) ;
00081    }
00082 
00083    /*- do the work -*/
00084 
00085    for( krep=0 ; krep < nrep ; krep++ ){
00086 
00087       kdiff = 0 ;  /* count of how many voxels were changed */
00088 
00089       for( kk=0 ; kk < nz ; kk++ ){        /* loops over 3D voxel indices */
00090          for( jj=0 ; jj < ny ; jj++ ){
00091             ijk = jj*nx+kk*nxy ;
00092             for( ii=0 ; ii < nx ; ii++,ijk++ ){
00093 
00094                if( mask != NULL && !mask[ijk] ){ shout[ijk]=shin[ijk]; continue; }
00095 
00096                val = shin[ijk] ;                        /* current voxel */
00097 
00098                if( clipval > 0 && val <= clipval )      /* 19 Oct 2001 */
00099                   val = shout[ijk] = 0 ;
00100 
00101                if( keep_zero && val == 0 ) continue ;   /* don't filter 0 */
00102 
00103                for( dd=0 ; dd < nd ; dd++ ){            /* loop over nbhd */
00104                   ip = ii+di[dd] ; ip = LIM(ip,0,nx1) ;
00105                   jp = jj+dj[dd] ; jp = LIM(jp,0,ny1) ;
00106                   kp = kk+dk[dd] ; kp = LIM(kp,0,nz1) ;
00107                   tmp[dd] = shin[ip+jp*nx+kp*nxy] ;
00108                }
00109 
00110                qsort_sh( nd , tmp ) ;
00111 
00112                shout[ijk] = nval = LIM(val,tmp[cbot],tmp[ctop]) ;
00113 
00114                if( nval != val ) kdiff++ ;
00115             }
00116          }
00117       }
00118 
00119       /* prepare for next iteration */
00120 
00121       if( verb ) fprintf(stderr,"+++ WINsorize iter%2d: # changed=%d\n",krep+1,kdiff) ;
00122 
00123       if( kdiff < nrep_until ) break ;
00124 
00125       if( krep < nrep-1 )
00126          memcpy( shin , shout , sizeof(short)*nxyz ) ;
00127    }
00128 
00129    /*- toss the trashola */
00130 
00131    KILL_CLUSTER(cl) ;
00132    free(tmp) ;
00133    if( shin != DSET_BRICK_ARRAY(inset,0) ) free(shin) ;
00134 
00135    /*- make output dataset */
00136 
00137    outset = EDIT_empty_copy( inset )  ;
00138 
00139    EDIT_dset_items( outset ,
00140                        ADN_prefix , prefix ,
00141                        ADN_nvals  , 1 ,
00142                        ADN_ntt    , 0 ,
00143                     ADN_none ) ;
00144 
00145    EDIT_substitute_brick( outset , 0 , MRI_short , shout ) ;
00146 
00147    return outset ;
00148 }
00149 
00150 /*****************************************************************************/
00151 
00152 void isort_sh( int n , short * ar )
00153 {
00154    register int  j , p ;        /* array indices */
00155    register short temp ;        /* a[j] holding place */
00156    register short * a = ar ;
00157 
00158    if( n < 2 ) return ;
00159 
00160    for( j=1 ; j < n ; j++ ){
00161 
00162      if( a[j] < a[j-1] ){       /* out of order */
00163         p    = j ;
00164         temp = a[j] ;
00165 
00166         do{
00167            a[p] = a[p-1] ;      /* now have a[p-1] > temp, so move it up */
00168            p-- ;
00169         } while( p > 0 && temp < a[p-1] ) ;
00170 
00171         a[p] = temp ;           /* finally, put temp in its place */
00172      }
00173    }
00174 }
00175 
00176 
00177 /*****************************************************************************/
00178 /* qsrec : recursive part of quicksort (stack implementation) */
00179 
00180 #define QS_STACK  1024  /* stack size */
00181 #define QS_SWAP(x,y) (temp=(x),(x)=(y),(y)=temp)
00182 
00183 static void qsrec_sh( int n , short * ar , int cutoff )
00184 {
00185    register int i , j ;           /* scanning indices */
00186    register short temp , pivot ;  /* holding places */
00187    register short * a = ar ;
00188 
00189    int left , right , mst , stack[QS_STACK] ;
00190 
00191    /* return if too short (insertion sort will clean up) */
00192 
00193    if( cutoff < 3 ) cutoff = 3 ;
00194    if( n < cutoff ) return ;
00195 
00196    /* initialize stack to start with whole array */
00197 
00198    stack[0] = 0   ;
00199    stack[1] = n-1 ;
00200    mst      = 2   ;
00201 
00202    /* loop while the stack is nonempty */
00203 
00204    while( mst > 0 ){
00205       right = stack[--mst] ;         /* subarray from left -> right */
00206       left  = stack[--mst] ;
00207 
00208       i = ( left + right ) / 2 ;     /* middle of subarray */
00209 
00210       /* sort the left, middle, and right a[]'s */
00211 
00212       if( a[left] > a[i]     ) QS_SWAP( a[left]  , a[i]     ) ;
00213       if( a[left] > a[right] ) QS_SWAP( a[left]  , a[right] ) ;
00214       if( a[i] > a[right]    ) QS_SWAP( a[right] , a[i]     ) ;
00215 
00216       pivot = a[i] ;                 /* a[i] is the median-of-3 pivot! */
00217       a[i]  = a[right] ;
00218 
00219       i = left ;                     /* initialize scanning */
00220       j = right ;
00221 
00222       /*----- partition:  move elements bigger than pivot up and elements
00223                           smaller than pivot down, scanning in from ends -----*/
00224 
00225       do{
00226         for( ; a[++i] < pivot ; ) ;  /* scan i up,   until a[i] >= pivot */
00227         for( ; a[--j] > pivot ; ) ;  /* scan j down, until a[j] <= pivot */
00228 
00229         if( j <= i ) break ;         /* if j meets i, quit */
00230 
00231         QS_SWAP( a[i] , a[j] ) ;
00232       } while( 1 ) ;
00233 
00234       /*----- at this point, the array is partitioned -----*/
00235 
00236       a[right] = a[i] ;              /* restore the pivot */
00237       a[i]     = pivot ;
00238 
00239       if( (i-left)  > cutoff ){ stack[mst++] = left ; stack[mst++] = i-1   ; }
00240       if( (right-i) > cutoff ){ stack[mst++] = i+1  ; stack[mst++] = right ; }
00241 
00242    }  /* end of while stack is non-empty */
00243 }
00244 
00245 /******************************************************************************/
00246 /* quick_sort :  sort an array partially recursively, and partially insertion */
00247 
00248 #ifndef QS_CUTOFF
00249 #define QS_CUTOFF 40
00250 #endif
00251 
00252 void qsort_sh( int n , short * a )
00253 {
00254    if( n < QS_CUTOFF )
00255       qsrec_sh( n , a , QS_CUTOFF ) ;
00256 
00257    isort_sh( n , a ) ;
00258 }