1 #include "mrilib.h"
2
3 #define SKIPVOX(i) ( mask != NULL && !mask[i] )
4
5 static int use_dxyz = 0 ;
mri_medianfilter_usedxyz(int i)6 void mri_medianfilter_usedxyz( int i ){ use_dxyz = i; }
mri_flatfilter_usedxyz(int i)7 void mri_flatfilter_usedxyz ( int i ){ use_dxyz = i; }
8
9 /*-----------------------------------------------------------------------*/
10 /*! Compute the median filter of an input image.
11 Output image is always in float format.
12 -------------------------------------------------------------------------*/
13
mri_medianfilter(MRI_IMAGE * imin,float irad,byte * mask,int verb)14 MRI_IMAGE *mri_medianfilter( MRI_IMAGE *imin, float irad, byte *mask, int verb )
15 {
16 MRI_IMAGE *imout ;
17 float *fin=NULL , *fout , *tmp ;
18 short *sin=NULL ; byte *bin=NULL ; void *vin ;
19 short *di , *dj , *dk ;
20 int nd, ii,jj,kk, ip,jp,kp, nx,ny,nz, nxy, ijk, dd,nt=0,pjk, kd=0 ;
21 MCW_cluster *cl ;
22 float dz ;
23
24 ENTRY("mri_medianfilter") ;
25
26 if( imin == NULL || irad <= 0.0f ) RETURN(NULL) ;
27
28 /** deal with vector-valued images [15 Dec 2008] -- see mrilib.h **/
29
30 #undef CALLME
31 #define CALLME(inn,out) (out) = mri_medianfilter( (inn), irad,mask,verb )
32 if( ISVECTIM(imin) ){ VECTORME(imin,imout) ; RETURN(imout) ; }
33
34 /** if not a good input data type, floatize and try again **/
35
36 if( imin->kind != MRI_float &&
37 imin->kind != MRI_short &&
38 imin->kind != MRI_byte ){
39
40 MRI_IMAGE *qim ;
41 qim = mri_to_float( imin ) ;
42 imout = mri_medianfilter( qim , irad , mask , verb ) ;
43 mri_free( qim ) ;
44 RETURN(imout) ;
45 }
46
47 /** build cluster of points for median-izing **/
48
49 if( !use_dxyz ){
50 if( irad < 1.01f ) irad = 1.01f ;
51 dz = (imin->nz == 1) ? 6666.0f : 1.0f ;
52 cl = MCW_build_mask( 1.0f,1.0f,dz , irad ) ;
53 } else {
54 float dm ;
55 dz = (imin->nz == 1) ? 6666.0f : imin->dz ;
56 dm = MIN(imin->dx,imin->dy) ; dm = MIN(dm,dz) ;
57 dm *= 1.01f ; if( irad < dm ) irad = dm ;
58 cl = MCW_build_mask( imin->dx,imin->dy,dz , irad ) ;
59 }
60
61 if( cl == NULL || cl->num_pt < 6 ){ KILL_CLUSTER(cl); RETURN(NULL); }
62
63 ADDTO_CLUSTER(cl,0,0,0,0) ;
64
65 di = cl->i ; dj = cl->j ; dk = cl->k ; nd = cl->num_pt ;
66 nx = imin->nx ; ny = imin->ny ; nz = imin->nz ; nxy = nx*ny ;
67
68 if( verb ){
69 INFO_message("Median filter mask has %d voxels",nd) ;
70 if( mask != NULL )
71 ININFO_message(" Data mask has %d voxels",THD_countmask(nxy*nz,mask)) ;
72 }
73
74 imout = mri_new_conforming( imin , MRI_float ) ;
75 fout = MRI_FLOAT_PTR( imout ) ;
76 vin = mri_data_pointer( imin ) ;
77 tmp = (float *) malloc(sizeof(float)*nd) ;
78 switch( imin->kind ){
79 default: break ;
80 case MRI_float: fin = (float *)vin ; break ;
81 case MRI_short: sin = (short *)vin ; break ;
82 case MRI_byte : bin = (byte *)vin ; break ;
83 }
84
85 if( verb ){
86 kd = (int)rint(0.05*nz); if( kd < 1 ) kd = 1;
87 fprintf(stderr," + Median filter loop") ;
88 }
89
90 for( kk=0 ; kk < nz ; kk++ ){
91 if( verb && kk%kd == 0 ) fprintf(stderr,".") ;
92 for( jj=0 ; jj < ny ; jj++ ){
93 for( ii=0 ; ii < nx ; ii++ ){
94 ijk = ii + jj*nx + kk*nxy ;
95 if( SKIPVOX(ijk) ) continue ;
96
97 /* extract neighborhood values */
98
99 switch( imin->kind ){
100 default: break ;
101 case MRI_float:
102 for( nt=dd=0 ; dd < nd ; dd++ ){
103 ip = ii+di[dd] ; if( ip < 0 || ip >= nx ) continue ;
104 jp = jj+dj[dd] ; if( jp < 0 || jp >= ny ) continue ;
105 kp = kk+dk[dd] ; if( kp < 0 || kp >= nz ) continue ;
106 pjk = ip+jp*nx+kp*nxy ;
107 if( SKIPVOX(pjk) ) continue ;
108 tmp[nt++] = fin[pjk] ;
109 }
110 break ;
111 case MRI_short:
112 for( nt=dd=0 ; dd < nd ; dd++ ){
113 ip = ii+di[dd] ; if( ip < 0 || ip >= nx ) continue ;
114 jp = jj+dj[dd] ; if( jp < 0 || jp >= ny ) continue ;
115 kp = kk+dk[dd] ; if( kp < 0 || kp >= nz ) continue ;
116 pjk = ip+jp*nx+kp*nxy ;
117 if( SKIPVOX(pjk) ) continue ;
118 tmp[nt++] = sin[pjk] ;
119 }
120 break ;
121 case MRI_byte:
122 for( nt=dd=0 ; dd < nd ; dd++ ){
123 ip = ii+di[dd] ; if( ip < 0 || ip >= nx ) continue ;
124 jp = jj+dj[dd] ; if( jp < 0 || jp >= ny ) continue ;
125 kp = kk+dk[dd] ; if( kp < 0 || kp >= nz ) continue ;
126 pjk = ip+jp*nx+kp*nxy ;
127 if( SKIPVOX(pjk) ) continue ;
128 tmp[nt++] = bin[pjk] ;
129 }
130 break ;
131 }
132
133 fout[ijk] = qmed_float( nt , tmp ) ; /* the actual median-izing */
134 }}}
135 if( verb ) fprintf(stderr,"\n") ;
136
137 KILL_CLUSTER(cl); free((void *)tmp); /* toss the trash */
138 RETURN(imout) ;
139 }
140
141 /*-----------------------------------------------------------------------*/
142 /*! Compute the flat (local average) filter of an input image.
143 Output image is always in float format.
144 -------------------------------------------------------------------------*/
145
mri_flatfilter(MRI_IMAGE * imin,float irad,byte * mask,int verb)146 MRI_IMAGE *mri_flatfilter( MRI_IMAGE *imin, float irad, byte *mask, int verb )
147 {
148 MRI_IMAGE *imout ;
149 float *fin=NULL , *fout , *tmp ;
150 short *sin=NULL ; byte *bin=NULL ; void *vin ;
151 short *di , *dj , *dk ;
152 int nd, ii,jj,kk, ip,jp,kp, nx,ny,nz, nxy, ijk, dd,nt=0,pjk, kd=0 ;
153 MCW_cluster *cl ;
154 float dz , sum ;
155
156 ENTRY("mri_flatfilter") ;
157
158 if( imin == NULL || irad <= 0.0f ) RETURN(NULL) ;
159
160 /** deal with vector-valued images [15 Dec 2008] -- see mrilib.h **/
161
162 #undef CALLME
163 #define CALLME(inn,out) (out) = mri_flatfilter( (inn), irad,mask,verb )
164 if( ISVECTIM(imin) ){ VECTORME(imin,imout) ; RETURN(imout) ; }
165
166 /** if not a good input data type, floatize and try again **/
167
168 if( imin->kind != MRI_float &&
169 imin->kind != MRI_short &&
170 imin->kind != MRI_byte ){
171
172 MRI_IMAGE *qim ;
173 qim = mri_to_float( imin ) ;
174 imout = mri_flatfilter( qim , irad , mask , verb ) ;
175 mri_free( qim ) ;
176 RETURN(imout) ;
177 }
178
179 /** build cluster of points for filter-izing **/
180
181 if( !use_dxyz ){
182 if( irad < 1.01f ) irad = 1.01f ;
183 dz = (imin->nz == 1) ? 6666.0f : 1.0f ;
184 cl = MCW_build_mask( 1.0f,1.0f,dz , irad ) ;
185 } else {
186 float dm ;
187 dz = (imin->nz == 1) ? 6666.0f : imin->dz ;
188 dm = MIN(imin->dx,imin->dy) ; dm = MIN(dm,dz) ;
189 dm *= 1.01f ; if( irad < dm ) irad = dm ;
190 cl = MCW_build_mask( imin->dx,imin->dy,dz , irad ) ;
191 }
192
193 if( cl == NULL || cl->num_pt < 6 ){ KILL_CLUSTER(cl); RETURN(NULL); }
194
195 ADDTO_CLUSTER(cl,0,0,0,0) ; /* add central point to cluster */
196
197 di = cl->i ; dj = cl->j ; dk = cl->k ; nd = cl->num_pt ;
198 nx = imin->nx ; ny = imin->ny ; nz = imin->nz ; nxy = nx*ny ;
199
200 if( verb ){
201 fprintf(stderr,"++ Flat filter mask size=%d voxels",nd) ;
202 if( mask != NULL )
203 fprintf(stderr," Data mask=%d",THD_countmask(nxy*nz,mask)) ;
204 }
205
206 imout = mri_new_conforming( imin , MRI_float ) ; /* zero filled */
207 fout = MRI_FLOAT_PTR( imout ) ;
208 vin = mri_data_pointer( imin ) ;
209 tmp = (float *) malloc(sizeof(float)*nd) ;
210 switch( imin->kind ){
211 default: break ;
212 case MRI_float: fin = (float *)vin ; break ;
213 case MRI_short: sin = (short *)vin ; break ;
214 case MRI_byte : bin = (byte *)vin ; break ;
215 }
216
217 if( verb ){ kd = (int)rint(0.05*nz); if( kd < 1 ) kd = 1; }
218
219 for( kk=0 ; kk < nz ; kk++ ){
220 if( verb && kk%kd == 0 ) fprintf(stderr,".") ;
221 for( jj=0 ; jj < ny ; jj++ ){
222 for( ii=0 ; ii < nx ; ii++ ){
223 ijk = ii + jj*nx + kk*nxy ;
224 if( SKIPVOX(ijk) ) continue ;
225
226 /* extract neighborhood values */
227
228 switch( imin->kind ){
229 default: break ;
230 case MRI_float:
231 for( nt=dd=0 ; dd < nd ; dd++ ){
232 ip = ii+di[dd] ; if( ip < 0 || ip >= nx ) continue ;
233 jp = jj+dj[dd] ; if( jp < 0 || jp >= ny ) continue ;
234 kp = kk+dk[dd] ; if( kp < 0 || kp >= nz ) continue ;
235 pjk = ip+jp*nx+kp*nxy ;
236 if( SKIPVOX(pjk) ) continue ;
237 tmp[nt++] = fin[pjk] ;
238 }
239 break ;
240 case MRI_short:
241 for( nt=dd=0 ; dd < nd ; dd++ ){
242 ip = ii+di[dd] ; if( ip < 0 || ip >= nx ) continue ;
243 jp = jj+dj[dd] ; if( jp < 0 || jp >= ny ) continue ;
244 kp = kk+dk[dd] ; if( kp < 0 || kp >= nz ) continue ;
245 pjk = ip+jp*nx+kp*nxy ;
246 if( SKIPVOX(pjk) ) continue ;
247 tmp[nt++] = sin[pjk] ;
248 }
249 break ;
250 case MRI_byte:
251 for( nt=dd=0 ; dd < nd ; dd++ ){
252 ip = ii+di[dd] ; if( ip < 0 || ip >= nx ) continue ;
253 jp = jj+dj[dd] ; if( jp < 0 || jp >= ny ) continue ;
254 kp = kk+dk[dd] ; if( kp < 0 || kp >= nz ) continue ;
255 pjk = ip+jp*nx+kp*nxy ;
256 if( SKIPVOX(pjk) ) continue ;
257 tmp[nt++] = bin[pjk] ;
258 }
259 break ;
260 }
261
262 sum = 0.0f ;
263 for( dd=0 ; dd < nt ; dd++ ) sum += tmp[dd] ;
264 fout[ijk] = sum / nt ; /* the actual filtering */
265 }}}
266 if( verb ) fprintf(stderr,"\n") ;
267
268 KILL_CLUSTER(cl); free((void *)tmp); /* toss the trash */
269 RETURN(imout) ;
270 }
271