1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "mrilib.h"
8
9 #define CFRAC 0.20
10 #define LIM(x,b,t) ((x<b) ? b : (x>t) ? t : x)
11
12 void qsort_sh( int n , short * a ) ; /* at end of file */
13
14 /*-------- 06 Jul 2000 - RWCox ----------------------------------------*/
15
WINsorize(THD_3dim_dataset * inset,int nrep,int cbot,int ctop,float irad,char * prefix,int keep_zero,int clipval,byte * mask)16 THD_3dim_dataset * WINsorize( THD_3dim_dataset *inset ,
17 int nrep , int cbot , int ctop ,
18 float irad , char *prefix ,
19 int keep_zero , int clipval , byte *mask )
20 {
21 THD_3dim_dataset *outset ;
22 short *shin , *shout , *di,*dj,*dk , *tmp , val,nval ;
23 MCW_cluster *cl ;
24 int jj,kk , krep,kdiff, nx,ny,nz,nxy,nxyz , nd,dd ;
25 int ip,jp,kp , nx1,ny1,nz1 , verb=1 ;
26 int nrep_until ;
27 register int ii,ijk ;
28
29 /*- check inputs -*/
30
31 if( inset == NULL || DSET_BRICK_TYPE(inset,0) != MRI_short ) return NULL ;
32 DSET_load(inset) ;
33 if( DSET_BRICK_ARRAY(inset,0) == NULL ) return NULL ;
34
35 if( nrep == 0 ) return NULL ;
36
37 if( nrep < 0 ){ nrep_until = abs(nrep) ; nrep = 999 ; }
38 else { nrep_until = 2 ; }
39
40 if( irad < 0.0 ){ verb=0 ; irad = -irad ; }
41
42 if( irad < 1.01 ) irad = 1.01 ;
43 if( !THD_filename_ok(prefix) ) prefix = "Winsor" ;
44
45 /*- build list of points to use -*/
46
47 cl = MCW_build_mask( 1.0,1.0,1.0 , irad ) ;
48
49 if( cl == NULL || cl->num_pt < 6 ){ KILL_CLUSTER(cl); return NULL; }
50
51 ADDTO_CLUSTER(cl,0,0,0,0) ;
52
53 di = cl->i ; dj = cl->j ; dk = cl->k ; nd = cl->num_pt ;
54
55 if( verb ) fprintf(stderr,"+++ WINsorize irad=%f nbhd=%d\n",irad,nd) ;
56
57 /*- make output array -*/
58
59 nx = DSET_NX(inset) ; nx1 = nx-1 ;
60 ny = DSET_NY(inset) ; ny1 = ny-1 ;
61 nz = DSET_NZ(inset) ; nz1 = nz-1 ; nxy = nx*ny ; nxyz = nxy*nz ;
62
63 shout = (short *) malloc(sizeof(short)*nxyz) ;
64 tmp = (short *) malloc(sizeof(short)*nd) ;
65
66 if( nrep > 1 ){
67 shin = (short *) malloc(sizeof(short)*nxyz) ;
68 memcpy( shin , DSET_BRICK_ARRAY(inset,0) , sizeof(short)*nxyz ) ;
69 } else {
70 shin = DSET_BRICK_ARRAY(inset,0) ;
71 }
72
73 if( cbot <= 0 || cbot >= nd-1 ){
74 cbot = rint( CFRAC*nd ) ;
75 if( cbot <= 0 ) cbot = 1 ;
76 if( verb ) fprintf(stderr,"+++ WINsorize cbot=%d\n",cbot) ;
77 }
78 if( ctop <= cbot || cbot >= nd-1 ){
79 ctop = nd-1-cbot ;
80 if( verb ) fprintf(stderr,"+++ WINsorize ctop=%d\n",ctop) ;
81 }
82
83 /*- do the work -*/
84
85 for( krep=0 ; krep < nrep ; krep++ ){
86
87 kdiff = 0 ; /* count of how many voxels were changed */
88
89 for( kk=0 ; kk < nz ; kk++ ){ /* loops over 3D voxel indices */
90 for( jj=0 ; jj < ny ; jj++ ){
91 ijk = jj*nx+kk*nxy ;
92 for( ii=0 ; ii < nx ; ii++,ijk++ ){
93
94 if( mask != NULL && !mask[ijk] ){ shout[ijk]=shin[ijk]; continue; }
95
96 val = shin[ijk] ; /* current voxel */
97
98 if( clipval > 0 && val <= clipval ) /* 19 Oct 2001 */
99 val = shout[ijk] = 0 ;
100
101 if( keep_zero && val == 0 ) continue ; /* don't filter 0 */
102
103 for( dd=0 ; dd < nd ; dd++ ){ /* loop over nbhd */
104 ip = ii+di[dd] ; ip = LIM(ip,0,nx1) ;
105 jp = jj+dj[dd] ; jp = LIM(jp,0,ny1) ;
106 kp = kk+dk[dd] ; kp = LIM(kp,0,nz1) ;
107 tmp[dd] = shin[ip+jp*nx+kp*nxy] ;
108 }
109
110 qsort_sh( nd , tmp ) ;
111
112 shout[ijk] = nval = LIM(val,tmp[cbot],tmp[ctop]) ;
113
114 if( nval != val ) kdiff++ ;
115 }
116 }
117 }
118
119 /* prepare for next iteration */
120
121 if( verb ) fprintf(stderr,"+++ WINsorize iter%2d: # changed=%d\n",krep+1,kdiff) ;
122
123 if( kdiff < nrep_until ) break ;
124
125 if( krep < nrep-1 )
126 memcpy( shin , shout , sizeof(short)*nxyz ) ;
127 }
128
129 /*- toss the trashola */
130
131 KILL_CLUSTER(cl) ;
132 free(tmp) ;
133 if( shin != DSET_BRICK_ARRAY(inset,0) ) free(shin) ;
134
135 /*- make output dataset */
136
137 outset = EDIT_empty_copy( inset ) ;
138
139 EDIT_dset_items( outset ,
140 ADN_prefix , prefix ,
141 ADN_nvals , 1 ,
142 ADN_ntt , 0 ,
143 ADN_none ) ;
144
145 EDIT_substitute_brick( outset , 0 , MRI_short , shout ) ;
146
147 return outset ;
148 }
149
150 /*****************************************************************************/
151
isort_sh(int n,short * ar)152 void isort_sh( int n , short * ar )
153 {
154 register int j , p ; /* array indices */
155 register short temp ; /* a[j] holding place */
156 register short * a = ar ;
157
158 if( n < 2 ) return ;
159
160 for( j=1 ; j < n ; j++ ){
161
162 if( a[j] < a[j-1] ){ /* out of order */
163 p = j ;
164 temp = a[j] ;
165
166 do{
167 a[p] = a[p-1] ; /* now have a[p-1] > temp, so move it up */
168 p-- ;
169 } while( p > 0 && temp < a[p-1] ) ;
170
171 a[p] = temp ; /* finally, put temp in its place */
172 }
173 }
174 }
175
176
177 /*****************************************************************************/
178 /* qsrec : recursive part of quicksort (stack implementation) */
179
180 #define QS_SWAP(x,y) (temp=(x),(x)=(y),(y)=temp)
181 #ifndef QS_STACK
182 # define QS_STACK 9999
183 #endif
184
qsrec_sh(int n,short * ar,int cutoff)185 static void qsrec_sh( int n , short * ar , int cutoff )
186 {
187 register int i , j ; /* scanning indices */
188 register short temp , pivot ; /* holding places */
189 register short * a = ar ;
190
191 int left , right , mst , stack[QS_STACK] ;
192
193 /* return if too short (insertion sort will clean up) */
194
195 if( cutoff < 3 ) cutoff = 3 ;
196 if( n < cutoff ) return ;
197
198 /* initialize stack to start with whole array */
199
200 stack[0] = 0 ;
201 stack[1] = n-1 ;
202 mst = 2 ;
203
204 /* loop while the stack is nonempty */
205
206 while( mst > 0 ){
207 right = stack[--mst] ; /* subarray from left -> right */
208 left = stack[--mst] ;
209
210 i = ( left + right ) / 2 ; /* middle of subarray */
211
212 /* sort the left, middle, and right a[]'s */
213
214 if( a[left] > a[i] ) QS_SWAP( a[left] , a[i] ) ;
215 if( a[left] > a[right] ) QS_SWAP( a[left] , a[right] ) ;
216 if( a[i] > a[right] ) QS_SWAP( a[right] , a[i] ) ;
217
218 pivot = a[i] ; /* a[i] is the median-of-3 pivot! */
219 a[i] = a[right] ;
220
221 i = left ; /* initialize scanning */
222 j = right ;
223
224 /*----- partition: move elements bigger than pivot up and elements
225 smaller than pivot down, scanning in from ends -----*/
226
227 do{
228 for( ; a[++i] < pivot ; ) ; /* scan i up, until a[i] >= pivot */
229 for( ; a[--j] > pivot ; ) ; /* scan j down, until a[j] <= pivot */
230
231 if( j <= i ) break ; /* if j meets i, quit */
232
233 QS_SWAP( a[i] , a[j] ) ;
234 } while( 1 ) ;
235
236 /*----- at this point, the array is partitioned -----*/
237
238 a[right] = a[i] ; /* restore the pivot */
239 a[i] = pivot ;
240
241 if( (i-left) > cutoff ){ stack[mst++] = left ; stack[mst++] = i-1 ; }
242 if( (right-i) > cutoff ){ stack[mst++] = i+1 ; stack[mst++] = right ; }
243
244 } /* end of while stack is non-empty */
245 }
246
247 /******************************************************************************/
248 /* quick_sort : sort an array partially recursively, and partially insertion */
249
250 #ifndef QS_CUTOFF
251 #define QS_CUTOFF 40
252 #endif
253
qsort_sh(int n,short * a)254 void qsort_sh( int n , short * a )
255 {
256 if( n < QS_CUTOFF )
257 qsrec_sh( n , a , QS_CUTOFF ) ;
258
259 isort_sh( n , a ) ;
260 }
261