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    /*-- 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 }