1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include <string.h>
8 #include "mrilib.h"
9
10 #define DFILT_SIGMA (4.0*0.42466090) /* 4.0 = FWHM */
11 #define MAX_ITER 5
12 #define THRESH 0.10
13
main(int argc,char * argv[])14 int main( int argc , char * argv[] )
15 {
16 MRI_IMAGE * im1 , * im2 , *bim,*xim,*yim,*tim , *bim2 ;
17 MRI_IMARR * fitim ;
18 float * fit , *tar,*xar,*yar , *fit2 ;
19 int nx,ny , ii,jj , joff , iter , good ;
20 float hnx,hny,fac ;
21
22 if( argc < 3 || strncmp(argv[1],"-help",4) == 0 ){
23 printf("Usage: imfit image1 image2 [output_file]\n"
24 " Tries to find shift and rotation parameters to fit\n"
25 " image2 to image1.\n"
26 " N.B.: the method used is valid only for 1-2 pixel\n"
27 " shifts and 1-2 degree rotations.\n"
28 " If 'output_file' is given, image2 will be transformed\n"
29 " accordingly and written into this file.\n"
30 ) ;
31 exit(0) ;
32 }
33
34 machdep() ;
35
36 im1 = mri_read_just_one( argv[1] ) ;
37 im2 = mri_read_just_one( argv[2] ) ;
38 if( im1 == NULL || im2 == NULL ) exit(1) ;
39
40 if( im1->nx != im2->nx || im1->ny != im2->ny ){
41 fprintf(stderr,"*** Input image mismatch: fatal error!\a\n");
42 exit(1);
43 }
44 if( ! MRI_IS_2D(im1) || ! MRI_IS_2D(im2) ){
45 fprintf(stderr,"*** Input images are not 2D!\a\n") ;
46 exit(1) ;
47 }
48
49 im1->dx = im1->dy = 1.0 ;
50 nx = im1->nx ; hnx = 0.5 * nx ;
51 ny = im1->ny ; hny = 0.5 * ny ;
52 fac = (PI/180.0) ;
53
54 bim = mri_filt_fft( im1 , DFILT_SIGMA , 0 , 0 , 0 ) ; /* smooth */
55 xim = mri_filt_fft( im1 , DFILT_SIGMA , 1 , 0 , 0 ) ; /* smooth and d/dx */
56 yim = mri_filt_fft( im1 , DFILT_SIGMA , 0 , 1 , 0 ) ; /* smooth and d/dy */
57
58 tim = mri_new( nx , ny , MRI_float ) ; /* x * d/dy - y * d/dx */
59 tar = mri_data_pointer( tim ) ; /* which is d/d(theta) */
60 xar = mri_data_pointer( xim ) ;
61 yar = mri_data_pointer( yim ) ;
62 for( jj=0 ; jj < ny ; jj++ ){
63 joff = jj * nx ;
64 for( ii=0 ; ii < nx ; ii++ ){
65 tar[ii+joff] = fac * ( (ii-hnx) * yar[ii+joff]
66 - (jj-hny) * xar[ii+joff] ) ;
67 }
68 }
69 INIT_IMARR ( fitim ) ;
70 ADDTO_IMARR( fitim , bim ) ;
71 ADDTO_IMARR( fitim , xim ) ;
72 ADDTO_IMARR( fitim , yim ) ;
73 ADDTO_IMARR( fitim , tim ) ;
74
75
76 bim2 = mri_filt_fft( im2 , DFILT_SIGMA , 0 , 0 , 0 ) ; /* smooth input image */
77 fit = mri_lsqfit( bim2 , fitim , bim2 ) ; /* fit input image to ref images */
78 mri_free( bim2 ) ;
79
80 #if 0
81 printf("Command that transforms second image to first would be\n"
82 " imrotate %g %g %g %s ...\n" ,
83 fit[1] , fit[2] , fit[3] , argv[2] ) ;
84 #endif
85
86 iter = 0 ;
87 do {
88
89 tim = mri_rota( im2 , fit[1] , fit[2] , fit[3]*(PI/180.0) ) ;
90 bim2 = mri_filt_fft( tim , DFILT_SIGMA , 0 , 0 , 0 ) ; /* smooth input image */
91 fit2 = mri_lsqfit( bim2 , fitim , bim2 ) ; /* fit input image to ref images */
92 mri_free( bim2 ) ; mri_free( tim ) ;
93
94 fit[1] += fit2[1] ;
95 fit[2] += fit2[2] ;
96 fit[3] += fit2[3] ;
97
98 #if 0
99 printf("Command that transforms second image to first would be\n"
100 " imrotate %g %g %g %s ...\n" ,
101 fit[1] , fit[2] , fit[3] , argv[2] ) ;
102 #endif
103
104 iter++ ;
105 good = (iter < MAX_ITER) &&
106 ( fabs(fit2[1]) > THRESH ||
107 fabs(fit2[2]) > THRESH || fabs(fit2[3]) > THRESH ) ;
108 } while(good) ;
109
110 if( argc < 4 ){
111 printf(" Command that transforms second image to first would be\n"
112 " imrotate %g %g %g %s ...\n" ,
113 fit[1] , fit[2] , fit[3] , argv[2] ) ;
114 } else {
115 tim = mri_rota( im2 , fit[1] , fit[2] , fit[3]*(PI/180.0) ) ;
116 switch( im2->kind ){
117 default: bim2 = tim ; break ;
118 case MRI_short: bim2 = mri_to_short( 1.0 , tim ) ; break ;
119 case MRI_byte: bim2 = mri_to_byte( tim ) ;
120 }
121 mri_write( argv[3] , bim2 ) ;
122 printf(" imrotate %g %g %g %s %s\n",
123 fit[1] , fit[2] , fit[3] , argv[2] , argv[3] ) ;
124 }
125
126 exit(0) ;
127 }
128