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