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 /*** NOT 7D SAFE ***/
10
11 /*** prototype ***/
12
13 void nonmax_kill( int , int , float * ) ;
14
15 /*-------------------------------------------------------------------------
16 Do Sobel edge enhancement on an image. Output image will be floats.
17 nloc = distance over which to suppress filter outputs which are not
18 local maxima: _individual applies to each separate filter,
19 _collective applies to the merger of all 4 filters;
20 (either or both may be zero).
21 ---------------------------------------------------------------------------*/
22
mri_sobel(int nloc_individual,int nloc_collective,MRI_IMAGE * imin)23 MRI_IMAGE * mri_sobel( int nloc_individual , int nloc_collective , MRI_IMAGE *imin )
24 {
25 MRI_IMAGE *imfl , *imout ;
26 float *flin , *flout , *fjj,*fjp,*fjm ;
27 int ii , jj , joff , nx,ny , ij , nij , nloc ;
28 register float sxx,syy,sdd,see ;
29
30 MRI_IMAGE *imxx , *imyy , *imdd , *imee ;
31 float *flxx , *flyy , *fldd , *flee ;
32
33 #ifdef DEBUG
34 printf("Entry: mri_sobel\n") ;
35 #endif
36
37 /*** setup input and output images ***/
38
39 if( imin == NULL || ! MRI_IS_2D(imin) ){
40 fprintf(stderr,"\n*** mri_sobel only works on 2D images!\n") ;
41 EXIT(1) ;
42 }
43
44 nx = imin->nx ;
45 ny = imin->ny ;
46 if( imin->kind == MRI_float ) imfl = imin ;
47 else imfl = mri_to_float( imin ) ;
48
49 imout = mri_new( nx , ny , MRI_float ) ;
50 flout = mri_data_pointer( imout ) ;
51 flin = mri_data_pointer( imfl ) ;
52
53 /*** setup arrays for each direction ***/
54
55 imxx = mri_new( nx , ny , MRI_float ) ; flxx = mri_data_pointer( imxx ) ;
56 imyy = mri_new( nx , ny , MRI_float ) ; flyy = mri_data_pointer( imyy ) ;
57 imdd = mri_new( nx , ny , MRI_float ) ; fldd = mri_data_pointer( imdd ) ;
58 imee = mri_new( nx , ny , MRI_float ) ; flee = mri_data_pointer( imee ) ;
59
60 /*** initialize edges to be zeros ***/
61
62 for( jj=0 ; jj < ny ; jj++ ){
63 joff = nx * jj ;
64 flxx[joff] = flyy[joff] = fldd[joff] = flee[joff] = 0;
65 flxx[joff+nx-1] = flyy[joff+nx-1] = fldd[joff+nx-1] = flee[joff+nx-1]= 0;
66 }
67 for( ii=0 ; ii < nx ; ii++ ){
68 flxx[ii] = flyy[ii] = fldd[ii] = flee[ii] = 0;
69 flxx[joff+ii] = flyy[joff+ii] = fldd[joff+ii] = flee[joff+ii] = 0;
70 }
71
72 /*** do scan over interior points to produce images for each direction ***/
73
74 for( jj=1 ; jj < ny-1 ; jj++ ){
75
76 joff = nx * jj ;
77 fjj = &flin[joff] ;
78 fjm = &flin[joff-nx] ;
79 fjp = &flin[joff+nx] ;
80
81 for( ii=1 ; ii < nx-1 ; ii++ ){
82
83 sxx = 2.0 * ( fjj[ii+1] - fjj[ii-1] )
84 + ( fjp[ii+1] + fjm[ii+1] - fjp[ii-1] - fjm[ii-1] ) ;
85
86 syy = 2.0 * ( fjp[ii] - fjm[ii] )
87 + ( fjp[ii+1] + fjp[ii-1] - fjm[ii+1] - fjm[ii-1] ) ;
88
89 sdd = 2.0 * ( fjp[ii+1] - fjm[ii-1] )
90 + ( fjj[ii+1] + fjp[ii] - fjj[ii-1] - fjm[ii] ) ;
91
92 see = 2.0 * ( fjp[ii-1] - fjm[ii+1] )
93 + ( fjp[ii] + fjj[ii-1] - fjm[ii] - fjj[ii+1] ) ;
94
95 flxx[ii+joff] = fabs(sxx) ; flyy[ii+joff] = fabs(syy) ;
96 fldd[ii+joff] = fabs(sdd) ; flee[ii+joff] = fabs(see) ;
97 }
98 }
99
100 if( imfl != imin ) mri_free( imfl ) ; /* don't need anymore */
101
102 /*.......................................................................*/
103 /*** if nloc > 0, scan over +/- nloc points and kill non maxima points ***/
104
105 nloc = nloc_individual ;
106 if( nloc > 0 ){
107
108 /*** do xx ***/
109
110 for( jj=1 ; jj < ny-1 ; jj++ ){
111 joff = jj * nx ;
112 nonmax_kill( nloc , nx , &flxx[joff] ) ;
113 }
114
115 /*** do yy ***/
116
117 for( ii=1 ; ii < nx-1 ; ii++ ){
118 for( jj=0 ; jj < ny ; jj++ ) flout[jj] = flyy[jj*nx+ii] ;
119 nonmax_kill( nloc , ny , flout ) ;
120 for( jj=0 ; jj < ny ; jj++ ) flyy[jj*nx+ii] = flout[jj] ;
121 }
122
123 /*** do dd: here, ij = ii-jj ***/
124
125 for( ij = -(ny-3) ; ij < nx-2 ; ij++ ){
126
127 ii = MAX( 0 , ij ) ; jj = ii - ij ;
128
129 for( nij=0 ; (ii<nx) && (jj<ny) ; nij++ , ii++ , jj++ ){
130 flout[nij] = fldd[ii+nx*jj] ;
131 }
132
133 nonmax_kill( nloc , nij , flout ) ;
134
135 ii = MAX( 0 , ij ) ; jj = ii - ij ;
136
137 for( nij=0 ; (ii<nx) && (jj<ny) ; nij++ , ii++ , jj++ ){
138 fldd[ii+nx*jj] = flout[nij] ;
139 }
140 }
141
142 /*** do ee: here, ij = ii+jj ***/
143
144 for( ij = 2 ; ij < (nx+ny-4) ; ij++ ){
145
146 jj = MIN( ij , ny-1 ) ; ii = ij - jj ;
147
148 for( nij=0 ; (ii<nx) && (jj>=0) ; nij++ , ii++ , jj-- ){
149 flout[nij] = flee[ii+nx*jj] ;
150 }
151
152 nonmax_kill( nloc , nij , flout ) ;
153
154 jj = MIN( ij , ny-1 ) ; ii = ij - jj ;
155
156 for( nij=0 ; (ii<nx) && (jj>=0) ; nij++ , ii++ , jj-- ){
157 flee[ii+nx*jj] = flout[nij] ;
158 }
159 }
160
161 } /* end if( nloc > 0 ) */
162
163 /*.......................................................................*/
164 /*** assign maximum of outputs at a given location to the final result ***/
165
166 for( jj=0 ; jj < ny ; jj++ ){ /* clear edges */
167 joff = nx * jj ;
168 flout[joff] = flout[joff+nx-1] = 0.0 ;
169 }
170 for( ii=0 ; ii < nx ; ii++ ){
171 flout[ii] = flout[joff+ii] = 0.0 ;
172 }
173
174 for( jj=1 ; jj < ny-1 ; jj++ ){
175 joff = nx * jj ;
176 for( ii=1 ; ii < nx-1 ; ii++ ){
177 sxx = MAX( flxx[ii+joff] , flyy[ii+joff] ) ;
178 sdd = MAX( fldd[ii+joff] , flee[ii+joff] ) ;
179 flout[ii+joff] = MAX( sxx , sdd ) ;
180 }
181 }
182
183 /*.......................................................................*/
184
185 nloc = nloc_collective ;
186 if( nloc > 0 ){
187 int xx,yy , xbot,xtop , ybot,ytop , dx,dy ;
188 float val ;
189
190 for( jj=1 ; jj < ny-1 ; jj++ ){
191 joff = nx * jj ;
192 ybot = MAX(0 ,jj-nloc) ;
193 ytop = MIN(ny-1,jj+nloc) ;
194
195 for( ii=1 ; ii < nx-1 ; ii++ ){
196 xbot = MAX(0 ,ii-nloc) ;
197 xtop = MIN(nx-1,ii+nloc) ;
198
199 sxx = flxx[ii+joff] ; syy = flyy[ii+joff] ;
200 sdd = fldd[ii+joff] ; see = flee[ii+joff] ; val = flout[ii+joff] ;
201
202 if( val == sxx ){ dx = 1 ; dy = 0 ; }
203 else if( val == syy ){ dx = 0 ; dy = 1 ; }
204 else if( val == sdd ){ dx = 1 ; dy = 1 ; }
205 else { dx = 1 ; dy = -1 ; }
206
207 for( xx=ii+dx , yy=jj+dy ;
208 (xx <= xtop) && (yy >= ybot) && (yy <= ytop) ;
209 xx += dx , yy+= dy ) val = MAX(val,flout[xx+nx*yy]) ;
210
211 for( xx=ii-dx , yy=jj-dy ;
212 (xx >= xbot) && (yy >= ybot) && (yy <= ytop) ;
213 xx -= dx , yy-= dy ) val = MAX(val,flout[xx+nx*yy]) ;
214
215 if( flout[ii+joff] < val ) flout[ii+joff] = 0.0 ;
216
217 } /* end for ii */
218 } /* end for jj */
219
220 /*** now, destroy points in flout with no neighbors ***/
221
222 for( ii=0 ; ii < nx*ny ; ii++ ) flee[ii] = 0.0 ;
223
224 for( jj=1 ; jj < ny-1 ; jj++ ){
225 joff = nx * jj ;
226 for( ii=1 ; ii < nx-1 ; ii++ ){
227 xx = ii+joff ;
228 if( flout[xx] == 0.0 ) continue ;
229
230 if( flout[xx-1] != 0.0 || flout[xx+1] != 0.0 ||
231 flout[xx+nx] != 0.0 || flout[xx-nx] != 0.0 ||
232 flout[xx+nx+1] != 0.0 || flout[xx+nx-1] != 0.0 ||
233 flout[xx-nx+1] != 0.0 || flout[xx-nx-1] != 0.0 ) flee[xx] = flout[xx] ;
234 }
235 }
236
237 for( ii=0 ; ii < nx*ny ; ii++ ) flout[ii] = flee[ii] ;
238
239 } /* end if nloc */
240
241 /*........................................................................*/
242
243 mri_free(imxx) ; mri_free(imyy) ; mri_free(imdd) ; mri_free(imee) ;
244
245 MRI_COPY_AUX(imout,imin) ;
246 return imout ;
247 }
248
249 /*--------------------------------------------------------------------------*/
250
nonmax_kill(int nloc,int npt,float * ar)251 void nonmax_kill( int nloc , int npt , float *ar )
252 {
253 int ii ;
254 register int jj , jbot,jtop ;
255 register float val ;
256
257 for( ii=0 ; ii < npt ; ii++ ){
258
259 jbot = MAX( 0 , ii-nloc ) ;
260 jtop = MIN( npt , ii+nloc ) ;
261 val = ar[jbot++] ;
262 for( jj=jbot ; jj < jtop ; jj++ ) val = MAX(val,ar[jj]) ;
263
264 if( ar[ii] != val ) ar[ii] = 0.0 ;
265 }
266
267 return ;
268 }
269
270 /*----------------------------------------------------------------
271 Sharpen an image; method same as "pgmenhance".
272 phi is a parameter between 0.0 and 1.0: larger --> more sharp.
273 logify is a flag: nonzero --> take log, sharpen, then exp
274 (e.g., homomorphic filtering).
275 ------------------------------------------------------------------*/
276
mri_sharpen(float phi,int logify,MRI_IMAGE * im)277 MRI_IMAGE * mri_sharpen( float phi , int logify , MRI_IMAGE *im )
278 {
279 int ii,jj , nx , ny , joff,ijoff , npix ;
280 MRI_IMAGE *flim , *outim ;
281 float *flar , *outar ;
282 float nphi , omphi , sum , bot,top ;
283
284 ENTRY("mri_sharpen") ;
285
286 if( phi <= 0.0 || phi >= 1.0 ){
287 ERROR_message("mri_sharpen: illegal phi=%g\n",phi) ;
288 phi = (phi <= 0.0) ? 0.1 : 0.9 ;
289 }
290
291 if( im->kind == MRI_float && !logify ) flim = im ;
292 else flim = mri_to_float( im ) ;
293 flar = mri_data_pointer( flim ) ;
294
295 nx = flim->nx ; ny = flim->ny ; npix = nx*ny ;
296 outim = mri_new( nx , ny , MRI_float ) ;
297 outar = mri_data_pointer( outim ) ;
298
299 if( logify ){
300 for( ii=0 ; ii < npix ; ii++ ) flar[ii] = log( fabs(flar[ii])+1.0 ) ;
301 }
302
303 for( ii=0 ; ii < nx ; ii++ ) outar[ii] = flar[ii] ; /* copy 1st row */
304
305 nphi = phi / 9.0 ;
306 omphi = 1.0/(1.0-phi) ;
307 bot = mri_min(flim) ;
308 top = mri_max(flim) ;
309
310 for( jj=1 ; jj < ny-1 ; jj++ ){
311 joff = jj * nx ;
312
313 outar[joff] = flar[joff] ; /* copy 1st and last columns */
314 outar[joff+nx-1] = flar[joff+nx-1] ;
315
316 for( ii=1 ; ii < nx-1 ; ii++ ){ /* filter all intermediate points */
317 ijoff = joff + ii ;
318
319 sum = flar[ijoff-nx-1] + flar[ijoff-nx] + flar[ijoff-nx+1]
320 +flar[ijoff-1] + flar[ijoff] + flar[ijoff+1]
321 +flar[ijoff+nx-1] + flar[ijoff+nx] + flar[ijoff+nx+1] ;
322
323 outar[ijoff] = (flar[ijoff] - nphi*sum) * omphi ;
324
325 if( outar[ijoff] < bot ) outar[ijoff] = bot ;
326 else if( outar[ijoff] > top ) outar[ijoff] = top ;
327 }
328 }
329
330 joff = (ny-1) * nx ;
331 for( ii=0 ; ii < nx ; ii++ ) outar[joff+ii] = flar[joff+ii] ; /* copy last row */
332
333 if( logify ){
334 for( ii=0 ; ii < npix ; ii++ ) outar[ii] = exp(outar[ii]) ;
335 }
336
337 if( flim != im ) mri_free( flim ) ;
338 RETURN(outim) ;
339 }
340