00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #include "mrilib.h"
00008 
00009 /*===========================================================================
00010    Function to extract a plane of shifted bytes from a 3D volume.
00011      nx, ny, nz = dimensions of vol
00012      vol        = input 3D volume of bytes
00013 
00014      fixdir = fixed direction (1=x, 2=y, 3=z)
00015      fixijk = fixed index
00016      da, db = shift in planar coordinaes (non-fixed directions)
00017      ma, mb = dimensions of im
00018      im     = output 2D image
00019 
00020    Goal is im[a,b] = vol[ P(a-da,b-db,c=fixijk) ] for a=0..ma-1, b=0..mb-1,
00021    where P(a,b,c) is the permutation of (a,b,c) that goes with fixdir:
00022      P(x,y,z) = (y,z,x) for fixdir == 1
00023      P(x,y,z) = (z,x,y) for fixdir == 2
00024      P(x,y,z) = (x,y,z) for fixdir == 3
00025    For values outside the range of vol[], im[] is set to 0.
00026 =============================================================================*/
00027 
00028   /* macros for offsets in vol[] to corners of the interpolation square */
00029 
00030 #undef LL
00031 #undef LR
00032 #undef UL
00033 #undef UR
00034 
00035 #define LL 0                /* voxel offset to lower left  */
00036 #define LR astep            /* voxel offset to lower right */
00037 #define UL bstep            /* voxel offset to upper left  */
00038 #define UR (astep+bstep)    /* voxel offset to upper right */
00039 
00040 #define ASSIGN_DIRECTIONS                                       \
00041  do{ switch( fixdir ){                                          \
00042       default:                                                  \
00043       case 1:            /* x-direction: (a,b,c) = (y,z,x) */   \
00044          astep = nx ; bstep = nxy ; cstep = 1  ;                \
00045          na    = ny ; nb    = nz  ; nc    = nx ;                \
00046       break ;                                                   \
00047                                                                 \
00048       case 2:            /* y-direction: (a,b,c) = (z,x,y) */   \
00049          astep = nxy ; bstep = 1  ; cstep = nx ;                \
00050          na    = nz  ; nb    = nx ; nc    = ny ;                \
00051       break ;                                                   \
00052                                                                 \
00053       case 3:            /* z-direction: (a,b,c) = (x,y,z) */   \
00054          astep = 1  ; bstep = nx ; cstep = nxy ;                \
00055          na    = nx ; nb    = ny ; nc    = nz  ;                \
00056       break ;                                                   \
00057     } } while(0)
00058 
00059 /*-----------------------------------------------------------------------*/
00060 
00061 void extract_assign_directions( int nx, int ny, int nz, int fixdir ,
00062                                 int *Astep, int *Bstep, int *Cstep ,
00063                                 int *Na   , int *Nb   , int *Nc     )
00064 {
00065    int astep,bstep,cstep , na,nb,nc , nxy=nx*ny ;
00066 
00067    ASSIGN_DIRECTIONS ;
00068 
00069    *Astep = astep ; *Bstep = bstep ; *Cstep = cstep ;
00070    *Na    = na    ; *Nb    = nb    ; *Nc    = nc    ; return ;
00071 }
00072 
00073 /*-----------------------------------------------------------------------
00074    NN "interpolation"
00075 -------------------------------------------------------------------------*/
00076 
00077 void extract_byte_nn( int nx , int ny , int nz , byte * vol ,
00078                       Tmask * tm ,
00079                       int fixdir , int fixijk , float da , float db ,
00080                       int ma , int mb , byte * im )
00081 {
00082    int adel,bdel , abot,atop , bb,bbot,btop , nxy=nx*ny ;
00083    register int aa , ijkoff , aoff,boff ;
00084    int astep,bstep,cstep , na,nb,nc ;
00085    byte * mask ;
00086 
00087    memset( im , 0 , ma*mb ) ;  /* initialize output to zero */
00088 
00089    if( fixijk < 0 ) return ;
00090 
00091    ASSIGN_DIRECTIONS ;
00092 
00093    if( fixijk >= nc ) return ;
00094 
00095    da += 0.5 ; adel = (int) da ; if( da < 0.0 ) adel-- ;  /* floor(da+0.5) */
00096    db += 0.5 ; bdel = (int) db ; if( db < 0.0 ) bdel-- ;  /* floor(db+0.5) */
00097 
00098    abot = 0       ; if( abot < adel ) abot = adel ;       /* range in im[] */
00099    atop = na+adel ; if( atop > ma   ) atop = ma ;
00100 
00101    bbot = 0       ; if( bbot < bdel ) bbot = bdel ;
00102    btop = nb+bdel ; if( btop > mb   ) btop = mb ;
00103 
00104    ijkoff = fixijk*cstep + (abot-adel)*astep + (bbot-bdel)*bstep ;
00105    boff   = bbot * ma ;
00106 
00107    mask = (tm == NULL) ? NULL
00108                        : tm->mask[fixdir%3] + (fixijk*nb - bdel) ;
00109 
00110    for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
00111       if( mask == NULL || mask[bb] )
00112          for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
00113             im[aa+boff] = vol[aoff+ijkoff] ;  /* im(aa,bb) = vol(aa-adel,bb-bdel,fixijk) */
00114                                               /*           = vol[ (aa-adel)*astep +
00115                                                                   (bb-bdel)*bstep +
00116                                                                   fixijk   *cstep   ]    */
00117 
00118    return ;
00119 }
00120 
00121 /*---------------------------------------------------------------------------
00122     Two-step interpolation
00123 -----------------------------------------------------------------------------*/
00124 
00125 #if 1
00126 # define TSBOT 0.3
00127 # define TSTOP 0.7
00128 #else
00129 # define TSBOT 0.25
00130 # define TSTOP 0.75
00131 #endif
00132 
00133 void extract_byte_ts( int nx , int ny , int nz , byte * vol ,
00134                       Tmask * tm ,
00135                       int fixdir , int fixijk , float da , float db ,
00136                       int ma , int mb , byte * im )
00137 {
00138    int adel,bdel , abot,atop , bb,bbot,btop , nxy=nx*ny ;
00139    register int aa , ijkoff , aoff,boff ;
00140    int astep,bstep,cstep , na,nb,nc , nts,dts1,dts2 ;
00141    float fa , fb ;
00142    byte * mask ;
00143 
00144    memset( im , 0 , ma*mb ) ;  /* initialize output to zero */
00145 
00146    if( fixijk < 0 ) return ;
00147 
00148    ASSIGN_DIRECTIONS ;
00149 
00150    if( fixijk >= nc ) return ;
00151 
00152    adel = (int) da ; if( da < 0.0 ) adel-- ;  /* floor(da) */
00153    bdel = (int) db ; if( db < 0.0 ) bdel-- ;  /* floor(db) */
00154 
00155    fa = da - adel ;               /* fractional part of dj */
00156    fb = db - bdel ;               /* fractional part of dk */
00157 
00158    fa = 1.0-fa ; fb = 1.0-fb ;
00159 
00160    if( fa < TSBOT ){                      /*- Left 30% -*/
00161       if( fb < TSBOT ){                   /*- Lower 30% -*/
00162         nts = 1 ; dts1 = LL ;               /* [0,0] */
00163       } else if( fb > TSTOP ){            /*- Upper 30% -*/
00164         nts = 1 ; dts1 = UL ;               /* [0,1] */
00165       } else {                            /*- Middle 40% -*/
00166         nts = 2 ; dts1 = LL ; dts2 = UL ;   /* mid of [0,0] and [0,1] */
00167       }
00168    } else if( fa > TSTOP ){               /*- Right 30% -*/
00169       if( fb < TSBOT ){                   /*- Lower 30% -*/
00170         nts = 1 ; dts1 = LR ;               /* [1,0] */
00171       } else if( fb > TSTOP ){            /*- Upper 30% -*/
00172         nts = 1 ; dts1 = UR ;               /* [1,1] */
00173       } else {                            /*- Middle 40% -*/
00174         nts = 2 ; dts1 = LR ; dts2 = UR ;   /* mid of [1,0] and [1,1] */
00175       }
00176    } else {                               /*- Middle 40% -*/
00177       if( fb < TSBOT ){                   /*- Lower 30% -*/
00178         nts = 2 ; dts1 = LL ; dts2 = LR ;   /* mid of [0,0] and [1,0] */
00179       } else if( fb > TSTOP ){            /*- Upper 30% -*/
00180         nts = 2 ; dts1 = UL ; dts2 = UR ;   /* mid of [0,1] and [1,1] */
00181       } else {                            /*- Middle 40% -*/
00182         nts = 4 ;                           /* mid of all 4 points */
00183       }
00184    }
00185 
00186    adel++ ; bdel++ ;
00187 
00188    abot = 0         ; if( abot < adel ) abot = adel ;       /* range in im[] */
00189    atop = na+adel-1 ; if( atop > ma   ) atop = ma ;
00190 
00191    bbot = 0         ; if( bbot < bdel ) bbot = bdel ;
00192    btop = nb+bdel-1 ; if( btop > mb   ) btop = mb ;
00193 
00194    ijkoff = fixijk*cstep + (abot-adel)*astep + (bbot-bdel)*bstep ;
00195    boff   = bbot * ma ;
00196 
00197    mask = (tm == NULL) ? NULL
00198                        : tm->mask[fixdir%3] + (fixijk*nb - bdel) ;
00199 
00200    switch( nts ){
00201 
00202       case 1:
00203          ijkoff += dts1 ;
00204          for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
00205            if( mask == NULL || mask[bb] || mask[bb+1] )
00206              for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
00207                im[aa+boff] = vol[aoff+ijkoff] ;
00208       break ;
00209 
00210       case 2:
00211          ijkoff += dts1 ; dts2 -= dts1 ;
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] = (vol[aoff+ijkoff] + vol[aoff+(ijkoff+dts2)]) >> 1;
00216       break ;
00217 
00218       case 4:
00219          for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
00220            if( mask == NULL || mask[bb] || mask[bb+1] )
00221              for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
00222                im[aa+boff] = ( vol[aoff+ijkoff]     +vol[aoff+(ijkoff+LR)]
00223                               +vol[aoff+(ijkoff+UL)]+vol[aoff+(ijkoff+UR)]) >> 2;
00224       break ;
00225    }
00226 
00227    return ;
00228 }
00229 
00230 /*----------------------------------------------------------------------------*/
00231 
00232 #undef U
00233 #define U(i,j) uu.mat[i][j]
00234 
00235 MRI_IMAGE * project_byte_mip( int nx, int ny, int nz, byte * vol, Tmask * tm,
00236                               THD_mat33 uu )
00237 {
00238    int ii,jj,kk , ni,nj,nk,pk , ma,mb,mab , pij , nnn[3] ;
00239    float utop,uabs , a,b , aii,aij,aji,ajj , hnk , ba,bb ;
00240    byte *im , *sl ;
00241    MRI_IMAGE *bim , *qim ;
00242 
00243 #if 1
00244 DUMP_MAT33("rotation",uu) ;
00245 #endif
00246 
00247    /*-- find element U(kk,2) that is largest --*/
00248 
00249    nnn[0] = nx ; nnn[1] = ny ; nnn[2] = nz ;
00250 
00251    kk = 0 ; utop = fabs(U(0,2)) ;
00252    uabs = fabs(U(1,2)) ; if( uabs > utop ){ utop = uabs; kk = 1; }
00253    uabs = fabs(U(2,2)) ; if( uabs > utop ){ utop = uabs; kk = 2; }
00254 
00255    if( utop == 0.0 ) return ;   /* bad matrix */
00256 
00257    ii = (kk+1) % 3 ;  /* image axes */
00258    jj = (kk+2) % 3 ;
00259 
00260    a = U(ii,2) / U(kk,2) ;  /* shearing parameters */
00261    b = U(jj,2) / U(kk,2) ;
00262 
00263 #if 0
00264 fprintf(stderr,"kk=%d a=%g b=%g\n",kk,a,b) ;
00265 #endif
00266 
00267    aii = U(ii,0) - a * U(kk,0) ;  /* warping parameters */
00268    aij = U(ii,1) - a * U(kk,1) ;  /* [not used just yet] */
00269    aji = U(jj,0) - b * U(kk,0) ;
00270    ajj = U(jj,1) - b * U(kk,1) ;
00271 
00272 #if 0
00273 fprintf(stderr,"warp: aii=%g  aij=%g\n"
00274                "      aji=%g  ajj=%g\n" , aii,aij,aji,ajj ) ;
00275 #endif
00276 
00277    ni  = nnn[ii] ; nj = nnn[jj] ; nk = nnn[kk] ; hnk = 0.5*nk ;
00278    ma  = MAX(ni,nj) ; ma = MAX(ma,nk) ; ma *= 1.2 ;
00279    mb  = ma ; mab = ma * mb ; ba = 0.5*(ma-ni) ; bb = 0.5*(mb-nj) ;
00280    sl  = (byte *) malloc(mab) ;
00281    bim = mri_new(ma,mb,MRI_byte) ; im = MRI_BYTE_PTR(bim) ; memset(im,0,mab) ;
00282    for( pk=0 ; pk < nk ; pk++ ){
00283       extract_byte_ts( nx,ny,nz , vol , tm ,
00284                        kk+1 , pk , ba-a*(pk-hnk) , bb-b*(pk-hnk) ,
00285                        ma,mb , sl ) ;
00286       for( pij=0 ; pij < mab ; pij++ )
00287          if( sl[pij] > im[pij] ) im[pij] = sl[pij] ;
00288    }
00289 
00290    free(sl) ;
00291 
00292 #if 1
00293    qim = mri_aff2d_byte( bim , 1 , aii,aij,aji,ajj ) ;
00294    mri_free(bim) ; bim = qim ;
00295 #endif
00296 
00297    return bim ;
00298 }
00299 
00300 /*$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$*/
00301 
00302 static THD_mat33 rotmatrix( int ax1,float th1 ,
00303                             int ax2,float th2 , int ax3,float th3  )
00304 {
00305    THD_mat33 q , p ;
00306 
00307    LOAD_ROT_MAT( q , th1 , ax1 ) ;
00308    LOAD_ROT_MAT( p , th2 , ax2 ) ; q = MAT_MUL( p , q ) ;
00309    LOAD_ROT_MAT( p , th3 , ax3 ) ; q = MAT_MUL( p , q ) ;
00310 
00311    return q ;
00312 }
00313 
00314 int main( int argc , char * argv[] )
00315 {
00316    THD_3dim_dataset *dset ;
00317    int iarg=1 ;
00318    char *cc1="x",*cc2="y",*cc3="z" ;
00319    float th1=0.0, th2=0.0, th3=0.0 ;
00320    float thx,thy,thz ;
00321    int   axx,ayy,azz ;
00322    char *fname="testcox.ppm" ;
00323    void * rhand ;
00324    int bot=1 , ii ;
00325    float omap[128] , bfac ;
00326    MRI_IMAGE * im , * brim ;
00327    int hbr[256] , nperc,ibot,itop,sum ;
00328    byte * bar ;
00329    THD_mat33 rmat ;
00330 
00331    if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00332       printf("Usage: testcox [-rotate a b c] [-out f] [-bot b] dset\n") ;
00333       exit(0) ;
00334    }
00335 
00336    while( iarg < argc && argv[iarg][0] == '-' ){
00337 
00338       if( strcmp(argv[iarg],"-bot") == 0 ){
00339          bot = strtod( argv[++iarg] , NULL ) ;
00340          iarg++ ; continue ;
00341       }
00342 
00343       if( strcmp(argv[iarg],"-rotate") == 0 ){
00344          th1 = (PI/180.0) * strtod( argv[++iarg] , &cc1 ) ;
00345          th2 = (PI/180.0) * strtod( argv[++iarg] , &cc2 ) ;
00346          th3 = (PI/180.0) * strtod( argv[++iarg] , &cc3 ) ;
00347 
00348          iarg++ ; continue ;
00349       }
00350 
00351       if( strcmp(argv[iarg],"-out") == 0 ){
00352          fname = argv[++iarg] ;
00353          iarg++ ; continue ;
00354       }
00355 
00356       fprintf(stderr,"Illegal option: %s\n",argv[iarg]); exit(1);
00357    }
00358 
00359    if( iarg >= argc ){fprintf(stderr,"No dataset?\n"); exit(1); }
00360 
00361    dset = THD_open_dataset( argv[iarg] ) ;
00362    if( dset == NULL ){fprintf(stderr,"Can't open dataset!\n");exit(1);}
00363    if( DSET_BRICK_TYPE(dset,0) != MRI_byte ){
00364       fprintf(stderr,"Non-byte dataset input!\n");exit(1);
00365    }
00366    DSET_mallocize(dset) ; DSET_load(dset) ;
00367    if( !DSET_LOADED(dset) ){
00368       fprintf(stderr,"Can't load dataset!\n");exit(1);
00369    }
00370 
00371    /* correct angles to go with rendering plugin are
00372          <yaw>A <pitch>R <roll>I
00373       in that order                                 */
00374 
00375    THD_rotangle_user_to_dset( dset ,
00376                               th1,*cc1  , th2,*cc2  , th3,*cc3 ,
00377                               &thx,&axx , &thy,&ayy , &thz,&azz ) ;
00378    rmat = rotmatrix( axx,thx , ayy,thy , azz,thz ) ;
00379 
00380    im = project_byte_mip( DSET_NX(dset) , DSET_NY(dset) , DSET_NZ(dset) ,
00381                           DSET_ARRAY(dset,0) , NULL , rmat ) ;
00382 
00383    mri_write_pnm( fname , im ) ;
00384    fprintf(stderr,"+++ Output to file %s\n",fname);
00385    exit(0) ;
00386 }