00001
00002
00003
00004
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 ) ;
00013
00014
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
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
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
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
00084
00085 for( krep=0 ; krep < nrep ; krep++ ){
00086
00087 kdiff = 0 ;
00088
00089 for( kk=0 ; kk < nz ; kk++ ){
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] ;
00097
00098 if( clipval > 0 && val <= clipval )
00099 val = shout[ijk] = 0 ;
00100
00101 if( keep_zero && val == 0 ) continue ;
00102
00103 for( dd=0 ; dd < nd ; dd++ ){
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
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
00130
00131 KILL_CLUSTER(cl) ;
00132 free(tmp) ;
00133 if( shin != DSET_BRICK_ARRAY(inset,0) ) free(shin) ;
00134
00135
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 ;
00155 register short temp ;
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] ){
00163 p = j ;
00164 temp = a[j] ;
00165
00166 do{
00167 a[p] = a[p-1] ;
00168 p-- ;
00169 } while( p > 0 && temp < a[p-1] ) ;
00170
00171 a[p] = temp ;
00172 }
00173 }
00174 }
00175
00176
00177
00178
00179
00180 #define QS_STACK 1024
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 ;
00186 register short temp , pivot ;
00187 register short * a = ar ;
00188
00189 int left , right , mst , stack[QS_STACK] ;
00190
00191
00192
00193 if( cutoff < 3 ) cutoff = 3 ;
00194 if( n < cutoff ) return ;
00195
00196
00197
00198 stack[0] = 0 ;
00199 stack[1] = n-1 ;
00200 mst = 2 ;
00201
00202
00203
00204 while( mst > 0 ){
00205 right = stack[--mst] ;
00206 left = stack[--mst] ;
00207
00208 i = ( left + right ) / 2 ;
00209
00210
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] ;
00217 a[i] = a[right] ;
00218
00219 i = left ;
00220 j = right ;
00221
00222
00223
00224
00225 do{
00226 for( ; a[++i] < pivot ; ) ;
00227 for( ; a[--j] > pivot ; ) ;
00228
00229 if( j <= i ) break ;
00230
00231 QS_SWAP( a[i] , a[j] ) ;
00232 } while( 1 ) ;
00233
00234
00235
00236 a[right] = a[i] ;
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 }
00243 }
00244
00245
00246
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 }