Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
imfit.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007 #include <string.h>
00008 #include "mrilib.h"
00009
00010 #define DFILT_SIGMA (4.0*0.42466090)
00011 #define MAX_ITER 5
00012 #define THRESH 0.10
00013
00014 int main( int argc , char * argv[] )
00015 {
00016 MRI_IMAGE * im1 , * im2 , *bim,*xim,*yim,*tim , *bim2 ;
00017 MRI_IMARR * fitim ;
00018 float * fit , *tar,*xar,*yar , *fit2 ;
00019 int nx,ny , ii,jj , joff , iter , good ;
00020 float hnx,hny,fac ;
00021
00022 if( argc < 3 || strncmp(argv[1],"-help",4) == 0 ){
00023 printf("Usage: imfit image1 image2 [output_file]\n"
00024 " Tries to find shift and rotation parameters to fit\n"
00025 " image2 to image1.\n"
00026 " N.B.: the method used is valid only for 1-2 pixel\n"
00027 " shifts and 1-2 degree rotations.\n"
00028 " If 'output_file' is given, image2 will be transformed\n"
00029 " accordingly and written into this file.\n"
00030 ) ;
00031 exit(0) ;
00032 }
00033
00034 machdep() ;
00035
00036 im1 = mri_read_just_one( argv[1] ) ;
00037 im2 = mri_read_just_one( argv[2] ) ;
00038 if( im1 == NULL || im2 == NULL ) exit(1) ;
00039
00040 if( im1->nx != im2->nx || im1->ny != im2->ny ){
00041 fprintf(stderr,"*** Input image mismatch: fatal error!\a\n");
00042 exit(1);
00043 }
00044 if( ! MRI_IS_2D(im1) || ! MRI_IS_2D(im2) ){
00045 fprintf(stderr,"*** Input images are not 2D!\a\n") ;
00046 exit(1) ;
00047 }
00048
00049 im1->dx = im1->dy = 1.0 ;
00050 nx = im1->nx ; hnx = 0.5 * nx ;
00051 ny = im1->ny ; hny = 0.5 * ny ;
00052 fac = (PI/180.0) ;
00053
00054 bim = mri_filt_fft( im1 , DFILT_SIGMA , 0 , 0 , 0 ) ;
00055 xim = mri_filt_fft( im1 , DFILT_SIGMA , 1 , 0 , 0 ) ;
00056 yim = mri_filt_fft( im1 , DFILT_SIGMA , 0 , 1 , 0 ) ;
00057
00058 tim = mri_new( nx , ny , MRI_float ) ;
00059 tar = mri_data_pointer( tim ) ;
00060 xar = mri_data_pointer( xim ) ;
00061 yar = mri_data_pointer( yim ) ;
00062 for( jj=0 ; jj < ny ; jj++ ){
00063 joff = jj * nx ;
00064 for( ii=0 ; ii < nx ; ii++ ){
00065 tar[ii+joff] = fac * ( (ii-hnx) * yar[ii+joff]
00066 - (jj-hny) * xar[ii+joff] ) ;
00067 }
00068 }
00069 INIT_IMARR ( fitim ) ;
00070 ADDTO_IMARR( fitim , bim ) ;
00071 ADDTO_IMARR( fitim , xim ) ;
00072 ADDTO_IMARR( fitim , yim ) ;
00073 ADDTO_IMARR( fitim , tim ) ;
00074
00075
00076 bim2 = mri_filt_fft( im2 , DFILT_SIGMA , 0 , 0 , 0 ) ;
00077 fit = mri_lsqfit( bim2 , fitim , bim2 ) ;
00078 mri_free( bim2 ) ;
00079
00080 #if 0
00081 printf("Command that transforms second image to first would be\n"
00082 " imrotate %g %g %g %s ...\n" ,
00083 fit[1] , fit[2] , fit[3] , argv[2] ) ;
00084 #endif
00085
00086 iter = 0 ;
00087 do {
00088
00089 tim = mri_rota( im2 , fit[1] , fit[2] , fit[3]*(PI/180.0) ) ;
00090 bim2 = mri_filt_fft( tim , DFILT_SIGMA , 0 , 0 , 0 ) ;
00091 fit2 = mri_lsqfit( bim2 , fitim , bim2 ) ;
00092 mri_free( bim2 ) ; mri_free( tim ) ;
00093
00094 fit[1] += fit2[1] ;
00095 fit[2] += fit2[2] ;
00096 fit[3] += fit2[3] ;
00097
00098 #if 0
00099 printf("Command that transforms second image to first would be\n"
00100 " imrotate %g %g %g %s ...\n" ,
00101 fit[1] , fit[2] , fit[3] , argv[2] ) ;
00102 #endif
00103
00104 iter++ ;
00105 good = (iter < MAX_ITER) &&
00106 ( fabs(fit2[1]) > THRESH ||
00107 fabs(fit2[2]) > THRESH || fabs(fit2[3]) > THRESH ) ;
00108 } while(good) ;
00109
00110 if( argc < 4 ){
00111 printf(" Command that transforms second image to first would be\n"
00112 " imrotate %g %g %g %s ...\n" ,
00113 fit[1] , fit[2] , fit[3] , argv[2] ) ;
00114 } else {
00115 tim = mri_rota( im2 , fit[1] , fit[2] , fit[3]*(PI/180.0) ) ;
00116 switch( im2->kind ){
00117 default: bim2 = tim ; break ;
00118 case MRI_short: bim2 = mri_to_short( 1.0 , tim ) ; break ;
00119 case MRI_byte: bim2 = mri_to_byte( tim ) ;
00120 }
00121 mri_write( argv[3] , bim2 ) ;
00122 printf(" imrotate %g %g %g %s %s\n",
00123 fit[1] , fit[2] , fit[3] , argv[2] , argv[3] ) ;
00124 }
00125
00126 exit(0) ;
00127 }