00001 #include "mrilib.h"
00002 #include <string.h>
00003 
00004 float * mri_align_crao( float filt_fwhm , MRI_IMARR * ims ) ;
00005 
00006 int main( int argc , char * argv[] )
00007 {
00008    MRI_IMAGE *imin ;
00009    MRI_IMARR *ims ;
00010    float * dsig ;
00011    float fwhm ;
00012    int ii , iarg ;
00013 
00014    if( argc < 4 || strncmp(argv[1],"-help",4) == 0 ){
00015       printf("Usage: crao fwhm image_files ...\n") ;
00016       exit(0) ;
00017    }
00018 
00019    fwhm = strtod( argv[1] , NULL ) ;
00020    if( fwhm <= 0.0 ){fprintf(stderr,"Illegal fwhm!\a\n");exit(1);}
00021 
00022    INIT_IMARR(ims) ;
00023    for( iarg=2 ; iarg < argc ; iarg++ ){
00024       imin  = mri_read( argv[iarg] ) ;
00025       ADDTO_IMARR(ims,imin) ;
00026    }
00027 
00028    dsig = mri_align_crao( fwhm , ims ) ;
00029    DESTROY_IMARR(ims) ;
00030 
00031    printf("fwhm = %g  CR: dx = %g   dy = %g  phi = %g\n",
00032           fwhm,dsig[0],dsig[1],dsig[2] ) ;
00033 
00034    exit(0) ;
00035 }
00036 
00037 #define DFAC         (PI/180.0)
00038 
00039 /***---------------------------------------------------------------------
00040   Compute the Cramer-Rao bounds on registration accuracy from
00041   a sequence of images.  Returns an array of length 3: dx,dy,phi (degrees).
00042 -------------------------------------------------------------------------***/
00043 
00044 float * mri_align_crao( float filt_fwhm , MRI_IMARR * ims )
00045 {
00046    MRI_IMAGE **imstat , *imbar , *imsig , *imdx,*imdy,*imphi ;
00047    float sthr,fac , hnx,hny , filt_rms = filt_fwhm*0.42466090 ;
00048    float *xar , *yar , *par , *sar , *bar , *crao ;
00049    int ii , npix , nx,ny , jj , joff , nzero=0 ;
00050    double vdx,vdy,vphi ;
00051 
00052    for( ii=1 ; ii < ims->num ; ii++ )
00053       (void) mri_stat_seq( ims->imarr[ii] ) ;
00054 
00055    imstat = mri_stat_seq( NULL ) ;
00056    imbar  = imstat[0] ;
00057    imsig  = imstat[1] ;
00058 
00059    nx = imbar->nx ; ny = imbar->ny ; npix = nx * ny ;
00060    hnx = 0.5*nx ; hny = 0.5*ny ;
00061 
00062    imdx   = mri_filt_fft( imbar , filt_rms , 1 , 0 , FILT_FFT_WRAPAROUND ) ;  /* d/dx */
00063    imdy   = mri_filt_fft( imbar , filt_rms , 0 , 1 , FILT_FFT_WRAPAROUND ) ;  /* d/dy */
00064    imphi  = mri_new( nx , ny , MRI_float ) ;
00065 
00066    xar = MRI_FLOAT_PTR(imdx)  ; yar = MRI_FLOAT_PTR(imdy) ;
00067    par = MRI_FLOAT_PTR(imphi) ; sar = MRI_FLOAT_PTR(imsig) ;
00068 
00069    for( jj=0 ; jj < ny ; jj++ ){
00070       joff = jj * nx ;
00071       for( ii=0 ; ii < nx ; ii++ ){
00072          par[ii+joff] = DFAC * (  (ii-hnx) * yar[ii+joff]
00073                                 - (jj-hny) * xar[ii+joff] ) ;
00074       }
00075    }
00076 
00077    sthr = 0.01 * mri_max( imsig ) ;
00078    for( ii=0 ; ii < npix ; ii++ )
00079       if( sar[ii] < sthr ){
00080          sar[ii] = 0.0 ;
00081          nzero++ ;
00082       }
00083    if( nzero > 0 ) printf("set %d sigmas to zero\n",nzero) ;
00084 
00085    for( ii=0 ; ii < npix ; ii++ ){
00086       fac      = (sar[ii] > 0.0) ? (1.0/SQR(sar[ii])) : 0.0 ;
00087       xar[ii] *= fac ;
00088       yar[ii] *= fac ;
00089       par[ii] *= fac ;
00090    }
00091 
00092    mri_free(imbar) ;
00093 
00094    imbar = mri_filt_fft( imdx , filt_rms , 0 , 0 , FILT_FFT_WRAPAROUND ) ;
00095    mri_free(imdx) ; imdx = imbar ; xar = MRI_FLOAT_PTR(imdx) ;
00096 
00097    imbar = mri_filt_fft( imdy , filt_rms , 0 , 0 , FILT_FFT_WRAPAROUND ) ;
00098    mri_free(imdy) ; imdy = imbar ; yar = MRI_FLOAT_PTR(imdy) ;
00099 
00100    imbar = mri_filt_fft( imphi , filt_rms , 0 , 0 , FILT_FFT_WRAPAROUND ) ;
00101    mri_free(imphi) ; imphi = imbar ; par = MRI_FLOAT_PTR(imphi) ;
00102 
00103    for( ii=0 ; ii < npix ; ii++ ){
00104       xar[ii] *= sar[ii] ; yar[ii] *= sar[ii] ; par[ii] *= sar[ii] ;
00105    }
00106 
00107    vdx = vdy = vphi = 0.0 ;
00108 
00109    for( ii=0 ; ii < npix ; ii++ ){
00110       vdx  += SQR( xar[ii] ) ;
00111       vdy  += SQR( yar[ii] ) ;
00112       vphi += SQR( par[ii] ) ;
00113    }
00114 
00115    mri_free(imsig) ;
00116    mri_free(imdx)  ; mri_free(imdy)  ; mri_free(imphi) ;
00117 
00118    vdx  = (vdx  > 0.0) ? (1.0/sqrt(vdx))  : 0.0 ;
00119    vdy  = (vdy  > 0.0) ? (1.0/sqrt(vdy))  : 0.0 ;
00120    vphi = (vphi > 0.0) ? (1.0/sqrt(vphi)) : 0.0 ;
00121 
00122    crao = (float *) malloc( sizeof(float) * 3 ) ;
00123    if( crao == NULL ){fprintf(stderr,"malloc(3) fails in mri_align_crao!\a\n");exit(1);}
00124 
00125    crao[0] = vdx ;
00126    crao[1] = vdy ;
00127    crao[2] = vphi ;
00128 
00129    return crao ;
00130 }