/*****************************************************************************
Major portions of this software are copyrighted by the Medical College
of Wisconsin, 1994-2000, and are released under the Gnu General Public
License, Version 2. See the file README.Copyright for details.
******************************************************************************/
#include "mrilib.h"
/*** NOT 7D SAFE ***/
/*************************************************************************
** Adapted from mri_align.c, for the purpose of real-time registration **
*************************************************************************/
#define DFILT_SIGMA (4.0*0.42466090) /* 4.0 = FWHM */
#define MAX_ITER 5
#define DXY_THRESH 0.15 /* pixels */
#define PHI_THRESH 0.45 /* degrees */
#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
/*********************************************************************
05 Nov 1997: make the parameters that control the iterations
be alterable.
**********************************************************************/
static float dfilt_sigma = DFILT_SIGMA ,
dxy_thresh = DXY_THRESH ,
phi_thresh = PHI_THRESH ,
fine_sigma = FINE_SIGMA ,
fine_dxy_thresh = FINE_DXY_THRESH ,
fine_phi_thresh = FINE_PHI_THRESH ;
static int max_iter = MAX_ITER ;
void mri_2dalign_params( int maxite,
float sig , float dxy , float dph ,
float fsig, float fdxy, float fdph )
{
if( maxite > 0 ) max_iter = maxite ; else max_iter = MAX_ITER ;
if( sig > 0.0 ) dfilt_sigma = sig ; else dfilt_sigma = DFILT_SIGMA ;
if( dxy > 0.0 ) dxy_thresh = dxy ; else dxy_thresh = DXY_THRESH ;
if( dph > 0.0 ) phi_thresh = dph ; else phi_thresh = PHI_THRESH ;
fine_sigma = fsig ;
if( fdxy > 0.0 ) fine_dxy_thresh = fdxy ; else fine_dxy_thresh = FINE_DXY_THRESH ;
if( fdph > 0.0 ) fine_phi_thresh = fdph ; else fine_phi_thresh = FINE_PHI_THRESH ;
return ;
}
/*-------------------------------------------------------------------------------*/
static int almode_coarse = MRI_BICUBIC ; /* 1 Oct 1998 */
static int almode_fine = MRI_BICUBIC ;
static int almode_reg = MRI_BICUBIC ;
#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))
void mri_2dalign_method( int coarse , int fine , int reg ) /* 1 Oct 1998 */
{
if( coarse > 0 ) almode_coarse = coarse ;
if( fine > 0 ) almode_fine = fine ;
if( reg > 0 ) almode_reg = reg ;
return ;
}
/*-------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------
Inputs: imbase = base image for alignment
imwt = image of weight factors to align to
(if NULL, will generate one internally)
Output: pointer to a MRI_2dalign_basis struct, for later use.
The malloc-ed data in there can be freed using
routine mri_2dalign_cleanup.
----------------------------------------------------------------------*/
MRI_2dalign_basis * mri_2dalign_setup( MRI_IMAGE * imbase , MRI_IMAGE * imwt )
{
MRI_IMAGE * im1 , *bim,*xim,*yim,*tim , *bim2 , * im2 , *imww ;
float *tar,*xar,*yar ;
int nx,ny , ii,jj , joff ;
int use_fine_fit = (fine_sigma > 0.0) ;
float hnx,hny ;
MRI_IMARR * fitim =NULL;
double * chol_fitim=NULL ;
MRI_IMARR * fine_fitim =NULL ;
double * chol_fine_fitim=NULL ;
MRI_2dalign_basis * bout = NULL ;
if( ! MRI_IS_2D(imbase) ){
fprintf(stderr,"\n*** mri_2dalign_setup: cannot use nD images!\a\n") ;
return NULL ;
}
im1 = mri_to_float( imbase ) ;
nx = im1->nx ; hnx = 0.5 * nx ;
ny = im1->ny ; hny = 0.5 * ny ;
bim = mri_filt_fft( im1, dfilt_sigma, 0, 0, FILT_FFT_WRAPAROUND ) ; /* smooth */
xim = mri_filt_fft( im1, dfilt_sigma, 1, 0, FILT_FFT_WRAPAROUND ) ; /* d/dx */
yim = mri_filt_fft( im1, dfilt_sigma, 0, 1, FILT_FFT_WRAPAROUND ) ; /* d/dy */
tim = mri_new( nx , ny , MRI_float ) ; /* x * d/dy - y * d/dx */
tar = mri_data_pointer( tim ) ; /* which is d/d(theta) */
xar = mri_data_pointer( xim ) ;
yar = mri_data_pointer( yim ) ;
for( jj=0 ; jj < ny ; jj++ ){
joff = jj * nx ;
for( ii=0 ; ii < nx ; ii++ ){
tar[ii+joff] = DFAC * ( (ii-hnx) * yar[ii+joff]
- (jj-hny) * xar[ii+joff] ) ;
}
}
INIT_IMARR ( fitim ) ;
ADDTO_IMARR( fitim , bim ) ;
ADDTO_IMARR( fitim , xim ) ;
ADDTO_IMARR( fitim , yim ) ;
ADDTO_IMARR( fitim , tim ) ;
if( imwt == NULL ) imww = mri_to_float( bim ) ; /* 28 Oct 1996 */
else imww = mri_to_float( imwt ) ;
tar = MRI_FLOAT_PTR(imww) ;
for( ii=0 ; ii < nx*ny ; ii++ ) tar[ii] = fabs(tar[ii]) ; /* 16 Nov 1998 */
chol_fitim = mri_startup_lsqfit( fitim , imww ) ;
mri_free(imww) ;
if( use_fine_fit ){
bim = mri_filt_fft( im1, fine_sigma, 0, 0, FILT_FFT_WRAPAROUND ) ; /* smooth */
xim = mri_filt_fft( im1, fine_sigma, 1, 0, FILT_FFT_WRAPAROUND ) ; /* d/dx */
yim = mri_filt_fft( im1, fine_sigma, 0, 1, FILT_FFT_WRAPAROUND ) ; /* d/dy */
tim = mri_new( nx , ny , MRI_float ) ; /* x * d/dy - y * d/dx */
tar = mri_data_pointer( tim ) ; /* which is d/d(theta) */
xar = mri_data_pointer( xim ) ;
yar = mri_data_pointer( yim ) ;
for( jj=0 ; jj < ny ; jj++ ){
joff = jj * nx ;
for( ii=0 ; ii < nx ; ii++ ){
tar[ii+joff] = DFAC * ( (ii-hnx) * yar[ii+joff]
- (jj-hny) * xar[ii+joff] ) ;
}
}
INIT_IMARR ( fine_fitim ) ;
ADDTO_IMARR( fine_fitim , bim ) ;
ADDTO_IMARR( fine_fitim , xim ) ;
ADDTO_IMARR( fine_fitim , yim ) ;
ADDTO_IMARR( fine_fitim , tim ) ;
if( imwt == NULL ) imww = mri_to_float( bim ) ; /* 03 Oct 1997 */
else imww = mri_to_float( imwt ) ;
tar = MRI_FLOAT_PTR(imww) ;
for( ii=0 ; ii < nx*ny ; ii++ ) tar[ii] = fabs(tar[ii]) ;
chol_fine_fitim = mri_startup_lsqfit( fine_fitim , imww ) ;
mri_free(imww) ;
}
mri_free(im1) ;
bout = (MRI_2dalign_basis *) malloc( sizeof(MRI_2dalign_basis) ) ;
bout->fitim = fitim ;
bout->chol_fitim = chol_fitim ;
bout->fine_fitim = fine_fitim ;
bout->chol_fine_fitim = chol_fine_fitim ;
return bout ;
}
/*-----------------------------------------------------------------------
Input: basis = MRI_2dalign_basis * return from setup routine above.
im = MRI_IMAGE * to align to base image
Output: return value is aligned image;
*dx , *dy, and *phi are set to rotation parameters.
Note that the returned image is floats.
-------------------------------------------------------------------------*/
MRI_IMAGE * mri_2dalign_one( MRI_2dalign_basis * basis , MRI_IMAGE * im ,
float * dx , float * dy , float * phi )
{
MRI_IMARR * fitim ;
double * chol_fitim=NULL ;
MRI_IMARR * fine_fitim =NULL ;
double * chol_fine_fitim=NULL ;
float * fit , *dfit ;
int nx,ny , ii,jj , joff , iter , good ;
int use_fine_fit = (fine_sigma > 0.0) ;
MRI_IMAGE * im2 , * bim2 , * tim ;
nx = im->nx ; ny = im->ny ;
fitim = basis->fitim ;
chol_fitim = basis->chol_fitim ;
fine_fitim = basis->fine_fitim ;
chol_fine_fitim = basis->chol_fine_fitim ;
/*-- register the image: coarsely --*/
im2 = mri_to_float( im ) ;
bim2 = mri_filt_fft( im2, dfilt_sigma, 0,0, FILT_FFT_WRAPAROUND ) ;
fit = mri_delayed_lsqfit( bim2 , fitim , chol_fitim ) ;
mri_free( bim2 ) ;
iter = 0 ;
good = ( fabs(fit[1]) > dxy_thresh ||
fabs(fit[2]) > dxy_thresh || fabs(fit[3]) > phi_thresh ) ;
/*-- iterate coarse fit --*/
while( good ){
tim = MRI_ROTA_COARSE( im2 , fit[1] , fit[2] , fit[3]*DFAC ) ;
bim2 = mri_filt_fft( tim, dfilt_sigma, 0,0, FILT_FFT_WRAPAROUND ) ;
dfit = mri_delayed_lsqfit( bim2 , fitim , chol_fitim ) ;
mri_free( bim2 ) ; mri_free( tim ) ;
fit[1] += dfit[1] ;
fit[2] += dfit[2] ;
fit[3] += dfit[3] ;
good = (++iter < max_iter) &&
( fabs(dfit[1]) > dxy_thresh ||
fabs(dfit[2]) > dxy_thresh || fabs(dfit[3]) > phi_thresh ) ;
free(dfit) ; dfit = NULL ;
} /* end while */
/*-- perform fine adjustments --*/
if( use_fine_fit ){
good = 1 ;
while( good ){
tim = MRI_ROTA_FINE( im2 , fit[1] , fit[2] , fit[3]*DFAC ) ;
bim2 = mri_filt_fft( tim, fine_sigma, 0,0, FILT_FFT_WRAPAROUND ) ;
dfit = mri_delayed_lsqfit( bim2 , fine_fitim , chol_fine_fitim ) ;
mri_free( bim2 ) ; mri_free( tim ) ;
fit[1] += dfit[1] ;
fit[2] += dfit[2] ;
fit[3] += dfit[3] ;
good = (++iter < max_iter) &&
( fabs(dfit[1]) > fine_dxy_thresh ||
fabs(dfit[2]) > fine_dxy_thresh || fabs(dfit[3]) > fine_phi_thresh ) ;
free(dfit) ; dfit = NULL ;
} /* end while */
}
/*-- save final alignment parameters --*/
if( dx != NULL ) *dx = fit[1] ;
if( dy != NULL ) *dy = fit[2] ;
if( phi!= NULL ) *phi= fit[3]*DFAC ;
/*-- do the actual realignment --*/
tim = MRI_ROTA_REG( im2 , fit[1],fit[2],fit[3]*DFAC ) ;
mri_free( im2 ) ;
return tim ;
}
MRI_IMARR * mri_2dalign_many( MRI_IMAGE * im , MRI_IMAGE * imwt , MRI_IMARR * ims ,
float * dx , float * dy , float * phi )
{
int kim ;
MRI_IMAGE * tim ;
MRI_IMARR * alim ;
MRI_2dalign_basis * basis ;
basis = mri_2dalign_setup( im , imwt ) ;
if( basis == NULL ) return NULL ;
INIT_IMARR( alim ) ;
for( kim=0 ; kim < ims->num ; kim++ ){
tim = mri_2dalign_one( basis , ims->imarr[kim] , dx+kim , dy+kim , phi+kim ) ;
ADDTO_IMARR(alim,tim) ;
}
mri_2dalign_cleanup( basis ) ;
return alim ;
}
void mri_2dalign_cleanup( MRI_2dalign_basis * basis )
{
if( basis == NULL ) return ;
if( basis->fitim != NULL ){ DESTROY_IMARR( basis->fitim ) ; }
if( basis->chol_fitim != NULL ){ free(basis->chol_fitim) ; }
if( basis->fine_fitim != NULL ){ DESTROY_IMARR( basis->fine_fitim ) ; }
if( basis->chol_fine_fitim != NULL ){ free(basis->chol_fine_fitim) ; }
free(basis) ; return ;
}
syntax highlighted by Code2HTML, v. 0.9.1