00001 #include "mrilib.h"
00002
00003
00004
00005 #undef OV
00006 #undef EP
00007 #undef AN
00008
00009 #define OV(i,j,k) ov[(i)+(j)*ovx+(k)*ovxy]
00010 #define EP(i,j,k) epiar[(i)+(j)*nxepi+(k)*nxyepi]
00011 #define AN(i,j,k) antar[(i)+(j)*nxant+(k)*nxyant]
00012
00013 THD_fvec3 THD_autonudge( THD_3dim_dataset *dsepi, int ivepi,
00014 THD_3dim_dataset *dsant, int ivant,
00015 float step, int xstep, int ystep, int zstep, int code )
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
00033
00034 ENTRY("THD_autonudge") ;
00035
00036
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
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
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 ) ;
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++ )
00092 epiar[ii] = (tar[ii] > epiclip) ;
00093
00094 free(tar) ;
00095
00096
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++ )
00127 antar[ii] = (tar[ii] > 0.0) ;
00128
00129 free(tar) ;
00130
00131
00132
00133 LOAD_FVEC3(fv1,0,0,0) ;
00134 fv1 = THD_3dfind_to_3dmm( dsepi, fv1 ) ;
00135 fv1 = THD_3dmm_to_dicomm( dsepi, fv1 ) ;
00136 fv1 = THD_dicomm_to_3dmm( dsant, fv1 ) ;
00137 iv1 = THD_3dmm_to_3dind ( dsant, fv1 ) ;
00138
00139 LOAD_FVEC3(fv2,nxepi-1,0,0) ;
00140 fv2 = THD_3dfind_to_3dmm( dsepi , fv2 ) ;
00141 fv2 = THD_3dmm_to_dicomm( dsepi , fv2 ) ;
00142 fv2 = THD_dicomm_to_3dmm( dsant, fv2 ) ;
00143 iv2 = THD_3dmm_to_3dind ( dsant, fv2 ) ;
00144
00145 if( iv1.ijk[0] != iv2.ijk[0] ) xant = 0 ;
00146 else if( iv1.ijk[1] != iv2.ijk[1] ) xant = 1 ;
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
00154
00155 LOAD_FVEC3(fv2,0,nyepi-1,0) ;
00156 fv2 = THD_3dfind_to_3dmm( dsepi, fv2 ) ;
00157 fv2 = THD_3dmm_to_dicomm( dsepi, fv2 ) ;
00158 fv2 = THD_dicomm_to_3dmm( dsant, fv2 ) ;
00159 iv2 = THD_3dmm_to_3dind ( dsant, fv2 ) ;
00160
00161 if( iv1.ijk[0] != iv2.ijk[0] ) yant = 0 ;
00162 else if( iv1.ijk[1] != iv2.ijk[1] ) yant = 1 ;
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
00170
00171 LOAD_FVEC3(fv2,0,0,nzepi-1) ;
00172 fv2 = THD_3dfind_to_3dmm( dsepi, fv2 ) ;
00173 fv2 = THD_3dmm_to_dicomm( dsepi, fv2 ) ;
00174 fv2 = THD_dicomm_to_3dmm( dsant, fv2 ) ;
00175 iv2 = THD_3dmm_to_3dind ( dsant, fv2 ) ;
00176
00177 if( iv1.ijk[0] != iv2.ijk[0] ) zant = 0 ;
00178 else if( iv1.ijk[1] != iv2.ijk[1] ) zant = 1 ;
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
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
00207
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++ ){
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 ;
00238
00239 fvorg_new.xyz[xant] += (pp-xstep)*step*dxorg.xyz[xant];
00240 fvorg_new.xyz[yant] += (qq-ystep)*step*dxorg.xyz[yant];
00241 fvorg_new.xyz[zant] += (rr-zstep)*step*dxorg.xyz[zant];
00242
00243 xorgant = fvorg_new.xyz[0] ;
00244 yorgant = fvorg_new.xyz[1] ;
00245 zorgant = fvorg_new.xyz[2] ;
00246
00247 ovp = ov + (pp + qq*ovx + rr*ovxy) ;
00248
00249
00250
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 ;
00262
00263 x1 = xorgepi + dxepi*(ii-0.5) ; x2 = xorgepi + dxepi*(ii+0.49999) ;
00264
00265
00266
00267
00268
00269 LOAD_FVEC3(fv1,x1,y1,z1) ;
00270 fv1 = THD_3dmm_to_dicomm( dsepi, fv1 ) ;
00271 fv1 = THD_dicomm_to_3dmm( dsant, fv1 ) ;
00272 UNLOAD_FVEC3(fv1,xx1,yy1,zz1) ;
00273
00274 LOAD_FVEC3(fv2,x2,y2,z2) ;
00275 fv2 = THD_3dmm_to_dicomm( dsepi, fv2 ) ;
00276 fv2 = THD_dicomm_to_3dmm( dsant, fv2 ) ;
00277 UNLOAD_FVEC3(fv2,xx2,yy2,zz2) ;
00278
00279
00280
00281
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
00299
00300
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 }}}
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
00318
00319 if( vsum > vsum_thresh ) (*ovp)++ ;
00320
00321 }}}
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 }}}
00330
00331
00332
00333 ii = jj = kk = ip = 0 ;
00334
00335 for( pp=0 ; pp < ovx ; pp++ ){
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];
00344 fvorg_new.xyz[yant] = (jj-ystep)*step*dxorg.xyz[yant];
00345 fvorg_new.xyz[zant] = (kk-zstep)*step*dxorg.xyz[zant];
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 }