Doxygen Source Code Documentation
mri_sobel.c File Reference
#include "mrilib.h"Go to the source code of this file.
Functions | |
| void | nonmax_kill (int, int, float *) |
| MRI_IMAGE * | mri_sobel (int nloc_individual, int nloc_collective, MRI_IMAGE *imin) |
| MRI_IMAGE * | mri_sharpen (float phi, int logify, MRI_IMAGE *im) |
Function Documentation
|
||||||||||||||||
|
Definition at line 277 of file mri_sobel.c. References MRI_IMAGE::kind, mri_data_pointer(), mri_free(), mri_max(), mri_min(), mri_new(), mri_to_float(), MRI_IMAGE::nx, MRI_IMAGE::ny, and top. Referenced by ISQ_process_mri(), and mri_sharpen_rgb().
00278 {
00279 int ii,jj , nx , ny , joff,ijoff , npix ;
00280 MRI_IMAGE * flim , * outim ;
00281 float * flar , * outar ;
00282 float nphi , omphi , sum , bot,top ;
00283
00284 #ifdef DEBUG
00285 printf("Entry: mri_sharpen\n") ;
00286 #endif
00287
00288 if( phi <= 0.0 || phi >= 1.0 ){
00289 fprintf(stderr,"*** mri_sharpen: illegal phi=%g\n",phi) ;
00290 return NULL ;
00291 }
00292
00293 if( im->kind == MRI_float && !logify ){
00294 flim = im ;
00295 } else {
00296 flim = mri_to_float( im ) ;
00297 }
00298 flar = mri_data_pointer( flim ) ;
00299
00300 nx = flim->nx ; ny = flim->ny ; npix = nx*ny ;
00301 outim = mri_new( nx , ny , MRI_float ) ;
00302 outar = mri_data_pointer( outim ) ;
00303
00304 if( logify ){
00305 for( ii=0 ; ii < npix ; ii++ ) flar[ii] = log( fabs(flar[ii])+1.0 ) ;
00306 }
00307
00308 for( ii=0 ; ii < nx ; ii++ ) outar[ii] = flar[ii] ; /* copy 1st row */
00309
00310 nphi = phi / 9.0 ;
00311 omphi = 1.0/(1.0-phi) ;
00312 bot = mri_min(flim) ;
00313 top = mri_max(flim) ;
00314
00315 for( jj=1 ; jj < ny-1 ; jj++ ){
00316 joff = jj * nx ;
00317
00318 outar[joff] = flar[joff] ; /* copy 1st and last columns */
00319 outar[joff+nx-1] = flar[joff+nx-1] ;
00320
00321 for( ii=1 ; ii < nx-1 ; ii++ ){ /* filter all intermediate points */
00322 ijoff = joff + ii ;
00323
00324 sum = flar[ijoff-nx-1] + flar[ijoff-nx] + flar[ijoff-nx+1]
00325 +flar[ijoff-1] + flar[ijoff] + flar[ijoff+1]
00326 +flar[ijoff+nx-1] + flar[ijoff+nx] + flar[ijoff+nx+1] ;
00327
00328 outar[ijoff] = (flar[ijoff] - nphi*sum) * omphi ;
00329
00330 if( outar[ijoff] < bot ) outar[ijoff] = bot ;
00331 else if( outar[ijoff] > top ) outar[ijoff] = top ;
00332 }
00333 }
00334
00335 joff = (ny-1) * nx ;
00336 for( ii=0 ; ii < nx ; ii++ ) outar[joff+ii] = flar[joff+ii] ; /* copy last row */
00337
00338 if( logify ){
00339 for( ii=0 ; ii < npix ; ii++ ) outar[ii] = exp(outar[ii]) ;
00340 }
00341
00342
00343 if( flim != im ) mri_free( flim ) ;
00344 return outim ;
00345 }
|
|
||||||||||||||||
|
Definition at line 23 of file mri_sobel.c. References EXIT, MRI_IMAGE::kind, MAX, MIN, MRI_COPY_AUX, mri_data_pointer(), mri_free(), MRI_IS_2D, mri_new(), mri_to_float(), nonmax_kill(), MRI_IMAGE::nx, and MRI_IMAGE::ny. Referenced by ISQ_process_mri().
00024 {
00025 MRI_IMAGE * imfl , * imout ;
00026 float * flin , * flout , * fjj,*fjp,*fjm ;
00027 int ii , jj , joff , nx,ny , ij , nij , nloc ;
00028 register float sxx,syy,sdd,see ;
00029
00030 MRI_IMAGE * imxx , * imyy , * imdd , * imee ;
00031 float * flxx , * flyy , * fldd , * flee ;
00032
00033 #ifdef DEBUG
00034 printf("Entry: mri_sobel\n") ;
00035 #endif
00036
00037 /*** setup input and output images ***/
00038
00039 if( imin == NULL || ! MRI_IS_2D(imin) ){
00040 fprintf(stderr,"\n*** mri_sobel only works on 2D images!\n") ;
00041 EXIT(1) ;
00042 }
00043
00044 nx = imin->nx ;
00045 ny = imin->ny ;
00046 if( imin->kind == MRI_float ) imfl = imin ;
00047 else imfl = mri_to_float( imin ) ;
00048
00049 imout = mri_new( nx , ny , MRI_float ) ;
00050 flout = mri_data_pointer( imout ) ;
00051 flin = mri_data_pointer( imfl ) ;
00052
00053 /*** setup arrays for each direction ***/
00054
00055 imxx = mri_new( nx , ny , MRI_float ) ; flxx = mri_data_pointer( imxx ) ;
00056 imyy = mri_new( nx , ny , MRI_float ) ; flyy = mri_data_pointer( imyy ) ;
00057 imdd = mri_new( nx , ny , MRI_float ) ; fldd = mri_data_pointer( imdd ) ;
00058 imee = mri_new( nx , ny , MRI_float ) ; flee = mri_data_pointer( imee ) ;
00059
00060 /*** initialize edges to be zeros ***/
00061
00062 for( jj=0 ; jj < ny ; jj++ ){
00063 joff = nx * jj ;
00064 flxx[joff] = flyy[joff] = fldd[joff] = flee[joff] = 0;
00065 flxx[joff+nx-1] = flyy[joff+nx-1] = fldd[joff+nx-1] = flee[joff+nx-1]= 0;
00066 }
00067 for( ii=0 ; ii < nx ; ii++ ){
00068 flxx[ii] = flyy[ii] = fldd[ii] = flee[ii] = 0;
00069 flxx[joff+ii] = flyy[joff+ii] = fldd[joff+ii] = flee[joff+ii] = 0;
00070 }
00071
00072 /*** do scan over interior points to produce images for each direction ***/
00073
00074 for( jj=1 ; jj < ny-1 ; jj++ ){
00075
00076 joff = nx * jj ;
00077 fjj = &flin[joff] ;
00078 fjm = &flin[joff-nx] ;
00079 fjp = &flin[joff+nx] ;
00080
00081 for( ii=1 ; ii < nx-1 ; ii++ ){
00082
00083 sxx = 2.0 * ( fjj[ii+1] - fjj[ii-1] )
00084 + ( fjp[ii+1] + fjm[ii+1] - fjp[ii-1] - fjm[ii-1] ) ;
00085
00086 syy = 2.0 * ( fjp[ii] - fjm[ii] )
00087 + ( fjp[ii+1] + fjp[ii-1] - fjm[ii+1] - fjm[ii-1] ) ;
00088
00089 sdd = 2.0 * ( fjp[ii+1] - fjm[ii-1] )
00090 + ( fjj[ii+1] + fjp[ii] - fjj[ii-1] - fjm[ii] ) ;
00091
00092 see = 2.0 * ( fjp[ii-1] - fjm[ii+1] )
00093 + ( fjp[ii] + fjj[ii-1] - fjm[ii] - fjj[ii+1] ) ;
00094
00095 flxx[ii+joff] = fabs(sxx) ; flyy[ii+joff] = fabs(syy) ;
00096 fldd[ii+joff] = fabs(sdd) ; flee[ii+joff] = fabs(see) ;
00097 }
00098 }
00099
00100 if( imfl != imin ) mri_free( imfl ) ; /* don't need anymore */
00101
00102 /*.......................................................................*/
00103 /*** if nloc > 0, scan over +/- nloc points and kill non maxima points ***/
00104
00105 nloc = nloc_individual ;
00106 if( nloc > 0 ){
00107
00108 /*** do xx ***/
00109
00110 for( jj=1 ; jj < ny-1 ; jj++ ){
00111 joff = jj * nx ;
00112 nonmax_kill( nloc , nx , &flxx[joff] ) ;
00113 }
00114
00115 /*** do yy ***/
00116
00117 for( ii=1 ; ii < nx-1 ; ii++ ){
00118 for( jj=0 ; jj < ny ; jj++ ) flout[jj] = flyy[jj*nx+ii] ;
00119 nonmax_kill( nloc , ny , flout ) ;
00120 for( jj=0 ; jj < ny ; jj++ ) flyy[jj*nx+ii] = flout[jj] ;
00121 }
00122
00123 /*** do dd: here, ij = ii-jj ***/
00124
00125 for( ij = -(ny-3) ; ij < nx-2 ; ij++ ){
00126
00127 ii = MAX( 0 , ij ) ; jj = ii - ij ;
00128
00129 for( nij=0 ; (ii<nx) && (jj<ny) ; nij++ , ii++ , jj++ ){
00130 flout[nij] = fldd[ii+nx*jj] ;
00131 }
00132
00133 nonmax_kill( nloc , nij , flout ) ;
00134
00135 ii = MAX( 0 , ij ) ; jj = ii - ij ;
00136
00137 for( nij=0 ; (ii<nx) && (jj<ny) ; nij++ , ii++ , jj++ ){
00138 fldd[ii+nx*jj] = flout[nij] ;
00139 }
00140 }
00141
00142 /*** do ee: here, ij = ii+jj ***/
00143
00144 for( ij = 2 ; ij < (nx+ny-4) ; ij++ ){
00145
00146 jj = MIN( ij , ny-1 ) ; ii = ij - jj ;
00147
00148 for( nij=0 ; (ii<nx) && (jj>=0) ; nij++ , ii++ , jj-- ){
00149 flout[nij] = flee[ii+nx*jj] ;
00150 }
00151
00152 nonmax_kill( nloc , nij , flout ) ;
00153
00154 jj = MIN( ij , ny-1 ) ; ii = ij - jj ;
00155
00156 for( nij=0 ; (ii<nx) && (jj>=0) ; nij++ , ii++ , jj-- ){
00157 flee[ii+nx*jj] = flout[nij] ;
00158 }
00159 }
00160
00161 } /* end if( nloc > 0 ) */
00162
00163 /*.......................................................................*/
00164 /*** assign maximum of outputs at a given location to the final result ***/
00165
00166 for( jj=0 ; jj < ny ; jj++ ){ /* clear edges */
00167 joff = nx * jj ;
00168 flout[joff] = flout[joff+nx-1] = 0.0 ;
00169 }
00170 for( ii=0 ; ii < nx ; ii++ ){
00171 flout[ii] = flout[joff+ii] = 0.0 ;
00172 }
00173
00174 for( jj=1 ; jj < ny-1 ; jj++ ){
00175 joff = nx * jj ;
00176 for( ii=1 ; ii < nx-1 ; ii++ ){
00177 sxx = MAX( flxx[ii+joff] , flyy[ii+joff] ) ;
00178 sdd = MAX( fldd[ii+joff] , flee[ii+joff] ) ;
00179 flout[ii+joff] = MAX( sxx , sdd ) ;
00180 }
00181 }
00182
00183 /*.......................................................................*/
00184
00185 nloc = nloc_collective ;
00186 if( nloc > 0 ){
00187 int xx,yy , xbot,xtop , ybot,ytop , dx,dy ;
00188 float val ;
00189
00190 for( jj=1 ; jj < ny-1 ; jj++ ){
00191 joff = nx * jj ;
00192 ybot = MAX(0 ,jj-nloc) ;
00193 ytop = MIN(ny-1,jj+nloc) ;
00194
00195 for( ii=1 ; ii < nx-1 ; ii++ ){
00196 xbot = MAX(0 ,ii-nloc) ;
00197 xtop = MIN(nx-1,ii+nloc) ;
00198
00199 sxx = flxx[ii+joff] ; syy = flyy[ii+joff] ;
00200 sdd = fldd[ii+joff] ; see = flee[ii+joff] ; val = flout[ii+joff] ;
00201
00202 if( val == sxx ){ dx = 1 ; dy = 0 ; }
00203 else if( val == syy ){ dx = 0 ; dy = 1 ; }
00204 else if( val == sdd ){ dx = 1 ; dy = 1 ; }
00205 else { dx = 1 ; dy = -1 ; }
00206
00207 for( xx=ii+dx , yy=jj+dy ;
00208 (xx <= xtop) && (yy >= ybot) && (yy <= ytop) ;
00209 xx += dx , yy+= dy ) val = MAX(val,flout[xx+nx*yy]) ;
00210
00211 for( xx=ii-dx , yy=jj-dy ;
00212 (xx >= xbot) && (yy >= ybot) && (yy <= ytop) ;
00213 xx -= dx , yy-= dy ) val = MAX(val,flout[xx+nx*yy]) ;
00214
00215 if( flout[ii+joff] < val ) flout[ii+joff] = 0.0 ;
00216
00217 } /* end for ii */
00218 } /* end for jj */
00219
00220 /*** now, destroy points in flout with no neighbors ***/
00221
00222 for( ii=0 ; ii < nx*ny ; ii++ ) flee[ii] = 0.0 ;
00223
00224 for( jj=1 ; jj < ny-1 ; jj++ ){
00225 joff = nx * jj ;
00226 for( ii=1 ; ii < nx-1 ; ii++ ){
00227 xx = ii+joff ;
00228 if( flout[xx] == 0.0 ) continue ;
00229
00230 if( flout[xx-1] != 0.0 || flout[xx+1] != 0.0 ||
00231 flout[xx+nx] != 0.0 || flout[xx-nx] != 0.0 ||
00232 flout[xx+nx+1] != 0.0 || flout[xx+nx-1] != 0.0 ||
00233 flout[xx-nx+1] != 0.0 || flout[xx-nx-1] != 0.0 ) flee[xx] = flout[xx] ;
00234 }
00235 }
00236
00237 for( ii=0 ; ii < nx*ny ; ii++ ) flout[ii] = flee[ii] ;
00238
00239 } /* end if nloc */
00240
00241 /*........................................................................*/
00242
00243 mri_free(imxx) ; mri_free(imyy) ; mri_free(imdd) ; mri_free(imee) ;
00244
00245 MRI_COPY_AUX(imout,imin) ;
00246 return imout ;
00247 }
|
|
||||||||||||||||
|
Definition at line 251 of file mri_sobel.c. Referenced by mri_sobel().
00252 {
00253 int ii ;
00254 register int jj , jbot,jtop ;
00255 register float val ;
00256
00257 for( ii=0 ; ii < npt ; ii++ ){
00258
00259 jbot = MAX( 0 , ii-nloc ) ;
00260 jtop = MIN( npt , ii+nloc ) ;
00261 val = ar[jbot++] ;
00262 for( jj=jbot ; jj < jtop ; jj++ ) val = MAX(val,ar[jj]) ;
00263
00264 if( ar[ii] != val ) ar[ii] = 0.0 ;
00265 }
00266
00267 return ;
00268 }
|
