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 "mri_render.h"
00008 
00009 /*============================================================================
00010   MREN = MRI Renderer, a set of routines for volume rendering 3D bricks.
00011   Built on top of VolPack, by Philippe Lacroute.
00012 
00013   (1) Use new_MREN_renderer() to create a renderer.
00014   (2) Load grayscale or color bricks into it with functions
00015       MREN_set_graybytes(), MREN_set_rgbbytes(), or MREN_set_rgbshorts.
00016   (3) Load the opacity brick with MREN_set_opabrick().
00017   (4) Set the viewing angles with MREN_set_viewpoint().
00018   (5) Create an image with MREN_render().
00019       Loop back through (2-4) as needed to create new images.
00020   (6) When finished, use destroy_MREN_renderer() to clean up.
00021 ==============================================================================*/
00022 
00023 /*--------------------------------------------------------------------------
00024   Initialize the short -> float color table, and
00025              the byte  -> float gray/opacity tables for VolPack.
00026 
00027   N.B.: The documentation for VolPack is incorrect -- the color values
00028         should be in the range 0..255, NOT 0..1.  The opacity values
00029         should still be in the range 0..1.
00030 ----------------------------------------------------------------------------*/
00031 
00032 static float * MREN_colorshorts = NULL ;
00033 static float * MREN_graytable   = NULL ;
00034 static float * MREN_opatable    = NULL ;
00035 static float * MREN_colorbytes  = NULL ;
00036 
00037 void init_MREN_colortable( void )
00038 {
00039    int ii , rr,gg,bb , ss ;
00040 
00041    if( MREN_colorshorts != NULL ) return ;  /* been here already */
00042 
00043    MREN_colorshorts = (float *) malloc( sizeof(float) * TOT_COLORS * 3 ) ;
00044    MREN_graytable   = (float *) malloc( sizeof(float) * MREN_MAX_GRAYS ) ;
00045    MREN_opatable    = (float *) malloc( sizeof(float) * MREN_MAX_GRAYS ) ;
00046    MREN_colorbytes  = (float *) malloc( sizeof(float) * MREN_MAX_GRAYS * 3 ) ;
00047 
00048    /*-- load linear ramp for grayscale and opacity --*/
00049 
00050    for( ii=0 ; ii < MREN_MAX_GRAYS ; ii++ ){
00051       MREN_graytable[ii] = ii ;
00052       MREN_opatable[ii]  = ii / 255.0 ;
00053    }
00054 
00055    /*-- load a 32 x 32 x 32 color cube [indexed by unsigned shorts] --*/
00056 
00057    for( rr=0 ; rr < MREN_MAX_CDIM ; rr++ ){
00058       for( gg=0 ; gg < MREN_MAX_CDIM ; gg++ ){
00059          for( bb=0 ; bb < MREN_MAX_CDIM ; bb++ ){
00060 
00061             ss = FIVE_TO_SHORT(rr,gg,bb) ;   /* color index */
00062 
00063             MREN_colorshorts[3*ss  ] = (rr * 255.0) / 31.0 ;
00064             MREN_colorshorts[3*ss+1] = (gg * 255.0) / 31.0 ;
00065             MREN_colorshorts[3*ss+2] = (bb * 255.0) / 31.0 ;
00066          }
00067       }
00068    }
00069 
00070    /*-- at the end, add the pure grays (at a higher resolution) --*/
00071 
00072    ss = 3 * MREN_MAX_COLORS ;
00073    for( ii=0 ; ii < MREN_MAX_GRAYS ; ii++ ){
00074       MREN_colorshorts[ss++] = ii ;
00075       MREN_colorshorts[ss++] = ii ;
00076       MREN_colorshorts[ss++] = ii ;
00077    }
00078 
00079    /*-- load a short 8 x 8 x 4 color cube [indexed by unsigned chars] --*/
00080 
00081    for( rr=0 ; rr < 8 ; rr++ ){
00082       for( gg=0 ; gg < 8 ; gg++ ){
00083          for( bb=0 ; bb < 4 ; bb++ ){
00084 
00085             ss = (rr << 5) | (gg << 2) || (bb) ;  /* color index */
00086 
00087             MREN_colorbytes[3*ss  ] = (rr * 255.0) / 8.0 ;
00088             MREN_colorbytes[3*ss+1] = (gg * 255.0) / 8.0 ;
00089             MREN_colorbytes[3*ss+2] = (bb * 255.0) / 4.0 ;
00090          }
00091       }
00092    }
00093 
00094    return ;
00095 }
00096 
00097 void destroy_MREN_colortable( void )
00098 {
00099    if( MREN_colorshorts == NULL ) return ;
00100    free( MREN_colorshorts ); MREN_colorshorts = NULL ;
00101    free( MREN_graytable   ); MREN_graytable   = NULL ;
00102    free( MREN_opatable    ); MREN_opatable    = NULL ;
00103    free( MREN_colorbytes  ); MREN_colorbytes  = NULL ;
00104    return ;
00105 }
00106 
00107 static int num_renderers = 0 ;  /* global count of how many are open */
00108 
00109 /*--------------------------------------------------------------------------
00110   Create a new MRI renderer.
00111   The return pointer is passed to all other MREN routines.
00112 ----------------------------------------------------------------------------*/
00113 
00114 #define DEFAULT_THETA 130.0
00115 #define DEFAULT_PHI   285.0
00116 #define DEFAULT_PSI     0.0
00117 
00118 void * new_MREN_renderer( void )
00119 {
00120    MREN_stuff * ar ;
00121 
00122    ar = (MREN_stuff *) malloc( sizeof(MREN_stuff) ) ;
00123    ar->type = MREN_TYPE ;
00124 
00125    init_MREN_colortable() ;  /* in case it's not already setup */
00126 
00127    /*-- initialize VolPack --*/
00128 
00129    ar->vpc = vpCreateContext() ;
00130 
00131    vpSeti( ar->vpc , VP_CONCAT_MODE , VP_CONCAT_LEFT ) ;
00132 
00133    vpCurrentMatrix( ar->vpc , VP_MODEL ) ;
00134    vpIdentityMatrix( ar->vpc ) ;
00135 
00136    vpCurrentMatrix( ar->vpc , VP_VIEW ) ;
00137    vpIdentityMatrix( ar->vpc ) ;
00138    vpRotate( ar->vpc , VP_X_AXIS , DEFAULT_PHI   ) ;
00139    vpRotate( ar->vpc , VP_Y_AXIS , DEFAULT_THETA ) ;
00140 
00141 #undef USE_CUEING
00142 #ifdef USE_CUEING
00143    vpSetDepthCueing( ar->vpc , 1.0 , 0.5 ) ;
00144    vpEnable( ar->vpc , VP_DEPTH_CUE , 1 ) ;
00145 #endif
00146 
00147    vpCurrentMatrix( ar->vpc , VP_PROJECT ) ;
00148    vpIdentityMatrix( ar->vpc ) ;
00149    vpWindow( ar->vpc , VP_PARALLEL , -0.55,0.55 , -0.55,0.55 , -0.55,0.55 ) ;
00150 
00151    /*-- initialize the rest of the data --*/
00152 
00153    ar->nx = ar->ny = ar->nz = ar->verbose = ar->newopac = ar->newvox = 0 ;
00154    ar->sx = ar->sy = ar->sz = 1.0 ;
00155 
00156    ar->theta = DEFAULT_THETA ;
00157    ar->phi   = DEFAULT_PHI ;
00158    ar->psi   = DEFAULT_PSI ;
00159    ar->shim  = ar->opim = NULL ;
00160    ar->vox   = NULL ;
00161    ar->pmode = PMODE_LOW ;
00162 
00163    ar->grayset = ar->rgbset = ar->opaset = 0 ;  /* nothing set yet */
00164 
00165    ar->ncmap = ar->newcmap = 0 ;
00166    ar->cmap  = NULL ;
00167 
00168    ar->min_opacity = 0.05 ;
00169 
00170    num_renderers ++ ;
00171    return (void *) ar ;
00172 }
00173 
00174 /*-----------------------------------------------------------------------------
00175   Turn depth cueing on or off -- 11 Sep 2001
00176 -------------------------------------------------------------------------------*/
00177 
00178 void MREN_depth_cue( void *ah , int onoff )
00179 {
00180    MREN_stuff * ar = (MREN_stuff *) ah ;
00181 
00182    if( !ISVALID_MREN(ar) ) return ;
00183 
00184    vpSetDepthCueing( ar->vpc , 2.0 , 1.3863 ) ;
00185    vpEnable( ar->vpc , VP_DEPTH_CUE , onoff ) ;
00186    return ;
00187 }
00188 
00189 /*-----------------------------------------------------------------------------
00190    Get rid of an MRI renderer.
00191 -------------------------------------------------------------------------------*/
00192 
00193 void destroy_MREN_renderer( void * ah )
00194 {
00195    MREN_stuff * ar = (MREN_stuff *) ah ;
00196 
00197    if( !ISVALID_MREN(ar) ) return ;
00198 
00199    if( ar->vox  != NULL ) free(ar->vox) ;
00200    if( ar->cmap != NULL ) free(ar->cmap) ;
00201    vpDestroyContext( ar->vpc ) ;
00202    free(ar) ;
00203 
00204    num_renderers-- ; if( num_renderers == 0 ) destroy_MREN_colortable() ;
00205    return ;
00206 }
00207 
00208 /*-----------------------------------------------------------------------------
00209    For debugging purposes
00210 -------------------------------------------------------------------------------*/
00211 
00212 void MREN_be_verbose( void * ah )
00213 {
00214    MREN_stuff * ar = (MREN_stuff *) ah ;
00215    if( !ISVALID_MREN(ar) ) return ;
00216    ar->verbose = 1 ; return ;
00217 }
00218 
00219 void MREN_be_quiet( void * ah )
00220 {
00221    MREN_stuff * ar = (MREN_stuff *) ah ;
00222    if( !ISVALID_MREN(ar) ) return ;
00223    ar->verbose = 0 ; return ;
00224 }
00225 
00226 /*-----------------------------------------------------------------------------
00227    Set the minimum voxel opacity to consider
00228 -------------------------------------------------------------------------------*/
00229 
00230 void MREN_set_min_opacity( void * ah , float opm )
00231 {
00232    MREN_stuff * ar = (MREN_stuff *) ah ;
00233 
00234    if( !ISVALID_MREN(ar) ) return ;
00235    if( opm <= 0.0 || opm >= 1.0 ) opm = 0.05 ;
00236    ar->min_opacity = opm ;
00237 
00238    if( ar->verbose ) fprintf(stderr,"--MREN: min_opacity = %f\n",opm) ;
00239    return ;
00240 }
00241 
00242 /*-----------------------------------------------------------------------------
00243   Load an user defined colormap.
00244     ncol    = number of colors (2..65535)
00245     rmap[i] = red map for index i=0..ncol-1 (values 0..255), etc.
00246 -------------------------------------------------------------------------------*/
00247 
00248 void MREN_set_rgbmap( void * ah, int ncol, byte * rmap, byte * gmap, byte * bmap )
00249 {
00250    MREN_stuff * ar = (MREN_stuff *) ah ;
00251    int ii ;
00252 
00253    if( !ISVALID_MREN(ar) ) return ;
00254    if( ncol < 2 || ncol > 65535 || rmap==NULL || gmap==NULL || bmap==NULL ) return ;
00255 
00256    if( ar->cmap != NULL ) free(ar->cmap) ;
00257 
00258    ar->cmap  = (float *) malloc( sizeof(float) * (3*ncol) ) ;
00259    ar->ncmap = ncol ;
00260 
00261    for( ii=0 ; ii < ncol ; ii++ ){
00262       ar->cmap[3*ii  ] = rmap[ii] ;
00263       ar->cmap[3*ii+1] = gmap[ii] ;
00264       ar->cmap[3*ii+2] = bmap[ii] ;
00265    }
00266 
00267    ar->newcmap = 1 ;
00268 
00269    if( ar->verbose ){
00270       fprintf(stderr,"--MREN: new colormap\n") ;
00271       for( ii=0 ; ii < ncol ; ii++ ){
00272          fprintf(stderr,"#%3d: %5.1f %5.1f %5.1f",
00273                  ii , ar->cmap[3*ii],ar->cmap[3*ii+1],ar->cmap[3*ii+2]) ;
00274          ii++ ;
00275          if( ii < ncol )
00276             fprintf(stderr,"  #%3d: %5.1f %5.1f %5.1f",
00277                     ii , ar->cmap[3*ii],ar->cmap[3*ii+1],ar->cmap[3*ii+2]) ;
00278          ii++ ;
00279          if( ii < ncol )
00280             fprintf(stderr,"  #%3d: %5.1f %5.1f %5.1f",
00281                     ii , ar->cmap[3*ii],ar->cmap[3*ii+1],ar->cmap[3*ii+2]) ;
00282          fprintf(stderr,"\n") ;
00283       }
00284    }
00285    return ;
00286 }
00287 
00288 void MREN_unset_rgbmap( void * ah )
00289 {
00290    MREN_stuff * ar = (MREN_stuff *) ah ;
00291 
00292    if( !ISVALID_MREN(ar) || ar->cmap == NULL ) return ;
00293    if( ar->cmap != NULL ){ free(ar->cmap) ; ar->cmap = NULL ; }
00294    ar->ncmap = 0 ; ar->newcmap = 1 ;
00295 
00296    if( ar->verbose ) fprintf(stderr,"--MREN: delete colormap\n") ;
00297    return ;
00298 }
00299 
00300 /*-----------------------------------------------------------------------------
00301    Set the grayscale brick in a renderer;
00302    Returns -1 if an error occurs, otherwise returns 0.
00303    Note that this brick is NOT free-d by MREN at any point -- that is
00304    the user's responsibility.
00305 -------------------------------------------------------------------------------*/
00306 
00307 int MREN_set_graybytes( void * ah , MRI_IMAGE * grim )
00308 {
00309    MREN_stuff * ar = (MREN_stuff *) ah ;
00310    int newvox=0 , nvox,ii ;
00311    byte    * gar ;
00312    rgbvox  * rvox ;
00313 
00314    /*-- sanity checks --*/
00315 
00316    if( !ISVALID_MREN(ar) || grim == NULL || grim->kind != MRI_byte ) return -1 ;
00317 
00318    if( grim->nx < 3 || grim->ny < 3 || grim->nz < 3 ){
00319       fprintf(stderr,"**MREN: illegal dimensions for a gray brick\n") ;
00320       return -1 ;
00321    }
00322 
00323    if( ar->verbose ){
00324       if( ar->rgbset ) fprintf(stderr,"--MREN: switching from rgb to gray brick\n") ;
00325       else             fprintf(stderr,"--MREN: input a new gray brick\n") ;
00326    }
00327 
00328    /*-- if have new dimensions, have to invalidate pre-existing opacity --*/
00329 
00330    if( ar->nx > 0 &&
00331        ( ar->nx != grim->nx || ar->ny != grim->ny || ar->nz != grim->nz ) ){
00332 
00333       ar->opim = NULL ; ar->opaset = 0 ;
00334 
00335       if( ar->vox != NULL ){ free(ar->vox) ; ar->vox = NULL ; }
00336 
00337       if( ar->verbose )
00338          fprintf(stderr,"--MREN: new gray brick changes volume dimensions\n"
00339                         "        nx:%d->%d  ny:%d->%d  nz:%d->%d\n",
00340                         ar->nx,grim->nx , ar->ny,grim->ny , ar->nz,grim->nz ) ;
00341    }
00342 
00343    /*-- set dimensions --*/
00344 
00345    ar->shim = grim ;
00346    ar->nx   = grim->nx ;
00347    ar->ny   = grim->ny ;
00348    ar->nz   = grim->nz ; nvox = ar->nx * ar->ny * ar->nz ;
00349 
00350    /*-- if need be, allocate a voxel array to hold the data --*/
00351 
00352    if( ar->vox == NULL ){
00353       ar->newvox = newvox = 1 ;
00354       ar->vox = (rgbvox *) malloc( sizeof(rgbvox) * nvox ) ;
00355       if( ar->vox == NULL ){
00356          fprintf(stderr,"**MREN: can't malloc workspace with new gray brick\n") ;
00357          return -1 ;
00358       } else if( ar->verbose ){
00359          fprintf(stderr,"--MREN: allocated new voxel array\n") ;
00360       }
00361    }
00362 
00363    /*-- copy grayscale data into voxel array --*/
00364 
00365    rvox = ar->vox ;
00366    gar  = MRI_BYTE_PTR(grim) ;
00367    for( ii=0 ; ii < nvox ; ii++ ) rvox[ii].rgb = (unsigned short) gar[ii] ;
00368 
00369    if( ar->rgbset ) ar->newvox = 1 ;  /* changed from color to gray */
00370 
00371    ar->grayset = 1 ; ar->rgbset = 0 ;
00372    return 0 ;
00373 }
00374 
00375 /*-----------------------------------------------------------------------------
00376    Set the opacity brick in a renderer.
00377    Returns -1 if an error occurs, otherwise returns 0.
00378    Note that this brick is NOT free-d by MREN at any point -- that is
00379    the user's responsibility.
00380 -------------------------------------------------------------------------------*/
00381 
00382 int MREN_set_opabytes( void * ah , MRI_IMAGE * opim )
00383 {
00384    MREN_stuff * ar = (MREN_stuff *) ah ;
00385    int nvox,ii , newvox=0 ;
00386    byte    * gar ;
00387    rgbvox  * rvox ;
00388 
00389    /*-- sanity checks --*/
00390 
00391    if( !ISVALID_MREN(ar) || opim == NULL || opim->kind != MRI_byte ) return -1 ;
00392 
00393    if( opim->nx < 3 || opim->ny < 3 || opim->nz < 3 ){
00394       fprintf(stderr,"**MREN: illegal dimensions for an opacity brick\n") ;
00395       return -1 ;
00396    }
00397 
00398    /*-- if have new dimensions, toss old stuff that doesn't match --*/
00399 
00400    if( ar->nx > 0 &&
00401        ( ar->nx != opim->nx || ar->ny != opim->ny || ar->nz != opim->nz ) ){
00402 
00403       ar->shim = NULL ; ar->grayset = ar->rgbset = 0 ;
00404 
00405       if( ar->vox != NULL ){ free(ar->vox) ; ar->vox = NULL ; }
00406 
00407       if( ar->verbose )
00408          fprintf(stderr,"--MREN: new opacity brick changes volume dimensions\n"
00409                         "        nx:%d->%d  ny:%d->%d  nz:%d->%d\n",
00410                         ar->nx,opim->nx , ar->ny,opim->ny , ar->nz,opim->nz ) ;
00411    } else {
00412       if( ar->verbose ) fprintf(stderr,"--MREN: new opacity brick\n") ;
00413    }
00414 
00415    /*-- set dimensions --*/
00416 
00417    ar->opim = opim ;
00418    ar->nx   = opim->nx ;
00419    ar->ny   = opim->ny ;
00420    ar->nz   = opim->nz ; nvox = ar->nx * ar->ny * ar->nz ;
00421 
00422    /*-- if need be, allocate a voxel array to hold the data --*/
00423 
00424    if( ar->vox == NULL ){
00425       ar->newvox = newvox = 1 ;
00426       ar->vox = (rgbvox *) malloc( sizeof(rgbvox) * nvox ) ;
00427       if( ar->vox == NULL ){
00428          fprintf(stderr,"**MREN: can't malloc workspace with new opacity brick\n") ;
00429          return -1 ;
00430       } else if( ar->verbose ){
00431          fprintf(stderr,"--MREN: allocated new voxel array\n") ;
00432       }
00433    }
00434 
00435    /*-- load the opacity into voxel array --*/
00436 
00437    gar  = MRI_BYTE_PTR(ar->opim) ;
00438    rvox = ar->vox ;
00439    for( ii=0 ; ii < nvox ; ii++ ) rvox[ii].alpha = (unsigned short) gar[ii] ;
00440 
00441    ar->newopac = 1 ; ar->opaset = 1 ;
00442    return 0 ;
00443 }
00444 
00445 /*-----------------------------------------------------------------------------
00446    Convert an RGB image to a standard 8x8x4 (MREN_colorbytes)
00447    color-mapped image, suitable for use in MREN_set_rgbbytes.
00448    The output image is composed of bytes.
00449 -------------------------------------------------------------------------------*/
00450 
00451 MRI_IMAGE * MREN_rgb_to_colorbytes( MRI_IMAGE * rgbim )
00452 {
00453    byte * rgbar , rb,gb,bb ;
00454    byte * shar ;
00455    MRI_IMAGE * shim ;
00456    int ii ;
00457 
00458    if( rgbim == NULL || rgbim->kind != MRI_rgb ) return NULL ;
00459 
00460    shim  = mri_new_conforming( rgbim , MRI_byte ) ;
00461    shar  = MRI_BYTE_PTR(shim) ;
00462    rgbar = MRI_RGB_PTR(rgbim) ;
00463 
00464    for( ii=0 ; ii < shim->nvox ; ii++ ){
00465       rb = rgbar[3*ii  ] >> 5 ;
00466       gb = rgbar[3*ii+1] >> 5 ;
00467       bb = rgbar[3*ii+2] >> 6 ;
00468 
00469       shar[ii] = (rb << 5) | (gb << 2) | bb ;  /* index into colorbytes */
00470    }
00471 
00472    return shim ;
00473 }
00474 
00475 /*-----------------------------------------------------------------------------
00476    Convert an RGB image to a standard 32x32x32+256 (MREN_colorshorts)
00477    color-mapped image, suitable for use in MREN_set_rgbshorts.
00478    The output image is composed of shorts.
00479 -------------------------------------------------------------------------------*/
00480 
00481 MRI_IMAGE * MREN_rgb_to_colorshorts( MRI_IMAGE * rgbim )
00482 {
00483    byte * rgbar , rb,gb,bb ;
00484    unsigned short * shar ;
00485    MRI_IMAGE * shim ;
00486    int ii ;
00487 
00488    if( rgbim == NULL || rgbim->kind != MRI_rgb ) return NULL ;
00489 
00490    shim  = mri_new_conforming( rgbim , MRI_short ) ;
00491    shar  = (unsigned short *) MRI_SHORT_PTR(shim) ;
00492    rgbar = MRI_RGB_PTR(rgbim) ;
00493 
00494    for( ii=0 ; ii < shim->nvox ; ii++ ){
00495       rb = EIGHT_TO_FIVE(rgbar[3*ii  ]) ;
00496       gb = EIGHT_TO_FIVE(rgbar[3*ii+1]) ;
00497       bb = EIGHT_TO_FIVE(rgbar[3*ii+2]) ;
00498 
00499       if( rb == gb && rb == bb ){
00500          shar[ii] = MREN_MAX_COLORS + rgbar[3*ii] ; /* index into grayscale */
00501       } else {
00502          shar[ii] = FIVE_TO_SHORT( rb , gb, bb ) ;  /* index into color cube */
00503       }
00504    }
00505 
00506    return shim ;
00507 }
00508 
00509 /*-----------------------------------------------------------------------------
00510    Set the color brick in a renderer -- values in the input image are
00511    indices into a colormap of length <= 256 -- by default this will
00512    be MREN_colorbytes, but can be replaced by the user.
00513    Returns -1 if an error occurs, otherwise returns 0.
00514    Note that these bricks are NOT free-d by MREN at any point -- that is
00515    the user's responsibility.
00516 -------------------------------------------------------------------------------*/
00517 
00518 int MREN_set_rgbbytes( void * ah , MRI_IMAGE * rgbim )
00519 {
00520    MREN_stuff * ar = (MREN_stuff *) ah ;
00521    int newvox=0 , nvox,ii ;
00522    byte    * gar ;
00523    rgbvox  * rvox ;
00524 
00525    /*-- sanity checks --*/
00526 
00527    if( !ISVALID_MREN(ar) || rgbim == NULL || rgbim->kind != MRI_byte ) return -1 ;
00528 
00529    if( rgbim->nx < 3 || rgbim->ny < 3 || rgbim->nz < 3 ){
00530       fprintf(stderr,"**MREN: illegal dimensions for a color brick\n") ; return -1 ;
00531    }
00532 
00533    /*-- if had an old gray brick, toss it (or at least its pointer) --*/
00534 
00535    if( ar->verbose ){
00536       if( ar->grayset ) fprintf(stderr,"--MREN: switching from gray to rgb brick\n") ;
00537       else              fprintf(stderr,"--MREN: input new rgb brick of bytes\n") ;
00538    }
00539 
00540    /*-- if have new dimensions, toss old stuff that doesn't match --*/
00541 
00542    if( ar->nx > 0 &&
00543        ( ar->nx != rgbim->nx || ar->ny != rgbim->ny || ar->nz != rgbim->nz ) ){
00544 
00545       ar->opim = NULL ; ar->opaset = 0 ;
00546 
00547       if( ar->vox != NULL ){ free(ar->vox) ; ar->vox = NULL ; }
00548 
00549       if( ar->verbose )
00550          fprintf(stderr,"--MREN: new rgb brick changes volume dimensions\n"
00551                         "        nx:%d->%d  ny:%d->%d  nz:%d->%d\n",
00552                         ar->nx,rgbim->nx , ar->ny,rgbim->ny , ar->nz,rgbim->nz ) ;
00553    }
00554 
00555    /*-- set dimensions --*/
00556 
00557    ar->shim = rgbim ;
00558    ar->nx   = rgbim->nx ;
00559    ar->ny   = rgbim->ny ;
00560    ar->nz   = rgbim->nz ; nvox = ar->nx * ar->ny * ar->nz ;
00561 
00562    /*-- if need be, allocate a voxel array to hold the data --*/
00563 
00564    if( ar->vox == NULL ){
00565       ar->newvox = newvox = 1 ;
00566       ar->vox = (rgbvox *) malloc( sizeof(rgbvox) * nvox ) ;
00567       if( ar->vox == NULL ){
00568          fprintf(stderr,"**MREN: can't malloc workspace with new color bricks\n") ;
00569          return -1 ;
00570       } else if( ar->verbose ){
00571          fprintf(stderr,"--MREN: allocated new voxel array\n") ;
00572       }
00573    }
00574 
00575    /*-- copy color data into voxel array --*/
00576 
00577    rvox = ar->vox ;
00578    gar  = MRI_BYTE_PTR(rgbim) ;
00579    for( ii=0 ; ii < nvox ; ii++ ) rvox[ii].rgb = (unsigned short) gar[ii] ;
00580 
00581    if( ar->grayset ) ar->newvox = 1 ;  /* changed from gray to color */
00582 
00583    ar->rgbset = 1 ; ar->grayset = 0 ;
00584    return 0 ;
00585 }
00586 
00587 /*-----------------------------------------------------------------------------
00588   Similar, but where the colormap indexes are unsigned shorts.
00589   They can refer to a user-supplied colormap, or to the standard
00590   32x32x32+256 colormap (MREN_colorshorts).
00591 -------------------------------------------------------------------------------*/
00592 
00593 int MREN_set_rgbshorts( void * ah , MRI_IMAGE * rgbim )
00594 {
00595    MREN_stuff * ar = (MREN_stuff *) ah ;
00596    int newvox=0 , nvox,ii ;
00597    unsigned short * gar ;
00598    rgbvox * rvox ;
00599 
00600    /*-- sanity checks --*/
00601 
00602    if( !ISVALID_MREN(ar) || rgbim == NULL || rgbim->kind != MRI_short ) return -1 ;
00603 
00604    if( rgbim->nx < 3 || rgbim->ny < 3 || rgbim->nz < 3 ){
00605       fprintf(stderr,"**MREN: illegal dimensions for a color brick\n") ; return -1 ;
00606    }
00607 
00608    if( ar->verbose ){
00609       if( ar->grayset ) fprintf(stderr,"--MREN: switching from gray to rgb brick\n") ;
00610       else              fprintf(stderr,"--MREN: input new rgb brick of shorts\n") ;
00611    }
00612 
00613    /*-- if have new dimensions, toss old stuff that doesn't match --*/
00614 
00615    if( ar->nx > 0 &&
00616        ( ar->nx != rgbim->nx || ar->ny != rgbim->ny || ar->nz != rgbim->nz ) ){
00617 
00618       ar->opim = NULL ; ar->opaset = 0 ;
00619 
00620       if( ar->vox != NULL ){ free(ar->vox) ; ar->vox = NULL ; }
00621 
00622       if( ar->verbose )
00623          fprintf(stderr,"--MREN: new rgb brick changes volume dimensions\n"
00624                         "        nx:%d->%d  ny:%d->%d  nz:%d->%d\n",
00625                         ar->nx,rgbim->nx , ar->ny,rgbim->ny , ar->nz,rgbim->nz ) ;
00626    }
00627 
00628    /*-- set dimensions --*/
00629 
00630    ar->shim = rgbim ;
00631    ar->nx   = rgbim->nx ;
00632    ar->ny   = rgbim->ny ;
00633    ar->nz   = rgbim->nz ; nvox = ar->nx * ar->ny * ar->nz ;
00634 
00635    /*-- if need be, allocate a voxel array to hold the data --*/
00636 
00637    if( ar->vox == NULL ){
00638       ar->newvox = newvox = 1 ;
00639       ar->vox = (rgbvox *) malloc( sizeof(rgbvox) * nvox ) ;
00640       if( ar->vox == NULL ){
00641          fprintf(stderr,"**MREN: can't malloc workspace with new color bricks\n") ;
00642          return -1 ;
00643       } else if( ar->verbose ){
00644          fprintf(stderr,"--MREN: allocated new voxel array\n") ;
00645       }
00646    }
00647 
00648    /*-- copy color data into voxel array --*/
00649 
00650    rvox = ar->vox ;
00651    gar  = (unsigned short *) MRI_SHORT_PTR(rgbim) ;
00652    for( ii=0 ; ii < nvox ; ii++ ) rvox[ii].rgb = gar[ii] ;
00653 
00654    if( ar->grayset ) ar->newvox = 1 ;  /* changed from gray to color */
00655 
00656    ar->rgbset = 2 ; ar->grayset = 0 ;
00657    return 0 ;
00658 }
00659 
00660 /*-------------------------------------------------------------------------------
00661    Set the viewpoint of the user (in polar angles)
00662 ---------------------------------------------------------------------------------*/
00663 
00664 void MREN_set_viewpoint( void * ah , float theta , float phi , float psi )
00665 {
00666    MREN_stuff * ar = (MREN_stuff *) ah ;
00667 
00668    if( !ISVALID_MREN(ar) ) return ;
00669 
00670    ar->theta = theta ; ar->phi = phi ; ar->psi = psi ;
00671 
00672    vpCurrentMatrix( ar->vpc , VP_VIEW ) ;
00673    vpIdentityMatrix( ar->vpc ) ;
00674    vpRotate( ar->vpc , VP_Z_AXIS , psi   ) ;  /* roll  */
00675    vpRotate( ar->vpc , VP_X_AXIS , phi   ) ;  /* pitch */
00676    vpRotate( ar->vpc , VP_Y_AXIS , theta ) ;  /* yaw   */
00677 
00678    if( ar->verbose ){
00679       vpMatrix4 vpm ;
00680 
00681       fprintf(stderr,"--MREN: set theta=%f  phi=%f  psi=%f\n",theta,phi,psi) ;
00682 
00683       vpGetMatrix( ar->vpc , VP_VIEW , vpm ) ;
00684       fprintf(stderr,"--matrix: %8.5f %8.5f %8.5f %8.5f\n"
00685                      "          %8.5f %8.5f %8.5f %8.5f\n"
00686                      "          %8.5f %8.5f %8.5f %8.5f\n"
00687                      "          %8.5f %8.5f %8.5f %8.5f\n" ,
00688               vpm[0][0] , vpm[0][1] , vpm[0][2] , vpm[0][3] ,
00689               vpm[1][0] , vpm[1][1] , vpm[1][2] , vpm[1][3] ,
00690               vpm[2][0] , vpm[2][1] , vpm[2][2] , vpm[2][3] ,
00691               vpm[3][0] , vpm[3][1] , vpm[3][2] , vpm[3][3]  ) ;
00692    }
00693 
00694    return ;
00695 }
00696 
00697 /*-------------------------------------------------------------------------------
00698    Set the precalculation mode of the renderer
00699 ---------------------------------------------------------------------------------*/
00700 
00701 void MREN_set_precalculation( void * ah , int mode )
00702 {
00703    MREN_stuff * ar = (MREN_stuff *) ah ;
00704 
00705    if( !ISVALID_MREN(ar) || mode < PMODE_LOW || mode > PMODE_HIGH ) return ;
00706 
00707    if( ar->pmode != mode ){ ar->pmode = mode ; ar->newopac = 1 ; }
00708    return ;
00709 }
00710 
00711 /*-------------------------------------------------------------------------------
00712    Set the scale factors for each axis (default = 1).  The inputs should
00713    be the size of the brick along each axes (e.g., sx = nx * dx).  This is
00714    needed because VolPack assumes the input data is in a unit cube, but
00715    our data may not be so uniform.
00716 
00717    N.B.: This is disabled, since VolPack doesn't seem to work in this regards.
00718 ---------------------------------------------------------------------------------*/
00719 
00720 void MREN_set_size( void * ah , float sx , float sy , float sz )
00721 {
00722 #if 0
00723    MREN_stuff * ar = (MREN_stuff *) ah ;
00724    float mmm ;
00725 
00726    if( !ISVALID_MREN(ar) ) return ;
00727 
00728    sx = fabs(sx) ; if( sx == 0.0 ) sx = 1.0 ;  /* don't allow non-positive sizes */
00729    sy = fabs(sy) ; if( sy == 0.0 ) sy = 1.0 ;
00730    sz = fabs(sz) ; if( sz == 0.0 ) sz = 1.0 ;
00731 
00732    mmm = sx ;
00733    if( mmm < sy ) mmm = sy ;
00734    if( mmm < sz ) mmm = sz ;  /* mmm = maximum size */
00735 
00736    ar->sx = sx / mmm ;        /* scale factors are <= 1.0 */
00737    ar->sy = sy / mmm ;
00738    ar->sz = sz / mmm ;
00739 
00740    vpCurrentMatrix( ar->vpc , VP_MODEL ) ;  /* scale model to world */
00741    vpIdentityMatrix( ar->vpc ) ;
00742    vpScale( ar->vpc , sx , sy , sz ) ;
00743 
00744    if( ar->verbose )
00745       fprintf(stderr,"--MREN: set scale factors = %f %f %f\n",ar->sx,ar->sy,ar->sz) ;
00746 #endif
00747 
00748    return ;
00749 }
00750 
00751 /*-------------------------------------------------------------------------------
00752    Find out if a renderer needs data
00753 ---------------------------------------------------------------------------------*/
00754 
00755 int MREN_needs_data( void * ah )
00756 {
00757    MREN_stuff * ar = (MREN_stuff *) ah ;
00758 
00759    return (ar->vox == NULL) ;
00760 }
00761 
00762 /*-------------------------------------------------------------------------------
00763    Actually render an image.  Returns NULL if an error occurs.
00764    Input npix = number of pixels on a side of the image (always will be square).
00765    If rendering a grayscale brick, returns an image of kind MRI_byte.
00766    If rendering RGB bricks, returns an image of kind MRI_rgb.
00767 ---------------------------------------------------------------------------------*/
00768 
00769 MRI_IMAGE * MREN_render( void * ah , int npix )
00770 {
00771    MREN_stuff * ar = (MREN_stuff *) ah ;
00772    int isgray , isrgb ;
00773    MRI_IMAGE * im ;
00774    byte * imar ;
00775    vpResult fred ;
00776 
00777    /*-- sanity checks --*/
00778 
00779    if( !ISVALID_MREN(ar) ) return NULL ;
00780 
00781    if( npix < 16 ){
00782       fprintf(stderr,"**MREN: attempt to render with less than 16 pixels!\n") ;
00783       return NULL ;
00784    }
00785 
00786    isgray = (ar->grayset > 0) ;
00787    isrgb  = (ar->rgbset  > 0) ;
00788 
00789    if( isgray && isrgb ){
00790       fprintf(stderr,"**MREN: attempt to render gray and color simultaneously?\n");
00791       return NULL ;
00792    }
00793 
00794    if( (!isgray && !isrgb) || ar->vox == NULL ){
00795       fprintf(stderr,"**MREN: attempt to render without data being loaded!\n") ;
00796       return NULL ;
00797    }
00798 
00799    if( ar->opaset == 0 ){
00800       fprintf(stderr,"**MREN: attempt to render without opacity being loaded!\n") ;
00801       return NULL ;
00802    }
00803 
00804    if( ar->nx < 3 || ar->ny < 3 || ar->nz < 3 ){
00805       fprintf(stderr,"**MREN: attempt to render without initialization!\n") ;
00806       return NULL ;
00807    }
00808 
00809    if( ar->verbose ) fprintf(stderr,"--MREN: setup for rendering\n") ;
00810 
00811    /*-- if have new voxel array, must tell VolPack all about it --*/
00812 
00813    if( ar->newvox || (isrgb && ar->newcmap) ){
00814       int nvox = ar->nx * ar->ny * ar->nz ;
00815       int oct_range ;
00816       rgbvox vv , *rv = &vv ;
00817 
00818       /* 3D dimensions */
00819 
00820       if( ar->verbose ) fprintf(stderr,"  call vpSetVolumeSize\n") ;
00821 
00822       fred = vpSetVolumeSize( ar->vpc , ar->nx , ar->ny , ar->nz ) ;
00823       if( fred != VP_OK ){
00824          fprintf(stderr,"**MREN: vpSetVolumeSize failed: code=%d\n",(int)fred) ;
00825          return NULL ;
00826       }
00827 
00828       /* each voxel has 2 data fields; 1 for shading and 1 for opacity */
00829 
00830       if( ar->verbose ) fprintf(stderr,"  call vpSetVoxelSize\n") ;
00831 
00832       fred = vpSetVoxelSize( ar->vpc , sizeof(rgbvox) , 2 , 1 , 1 ) ;
00833       if( fred != VP_OK ){
00834          fprintf(stderr,"**MREN: vpSetVoxelSize failed: code=%d\n",(int)fred) ;
00835          return NULL ;
00836       }
00837 
00838       /* voxel field 1 (alpha) is an index into MREN_opatable */
00839 
00840       if( ar->verbose ) fprintf(stderr,"  call vpSetVoxelField(1)\n") ;
00841 
00842       fred = vpSetVoxelField( ar->vpc, 1, sizeof(short),
00843                               vpFieldOffset(rv,alpha), MREN_MAX_GRAYS-1 );
00844       if( fred != VP_OK ){
00845          fprintf(stderr,"**MREN: vpSetVoxelField(1) failed: code=%d\n",(int)fred) ;
00846          return NULL ;
00847       }
00848 
00849       /* tell VolPack where the voxels are */
00850 
00851       if( ar->verbose ) fprintf(stderr,"  call vpSetRawVoxels\n") ;
00852 
00853       fred = vpSetRawVoxels( ar->vpc , ar->vox ,
00854                              sizeof(rgbvox)*nvox ,
00855                              sizeof(rgbvox) ,
00856                              sizeof(rgbvox)*(ar->nx) ,
00857                              sizeof(rgbvox)*(ar->nx * ar->ny) ) ;
00858       if( fred != VP_OK ){
00859          fprintf(stderr,"**MREN: vpSetRawVoxels failed: code=%d\n",(int)fred) ;
00860          return NULL ;
00861       }
00862 
00863       /*--- gray scale input ---*/
00864 
00865       if( isgray ){
00866 
00867          /* voxel field 0 (rgb) is an index into MREN_graytable */
00868 
00869          if( ar->verbose ) fprintf(stderr,"  call vpSetVoxelField(grays)\n") ;
00870 
00871          fred = vpSetVoxelField( ar->vpc, 0, sizeof(short),
00872                                  vpFieldOffset(rv,rgb), MREN_MAX_GRAYS-1 );
00873          if( fred != VP_OK ){
00874             fprintf(stderr,"**MREN: vpSetVoxelField(0) failed: code=%d\n",(int)fred) ;
00875             return NULL ;
00876          }
00877 
00878          /* setup MREN_graytable to hold the colormap */
00879 
00880          if( ar->verbose ) fprintf(stderr,"  call vpSetLookupShader(graytable)\n") ;
00881 
00882          fred = vpSetLookupShader( ar->vpc , 1 , 1 , 0 ,
00883                                    MREN_graytable , sizeof(float)*MREN_MAX_GRAYS ,
00884                                    0 , NULL , 0 ) ;
00885          if( fred != VP_OK ){
00886             fprintf(stderr,"**MREN: vpSetLookupShader failed: code=%d\n",(int)fred) ;
00887             return NULL ;
00888          }
00889 
00890       /*--- color input ---*/
00891 
00892       } else if( isrgb ){
00893 
00894          /* There are 3 possible cases for the colormap:
00895               a) user supplied colormap
00896               b) standard 16 bit index colormap MREN_colorshorts
00897               c) standard 8 bit index colormap MREN_colorbytes    */
00898 
00899          if( ar->cmap != NULL && ar->ncmap > 1 ){   /* user supplied colormap */
00900 
00901             if( ar->verbose ) fprintf(stderr,"  call vpSetVoxelField(cmap)\n") ;
00902 
00903             fred = vpSetVoxelField( ar->vpc, 0, sizeof(short),
00904                                     vpFieldOffset(rv,rgb), ar->ncmap-1 );
00905 
00906             if( ar->verbose ) fprintf(stderr,"  call vpSetLookupShader(cmap)\n") ;
00907 
00908             fred = vpSetLookupShader( ar->vpc , 3 , 1 , 0 ,
00909                                       ar->cmap , sizeof(float)*3*ar->ncmap ,
00910                                       0 , NULL , 0 ) ;
00911 
00912          } else if( ar->rgbset == 2 ){          /* MREN_colorshorts */
00913 
00914             if( ar->verbose ) fprintf(stderr,"  call vpSetVoxelField(rgb shorts)\n") ;
00915 
00916             fred = vpSetVoxelField( ar->vpc, 0, sizeof(short),
00917                                     vpFieldOffset(rv,rgb), TOT_COLORS-1 );
00918 
00919             if( ar->verbose ) fprintf(stderr,"  call vpSetLookupShader(colorshorts)\n") ;
00920 
00921             fred = vpSetLookupShader( ar->vpc , 3 , 1 , 0 ,
00922                                       MREN_colorshorts , sizeof(float)*TOT_COLORS*3 ,
00923                                       0 , NULL , 0 ) ;
00924 
00925          } else {                               /* MREN_colorbytes */
00926 
00927             if( ar->verbose ) fprintf(stderr,"  call vpSetVoxelField(rgb bytes)\n") ;
00928 
00929             fred = vpSetVoxelField( ar->vpc, 0, sizeof(short),
00930                                     vpFieldOffset(rv,rgb), MREN_MAX_GRAYS-1 );
00931 
00932             if( ar->verbose ) fprintf(stderr,"  call vpSetLookupShader(colorbytes)\n") ;
00933 
00934             fred = vpSetLookupShader( ar->vpc , 3 , 1 , 0 ,
00935                                       MREN_colorbytes , sizeof(float)*MREN_MAX_GRAYS*3 ,
00936                                       0 , NULL , 0 ) ;
00937          }
00938 
00939          if( fred != VP_OK ){
00940             fprintf(stderr,"**MREN: vpSetLookupShader failed: code=%d\n",(int)fred) ;
00941             return NULL ;
00942          }
00943       }
00944 
00945       /*-- in all cases, voxel field 1 (alpha) is an index into MREN_opatable --*/
00946 
00947       if( ar->verbose ) fprintf(stderr,"  call vpSetClassifierTable\n") ;
00948 
00949       fred = vpSetClassifierTable( ar->vpc, 0, 1, MREN_opatable, sizeof(float)*MREN_MAX_GRAYS ) ;
00950       if( fred != VP_OK ){
00951          fprintf(stderr,"**MREN: vpSetClassifierTable failed: code=%d\n",(int)fred) ;
00952          return NULL ;
00953       }
00954 
00955       /* threshold for octree bins: 12 = 5% of possible opacity range */
00956 
00957       if( ar->verbose ) fprintf(stderr,"  call vpMinMaxOctreeThreshold\n") ;
00958 
00959       fred = vpMinMaxOctreeThreshold( ar->vpc , 0 , 12 ) ;
00960       if( fred != VP_OK ){
00961          fprintf(stderr,"**MREN: vpMinMaxOctreeThreshold failed: code=%d\n",(int)fred) ;
00962          return NULL ;
00963       }
00964 
00965       ar->newopac = 1 ;  /* for the precalculations below */
00966       ar->newvox  = 0 ;
00967       ar->newcmap = 0 ;
00968    }
00969 
00970    /*-- if have new data in the voxel array, must do precalculations --*/
00971 
00972    vpSetd( ar->vpc , VP_MAX_RAY_OPACITY   , 0.95 ) ;
00973    vpSetd( ar->vpc , VP_MIN_VOXEL_OPACITY , ar->min_opacity ) ;
00974 
00975    if( ar->newopac ){
00976 
00977       (void) vpDestroyMinMaxOctree( ar->vpc ) ;      /* toss previous work */
00978       (void) vpDestroyClassifiedVolume( ar->vpc ) ;
00979 
00980       if( ar->pmode == PMODE_MEDIUM ){
00981 
00982          if( ar->verbose ) fprintf(stderr,"--MREN: computing octree\n") ;
00983 
00984          /* make octree, down to 4x4x4 voxel bins */
00985 
00986          if( ar->verbose ) fprintf(stderr,"  call vpCreateMinMaxOctree\n") ;
00987 
00988          fred = vpCreateMinMaxOctree( ar->vpc , 0 , 4 ) ;
00989          if( fred != VP_OK ){
00990             fprintf(stderr,"**MREN: vpCreateMinMaxOctree failed: code=%d\n",(int)fred) ;
00991             return NULL ;
00992          }
00993 
00994       } else if( ar->pmode == PMODE_HIGH ){
00995 
00996          if( ar->verbose ) fprintf(stderr,"--MREN: computing classified volume\n") ;
00997 
00998          /* classify volume (slower than octree, but may do faster rendering) */
00999 
01000          if( ar->verbose ) fprintf(stderr,"  call vpClassifyVolume\n") ;
01001 
01002          fred = vpClassifyVolume( ar->vpc ) ;
01003          if( fred != VP_OK ){
01004             fprintf(stderr,"**MREN: vpClassifyVolume failed: code=%d\n",(int)fred) ;
01005             return NULL ;
01006          }
01007       }
01008 
01009       ar->newopac = 0 ;
01010    }
01011 
01012    /*-- create the output image --*/
01013 
01014 #undef GET_ALPHA  /* for debugging: compute the opacity image */
01015 
01016    if( isgray ){
01017       im   = mri_new( npix , npix , MRI_byte ) ;
01018       imar = MRI_BYTE_PTR(im) ;
01019 #ifndef GET_ALPHA
01020       if( ar->verbose ) fprintf(stderr,"  call vpSetImage(LUMINANCE)\n") ;
01021       vpSetImage( ar->vpc , imar , npix,npix,npix , VP_LUMINANCE ) ;
01022 #else
01023       if( ar->verbose ) fprintf(stderr,"  call vpSetImage(ALPHA)\n") ;
01024       vpSetImage( ar->vpc , imar , npix,npix,npix , VP_ALPHA ) ;
01025 #endif
01026    } else if( isrgb ){
01027 #ifndef GET_ALPHA
01028       im   = mri_new( npix , npix , MRI_rgb ) ;
01029       imar = MRI_RGB_PTR(im) ;
01030       if( ar->verbose ) fprintf(stderr,"  call vpSetImage(RGB)\n") ;
01031       vpSetImage( ar->vpc , imar , npix,npix,3*npix , VP_RGB ) ;
01032 #else
01033       im   = mri_new( npix , npix , MRI_byte ) ;
01034       imar = MRI_BYTE_PTR(im) ;
01035       if( ar->verbose ) fprintf(stderr,"  call vpSetImage(ALPHA)\n") ;
01036       vpSetImage( ar->vpc , imar , npix,npix,npix , VP_ALPHA ) ;
01037 #endif
01038    }
01039 
01040    if( ar->verbose ) fprintf(stderr,"--MREN: rendering image\n") ;
01041 
01042    if( ar->pmode == PMODE_HIGH ){
01043       if( ar->verbose ) fprintf(stderr,"  call vpRenderClassifiedVolume\n") ;
01044       fred = vpRenderClassifiedVolume(ar->vpc) ;
01045       if( fred != VP_OK ){
01046          fprintf(stderr,"**MREN: vpRenderClassifiedVolume failed: code=%d\n",(int)fred) ;
01047          mri_free(im) ; return NULL ;
01048       }
01049    } else if( ar->pmode == PMODE_MEDIUM ){
01050       if( ar->verbose ) fprintf(stderr,"  call vpRenderRawVolume\n") ;
01051       fred = vpRenderRawVolume(ar->vpc) ;
01052       if( fred != VP_OK ){
01053          fprintf(stderr,"**MREN: vpRenderRawVolume failed: code=%d\n",(int)fred) ;
01054          mri_free(im) ; return NULL ;
01055       }
01056    } else {
01057       if( ar->verbose ) fprintf(stderr,"  call vpBruteForceRender\n") ;
01058       fred = vpBruteForceRender(ar->vpc) ;
01059       if( fred != VP_OK ){
01060          fprintf(stderr,"**MREN: vpBruteForceRender failed: code=%d\n",(int)fred) ;
01061          mri_free(im) ; return NULL ;
01062       }
01063    }
01064 
01065    return im ;
01066 }