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 }