00001 #include "mrilib.h"
00002 
00003 typedef struct {
00004    int   nmask[3] ;
00005    byte * mask[3] ;
00006 } Tmask ;
00007 
00008 void free_Tmask( Tmask * tm )
00009 {
00010    if( tm != NULL ){
00011       free(tm->mask[0]) ; free(tm->mask[1]) ; free(tm->mask[2]) ; free(tm) ;
00012    }
00013    return ;
00014 }
00015 
00016 #define IXY 2  /* fixdir-1 for each plane */
00017 #define IYZ 0
00018 #define IZX 1
00019 
00020 Tmask * create_Tmask( int nx, int ny, int nz, byte * vol )
00021 {
00022    Tmask * tm ;
00023    int ii,jj,kk,vv , nxy,nyz,nzx ;
00024    byte * bz , *xym,*yzm,*zxm , *bxy,*byz,*bzx ;
00025 
00026    tm = (Tmask *) malloc(sizeof(Tmask)) ;
00027    tm->nmask[IXY] = nxy = nx*ny ;
00028    tm->nmask[IYZ] = nyz = ny*nz ;
00029    tm->nmask[IZX] = nzx = nz*nx ;
00030 
00031    tm->mask[IXY] = xym = (byte *) calloc(1,sizeof(byte)*nxy) ;
00032    tm->mask[IYZ] = yzm = (byte *) calloc(1,sizeof(byte)*nyz) ;
00033    tm->mask[IZX] = zxm = (byte *) calloc(1,sizeof(byte)*nzx) ;
00034 
00035    for( byz=yzm,kk=0 ; kk < nz ; kk++,byz+=ny ){
00036       bz = vol + kk*nxy ;
00037       for( bxy=xym,jj=0 ; jj < ny ; jj++,bz+=nx,bxy+=nx ){
00038          for( bzx=zxm,ii=0 ; ii < nx ; ii++,bzx+=nz ){
00039             if( bz[ii] ){ bxy[ii] = byz[jj] = bzx[kk] = 1 ; }
00040          }
00041       }
00042    }
00043 
00044    return tm ;
00045 }
00046 
00047 /*===========================================================================
00048    Functions to extract a plane of shifted bytes from a 3D volume.
00049      nx, ny, nz = dimensions of vol
00050      vol        = input 3D volume of bytes
00051 
00052      fixdir = fixed direction (1=x, 2=y, 3=z)
00053      fixijk = fixed index
00054      da, db = shift in planar coordinaes (non-fixed directions)
00055      ma, mb = dimensions of im
00056      im     = output 2D image
00057 
00058    Goal is im[a,b] = vol[ P(a-da,b-db,c=fixijk) ] for a=0..ma-1, b=0..mb-1,
00059    where P(a,b,c) is the permutation of (a,b,c) that goes with fixdir:
00060      P(x,y,z) = (y,z,x) for fixdir == 1
00061      P(x,y,z) = (z,x,y) for fixdir == 2
00062      P(x,y,z) = (x,y,z) for fixdir == 3
00063    For values outside the range of vol[], im[] is set to 0.
00064 
00065    The five interpolation routines that follow are:
00066      _nn   = nearest neigbhor "interpolation"
00067      _lifl = linear interpolation, with floating point arithmetic
00068      _liby = linear interpolation, with byte arithmetic
00069      _ts   = two-step interpolation
00070      _fs   = four-step interpolation
00071 =============================================================================*/
00072 
00073   /* macros for offsets in vol[] to corners of the interpolation square */
00074 
00075 #undef LL
00076 #undef LR
00077 #undef UL
00078 #undef UR
00079 
00080 #define LL 0                /* lower left  */
00081 #define LR astep            /* lower right */
00082 #define UL bstep            /* upper left  */
00083 #define UR (astep+bstep)    /* upper right */
00084 
00085 #define ASSIGN_DIRECTIONS                                       \
00086  do{ switch( fixdir ){                                          \
00087       default:                                                  \
00088       case 1:            /* x-direction: (a,b,c) = (y,z,x) */   \
00089          astep = nx ; bstep = nxy ; cstep = 1  ;                \
00090          na    = ny ; nb    = nz  ; nc    = nx ;                \
00091       break ;                                                   \
00092                                                                 \
00093       case 2:            /* y-direction: (a,b,c) = (z,x,y) */   \
00094          astep = nxy ; bstep = 1  ; cstep = nx ;                \
00095          na    = nz  ; nb    = nx ; nc    = ny ;                \
00096       break ;                                                   \
00097                                                                 \
00098       case 3:            /* z-direction: (a,b,c) = (x,y,z) */   \
00099          astep = 1  ; bstep = nx ; cstep = nxy ;                \
00100          na    = nx ; nb    = ny ; nc    = nz  ;                \
00101       break ;                                                   \
00102     } } while(0)
00103 
00104 /*-----------------------------------------------------------------------*/
00105 
00106 void extract_assign_directions( int nx, int ny, int nz, int fixdir ,
00107                                 int *Astep, int *Bstep, int *Cstep ,
00108                                 int *Na   , int *Nb   , int *Nc     )
00109 {
00110    int astep,bstep,cstep , na,nb,nc , nxy=nx*ny ;
00111 
00112    ASSIGN_DIRECTIONS ;
00113 
00114    *Astep = astep ; *Bstep = bstep ; *Cstep = cstep ;
00115    *Na    = na    ; *Nb    = nb    ; *Nc    = nc    ; return ;
00116 }
00117 
00118 /*-----------------------------------------------------------------------
00119    NN "interpolation"
00120 -------------------------------------------------------------------------*/
00121 
00122 void extract_byte_nn( int nx , int ny , int nz , byte * vol ,
00123                       Tmask * tm ,
00124                       int fixdir , int fixijk , float da , float db ,
00125                       int ma , int mb , byte * im )
00126 {
00127    int adel,bdel , abot,atop , bb,bbot,btop , nxy=nx*ny ;
00128    register int aa , ijkoff , aoff,boff ;
00129    int astep,bstep,cstep , na,nb,nc ;
00130    byte * mask ;
00131 
00132    memset( im , 0 , ma*mb ) ;  /* initialize output to zero */
00133 
00134    if( fixijk < 0 ) return ;
00135 
00136    ASSIGN_DIRECTIONS ;
00137 
00138    if( fixijk >= nc ) return ;
00139 
00140    da += 0.5 ; adel = (int) da ; if( da < 0.0 ) adel-- ;  /* floor(da+0.5) */
00141    db += 0.5 ; bdel = (int) db ; if( db < 0.0 ) bdel-- ;  /* floor(db+0.5) */
00142 
00143    abot = 0       ; if( abot < adel ) abot = adel ;       /* range in im[] */
00144    atop = na+adel ; if( atop > ma   ) atop = ma ;
00145 
00146    bbot = 0       ; if( bbot < bdel ) bbot = bdel ;
00147    btop = nb+bdel ; if( btop > mb   ) btop = mb ;
00148 
00149    ijkoff = fixijk*cstep + (abot-adel)*astep + (bbot-bdel)*bstep ;
00150    boff   = bbot * ma ;
00151 
00152    mask = (tm == NULL) ? NULL
00153                        : tm->mask[fixdir%3] + (fixijk*nb - bdel) ;
00154 
00155    for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
00156       if( mask == NULL || mask[bb] )
00157          for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
00158             im[aa+boff] = vol[aoff+ijkoff] ;
00159 
00160    return ;
00161 }
00162 
00163 /*---------------------------------------------------------------------------
00164     Linear interpolation with floating point arithmetic
00165 -----------------------------------------------------------------------------*/
00166 
00167 void extract_byte_lifl( int nx , int ny , int nz , byte * vol ,
00168                         Tmask * tm ,
00169                         int fixdir , int fixijk , float da , float db ,
00170                         int ma , int mb , byte * im )
00171 {
00172    int adel,bdel , abot,atop , bb,bbot,btop , nxy=nx*ny ;
00173    register int aa , ijkoff , aoff,boff ;
00174    int astep,bstep,cstep , na,nb,nc ;
00175    float fa , fb ;
00176    float f_a_b , f_ap_b , f_a_bp , f_ap_bp ;
00177    byte * mask ;
00178 
00179    memset( im , 0 , ma*mb ) ;  /* initialize output to zero */
00180 
00181    if( fixijk < 0 ) return ;
00182 
00183    ASSIGN_DIRECTIONS ;
00184 
00185    if( fixijk >= nc ) return ;
00186 
00187    adel = (int) da ; if( da < 0.0 ) adel-- ;  /* floor(da) */
00188    bdel = (int) db ; if( db < 0.0 ) bdel-- ;  /* floor(db) */
00189 
00190    fa = da - adel ;               /* fractional part of dj */
00191    fb = db - bdel ;               /* fractional part of dk */
00192 
00193    adel++ ; bdel++ ;
00194 
00195    f_a_b   = fa      * fb      ;
00196    f_ap_b  = (1.0-fa)* fb      ;
00197    f_a_bp  = fa      *(1.0-fb) ;
00198    f_ap_bp = (1.0-fa)*(1.0-fb) ;
00199 
00200    abot = 0         ; if( abot < adel ) abot = adel ;       /* range in im[] */
00201    atop = na+adel-1 ; if( atop > ma   ) atop = ma ;
00202 
00203    bbot = 0         ; if( bbot < bdel ) bbot = bdel ;
00204    btop = nb+bdel-1 ; if( btop > mb   ) btop = mb ;
00205 
00206    ijkoff = fixijk*cstep + (abot-adel)*astep + (bbot-bdel)*bstep ;
00207    boff   = bbot * ma ;
00208 
00209    mask = (tm == NULL) ? NULL
00210                        : tm->mask[fixdir%3] + (fixijk*nb - bdel) ;
00211 
00212    for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
00213       if( mask == NULL || mask[bb] || mask[bb+1] )
00214          for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
00215             im[aa+boff] = (byte)(  f_a_b   * vol[aoff+ijkoff]
00216                                  + f_ap_b  * vol[aoff+(ijkoff+LR)]
00217                                  + f_a_bp  * vol[aoff+(ijkoff+UL)]
00218                                  + f_ap_bp * vol[aoff+(ijkoff+UR)] ) ;
00219    return ;
00220 }
00221 
00222 /*---------------------------------------------------------------------------
00223     Linear interpolation with fixed point arithmetic
00224 -----------------------------------------------------------------------------*/
00225 
00226 void extract_byte_liby( int nx , int ny , int nz , byte * vol ,
00227                         Tmask * tm ,
00228                         int fixdir , int fixijk , float da , float db ,
00229                         int ma , int mb , byte * im )
00230 {
00231    int adel,bdel , abot,atop , bb,bbot,btop , nxy=nx*ny ;
00232    register int aa , ijkoff , aoff,boff ;
00233    int astep,bstep,cstep , na,nb,nc ;
00234    float fa , fb ;
00235    float f_a_b , f_ap_b , f_a_bp , f_ap_bp ;
00236    byte  b_a_b , b_ap_b , b_a_bp , b_ap_bp ;
00237    byte * mask ;
00238 
00239    memset( im , 0 , ma*mb ) ;  /* initialize output to zero */
00240 
00241    if( fixijk < 0 ) return ;
00242 
00243    ASSIGN_DIRECTIONS ;
00244 
00245    if( fixijk >= nc ) return ;
00246 
00247    adel = (int) da ; if( da < 0.0 ) adel-- ;  /* floor(da) */
00248    bdel = (int) db ; if( db < 0.0 ) bdel-- ;  /* floor(db) */
00249 
00250    fa = da - adel ;               /* fractional part of dj */
00251    fb = db - bdel ;               /* fractional part of dk */
00252 
00253    adel++ ; bdel++ ;
00254 
00255    f_a_b   = fa      * fb      ;
00256    f_ap_b  = (1.0-fa)* fb      ;
00257    f_a_bp  = fa      *(1.0-fb) ;
00258    f_ap_bp = (1.0-fa)*(1.0-fb) ;
00259 
00260    bb = (int)(256*f_a_b  + 0.499) ; if( bb == 256 ) bb-- ; b_a_b  = (byte) bb ;
00261    bb = (int)(256*f_ap_b + 0.499) ; if( bb == 256 ) bb-- ; b_ap_b = (byte) bb ;
00262    bb = (int)(256*f_a_bp + 0.499) ; if( bb == 256 ) bb-- ; b_a_bp = (byte) bb ;
00263    bb = (int)(256*f_ap_bp+ 0.499) ; if( bb == 256 ) bb-- ; b_ap_bp= (byte) bb ;
00264 
00265    abot = 0         ; if( abot < adel ) abot = adel ;       /* range in im[] */
00266    atop = na+adel-1 ; if( atop > ma   ) atop = ma ;
00267 
00268    bbot = 0         ; if( bbot < bdel ) bbot = bdel ;
00269    btop = nb+bdel-1 ; if( btop > mb   ) btop = mb ;
00270 
00271    ijkoff = fixijk*cstep + (abot-adel)*astep + (bbot-bdel)*bstep ;
00272    boff   = bbot * ma ;
00273 
00274    mask = (tm == NULL) ? NULL
00275                        : tm->mask[fixdir%3] + (fixijk*nb - bdel) ;
00276 
00277    for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
00278       if( mask == NULL || mask[bb] || mask[bb+1] )
00279         for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
00280            im[aa+boff] = (byte)((  b_a_b   * vol[aoff+ijkoff]
00281                                  + b_ap_b  * vol[aoff+(ijkoff+LR)]
00282                                  + b_a_bp  * vol[aoff+(ijkoff+UL)]
00283                                  + b_ap_bp * vol[aoff+(ijkoff+UR)] ) >> 8 ) ;
00284    return ;
00285 }
00286 
00287 /*---------------------------------------------------------------------------
00288     Two-step interpolation
00289 -----------------------------------------------------------------------------*/
00290 
00291 #if 0
00292 # define TSBOT 0.3
00293 # define TSTOP 0.7
00294 #else
00295 # define TSBOT 0.25
00296 # define TSTOP 0.75
00297 #endif
00298 
00299 void extract_byte_ts( int nx , int ny , int nz , byte * vol ,
00300                       Tmask * tm ,
00301                       int fixdir , int fixijk , float da , float db ,
00302                       int ma , int mb , byte * im )
00303 {
00304    int adel,bdel , abot,atop , bb,bbot,btop , nxy=nx*ny ;
00305    register int aa , ijkoff , aoff,boff ;
00306    int astep,bstep,cstep , na,nb,nc , nts,dts1,dts2 ;
00307    float fa , fb ;
00308    byte * mask ;
00309 
00310    memset( im , 0 , ma*mb ) ;  /* initialize output to zero */
00311 
00312    if( fixijk < 0 ) return ;
00313 
00314    ASSIGN_DIRECTIONS ;
00315 
00316    if( fixijk >= nc ) return ;
00317 
00318    adel = (int) da ; if( da < 0.0 ) adel-- ;  /* floor(da) */
00319    bdel = (int) db ; if( db < 0.0 ) bdel-- ;  /* floor(db) */
00320 
00321    fa = da - adel ;               /* fractional part of dj */
00322    fb = db - bdel ;               /* fractional part of dk */
00323 
00324    fa = 1.0-fa ; fb = 1.0-fb ;
00325 
00326    if( fa < TSBOT ){                      /*- Left 30% -*/
00327       if( fb < TSBOT ){                   /*- Lower 30% -*/
00328         nts = 1 ; dts1 = LL ;               /* [0,0] */
00329       } else if( fb > TSTOP ){            /*- Upper 30% -*/
00330         nts = 1 ; dts1 = UL ;               /* [0,1] */
00331       } else {                            /*- Middle 40% -*/
00332         nts = 2 ; dts1 = LL ; dts2 = UL ;   /* mid of [0,0] and [0,1] */
00333       }
00334    } else if( fa > TSTOP ){               /*- Right 30% -*/
00335       if( fb < TSBOT ){                   /*- Lower 30% -*/
00336         nts = 1 ; dts1 = LR ;               /* [1,0] */
00337       } else if( fb > TSTOP ){            /*- Upper 30% -*/
00338         nts = 1 ; dts1 = UR ;               /* [1,1] */
00339       } else {                            /*- Middle 40% -*/
00340         nts = 2 ; dts1 = LR ; dts2 = UR ;   /* mid of [1,0] and [1,1] */
00341       }
00342    } else {                               /*- Middle 40% -*/
00343       if( fb < TSBOT ){                   /*- Lower 30% -*/
00344         nts = 2 ; dts1 = LL ; dts2 = LR ;   /* mid of [0,0] and [1,0] */
00345       } else if( fb > TSTOP ){            /*- Upper 30% -*/
00346         nts = 2 ; dts1 = UL ; dts2 = UR ;   /* mid of [0,1] and [1,1] */
00347       } else {                            /*- Middle 40% -*/
00348         nts = 4 ;                           /* mid of all 4 points */
00349       }
00350    }
00351 
00352    adel++ ; bdel++ ;
00353 
00354    abot = 0         ; if( abot < adel ) abot = adel ;       /* range in im[] */
00355    atop = na+adel-1 ; if( atop > ma   ) atop = ma ;
00356 
00357    bbot = 0         ; if( bbot < bdel ) bbot = bdel ;
00358    btop = nb+bdel-1 ; if( btop > mb   ) btop = mb ;
00359 
00360    ijkoff = fixijk*cstep + (abot-adel)*astep + (bbot-bdel)*bstep ;
00361    boff   = bbot * ma ;
00362 
00363    mask = (tm == NULL) ? NULL
00364                        : tm->mask[fixdir%3] + (fixijk*nb - bdel) ;
00365 
00366    switch( nts ){
00367 
00368       case 1:
00369          ijkoff += dts1 ;
00370          for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
00371            if( mask == NULL || mask[bb] || mask[bb+1] )
00372              for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
00373                im[aa+boff] = vol[aoff+ijkoff] ;
00374       break ;
00375 
00376       case 2:
00377          ijkoff += dts1 ; dts2 -= dts1 ;
00378          for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
00379            if( mask == NULL || mask[bb] || mask[bb+1] )
00380              for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
00381                im[aa+boff] = (vol[aoff+ijkoff] + vol[aoff+(ijkoff+dts2)]) >> 1;
00382       break ;
00383 
00384       case 4:
00385          for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
00386            if( mask == NULL || mask[bb] || mask[bb+1] )
00387              for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
00388                im[aa+boff] = ( vol[aoff+ijkoff]     +vol[aoff+(ijkoff+LR)]
00389                               +vol[aoff+(ijkoff+UL)]+vol[aoff+(ijkoff+UR)]) >> 2;
00390       break ;
00391    }
00392 
00393    return ;
00394 }
00395 
00396 /*---------------------------------------------------------------------------
00397     Four-step interpolation
00398 -----------------------------------------------------------------------------*/
00399 
00400 #if 0
00401 # define FSA 0.175
00402 # define FSB 0.400
00403 # define FSC 0.600
00404 # define FSD 0.825
00405 #else
00406 # define FSA 0.125
00407 # define FSB 0.375
00408 # define FSC 0.625
00409 # define FSD 0.875
00410 #endif
00411 
00412 void extract_byte_fs( int nx , int ny , int nz , byte * vol ,
00413                       Tmask * tm ,
00414                       int fixdir , int fixijk , float da , float db ,
00415                       int ma , int mb , byte * im )
00416 {
00417    int adel,bdel , abot,atop , bb,bbot,btop , nxy=nx*ny ;
00418    register int aa , ijkoff , aoff,boff ;
00419    int astep,bstep,cstep , na,nb,nc , nfs,dfs1,dfs2,dfs3,dfs4 , ap,bp ;
00420    float fa , fb ;
00421    byte * mask ;
00422 
00423    memset( im , 0 , ma*mb ) ;  /* initialize output to zero */
00424 
00425    if( fixijk < 0 ) return ;
00426 
00427    ASSIGN_DIRECTIONS ;
00428 
00429    if( fixijk >= nc ) return ;
00430 
00431    adel = (int) da ; if( da < 0.0 ) adel-- ;  /* floor(da) */
00432    bdel = (int) db ; if( db < 0.0 ) bdel-- ;  /* floor(db) */
00433 
00434    fa = da - adel ;               /* fractional part of dj */
00435    fb = db - bdel ;               /* fractional part of dk */
00436 
00437    fa = 1.0-fa ; fb = 1.0-fb ;   /* weights for right/upper sides */
00438 
00439         if( fa < FSA ) ap = 0 ;  /* left-right position */
00440    else if( fa < FSB ) ap = 1 ;
00441    else if( fa < FSC ) ap = 2 ;
00442    else if( fa < FSD ) ap = 3 ;
00443    else                ap = 4 ;
00444 
00445         if( fb < FSA ) bp = 0 ;  /* down-up position */
00446    else if( fb < FSB ) bp = 1 ;
00447    else if( fb < FSC ) bp = 2 ;
00448    else if( fb < FSD ) bp = 3 ;
00449    else                bp = 4 ;
00450 
00451    /*----- 5x5 grid of possible interpolation cases (nfs): -----------------
00452 
00453                    bp = 4|  1 3 2 3 1     04 14 24 34 44 <- grid of
00454                         3|  3 4 5 4 3     03 13 23 33 43 <- 10*ap + bp
00455                         2|  2 5 6 5 2     02 12 22 32 42 <- values
00456                         1|  3 4 5 4 3     01 11 21 31 41
00457                         0|  1 3 2 3 1     00 10 20 30 40
00458                            -----------
00459                        ap = 0 1 2 3 4
00460 
00461      ----- The indices and nfs cases are assigned in the switch below. -----*/
00462 
00463 
00464    dfs2=dfs3=dfs4=-1 ;
00465    switch( 10*ap + bp ){
00466 
00467       default: return ;  /* should never be executed */
00468 
00469       case 00: nfs = 1 ; dfs1 = LL ; break ;              /* 1 point */
00470       case 04: nfs = 1 ; dfs1 = UL ; break ;
00471       case 40: nfs = 1 ; dfs1 = LR ; break ;
00472       case 44: nfs = 1 ; dfs1 = UR ; break ;
00473 
00474       case 20: nfs = 2 ; dfs1 = LL ; dfs2 = LR ; break ;  /* 2 points:  */
00475       case 02: nfs = 2 ; dfs1 = LL ; dfs2 = UL ; break ;  /* 1/2 = dfs1 */
00476       case 24: nfs = 2 ; dfs1 = UL ; dfs2 = UR ; break ;  /* 1/2 = dfs2 */
00477       case 42: nfs = 2 ; dfs1 = LR ; dfs2 = UR ; break ;
00478 
00479       case 10: nfs = 3 ; dfs1 = LL ; dfs2 = LR ; break ;  /* 2 points:  */
00480       case 30: nfs = 3 ; dfs1 = LR ; dfs2 = LL ; break ;  /* 3/4 = dfs1 */
00481       case 01: nfs = 3 ; dfs1 = LL ; dfs2 = UL ; break ;  /* 1/4 = dfs2 */
00482       case 03: nfs = 3 ; dfs1 = UL ; dfs2 = LL ; break ;
00483       case 14: nfs = 3 ; dfs1 = UL ; dfs2 = UR ; break ;
00484       case 34: nfs = 3 ; dfs1 = UR ; dfs2 = UL ; break ;
00485       case 41: nfs = 3 ; dfs1 = LR ; dfs2 = UR ; break ;
00486       case 43: nfs = 3 ; dfs1 = UR ; dfs2 = LR ; break ;
00487 
00488       case 11: nfs = 4 ; dfs1 = LL ; dfs2 = LR ;          /* 4 points:   */
00489                          dfs3 = UL ; dfs4 = UR ; break ;  /* 9/16 = dfs1 */
00490       case 13: nfs = 4 ; dfs1 = UL ; dfs2 = UR ;          /* 3/16 = dfs2 */
00491                          dfs3 = LL ; dfs4 = LR ; break ;  /* 3/16 = dfs3 */
00492       case 31: nfs = 4 ; dfs1 = LR ; dfs2 = LL ;          /* 1/16 = dfs4 */
00493                          dfs3 = UR ; dfs4 = UL ; break ;
00494       case 33: nfs = 4 ; dfs1 = UR ; dfs2 = UL ;
00495                          dfs3 = LR ; dfs4 = LL ; break ;
00496 
00497       case 12: nfs = 5 ; dfs1 = LL ; dfs2 = UL ;          /* 4 points:  */
00498                          dfs3 = LR ; dfs4 = UR ; break ;  /* 3/8 = dfs1 */
00499       case 21: nfs = 5 ; dfs1 = LL ; dfs2 = LR ;          /* 3/8 = dfs2 */
00500                          dfs3 = UL ; dfs4 = UR ; break ;  /* 1/8 = dfs3 */
00501       case 23: nfs = 5 ; dfs1 = UL ; dfs2 = UR ;          /* 1/8 = dfs4 */
00502                          dfs3 = LL ; dfs4 = LR ; break ;
00503       case 32: nfs = 5 ; dfs1 = LR ; dfs2 = UR ;
00504                          dfs3 = LL ; dfs4 = UL ; break ;
00505 
00506       case 22: nfs = 6 ; dfs1 = LL ; dfs2 = LR ;          /* 4 points: */
00507                          dfs3 = UL ; dfs4 = UR ; break ;  /* 1/4 = all */
00508    }
00509 
00510    adel++ ; bdel++ ;
00511 
00512    abot = 0         ; if( abot < adel ) abot = adel ;       /* range in im[] */
00513    atop = na+adel-1 ; if( atop > ma   ) atop = ma ;
00514 
00515    bbot = 0         ; if( bbot < bdel ) bbot = bdel ;
00516    btop = nb+bdel-1 ; if( btop > mb   ) btop = mb ;
00517 
00518    ijkoff = fixijk*cstep + (abot-adel)*astep + (bbot-bdel)*bstep ;
00519    boff   = bbot * ma ;
00520 
00521 #if 0
00522 printf("fixijk=%3d  nfs=%d  dfs1=%d  dfs2=%d  dfs3=%d  dfs4=%d\n",
00523         fixijk,nfs,dfs1,dfs2,dfs3,dfs4);
00524 #endif
00525 
00526    mask = (tm == NULL) ? NULL
00527                        : tm->mask[fixdir%3] + (fixijk*nb - bdel) ;
00528 
00529    switch( nfs ){
00530 
00531       case 1:                                          /* 1 point (NN copy) */
00532          ijkoff += dfs1 ;
00533          for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
00534            if( mask == NULL || mask[bb] || mask[bb+1] )
00535             for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
00536                im[aa+boff] = vol[aoff+ijkoff] ;
00537       break ;
00538 
00539       case 2:                                          /* 2 points (1/2+1/2) */
00540          ijkoff += dfs1 ; dfs2 -= dfs1 ;
00541          for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
00542            if( mask == NULL || mask[bb] || mask[bb+1] )
00543             for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
00544                im[aa+boff] = (vol[aoff+ijkoff] + vol[aoff+(ijkoff+dfs2)]) >> 1 ;
00545       break ;
00546 
00547       case 3:                                          /* 2 points (3/4+1/4) */
00548          ijkoff += dfs1 ; dfs2 -= dfs1 ;
00549          for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
00550            if( mask == NULL || mask[bb] || mask[bb+1] )
00551             for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
00552                im[aa+boff] = ( (vol[aoff+ijkoff] << 1) + vol[aoff+ijkoff]
00553                               + vol[aoff+(ijkoff+dfs2)]                  ) >> 2 ;
00554       break ;
00555 
00556       case 4:                                          /* 4 points (9/16+3/16+3/16+1/16) */
00557          ijkoff += dfs1 ; dfs2 -= dfs1 ; dfs3 -= dfs1 ; dfs4 -= dfs1 ;
00558          for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
00559            if( mask == NULL || mask[bb] || mask[bb+1] )
00560             for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
00561                im[aa+boff] = ( (vol[aoff+ijkoff] << 3)
00562                               + vol[aoff+ijkoff]
00563                               +((vol[aoff+(ijkoff+dfs2)] + vol[aoff+(ijkoff+dfs3)]) << 1)
00564                               + (vol[aoff+(ijkoff+dfs2)] + vol[aoff+(ijkoff+dfs3)])
00565                               + vol[aoff+(ijkoff+dfs4)]                                ) >> 4 ;
00566       break ;
00567 
00568       case 5:                                          /* 4 points (3/8+3/8+1/8+1/8) */
00569          ijkoff += dfs1 ; dfs2 -= dfs1 ; dfs3 -= dfs1 ; dfs4 -= dfs1 ;
00570          for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
00571            if( mask == NULL || mask[bb] || mask[bb+1] )
00572             for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
00573                im[aa+boff] = ( ((vol[aoff+ijkoff] + vol[aoff+(ijkoff+dfs2)]) << 1)
00574                               + (vol[aoff+ijkoff] + vol[aoff+(ijkoff+dfs2)])
00575                               + vol[aoff+(ijkoff+dfs3)] + vol[aoff+(ijkoff+dfs4)] ) >> 3 ;
00576       break;
00577 
00578       case 6:                                          /* 4 points (1/4+1/4+1/4+1/4) */
00579          ijkoff += dfs1 ; dfs2 -= dfs1 ; dfs3 -= dfs1 ; dfs4 -= dfs1 ;
00580          for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
00581            if( mask == NULL || mask[bb] || mask[bb+1] )
00582             for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
00583                im[aa+boff] = (  vol[aoff+ijkoff]        + vol[aoff+(ijkoff+dfs2)]
00584                               + vol[aoff+(ijkoff+dfs3)] + vol[aoff+(ijkoff+dfs4)] ) >> 2 ;
00585       break;
00586    }
00587 
00588    return ;
00589 }
00590 
00591 /*---------------------------------------------------------------------------
00592     Test the speeds of the above routines:
00593       nrep = number of repetitions to execute
00594       ct   = float [5] array (must be allocated by caller)
00595              ct[0] = CPU time for _nn
00596              ct[1] = CPU time for _lifl
00597              ct[2] = CPU time for _liby
00598              ct[3] = CPU time for _ts
00599              ct[4] = CPU time for _fs
00600 -----------------------------------------------------------------------------*/
00601 
00602 void extract_byte_speedtest( int nrep , int fixdir , float * ct )
00603 {
00604    double cputim ;
00605    int pp , nx=161,ny=191,nz=141,nxy=nx*ny ,
00606        kk , ma,mb,mab , apad,bpad ;
00607    float aa=0.347 , bb=-0.521 , da,db ;
00608    byte * vin , * vout ;
00609    int astep,bstep,cstep , na,nb,nc ;
00610 
00611    ASSIGN_DIRECTIONS ;
00612 
00613    /* setup bricks */
00614 
00615    da = fabs( 0.5*aa*(nc-1.0) ) ; db = fabs( 0.5*bb*(nc-1.0) ) ;
00616    apad = (int)(2.0+da)         ; bpad = (int)(2.0+db) ;
00617    ma   = na + 2*apad           ; mb   = nb + 2*bpad   ; mab = ma*mb ;
00618 
00619    vin = (byte *) malloc( sizeof(byte) * (na*nb*nc) ) ;
00620    if( vin == NULL ) return ;
00621 
00622    vout = (byte *) malloc( sizeof(byte) * (ma*mb*nc) ) ;
00623    if( vout == NULL ){ free(vin) ; return ; }
00624 
00625    vin[0] = 1 ;
00626    for( kk=1 ; kk < na*nb*nc ; kk++ ) vin[kk] = (byte)((3*vin[kk-1]+7) % 256) ;
00627 
00628 #undef BTEST
00629 #define BTEST(func) do{ cputim = COX_cpu_time() ;                    \
00630                         for( pp=0 ; pp < nrep ; pp++ ){              \
00631                           for( kk=0 ; kk < nc ; kk++ ){              \
00632                              da = aa*(kk - 0.5*(nc-1.0)) + apad ;    \
00633                              db = bb*(kk - 0.5*(nc-1.0)) + bpad ;    \
00634                              func( nx,ny,nz , vin ,                  \
00635                                    NULL ,                            \
00636                                    fixdir , kk , da , db ,           \
00637                                    ma , mb , vout + kk*mab ) ;       \
00638                           }                                          \
00639                         }                                            \
00640                         cputim = COX_cpu_time() - cputim ; } while(0)
00641 
00642    BTEST(extract_byte_nn)   ; ct[0] = cputim ;
00643    BTEST(extract_byte_lifl) ; ct[1] = cputim ;
00644    BTEST(extract_byte_liby) ; ct[2] = cputim ;
00645    BTEST(extract_byte_ts)   ; ct[3] = cputim ;
00646    BTEST(extract_byte_fs)   ; ct[4] = cputim ;
00647 
00648 #undef BTEST
00649 
00650    free(vin) ; free(vout) ; return ;
00651 }
00652 
00653 /*-----------------------------------------------------------------------
00654    Simple get/put of a fixed plane (no shifting, zero padding).
00655 -------------------------------------------------------------------------*/
00656 
00657 void getplane_byte( int nx , int ny , int nz , byte * vol ,
00658                     int fixdir , int fixijk , byte * im )
00659 {
00660    int bb , nxy=nx*ny ;
00661    register int aa , ijkoff , aoff,boff ;
00662    int astep,bstep,cstep , na,nb,nc ;
00663 
00664    if( fixijk < 0 ) return ;
00665 
00666    ASSIGN_DIRECTIONS ;
00667 
00668    if( fixijk >= nc ) return ;
00669 
00670    ijkoff = fixijk*cstep ;
00671 
00672    for( bb=0,boff=0 ; bb < nb ; bb++,boff+=na,ijkoff+=bstep )
00673       for( aa=0,aoff=0 ; aa < na ; aa++,aoff+=astep )
00674          im[aa+boff] = vol[aoff+ijkoff] ;
00675 
00676    return ;
00677 }
00678 
00679 void putplane_byte( int nx , int ny , int nz , byte * vol ,
00680                     int fixdir , int fixijk , byte * im )
00681 {
00682    int bb , nxy=nx*ny ;
00683    register int aa , ijkoff , aoff,boff ;
00684    int astep,bstep,cstep , na,nb,nc ;
00685 
00686    if( fixijk < 0 ) return ;
00687 
00688    ASSIGN_DIRECTIONS ;
00689 
00690    if( fixijk >= nc ) return ;
00691 
00692    ijkoff = fixijk*cstep ;
00693 
00694    for( bb=0,boff=0 ; bb < nb ; bb++,boff+=na,ijkoff+=bstep )
00695       for( aa=0,aoff=0 ; aa < na ; aa++,aoff+=astep )
00696          vol[aoff+ijkoff] = im[aa+boff] ;
00697 
00698    return ;
00699 }
00700 
00701 /******************************************************************************
00702  ******************************************************************************
00703  ******************************************************************************/
00704 
00705 typedef void gfun( int , int , int , byte * , Tmask * ,
00706                    int , int , float , float , int , int , byte * ) ;
00707 
00708 int main( int argc , char * argv[] )
00709 {
00710    THD_3dim_dataset * in_dset ;
00711    Tmask * tmask ;
00712    int nx,ny,nz,nxy , kk,ii , ma,mb,mab ,
00713        apad,bpad , pp,ploop=1,fixdir;
00714    float aa , bb , da,db ;
00715    THD_ivec3 iv ;
00716    byte * vin , * vout , * vmax ;
00717    MRI_IMAGE * imout , * immax ;
00718    double cputim ;
00719    gfun * func = extract_byte_nn ;
00720    char * cfun = "nn" ;
00721    int astep,bstep,cstep , na,nb,nc , use_tmask ;
00722 
00723    if( argc < 3 ){
00724       printf("Usage 1: extor fixdir A B bytedset [loops [suffix]]\n") ;
00725       printf("Usage 2: extor fixdir loops\n") ;
00726       exit(0) ;
00727    }
00728 
00729    fixdir = strtol(argv[1],NULL,10) ;
00730    if( fixdir < 0 ){ use_tmask = 1 ; fixdir = -fixdir ; }
00731    if( fixdir<1 || fixdir>3 ){fprintf(stderr,"fixdir=%d?\n",fixdir);exit(1);}
00732 
00733    if( argc == 3 ){
00734       float ct[5] ;
00735       ploop = strtol(argv[2],NULL,10) ;
00736       if( ploop < 1 ){ fprintf(stderr,"loop=%d?\n",ploop);exit(1);}
00737       extract_byte_speedtest( ploop , fixdir , ct ) ;
00738       printf("Speed test with fixdir=%d\n"
00739              "_nn   = %g (%g/rep)\n"
00740              "_lifl = %g (%g/rep)\n"
00741              "_liby = %g (%g/rep)\n"
00742              "_ts   = %g (%g/rep)\n"
00743              "_fs   = %g (%g/rep)\n" ,
00744              fixdir ,
00745              ct[0],ct[0]/ploop, ct[1],ct[1]/ploop,
00746              ct[2],ct[2]/ploop, ct[3],ct[3]/ploop, ct[4],ct[4]/ploop ) ;
00747       exit(1) ;
00748    }
00749 
00750    aa = strtod(argv[2],NULL) ;
00751    bb = strtod(argv[3],NULL) ;
00752    if( aa == 0.0 && bb == 0.0 ){fprintf(stderr,"A=B=0?\n");exit(1);}
00753 
00754    if( argc > 5 ){
00755       ploop = strtol(argv[5],NULL,10) ;
00756       if( ploop < 1 ){ fprintf(stderr,"loop=%d?\n",ploop);exit(1); }
00757    }
00758 
00759    if( argc > 6 ){
00760       cfun = argv[6] ;
00761       if( strstr(argv[6],"nn") != NULL )
00762          func = extract_byte_nn ;
00763       else if( strstr(argv[6],"lifl") != NULL )
00764          func = extract_byte_lifl ;
00765       else if( strstr(argv[6],"liby") != NULL )
00766          func = extract_byte_liby ;
00767       else if( strstr(argv[6],"ts") != NULL )
00768          func = extract_byte_ts ;
00769       else if( strstr(argv[6],"fs") != NULL )
00770          func = extract_byte_fs ;
00771       else {
00772          fprintf(stderr,"Unknown func suffix\n");exit(1);
00773       }
00774    }
00775 
00776    in_dset = THD_open_dataset( argv[4] ) ;
00777    if( in_dset == NULL ){fprintf(stderr,"can't open dataset?\n");exit(1);}
00778    if( DSET_NVALS(in_dset) > 1 ){fprintf(stderr,"nvals > 1?\n");exit(1);}
00779    if( DSET_BRICK_TYPE(in_dset,0) != MRI_byte ){fprintf(stderr,"not byte?\n");exit(1);}
00780 
00781    nx = DSET_NX(in_dset) ;
00782    ny = DSET_NY(in_dset) ;
00783    nz = DSET_NZ(in_dset) ; nxy = nx*ny ;
00784 
00785    ASSIGN_DIRECTIONS ;
00786 
00787    da = fabs( 0.5*aa*(nc-1.0) ) ; db = fabs( 0.5*bb*(nc-1.0) ) ;
00788    if( da < 1.0 && db < 1.0 ){fprintf(stderr,"da=%g db=%g ?\n",da,db);exit(1);}
00789 
00790    apad = (int)(2.0+da) ; bpad = (int)(2.0+db) ;
00791    ma   = na + 2*apad   ; mb   = nb + 2*bpad   ; mab = ma*mb ;
00792 
00793    DSET_load(in_dset) ;
00794    vin = DSET_BRICK_ARRAY(in_dset,0) ;
00795 
00796    imout = mri_new( ma,mb , MRI_byte ) ; vout = MRI_BYTE_PTR(imout) ;
00797    immax = mri_new( ma,mb , MRI_byte ) ; vmax = MRI_BYTE_PTR(immax) ;
00798 
00799    tmask = (use_tmask) ? create_Tmask(nx,ny,nz,vin) : NULL ;
00800 
00801    cputim = COX_cpu_time() ;
00802 
00803    for( pp=0 ; pp < ploop ; pp++ ){
00804      memset( vmax , 0 , mab ) ;
00805      for( kk=0 ; kk < nc ; kk++ ){
00806         da = aa*(kk - 0.5*(nc-1.0)) + apad ;
00807         db = bb*(kk - 0.5*(nc-1.0)) + bpad ;
00808 
00809         func( nx,ny,nz , vin,tmask , fixdir,kk , da,db , ma,mb , vout ) ;
00810 
00811         for( ii=0 ; ii < mab ; ii++ )
00812           if( vout[ii] > vmax[ii] ) vmax[ii] = vout[ii] ;
00813      }
00814    }
00815 
00816    cputim = (COX_cpu_time() - cputim)/ploop ;
00817    fprintf(stderr,"CPU time per loop = %g [%s]\n",cputim,cfun) ;
00818 
00819    { char fname[128] = "exim_" ;
00820      strcat(fname,cfun) ; strcat(fname,".pgm") ;
00821      mri_write( fname , immax ) ;
00822    }
00823 
00824    exit(0) ;
00825 }