Doxygen Source Code Documentation
thd_autonudge.c File Reference
#include "mrilib.h"
Go to the source code of this file.
Defines | |
#define | OV(i, j, k) ov[(i)+(j)*ovx+(k)*ovxy] |
#define | EP(i, j, k) epiar[(i)+(j)*nxepi+(k)*nxyepi] |
#define | AN(i, j, k) antar[(i)+(j)*nxant+(k)*nxyant] |
Functions | |
THD_fvec3 | THD_autonudge (THD_3dim_dataset *dsepi, int ivepi, THD_3dim_dataset *dsant, int ivant, float step, int xstep, int ystep, int zstep, int code) |
Define Documentation
|
Definition at line 11 of file thd_autonudge.c. Referenced by THD_autonudge(). |
|
Definition at line 10 of file thd_autonudge.c. Referenced by THD_autonudge(). |
|
Definition at line 9 of file thd_autonudge.c. |
Function Documentation
|
Definition at line 13 of file thd_autonudge.c. References AN, calloc, DSET_ARRAY, DSET_BRICK_FACTOR, DSET_BRICK_TYPE, DSET_DX, DSET_DY, DSET_DZ, DSET_HEADNAME, DSET_load, DSET_LOADED, DSET_NVALS, DSET_NX, DSET_NY, DSET_NZ, DSET_unload, DSET_XORG, DSET_YORG, DSET_ZORG, DUMP_FVEC3, DUMP_IVEC3, EDIT_coerce_scale_type(), ENTRY, EP, EXIT, free, i, THD_ivec3::ijk, ISVALID_DSET, LOAD_FVEC3, malloc, MIN, mri_clear_data_pointer, mri_fix_data_pointer(), mri_free(), mri_new_vol_empty(), OV, RETURN, THD_3dfind_to_3dmm(), THD_3dmm_to_3dind(), THD_3dmm_to_dicomm(), THD_cliplevel(), THD_dicomm_to_3dmm(), UNLOAD_FVEC3, x2, THD_fvec3::xyz, y1, and z1. Referenced by main().
00016 { 00017 THD_fvec3 fv1,fv2 , fvorg_old,fvorg_new , dxorg ; 00018 THD_ivec3 iv1,iv2 ; 00019 float *tar ; 00020 byte *epiar , *antar ; 00021 int ii,jj,kk , nxepi,nyepi,nzepi , nxyepi,nxyzepi ; 00022 int nxant,nyant,nzant , nxyant,nxyzant ; 00023 MRI_IMAGE *tim ; 00024 float xorgepi , yorgepi , zorgepi , xx1,xx2,yy1,yy2,zz1,zz2 ; 00025 float epiclip , xorgant,yorgant,zorgant , f1,f2,g1,g2,h1,h2 , f,g,h ; 00026 float tx,ty,tz , dxepi,dyepi,dzepi , dxant,dyant,dzant , z1,z2,y1,y2,x1,x2 ; 00027 int *ov,*ovp , ovx=2*xstep+1 , ovy=2*ystep+1 , ovz=2*zstep+1 , ovxy=ovx*ovy ; 00028 int xant,yant,zant , pp,qq,rr , i,j,k , ip,jp,kp , ovtop , kstep ; 00029 float dxyz_ratio , vsum_thresh , vsum , sx,sy,sz ; 00030 int verb = ((code & 1) != 0) ; 00031 00032 /*-- start the action! --*/ 00033 00034 ENTRY("THD_autonudge") ; 00035 00036 /*-- sanity checks --*/ 00037 00038 if( !ISVALID_DSET(dsepi) || 00039 !ISVALID_DSET(dsant) || 00040 ivepi < 0 || ivepi >= DSET_NVALS(dsepi) || 00041 ivant < 0 || ivant >= DSET_NVALS(dsant) || 00042 step <= 0.0 || ovx < 1 || ovy < 1 || ovz < 1 || ovx*ovy*ovz < 3 ){ 00043 00044 fprintf(stderr,"THD_autonudge: bad inputs!\n") ; EXIT(1) ; 00045 } 00046 00047 /*-- load chosen sub-brick of epi into local float array --*/ 00048 00049 if( DSET_ARRAY(dsepi,ivepi) == NULL ){ 00050 DSET_load(dsepi) ; 00051 if( !DSET_LOADED(dsepi) ){ 00052 fprintf(stderr,"THD_autonudge: can't load %s\n",DSET_HEADNAME(dsepi)); 00053 EXIT(1) ; 00054 } 00055 } 00056 00057 nxepi = DSET_NX(dsepi) ; 00058 nyepi = DSET_NY(dsepi) ; nxyepi = nxepi * nyepi ; 00059 nzepi = DSET_NZ(dsepi) ; nxyzepi = nxyepi * nzepi ; 00060 00061 tar = (float *) malloc( sizeof(float) * nxyzepi ) ; 00062 if( tar == NULL ){ 00063 fprintf(stderr,"THD_autonudge: malloc failure for epiar\n"); EXIT(1); 00064 } 00065 00066 EDIT_coerce_scale_type( nxyzepi , 00067 DSET_BRICK_FACTOR(dsepi,ivepi) , 00068 DSET_BRICK_TYPE(dsepi,ivepi) , 00069 DSET_ARRAY(dsepi,ivepi) , MRI_float , tar ) ; 00070 DSET_unload(dsepi) ; 00071 00072 /*-- clip epi array values --*/ 00073 00074 tim = mri_new_vol_empty( nxepi , nyepi , nzepi , MRI_float ) ; 00075 mri_fix_data_pointer( tar , tim ) ; 00076 epiclip = THD_cliplevel( tim , 0.5 ) ; /* get clip value */ 00077 mri_clear_data_pointer(tim) ; mri_free(tim) ; 00078 00079 if( epiclip <= 0.0 ){ 00080 fprintf(stderr,"THD_autonudge: can't compute epiclip\n"); EXIT(1); 00081 } 00082 00083 if( verb ) 00084 fprintf(stderr,"THD_autonudge: epi clip level=%g\n",epiclip) ; 00085 00086 epiar = (byte *) malloc(sizeof(byte)*nxyzepi) ; 00087 if( epiar == 0 ){ 00088 fprintf(stderr,"THD_autonudge: malloc failed for epiar\n"); EXIT(1); 00089 } 00090 00091 for( ii=0 ; ii < nxyzepi ; ii++ ) /* mask of supra-clip voxels */ 00092 epiar[ii] = (tar[ii] > epiclip) ; 00093 00094 free(tar) ; 00095 00096 /*-- load chosen sub-brick of ant into local float array --*/ 00097 00098 if( DSET_ARRAY(dsant,ivant) == NULL ){ 00099 DSET_load(dsant) ; 00100 if( !DSET_LOADED(dsant) ){ 00101 fprintf(stderr,"THD_autonudge: can't load %s\n",DSET_HEADNAME(dsant)); 00102 EXIT(1) ; 00103 } 00104 } 00105 00106 nxant = DSET_NX(dsant) ; 00107 nyant = DSET_NY(dsant) ; nxyant = nxant * nyant ; 00108 nzant = DSET_NZ(dsant) ; nxyzant = nxyant * nzant ; 00109 00110 tar = (float *) malloc( sizeof(float) * nxyzant ) ; 00111 if( tar == NULL ){ 00112 fprintf(stderr,"THD_autonudge: malloc failure for antar\n"); EXIT(1); 00113 } 00114 00115 EDIT_coerce_scale_type( nxyzant , 00116 DSET_BRICK_FACTOR(dsant,ivant) , 00117 DSET_BRICK_TYPE(dsant,ivant) , 00118 DSET_ARRAY(dsant,ivant) , MRI_float , tar ) ; 00119 DSET_unload(dsant) ; 00120 00121 antar = (byte *) malloc(sizeof(byte)*nxyzant) ; 00122 if( antar == NULL ){ 00123 fprintf(stderr,"THD_autonudge: malloc failure for antar\n"); EXIT(1); 00124 } 00125 00126 for( ii=0 ; ii < nxyzant ; ii++ ) /* make mask */ 00127 antar[ii] = (tar[ii] > 0.0) ; 00128 00129 free(tar) ; 00130 00131 /*-- find axis in ant that corresponds to x-axis in epi --*/ 00132 00133 LOAD_FVEC3(fv1,0,0,0) ; 00134 fv1 = THD_3dfind_to_3dmm( dsepi, fv1 ) ; /* coords in dsepi */ 00135 fv1 = THD_3dmm_to_dicomm( dsepi, fv1 ) ; /* DICOM in dsepi */ 00136 fv1 = THD_dicomm_to_3dmm( dsant, fv1 ) ; /* coords in dsant */ 00137 iv1 = THD_3dmm_to_3dind ( dsant, fv1 ) ; /* index in dsant */ 00138 00139 LOAD_FVEC3(fv2,nxepi-1,0,0) ; 00140 fv2 = THD_3dfind_to_3dmm( dsepi , fv2 ) ; /* coords in dsepi */ 00141 fv2 = THD_3dmm_to_dicomm( dsepi , fv2 ) ; /* DICOM in dsepi */ 00142 fv2 = THD_dicomm_to_3dmm( dsant, fv2 ) ; /* coords in dsant */ 00143 iv2 = THD_3dmm_to_3dind ( dsant, fv2 ) ; /* index in dsant */ 00144 00145 if( iv1.ijk[0] != iv2.ijk[0] ) xant = 0 ; /* epi x-axis */ 00146 else if( iv1.ijk[1] != iv2.ijk[1] ) xant = 1 ; /* in ant */ 00147 else if( iv1.ijk[2] != iv2.ijk[2] ) xant = 2 ; 00148 else { 00149 fprintf(stderr,"THD_autonudge: incoherent x slicing?!\n"); 00150 DUMP_IVEC3("iv1",iv1); DUMP_IVEC3("iv2",iv2); EXIT(1); 00151 } 00152 00153 /*-- find axis in ant that corresponds to y-axis in epi --*/ 00154 00155 LOAD_FVEC3(fv2,0,nyepi-1,0) ; 00156 fv2 = THD_3dfind_to_3dmm( dsepi, fv2 ) ; /* coords in dsepi */ 00157 fv2 = THD_3dmm_to_dicomm( dsepi, fv2 ) ; /* DICOM in dsepi */ 00158 fv2 = THD_dicomm_to_3dmm( dsant, fv2 ) ; /* coords in dsant */ 00159 iv2 = THD_3dmm_to_3dind ( dsant, fv2 ) ; /* index in dsant */ 00160 00161 if( iv1.ijk[0] != iv2.ijk[0] ) yant = 0 ; /* epi y-axis */ 00162 else if( iv1.ijk[1] != iv2.ijk[1] ) yant = 1 ; /* in ant */ 00163 else if( iv1.ijk[2] != iv2.ijk[2] ) yant = 2 ; 00164 else { 00165 fprintf(stderr,"THD_autonudge: incoherent y slicing?!\n"); 00166 DUMP_IVEC3("iv1",iv1); DUMP_IVEC3("iv2",iv2); EXIT(1); 00167 } 00168 00169 /*-- find axis in ant that corresponds to z-axis in epi --*/ 00170 00171 LOAD_FVEC3(fv2,0,0,nzepi-1) ; 00172 fv2 = THD_3dfind_to_3dmm( dsepi, fv2 ) ; /* coords in dsepi */ 00173 fv2 = THD_3dmm_to_dicomm( dsepi, fv2 ) ; /* DICOM in dsepi */ 00174 fv2 = THD_dicomm_to_3dmm( dsant, fv2 ) ; /* coords in dsant */ 00175 iv2 = THD_3dmm_to_3dind ( dsant, fv2 ) ; /* index in dsant */ 00176 00177 if( iv1.ijk[0] != iv2.ijk[0] ) zant = 0 ; /* epi z-axis */ 00178 else if( iv1.ijk[1] != iv2.ijk[1] ) zant = 1 ; /* in ant */ 00179 else if( iv1.ijk[2] != iv2.ijk[2] ) zant = 2 ; 00180 else { 00181 fprintf(stderr,"THD_autonudge: incoherent z slicing?!\n"); 00182 DUMP_IVEC3("iv1",iv1); DUMP_IVEC3("iv2",iv2); EXIT(1); 00183 } 00184 00185 #if 0 00186 if( verb ) 00187 fprintf(stderr," xant=%d yant=%d zant=%d\n", 00188 xant,yant,zant) ; 00189 #endif 00190 00191 if( ((1<<xant) | (1<<yant) | (1<<zant)) != 7 ){ 00192 fprintf(stderr,"THD_autonudge: incoherent xyz slicing!\n" 00193 " xant=%d yant=%d zant=%d\n", 00194 xant,yant,zant) ; 00195 EXIT(1) ; 00196 } 00197 00198 /*-- allocate space for array of overlap counts --*/ 00199 00200 ov = (int *) calloc( sizeof(int) , ovx*ovy*ovz ) ; 00201 if( ov == NULL ){ 00202 fprintf(stderr,"THD_autonudge: can't malloc space for overlap counts!\n"); 00203 EXIT(1) ; 00204 } 00205 00206 /*-- for each origin shift in the xant,yant directions, 00207 shift origin of ant dataset, then compare to epi dataset --*/ 00208 00209 xorgepi = DSET_XORG(dsepi) ; dxepi = DSET_DX(dsepi) ; 00210 yorgepi = DSET_YORG(dsepi) ; dyepi = DSET_DY(dsepi) ; 00211 zorgepi = DSET_ZORG(dsepi) ; dzepi = DSET_DZ(dsepi) ; 00212 00213 xorgant = DSET_XORG(dsant) ; dxant = DSET_DX(dsant) ; 00214 yorgant = DSET_YORG(dsant) ; dyant = DSET_DY(dsant) ; 00215 zorgant = DSET_ZORG(dsant) ; dzant = DSET_DZ(dsant) ; 00216 00217 dxyz_ratio = fabs( (dxepi*dyepi*dzepi)/(dxant*dyant*dzant) ) ; 00218 vsum_thresh = 0.5*dxyz_ratio ; 00219 00220 LOAD_FVEC3( fvorg_old,xorgant,yorgant,zorgant) ; 00221 LOAD_FVEC3( dxorg ,dxant ,dyant ,dzant ) ; 00222 00223 ovtop = 0 ; 00224 if( nzepi < 11 ){ 00225 kstep = 2 ; 00226 } else { 00227 kstep = (int)(0.2*nzepi+0.5) ; 00228 } 00229 00230 for( pp=0 ; pp < ovx ; pp++ ){ /* loop over shifts in 3 directions */ 00231 for( qq=0 ; qq < ovy ; qq++ ){ 00232 for( rr=0 ; rr < ovz ; rr++ ){ 00233 00234 if( verb ) fprintf(stderr," starting shift %2d %2d %2d ", 00235 pp-xstep,qq-ystep,rr-zstep) ; 00236 00237 fvorg_new = fvorg_old ; /* old ant origin */ 00238 00239 fvorg_new.xyz[xant] += (pp-xstep)*step*dxorg.xyz[xant]; /* shift epi-x */ 00240 fvorg_new.xyz[yant] += (qq-ystep)*step*dxorg.xyz[yant]; /* shift epi-y */ 00241 fvorg_new.xyz[zant] += (rr-zstep)*step*dxorg.xyz[zant]; /* shift epi-z */ 00242 00243 xorgant = fvorg_new.xyz[0] ; /* load new ant origin */ 00244 yorgant = fvorg_new.xyz[1] ; 00245 zorgant = fvorg_new.xyz[2] ; 00246 00247 ovp = ov + (pp + qq*ovx + rr*ovxy) ; /* place to store result */ 00248 00249 /*-- foreach voxel in epi dataset, 00250 find how much of it is filled by nonzero ant voxels --*/ 00251 00252 for( kk=0 ; kk < nzepi ; kk++ ){ 00253 z1 = zorgepi + dzepi*(kk-0.5) ; z2 = zorgepi + dzepi*(kk+0.49999) ; 00254 00255 if( verb && kk%kstep == 0 ) fprintf(stderr,".") ; 00256 00257 for( jj=0 ; jj < nyepi ; jj++ ){ 00258 y1 = yorgepi + dyepi*(jj-0.5) ; y2 = yorgepi + dyepi*(jj+0.49999) ; 00259 00260 for( ii=0 ; ii < nxepi ; ii++ ){ 00261 if( EP(ii,jj,kk) == 0 ) continue ; /* skip voxel */ 00262 00263 x1 = xorgepi + dxepi*(ii-0.5) ; x2 = xorgepi + dxepi*(ii+0.49999) ; 00264 00265 /* epi voxel covers coords [x1..x2] X [y1..y2] X [z1..z2] */ 00266 00267 /* transform these to ant dataset coords */ 00268 00269 LOAD_FVEC3(fv1,x1,y1,z1) ; /* coords in epi */ 00270 fv1 = THD_3dmm_to_dicomm( dsepi, fv1 ) ; /* DICOM coords */ 00271 fv1 = THD_dicomm_to_3dmm( dsant, fv1 ) ; /* coords in ant */ 00272 UNLOAD_FVEC3(fv1,xx1,yy1,zz1) ; 00273 00274 LOAD_FVEC3(fv2,x2,y2,z2) ; /* coords in epi */ 00275 fv2 = THD_3dmm_to_dicomm( dsepi, fv2 ) ; /* DICOM coords */ 00276 fv2 = THD_dicomm_to_3dmm( dsant, fv2 ) ; /* coords in ant */ 00277 UNLOAD_FVEC3(fv2,xx2,yy2,zz2) ; 00278 00279 /* epi voxel spans ant coords [xx1..xx2] X [yy1..yy2] X [zz1..zz2] */ 00280 00281 /* compute indices into ant dataset voxels */ 00282 00283 f1 = (xx1-xorgant)/dxant + 0.49999 ; f2 = (xx2-xorgant)/dxant + 0.49999 ; 00284 if( f1 > f2 ){ tx = f1 ; f1 = f2 ; f2 = tx ; } 00285 if( f1 >= nxant || f2 <= 0.0 ) continue ; 00286 if( f1 < 0.0 ) f1 = 0.0 ; if( f2 >= nxant ) f2 = nxant - 0.001 ; 00287 00288 g1 = (yy1-yorgant)/dyant + 0.49999 ; g2 = (yy2-yorgant)/dyant + 0.49999 ; 00289 if( g1 > g2 ){ ty = g1 ; g1 = g2 ; g2 = ty ; } 00290 if( g1 >= nyant || g2 <= 0.0 ) continue ; 00291 if( g1 < 0.0 ) g1 = 0.0 ; if( g2 >= nyant ) g2 = nyant - 0.001 ; 00292 00293 h1 = (zz1-zorgant)/dzant + 0.49999 ; h2 = (zz2-zorgant)/dzant + 0.49999 ; 00294 if( h1 > h2 ){ tz = h1 ; h1 = h2 ; h2 = tz ; } 00295 if( h1 >= nzant || h2 <= 0.0 ) continue ; 00296 if( h1 < 0.0 ) h1 = 0.0 ; if( h2 >= nzant ) h2 = nzant - 0.001 ; 00297 00298 /* epi voxel covers voxels [f1..f2] X [g1..g2] X [h1..h2] in ant */ 00299 00300 /* loop over these, and count how much is nonzero */ 00301 00302 vsum = 0.0 ; 00303 for( f=f1 ; f < f2 ; f = ip ){ 00304 i = (int) f ; ip = i+1 ; tx = MIN(ip,f2) ; sx = tx - f ; 00305 for( g=g1 ; g < g2 ; g = jp ){ 00306 j = (int) g ; jp = j+1 ; ty = MIN(jp,g2) ; sy = ty - g ; 00307 for( h=h1 ; h < h2 ; h = kp ){ 00308 k = (int) h ; kp = k+1 ; tz = MIN(kp,h2) ; sz = tz - h ; 00309 if( AN(i,j,k) ) vsum += sx * sy * sz ; 00310 }}} /* end of loop over ant voxels */ 00311 00312 #if 0 00313 fprintf(stderr," epi=%d %d %d ant=%6.2f..%6.2f %6.2f..%6.2f %6.2f..%6.2f vsum=%6.2f\n", 00314 ii,jj,kk , f1,f2 , g1,g2 , h1,h2 , vsum) ; 00315 #endif 00316 00317 /* add to results for this shift */ 00318 00319 if( vsum > vsum_thresh ) (*ovp)++ ; 00320 00321 }}} /* end of loop over epi voxels */ 00322 00323 if( verb ) fprintf(stderr," overlap=%d",*ovp) ; 00324 if( *ovp > ovtop ){ 00325 ovtop = *ovp ; if( verb ) fprintf(stderr," *") ; 00326 } 00327 if( verb ) fprintf(stderr,"\n") ; 00328 00329 }}} /* end of loop over shifts */ 00330 00331 /*-- find best shift in list --*/ 00332 00333 ii = jj = kk = ip = 0 ; 00334 00335 for( pp=0 ; pp < ovx ; pp++ ){ /* loop over shifts in 3 directions */ 00336 for( qq=0 ; qq < ovy ; qq++ ){ 00337 for( rr=0 ; rr < ovz ; rr++ ){ 00338 if( OV(pp,qq,rr) > ip ){ 00339 ii = pp ; jj = qq ; kk = rr ; ip = OV(pp,qq,rr) ; 00340 } 00341 }}} 00342 00343 fvorg_new.xyz[xant] = (ii-xstep)*step*dxorg.xyz[xant]; /* shift epi-x */ 00344 fvorg_new.xyz[yant] = (jj-ystep)*step*dxorg.xyz[yant]; /* shift epi-y */ 00345 fvorg_new.xyz[zant] = (kk-zstep)*step*dxorg.xyz[zant]; /* shift epi-z */ 00346 00347 if( verb ){ 00348 fprintf(stderr," best shift: %d %d %d overlap=%d\n", 00349 ii-xstep,jj-ystep,kk-zstep,ip) ; 00350 DUMP_FVEC3(" best shift",fvorg_new) ; 00351 } 00352 00353 free(ov) ; free(antar) ; free(epiar) ; 00354 00355 RETURN( fvorg_new ) ; 00356 } |