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 }

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 }