Doxygen Source Code Documentation

Main Page Alphabetical List Data Structures File List Data Fields Globals Search
edt_blur.c

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 /*************************************************************************
00010   Routine to blur a 3D volume with a Gaussian, using FFTs.
00011 **************************************************************************/
00012 
00013 #define GET_AS_BIG(name,type,dim)                                       \
00014    do{ if( (dim) > name ## _size ){                                     \
00015           if( name != NULL ) free(name) ;                               \
00016           name = (type *) malloc( sizeof(type) * (dim) ) ;              \
00017           if( name == NULL ){                                           \
00018              fprintf(stderr,"\n*** cannot malloc EDIT workspace!\n") ;  \
00019              EXIT(1) ; }                                                \
00020           name ## _size = (dim) ; }                                     \
00021        break ; } while(1)
00022 
00023 
00024 void EDIT_blur_volume( int nx, int ny, int nz,
00025                        float dx, float dy, float dz,
00026                        int ftype , void * vfim , float sigma )
00027 {
00028   EDIT_blur_volume_3d (nx, ny, nz, dx, dy, dz, ftype, vfim,
00029                        sigma, sigma, sigma);
00030 }
00031 
00032 
00033 /*
00034   The following slightly modified version of EDIT_blur_volume allows
00035   independent specification of Gaussian filter widths along the three
00036   perpendicular axes.
00037   BDW  21 Feb 1997
00038 */
00039 
00040 void EDIT_blur_volume_3d( int nx, int ny, int nz,
00041                           float dx, float dy, float dz,
00042                           int ftype , void * vfim ,
00043                           float sigmax, float sigmay, float sigmaz )
00044 {
00045    int jj,kk , nxy , base,nby2 ;
00046    float  dk , aa , k , fac ;
00047    register int ii , nup ;
00048 
00049    static int cx_size  = 0 ;     /* workspaces */
00050    static int gg_size  = 0 ;
00051    static complex * cx = NULL ;
00052    static float   * gg = NULL ;
00053 
00054    byte *    bfim = NULL ;
00055    short *   sfim = NULL ;
00056    float *   ffim = NULL ;
00057    complex * cfim = NULL ;
00058 
00059    float fbot,ftop ;     /* 10 Jan 2003 */
00060    int nxyz ;
00061 
00062    /*** initialize ***/
00063 
00064 ENTRY("EDIT_blur_volume_3d") ;
00065 
00066    if( vfim == NULL ||
00067        sigmax <= 0.0 || sigmay <= 0.0 || sigmaz <= 0.0)  EXRETURN ;
00068 
00069    switch( ftype ){
00070       default: EXRETURN ;
00071       case MRI_short:   sfim = (short *)   vfim ; break ;
00072       case MRI_float:   ffim = (float *)   vfim ; break ;
00073       case MRI_byte:    bfim = (byte *)    vfim ; break ;
00074       case MRI_complex: cfim = (complex *) vfim ; break ;
00075    }
00076    nxy = nx * ny ; nxyz = nxy * nz ;
00077 
00078    /*** 10 Jan 2003: find bot and top of data input */
00079 
00080    switch( ftype ){
00081      default:
00082        fbot = ftop = 0.0 ;
00083      break ;
00084 
00085      case MRI_short:
00086        fbot = ftop = sfim[0] ;
00087        for( ii=1 ; ii < nxyz ; ii++ )
00088               if( sfim[ii] < fbot ) fbot = sfim[ii] ;
00089          else if( sfim[ii] > ftop ) ftop = sfim[ii] ;
00090      break ;
00091 
00092      case MRI_float:
00093        fbot = ftop = ffim[0] ;
00094        for( ii=1 ; ii < nxyz ; ii++ )
00095               if( ffim[ii] < fbot ) fbot = ffim[ii] ;
00096          else if( ffim[ii] > ftop ) ftop = ffim[ii] ;
00097      break ;
00098 
00099      case MRI_byte:
00100        fbot = ftop = bfim[0] ;
00101        for( ii=1 ; ii < nxyz ; ii++ )
00102               if( bfim[ii] < fbot ) fbot = bfim[ii] ;
00103          else if( bfim[ii] > ftop ) ftop = bfim[ii] ;
00104      break ;
00105    }
00106 
00107    /*** do x-direction ***/
00108 
00109 STATUS("start x FFTs") ;
00110 
00111    aa  = sigmax * sigmax * 0.5 ;
00112    nup = nx + (int)(3.0 * sigmax / dx) ;      /* min FFT length */
00113 #if 0
00114    ii  = 2 ; while( ii < nup ){ ii *= 2 ; }  /* next power of 2 larger */
00115 #else
00116    ii = csfft_nextup_one35(nup) ;
00117 #endif
00118    nup = ii ; nby2 = nup / 2 ;
00119 
00120    GET_AS_BIG(cx,complex,nup) ; GET_AS_BIG(gg,float,nup) ;
00121 
00122    dk    = (2.0*PI) / (nup * dx) ;
00123    fac   = 1.0 / nup ;
00124    gg[0] = fac ;
00125    for( ii=1 ; ii<=nby2 ; ii++ ){ k=ii*dk; gg[nup-ii]=gg[ii]=fac*exp(-aa*k*k); }
00126 
00127    /** July 20: double up on FFTs **/
00128    /** Feb  09: extend to other data types besides shorts;
00129                 doubling up does not apply to complex data! **/
00130 
00131    switch( ftype ){
00132       case MRI_short:{
00133          register short * qfim ;
00134          for( kk=0 ; kk < nz ; kk++ ){
00135             for( jj=0 ; jj < ny ; jj+=2 ){
00136                base = jj*nx + kk*nxy ;
00137                qfim = sfim + base ;
00138                if( jj == ny-1 )
00139                   for( ii=0 ; ii<nx ; ii++){ cx[ii].r = qfim[ii] ; cx[ii].i = 0.0 ; }
00140                else
00141                   for( ii=0 ; ii<nx ; ii++){ cx[ii].r = qfim[ii] ; cx[ii].i = qfim[ii+nx] ; }
00142                for( ii=nx; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
00143                csfft_cox( -1 , nup , cx ) ;
00144                for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
00145                csfft_cox(  1 , nup , cx ) ;
00146                if( jj == ny-1 )
00147                   for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii].r ; }
00148                else
00149                   for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii].r ; qfim[ii+nx] = cx[ii].i ; }
00150             }
00151          }
00152       }
00153       break ;
00154 
00155       case MRI_float:{
00156          register float * qfim ;
00157          for( kk=0 ; kk < nz ; kk++ ){
00158             for( jj=0 ; jj < ny ; jj+=2 ){
00159                base = jj*nx + kk*nxy ;
00160                qfim = ffim + base ;
00161                if( jj == ny-1 )
00162                   for( ii=0 ; ii<nx ; ii++){ cx[ii].r = qfim[ii] ; cx[ii].i = 0.0 ; }
00163                else
00164                   for( ii=0 ; ii<nx ; ii++){ cx[ii].r = qfim[ii] ; cx[ii].i = qfim[ii+nx] ; }
00165                for( ii=nx; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
00166                csfft_cox( -1 , nup , cx ) ;
00167                for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
00168                csfft_cox(  1 , nup , cx ) ;
00169                if( jj == ny-1 )
00170                   for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii].r ; }
00171                else
00172                   for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii].r ; qfim[ii+nx] = cx[ii].i ; }
00173             }
00174          }
00175       }
00176       break ;
00177 
00178       case MRI_byte:{
00179          register byte * qfim ;
00180          for( kk=0 ; kk < nz ; kk++ ){
00181             for( jj=0 ; jj < ny ; jj+=2 ){
00182                base = jj*nx + kk*nxy ;
00183                qfim = bfim + base ;
00184                if( jj == ny-1 )
00185                   for( ii=0 ; ii<nx ; ii++){ cx[ii].r = qfim[ii] ; cx[ii].i = 0.0 ; }
00186                else
00187                   for( ii=0 ; ii<nx ; ii++){ cx[ii].r = qfim[ii] ; cx[ii].i = qfim[ii+nx] ; }
00188                for( ii=nx; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
00189                csfft_cox( -1 , nup , cx ) ;
00190                for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
00191                csfft_cox(  1 , nup , cx ) ;
00192                if( jj == ny-1 )
00193                   for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii].r ; }
00194                else
00195                   for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii].r ; qfim[ii+nx] = cx[ii].i ; }
00196             }
00197          }
00198       }
00199       break ;
00200 
00201       case MRI_complex:{
00202          register complex * qfim ;
00203          for( kk=0 ; kk < nz ; kk++ ){
00204             for( jj=0 ; jj < ny ; jj++ ){
00205                base = jj*nx + kk*nxy ;
00206                qfim = cfim + base ;
00207                for( ii=0 ; ii<nx ; ii++) { cx[ii] = qfim[ii] ; }
00208                for( ii=nx; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
00209                csfft_cox( -1 , nup , cx ) ;
00210                for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
00211                csfft_cox(  1 , nup , cx ) ;
00212                for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii] ; }
00213             }
00214          }
00215       }
00216       break ;
00217    }
00218 
00219    /*** do y-direction ***/
00220 
00221 STATUS("start y FFTs") ;
00222 
00223    aa  = sigmay * sigmay * 0.5 ;
00224    nup = ny + (int)(3.0 * sigmay / dy) ;      /* min FFT length */
00225 #if 0
00226    ii  = 2 ; while( ii < nup ){ ii *= 2 ; }  /* next power of 2 larger */
00227 #else
00228    ii = csfft_nextup_one35(nup) ;
00229 #endif
00230    nup = ii ; nby2 = nup / 2 ;
00231 
00232    GET_AS_BIG(cx,complex,nup) ; GET_AS_BIG(gg,float,nup) ;
00233 
00234    dk    = (2.0*PI) / (nup * dy) ;
00235    fac   = 1.0 / nup ;
00236    gg[0] = fac ;
00237    for( ii=1 ; ii<=nby2 ; ii++ ){ k=ii*dk; gg[nup-ii]=gg[ii]=fac*exp(-aa*k*k); }
00238 
00239    switch( ftype ){
00240       case MRI_short:{
00241          register short * qfim ;
00242          for( kk=0 ; kk < nz ; kk++ ){
00243             for( jj=0 ; jj < nx ; jj+=2 ){
00244                base = jj + kk*nxy ;
00245                qfim = sfim + base ;
00246                if( jj == nx-1 )
00247                   for( ii=0 ; ii<ny ; ii++){ cx[ii].r = qfim[ii*nx] ; cx[ii].i = 0.0 ; }
00248                else
00249                   for( ii=0 ; ii<ny ; ii++){ cx[ii].r = qfim[ii*nx] ; cx[ii].i = qfim[ii*nx+1] ; }
00250                for( ii=ny; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
00251                csfft_cox( -1 , nup , cx ) ;
00252                for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
00253                csfft_cox(  1 , nup , cx ) ;
00254                if( jj == nx-1 )
00255                   for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii].r ; }
00256                else
00257                   for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii].r ; qfim[ii*nx+1] = cx[ii].i ; }
00258             }
00259          }
00260       }
00261       break ;
00262 
00263       case MRI_byte:{
00264          register byte * qfim ;
00265          for( kk=0 ; kk < nz ; kk++ ){
00266             for( jj=0 ; jj < nx ; jj+=2 ){
00267                base = jj + kk*nxy ;
00268                qfim = bfim + base ;
00269                if( jj == nx-1 )
00270                   for( ii=0 ; ii<ny ; ii++){ cx[ii].r = qfim[ii*nx] ; cx[ii].i = 0.0 ; }
00271                else
00272                   for( ii=0 ; ii<ny ; ii++){ cx[ii].r = qfim[ii*nx] ; cx[ii].i = qfim[ii*nx+1] ; }
00273                for( ii=ny; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
00274                csfft_cox( -1 , nup , cx ) ;
00275                for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
00276                csfft_cox(  1 , nup , cx ) ;
00277                if( jj == nx-1 )
00278                   for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii].r ; }
00279                else
00280                   for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii].r ; qfim[ii*nx+1] = cx[ii].i ; }
00281             }
00282          }
00283       }
00284       break ;
00285 
00286       case MRI_float:{
00287          register float * qfim ;
00288          for( kk=0 ; kk < nz ; kk++ ){
00289             for( jj=0 ; jj < nx ; jj+=2 ){
00290                base = jj + kk*nxy ;
00291                qfim = ffim + base ;
00292                if( jj == nx-1 )
00293                   for( ii=0 ; ii<ny ; ii++){ cx[ii].r = qfim[ii*nx] ; cx[ii].i = 0.0 ; }
00294                else
00295                   for( ii=0 ; ii<ny ; ii++){ cx[ii].r = qfim[ii*nx] ; cx[ii].i = qfim[ii*nx+1] ; }
00296                for( ii=ny; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
00297                csfft_cox( -1 , nup , cx ) ;
00298                for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
00299                csfft_cox(  1 , nup , cx ) ;
00300                if( jj == nx-1 )
00301                   for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii].r ; }
00302                else
00303                   for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii].r ; qfim[ii*nx+1] = cx[ii].i ; }
00304             }
00305          }
00306       }
00307       break ;
00308 
00309       case MRI_complex:{
00310          register complex * qfim ;
00311          for( kk=0 ; kk < nz ; kk++ ){
00312             for( jj=0 ; jj < nx ; jj++ ){
00313                base = jj + kk*nxy ;
00314                qfim = cfim + base ;
00315                for( ii=0 ; ii<ny ; ii++){ cx[ii] = qfim[ii*nx] ; }
00316                for( ii=ny; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
00317                csfft_cox( -1 , nup , cx ) ;
00318                for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
00319                csfft_cox(  1 , nup , cx ) ;
00320                for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii] ; }
00321             }
00322          }
00323       }
00324       break ;
00325    }
00326 
00327    /*** do z-direction ***/
00328 
00329 STATUS("start z FFTs") ;
00330 
00331    aa  = sigmaz * sigmaz * 0.5 ;
00332    nup = nz + (int)(3.0 * sigmaz / dz) ;      /* min FFT length */
00333 #if 0
00334    ii  = 2 ; while( ii < nup ){ ii *= 2 ; }  /* next power of 2 larger */
00335 #else
00336    ii = csfft_nextup_one35(nup) ;
00337 #endif
00338    nup = ii ; nby2 = nup / 2 ;
00339 
00340    GET_AS_BIG(cx,complex,nup) ; GET_AS_BIG(gg,float,nup) ;
00341 
00342    dk    = (2.0*PI) / (nup * dz) ;
00343    fac   = 1.0 / nup ;
00344    gg[0] = fac ;
00345    for( ii=1 ; ii<=nby2 ; ii++ ){ k=ii*dk; gg[nup-ii]=gg[ii]=fac*exp(-aa*k*k); }
00346 
00347    switch( ftype ){
00348       case MRI_short:{
00349          register short * qfim ;
00350          for( kk=0 ; kk < ny ; kk++ ){
00351             for( jj=0 ; jj < nx ; jj+=2 ){
00352                base = jj + kk*nx ;
00353                qfim = sfim + base ;
00354                if( jj == nx-1 )
00355                   for( ii=0 ; ii<nz ; ii++){ cx[ii].r = qfim[ii*nxy] ; cx[ii].i = 0.0 ; }
00356                else
00357                   for( ii=0 ; ii<nz ; ii++){ cx[ii].r = qfim[ii*nxy] ; cx[ii].i = qfim[ii*nxy+1] ; }
00358                for( ii=nz; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
00359                csfft_cox( -1 , nup , cx ) ;
00360                for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
00361                csfft_cox(  1 , nup , cx ) ;
00362                if( jj == nx-1 )
00363                   for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii].r ; }
00364                else
00365                   for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii].r ; qfim[ii*nxy+1] = cx[ii].i ; }
00366             }
00367          }
00368       }
00369       break ;
00370 
00371       case MRI_float:{
00372          register float * qfim ;
00373          for( kk=0 ; kk < ny ; kk++ ){
00374             for( jj=0 ; jj < nx ; jj+=2 ){
00375                base = jj + kk*nx ;
00376                qfim = ffim + base ;
00377                if( jj == nx-1 )
00378                   for( ii=0 ; ii<nz ; ii++){ cx[ii].r = qfim[ii*nxy] ; cx[ii].i = 0.0 ; }
00379                else
00380                   for( ii=0 ; ii<nz ; ii++){ cx[ii].r = qfim[ii*nxy] ; cx[ii].i = qfim[ii*nxy+1] ; }
00381                for( ii=nz; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
00382                csfft_cox( -1 , nup , cx ) ;
00383                for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
00384                csfft_cox(  1 , nup , cx ) ;
00385                if( jj == nx-1 )
00386                   for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii].r ; }
00387                else
00388                   for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii].r ; qfim[ii*nxy+1] = cx[ii].i ; }
00389             }
00390          }
00391       }
00392       break ;
00393 
00394       case MRI_byte:{
00395          register byte * qfim ;
00396          for( kk=0 ; kk < ny ; kk++ ){
00397             for( jj=0 ; jj < nx ; jj+=2 ){
00398                base = jj + kk*nx ;
00399                qfim = bfim + base ;
00400                if( jj == nx-1 )
00401                   for( ii=0 ; ii<nz ; ii++){ cx[ii].r = qfim[ii*nxy] ; cx[ii].i = 0.0 ; }
00402                else
00403                   for( ii=0 ; ii<nz ; ii++){ cx[ii].r = qfim[ii*nxy] ; cx[ii].i = qfim[ii*nxy+1] ; }
00404                for( ii=nz; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
00405                csfft_cox( -1 , nup , cx ) ;
00406                for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
00407                csfft_cox(  1 , nup , cx ) ;
00408                if( jj == nx-1 )
00409                   for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii].r ; }
00410                else
00411                   for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii].r ; qfim[ii*nxy+1] = cx[ii].i ; }
00412             }
00413          }
00414       }
00415       break ;
00416 
00417       case MRI_complex:{
00418          register complex * qfim ;
00419          for( kk=0 ; kk < ny ; kk++ ){
00420             for( jj=0 ; jj < nx ; jj++ ){
00421                base = jj + kk*nx ;
00422                qfim = cfim + base ;
00423                for( ii=0 ; ii<nz ; ii++){ cx[ii] = qfim[ii*nxy] ; }
00424                for( ii=nz; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
00425                csfft_cox( -1 , nup , cx ) ;
00426                for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
00427                csfft_cox(  1 , nup , cx ) ;
00428                for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii] ; }
00429             }
00430          }
00431       }
00432       break ;
00433    }
00434 
00435    /*** 10 Jan 2003: clip data to bot and top found above ***/
00436    /***              to minimize Gibbs ringing artifacts  ***/
00437 
00438    switch( ftype ){
00439 
00440      case MRI_short:
00441        for( ii=0 ; ii < nxyz ; ii++ )
00442               if( sfim[ii] < fbot ) sfim[ii] = fbot ;
00443          else if( sfim[ii] > ftop ) sfim[ii] = ftop ;
00444      break ;
00445 
00446      case MRI_float:
00447        for( ii=0 ; ii < nxyz ; ii++ )
00448               if( ffim[ii] < fbot ) ffim[ii] = fbot ;
00449          else if( ffim[ii] > ftop ) ffim[ii] = ftop ;
00450      break ;
00451 
00452      case MRI_byte:
00453        for( ii=0 ; ii < nxyz ; ii++ )
00454               if( bfim[ii] < fbot ) bfim[ii] = fbot ;
00455          else if( bfim[ii] > ftop ) bfim[ii] = ftop ;
00456      break ;
00457    }
00458 
00459    /*** done! ***/
00460 
00461    EXRETURN ;
00462 }
AFNI/NIfTI Server

Sections

Personal tools

Navigation

Doxygen Source Code Documentation

edt_blur.c
