Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
crao2.c
Go to the documentation of this file.00001 #include "mrilib.h"
00002 #include <string.h>
00003
00004 #include "vecmat.h"
00005
00006 MRI_IMAGE * mri_align_crao( float filt_fwhm , MRI_IMARR * ims ) ;
00007
00008 int main( int argc , char * argv[] )
00009 {
00010 MRI_IMAGE *imin ;
00011 MRI_IMARR *ims ;
00012 float fwhm ;
00013 int ii , iarg ;
00014
00015 if( argc < 4 || strncmp(argv[1],"-help",4) == 0 ){
00016 printf("Usage: crao fwhm image_files ...\n") ;
00017 exit(0) ;
00018 }
00019
00020 fwhm = strtod( argv[1] , NULL ) ;
00021 if( fwhm <= 0.0 ){fprintf(stderr,"Illegal fwhm!\a\n");exit(1);}
00022
00023 INIT_IMARR(ims) ;
00024 for( iarg=2 ; iarg < argc ; iarg++ ){
00025 imin = mri_read( argv[iarg] ) ;
00026 ADDTO_IMARR(ims,imin) ;
00027 }
00028
00029 imin = mri_align_crao( fwhm , ims ) ;
00030 mri_write( "crao.ent" , imin ) ;
00031
00032 exit(0) ;
00033 }
00034
00035 #define DFAC (PI/180.0)
00036
00037
00038
00039
00040
00041
00042 MRI_IMAGE * mri_align_crao( float filt_fwhm , MRI_IMARR * ims )
00043 {
00044 MRI_IMAGE **imstat , *imbar , *imsig , *imdx,*imdy,*imphi , *iment ;
00045 float sthr,fac , hnx,hny , filt_rms = filt_fwhm*0.42466090 ;
00046 float *xar , *yar , *par , *sar , *bar , *entar ;
00047 int ii , npix , nx,ny , jj , joff , nzero=0 ;
00048 THD_dmat33 cov_all , cov_one ;
00049 double vdxdx,vdydy,vphiphi , vdxdy,vdxphi,vdyphi ;
00050 double det_all , det_one ;
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 ) ;
00063 imdy = mri_filt_fft( imbar , filt_rms , 0 , 1 , FILT_FFT_WRAPAROUND ) ;
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 mri_free(imsig) ;
00107
00108 vdxdx = vdydy = vphiphi = vdxdy = vdxphi = vdyphi = 0.0 ;
00109
00110 for( ii=0 ; ii < npix ; ii++ ){
00111 vdxdx += SQR( xar[ii] ) ;
00112 vdydy += SQR( yar[ii] ) ;
00113 vphiphi += SQR( par[ii] ) ;
00114 vdxdy += xar[ii] * yar[ii] ;
00115 vdxphi += xar[ii] * par[ii] ;
00116 vdyphi += yar[ii] * par[ii] ;
00117 }
00118 cov_all.mat[0][0] = vdxdx ;
00119 cov_all.mat[1][1] = vdydy ;
00120 cov_all.mat[2][2] = vphiphi ;
00121 cov_all.mat[0][1] = vdxdy ;
00122 cov_all.mat[1][0] = vdxdy ;
00123 cov_all.mat[0][2] = vdxphi ;
00124 cov_all.mat[2][0] = vdxphi ;
00125 cov_all.mat[1][2] = vdyphi ;
00126 cov_all.mat[2][1] = vdyphi ;
00127
00128 det_all = DMAT_DET(cov_all) ;
00129
00130 iment = mri_new( nx , ny , MRI_float ) ;
00131 entar = MRI_FLOAT_PTR(iment) ;
00132
00133 for( ii=0 ; ii < npix ; ii++ ){
00134 vdxdx = SQR( xar[ii] ) ;
00135 vdydy = SQR( yar[ii] ) ;
00136 vphiphi = SQR( par[ii] ) ;
00137 vdxdy = xar[ii] * yar[ii] ;
00138 vdxphi = xar[ii] * par[ii] ;
00139 vdyphi = yar[ii] * par[ii] ;
00140
00141 cov_one = cov_all ;
00142
00143 cov_one.mat[0][0] -= vdxdx ;
00144 cov_one.mat[1][1] -= vdydy ;
00145 cov_one.mat[2][2] -= vphiphi ;
00146 cov_one.mat[0][1] -= vdxdy ;
00147 cov_one.mat[1][0] -= vdxdy ;
00148 cov_one.mat[0][2] -= vdxphi ;
00149 cov_one.mat[2][0] -= vdxphi ;
00150 cov_one.mat[1][2] -= vdyphi ;
00151 cov_one.mat[2][1] -= vdyphi ;
00152
00153 det_one = DMAT_DET(cov_one) ;
00154
00155 entar[ii] = log( det_all / det_one ) ;
00156 }
00157
00158 mri_free(imdx) ; mri_free(imdy) ; mri_free(imphi) ;
00159
00160 return iment ;
00161 }