00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006    
00007 #include "mrilib.h"
00008 
00009 /*** NOT 7D SAFE ***/
00010 
00011 /*** prototype ***/
00012 
00013 void nonmax_kill( int , int , float * ) ;
00014 
00015 /*-------------------------------------------------------------------------
00016   Do Sobel edge enhancement on an image.  Output image will be floats.
00017   nloc = distance over which to suppress filter outputs which are not
00018          local maxima: _individual applies to each separate filter,
00019                        _collective applies to the merger of all 4 filters;
00020          (either or both may be zero).
00021 ---------------------------------------------------------------------------*/
00022 
00023 MRI_IMAGE * mri_sobel( int nloc_individual , int nloc_collective , MRI_IMAGE * imin )
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 }
00248 
00249 /*--------------------------------------------------------------------------*/
00250 
00251 void nonmax_kill( int nloc , int npt , float * ar )
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 }
00269 
00270 /*----------------------------------------------------------------
00271   Sharpen an image;  method same as "pgmenhance".
00272   phi is a parameter between 0.0 and 1.0: larger --> more sharp.
00273   logify is a flag: nonzero --> take log, sharpen, then exp
00274    (e.g., homomorphic filtering).
00275 ------------------------------------------------------------------*/
00276 
00277 MRI_IMAGE * mri_sharpen( float phi , int logify , MRI_IMAGE * im )
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 }