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 /**************************************************************************/
00012 /*
00013     mode = -1 for forward transform
00014          = +1 for inverse (including scaling)
00015 
00016    taper = fraction of data to taper (0 to 1)
00017 */
00018 
00019 void mri_fft_complex( int mode , float taper , MRI_IMAGE *im )
00020 {
00021    float *rbuf , *ibuf , *xtap , *ytap ;
00022    complex *cxim ;
00023    int ii , jj , npix , jbase , nx,ny ;
00024 
00025    if( im->kind != MRI_complex ){
00026       fprintf( stderr , "mri_fft_complex only works on complex images!\n" ) ;
00027       MRI_FATAL_ERROR ;
00028    }
00029 
00030    if( ! MRI_IS_2D(im) ){
00031       fprintf(stderr,"mri_fft_complex only works on 2D images!\n") ;
00032       MRI_FATAL_ERROR ;
00033    }
00034 
00035    /*** set up buffers ***/
00036 
00037    npix = im->nx * im->ny ;                           /* number of pixels */
00038    rbuf = (float *)malloc( sizeof(float) * npix ) ;   /* real and imag buffs */
00039    ibuf = (float *)malloc( sizeof(float) * npix ) ;
00040    cxim = mri_data_pointer( im ) ;                    /* easy acces to im */
00041 
00042    for( ii=0 ; ii < npix ; ii++ ){
00043       rbuf[ii] = cxim[ii].r ;
00044       ibuf[ii] = cxim[ii].i ;
00045    }
00046 
00047    /*** taper buffers, if desired ***/
00048 
00049    if( taper > 0.0 && taper <= 1.0 ){
00050       nx   = im->nx ;
00051       ny   = im->ny ;
00052       xtap = mri_setup_taper( nx , taper ) ;
00053 
00054 /***
00055       printf( "taper" ) ;
00056       for( ii=0 ; ii < nx ; ii++ ){
00057          if( (ii%5) == 0 ) printf("\n") ;
00058          printf( "%12.4e " , xtap[ii] ) ;
00059       }
00060       printf("\n") ;
00061 ***/
00062 
00063       if( nx == ny ) ytap = xtap ;
00064       else           ytap = mri_setup_taper( ny , taper ) ;
00065 
00066       for( jj=0 ; jj < ny ; jj++ ){
00067          jbase = jj * nx ;
00068          for( ii=0 ; ii < nx ; ii++ ){
00069             rbuf[ii] *= xtap[ii] * ytap[jj] ;
00070             ibuf[ii] *= xtap[ii] * ytap[jj] ;
00071          }
00072       }
00073       free( xtap ) ;
00074       if( ytap != xtap ) free(ytap) ;
00075    }
00076 
00077    /*** FFT buffers and copy them back to original image ***/
00078 
00079    cfft2d( mode , im->nx , im->ny , rbuf,ibuf ) ;
00080 
00081    for( ii=0 ; ii < npix ; ii++ ){
00082       cxim[ii].r = rbuf[ii] ;
00083       cxim[ii].i = ibuf[ii] ;
00084    }
00085 
00086    return ;
00087 }
00088 
00089 /***********************************************************************/
00090 
00091 float *mri_setup_taper( int nx , float taper )
00092 {
00093    register int ii ;
00094    int ntap ;
00095    float *tap ;
00096    float phi ; 
00097 
00098    tap = (float *)malloc( sizeof(float) * nx ) ;   /* make array */
00099 
00100    for( ii=0 ; ii < nx ; ii++ ) tap[ii] = 1.0 ;    /* default 1's */
00101 
00102    ntap = (int) (nx * 0.5 * taper ) ;              /* # pts on each end */
00103 
00104    if( ntap == 0 ){             /* special case of no points: taper is tiny */
00105       tap[0] = tap[nx-1] = 0.5 ;
00106       return tap ;
00107    }
00108 
00109    phi = PI / ntap ;
00110    for( ii=0 ; ii < ntap ; ii++ ){
00111       tap[ii]      = 0.5 - 0.5 * cos( ii*phi ) ;
00112       tap[nx-1-ii] = tap[ii] ;
00113    }
00114 
00115    return tap ;
00116 }