Doxygen Source Code Documentation

Main Page Alphabetical List Data Structures File List Data Fields Globals Search

mri_align.c File Reference

#include "mrilib.h"
#include "pcor.h"


Defines
#define	REF_FLOAT_SINGLE
#define	DFILT_SIGMA (4.0*0.42466090)
#define	MAX_ITER 5
#define	DXY_THRESH 0.15
#define	PHI_THRESH 0.45
#define	DFAC (PI/180.0)
#define	FINE_FIT
#define	FINE_SIGMA (1.0*0.42466090)
#define	FINE_DXY_THRESH 0.07
#define	FINE_PHI_THRESH 0.21
#define	USE_DELAYED_FIT
#define	MRI_ROTA_COARSE(a, b, c, d) mri_rota_variable(almode_coarse,(a),(b),(c),(d))
#define	MRI_ROTA_FINE(a, b, c, d) mri_rota_variable(almode_fine ,(a),(b),(c),(d))
#define	MRI_ROTA_REG(a, b, c, d) mri_rota_variable(almode_reg ,(a),(b),(c),(d))
Functions
void	mri_align_params (int maxite, float sig, float dxy, float dph, float fsig, float fdxy, float fdph)
void	mri_align_method (int coarse, int fine, int reg)
MRI_IMARR *	mri_align_dfspace (MRI_IMAGE imbase, MRI_IMAGE imwt, MRI_IMARR *ims, int code, float dx[], float dy[], float phi[])
Variables
float	dfilt_sigma = DFILT_SIGMA
float	dxy_thresh = DXY_THRESH
float	phi_thresh = PHI_THRESH
float	fine_sigma = FINE_SIGMA
float	fine_dxy_thresh = FINE_DXY_THRESH
float	fine_phi_thresh = FINE_PHI_THRESH
int	max_iter = MAX_ITER
int	almode_coarse = MRI_BICUBIC
int	almode_fine = MRI_BICUBIC
int	almode_reg = MRI_BICUBIC

Define Documentation

#define DFAC (PI/180.0)

Definition at line 19 of file mri_align.c.

Referenced by mri_align_dfspace().

#define DFILT_SIGMA (4.0*0.42466090)

Definition at line 15 of file mri_align.c.

Referenced by mri_align_params().

#define DXY_THRESH 0.15

Definition at line 17 of file mri_align.c.

Referenced by mri_align_params().

#define FINE_DXY_THRESH 0.07

Definition at line 30 of file mri_align.c.

Referenced by mri_align_params().

#define FINE_FIT

FINE_FIT: experimental code for a less smoothed fit at the end RW Cox, 17 July 1995

Definition at line 26 of file mri_align.c.

#define FINE_PHI_THRESH 0.21

Definition at line 31 of file mri_align.c.

Referenced by mri_align_params().

#define FINE_SIGMA (1.0*0.42466090)

Definition at line 29 of file mri_align.c.

#define MAX_ITER 5

Definition at line 16 of file mri_align.c.

Referenced by mri_align_params().

#define MRI_ROTA_COARSE	(	a,
		b,
		c,
		d	)	mri_rota_variable(almode_coarse,(a),(b),(c),(d))

Definition at line 71 of file mri_align.c.

#define MRI_ROTA_FINE	(	a,
		b,
		c,
		d	)	mri_rota_variable(almode_fine ,(a),(b),(c),(d))

Definition at line 72 of file mri_align.c.

#define MRI_ROTA_REG	(	a,
		b,
		c,
		d	)	mri_rota_variable(almode_reg ,(a),(b),(c),(d))

Definition at line 73 of file mri_align.c.

#define PHI_THRESH 0.45

Definition at line 18 of file mri_align.c.

Referenced by mri_align_params().

#define REF_FLOAT_SINGLE

Definition at line 11 of file mri_align.c.

#define USE_DELAYED_FIT

Definition at line 34 of file mri_align.c.

Function Documentation

MRI_IMARR* mri_align_dfspace	(	MRI_IMAGE *	imbase,
		MRI_IMAGE *	imwt,
		MRI_IMARR *	ims,
		int	code,
		float	dx[],
		float	dy[],
		float	phi[]
	)

Definition at line 107 of file mri_align.c.

References ADDTO_IMARR, ALIGN_DEBUG_CODE, ALIGN_NOITER_CODE, ALIGN_REGISTER_CODE, ALIGN_VERBOSE_CODE, DESTROY_IMARR, DFAC, dfilt_sigma, dxy_thresh, EXIT, FILT_FFT_WRAPAROUND, fine_dxy_thresh, fine_phi_thresh, fine_sigma, fit, free, MRI_IMARR::imarr, imbase, INIT_IMARR, max_iter, mri_data_pointer(), mri_delayed_lsqfit(), mri_filt_fft(), MRI_FLOAT_PTR, mri_free(), MRI_IS_2D, mri_lsqfit(), mri_new(), MRI_ROTA_COARSE, MRI_ROTA_FINE, MRI_ROTA_REG, mri_startup_lsqfit(), mri_to_float(), MRI_IMAGE::name, MRI_IMARR::num, MRI_IMAGE::nx, MRI_IMAGE::ny, and phi_thresh.

00110 {
00111    MRI_IMAGE * im1 , *bim,*xim,*yim,*tim , *bim2 , * im2 , *imww ;
00112    MRI_IMARR * fitim ;
00113    float * fit , *tar,*xar,*yar , *dfit ;
00114    int nx,ny, ii,jj, joff, iter, good, kim, debug, verbose ;
00115    float hnx,hny ;
00116    double * chol_fitim=NULL ;
00117 
00118 #ifdef FINE_FIT
00119    MRI_IMARR * fine_fitim  =NULL ;
00120    MRI_IMAGE * fine_imww   =NULL ;
00121    double * chol_fine_fitim=NULL ;
00122    int use_fine_fit = (fine_sigma > 0.0) ;
00123 #endif
00124 
00125    if( ! MRI_IS_2D(imbase) ){
00126       fprintf(stderr,"\n*** mri_align_dfspace: cannot use nD images!\a\n") ;
00127       EXIT(1) ;
00128    }
00129 
00130    debug    = (code & ALIGN_DEBUG_CODE)    != 0 ;
00131    verbose  = (code & ALIGN_VERBOSE_CODE)  != 0 && !debug ;
00132 
00133    if( verbose ){printf("-- mri_align_dfspace");fflush(stdout);}
00134 
00135    /** create the fitting images **/
00136 
00137    if( debug ){
00138       printf("-- mri_align_dfspace: code=%d\n",code);
00139       if( imbase->name != NULL )
00140         printf("-- imbase name = %s\n",imbase->name);
00141    }
00142 
00143    im1 = mri_to_float( imbase ) ;
00144    nx  = im1->nx ;  hnx = 0.5 * nx ;
00145    ny  = im1->ny ;  hny = 0.5 * ny ;
00146 
00147    bim = mri_filt_fft( im1 , dfilt_sigma , 0 , 0 , FILT_FFT_WRAPAROUND ) ;  /* smooth */
00148    xim = mri_filt_fft( im1 , dfilt_sigma , 1 , 0 , FILT_FFT_WRAPAROUND ) ;  /* d/dx */
00149    yim = mri_filt_fft( im1 , dfilt_sigma , 0 , 1 , FILT_FFT_WRAPAROUND ) ;  /* d/dy */
00150 
00151    tim = mri_new( nx , ny , MRI_float ) ;    /* x * d/dy - y * d/dx */
00152    tar = mri_data_pointer( tim ) ;           /* which is d/d(theta) */
00153    xar = mri_data_pointer( xim ) ;
00154    yar = mri_data_pointer( yim ) ;
00155    for( jj=0 ; jj < ny ; jj++ ){
00156       joff = jj * nx ;
00157       for( ii=0 ; ii < nx ; ii++ ){
00158          tar[ii+joff] = DFAC * (  (ii-hnx) * yar[ii+joff]
00159                                 - (jj-hny) * xar[ii+joff] ) ;
00160       }
00161    }
00162    INIT_IMARR ( fitim ) ;
00163    ADDTO_IMARR( fitim , bim ) ;
00164    ADDTO_IMARR( fitim , xim ) ;
00165    ADDTO_IMARR( fitim , yim ) ;
00166    ADDTO_IMARR( fitim , tim ) ;
00167 
00168    if( imwt == NULL ) imww = mri_to_float( bim ) ;  /* 28 Oct 1996 */
00169    else               imww = mri_to_float( imwt ) ;
00170 
00171    tar = MRI_FLOAT_PTR(imww) ;
00172    for( ii=0 ; ii < nx*ny ; ii++ ) tar[ii] = fabs(tar[ii]) ;  /* 16 Nov 1998 */
00173 
00174 #ifdef USE_DELAYED_FIT
00175    chol_fitim = mri_startup_lsqfit( fitim , imww ) ;
00176 #endif
00177 
00178    /*** create the FINE_FIT analog of the above, if required ***/
00179 #ifdef FINE_FIT
00180    if( use_fine_fit ){
00181      bim = mri_filt_fft( im1 , fine_sigma , 0 , 0 , FILT_FFT_WRAPAROUND ) ;  /* smooth */
00182      xim = mri_filt_fft( im1 , fine_sigma , 1 , 0 , FILT_FFT_WRAPAROUND ) ;  /* d/dx */
00183      yim = mri_filt_fft( im1 , fine_sigma , 0 , 1 , FILT_FFT_WRAPAROUND ) ;  /* d/dy */
00184 
00185      tim = mri_new( nx , ny , MRI_float ) ;    /* x * d/dy - y * d/dx */
00186      tar = mri_data_pointer( tim ) ;           /* which is d/d(theta) */
00187      xar = mri_data_pointer( xim ) ;
00188      yar = mri_data_pointer( yim ) ;
00189      for( jj=0 ; jj < ny ; jj++ ){
00190         joff = jj * nx ;
00191         for( ii=0 ; ii < nx ; ii++ ){
00192            tar[ii+joff] = DFAC * (  (ii-hnx) * yar[ii+joff]
00193                                   - (jj-hny) * xar[ii+joff] ) ;
00194         }
00195      }
00196      INIT_IMARR ( fine_fitim ) ;
00197      ADDTO_IMARR( fine_fitim , bim ) ;
00198      ADDTO_IMARR( fine_fitim , xim ) ;
00199      ADDTO_IMARR( fine_fitim , yim ) ;
00200      ADDTO_IMARR( fine_fitim , tim ) ;
00201 
00202      if( imwt == NULL ) fine_imww = mri_to_float( bim ) ;  /* 03 Oct 1997 */
00203      else               fine_imww = mri_to_float( imwt ) ;
00204 
00205      tar = MRI_FLOAT_PTR(fine_imww) ;
00206      for( ii=0 ; ii < nx*ny ; ii++ ) tar[ii] = fabs(tar[ii]) ;
00207 
00208 #ifdef USE_DELAYED_FIT
00209      chol_fine_fitim = mri_startup_lsqfit( fine_fitim , fine_imww ) ;
00210 #endif
00211    }
00212 #endif  /* FINE_FIT */
00213 
00214    mri_free( im1 ) ;
00215 
00216    /** fit each image to the fitting images **/
00217 
00218    for( kim=0 ; kim < ims->num ; kim++ ){
00219 
00220       if( verbose && kim%5 == 0 ){printf(".");fflush(stdout);}
00221 
00222       if( debug ) printf("-- Start image %d: %s\n",kim,
00223                          (ims->imarr[kim]->name==NULL)?" ":ims->imarr[kim]->name);
00224 
00225       im2  = mri_to_float( ims->imarr[kim] ) ;
00226       bim2 = mri_filt_fft( im2 , dfilt_sigma , 0 , 0 , FILT_FFT_WRAPAROUND ) ;
00227 #ifdef USE_DELAYED_FIT
00228       fit  = mri_delayed_lsqfit( bim2 , fitim , chol_fitim ) ;
00229 #else
00230       fit  = mri_lsqfit( bim2 , fitim , imww ) ;
00231 #endif
00232       mri_free( bim2 ) ;
00233 
00234       if( debug ) printf("   fit = %13.6g %13.6g %13.6g\n",
00235                          fit[1],fit[2],fit[3] ) ;
00236 
00237       iter = 0 ;
00238       good = ( fabs(fit[1]) > dxy_thresh ||
00239                fabs(fit[2]) > dxy_thresh || fabs(fit[3]) > phi_thresh ) &&
00240              ( (code & ALIGN_NOITER_CODE) == 0 ) ;
00241 
00242       while( good ){
00243          tim  = MRI_ROTA_COARSE( im2 , fit[1] , fit[2] , fit[3]*DFAC ) ;
00244          bim2 = mri_filt_fft( tim , dfilt_sigma , 0 , 0 , FILT_FFT_WRAPAROUND ) ;
00245 #ifdef USE_DELAYED_FIT
00246          dfit = mri_delayed_lsqfit( bim2 , fitim , chol_fitim ) ;
00247 #else
00248          dfit = mri_lsqfit( bim2 , fitim , imww ) ;
00249 #endif
00250          mri_free( bim2 ) ; mri_free( tim ) ;
00251 
00252          if( debug )
00253             printf(" Cdfit = %13.6g %13.6g %13.6g\n",
00254                    dfit[1],dfit[2],dfit[3] ) ;
00255 
00256          fit[1] += dfit[1] ;
00257          fit[2] += dfit[2] ;
00258          fit[3] += dfit[3] ;
00259 
00260          good = (++iter < max_iter) &&
00261                   ( fabs(dfit[1]) > dxy_thresh ||
00262                     fabs(dfit[2]) > dxy_thresh || fabs(dfit[3]) > phi_thresh ) ;
00263 
00264          free(dfit) ; dfit = NULL ;
00265       } /* end while */
00266 
00267       /*** perform fine adjustments (always use bicubic interpolation) ***/
00268 #ifdef FINE_FIT
00269       if( use_fine_fit && iter < max_iter && (code & ALIGN_NOITER_CODE) == 0 ){
00270          good = 1 ;
00271          while( good ){
00272             tim  = MRI_ROTA_FINE( im2 , fit[1] , fit[2] , fit[3]*DFAC ) ;
00273             bim2 = mri_filt_fft( tim , fine_sigma , 0 , 0 , FILT_FFT_WRAPAROUND ) ;
00274 #ifdef USE_DELAYED_FIT
00275             dfit = mri_delayed_lsqfit( bim2 , fine_fitim , chol_fine_fitim ) ;
00276 #else
00277             dfit = mri_lsqfit( bim2 , fine_fitim , fine_imww ) ;
00278 #endif
00279             mri_free( bim2 ) ; mri_free( tim ) ;
00280 
00281             if( debug )
00282                printf(" Fdfit = %13.6g %13.6g %13.6g\n",
00283                       dfit[1],dfit[2],dfit[3] ) ;
00284 
00285             fit[1] += dfit[1] ;
00286             fit[2] += dfit[2] ;
00287             fit[3] += dfit[3] ;
00288 
00289             good = (++iter < max_iter) &&
00290                      ( fabs(dfit[1]) > fine_dxy_thresh ||
00291                        fabs(dfit[2]) > fine_dxy_thresh || fabs(dfit[3]) > fine_phi_thresh ) ;
00292 
00293             free(dfit) ; dfit = NULL ;
00294          } /* end while */
00295       }
00296 #endif
00297 
00298       dx[kim]  = fit[1] ;
00299       dy[kim]  = fit[2] ;
00300       phi[kim] = fit[3]*DFAC ;
00301 
00302       if( debug )
00303          printf("   FIT = %13.6g %13.6g %13.6g\n",
00304                 fit[1],fit[2],fit[3] ) ;
00305 
00306       free(fit) ; fit = NULL ; mri_free( im2 ) ;
00307    }
00308 
00309    if( verbose ){printf("\n");fflush(stdout);}
00310 
00311    DESTROY_IMARR( fitim ) ;
00312 #ifdef USE_DELAYED_FIT
00313    free(chol_fitim) ; chol_fitim = NULL ;
00314 #endif
00315 
00316 #ifdef FINE_FIT
00317    if( use_fine_fit ){
00318      DESTROY_IMARR( fine_fitim ) ;
00319      mri_free( fine_imww ) ;
00320 #ifdef USE_DELAYED_FIT
00321      free(chol_fine_fitim) ; chol_fine_fitim = NULL ;
00322 #endif
00323    }
00324 #endif
00325 
00326    mri_free( imww ) ;
00327 
00328    if( (code & ALIGN_REGISTER_CODE) == 0 ) return NULL ;
00329 
00330    /** do the actual registration, if ordered (is always bicubic) **/
00331 
00332    INIT_IMARR( fitim ) ;
00333 
00334    if( verbose ){printf("-- registering");fflush(stdout);}
00335 
00336    for( kim=0 ; kim < ims->num ; kim++ ){
00337       tim = MRI_ROTA_REG( ims->imarr[kim] , dx[kim],dy[kim],phi[kim] ) ;
00338       ADDTO_IMARR( fitim , tim ) ;
00339       if( verbose && kim%5 == 0 ){printf(".");fflush(stdout);}
00340    }
00341 
00342    if( verbose ){printf("\n");fflush(stdout);}
00343 
00344    return fitim ;
00345 }

void mri_align_method	(	int	coarse,
		int	fine,
		int	reg
	)

Definition at line 75 of file mri_align.c.

References almode_coarse, almode_fine, and almode_reg.

Referenced by REG_command_line().

00076 {
00077    if( coarse > 0 ) almode_coarse = coarse ;
00078    if( fine   > 0 ) almode_fine   = fine   ;
00079    if( reg    > 0 ) almode_reg    = reg    ;
00080    return ;
00081 }

void mri_align_params	(	int	maxite,
		float	sig,
		float	dxy,
		float	dph,
		float	fsig,
		float	fdxy,
		float	fdph
	)

Definition at line 50 of file mri_align.c.

References dfilt_sigma, DFILT_SIGMA, dxy_thresh, DXY_THRESH, fine_dxy_thresh, FINE_DXY_THRESH, fine_phi_thresh, FINE_PHI_THRESH, fine_sigma, max_iter, MAX_ITER, phi_thresh, and PHI_THRESH.

Referenced by IMREG_main(), and REG_command_line().

00053 {
00054    if( maxite > 0   ) max_iter    = maxite ; else max_iter    = MAX_ITER    ;
00055    if( sig    > 0.0 ) dfilt_sigma = sig    ; else dfilt_sigma = DFILT_SIGMA ;
00056    if( dxy    > 0.0 ) dxy_thresh  = dxy    ; else dxy_thresh  = DXY_THRESH  ;
00057    if( dph    > 0.0 ) phi_thresh  = dph    ; else phi_thresh  = PHI_THRESH  ;
00058 
00059    fine_sigma = fsig ;
00060    if( fdxy > 0.0 ) fine_dxy_thresh = fdxy ; else fine_dxy_thresh = FINE_DXY_THRESH ;
00061    if( fdph > 0.0 ) fine_phi_thresh = fdph ; else fine_phi_thresh = FINE_PHI_THRESH ;
00062 
00063    return ;
00064 }