Doxygen Source Code Documentation
thd_winsor.c File Reference
#include "mrilib.h"Go to the source code of this file.
Defines | |
| #define | CFRAC 0.20 |
| #define | LIM(x, b, t) ((x<b) ? b : (x>t) ? t : x) |
| #define | QS_STACK 1024 |
| #define | QS_SWAP(x, y) (temp=(x),(x)=(y),(y)=temp) |
| #define | QS_CUTOFF 40 |
Functions | |
| void | qsort_sh (int n, short *a) |
| THD_3dim_dataset * | WINsorize (THD_3dim_dataset *inset, int nrep, int cbot, int ctop, float irad, char *prefix, int keep_zero, int clipval, byte *mask) |
| void | isort_sh (int n, short *ar) |
| void | qsrec_sh (int n, short *ar, int cutoff) |
Define Documentation
|
|
Definition at line 9 of file thd_winsor.c. Referenced by WINsorize(). |
|
|
Definition at line 10 of file thd_winsor.c. Referenced by WINsorize(). |
|
|
Definition at line 249 of file thd_winsor.c. Referenced by qsort_sh(). |
|
|
Definition at line 180 of file thd_winsor.c. Referenced by qsrec_sh(). |
|
|
Definition at line 181 of file thd_winsor.c. |
Function Documentation
|
||||||||||||
|
Definition at line 152 of file thd_winsor.c. Referenced by qsort_partial_sh(), and qsort_sh().
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 }
|
|
||||||||||||
|
Definition at line 252 of file thd_winsor.c. References a, isort_sh(), QS_CUTOFF, and qsrec_sh(). Referenced by WINsorize().
|
|
||||||||||||||||
|
Definition at line 183 of file thd_winsor.c. References a, i, left, QS_STACK, QS_SWAP, and right. Referenced by qsort_partial_sh(), and qsort_sh().
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 }
|
|
||||||||||||||||||||||||||||||||||||||||
|
Definition at line 16 of file thd_winsor.c. References abs, ADDTO_CLUSTER, ADN_none, ADN_ntt, ADN_nvals, ADN_prefix, CFRAC, DSET_BRICK_ARRAY, DSET_BRICK_TYPE, DSET_load, DSET_NX, DSET_NY, DSET_NZ, EDIT_dset_items(), EDIT_empty_copy(), EDIT_substitute_brick(), free, MCW_cluster::i, MCW_cluster::j, MCW_cluster::k, KILL_CLUSTER, LIM, malloc, MCW_build_mask(), MCW_cluster::num_pt, nz, qsort_sh(), and THD_filename_ok(). Referenced by main().
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 }
|
