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  ** Adapted from mri_align.c, for the purpose of real-time registration **
00013  *************************************************************************/
00014 
00015 #define DFILT_SIGMA  (4.0*0.42466090)  /* 4.0 = FWHM */
00016 #define MAX_ITER     5
00017 #define DXY_THRESH   0.15         /* pixels */
00018 #define PHI_THRESH   0.45         /* degrees */
00019 #define DFAC         (PI/180.0)
00020 
00021 #define FINE_FIT
00022 
00023 #define FINE_SIGMA   (1.0*0.42466090)
00024 #define FINE_DXY_THRESH  0.07
00025 #define FINE_PHI_THRESH  0.21
00026 
00027 /*********************************************************************
00028   05 Nov 1997: make the parameters that control the iterations
00029                be alterable.
00030 **********************************************************************/
00031 
00032 static float dfilt_sigma     = DFILT_SIGMA ,
00033              dxy_thresh      = DXY_THRESH ,
00034              phi_thresh      = PHI_THRESH ,
00035              fine_sigma      = FINE_SIGMA ,
00036              fine_dxy_thresh = FINE_DXY_THRESH ,
00037              fine_phi_thresh = FINE_PHI_THRESH  ;
00038 
00039 static int max_iter = MAX_ITER ;
00040 
00041 void mri_2dalign_params( int maxite,
00042                          float sig , float dxy , float dph ,
00043                          float fsig, float fdxy, float fdph )
00044 {
00045    if( maxite > 0   ) max_iter    = maxite ; else max_iter    = MAX_ITER    ;
00046    if( sig    > 0.0 ) dfilt_sigma = sig    ; else dfilt_sigma = DFILT_SIGMA ;
00047    if( dxy    > 0.0 ) dxy_thresh  = dxy    ; else dxy_thresh  = DXY_THRESH  ;
00048    if( dph    > 0.0 ) phi_thresh  = dph    ; else phi_thresh  = PHI_THRESH  ;
00049 
00050    fine_sigma = fsig ;
00051    if( fdxy > 0.0 ) fine_dxy_thresh = fdxy ; else fine_dxy_thresh = FINE_DXY_THRESH ;
00052    if( fdph > 0.0 ) fine_phi_thresh = fdph ; else fine_phi_thresh = FINE_PHI_THRESH ;
00053 
00054    return ;
00055 }
00056 
00057 /*-------------------------------------------------------------------------------*/
00058 static int almode_coarse = MRI_BICUBIC ;  /* 1 Oct 1998 */
00059 static int almode_fine   = MRI_BICUBIC ;
00060 static int almode_reg    = MRI_BICUBIC ;
00061 
00062 #define MRI_ROTA_COARSE(a,b,c,d) mri_rota_variable(almode_coarse,(a),(b),(c),(d))
00063 #define MRI_ROTA_FINE(a,b,c,d)   mri_rota_variable(almode_fine  ,(a),(b),(c),(d))
00064 #define MRI_ROTA_REG(a,b,c,d)    mri_rota_variable(almode_reg   ,(a),(b),(c),(d))
00065 
00066 void mri_2dalign_method( int coarse , int fine , int reg )  /* 1 Oct 1998 */
00067 {
00068    if( coarse > 0 ) almode_coarse = coarse ;
00069    if( fine   > 0 ) almode_fine   = fine   ;
00070    if( reg    > 0 ) almode_reg    = reg    ;
00071    return ;
00072 }
00073 /*-------------------------------------------------------------------------------*/
00074 
00075 /*--------------------------------------------------------------------
00076    Inputs: imbase = base image for alignment
00077            imwt   = image of weight factors to align to
00078                       (if NULL, will generate one internally)
00079 
00080    Output: pointer to a MRI_2dalign_basis struct, for later use.
00081            The malloc-ed data in there can be freed using
00082            routine mri_2dalign_cleanup.
00083 ----------------------------------------------------------------------*/
00084 
00085 MRI_2dalign_basis * mri_2dalign_setup( MRI_IMAGE * imbase , MRI_IMAGE * imwt )
00086 {
00087    MRI_IMAGE * im1 , *bim,*xim,*yim,*tim , *bim2 , * im2 , *imww ;
00088    float *tar,*xar,*yar ;
00089    int nx,ny , ii,jj , joff ;
00090    int use_fine_fit = (fine_sigma > 0.0) ;
00091    float hnx,hny ;
00092 
00093    MRI_IMARR * fitim  =NULL;
00094    double * chol_fitim=NULL ;
00095    MRI_IMARR * fine_fitim  =NULL ;
00096    double * chol_fine_fitim=NULL ;
00097    MRI_2dalign_basis * bout = NULL ;
00098 
00099    if( ! MRI_IS_2D(imbase) ){
00100       fprintf(stderr,"\n*** mri_2dalign_setup: cannot use nD images!\a\n") ;
00101       return NULL ;
00102    }
00103 
00104    im1 = mri_to_float( imbase ) ;
00105    nx  = im1->nx ;  hnx = 0.5 * nx ;
00106    ny  = im1->ny ;  hny = 0.5 * ny ;
00107 
00108    bim = mri_filt_fft( im1, dfilt_sigma, 0, 0, FILT_FFT_WRAPAROUND ) ; /* smooth */
00109    xim = mri_filt_fft( im1, dfilt_sigma, 1, 0, FILT_FFT_WRAPAROUND ) ; /* d/dx */
00110    yim = mri_filt_fft( im1, dfilt_sigma, 0, 1, FILT_FFT_WRAPAROUND ) ; /* d/dy */
00111 
00112    tim = mri_new( nx , ny , MRI_float ) ;  /* x * d/dy - y * d/dx */
00113    tar = mri_data_pointer( tim ) ;         /* which is d/d(theta) */
00114    xar = mri_data_pointer( xim ) ;
00115    yar = mri_data_pointer( yim ) ;
00116    for( jj=0 ; jj < ny ; jj++ ){
00117       joff = jj * nx ;
00118       for( ii=0 ; ii < nx ; ii++ ){
00119          tar[ii+joff] = DFAC * (  (ii-hnx) * yar[ii+joff]
00120                                 - (jj-hny) * xar[ii+joff] ) ;
00121       }
00122    }
00123    INIT_IMARR ( fitim ) ;
00124    ADDTO_IMARR( fitim , bim ) ;
00125    ADDTO_IMARR( fitim , xim ) ;
00126    ADDTO_IMARR( fitim , yim ) ;
00127    ADDTO_IMARR( fitim , tim ) ;
00128 
00129    if( imwt == NULL ) imww = mri_to_float( bim ) ;  /* 28 Oct 1996 */
00130    else               imww = mri_to_float( imwt ) ;
00131 
00132    tar = MRI_FLOAT_PTR(imww) ;       
00133    for( ii=0 ; ii < nx*ny ; ii++ ) tar[ii] = fabs(tar[ii]) ;  /* 16 Nov 1998 */
00134 
00135    chol_fitim = mri_startup_lsqfit( fitim , imww ) ;
00136    mri_free(imww) ;
00137 
00138    if( use_fine_fit ){
00139      bim = mri_filt_fft( im1, fine_sigma, 0, 0, FILT_FFT_WRAPAROUND ) ; /* smooth */
00140      xim = mri_filt_fft( im1, fine_sigma, 1, 0, FILT_FFT_WRAPAROUND ) ; /* d/dx */
00141      yim = mri_filt_fft( im1, fine_sigma, 0, 1, FILT_FFT_WRAPAROUND ) ; /* d/dy */
00142 
00143      tim = mri_new( nx , ny , MRI_float ) ;  /* x * d/dy - y * d/dx */
00144      tar = mri_data_pointer( tim ) ;         /* which is d/d(theta) */
00145      xar = mri_data_pointer( xim ) ;
00146      yar = mri_data_pointer( yim ) ;
00147      for( jj=0 ; jj < ny ; jj++ ){
00148         joff = jj * nx ;
00149         for( ii=0 ; ii < nx ; ii++ ){
00150            tar[ii+joff] = DFAC * (  (ii-hnx) * yar[ii+joff]
00151                                   - (jj-hny) * xar[ii+joff] ) ;
00152         }
00153      }
00154      INIT_IMARR ( fine_fitim ) ;
00155      ADDTO_IMARR( fine_fitim , bim ) ;
00156      ADDTO_IMARR( fine_fitim , xim ) ;
00157      ADDTO_IMARR( fine_fitim , yim ) ;
00158      ADDTO_IMARR( fine_fitim , tim ) ;
00159 
00160      if( imwt == NULL ) imww = mri_to_float( bim ) ;  /* 03 Oct 1997 */
00161      else               imww = mri_to_float( imwt ) ;
00162 
00163      tar = MRI_FLOAT_PTR(imww) ;       
00164      for( ii=0 ; ii < nx*ny ; ii++ ) tar[ii] = fabs(tar[ii]) ;
00165 
00166      chol_fine_fitim = mri_startup_lsqfit( fine_fitim , imww ) ;
00167      mri_free(imww) ;
00168    }
00169 
00170    mri_free(im1) ;
00171 
00172    bout = (MRI_2dalign_basis *) malloc( sizeof(MRI_2dalign_basis) ) ;
00173    bout->fitim           = fitim ;
00174    bout->chol_fitim      = chol_fitim ;
00175    bout->fine_fitim      = fine_fitim ;
00176    bout->chol_fine_fitim = chol_fine_fitim ;
00177    return bout ;
00178 }
00179 
00180 /*-----------------------------------------------------------------------
00181    Input:   basis  = MRI_2dalign_basis * return from setup routine above.
00182             im     = MRI_IMAGE * to align to base image
00183 
00184    Output:  return value is aligned image;
00185             *dx , *dy, and *phi are set to rotation parameters.
00186             Note that the returned image is floats.
00187 -------------------------------------------------------------------------*/
00188 
00189 MRI_IMAGE * mri_2dalign_one( MRI_2dalign_basis * basis , MRI_IMAGE * im ,
00190                              float * dx , float * dy , float * phi )
00191 {
00192    MRI_IMARR * fitim ;
00193    double * chol_fitim=NULL ;
00194    MRI_IMARR * fine_fitim  =NULL ;
00195    double * chol_fine_fitim=NULL ;
00196 
00197    float * fit , *dfit ;
00198    int nx,ny , ii,jj , joff , iter , good ;
00199    int use_fine_fit = (fine_sigma > 0.0) ;
00200    MRI_IMAGE * im2 , * bim2 , * tim ;
00201 
00202    nx  = im->nx ; ny  = im->ny ;
00203 
00204    fitim           = basis->fitim ;
00205    chol_fitim      = basis->chol_fitim ;
00206    fine_fitim      = basis->fine_fitim ;
00207    chol_fine_fitim = basis->chol_fine_fitim ;
00208 
00209    /*-- register the image: coarsely --*/
00210 
00211       im2  = mri_to_float( im ) ;
00212       bim2 = mri_filt_fft( im2, dfilt_sigma, 0,0, FILT_FFT_WRAPAROUND ) ;
00213       fit  = mri_delayed_lsqfit( bim2 , fitim , chol_fitim ) ;
00214       mri_free( bim2 ) ;
00215 
00216       iter = 0 ;
00217       good = ( fabs(fit[1]) > dxy_thresh ||
00218                fabs(fit[2]) > dxy_thresh || fabs(fit[3]) > phi_thresh ) ;
00219 
00220       /*-- iterate coarse fit --*/
00221 
00222       while( good ){
00223          tim  = MRI_ROTA_COARSE( im2 , fit[1] , fit[2] , fit[3]*DFAC ) ;
00224          bim2 = mri_filt_fft( tim, dfilt_sigma, 0,0, FILT_FFT_WRAPAROUND ) ;
00225          dfit = mri_delayed_lsqfit( bim2 , fitim , chol_fitim ) ;
00226          mri_free( bim2 ) ; mri_free( tim ) ;
00227 
00228          fit[1] += dfit[1] ;
00229          fit[2] += dfit[2] ;
00230          fit[3] += dfit[3] ;
00231 
00232          good = (++iter < max_iter) &&
00233                   ( fabs(dfit[1]) > dxy_thresh ||
00234                     fabs(dfit[2]) > dxy_thresh || fabs(dfit[3]) > phi_thresh ) ;
00235 
00236          free(dfit) ; dfit = NULL ;
00237       } /* end while */
00238 
00239       /*-- perform fine adjustments  --*/
00240 
00241       if( use_fine_fit ){
00242          good = 1 ;
00243          while( good ){
00244             tim  = MRI_ROTA_FINE( im2 , fit[1] , fit[2] , fit[3]*DFAC ) ;
00245             bim2 = mri_filt_fft( tim, fine_sigma, 0,0, FILT_FFT_WRAPAROUND ) ;
00246             dfit = mri_delayed_lsqfit( bim2 , fine_fitim , chol_fine_fitim ) ;
00247             mri_free( bim2 ) ; mri_free( tim ) ;
00248 
00249             fit[1] += dfit[1] ;
00250             fit[2] += dfit[2] ;
00251             fit[3] += dfit[3] ;
00252 
00253             good = (++iter < max_iter) &&
00254                      ( fabs(dfit[1]) > fine_dxy_thresh ||
00255                        fabs(dfit[2]) > fine_dxy_thresh || fabs(dfit[3]) > fine_phi_thresh ) ;
00256 
00257             free(dfit) ; dfit = NULL ;
00258          } /* end while */
00259       }
00260 
00261       /*-- save final alignment parameters --*/
00262 
00263       if( dx != NULL ) *dx = fit[1] ;
00264       if( dy != NULL ) *dy = fit[2] ;
00265       if( phi!= NULL ) *phi= fit[3]*DFAC ;
00266 
00267    /*-- do the actual realignment --*/
00268 
00269    tim = MRI_ROTA_REG( im2 , fit[1],fit[2],fit[3]*DFAC ) ;
00270    mri_free( im2 ) ;
00271    return tim ;
00272 }
00273 
00274 MRI_IMARR * mri_2dalign_many( MRI_IMAGE * im , MRI_IMAGE * imwt , MRI_IMARR * ims ,
00275                               float * dx , float * dy , float * phi )
00276 {
00277    int kim ;
00278    MRI_IMAGE * tim ;
00279    MRI_IMARR * alim ;
00280    MRI_2dalign_basis * basis ;
00281 
00282    basis = mri_2dalign_setup( im , imwt ) ;
00283    if( basis == NULL ) return NULL ;
00284 
00285    INIT_IMARR( alim ) ;
00286 
00287    for( kim=0 ; kim < ims->num ; kim++ ){
00288       tim = mri_2dalign_one( basis , ims->imarr[kim] , dx+kim , dy+kim , phi+kim ) ;
00289       ADDTO_IMARR(alim,tim) ;
00290    }
00291 
00292    mri_2dalign_cleanup( basis ) ;
00293    return alim ;
00294 }
00295 
00296 void mri_2dalign_cleanup( MRI_2dalign_basis * basis )
00297 {
00298    if( basis == NULL ) return ;
00299 
00300    if( basis->fitim           != NULL ){ DESTROY_IMARR( basis->fitim ) ; }
00301    if( basis->chol_fitim      != NULL ){ free(basis->chol_fitim) ; }
00302 
00303    if( basis->fine_fitim      != NULL ){ DESTROY_IMARR( basis->fine_fitim ) ; }
00304    if( basis->chol_fine_fitim != NULL ){ free(basis->chol_fine_fitim) ; }
00305 
00306    free(basis) ; return ;
00307 }