Skip to content

AFNI/NIfTI Server

Sections
Personal tools
You are here: Home » AFNI » Documentation
 

Doxygen Source Code Documentation

Main Page   Alphabetical List   Data Structures   File List   Data Fields   Globals   Search  

cox_render.c

Go to the documentation of this file.
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 "cox_render.h"
00008 
00009 #undef CREN_DEBUG
00010 
00011 /*============================================================================
00012   CREN = Cox Renderer, a set of routines for volume rendering 3D bricks.
00013 
00014   (1) Use new_CREN_renderer() to create a renderer.
00015   (2) Load grayscale+color bricks into it with function CREN_set_databytes()
00016   (3) Set the viewing angles with CREN_set_viewpoint().
00017       (Optional: set axes ordering with CREN_set_axes().)
00018   (4) Create an image with CREN_render().
00019       Loop back through (2-4) as needed to create new images.
00020   (5) When finished, use destroy_CREN_renderer() to clean up.
00021 ==============================================================================*/
00022 
00023 static int num_renderers = 0 ;  /* global count of how many are open */
00024 
00025 THD_mat33 rotmatrix( int ax1,float th1 ,
00026                      int ax2,float th2 , int ax3,float th3  ) ;
00027 
00028 void extract_byte_nn( int nx , int ny , int nz , byte * vol ,
00029                       Tmask * tm ,
00030                       int fixdir , int fixijk , float da , float db ,
00031                       int ma , int mb , byte * im ) ;
00032 
00033 void extract_rgba_nn( int nx , int ny , int nz , rgba * vol ,
00034                       Tmask * tm ,
00035                       int fixdir , int fixijk , float da , float db ,
00036                       int ma , int mb , rgba * im ) ;
00037 
00038 void extract_byte_tsx( int nx , int ny , int nz , byte * vol ,
00039                        Tmask * tm ,
00040                        int fixdir , int fixijk , float da , float db ,
00041                        int ma , int mb , byte * im ) ;
00042 
00043 void extract_byte_lix( int nx , int ny , int nz , byte * vol ,
00044                        Tmask * tm ,
00045                        int fixdir , int fixijk , float da , float db ,
00046                        int ma , int mb , byte * im ) ;
00047 
00048 void extract_byte_lixx( int nx , int ny , int nz , byte * vol ,
00049                         Tmask * tm ,
00050                         int fixdir , int fixijk , float da , float db ,
00051                         int ma , int mb , byte * im ) ;
00052 
00053 typedef void exfunc(int,int,int,byte*,Tmask*,int,int,float,float,int,int,byte*);
00054 
00055 #if 0
00056 static void extract_byte_ts( int nx , int ny , int nz , byte * vol ,
00057                              Tmask * tm ,
00058                              int fixdir , int fixijk , float da , float db ,
00059                              int ma , int mb , byte * im ) ;
00060 #endif
00061 
00062 /*--------------------------------------------------------------------------
00063   Create a new renderer.
00064   The return pointer is passed to all other CREN routines.
00065 ----------------------------------------------------------------------------*/
00066 
00067 void * new_CREN_renderer( void )
00068 {
00069    CREN_stuff * ar ;
00070    int ii ;
00071 
00072    /*-- make storage for rendering struct --*/
00073 
00074    ar = (CREN_stuff *) malloc( sizeof(CREN_stuff) ) ;
00075    ar->type = CREN_TYPE ;
00076 
00077    /*-- initialize rendering struct to somewhat random values --*/
00078 
00079    ar->nx = ar->ny = ar->nz = ar->newvox = 0 ;
00080    ar->dx = ar->dy = ar->dz = 1.0 ;
00081 
00082    /* default axis rotation order is as in plug_render.c */
00083 
00084    ar->ax1 = 1   ; ar->ax2 = 0   ; ar->ax3 = 2   ;
00085    ar->th1 = 0.0 ; ar->th2 = 0.0 ; ar->th3 = 0.0 ; ar->newangles = 1 ;
00086 
00087    ar->vox = NULL ;  /* no data yet */
00088    ar->vtm = NULL ;  /* no Tmask yet */
00089 
00090    ar->vox_is_gray = 0 ;
00091 
00092    ar->newopa = 0 ;
00093    ar->opargb = 1.0 ;             /* colored voxels are opaque */
00094    for( ii=0 ; ii < 128 ; ii++ )  /* linear map for gray opacity */
00095       ar->opamap[ii] = ii/127.0 ;
00096 
00097    ar->nrgb   = 0 ;               /* no color map set */
00098    memset( ar->rmap , 0 , 128 ) ; memset( ar->gmap , 0 , 128 ) ;
00099    memset( ar->bmap , 0 , 128 ) ; memset( ar->imap , 0 , 128 ) ;
00100 
00101    ar->min_opacity = 0.05 ;
00102    ar->renmode     = CREN_SUM_VOX ;
00103    ar->intmode     = CREN_TWOSTEP ;
00104 
00105    LOAD_DIAG_MAT( ar->skewmat , 1.0,1.0,1.0 ) ;
00106 
00107    num_renderers ++ ; return (void *) ar ;
00108 }
00109 
00110 /*-----------------------------------------------------------------------------
00111    Get rid of a renderer and all its contents
00112 -------------------------------------------------------------------------------*/
00113 
00114 void destroy_CREN_renderer( void * ah )
00115 {
00116    CREN_stuff * ar = (CREN_stuff *) ah ;
00117 
00118    if( !ISVALID_CREN(ar) ) return ;
00119 
00120    if( ar->vox != NULL ) free(ar->vox) ;
00121    if( ar->vtm != NULL ) free_Tmask(ar->vtm) ;
00122    free(ar) ;
00123 
00124    num_renderers -- ; return ;
00125 }
00126 
00127 
00128 /*-----------------------------------------------------------------------------
00129    Set the matrix to skew ijk indices to xyz coordinates (not normally needed)
00130 -------------------------------------------------------------------------------*/
00131 
00132 void CREN_set_skewmat( void * ah , THD_mat33 sm )
00133 {
00134    CREN_stuff * ar = (CREN_stuff *) ah ;
00135 
00136    if( !ISVALID_CREN(ar) ) return ;
00137 
00138    ar->skewmat = sm ; return ;
00139 }
00140 
00141 /*-----------------------------------------------------------------------------
00142    Set the minimum voxel opacity to consider
00143 -------------------------------------------------------------------------------*/
00144 
00145 void CREN_set_min_opacity( void * ah , float opm )
00146 {
00147    CREN_stuff * ar = (CREN_stuff *) ah ;
00148 
00149    if( !ISVALID_CREN(ar) ) return ;
00150    if( opm <= 0.0 || opm >= 1.0 ) opm = 0.05 ;
00151    ar->min_opacity = opm ;
00152    return ;
00153 }
00154 
00155 /*-----------------------------------------------------------------------------
00156    Set the rendering mode
00157      CREN_SUM_VOX = integral of voxel data times opacity
00158      CREN_MIP_VOX = maximum voxel intensity
00159      CREN_MIP_OPA = maximum opacity
00160 -------------------------------------------------------------------------------*/
00161 
00162 void CREN_set_render_mode( void * ah , int mmm )
00163 {
00164    CREN_stuff * ar = (CREN_stuff *) ah ;
00165    if( !ISVALID_CREN(ar) ) return ;
00166 
00167    if( mmm < 0 || mmm > CREN_LAST_MODE ) mmm = CREN_SUM_VOX ;
00168    ar->renmode = mmm ;
00169    return ;
00170 }
00171 
00172 /*-----------------------------------------------------------------------------
00173    Set the interpolation mode
00174      CREN_NN      = nearest neighbor
00175      CREN_TWOSTEP = two-step
00176      CREN_LINEAR  = linear
00177 -------------------------------------------------------------------------------*/
00178 
00179 void CREN_set_interp( void * ah , int mmm )
00180 {
00181    CREN_stuff * ar = (CREN_stuff *) ah ;
00182    if( !ISVALID_CREN(ar) ) return ;
00183 
00184    if( mmm < 0 || mmm > CREN_LINEAR ) mmm = CREN_TWOSTEP ;
00185    ar->intmode = mmm ;
00186    return ;
00187 }
00188 
00189 /*-----------------------------------------------------------------------------
00190   Load the user defined colormap.
00191     ncol    = number of colors (must be from 1 to 128, inclusive)
00192     rmap[i] = red   map for index i=0..ncol-1 (for voxel values i+128)
00193     gmap[i] = green map for index i=0..ncol-1 (for voxel values i+128)
00194     bmap[i] = blue  map for index i=0..ncol-1 (for voxel values i+128)
00195 -------------------------------------------------------------------------------*/
00196 
00197 void CREN_set_rgbmap( void *ah, int ncol, byte *rmap, byte *gmap, byte *bmap )
00198 {
00199    CREN_stuff * ar = (CREN_stuff *) ah ;
00200    int ii ;
00201 
00202    if( !ISVALID_CREN(ar) ) return ;
00203    if( ncol<1 || ncol>128 || rmap==NULL || gmap==NULL || bmap==NULL ) return ;
00204 
00205    ar->nrgb = ncol ;
00206 
00207    /* copy into rendering struct, and compute intensity of each color */
00208 
00209    for( ii=0 ; ii < ncol ; ii++ ){
00210       ar->rmap[ii] = rmap[ii] ;
00211       ar->gmap[ii] = gmap[ii] ;
00212       ar->bmap[ii] = bmap[ii] ;
00213       ar->imap[ii] = (byte)(0.299*rmap[ii]+0.587*gmap[ii]+0.114*bmap[ii]) ;
00214    }
00215 
00216    /* set leftovers to 0 */
00217 
00218    for( ii=ncol ; ii < 128 ; ii++ )
00219       ar->rmap[ii] = ar->gmap[ii] = ar->bmap[ii] = ar->imap[ii] = 0 ;
00220 
00221    return ;
00222 }
00223 
00224 /*----------------------------------------------------------------------
00225   Set the mapping from voxel value to opacity:
00226     opm[vv] = opacity of voxel value vv, for vv=0..127
00227     opgrb   = (fixed) opacity of colored voxels, with vv=128..255
00228   Opacities range from 0.0 to 1.0
00229 ------------------------------------------------------------------------*/
00230 
00231 void CREN_set_opamap( void *ah, float *opm , float oprgb )
00232 {
00233    CREN_stuff * ar = (CREN_stuff *) ah ;
00234 
00235    if( !ISVALID_CREN(ar) ) return ;
00236 
00237    if( opm != NULL )
00238       memcpy( ar->opamap , opm , sizeof(float)*128 ) ;
00239 
00240    if( oprgb >= 0.0 && oprgb <= 1.0 )
00241       ar->opargb = oprgb ;
00242 
00243 #ifdef CREN_DEBUG
00244 fprintf(stderr,"CREN_set_opamap: opm[0]=%g opm[127]=%g oprgb=%g\n",
00245 opm[0],opm[127],oprgb) ;
00246 #endif
00247 
00248    ar->newopa = 1 ; return ;
00249 }
00250 
00251 /*-----------------------------------------------------------------------
00252   See if the thing needs data bytes
00253 -------------------------------------------------------------------------*/
00254 
00255 int CREN_needs_data( void * ah )
00256 {
00257    CREN_stuff * ar = (CREN_stuff *) ah ;
00258 
00259    if( !ISVALID_CREN(ar) ) return -1 ;
00260 
00261    return (ar->vox == NULL) ;
00262 }
00263 
00264 /*-----------------------------------------------------------------------------
00265   Set the data axes stuff.  The volume has 3 axes, indexed by (i,j,k).
00266   aii is set to denote the spatial direction of the increasing i direction:
00267       1 == +x direction
00268      -1 == -x direction
00269       2 == +y direction
00270      -2 == -y direction
00271       3 == +z direction
00272      -3 == -z direction
00273   and similarly for ajj, akk for the j and k directions.
00274   di, dj, dk are the magnitude of the stepsizes in each direction.
00275   (The number of points along each axis is set via CREN_set_databytes.)
00276 -------------------------------------------------------------------------------*/
00277 
00278 void CREN_set_axes( void * ah , int  aii , int  ajj , int  akk ,
00279                                 float di , float dj , float dk  )
00280 {
00281    CREN_stuff * ar = (CREN_stuff *) ah ;
00282    int abii=abs(aii) , abjj=abs(ajj) , abkk=abs(akk) ;
00283 
00284    /*-- sanity checks --*/
00285 
00286    if( !ISVALID_CREN(ar) ) return ;
00287 
00288    if( abii < 1 || abii > 3 ||
00289        abjj < 1 || abjj > 3 ||
00290        abkk < 1 || abkk > 3 || abii+abjj+abkk != 6 ) return ;
00291 
00292    /*-- load stepsizes --*/
00293 
00294    ar->dx = fabs(di) ; if( ar->dx == 0.0 ) ar->dx = 1.0 ;
00295    ar->dy = fabs(dj) ; if( ar->dy == 0.0 ) ar->dy = 1.0 ;
00296    ar->dz = fabs(dk) ; if( ar->dz == 0.0 ) ar->dz = 1.0 ;
00297 
00298    /*-- construct skewmat --*/
00299 
00300    LOAD_ZERO_MAT(ar->skewmat) ;
00301 
00302 #if 0
00303    ar->skewmat.mat[abii-1][0] = (aii > 0) ? 1.0 : -1.0 ;
00304    ar->skewmat.mat[abjj-1][1] = (ajj > 0) ? 1.0 : -1.0 ;
00305    ar->skewmat.mat[abkk-1][2] = (akk > 0) ? 1.0 : -1.0 ;
00306 #else
00307    ar->skewmat.mat[0][abii-1] = (aii > 0) ? 1.0 : -1.0 ;
00308    ar->skewmat.mat[1][abjj-1] = (ajj > 0) ? 1.0 : -1.0 ;
00309    ar->skewmat.mat[2][abkk-1] = (akk > 0) ? 1.0 : -1.0 ;
00310 #endif
00311 
00312 #ifdef CREN_DEBUG
00313 DUMP_MAT33("skewmat",ar->skewmat) ;
00314 #endif
00315 
00316    ar->newangles = 1 ; return ;
00317 }
00318 
00319 /*---------------------------------------------------------------*/
00320 
00321 void CREN_dset_axes( void * ah , THD_3dim_dataset * dset )
00322 {
00323    CREN_stuff * ar = (CREN_stuff *) ah ;
00324    int aii=1 , ajj=2 , akk=3 ;
00325    float dx , dy , dz ;
00326 
00327    if( !ISVALID_CREN(ar) || !ISVALID_DSET(dset) ) return ;
00328 
00329    switch( dset->daxes->xxorient ){
00330       case ORI_R2L_TYPE: aii =  1 ; break ;
00331       case ORI_L2R_TYPE: aii = -1 ; break ;
00332       case ORI_A2P_TYPE: aii =  2 ; break ;
00333       case ORI_P2A_TYPE: aii = -2 ; break ;
00334       case ORI_I2S_TYPE: aii =  3 ; break ;
00335       case ORI_S2I_TYPE: aii = -3 ; break ;
00336       default:
00337         fprintf(stderr,"** CREN_dset_axes: illegal xxorient code!\n") ;
00338    }
00339    switch( dset->daxes->yyorient ){
00340       case ORI_R2L_TYPE: ajj =  1 ; break ;
00341       case ORI_L2R_TYPE: ajj = -1 ; break ;
00342       case ORI_A2P_TYPE: ajj =  2 ; break ;
00343       case ORI_P2A_TYPE: ajj = -2 ; break ;
00344       case ORI_I2S_TYPE: ajj =  3 ; break ;
00345       case ORI_S2I_TYPE: ajj = -3 ; break ;
00346       default:
00347         fprintf(stderr,"** CREN_dset_axes: illegal yyorient code!\n") ;
00348    }
00349    switch( dset->daxes->zzorient ){
00350       case ORI_R2L_TYPE: akk =  1 ; break ;
00351       case ORI_L2R_TYPE: akk = -1 ; break ;
00352       case ORI_A2P_TYPE: akk =  2 ; break ;
00353       case ORI_P2A_TYPE: akk = -2 ; break ;
00354       case ORI_I2S_TYPE: akk =  3 ; break ;
00355       case ORI_S2I_TYPE: akk = -3 ; break ;
00356       default:
00357         fprintf(stderr,"** CREN_dset_axes: illegal zzorient code!\n") ;
00358    }
00359 
00360    dx = fabs(dset->daxes->xxdel) ;
00361    dy = fabs(dset->daxes->yydel) ;
00362    dz = fabs(dset->daxes->zzdel) ;
00363 
00364    CREN_set_axes( ah , aii,ajj,akk , dx,dy,dz ) ;
00365    return ;
00366 }
00367 
00368 /*-----------------------------------------------------------------------------
00369    Set the data brick in a renderer;
00370    Returns -1 if an error occurs, otherwise returns 0.
00371 
00372    Data bytes:
00373      values vv=0..127   correspond to grayscale 0..254 (vv << 1)
00374      values vv=128..255 correspond to colormap index vv-128.
00375 -------------------------------------------------------------------------------*/
00376 
00377 void CREN_set_databytes( void * ah , int ni,int nj,int nk, byte * grim )
00378 {
00379    CREN_stuff * ar = (CREN_stuff *) ah ;
00380    int nvox , ii ;
00381 
00382    /*-- sanity checks --*/
00383 
00384    if( !ISVALID_CREN(ar) || grim == NULL ) return ;
00385    if( ni < 3 || nj < 3 || nk < 3 ) return ;
00386 
00387    /*-- free old data, if any --*/
00388 
00389    if( ar->vox != NULL ){ free(ar->vox)      ; ar->vox = NULL; }
00390    if( ar->vtm != NULL ){ free_Tmask(ar->vtm); ar->vtm = NULL; }
00391 
00392    /*-- set size of each axis--*/
00393 
00394    ar->nx = ni ; ar->ny = nj ; ar->nz = nk ;
00395 
00396    /*-- signal we have new voxel data --*/
00397 
00398    ar->newvox = 1 ;
00399 
00400    /*-- copy data from grim into internal storage --*/
00401 
00402    nvox = ni * nj * nk ;
00403    ar->vox = (byte *) malloc(nvox) ;
00404    memcpy( ar->vox , grim , nvox ) ;
00405 
00406    for( ii=0 ; ii < nvox && grim[ii] < 128 ; ii++ ) ; /* nada */
00407    ar->vox_is_gray = (ii == nvox) ;
00408 
00409    return ;
00410 }
00411 
00412 /*----------------------------------------------------------------------------
00413    Set the viewpoint of the user
00414      axI = axis about which rotation #I is to take place
00415             0 == x , 1 == y , 2 == z
00416      thI = angle of rotation (radians)
00417 ------------------------------------------------------------------------------*/
00418 
00419 void CREN_set_viewpoint( void * ah , int ax1,float th1 ,
00420                                      int ax2,float th2 , int ax3,float th3  )
00421 {
00422    CREN_stuff * ar = (CREN_stuff *) ah ;
00423 
00424    if( !ISVALID_CREN(ar) ) return ;
00425 
00426    ar->ax1 = ax1 ; ar->ax2 = ax2 ; ar->ax3 = ax3 ;
00427    ar->th1 = th1 ; ar->th2 = th2 ; ar->th3 = th3 ;
00428 
00429 #ifdef CREN_DEBUG
00430 fprintf(stderr,"CREN_set_viewpoint: ax1=%d th1=%g ax2=%d th2=%g ax3=%d th3=%g\n",
00431 ax1,th1,ax2,th2,ax3,th3) ;
00432 #endif
00433 
00434    ar->newangles = 1 ; return ;
00435 }
00436 
00437 /*-----------------------------------------------------------------------------*/
00438 
00439 void CREN_set_rotaxes( void * ah , int ax1 , int ax2 , int ax3 )
00440 {
00441    CREN_stuff * ar = (CREN_stuff *) ah ;
00442 
00443    if( !ISVALID_CREN(ar) ) return ;
00444    ar->ax1 = ax1 ; ar->ax2 = ax2 ; ar->ax3 = ax3 ;
00445    ar->newangles = 1 ; return ;
00446 }
00447 
00448 /*-----------------------------------------------------------------------------*/
00449 
00450 void CREN_set_angles( void * ah , float th1 , float th2 , float th3 )
00451 {
00452    CREN_stuff * ar = (CREN_stuff *) ah ;
00453 
00454    if( !ISVALID_CREN(ar) ) return ;
00455    ar->th1 = th1 ; ar->th2 = th2 ; ar->th3 = th3 ;
00456    ar->newangles = 1 ; return ;
00457 }
00458 
00459 /*-----------------------------------------------------------------------------
00460    Actually render an image.  Returns NULL if an error occurs.
00461    The image is always in MRI_rgb format.
00462 
00463    If rotm is valid, copy the rotation matrix.     2002.08.28 - rickr
00464 -------------------------------------------------------------------------------*/
00465 
00466 MRI_IMAGE * CREN_render( void * ah, THD_mat33 * rotm )
00467 {
00468    CREN_stuff * ar = (CREN_stuff *) ah ;
00469    MRI_IMAGE * bim , * qim ;
00470    byte * rgb , * var , * sl ;
00471    int nvox , nxy , vtop=128+ar->nrgb ;
00472    float *used=NULL , obot=ar->min_opacity ;
00473    THD_mat33 uu ;
00474    int ii,jj,kk , ni,nj,nk,pk , ma,mb,mab , nnn[3] ;
00475    float utop,uabs , a,b , aii,aij,aji,ajj , hnk , ba,bb ;
00476    int pk_reverse , ppk ;
00477    float dab ;
00478    register int vv , pij , p3 ;
00479    register float *omap , orgb ;
00480    exfunc *ifunc ;
00481 
00482    /*-- sanity checks --*/
00483 
00484    if( !ISVALID_CREN(ar) ) return NULL ;
00485 
00486    var = ar->vox ;
00487    if( var == NULL ) return NULL ;
00488    if( ar->nx < 3 || ar->ny < 3 || ar->nz < 3 ) return NULL ;
00489 
00490 #ifdef CREN_DEBUG
00491 fprintf(stderr,"CREN_render: nx=%d ny=%d nz=%d\n",ar->nx,ar->ny,ar->nz);
00492 #endif
00493 
00494    nxy  = ar->nx * ar->ny ;  /* number of voxels in one plane */
00495    nvox = nxy * ar->nz ;     /* number of voxels in whole volume */
00496 
00497    /*-- make a Tmask to show which 1D rows contain opacity --*/
00498 
00499    omap = ar->opamap ; orgb = ar->opargb ;
00500 
00501 #if 1
00502    if( ar->newvox || ar->newopa || ar->vtm == NULL ){
00503       register byte *tar, qrgb ;
00504 
00505       free_Tmask(ar->vtm) ;
00506       tar  = (byte *) malloc(nvox) ;
00507       qrgb = (orgb >= obot) ;  /* is color opaque? */
00508       for( pij=0 ; pij < nvox ; pij++ ){
00509          vv = var[pij] ; if( vv == 0 ) continue ;
00510          if( vv < 128 ) tar[pij] = (omap[vv] >= obot) ;
00511          else           tar[pij] = qrgb ;
00512       }
00513 
00514       ar->vtm = create_Tmask_byte( ar->nx,ar->ny,ar->nz , tar ) ;
00515       free(tar) ;
00516    }
00517 #endif
00518 
00519    /*-- compute rotation matrix uu --*/
00520 
00521    uu = rotmatrix( ar->ax1,ar->th1 , ar->ax2,ar->th2 , ar->ax3,ar->th3 ) ;
00522 
00523    if ( rotm != NULL ) /* if requested, copy the rotation matrix, angles only */
00524       *rotm = uu;
00525 
00526    dab = MIN(ar->dx,ar->dy) ; dab = MIN(dab,ar->dz) ;  /* output grid spacing */
00527 
00528    /*-- we want
00529            diag{1/dx,1/dy,1/dz} [skewmat] [uu] diag{dab,dab,dab}
00530         which is the matrix which transforms from
00531         image index coords to volume index coords  --*/
00532 
00533    uu = MAT_MUL( ar->skewmat , uu ) ; /* a signed permutation? */
00534 
00535    uu.mat[0][0] *= (dab / ar->dx ) ;  /* multiply columns by dab  */
00536    uu.mat[0][1] *= (dab / ar->dx ) ;  /* rows by {1/dx,1/dy,1/dz} */
00537    uu.mat[0][2] *= (dab / ar->dx ) ;
00538 
00539    uu.mat[1][0] *= (dab / ar->dy ) ;
00540    uu.mat[1][1] *= (dab / ar->dy ) ;
00541    uu.mat[1][2] *= (dab / ar->dy ) ;
00542 
00543    uu.mat[2][0] *= (dab / ar->dz ) ;
00544    uu.mat[2][1] *= (dab / ar->dz ) ;
00545    uu.mat[2][2] *= (dab / ar->dz ) ;
00546 
00547 #ifdef CREN_DEBUG
00548 DUMP_MAT33("uu matrix",uu) ;
00549 #endif
00550 
00551    /*-- find element uu(kk,2) that is largest: kk = projection axis --*/
00552 
00553    nnn[0] = ar->nx ; nnn[1] = ar->ny ; nnn[2] = ar->nz ;
00554 
00555 #undef U
00556 #define U(i,j) uu.mat[i][j]
00557 
00558    kk = 0 ; utop = fabs(U(0,2)) ;
00559    uabs = fabs(U(1,2)) ; if( uabs > utop ){ utop = uabs; kk = 1; }
00560    uabs = fabs(U(2,2)) ; if( uabs > utop ){ utop = uabs; kk = 2; }
00561 
00562    if( utop == 0.0 ) return NULL ;   /* bad matrix */
00563 
00564    ii = (kk+1) % 3 ;  /* image axes */
00565    jj = (kk+2) % 3 ;
00566 
00567    a = U(ii,2) / U(kk,2) ;  /* shearing parameters */
00568    b = U(jj,2) / U(kk,2) ;
00569 
00570 #ifdef CREN_DEBUG
00571 fprintf(stderr,"kk=%d a=%g b=%g\n",kk,a,b) ;
00572 #endif
00573 
00574    aii = U(ii,0) - a * U(kk,0) ;  /* warping parameters   */
00575    aij = U(ii,1) - a * U(kk,1) ;  /* [to make projection] */
00576    aji = U(jj,0) - b * U(kk,0) ;  /* [image look correct] */
00577    ajj = U(jj,1) - b * U(kk,1) ;
00578 
00579 #ifdef CREN_DEBUG
00580 fprintf(stderr,"warp: aii=%g  aij=%g\n"
00581                "      aji=%g  ajj=%g\n" , aii,aij,aji,ajj ) ;
00582 #endif
00583 
00584    /* n{ijk} = dimension of volume along axis {ijk} */
00585 
00586    ni = nnn[ii] ; nj = nnn[jj] ; nk = nnn[kk] ; hnk = 0.5*nk ;
00587 
00588    pk_reverse = ( U(kk,2) < 0.0 ) ;
00589 
00590    /* output image will be ma x mb pixels */
00591 
00592    ma = MAX(ni,nj) ; ma = MAX(ma,nk) ; ma *= 1.2 ;
00593    mb = ma ; mab = ma * mb ; ba = 0.5*(ma-ni) ; bb = 0.5*(mb-nj) ;
00594 
00595    sl = (byte *) malloc(mab) ;  /* will hold extracted sheared slice */
00596 
00597    /* create all zero output image (now done by mri_new) */
00598 
00599    bim = mri_new(ma,mb,MRI_rgb); rgb = MRI_RGB_PTR(bim);
00600 
00601    /* prepare for projections of different types */
00602 
00603    switch( ar->renmode ){
00604      default:
00605      case CREN_SUM_VOX:
00606         used = (float *) malloc(sizeof(float)*mab) ;       /* how much opacity */
00607         for( pij=0 ; pij < mab ; pij++ ) used[pij] = 0.0 ; /* has been used up */
00608      break ;
00609 
00610      case CREN_MIP_VOX:
00611      break ;
00612 
00613      case CREN_MIP_OPA:
00614      break ;
00615    }
00616 
00617    /* extract sheared slices, then project them */
00618 
00619    switch( ar->intmode ){
00620       default:
00621       case CREN_TWOSTEP: ifunc = extract_byte_tsx ; break ;
00622       case CREN_NN:      ifunc = extract_byte_nn  ; break ;
00623       case CREN_LINEAR:  ifunc = (ar->vox_is_gray) ? extract_byte_lixx
00624                                                    : extract_byte_lix ;
00625                          break ;
00626    }
00627 
00628    for( pk=0 ; pk < nk ; pk++ ){  /* loop over slices */
00629 
00630       ppk = (pk_reverse) ? (nk-1-pk) : pk ;  /* maybe going backwards? */
00631 
00632       ifunc( ar->nx,ar->ny,ar->nz , var , ar->vtm ,
00633              kk+1 , ppk , ba-a*(ppk-hnk) , bb-b*(ppk-hnk) , ma,mb , sl ) ;
00634 
00635       switch( ar->renmode ){
00636 
00637 #undef  MAX_OPACITY
00638 #define MAX_OPACITY 0.95
00639 
00640         default:
00641         case CREN_SUM_VOX:{  /* integrate along the axis */
00642           register float opa ;
00643 
00644           for( p3=pij=0 ; pij < mab ; pij++,p3+=3 ){ /* loop over pixels */
00645 
00646              vv = sl[pij] ; if( vv == 0 ) continue ; /* skip voxel */
00647 
00648              if( used[pij] > MAX_OPACITY ) continue; /* skip voxel */
00649 
00650                   if( vv < 128  ) opa = omap[vv] ;   /* gray  */
00651              else if( vv < vtop ) opa = orgb ;       /* color */
00652              else                 continue ;         /* skip voxel */
00653              if( opa < obot ) continue ;             /* skip voxel */
00654 
00655              opa *= (1.0-used[pij]) ; used[pij] += opa ;
00656 
00657              if( vv < 128 ){                         /* gray */
00658                 vv = (byte)( opa * (vv << 1) ) ;
00659                 rgb[p3  ] += vv ;
00660                 rgb[p3+1] += vv ;
00661                 rgb[p3+2] += vv ;
00662              } else if( vv < vtop ){                 /* color */
00663                 rgb[p3  ] += (byte)( opa * ar->rmap[vv-128] ) ;
00664                 rgb[p3+1] += (byte)( opa * ar->gmap[vv-128] ) ;
00665                 rgb[p3+2] += (byte)( opa * ar->bmap[vv-128] ) ;
00666              }
00667 
00668           }
00669         }
00670         break ;
00671 
00672         /* the MIP are grayscale projections:
00673            only fill in the R value now, and fix the GB values at the end */
00674 
00675         case CREN_MIP_VOX:   /* MIP on the signal intensity */
00676           for( p3=pij=0 ; pij < mab ; pij++,p3+=3 ){
00677              vv = sl[pij] ; if( vv == 0 ) continue ;      /* skip */
00678                   if( vv < 128  ) vv = vv << 1 ;          /* gray  */
00679              else if( vv < vtop ) vv = ar->imap[vv-128] ; /* color */
00680              else                 continue ;              /* skip  */
00681 
00682              if( vv > rgb[p3] ) rgb[p3] = vv ;      /* MIP   */
00683           }
00684         break ;
00685 
00686         case CREN_MIP_OPA:{  /* MIP on the opacity */
00687           float opa ;
00688           for( p3=pij=0 ; pij < mab ; pij++,p3+=3 ){
00689              vv = sl[pij] ; if( vv == 0 ) continue ;      /* skip */
00690                   if( vv < 128  ) opa = omap[vv] ;        /* gray  */
00691              else if( vv < vtop ) opa = orgb ;            /* color */
00692              else                 continue ;              /* skip  */
00693 
00694              vv = (byte)(255.9*opa) ;                     /* scale */
00695              if( vv > rgb[p3] ) rgb[p3] = vv ;      /* MIP   */
00696           }
00697         }
00698         break ;
00699 
00700      } /* end of switch over rendering mode */
00701 
00702    } /* end of loop over slices */
00703 
00704    free(sl) ;
00705 
00706    /*-- finalization --*/
00707 
00708    switch( ar->renmode ){
00709      default:
00710      case CREN_SUM_VOX:
00711         free(used) ;
00712      break ;
00713 
00714      case CREN_MIP_VOX:  /* fill in missing GB values */
00715      case CREN_MIP_OPA:
00716         for( p3=pij=0 ; pij < mab ; pij++,p3+=3 )
00717            rgb[p3+1] = rgb[p3+2] = rgb[p3] ;
00718      break ;
00719    }
00720 
00721    /* warp projection to final display coordinates */
00722 
00723    qim = mri_aff2d_rgb( bim , 1 , aii,aij,aji,ajj ) ;
00724    mri_free(bim) ; return qim ;
00725 }
00726 
00727 /*===========================================================================
00728    Functions to extract a plane of shifted bytes from a 3D volume.
00729      nx, ny, nz = dimensions of vol
00730      vol        = input 3D volume of bytes
00731 
00732      fixdir = fixed direction (1=x, 2=y, 3=z)
00733      fixijk = fixed index
00734      da, db = shift in planar coordinaes (non-fixed directions)
00735      ma, mb = dimensions of im
00736      im     = output 2D image
00737 
00738    Goal is im[a,b] = vol[ P(a-da,b-db,c=fixijk) ] for a=0..ma-1, b=0..mb-1,
00739    where P(a,b,c) is the permutation of (a,b,c) that goes with fixdir:
00740      P(x,y,z) = (y,z,x) for fixdir == 1
00741      P(x,y,z) = (z,x,y) for fixdir == 2
00742      P(x,y,z) = (x,y,z) for fixdir == 3
00743    For values outside the range of vol[], im[] is set to 0.
00744 =============================================================================*/
00745 
00746   /* macros for offsets in vol[] to corners of the interpolation square */
00747 
00748 #undef LL
00749 #undef LR
00750 #undef UL
00751 #undef UR
00752 
00753 #define LL 0                /* voxel offset to lower left  */
00754 #define LR astep            /* voxel offset to lower right */
00755 #define UL bstep            /* voxel offset to upper left  */
00756 #define UR (astep+bstep)    /* voxel offset to upper right */
00757 
00758   /* macro to compute the stepsizes and dimensions in the
00759      3D array for each direction, given the fixed direction */
00760 
00761 #undef  ASSIGN_DIRECTIONS
00762 #define ASSIGN_DIRECTIONS                                       \
00763  do{ switch( fixdir ){                                          \
00764       default:                                                  \
00765       case 1:            /* x-direction: (a,b,c) = (y,z,x) */   \
00766          astep = nx ; bstep = nxy ; cstep = 1  ;                \
00767          na    = ny ; nb    = nz  ; nc    = nx ;                \
00768       break ;                                                   \
00769                                                                 \
00770       case 2:            /* y-direction: (a,b,c) = (z,x,y) */   \
00771          astep = nxy ; bstep = 1  ; cstep = nx ;                \
00772          na    = nz  ; nb    = nx ; nc    = ny ;                \
00773       break ;                                                   \
00774                                                                 \
00775       case 3:            /* z-direction: (a,b,c) = (x,y,z) */   \
00776          astep = 1  ; bstep = nx ; cstep = nxy ;                \
00777          na    = nx ; nb    = ny ; nc    = nz  ;                \
00778       break ;                                                   \
00779     } } while(0)
00780 
00781 /*-----------------------------------------------------------------------*/
00782 
00783 void extract_assign_directions( int nx, int ny, int nz, int fixdir ,
00784                                 int *Astep, int *Bstep, int *Cstep ,
00785                                 int *Na   , int *Nb   , int *Nc     )
00786 {
00787    int astep,bstep,cstep , na,nb,nc , nxy=nx*ny ;
00788 
00789    ASSIGN_DIRECTIONS ;
00790 
00791    *Astep = astep ; *Bstep = bstep ; *Cstep = cstep ;
00792    *Na    = na    ; *Nb    = nb    ; *Nc    = nc    ; return ;
00793 }
00794 
00795 /*-----------------------------------------------------------------------
00796    NN "interpolation"
00797 -------------------------------------------------------------------------*/
00798 
00799 void extract_byte_nn( int nx , int ny , int nz , byte * vol ,
00800                       Tmask * tm ,
00801                       int fixdir , int fixijk , float da , float db ,
00802                       int ma , int mb , byte * im )
00803 {
00804    int adel,bdel , abot,atop , bb,bbot,btop , nxy=nx*ny ;
00805    register int aa , ijkoff , aoff,boff ;
00806    int astep,bstep,cstep , na,nb,nc ;
00807    byte * mask ;
00808 
00809    memset( im , 0 , ma*mb ) ;  /* initialize output to zero */
00810 
00811    if( fixijk < 0 ) return ;
00812 
00813    ASSIGN_DIRECTIONS ;
00814 
00815    if( fixijk >= nc ) return ;
00816 
00817    da += 0.5 ; adel = (int) da ; if( da < 0.0 ) adel-- ;  /* floor(da+0.5) */
00818    db += 0.5 ; bdel = (int) db ; if( db < 0.0 ) bdel-- ;  /* floor(db+0.5) */
00819 
00820    abot = 0       ; if( abot < adel ) abot = adel ;       /* range in im[] */
00821    atop = na+adel ; if( atop > ma   ) atop = ma ;
00822 
00823    bbot = 0       ; if( bbot < bdel ) bbot = bdel ;
00824    btop = nb+bdel ; if( btop > mb   ) btop = mb ;
00825 
00826    ijkoff = fixijk*cstep + (abot-adel)*astep + (bbot-bdel)*bstep ;
00827    boff   = bbot * ma ;
00828 
00829    if( atop <= abot || btop <= bbot ) return ;  /* nothing to do */
00830 
00831    mask = (tm == NULL) ? NULL
00832                        : tm->mask[fixdir%3] + (fixijk*nb - bdel) ;
00833 
00834    if( astep != 1 ){
00835     for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
00836      if( mask == NULL || mask[bb] )
00837       for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
00838        im[aa+boff] = vol[aoff+ijkoff]; /* im(aa,bb)=vol(aa-adel,bb-bdel,fixijk) */
00839                                        /*          =vol[ (aa-adel)*astep +
00840                                                          (bb-bdel)*bstep +
00841                                                          fixijk   *cstep   ]    */
00842    } else { /* 05 Nov 2003 */
00843      for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
00844       if( mask == NULL || mask[bb] )
00845         memcpy( im+(abot+boff) , vol+ijkoff , atop-abot ) ;
00846    }
00847 
00848    return ;
00849 }
00850 
00851 /*-----------------------------------------------------------------------
00852    NN "interpolation" of rgba data - 30 Jan 2003
00853 -------------------------------------------------------------------------*/
00854 
00855 void extract_rgba_nn( int nx , int ny , int nz , rgba * vol ,
00856                       Tmask * tm ,
00857                       int fixdir , int fixijk , float da , float db ,
00858                       int ma , int mb , rgba * im )
00859 {
00860    int adel,bdel , abot,atop , bb,bbot,btop , nxy=nx*ny ;
00861    register int aa , ijkoff , aoff,boff ;
00862    int astep,bstep,cstep , na,nb,nc ;
00863    byte * mask ;
00864 
00865    memset( im , 0 , sizeof(rgba)*ma*mb ) ;  /* initialize output to zero */
00866 
00867    if( fixijk < 0 ) return ;
00868 
00869    ASSIGN_DIRECTIONS ;
00870 
00871    if( fixijk >= nc ) return ;
00872 
00873    da += 0.5 ; adel = (int) da ; if( da < 0.0 ) adel-- ;  /* floor(da+0.5) */
00874    db += 0.5 ; bdel = (int) db ; if( db < 0.0 ) bdel-- ;  /* floor(db+0.5) */
00875 
00876    abot = 0       ; if( abot < adel ) abot = adel ;       /* range in im[] */
00877    atop = na+adel ; if( atop > ma   ) atop = ma ;
00878 
00879    bbot = 0       ; if( bbot < bdel ) bbot = bdel ;
00880    btop = nb+bdel ; if( btop > mb   ) btop = mb ;
00881 
00882    ijkoff = fixijk*cstep + (abot-adel)*astep + (bbot-bdel)*bstep ;
00883    boff   = bbot * ma ;
00884 
00885    mask = (tm == NULL) ? NULL
00886                        : tm->mask[fixdir%3] + (fixijk*nb - bdel) ;
00887 
00888    for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
00889     if( mask == NULL || mask[bb] )
00890      for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
00891       im[aa+boff] = vol[aoff+ijkoff]; /* im(aa,bb)=vol(aa-adel,bb-bdel,fixijk) */
00892                                       /*          =vol[ (aa-adel)*astep +
00893                                                         (bb-bdel)*bstep +
00894                                                         fixijk   *cstep   ]    */
00895 
00896    return ;
00897 }
00898 
00899 /*---------------------------------------------------------------------------
00900     Two-step interpolation
00901 -----------------------------------------------------------------------------*/
00902 
00903 #undef TSBOT
00904 #undef TSTOP
00905 
00906 #if 1
00907 # define TSBOT 0.3
00908 # define TSTOP 0.7
00909 #else
00910 # define TSBOT 0.25
00911 # define TSTOP 0.75
00912 #endif
00913 
00914 /*---------------------------------------------------------------------------*/
00915 #if 0
00916 static void extract_byte_ts( int nx , int ny , int nz , byte * vol ,
00917                              Tmask * tm ,
00918                              int fixdir , int fixijk , float da , float db ,
00919                              int ma , int mb , byte * im )
00920 {
00921    int adel,bdel , abot,atop , bb,bbot,btop , nxy=nx*ny ;
00922    register int aa , ijkoff , aoff,boff ;
00923    int astep,bstep,cstep , na,nb,nc , nts,dts1=0,dts2=0 ;
00924    float fa , fb ;
00925    byte * mask ;
00926 
00927    memset( im , 0 , ma*mb ) ; /* initialize output to zero */
00928 
00929    if( fixijk < 0 ) return ;
00930 
00931    ASSIGN_DIRECTIONS ;
00932 
00933    if( fixijk >= nc ) return ;
00934 
00935    adel = (int) da ; if( da < 0.0 ) adel-- ;  /* floor(da) */
00936    bdel = (int) db ; if( db < 0.0 ) bdel-- ;  /* floor(db) */
00937 
00938    fa = da - adel ;               /* fractional part of da */
00939    fb = db - bdel ;               /* fractional part of db */
00940 
00941    fa = 1.0-fa ; fb = 1.0-fb ;    /* since im[a,b] = vol[a-da,b-db] */
00942 
00943    if( fa < TSBOT ){                      /*- Left 30% -*/
00944       if( fb < TSBOT ){                   /*- Lower 30% -*/
00945         nts = 1 ; dts1 = LL ;               /* [0,0] */
00946       } else if( fb > TSTOP ){            /*- Upper 30% -*/
00947         nts = 1 ; dts1 = UL ;               /* [0,1] */
00948       } else {                            /*- Middle 40% -*/
00949         nts = 2 ; dts1 = LL ; dts2 = UL ;   /* mid of [0,0] and [0,1] */
00950       }
00951    } else if( fa > TSTOP ){               /*- Right 30% -*/
00952       if( fb < TSBOT ){                   /*- Lower 30% -*/
00953         nts = 1 ; dts1 = LR ;               /* [1,0] */
00954       } else if( fb > TSTOP ){            /*- Upper 30% -*/
00955         nts = 1 ; dts1 = UR ;               /* [1,1] */
00956       } else {                            /*- Middle 40% -*/
00957         nts = 2 ; dts1 = LR ; dts2 = UR ;   /* mid of [1,0] and [1,1] */
00958       }
00959    } else {                               /*- Middle 40% -*/
00960       if( fb < TSBOT ){                   /*- Lower 30% -*/
00961         nts = 2 ; dts1 = LL ; dts2 = LR ;   /* mid of [0,0] and [1,0] */
00962       } else if( fb > TSTOP ){            /*- Upper 30% -*/
00963         nts = 2 ; dts1 = UL ; dts2 = UR ;   /* mid of [0,1] and [1,1] */
00964       } else {                            /*- Middle 40% -*/
00965         nts = 4 ;                           /* mid of all 4 points */
00966       }
00967    }
00968 
00969    adel++ ; bdel++ ;
00970 
00971    abot = 0         ; if( abot < adel ) abot = adel ;       /* range in im[] */
00972    atop = na+adel-1 ; if( atop > ma   ) atop = ma ;
00973 
00974    bbot = 0         ; if( bbot < bdel ) bbot = bdel ;
00975    btop = nb+bdel-1 ; if( btop > mb   ) btop = mb ;
00976 
00977    ijkoff = fixijk*cstep + (abot-adel)*astep + (bbot-bdel)*bstep ;
00978    boff   = bbot * ma ;
00979 
00980    mask = (tm == NULL) ? NULL
00981                        : tm->mask[fixdir%3] + (fixijk*nb - bdel) ;
00982 
00983    switch( nts ){
00984 
00985       case 1:
00986          ijkoff += dts1 ;
00987          for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
00988            if( mask == NULL || mask[bb] || mask[bb+1] )
00989              for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
00990                im[aa+boff] = vol[aoff+ijkoff] ;
00991       break ;
00992 
00993       case 2:
00994          ijkoff += dts1 ; dts2 -= dts1 ;
00995          for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
00996            if( mask == NULL || mask[bb] || mask[bb+1] )
00997              for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
00998                im[aa+boff] = (vol[aoff+ijkoff] + vol[aoff+(ijkoff+dts2)]) >> 1;
00999       break ;
01000 
01001       case 4:
01002          for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
01003            if( mask == NULL || mask[bb] || mask[bb+1] )
01004              for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
01005                im[aa+boff] = ( vol[aoff+ijkoff]     +vol[aoff+(ijkoff+LR)]
01006                               +vol[aoff+(ijkoff+UL)]+vol[aoff+(ijkoff+UR)]) >> 2;
01007       break ;
01008    }
01009 
01010    return ;
01011 }
01012 #endif
01013 
01014 /*---------------------------------------------------------------------------*/
01015 
01016 void extract_byte_tsx( int nx , int ny , int nz , byte * vol ,
01017                        Tmask * tm ,
01018                        int fixdir , int fixijk , float da , float db ,
01019                        int ma , int mb , byte * im )
01020 {
01021    int adel,bdel , abot,atop , bb,bbot,btop , nxy=nx*ny ;
01022    register int aa , ijkoff , aoff,boff ;
01023    int astep,bstep,cstep , na,nb,nc , nts,dts1=0,dts2=0 , nn ;
01024    float fa , fb ;
01025    byte *mask , v1,v2,v3,v4 ;
01026 
01027    memset( im , 0 , ma*mb ) ; /* initialize output to zero */
01028 
01029    if( fixijk < 0 ) return ;
01030 
01031    ASSIGN_DIRECTIONS ;
01032 
01033    if( fixijk >= nc ) return ;
01034 
01035    adel = (int) da ; if( da < 0.0 ) adel-- ;  /* floor(da) */
01036    bdel = (int) db ; if( db < 0.0 ) bdel-- ;  /* floor(db) */
01037 
01038    fa = da - adel ;               /* fractional part of da */
01039    fb = db - bdel ;               /* fractional part of db */
01040 
01041    fa = 1.0-fa ; fb = 1.0-fb ;    /* since im[a,b] = vol[a-da,b-db] */
01042 
01043    if( fa < TSBOT ){                      /*- Left 30% -*/
01044       if( fb < TSBOT ){                   /*- Lower 30% -*/
01045         nts = 1 ; dts1 = LL ;               /* [0,0] */
01046       } else if( fb > TSTOP ){            /*- Upper 30% -*/
01047         nts = 1 ; dts1 = UL ;               /* [0,1] */
01048       } else {                            /*- Middle 40% -*/
01049         nts = 2 ; dts1 = LL ; dts2 = UL ;   /* mid of [0,0] and [0,1] */
01050       }
01051    } else if( fa > TSTOP ){               /*- Right 30% -*/
01052       if( fb < TSBOT ){                   /*- Lower 30% -*/
01053         nts = 1 ; dts1 = LR ;               /* [1,0] */
01054       } else if( fb > TSTOP ){            /*- Upper 30% -*/
01055         nts = 1 ; dts1 = UR ;               /* [1,1] */
01056       } else {                            /*- Middle 40% -*/
01057         nts = 2 ; dts1 = LR ; dts2 = UR ;   /* mid of [1,0] and [1,1] */
01058       }
01059    } else {                               /*- Middle 40% -*/
01060       if( fb < TSBOT ){                   /*- Lower 30% -*/
01061         nts = 2 ; dts1 = LL ; dts2 = LR ;   /* mid of [0,0] and [1,0] */
01062       } else if( fb > TSTOP ){            /*- Upper 30% -*/
01063         nts = 2 ; dts1 = UL ; dts2 = UR ;   /* mid of [0,1] and [1,1] */
01064       } else {                            /*- Middle 40% -*/
01065         nts = 4 ;                           /* mid of all 4 points */
01066       }
01067    }
01068 
01069    nn = (fa < 0.5) ? ( (fb < 0.5) ? LL : UL )   /* NN index point */
01070                    : ( (fb < 0.5) ? LR : UR ) ;
01071 
01072    adel++ ; bdel++ ;
01073 
01074    abot = 0         ; if( abot < adel ) abot = adel ;       /* range in im[] */
01075    atop = na+adel-1 ; if( atop > ma   ) atop = ma ;
01076 
01077    bbot = 0         ; if( bbot < bdel ) bbot = bdel ;
01078    btop = nb+bdel-1 ; if( btop > mb   ) btop = mb ;
01079 
01080    ijkoff = fixijk*cstep + (abot-adel)*astep + (bbot-bdel)*bstep ;
01081    boff   = bbot * ma ;
01082 
01083    mask = (tm == NULL) ? NULL
01084                        : tm->mask[fixdir%3] + (fixijk*nb - bdel) ;
01085 
01086    /*-- the following is intended to implement
01087 
01088          im(aa,bb) = vol(aa-adel,bb-bdel,fixijk)
01089                    = vol[ (aa-adel)*astep +
01090                           (bb-bdel)*bstep +
01091                           fixijk   *cstep   ]    --*/
01092 
01093 #define BECLEVER
01094 
01095    switch( nts ){
01096 
01097       case 1:
01098          ijkoff += dts1 ;
01099          for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
01100            if( mask == NULL || mask[bb] || mask[bb+1] )
01101              for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
01102                im[aa+boff] = vol[aoff+ijkoff] ;
01103       break ;
01104 
01105       case 2:
01106          ijkoff += dts1 ; dts2 -= dts1 ; nn -= dts1 ;
01107          for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
01108            if( mask == NULL || mask[bb] || mask[bb+1] )
01109              for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep ){
01110                v1 = vol[aoff+ijkoff] ;
01111                v2 = vol[aoff+(ijkoff+dts2)] ;
01112 #ifdef BECLEVER
01113                if( (v1|v2) & 128 != 0 )
01114 #else
01115                if( v1 < 128 && v2 < 128 )
01116 #endif
01117                   im[aa+boff] = (v1+v2) >> 1 ;           /* grayscale */
01118                else
01119                   im[aa+boff] = vol[aoff+(ijkoff+nn)] ;  /* color code */
01120              }
01121       break ;
01122 
01123       case 4:
01124          for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
01125            if( mask == NULL || mask[bb] || mask[bb+1] )
01126              for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep ){
01127                v1 = vol[aoff+ijkoff] ;
01128                v2 = vol[aoff+(ijkoff+LR)] ;
01129                v3 = vol[aoff+(ijkoff+UL)] ;
01130                v4 = vol[aoff+(ijkoff+UR)] ;
01131 #ifdef BECLEVER
01132                if( (v1|v2|v3|v4) & 128 != 0 )
01133 #else
01134                if( v1 < 128 && v2 < 128 && v3 < 128 && v4 < 128 )
01135 #endif
01136                   im[aa+boff] = (v1+v2+v3+v4) >> 2 ;     /* grayscale */
01137                else
01138                   im[aa+boff] = vol[aoff+(ijkoff+nn)] ;  /* color code */
01139             }
01140       break ;
01141    }
01142 
01143    return ;
01144 }
01145 
01146 /*---------------------------------------------------------------------------*/
01147 
01148 void extract_byte_lix( int nx , int ny , int nz , byte * vol ,
01149                        Tmask * tm ,
01150                        int fixdir , int fixijk , float da , float db ,
01151                        int ma , int mb , byte *im )
01152 {
01153    int adel,bdel , abot,atop , bb,bbot,btop , nxy=nx*ny ;
01154    register int aa , ijkoff , aoff,boff ;
01155    int astep,bstep,cstep , na,nb,nc , nn ;
01156    float fa , fb ;
01157    float f_a_b , f_ap_b , f_a_bp , f_ap_bp ;
01158    byte  b_a_b , b_ap_b , b_a_bp , b_ap_bp ;
01159    byte *mask , v1,v2,v3,v4 ;
01160 
01161    memset( im , 0 , ma*mb ) ; /* initialize output to zero */
01162 
01163    if( fixijk < 0 ) return ;
01164 
01165    ASSIGN_DIRECTIONS ;
01166 
01167    if( fixijk >= nc ) return ;
01168 
01169    adel = (int) da ; if( da < 0.0 ) adel-- ;  /* floor(da) */
01170    bdel = (int) db ; if( db < 0.0 ) bdel-- ;  /* floor(db) */
01171 
01172    fa = da - adel ;               /* fractional part of da */
01173    fb = db - bdel ;               /* fractional part of db */
01174 
01175    f_a_b   = fa      * fb      ;  /* bilinear interpolation */
01176    f_ap_b  = (1.0-fa)* fb      ;  /* coefficients */
01177    f_a_bp  = fa      *(1.0-fb) ;
01178    f_ap_bp = (1.0-fa)*(1.0-fb) ;
01179 
01180    bb = (int)(256*f_a_b  + 0.499); if( bb == 256 ) bb--; b_a_b  = (byte) bb;
01181    bb = (int)(256*f_ap_b + 0.499); if( bb == 256 ) bb--; b_ap_b = (byte) bb;
01182    bb = (int)(256*f_a_bp + 0.499); if( bb == 256 ) bb--; b_a_bp = (byte) bb;
01183    bb = (int)(256*f_ap_bp+ 0.499); if( bb == 256 ) bb--; b_ap_bp= (byte) bb;
01184 
01185    nn = (fa > 0.5) ? ( (fb > 0.5) ? LL : UL )   /* NN index point */
01186                    : ( (fb > 0.5) ? LR : UR ) ;
01187 
01188    adel++ ; bdel++ ;
01189 
01190    abot = 0         ; if( abot < adel ) abot = adel ;       /* range in im[] */
01191    atop = na+adel-1 ; if( atop > ma   ) atop = ma ;
01192 
01193    bbot = 0         ; if( bbot < bdel ) bbot = bdel ;
01194    btop = nb+bdel-1 ; if( btop > mb   ) btop = mb ;
01195 
01196    ijkoff = fixijk*cstep + (abot-adel)*astep + (bbot-bdel)*bstep ;
01197    boff   = bbot * ma ;
01198 
01199    mask = (tm == NULL) ? NULL
01200                        : tm->mask[fixdir%3] + (fixijk*nb - bdel) ;
01201 
01202    /*-- the following is intended to implement
01203 
01204          im(aa,bb) = vol(aa-adel,bb-bdel,fixijk)
01205                    = vol[ (aa-adel)*astep +
01206                           (bb-bdel)*bstep +
01207                           fixijk   *cstep   ]    --*/
01208 
01209    for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
01210      if( mask == NULL || mask[bb] || mask[bb+1] )
01211        for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep ){
01212          v1 = vol[aoff+ijkoff]      ;
01213          v2 = vol[aoff+(ijkoff+LR)] ;
01214          v3 = vol[aoff+(ijkoff+UL)] ;
01215          v4 = vol[aoff+(ijkoff+UR)] ;
01216 #ifdef BECLEVER
01217          if( (v1|v2|v3|v4) & 128 != 0 )
01218 #else
01219          if( v1 < 128 && v2 < 128 && v3 < 128 && v4 < 128 )  /* gray */
01220 #endif
01221 
01222 #if 0
01223            im[aa+boff] = (byte)(  f_a_b  * v1 + f_ap_b  * v2
01224                                 + f_a_bp * v3 + f_ap_bp * v4 ) ;
01225 #else
01226            im[aa+boff] = (byte)((  b_a_b  * v1 + b_ap_b  * v2
01227                                  + b_a_bp * v3 + b_ap_bp * v4 ) >> 8) ;
01228 #endif
01229          else
01230            im[aa+boff] = vol[aoff+(ijkoff+nn)] ;  /* color code */
01231        }
01232 
01233    return ;
01234 }
01235 
01236 /*---------------------------------------------------------------------------*/
01237 
01238 void extract_byte_lixx( int nx , int ny , int nz , byte * vol ,
01239                         Tmask * tm ,
01240                         int fixdir , int fixijk , float da , float db ,
01241                         int ma , int mb , byte * im )
01242 {
01243    int adel,bdel , abot,atop , bb,bbot,btop , nxy=nx*ny ;
01244    register int aa , ijkoff , aoff,boff ;
01245    int astep,bstep,cstep , na,nb,nc , nn ;
01246    float fa , fb ;
01247    float f_a_b , f_ap_b , f_a_bp , f_ap_bp ;
01248    byte  b_a_b , b_ap_b , b_a_bp , b_ap_bp ;
01249    byte *mask , v1,v2,v3,v4 ;
01250 
01251    memset( im , 0 , ma*mb ) ; /* initialize output to zero */
01252 
01253    if( fixijk < 0 ) return ;
01254 
01255    ASSIGN_DIRECTIONS ;
01256 
01257    if( fixijk >= nc ) return ;
01258 
01259    adel = (int) da ; if( da < 0.0 ) adel-- ;  /* floor(da) */
01260    bdel = (int) db ; if( db < 0.0 ) bdel-- ;  /* floor(db) */
01261 
01262    fa = da - adel ;               /* fractional part of da */
01263    fb = db - bdel ;               /* fractional part of db */
01264 
01265    f_a_b   = fa      * fb      ;  /* bilinear interpolation */
01266    f_ap_b  = (1.0-fa)* fb      ;  /* coefficients */
01267    f_a_bp  = fa      *(1.0-fb) ;
01268    f_ap_bp = (1.0-fa)*(1.0-fb) ;
01269 
01270    bb = (int)(256*f_a_b  + 0.499); if( bb == 256 ) bb--; b_a_b  = (byte) bb;
01271    bb = (int)(256*f_ap_b + 0.499); if( bb == 256 ) bb--; b_ap_b = (byte) bb;
01272    bb = (int)(256*f_a_bp + 0.499); if( bb == 256 ) bb--; b_a_bp = (byte) bb;
01273    bb = (int)(256*f_ap_bp+ 0.499); if( bb == 256 ) bb--; b_ap_bp= (byte) bb;
01274 
01275    nn = (fa > 0.5) ? ( (fb > 0.5) ? LL : UL )   /* NN index point */
01276                    : ( (fb > 0.5) ? LR : UR ) ;
01277 
01278    adel++ ; bdel++ ;
01279 
01280    abot = 0         ; if( abot < adel ) abot = adel ;       /* range in im[] */
01281    atop = na+adel-1 ; if( atop > ma   ) atop = ma ;
01282 
01283    bbot = 0         ; if( bbot < bdel ) bbot = bdel ;
01284    btop = nb+bdel-1 ; if( btop > mb   ) btop = mb ;
01285 
01286    if( atop <= abot || btop <= bbot ) return ;  /* nothing to do */
01287 
01288    ijkoff = fixijk*cstep + (abot-adel)*astep + (bbot-bdel)*bstep ;
01289    boff   = bbot * ma ;
01290 
01291    mask = (tm == NULL) ? NULL
01292                        : tm->mask[fixdir%3] + (fixijk*nb - bdel) ;
01293 
01294    /*-- the following is intended to implement
01295 
01296          im(aa,bb) = vol(aa-adel,bb-bdel,fixijk)
01297                    = vol[ (aa-adel)*astep +
01298                           (bb-bdel)*bstep +
01299                           fixijk   *cstep   ]    --*/
01300 
01301    if( astep != 1 ){
01302     for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
01303       if( mask == NULL || mask[bb] || mask[bb+1] )
01304         for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep ){
01305           v1 = vol[aoff+ijkoff]      ;
01306           v2 = vol[aoff+(ijkoff+LR)] ;
01307           v3 = vol[aoff+(ijkoff+UL)] ;
01308           v4 = vol[aoff+(ijkoff+UR)] ;
01309           im[aa+boff] = (byte)((  b_a_b  * v1 + b_ap_b  * v2
01310                                 + b_a_bp * v3 + b_ap_bp * v4 ) >> 8) ;
01311         }
01312    } else {
01313     ijkoff -= abot ;  /* aoff==aa-abot when astep==1 */
01314     for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
01315       if( mask == NULL || mask[bb] || mask[bb+1] )
01316         for( aa=abot ; aa < atop ; aa++ ){
01317           v1 = vol[aa+ijkoff]      ;
01318           v2 = vol[aa+(ijkoff+LR)] ;
01319           v3 = vol[aa+(ijkoff+UL)] ;
01320           v4 = vol[aa+(ijkoff+UR)] ;
01321           im[aa+boff] = (byte)((  b_a_b  * v1 + b_ap_b  * v2
01322                                 + b_a_bp * v3 + b_ap_bp * v4 ) >> 8) ;
01323         }
01324    }
01325 
01326    return ;
01327 }
01328 
01329 /*------------------------------------------------------------------------------*/
01330 
01331 THD_mat33 rotmatrix( int ax1,float th1 ,
01332                      int ax2,float th2 , int ax3,float th3  )
01333 {
01334    THD_mat33 q , p ;
01335 
01336    LOAD_ROT_MAT( q , th1 , ax1 ) ;
01337    LOAD_ROT_MAT( p , th2 , ax2 ) ; q = MAT_MUL( p , q ) ;
01338    LOAD_ROT_MAT( p , th3 , ax3 ) ; q = MAT_MUL( p , q ) ;
01339 
01340    return q ;
01341 }
 
Powered by Plone

This site conforms to the following standards: