Doxygen Source Code Documentation
imfit.c File Reference
#include <string.h>
#include "mrilib.h"
Go to the source code of this file.
Defines | |
#define | DFILT_SIGMA (4.0*0.42466090) |
#define | MAX_ITER 5 |
#define | THRESH 0.10 |
Functions | |
int | main (int argc, char *argv[]) |
Define Documentation
|
Definition at line 10 of file imfit.c. Referenced by main(). |
|
Definition at line 11 of file imfit.c. Referenced by main(). |
|
Definition at line 12 of file imfit.c. Referenced by main(). |
Function Documentation
|
\** File : SUMA.c
Input paramters :
Definition at line 14 of file imfit.c. References ADDTO_IMARR, argc, DFILT_SIGMA, MRI_IMAGE::dx, MRI_IMAGE::dy, fit, INIT_IMARR, MRI_IMAGE::kind, machdep(), MAX_ITER, mri_data_pointer(), mri_filt_fft(), mri_free(), MRI_IS_2D, mri_lsqfit(), mri_new(), mri_read_just_one(), mri_rota(), mri_to_byte(), mri_to_short(), mri_write(), MRI_IMAGE::nx, MRI_IMAGE::ny, and THRESH.
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 ) ; /* smooth */ 00055 xim = mri_filt_fft( im1 , DFILT_SIGMA , 1 , 0 , 0 ) ; /* smooth and d/dx */ 00056 yim = mri_filt_fft( im1 , DFILT_SIGMA , 0 , 1 , 0 ) ; /* smooth and d/dy */ 00057 00058 tim = mri_new( nx , ny , MRI_float ) ; /* x * d/dy - y * d/dx */ 00059 tar = mri_data_pointer( tim ) ; /* which is d/d(theta) */ 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 ) ; /* smooth input image */ 00077 fit = mri_lsqfit( bim2 , fitim , bim2 ) ; /* fit input image to ref images */ 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 ) ; /* smooth input image */ 00091 fit2 = mri_lsqfit( bim2 , fitim , bim2 ) ; /* fit input image to ref images */ 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 } |