1 #include "mrilib.h"
2 #include <string.h>
3
4 float * mri_align_crao( float filt_fwhm , MRI_IMARR * ims ) ;
5
main(int argc,char * argv[])6 int main( int argc , char * argv[] )
7 {
8 MRI_IMAGE *imin ;
9 MRI_IMARR *ims ;
10 float * dsig ;
11 float fwhm ;
12 int ii , iarg ;
13
14 if( argc < 4 || strncmp(argv[1],"-help",4) == 0 ){
15 printf("Usage: crao fwhm image_files ...\n") ;
16 exit(0) ;
17 }
18
19 fwhm = strtod( argv[1] , NULL ) ;
20 if( fwhm <= 0.0 ){fprintf(stderr,"Illegal fwhm!\a\n");exit(1);}
21
22 INIT_IMARR(ims) ;
23 for( iarg=2 ; iarg < argc ; iarg++ ){
24 imin = mri_read( argv[iarg] ) ;
25 ADDTO_IMARR(ims,imin) ;
26 }
27
28 dsig = mri_align_crao( fwhm , ims ) ;
29 DESTROY_IMARR(ims) ;
30
31 printf("fwhm = %g CR: dx = %g dy = %g phi = %g\n",
32 fwhm,dsig[0],dsig[1],dsig[2] ) ;
33
34 exit(0) ;
35 }
36
37 #define DFAC (PI/180.0)
38
39 /***---------------------------------------------------------------------
40 Compute the Cramer-Rao bounds on registration accuracy from
41 a sequence of images. Returns an array of length 3: dx,dy,phi (degrees).
42 -------------------------------------------------------------------------***/
43
mri_align_crao(float filt_fwhm,MRI_IMARR * ims)44 float * mri_align_crao( float filt_fwhm , MRI_IMARR * ims )
45 {
46 MRI_IMAGE **imstat , *imbar , *imsig , *imdx,*imdy,*imphi ;
47 float sthr,fac , hnx,hny , filt_rms = filt_fwhm*0.42466090 ;
48 float *xar , *yar , *par , *sar , *bar , *crao ;
49 int ii , npix , nx,ny , jj , joff , nzero=0 ;
50 double vdx,vdy,vphi ;
51
52 for( ii=1 ; ii < ims->num ; ii++ )
53 (void) mri_stat_seq( ims->imarr[ii] ) ;
54
55 imstat = mri_stat_seq( NULL ) ;
56 imbar = imstat[0] ;
57 imsig = imstat[1] ;
58
59 nx = imbar->nx ; ny = imbar->ny ; npix = nx * ny ;
60 hnx = 0.5*nx ; hny = 0.5*ny ;
61
62 imdx = mri_filt_fft( imbar , filt_rms , 1 , 0 , FILT_FFT_WRAPAROUND ) ; /* d/dx */
63 imdy = mri_filt_fft( imbar , filt_rms , 0 , 1 , FILT_FFT_WRAPAROUND ) ; /* d/dy */
64 imphi = mri_new( nx , ny , MRI_float ) ;
65
66 xar = MRI_FLOAT_PTR(imdx) ; yar = MRI_FLOAT_PTR(imdy) ;
67 par = MRI_FLOAT_PTR(imphi) ; sar = MRI_FLOAT_PTR(imsig) ;
68
69 for( jj=0 ; jj < ny ; jj++ ){
70 joff = jj * nx ;
71 for( ii=0 ; ii < nx ; ii++ ){
72 par[ii+joff] = DFAC * ( (ii-hnx) * yar[ii+joff]
73 - (jj-hny) * xar[ii+joff] ) ;
74 }
75 }
76
77 sthr = 0.01 * mri_max( imsig ) ;
78 for( ii=0 ; ii < npix ; ii++ )
79 if( sar[ii] < sthr ){
80 sar[ii] = 0.0 ;
81 nzero++ ;
82 }
83 if( nzero > 0 ) printf("set %d sigmas to zero\n",nzero) ;
84
85 for( ii=0 ; ii < npix ; ii++ ){
86 fac = (sar[ii] > 0.0) ? (1.0/SQR(sar[ii])) : 0.0 ;
87 xar[ii] *= fac ;
88 yar[ii] *= fac ;
89 par[ii] *= fac ;
90 }
91
92 mri_free(imbar) ;
93
94 imbar = mri_filt_fft( imdx , filt_rms , 0 , 0 , FILT_FFT_WRAPAROUND ) ;
95 mri_free(imdx) ; imdx = imbar ; xar = MRI_FLOAT_PTR(imdx) ;
96
97 imbar = mri_filt_fft( imdy , filt_rms , 0 , 0 , FILT_FFT_WRAPAROUND ) ;
98 mri_free(imdy) ; imdy = imbar ; yar = MRI_FLOAT_PTR(imdy) ;
99
100 imbar = mri_filt_fft( imphi , filt_rms , 0 , 0 , FILT_FFT_WRAPAROUND ) ;
101 mri_free(imphi) ; imphi = imbar ; par = MRI_FLOAT_PTR(imphi) ;
102
103 for( ii=0 ; ii < npix ; ii++ ){
104 xar[ii] *= sar[ii] ; yar[ii] *= sar[ii] ; par[ii] *= sar[ii] ;
105 }
106
107 vdx = vdy = vphi = 0.0 ;
108
109 for( ii=0 ; ii < npix ; ii++ ){
110 vdx += SQR( xar[ii] ) ;
111 vdy += SQR( yar[ii] ) ;
112 vphi += SQR( par[ii] ) ;
113 }
114
115 mri_free(imsig) ;
116 mri_free(imdx) ; mri_free(imdy) ; mri_free(imphi) ;
117
118 vdx = (vdx > 0.0) ? (1.0/sqrt(vdx)) : 0.0 ;
119 vdy = (vdy > 0.0) ? (1.0/sqrt(vdy)) : 0.0 ;
120 vphi = (vphi > 0.0) ? (1.0/sqrt(vphi)) : 0.0 ;
121
122 crao = (float *) malloc( sizeof(float) * 3 ) ;
123 if( crao == NULL ){fprintf(stderr,"malloc(3) fails in mri_align_crao!\a\n");exit(1);}
124
125 crao[0] = vdx ;
126 crao[1] = vdy ;
127 crao[2] = vphi ;
128
129 return crao ;
130 }
131