00001 #include "aren.h"
00002 
00003 /*****************************************************************************
00004   This software is copyrighted and owned by the Medical College of Wisconsin.
00005   See the file README.Copyright for details.
00006 ******************************************************************************/
00007 
00008 /*============================================================================
00009   AREN = AFNI Renderer, a set of routines for volume rendering 3D bricks.
00010   Built on top of VolPack, by Philippe Lacroute.
00011 
00012   (1) Use new_AREN_renderer() to create a renderer.
00013   (2) Load grayscale or color bricks into it with functions
00014       AREN_set_graybrick() or AREN_set_rgbbricks().
00015   (3) Load the opacity brick with AREN_set_opabrick().
00016   (4) Set the viewing angles with AREN_set_viewpoint().
00017       Set the brick dimensions with AREN_set_size().
00018   (5) Create an image with AREN_render().
00019       Loop back through (2-4) as needed to create new images.
00020   (6) When finished, use destroy_AREN_renderer() to clean up.
00021 ==============================================================================*/
00022 
00023 /*--------------------------------------------------------------------------
00024   Initialize the short -> float color table, and
00025              the byte  -> float gray/opacity table for VolPack.
00026   This routine only does something the first time it is called -- it
00027   creates space for the color tables and fills them.  They will never
00028   be deallocated.
00029 
00030   N.B.: The documentation for VolPack is incorrect -- the color values
00031         should be in the range 0..255, NOT 0..1.  The opacity values
00032         should still be in the range 0..1.
00033 ----------------------------------------------------------------------------*/
00034 
00035 static float * AREN_colortable = NULL ;
00036 static float * AREN_graytable  = NULL ;
00037 static float * AREN_opatable   = NULL ;
00038 
00039 void init_AREN_colortable( void )
00040 {
00041    static int been_here = 0 ;
00042    int ii , rr,gg,bb , ss ;
00043 
00044    if( been_here ) return ;  /* done that */
00045 
00046    AREN_colortable = (float *) malloc( sizeof(float) * TOT_COLORS * 3 ) ;
00047    AREN_graytable  = (float *) malloc( sizeof(float) * MAX_GRAYS ) ;
00048    AREN_opatable   = (float *) malloc( sizeof(float) * MAX_GRAYS ) ;
00049 
00050    /*-- load linear ramp for grayscale and opacity --*/
00051 
00052    for( ii=0 ; ii < MAX_GRAYS ; ii++ ){
00053       AREN_graytable[ii] = ii ;
00054       AREN_opatable[ii]  = ii / 255.0 ;
00055    }
00056 
00057    /*-- load a 32 x 32 x 32 color cube --*/
00058 
00059    for( rr=0 ; rr < MAX_CDIM ; rr++ ){
00060       for( gg=0 ; gg < MAX_CDIM ; gg++ ){
00061          for( bb=0 ; bb < MAX_CDIM ; bb++ ){
00062 
00063             ss = FIVE_to_short(rr,gg,bb) ;   /* color index */
00064 
00065             AREN_colortable[3*ss  ] = (rr * 255.0) / 31.0 ;
00066             AREN_colortable[3*ss+1] = (gg * 255.0) / 31.0 ;
00067             AREN_colortable[3*ss+2] = (bb * 255.0) / 31.0 ;
00068          }
00069       }
00070    }
00071 
00072    /*-- at the end, add the pure grays (at a higher resolution) --*/
00073 
00074    ss = 3 * MAX_COLORS ;
00075    for( ii=0 ; ii < MAX_GRAYS ; ii++ ){
00076       AREN_colortable[ss++] = ii ;
00077       AREN_colortable[ss++] = ii ;
00078       AREN_colortable[ss++] = ii ;
00079    }
00080 
00081    been_here = 1 ;
00082    return ;
00083 }
00084 
00085 /*--------------------------------------------------------------------------
00086   Create a new AFNI renderer.
00087   The return pointer is passed to all other AREN routines.
00088 ----------------------------------------------------------------------------*/
00089 
00090 #define DEFAULT_THETA 130.0
00091 #define DEFAULT_PHI   285.0
00092 
00093 void * new_AREN_renderer( void )
00094 {
00095    AREN_stuff * ar ;
00096 
00097    ar = (AREN_stuff *) malloc( sizeof(AREN_stuff) ) ;
00098    ar->type = AREN_TYPE ;
00099 
00100    init_AREN_colortable() ;  /* only does something 1st time in */
00101 
00102    /*-- initialize VolPack --*/
00103 
00104    ar->vpc = vpCreateContext() ;
00105 
00106    vpSeti( ar->vpc , VP_CONCAT_MODE , VP_CONCAT_LEFT ) ;
00107 
00108    vpCurrentMatrix( ar->vpc , VP_MODEL ) ;
00109    vpIdentityMatrix( ar->vpc ) ;
00110 
00111    vpCurrentMatrix( ar->vpc , VP_VIEW ) ;
00112    vpIdentityMatrix( ar->vpc ) ;
00113    vpRotate( ar->vpc , VP_Y_AXIS , DEFAULT_THETA ) ;
00114    vpRotate( ar->vpc , VP_X_AXIS , DEFAULT_PHI   ) ;
00115 
00116 #undef USE_CUEING
00117 #ifdef USE_CUEING
00118    vpSetDepthCueing( ar->vpc , 1.0 , 0.5 ) ;
00119    vpEnable( ar->vpc , VP_DEPTH_CUE , 1 ) ;
00120 #endif
00121 
00122    vpCurrentMatrix( ar->vpc , VP_PROJECT ) ;
00123    vpIdentityMatrix( ar->vpc ) ;
00124    vpWindow( ar->vpc , VP_PARALLEL , -0.5,0.5 , -0.5,0.5 , -0.5,0.5 ) ;
00125 
00126    /*-- initialize the rest of the data --*/
00127 
00128    ar->nx = ar->ny = ar->nz = ar->verbose = ar->newdata = ar->newvox = 0 ;
00129    ar->sx = ar->sy = ar->sz = 1.0 ;
00130 
00131    ar->theta = DEFAULT_THETA ;
00132    ar->phi   = DEFAULT_PHI ;
00133    ar->grim  = ar->rim = ar->gim = ar->bim = ar->opim = NULL ;
00134    ar->vox   = NULL ;
00135 
00136    return (void *) ar ;
00137 }
00138 
00139 /*-----------------------------------------------------------------------------
00140    Get rid of an AFNI renderer.
00141 -------------------------------------------------------------------------------*/
00142 
00143 void destroy_AREN_renderer( void * ah )
00144 {
00145    AREN_stuff * ar = (AREN_stuff *) ah ;
00146 
00147    if( !ISVALID_AREN(ar) ) return ;
00148 
00149    if( ar->vox != NULL ) free(ar->vox) ;
00150    vpDestroyContext( ar->vpc ) ;
00151    free(ar) ;
00152    return ;
00153 }
00154 
00155 /*-----------------------------------------------------------------------------
00156    For debugging purposes
00157 -------------------------------------------------------------------------------*/
00158 
00159 void AREN_be_verbose( void * ah )
00160 {
00161    AREN_stuff * ar = (AREN_stuff *) ah ;
00162    if( !ISVALID_AREN(ar) ) return ;
00163    ar->verbose = 1 ; return ;
00164 }
00165 
00166 void AREN_be_quiet( void * ah )
00167 {
00168    AREN_stuff * ar = (AREN_stuff *) ah ;
00169    if( !ISVALID_AREN(ar) ) return ;
00170    ar->verbose = 0 ; return ;
00171 }
00172 
00173 /*-----------------------------------------------------------------------------
00174    Set the grayscale brick in a renderer.
00175    Returns -1 if an error occurs, otherwise returns 0.
00176    Note that this brick is NOT free-d by AREN at any point -- that is
00177    the user's responsibility.
00178 -------------------------------------------------------------------------------*/
00179 
00180 int AREN_set_graybrick( void * ah , MRI_IMAGE * grim )
00181 {
00182    AREN_stuff * ar = (AREN_stuff *) ah ;
00183    int newvox=0 , nvox,ii ;
00184    byte    * gar ;
00185    rgbvox  * gvox ;
00186    byte      gb ;
00187 
00188    /*-- sanity checks --*/
00189 
00190    if( !ISVALID_AREN(ar) || grim == NULL || grim->kind != MRI_byte ) return -1 ;
00191 
00192    if( grim->nx < 3 || grim->ny < 3 || grim->nz < 3 ){
00193       fprintf(stderr,"**AREN: illegal dimensions for a gray brick\n") ;
00194       return -1 ;
00195    }
00196 
00197    /*-- if had old color bricks, toss them (or their pointers) --*/
00198 
00199    if( ar->rim != NULL ){
00200       ar->rim = ar->gim = ar->bim = NULL ;
00201       if( ar->verbose )
00202          fprintf(stderr,"--AREN: switching from color bricks to gray brick\n") ;
00203    } else {
00204       if( ar->verbose )
00205          fprintf(stderr,"--AREN: input a new gray brick\n") ;
00206    }
00207 
00208    /*-- if have new dimensions, toss old stuff that doesn't match --*/
00209 
00210    if( ar->nx > 0 &&
00211        ( ar->nx != grim->nx || ar->ny != grim->ny || ar->nz != grim->nz ) ){
00212 
00213       ar->opim = NULL ;
00214       if( ar->vox != NULL ){ free(ar->vox) ; ar->vox = NULL ; }
00215 
00216       if( ar->verbose )
00217          fprintf(stderr,"--AREN: new gray brick changes volume dimensions\n"
00218                         "        nx:%d->%d  ny:%d->%d  nz:%d->%d\n",
00219                         ar->nx,grim->nx , ar->ny,grim->ny , ar->nz,grim->nz ) ;
00220    }
00221 
00222    /*-- set dimensions --*/
00223 
00224    ar->grim = grim ;
00225    ar->nx   = grim->nx ;
00226    ar->ny   = grim->ny ;
00227    ar->nz   = grim->nz ; nvox = ar->nx * ar->ny * ar->nz ;
00228 
00229    /*-- if need be, allocate a voxel array to hold the data --*/
00230 
00231    if( ar->vox == NULL ){
00232       ar->newvox = newvox = 1 ;
00233       ar->vox = (rgbvox *) malloc( sizeof(rgbvox) * nvox ) ;
00234       if( ar->vox == NULL ){
00235          fprintf(stderr,"**AREN: can't malloc workspace with new gray brick\n") ;
00236          return -1 ;
00237       } else if( ar->verbose ){
00238          fprintf(stderr,"--AREN: allocated new voxel array\n") ;
00239       }
00240    }
00241 
00242    /*-- copy grayscale data into voxel array --*/
00243 
00244    gvox = ar->vox ;
00245    gar  = MRI_BYTE_PTR(grim) ;
00246    for( ii=0 ; ii < nvox ; ii++ ) gvox[ii].rgb = (unsigned short) gar[ii] ;
00247 
00248    ar->newdata = 1 ;
00249    return 0 ;
00250 }
00251 
00252 /*-----------------------------------------------------------------------------
00253    Set the opacity brick in a renderer.
00254    Returns -1 if an error occurs, otherwise returns 0.
00255    Note that this brick is NOT free-d by AREN at any point -- that is
00256    the user's responsibility.
00257 -------------------------------------------------------------------------------*/
00258 
00259 int AREN_set_opabrick( void * ah , MRI_IMAGE * opim )
00260 {
00261    AREN_stuff * ar = (AREN_stuff *) ah ;
00262    int nvox,ii , newvox ;
00263    byte    * gar ;
00264    rgbvox  * rvox ;
00265 
00266    /*-- sanity checks --*/
00267 
00268    if( !ISVALID_AREN(ar) || opim == NULL || opim->kind != MRI_byte ) return -1 ;
00269 
00270    if( opim->nx < 3 || opim->ny < 3 || opim->nz < 3 ){
00271       fprintf(stderr,"**AREN: illegal dimensions for an opacity brick\n") ;
00272       return -1 ;
00273    }
00274 
00275    /*-- if have new dimensions, toss old stuff that doesn't match --*/
00276 
00277    if( ar->nx > 0 &&
00278        ( ar->nx != opim->nx || ar->ny != opim->ny || ar->nz != opim->nz ) ){
00279 
00280       ar->grim = ar->rim = ar->gim = ar->bim = NULL ;
00281       if( ar->vox != NULL ){ free(ar->vox) ; ar->vox = NULL ; }
00282 
00283       if( ar->verbose )
00284          fprintf(stderr,"--AREN: new opacity brick changes volume dimensions\n"
00285                         "        nx:%d->%d  ny:%d->%d  nz:%d->%d\n",
00286                         ar->nx,opim->nx , ar->ny,opim->ny , ar->nz,opim->nz ) ;
00287    } else {
00288       if( ar->verbose )
00289          fprintf(stderr,"--AREN: new opacity brick\n") ;
00290    }
00291 
00292    /*-- set dimensions --*/
00293 
00294    ar->opim = opim ;
00295    ar->nx   = opim->nx ;
00296    ar->ny   = opim->ny ;
00297    ar->nz   = opim->nz ; nvox = ar->nx * ar->ny * ar->nz ;
00298 
00299    /*-- if need be, allocate a voxel array to hold the data --*/
00300 
00301    if( ar->vox == NULL ){
00302       ar->newvox = newvox = 1 ;
00303       ar->vox = (rgbvox *) malloc( sizeof(rgbvox) * nvox ) ;
00304       if( ar->vox == NULL ){
00305          fprintf(stderr,"**AREN: can't malloc workspace with new opacity brick\n") ;
00306          return -1 ;
00307       } else if( ar->verbose ){
00308          fprintf(stderr,"--AREN: allocated new voxel array\n") ;
00309       }
00310    }
00311 
00312    /*-- load the opacity into voxel array --*/
00313 
00314    gar  = MRI_BYTE_PTR(ar->opim) ;
00315    rvox = ar->vox ;
00316    for( ii=0 ; ii < nvox ; ii++ ) rvox[ii].alpha = (unsigned short) gar[ii] ;
00317 
00318    ar->newdata = 1 ;
00319    return 0 ;
00320 }
00321 
00322 /*-----------------------------------------------------------------------------
00323    Set the color bricks in a renderer.
00324    Returns -1 if an error occurs, otherwise returns 0.
00325    Note that these bricks are NOT free-d by AREN at any point -- that is
00326    the user's responsibility.
00327 -------------------------------------------------------------------------------*/
00328 
00329 int AREN_set_rgbbricks( void * ah , MRI_IMAGE *rim, MRI_IMAGE *gim, MRI_IMAGE *bim )
00330 {
00331    AREN_stuff * ar = (AREN_stuff *) ah ;
00332    int newvox=0 , nvox,ii ;
00333    byte    * rar , * gar , * bar ;
00334    rgbvox  * rvox ;
00335    byte      rb,gb,bb ;
00336 
00337    /*-- sanity checks --*/
00338 
00339    if( rim == NULL || rim->kind != MRI_byte ||
00340        gim == NULL || gim->kind != MRI_byte ||
00341        bim == NULL || bim->kind != MRI_byte || !ISVALID_AREN(ar) ) return -1 ;
00342 
00343    if( rim->nx < 3 || rim->ny < 3 || rim->nz < 3 ){
00344       fprintf(stderr,"**AREN: illegal dimensions for a color brick\n") ; return -1 ;
00345    }
00346 
00347    if( gim->nx != rim->nx || gim->ny != rim->ny || gim->nz != rim->nz ){
00348       fprintf(stderr,"**AREN: green & red rgb bricks don't match\n") ; return -1 ;
00349    }
00350 
00351    if( bim->nx != rim->nx || bim->ny != rim->ny || bim->nz != rim->nz ){
00352       fprintf(stderr,"**AREN: blue & red rgb bricks don't match\n") ; return -1 ;
00353    }
00354 
00355    /*-- if had an old gray brick, toss it (or at least its pointer) --*/
00356 
00357    if( ar->grim != NULL ){
00358       ar->grim = NULL ;
00359       if( ar->verbose )
00360          fprintf(stderr,"--AREN: switching from gray brick to color bricks\n") ;
00361    } else {
00362       if( ar->verbose )
00363          fprintf(stderr,"--AREN: input new color bricks\n") ;
00364    }
00365 
00366    /*-- if have new dimensions, toss old stuff that doesn't match --*/
00367 
00368    if( ar->nx > 0 &&
00369        ( ar->nx != rim->nx || ar->ny != rim->ny || ar->nz != rim->nz ) ){
00370 
00371       ar->opim = NULL ;
00372       if( ar->vox != NULL ){ free(ar->vox) ; ar->vox = NULL ; }
00373 
00374       if( ar->verbose )
00375          fprintf(stderr,"--AREN: new color bricks change volume dimensions\n"
00376                         "        nx:%d->%d  ny:%d->%d  nz:%d->%d\n",
00377                         ar->nx,rim->nx , ar->ny,rim->ny , ar->nz,rim->nz ) ;
00378    }
00379 
00380    /*-- set dimensions --*/
00381 
00382    ar->rim = rim ; ar->gim = gim ; ar->bim = bim ;
00383    ar->nx  = rim->nx ;
00384    ar->ny  = rim->ny ;
00385    ar->nz  = rim->nz ; nvox = ar->nx * ar->ny * ar->nz ;
00386 
00387    /*-- if need be, allocate a voxel array to hold the data --*/
00388 
00389    if( ar->vox == NULL ){
00390       ar->newvox = newvox = 1 ;
00391       ar->vox = (rgbvox *) malloc( sizeof(rgbvox) * nvox ) ;
00392       if( ar->vox == NULL ){
00393          fprintf(stderr,"**AREN: can't malloc workspace with new color bricks\n") ;
00394          return -1 ;
00395       } else if( ar->verbose ){
00396          fprintf(stderr,"--AREN: allocated new voxel array\n") ;
00397       }
00398    }
00399 
00400    /*-- copy color data into voxel array --*/
00401 
00402    rvox = ar->vox ;
00403    rar  = MRI_BYTE_PTR(rim) ;
00404    gar  = MRI_BYTE_PTR(gim) ;
00405    bar  = MRI_BYTE_PTR(bim) ;
00406    for( ii=0 ; ii < nvox ; ii++ ){
00407       rb = EIGHT_TO_FIVE(rar[ii]) ;
00408       gb = EIGHT_TO_FIVE(gar[ii]) ;
00409       bb = EIGHT_TO_FIVE(bar[ii]) ;
00410 
00411       if( rb == gb && rb == bb ){
00412          rvox[ii].rgb = MAX_COLORS + rar[ii] ;         /* index into grayscale */
00413       } else {
00414          rvox[ii].rgb = FIVE_to_short( rb , gb, bb ) ; /* index into color cube */
00415       }
00416    }
00417 
00418    ar->newdata = 1 ;
00419    return 0 ;
00420 }
00421 
00422 /*-------------------------------------------------------------------------------
00423    Set the viewpoint of the user (in polar angles)
00424 ---------------------------------------------------------------------------------*/
00425 
00426 void AREN_set_viewpoint( void * ah , float theta , float phi )
00427 {
00428    AREN_stuff * ar = (AREN_stuff *) ah ;
00429 
00430    if( !ISVALID_AREN(ar) ) return ;
00431 
00432    ar->theta = theta ; ar->phi = phi ;
00433 
00434    vpCurrentMatrix( ar->vpc , VP_VIEW ) ;
00435    vpIdentityMatrix( ar->vpc ) ;
00436    vpRotate( ar->vpc , VP_Y_AXIS , theta ) ;
00437    vpRotate( ar->vpc , VP_X_AXIS , phi   ) ;
00438 
00439    if( ar->verbose )
00440       fprintf(stderr,"--AREN: set theta=%f  phi=%f\n",theta,phi) ;
00441 
00442    return ;
00443 }
00444 
00445 /*-------------------------------------------------------------------------------
00446    Set the scale factors for each axis (default = 1).  The inputs should
00447    be the size of the brick along each axes (e.g., sx = nx * dx).  This is
00448    needed because VolPack assumes the input data is in a unit cube, but
00449    our data may not be so uniform.
00450 
00451    N.B.: This is disabled, since VolPack doesn't seem to work in this regards.
00452 ---------------------------------------------------------------------------------*/
00453 
00454 void AREN_set_size( void * ah , float sx , float sy , float sz )
00455 {
00456 #if 0
00457    AREN_stuff * ar = (AREN_stuff *) ah ;
00458    float mmm ;
00459 
00460    if( !ISVALID_AREN(ar) ) return ;
00461 
00462    sx = fabs(sx) ; if( sx == 0.0 ) sx = 1.0 ;  /* don't allow non-positive sizes */
00463    sy = fabs(sy) ; if( sy == 0.0 ) sy = 1.0 ;
00464    sz = fabs(sz) ; if( sz == 0.0 ) sz = 1.0 ;
00465 
00466    mmm = sx ;
00467    if( mmm < sy ) mmm = sy ;
00468    if( mmm < sz ) mmm = sz ;  /* mmm = maximum size */
00469 
00470    ar->sx = sx / mmm ;        /* scale factors are <= 1.0 */
00471    ar->sy = sy / mmm ;
00472    ar->sz = sz / mmm ;
00473 
00474    vpCurrentMatrix( ar->vpc , VP_MODEL ) ;  /* scale model to world */
00475    vpIdentityMatrix( ar->vpc ) ;
00476    vpScale( ar->vpc , sx , sy , sz ) ;
00477 
00478    if( ar->verbose )
00479       fprintf(stderr,"--AREN: set scale factors = %f %f %f\n",ar->sx,ar->sy,ar->sz) ;
00480 #endif
00481 
00482    return ;
00483 }
00484 
00485 #undef USE_CALLBACKS
00486 #ifdef USE_CALLBACKS
00487 /*------------------------------------------------------------------------------
00488    Routines that can be used for callback shading, vs. lookup table shading.
00489    Mostly for debugging purposes.
00490 --------------------------------------------------------------------------------*/
00491 
00492 static int ncall = 0 ;
00493 
00494 void AREN_grayfunc( void * vp , float * g , void * cd )
00495 {
00496    rgbvox * rv = (rgbvox *) vp ;
00497 
00498    *g = AREN_graytable[ rv->rgb ] ;
00499 
00500 #if 0
00501    ncall++ ;
00502    if( ncall%1000 == 0 && rv->rgb > 0 )
00503       fprintf(stderr,"Grayfunc: rgb=%d g=%f\n",(int)(rv->rgb) , *g ) ;
00504 #endif
00505 
00506    return ;
00507 }
00508 
00509 void AREN_colorfunc( void * vp , float * r , float * g , float * b , void * cd )
00510 {
00511    rgbvox * rv = (rgbvox *) vp ;
00512    *r = AREN_colortable[ 3*(rv->rgb)   ] ;
00513    *g = AREN_colortable[ 3*(rv->rgb)+1 ] ;
00514    *b = AREN_colortable[ 3*(rv->rgb)+2 ] ;
00515 
00516 #if 0
00517    ncall++ ;
00518    if( ncall%1000 == 0 && rv->rgb > 0 )
00519       fprintf(stderr,"Colorfunc: rgb=%d r=%f g=%f b=%f\n",
00520                      (int)(rv->rgb) , *r,*g,*b ) ;
00521 #endif
00522 
00523    return ;
00524 }
00525 #endif /* USE_CALLBACKS */
00526 
00527 /*-------------------------------------------------------------------------------
00528    Actually render an image.  Returns NULL if an error occurs.
00529    Input npix = number of pixels on a side of the image (always will be square).
00530    If rendering a grayscale brick, returns an image of kind MRI_byte.
00531    If rendering RGB bricks, returns an image of kind MRI_rgb.
00532 ---------------------------------------------------------------------------------*/
00533 
00534 MRI_IMAGE * AREN_render( void * ah , int npix )
00535 {
00536    AREN_stuff * ar = (AREN_stuff *) ah ;
00537    int isgray , isrgb ;
00538    MRI_IMAGE * im ;
00539    byte * imar ;
00540    vpResult fred ;
00541 
00542    /*-- sanity checks --*/
00543 
00544    if( !ISVALID_AREN(ar) ) return NULL ;
00545 
00546    if( npix < 16 ){
00547       fprintf(stderr,"**AREN: attempt to render with less than 16 pixels!\n") ;
00548       return NULL ;
00549    }
00550 
00551    isgray = (ar->grim != NULL) ;
00552    isrgb  = (ar->rim  != NULL) ;
00553 
00554    if( isgray && isrgb ){
00555       fprintf(stderr,"**AREN: attempt to render gray and color simultaneously?\n");
00556       return NULL ;
00557    }
00558 
00559    if( (!isgray && !isrgb) || ar->vox == NULL ){
00560       fprintf(stderr,"**AREN: attempt to render without data being loaded!\n") ;
00561       return NULL ;
00562    }
00563 
00564    if( ar->opim == NULL ){
00565       fprintf(stderr,"**AREN: attempt to render without opacity being loaded!\n") ;
00566       return NULL ;
00567    }
00568 
00569    if( ar->nx < 3 || ar->ny < 3 || ar->nz < 3 ){
00570       fprintf(stderr,"**AREN: attempt to render without initialization!\n") ;
00571       return NULL ;
00572    }
00573 
00574    /*-- if have new voxel array, must tell VolPack all about it --*/
00575 
00576    if( ar->newvox ){
00577       int nvox = ar->nx * ar->ny * ar->nz ;
00578       rgbvox vv , *rv = &vv ;
00579 
00580       if( ar->verbose )
00581          fprintf(stderr,"--AREN: setting up new voxel array\n") ;
00582 
00583       /* 3D dimensions */
00584 
00585       fred = vpSetVolumeSize( ar->vpc , ar->nx , ar->ny , ar->nz ) ;
00586       if( fred != VP_OK ){
00587          fprintf(stderr,"**AREN: vpSetVolumeSize failed: code=%d\n",(int)fred) ;
00588          return NULL ;
00589       }
00590 
00591       /* each voxel has 2 data fields; 1 for shading and 1 for opacity */
00592 
00593       fred = vpSetVoxelSize( ar->vpc , sizeof(rgbvox) , 2 , 1 , 1 ) ;
00594       if( fred != VP_OK ){
00595          fprintf(stderr,"**AREN: vpSetVoxelSize failed: code=%d\n",(int)fred) ;
00596          return NULL ;
00597       }
00598 
00599       /* gray scale input */
00600 
00601       if( isgray ){
00602 
00603          /* voxel field 0 (rgb) is an index into AREN_graytable */
00604 
00605          fred = vpSetVoxelField( ar->vpc, 0, sizeof(short),
00606                                  vpFieldOffset(rv,rgb), MAX_GRAYS-1 );
00607          if( fred != VP_OK ){
00608             fprintf(stderr,"**AREN: vpSetVoxelField(0) failed: code=%d\n",(int)fred) ;
00609             return NULL ;
00610          }
00611 
00612          /* voxel field 1 (alpha) is an index into AREN_opatable */
00613 
00614          fred = vpSetVoxelField( ar->vpc, 1, sizeof(short),
00615                                  vpFieldOffset(rv,alpha), MAX_GRAYS-1 );
00616          if( fred != VP_OK ){
00617             fprintf(stderr,"**AREN: vpSetVoxelField(1) failed: code=%d\n",(int)fred) ;
00618             return NULL ;
00619          }
00620 
00621          /* tell VolPack where the voxels are */
00622 
00623          fred = vpSetRawVoxels( ar->vpc , ar->vox ,
00624                                 sizeof(rgbvox)*nvox ,
00625                                 sizeof(rgbvox) ,
00626                                 sizeof(rgbvox)*(ar->nx) ,
00627                                 sizeof(rgbvox)*(ar->nx * ar->ny) ) ;
00628          if( fred != VP_OK ){
00629             fprintf(stderr,"**AREN: vpSetRawVoxels failed: code=%d\n",(int)fred) ;
00630             return NULL ;
00631          }
00632 
00633 #ifndef USE_CALLBACKS
00634 
00635          /* voxel field 0 (rgb) is an index into AREN_graytable */
00636 
00637          fred = vpSetLookupShader( ar->vpc , 1 , 1 , 0 ,
00638                                    AREN_graytable , sizeof(float)*MAX_GRAYS ,
00639                                    0 , NULL , 0 ) ;
00640          if( fred != VP_OK ){
00641             fprintf(stderr,"**AREN: vpSetLookupShader failed: code=%d\n",(int)fred) ;
00642             return NULL ;
00643          }
00644 #else
00645          fred = vpSetCallback( ar->vpc , VP_GRAY_SHADE_FUNC , AREN_grayfunc ) ;
00646          if( fred != VP_OK ){
00647             fprintf(stderr,"**AREN: vpSetCallback(GRAY) failed: code=%d\n",(int)fred) ;
00648             return NULL ;
00649          }
00650 #endif
00651 
00652       /* color input */
00653 
00654       } else if( isrgb ){
00655 
00656          /* voxel field 0 (rgb) is an index into AREN_colortable */
00657 
00658          fred = vpSetVoxelField( ar->vpc, 0, sizeof(short),
00659                                  vpFieldOffset(rv,rgb), TOT_COLORS-1 );
00660          if( fred != VP_OK ){
00661             fprintf(stderr,"**AREN: vpSetVoxelField(0) failed: code=%d\n",(int)fred) ;
00662             return NULL ;
00663          }
00664 
00665          /* voxel field 1 (alpha) is an index into AREN_opatable */
00666 
00667          fred = vpSetVoxelField( ar->vpc, 1, sizeof(short) ,
00668                                  vpFieldOffset(rv,alpha), MAX_GRAYS-1 );
00669          if( fred != VP_OK ){
00670             fprintf(stderr,"**AREN: vpSetVoxelField(1) failed: code=%d\n",(int)fred) ;
00671             return NULL ;
00672          }
00673 
00674          /* tell VolPack where the voxels are */
00675 
00676          fred = vpSetRawVoxels( ar->vpc , ar->vox ,
00677                                 sizeof(rgbvox)*nvox ,
00678                                 sizeof(rgbvox) ,
00679                                 sizeof(rgbvox)*(ar->nx) ,
00680                                 sizeof(rgbvox)*(ar->nx * ar->ny) ) ;
00681          if( fred != VP_OK ){
00682             fprintf(stderr,"**AREN: vpSetRawVoxels failed: code=%d\n",(int)fred) ;
00683             return NULL ;
00684          }
00685 
00686 #ifndef USE_CALLBACKS
00687 
00688          /* voxel field 0 (rgb) is an index into AREN_colortable */
00689 
00690          fred = vpSetLookupShader( ar->vpc , 3 , 1 , 0 ,
00691                                    AREN_colortable , sizeof(float)*TOT_COLORS*3 ,
00692                                    0 , NULL , 0 ) ;
00693          if( fred != VP_OK ){
00694             fprintf(stderr,"**AREN: vpSetLookupShader failed: code=%d\n",(int)fred) ;
00695             return NULL ;
00696          }
00697 #else
00698          fred = vpSetCallback( ar->vpc , VP_RGB_SHADE_FUNC , AREN_colorfunc ) ;
00699          if( fred != VP_OK ){
00700             fprintf(stderr,"**AREN: vpSetCallback(COLOR) failed: code=%d\n",(int)fred) ;
00701             return NULL ;
00702          }
00703 #endif
00704       }
00705 
00706       /* voxel field 1 (alpha) is an index into AREN_opatable */
00707 
00708       fred = vpSetClassifierTable( ar->vpc, 0, 1, AREN_opatable, sizeof(float)*MAX_GRAYS ) ;
00709       if( fred != VP_OK ){
00710          fprintf(stderr,"**AREN: vpSetClassifierTable failed: code=%d\n",(int)fred) ;
00711          return NULL ;
00712       }
00713 
00714       /* threshold for octree bins: 12 = 5% of opacity range */
00715 
00716       fred = vpMinMaxOctreeThreshold( ar->vpc , 0 , 12 ) ;
00717       if( fred != VP_OK ){
00718          fprintf(stderr,"**AREN: vpMinMaxOctreeThreshold failed: code=%d\n",(int)fred) ;
00719          return NULL ;
00720       }
00721 
00722       ar->newvox = 0 ;
00723    }
00724 
00725    /*-- if have new data in the voxel array, must recreate the octree --*/
00726 
00727 #if 1
00728    vpSetd( ar->vpc , VP_MAX_RAY_OPACITY   , 0.95 ) ;
00729    vpSetd( ar->vpc , VP_MIN_VOXEL_OPACITY , 0.05 ) ;
00730 #endif
00731 
00732    if( ar->newdata ){
00733 
00734 #undef USE_CLASSIFIED  /* for using the "classified" volume feature of VolPack */
00735 
00736 #ifndef USE_CLASSIFIED
00737       if( ar->verbose )
00738          fprintf(stderr,"--AREN: computing octree\n") ;
00739 
00740       /* make octree, down to 4x4x4 voxel bins */
00741 
00742       fred = vpCreateMinMaxOctree( ar->vpc , 0 , 4 ) ;
00743       if( fred != VP_OK ){
00744          fprintf(stderr,"**AREN: vpCreateMinMaxOctree failed: code=%d\n",(int)fred) ;
00745          return NULL ;
00746       }
00747 #else
00748       if( ar->verbose )
00749          fprintf(stderr,"--AREN: computing classified volume\n") ;
00750 
00751       /* classify volume (slower than octree-ifying) */
00752 
00753       fred = vpClassifyVolume( ar->vpc ) ;
00754       if( fred != VP_OK ){
00755          fprintf(stderr,"**AREN: vpClassifyVolume failed: code=%d\n",(int)fred) ;
00756          return NULL ;
00757       }
00758 #endif
00759 
00760       ar->newdata = 0 ;
00761    }
00762 
00763    /*-- create the output image --*/
00764 
00765 #undef GET_ALPHA  /* for debugging: compute the opacity image */
00766 
00767    if( isgray ){
00768       im   = mri_new( npix , npix , MRI_byte ) ;
00769       imar = MRI_BYTE_PTR(im) ;
00770 #ifndef GET_ALPHA
00771       vpSetImage( ar->vpc , imar , npix,npix,npix , VP_LUMINANCE ) ;
00772 #else
00773       vpSetImage( ar->vpc , imar , npix,npix,npix , VP_ALPHA ) ;
00774 #endif
00775    } else if( isrgb ){
00776 #ifndef GET_ALPHA
00777       im   = mri_new( npix , npix , MRI_rgb ) ;
00778       imar = MRI_RGB_PTR(im) ;
00779       vpSetImage( ar->vpc , imar , npix,npix,3*npix , VP_RGB ) ;
00780 #else
00781       im   = mri_new( npix , npix , MRI_byte ) ;
00782       imar = MRI_BYTE_PTR(im) ;
00783       vpSetImage( ar->vpc , imar , npix,npix,npix , VP_ALPHA ) ;
00784 #endif
00785    }
00786 
00787    if( ar->verbose )
00788       fprintf(stderr,"--AREN: rendering image\n") ;
00789 
00790 #ifdef USE_CLASSIFIED
00791    fred = vpRenderClassifiedVolume(ar->vpc) ;
00792    if( fred != VP_OK ){
00793       fprintf(stderr,"**AREN: vpRenderClassifiedVolume failed: code=%d\n",(int)fred) ;
00794       mri_free(im) ; return NULL ;
00795    }
00796 #else
00797    fred = vpRenderRawVolume(ar->vpc) ;
00798    if( fred != VP_OK ){
00799       fprintf(stderr,"**AREN: vpRenderRawVolume failed: code=%d\n",(int)fred) ;
00800       mri_free(im) ; return NULL ;
00801    }
00802 #endif
00803 
00804    return im ;
00805 }