do{ switch( fixdir ){                                          \
      default:                                                  \
      case 1:                \
         astep = nx ; bstep = nxy ; cstep = 1  ;                \
         na    = ny ; nb    = nz  ; nc    = nx ;                \
      break ;                                                   \
                                                                \
      case 2:                \
         astep = nxy ; bstep = 1  ; cstep = nx ;                \
         na    = nz  ; nb    = nx ; nc    = ny ;                \
      break ;                                                   \
                                                                \
      case 3:                \
         astep = 1  ; bstep = nx ; cstep = nxy ;                \
         na    = nx ; nb    = ny ; nc    = nz  ;                \
      break ;                                                   \
    } } while(0)

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 }

00256 {
00257    CREN_stuff * ar = (CREN_stuff *) ah ;
00258 
00259    if( !ISVALID_CREN(ar) ) return -1 ;
00260 
00261    return (ar->vox == NULL) ;
00262 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }

00133 {
00134    CREN_stuff * ar = (CREN_stuff *) ah ;
00135 
00136    if( !ISVALID_CREN(ar) ) return ;
00137 
00138    ar->skewmat = sm ; return ;
00139 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }