/*****************************************************************************
   Major portions of this software are copyrighted by the Medical College
   of Wisconsin, 1994-2000, and are released under the Gnu General Public
   License, Version 2.  See the file README.Copyright for details.
******************************************************************************/

#include "mri_render.h"

/*============================================================================
  MREN = MRI Renderer, a set of routines for volume rendering 3D bricks.
  Built on top of VolPack, by Philippe Lacroute.

  (1) Use new_MREN_renderer() to create a renderer.
  (2) Load grayscale or color bricks into it with functions
      MREN_set_graybytes(), MREN_set_rgbbytes(), or MREN_set_rgbshorts.
  (3) Load the opacity brick with MREN_set_opabrick().
  (4) Set the viewing angles with MREN_set_viewpoint().
  (5) Create an image with MREN_render().
      Loop back through (2-4) as needed to create new images.
  (6) When finished, use destroy_MREN_renderer() to clean up.
==============================================================================*/

/*--------------------------------------------------------------------------
  Initialize the short -> float color table, and
             the byte  -> float gray/opacity tables for VolPack.

  N.B.: The documentation for VolPack is incorrect -- the color values
        should be in the range 0..255, NOT 0..1.  The opacity values
        should still be in the range 0..1.
----------------------------------------------------------------------------*/

static float * MREN_colorshorts = NULL ;
static float * MREN_graytable   = NULL ;
static float * MREN_opatable    = NULL ;
static float * MREN_colorbytes  = NULL ;

void init_MREN_colortable( void )
{
   int ii , rr,gg,bb , ss ;

   if( MREN_colorshorts != NULL ) return ;  /* been here already */

   MREN_colorshorts = (float *) malloc( sizeof(float) * TOT_COLORS * 3 ) ;
   MREN_graytable   = (float *) malloc( sizeof(float) * MREN_MAX_GRAYS ) ;
   MREN_opatable    = (float *) malloc( sizeof(float) * MREN_MAX_GRAYS ) ;
   MREN_colorbytes  = (float *) malloc( sizeof(float) * MREN_MAX_GRAYS * 3 ) ;

   /*-- load linear ramp for grayscale and opacity --*/

   for( ii=0 ; ii < MREN_MAX_GRAYS ; ii++ ){
      MREN_graytable[ii] = ii ;
      MREN_opatable[ii]  = ii / 255.0 ;
   }

   /*-- load a 32 x 32 x 32 color cube [indexed by unsigned shorts] --*/

   for( rr=0 ; rr < MREN_MAX_CDIM ; rr++ ){
      for( gg=0 ; gg < MREN_MAX_CDIM ; gg++ ){
         for( bb=0 ; bb < MREN_MAX_CDIM ; bb++ ){

            ss = FIVE_TO_SHORT(rr,gg,bb) ;   /* color index */

            MREN_colorshorts[3*ss  ] = (rr * 255.0) / 31.0 ;
            MREN_colorshorts[3*ss+1] = (gg * 255.0) / 31.0 ;
            MREN_colorshorts[3*ss+2] = (bb * 255.0) / 31.0 ;
         }
      }
   }

   /*-- at the end, add the pure grays (at a higher resolution) --*/

   ss = 3 * MREN_MAX_COLORS ;
   for( ii=0 ; ii < MREN_MAX_GRAYS ; ii++ ){
      MREN_colorshorts[ss++] = ii ;
      MREN_colorshorts[ss++] = ii ;
      MREN_colorshorts[ss++] = ii ;
   }

   /*-- load a short 8 x 8 x 4 color cube [indexed by unsigned chars] --*/

   for( rr=0 ; rr < 8 ; rr++ ){
      for( gg=0 ; gg < 8 ; gg++ ){
         for( bb=0 ; bb < 4 ; bb++ ){

            ss = (rr << 5) | (gg << 2) || (bb) ;  /* color index */

            MREN_colorbytes[3*ss  ] = (rr * 255.0) / 8.0 ;
            MREN_colorbytes[3*ss+1] = (gg * 255.0) / 8.0 ;
            MREN_colorbytes[3*ss+2] = (bb * 255.0) / 4.0 ;
         }
      }
   }

   return ;
}

void destroy_MREN_colortable( void )
{
   if( MREN_colorshorts == NULL ) return ;
   free( MREN_colorshorts ); MREN_colorshorts = NULL ;
   free( MREN_graytable   ); MREN_graytable   = NULL ;
   free( MREN_opatable    ); MREN_opatable    = NULL ;
   free( MREN_colorbytes  ); MREN_colorbytes  = NULL ;
   return ;
}

static int num_renderers = 0 ;  /* global count of how many are open */

/*--------------------------------------------------------------------------
  Create a new MRI renderer.
  The return pointer is passed to all other MREN routines.
----------------------------------------------------------------------------*/

#define DEFAULT_THETA 130.0
#define DEFAULT_PHI   285.0
#define DEFAULT_PSI     0.0

void * new_MREN_renderer( void )
{
   MREN_stuff * ar ;

   ar = (MREN_stuff *) malloc( sizeof(MREN_stuff) ) ;
   ar->type = MREN_TYPE ;

   init_MREN_colortable() ;  /* in case it's not already setup */

   /*-- initialize VolPack --*/

   ar->vpc = vpCreateContext() ;

   vpSeti( ar->vpc , VP_CONCAT_MODE , VP_CONCAT_LEFT ) ;

   vpCurrentMatrix( ar->vpc , VP_MODEL ) ;
   vpIdentityMatrix( ar->vpc ) ;

   vpCurrentMatrix( ar->vpc , VP_VIEW ) ;
   vpIdentityMatrix( ar->vpc ) ;
   vpRotate( ar->vpc , VP_X_AXIS , DEFAULT_PHI   ) ;
   vpRotate( ar->vpc , VP_Y_AXIS , DEFAULT_THETA ) ;

#undef USE_CUEING
#ifdef USE_CUEING
   vpSetDepthCueing( ar->vpc , 1.0 , 0.5 ) ;
   vpEnable( ar->vpc , VP_DEPTH_CUE , 1 ) ;
#endif

   vpCurrentMatrix( ar->vpc , VP_PROJECT ) ;
   vpIdentityMatrix( ar->vpc ) ;
   vpWindow( ar->vpc , VP_PARALLEL , -0.55,0.55 , -0.55,0.55 , -0.55,0.55 ) ;

   /*-- initialize the rest of the data --*/

   ar->nx = ar->ny = ar->nz = ar->verbose = ar->newopac = ar->newvox = 0 ;
   ar->sx = ar->sy = ar->sz = 1.0 ;

   ar->theta = DEFAULT_THETA ;
   ar->phi   = DEFAULT_PHI ;
   ar->psi   = DEFAULT_PSI ;
   ar->shim  = ar->opim = NULL ;
   ar->vox   = NULL ;
   ar->pmode = PMODE_LOW ;

   ar->grayset = ar->rgbset = ar->opaset = 0 ;  /* nothing set yet */

   ar->ncmap = ar->newcmap = 0 ;
   ar->cmap  = NULL ;

   ar->min_opacity = 0.05 ;

   num_renderers ++ ;
   return (void *) ar ;
}

/*-----------------------------------------------------------------------------
  Turn depth cueing on or off -- 11 Sep 2001
-------------------------------------------------------------------------------*/

void MREN_depth_cue( void *ah , int onoff )
{
   MREN_stuff * ar = (MREN_stuff *) ah ;

   if( !ISVALID_MREN(ar) ) return ;

   vpSetDepthCueing( ar->vpc , 2.0 , 1.3863 ) ;
   vpEnable( ar->vpc , VP_DEPTH_CUE , onoff ) ;
   return ;
}

/*-----------------------------------------------------------------------------
   Get rid of an MRI renderer.
-------------------------------------------------------------------------------*/

void destroy_MREN_renderer( void * ah )
{
   MREN_stuff * ar = (MREN_stuff *) ah ;

   if( !ISVALID_MREN(ar) ) return ;

   if( ar->vox  != NULL ) free(ar->vox) ;
   if( ar->cmap != NULL ) free(ar->cmap) ;
   vpDestroyContext( ar->vpc ) ;
   free(ar) ;

   num_renderers-- ; if( num_renderers == 0 ) destroy_MREN_colortable() ;
   return ;
}

/*-----------------------------------------------------------------------------
   For debugging purposes
-------------------------------------------------------------------------------*/

void MREN_be_verbose( void * ah )
{
   MREN_stuff * ar = (MREN_stuff *) ah ;
   if( !ISVALID_MREN(ar) ) return ;
   ar->verbose = 1 ; return ;
}

void MREN_be_quiet( void * ah )
{
   MREN_stuff * ar = (MREN_stuff *) ah ;
   if( !ISVALID_MREN(ar) ) return ;
   ar->verbose = 0 ; return ;
}

/*-----------------------------------------------------------------------------
   Set the minimum voxel opacity to consider
-------------------------------------------------------------------------------*/

void MREN_set_min_opacity( void * ah , float opm )
{
   MREN_stuff * ar = (MREN_stuff *) ah ;

   if( !ISVALID_MREN(ar) ) return ;
   if( opm <= 0.0 || opm >= 1.0 ) opm = 0.05 ;
   ar->min_opacity = opm ;

   if( ar->verbose ) fprintf(stderr,"--MREN: min_opacity = %f\n",opm) ;
   return ;
}

/*-----------------------------------------------------------------------------
  Load an user defined colormap.
    ncol    = number of colors (2..65535)
    rmap[i] = red map for index i=0..ncol-1 (values 0..255), etc.
-------------------------------------------------------------------------------*/

void MREN_set_rgbmap( void * ah, int ncol, byte * rmap, byte * gmap, byte * bmap )
{
   MREN_stuff * ar = (MREN_stuff *) ah ;
   int ii ;

   if( !ISVALID_MREN(ar) ) return ;
   if( ncol < 2 || ncol > 65535 || rmap==NULL || gmap==NULL || bmap==NULL ) return ;

   if( ar->cmap != NULL ) free(ar->cmap) ;

   ar->cmap  = (float *) malloc( sizeof(float) * (3*ncol) ) ;
   ar->ncmap = ncol ;

   for( ii=0 ; ii < ncol ; ii++ ){
      ar->cmap[3*ii  ] = rmap[ii] ;
      ar->cmap[3*ii+1] = gmap[ii] ;
      ar->cmap[3*ii+2] = bmap[ii] ;
   }

   ar->newcmap = 1 ;

   if( ar->verbose ){
      fprintf(stderr,"--MREN: new colormap\n") ;
      for( ii=0 ; ii < ncol ; ii++ ){
         fprintf(stderr,"#%3d: %5.1f %5.1f %5.1f",
                 ii , ar->cmap[3*ii],ar->cmap[3*ii+1],ar->cmap[3*ii+2]) ;
         ii++ ;
         if( ii < ncol )
            fprintf(stderr,"  #%3d: %5.1f %5.1f %5.1f",
                    ii , ar->cmap[3*ii],ar->cmap[3*ii+1],ar->cmap[3*ii+2]) ;
         ii++ ;
         if( ii < ncol )
            fprintf(stderr,"  #%3d: %5.1f %5.1f %5.1f",
                    ii , ar->cmap[3*ii],ar->cmap[3*ii+1],ar->cmap[3*ii+2]) ;
         fprintf(stderr,"\n") ;
      }
   }
   return ;
}

void MREN_unset_rgbmap( void * ah )
{
   MREN_stuff * ar = (MREN_stuff *) ah ;

   if( !ISVALID_MREN(ar) || ar->cmap == NULL ) return ;
   if( ar->cmap != NULL ){ free(ar->cmap) ; ar->cmap = NULL ; }
   ar->ncmap = 0 ; ar->newcmap = 1 ;

   if( ar->verbose ) fprintf(stderr,"--MREN: delete colormap\n") ;
   return ;
}

/*-----------------------------------------------------------------------------
   Set the grayscale brick in a renderer;
   Returns -1 if an error occurs, otherwise returns 0.
   Note that this brick is NOT free-d by MREN at any point -- that is
   the user's responsibility.
-------------------------------------------------------------------------------*/

int MREN_set_graybytes( void * ah , MRI_IMAGE * grim )
{
   MREN_stuff * ar = (MREN_stuff *) ah ;
   int newvox=0 , nvox,ii ;
   byte    * gar ;
   rgbvox  * rvox ;

   /*-- sanity checks --*/

   if( !ISVALID_MREN(ar) || grim == NULL || grim->kind != MRI_byte ) return -1 ;

   if( grim->nx < 3 || grim->ny < 3 || grim->nz < 3 ){
      fprintf(stderr,"**MREN: illegal dimensions for a gray brick\n") ;
      return -1 ;
   }

   if( ar->verbose ){
      if( ar->rgbset ) fprintf(stderr,"--MREN: switching from rgb to gray brick\n") ;
      else             fprintf(stderr,"--MREN: input a new gray brick\n") ;
   }

   /*-- if have new dimensions, have to invalidate pre-existing opacity --*/

   if( ar->nx > 0 &&
       ( ar->nx != grim->nx || ar->ny != grim->ny || ar->nz != grim->nz ) ){

      ar->opim = NULL ; ar->opaset = 0 ;

      if( ar->vox != NULL ){ free(ar->vox) ; ar->vox = NULL ; }

      if( ar->verbose )
         fprintf(stderr,"--MREN: new gray brick changes volume dimensions\n"
                        "        nx:%d->%d  ny:%d->%d  nz:%d->%d\n",
                        ar->nx,grim->nx , ar->ny,grim->ny , ar->nz,grim->nz ) ;
   }

   /*-- set dimensions --*/

   ar->shim = grim ;
   ar->nx   = grim->nx ;
   ar->ny   = grim->ny ;
   ar->nz   = grim->nz ; nvox = ar->nx * ar->ny * ar->nz ;

   /*-- if need be, allocate a voxel array to hold the data --*/

   if( ar->vox == NULL ){
      ar->newvox = newvox = 1 ;
      ar->vox = (rgbvox *) malloc( sizeof(rgbvox) * nvox ) ;
      if( ar->vox == NULL ){
         fprintf(stderr,"**MREN: can't malloc workspace with new gray brick\n") ;
         return -1 ;
      } else if( ar->verbose ){
         fprintf(stderr,"--MREN: allocated new voxel array\n") ;
      }
   }

   /*-- copy grayscale data into voxel array --*/

   rvox = ar->vox ;
   gar  = MRI_BYTE_PTR(grim) ;
   for( ii=0 ; ii < nvox ; ii++ ) rvox[ii].rgb = (unsigned short) gar[ii] ;

   if( ar->rgbset ) ar->newvox = 1 ;  /* changed from color to gray */

   ar->grayset = 1 ; ar->rgbset = 0 ;
   return 0 ;
}

/*-----------------------------------------------------------------------------
   Set the opacity brick in a renderer.
   Returns -1 if an error occurs, otherwise returns 0.
   Note that this brick is NOT free-d by MREN at any point -- that is
   the user's responsibility.
-------------------------------------------------------------------------------*/

int MREN_set_opabytes( void * ah , MRI_IMAGE * opim )
{
   MREN_stuff * ar = (MREN_stuff *) ah ;
   int nvox,ii , newvox=0 ;
   byte    * gar ;
   rgbvox  * rvox ;

   /*-- sanity checks --*/

   if( !ISVALID_MREN(ar) || opim == NULL || opim->kind != MRI_byte ) return -1 ;

   if( opim->nx < 3 || opim->ny < 3 || opim->nz < 3 ){
      fprintf(stderr,"**MREN: illegal dimensions for an opacity brick\n") ;
      return -1 ;
   }

   /*-- if have new dimensions, toss old stuff that doesn't match --*/

   if( ar->nx > 0 &&
       ( ar->nx != opim->nx || ar->ny != opim->ny || ar->nz != opim->nz ) ){

      ar->shim = NULL ; ar->grayset = ar->rgbset = 0 ;

      if( ar->vox != NULL ){ free(ar->vox) ; ar->vox = NULL ; }

      if( ar->verbose )
         fprintf(stderr,"--MREN: new opacity brick changes volume dimensions\n"
                        "        nx:%d->%d  ny:%d->%d  nz:%d->%d\n",
                        ar->nx,opim->nx , ar->ny,opim->ny , ar->nz,opim->nz ) ;
   } else {
      if( ar->verbose ) fprintf(stderr,"--MREN: new opacity brick\n") ;
   }

   /*-- set dimensions --*/

   ar->opim = opim ;
   ar->nx   = opim->nx ;
   ar->ny   = opim->ny ;
   ar->nz   = opim->nz ; nvox = ar->nx * ar->ny * ar->nz ;

   /*-- if need be, allocate a voxel array to hold the data --*/

   if( ar->vox == NULL ){
      ar->newvox = newvox = 1 ;
      ar->vox = (rgbvox *) malloc( sizeof(rgbvox) * nvox ) ;
      if( ar->vox == NULL ){
         fprintf(stderr,"**MREN: can't malloc workspace with new opacity brick\n") ;
         return -1 ;
      } else if( ar->verbose ){
         fprintf(stderr,"--MREN: allocated new voxel array\n") ;
      }
   }

   /*-- load the opacity into voxel array --*/

   gar  = MRI_BYTE_PTR(ar->opim) ;
   rvox = ar->vox ;
   for( ii=0 ; ii < nvox ; ii++ ) rvox[ii].alpha = (unsigned short) gar[ii] ;

   ar->newopac = 1 ; ar->opaset = 1 ;
   return 0 ;
}

/*-----------------------------------------------------------------------------
   Convert an RGB image to a standard 8x8x4 (MREN_colorbytes)
   color-mapped image, suitable for use in MREN_set_rgbbytes.
   The output image is composed of bytes.
-------------------------------------------------------------------------------*/

MRI_IMAGE * MREN_rgb_to_colorbytes( MRI_IMAGE * rgbim )
{
   byte * rgbar , rb,gb,bb ;
   byte * shar ;
   MRI_IMAGE * shim ;
   int ii ;

   if( rgbim == NULL || rgbim->kind != MRI_rgb ) return NULL ;

   shim  = mri_new_conforming( rgbim , MRI_byte ) ;
   shar  = MRI_BYTE_PTR(shim) ;
   rgbar = MRI_RGB_PTR(rgbim) ;

   for( ii=0 ; ii < shim->nvox ; ii++ ){
      rb = rgbar[3*ii  ] >> 5 ;
      gb = rgbar[3*ii+1] >> 5 ;
      bb = rgbar[3*ii+2] >> 6 ;

      shar[ii] = (rb << 5) | (gb << 2) | bb ;  /* index into colorbytes */
   }

   return shim ;
}

/*-----------------------------------------------------------------------------
   Convert an RGB image to a standard 32x32x32+256 (MREN_colorshorts)
   color-mapped image, suitable for use in MREN_set_rgbshorts.
   The output image is composed of shorts.
-------------------------------------------------------------------------------*/

MRI_IMAGE * MREN_rgb_to_colorshorts( MRI_IMAGE * rgbim )
{
   byte * rgbar , rb,gb,bb ;
   unsigned short * shar ;
   MRI_IMAGE * shim ;
   int ii ;

   if( rgbim == NULL || rgbim->kind != MRI_rgb ) return NULL ;

   shim  = mri_new_conforming( rgbim , MRI_short ) ;
   shar  = (unsigned short *) MRI_SHORT_PTR(shim) ;
   rgbar = MRI_RGB_PTR(rgbim) ;

   for( ii=0 ; ii < shim->nvox ; ii++ ){
      rb = EIGHT_TO_FIVE(rgbar[3*ii  ]) ;
      gb = EIGHT_TO_FIVE(rgbar[3*ii+1]) ;
      bb = EIGHT_TO_FIVE(rgbar[3*ii+2]) ;

      if( rb == gb && rb == bb ){
         shar[ii] = MREN_MAX_COLORS + rgbar[3*ii] ; /* index into grayscale */
      } else {
         shar[ii] = FIVE_TO_SHORT( rb , gb, bb ) ;  /* index into color cube */
      }
   }

   return shim ;
}

/*-----------------------------------------------------------------------------
   Set the color brick in a renderer -- values in the input image are
   indices into a colormap of length <= 256 -- by default this will
   be MREN_colorbytes, but can be replaced by the user.
   Returns -1 if an error occurs, otherwise returns 0.
   Note that these bricks are NOT free-d by MREN at any point -- that is
   the user's responsibility.
-------------------------------------------------------------------------------*/

int MREN_set_rgbbytes( void * ah , MRI_IMAGE * rgbim )
{
   MREN_stuff * ar = (MREN_stuff *) ah ;
   int newvox=0 , nvox,ii ;
   byte    * gar ;
   rgbvox  * rvox ;

   /*-- sanity checks --*/

   if( !ISVALID_MREN(ar) || rgbim == NULL || rgbim->kind != MRI_byte ) return -1 ;

   if( rgbim->nx < 3 || rgbim->ny < 3 || rgbim->nz < 3 ){
      fprintf(stderr,"**MREN: illegal dimensions for a color brick\n") ; return -1 ;
   }

   /*-- if had an old gray brick, toss it (or at least its pointer) --*/

   if( ar->verbose ){
      if( ar->grayset ) fprintf(stderr,"--MREN: switching from gray to rgb brick\n") ;
      else              fprintf(stderr,"--MREN: input new rgb brick of bytes\n") ;
   }

   /*-- if have new dimensions, toss old stuff that doesn't match --*/

   if( ar->nx > 0 &&
       ( ar->nx != rgbim->nx || ar->ny != rgbim->ny || ar->nz != rgbim->nz ) ){

      ar->opim = NULL ; ar->opaset = 0 ;

      if( ar->vox != NULL ){ free(ar->vox) ; ar->vox = NULL ; }

      if( ar->verbose )
         fprintf(stderr,"--MREN: new rgb brick changes volume dimensions\n"
                        "        nx:%d->%d  ny:%d->%d  nz:%d->%d\n",
                        ar->nx,rgbim->nx , ar->ny,rgbim->ny , ar->nz,rgbim->nz ) ;
   }

   /*-- set dimensions --*/

   ar->shim = rgbim ;
   ar->nx   = rgbim->nx ;
   ar->ny   = rgbim->ny ;
   ar->nz   = rgbim->nz ; nvox = ar->nx * ar->ny * ar->nz ;

   /*-- if need be, allocate a voxel array to hold the data --*/

   if( ar->vox == NULL ){
      ar->newvox = newvox = 1 ;
      ar->vox = (rgbvox *) malloc( sizeof(rgbvox) * nvox ) ;
      if( ar->vox == NULL ){
         fprintf(stderr,"**MREN: can't malloc workspace with new color bricks\n") ;
         return -1 ;
      } else if( ar->verbose ){
         fprintf(stderr,"--MREN: allocated new voxel array\n") ;
      }
   }

   /*-- copy color data into voxel array --*/

   rvox = ar->vox ;
   gar  = MRI_BYTE_PTR(rgbim) ;
   for( ii=0 ; ii < nvox ; ii++ ) rvox[ii].rgb = (unsigned short) gar[ii] ;

   if( ar->grayset ) ar->newvox = 1 ;  /* changed from gray to color */

   ar->rgbset = 1 ; ar->grayset = 0 ;
   return 0 ;
}

/*-----------------------------------------------------------------------------
  Similar, but where the colormap indexes are unsigned shorts.
  They can refer to a user-supplied colormap, or to the standard
  32x32x32+256 colormap (MREN_colorshorts).
-------------------------------------------------------------------------------*/

int MREN_set_rgbshorts( void * ah , MRI_IMAGE * rgbim )
{
   MREN_stuff * ar = (MREN_stuff *) ah ;
   int newvox=0 , nvox,ii ;
   unsigned short * gar ;
   rgbvox * rvox ;

   /*-- sanity checks --*/

   if( !ISVALID_MREN(ar) || rgbim == NULL || rgbim->kind != MRI_short ) return -1 ;

   if( rgbim->nx < 3 || rgbim->ny < 3 || rgbim->nz < 3 ){
      fprintf(stderr,"**MREN: illegal dimensions for a color brick\n") ; return -1 ;
   }

   if( ar->verbose ){
      if( ar->grayset ) fprintf(stderr,"--MREN: switching from gray to rgb brick\n") ;
      else              fprintf(stderr,"--MREN: input new rgb brick of shorts\n") ;
   }

   /*-- if have new dimensions, toss old stuff that doesn't match --*/

   if( ar->nx > 0 &&
       ( ar->nx != rgbim->nx || ar->ny != rgbim->ny || ar->nz != rgbim->nz ) ){

      ar->opim = NULL ; ar->opaset = 0 ;

      if( ar->vox != NULL ){ free(ar->vox) ; ar->vox = NULL ; }

      if( ar->verbose )
         fprintf(stderr,"--MREN: new rgb brick changes volume dimensions\n"
                        "        nx:%d->%d  ny:%d->%d  nz:%d->%d\n",
                        ar->nx,rgbim->nx , ar->ny,rgbim->ny , ar->nz,rgbim->nz ) ;
   }

   /*-- set dimensions --*/

   ar->shim = rgbim ;
   ar->nx   = rgbim->nx ;
   ar->ny   = rgbim->ny ;
   ar->nz   = rgbim->nz ; nvox = ar->nx * ar->ny * ar->nz ;

   /*-- if need be, allocate a voxel array to hold the data --*/

   if( ar->vox == NULL ){
      ar->newvox = newvox = 1 ;
      ar->vox = (rgbvox *) malloc( sizeof(rgbvox) * nvox ) ;
      if( ar->vox == NULL ){
         fprintf(stderr,"**MREN: can't malloc workspace with new color bricks\n") ;
         return -1 ;
      } else if( ar->verbose ){
         fprintf(stderr,"--MREN: allocated new voxel array\n") ;
      }
   }

   /*-- copy color data into voxel array --*/

   rvox = ar->vox ;
   gar  = (unsigned short *) MRI_SHORT_PTR(rgbim) ;
   for( ii=0 ; ii < nvox ; ii++ ) rvox[ii].rgb = gar[ii] ;

   if( ar->grayset ) ar->newvox = 1 ;  /* changed from gray to color */

   ar->rgbset = 2 ; ar->grayset = 0 ;
   return 0 ;
}

/*-------------------------------------------------------------------------------
   Set the viewpoint of the user (in polar angles)
---------------------------------------------------------------------------------*/

void MREN_set_viewpoint( void * ah , float theta , float phi , float psi )
{
   MREN_stuff * ar = (MREN_stuff *) ah ;

   if( !ISVALID_MREN(ar) ) return ;

   ar->theta = theta ; ar->phi = phi ; ar->psi = psi ;

   vpCurrentMatrix( ar->vpc , VP_VIEW ) ;
   vpIdentityMatrix( ar->vpc ) ;
   vpRotate( ar->vpc , VP_Z_AXIS , psi   ) ;  /* roll  */
   vpRotate( ar->vpc , VP_X_AXIS , phi   ) ;  /* pitch */
   vpRotate( ar->vpc , VP_Y_AXIS , theta ) ;  /* yaw   */

   if( ar->verbose ){
      vpMatrix4 vpm ;

      fprintf(stderr,"--MREN: set theta=%f  phi=%f  psi=%f\n",theta,phi,psi) ;

      vpGetMatrix( ar->vpc , VP_VIEW , vpm ) ;
      fprintf(stderr,"--matrix: %8.5f %8.5f %8.5f %8.5f\n"
                     "          %8.5f %8.5f %8.5f %8.5f\n"
                     "          %8.5f %8.5f %8.5f %8.5f\n"
                     "          %8.5f %8.5f %8.5f %8.5f\n" ,
              vpm[0][0] , vpm[0][1] , vpm[0][2] , vpm[0][3] ,
              vpm[1][0] , vpm[1][1] , vpm[1][2] , vpm[1][3] ,
              vpm[2][0] , vpm[2][1] , vpm[2][2] , vpm[2][3] ,
              vpm[3][0] , vpm[3][1] , vpm[3][2] , vpm[3][3]  ) ;
   }

   return ;
}

/*-------------------------------------------------------------------------------
   Set the precalculation mode of the renderer
---------------------------------------------------------------------------------*/

void MREN_set_precalculation( void * ah , int mode )
{
   MREN_stuff * ar = (MREN_stuff *) ah ;

   if( !ISVALID_MREN(ar) || mode < PMODE_LOW || mode > PMODE_HIGH ) return ;

   if( ar->pmode != mode ){ ar->pmode = mode ; ar->newopac = 1 ; }
   return ;
}

/*-------------------------------------------------------------------------------
   Set the scale factors for each axis (default = 1).  The inputs should
   be the size of the brick along each axes (e.g., sx = nx * dx).  This is
   needed because VolPack assumes the input data is in a unit cube, but
   our data may not be so uniform.

   N.B.: This is disabled, since VolPack doesn't seem to work in this regards.
---------------------------------------------------------------------------------*/

void MREN_set_size( void * ah , float sx , float sy , float sz )
{
#if 0
   MREN_stuff * ar = (MREN_stuff *) ah ;
   float mmm ;

   if( !ISVALID_MREN(ar) ) return ;

   sx = fabs(sx) ; if( sx == 0.0 ) sx = 1.0 ;  /* don't allow non-positive sizes */
   sy = fabs(sy) ; if( sy == 0.0 ) sy = 1.0 ;
   sz = fabs(sz) ; if( sz == 0.0 ) sz = 1.0 ;

   mmm = sx ;
   if( mmm < sy ) mmm = sy ;
   if( mmm < sz ) mmm = sz ;  /* mmm = maximum size */

   ar->sx = sx / mmm ;        /* scale factors are <= 1.0 */
   ar->sy = sy / mmm ;
   ar->sz = sz / mmm ;

   vpCurrentMatrix( ar->vpc , VP_MODEL ) ;  /* scale model to world */
   vpIdentityMatrix( ar->vpc ) ;
   vpScale( ar->vpc , sx , sy , sz ) ;

   if( ar->verbose )
      fprintf(stderr,"--MREN: set scale factors = %f %f %f\n",ar->sx,ar->sy,ar->sz) ;
#endif

   return ;
}

/*-------------------------------------------------------------------------------
   Find out if a renderer needs data
---------------------------------------------------------------------------------*/

int MREN_needs_data( void * ah )
{
   MREN_stuff * ar = (MREN_stuff *) ah ;

   return (ar->vox == NULL) ;
}

/*-------------------------------------------------------------------------------
   Actually render an image.  Returns NULL if an error occurs.
   Input npix = number of pixels on a side of the image (always will be square).
   If rendering a grayscale brick, returns an image of kind MRI_byte.
   If rendering RGB bricks, returns an image of kind MRI_rgb.
---------------------------------------------------------------------------------*/

MRI_IMAGE * MREN_render( void * ah , int npix )
{
   MREN_stuff * ar = (MREN_stuff *) ah ;
   int isgray , isrgb ;
   MRI_IMAGE * im ;
   byte * imar ;
   vpResult fred ;

   /*-- sanity checks --*/

   if( !ISVALID_MREN(ar) ) return NULL ;

   if( npix < 16 ){
      fprintf(stderr,"**MREN: attempt to render with less than 16 pixels!\n") ;
      return NULL ;
   }

   isgray = (ar->grayset > 0) ;
   isrgb  = (ar->rgbset  > 0) ;

   if( isgray && isrgb ){
      fprintf(stderr,"**MREN: attempt to render gray and color simultaneously?\n");
      return NULL ;
   }

   if( (!isgray && !isrgb) || ar->vox == NULL ){
      fprintf(stderr,"**MREN: attempt to render without data being loaded!\n") ;
      return NULL ;
   }

   if( ar->opaset == 0 ){
      fprintf(stderr,"**MREN: attempt to render without opacity being loaded!\n") ;
      return NULL ;
   }

   if( ar->nx < 3 || ar->ny < 3 || ar->nz < 3 ){
      fprintf(stderr,"**MREN: attempt to render without initialization!\n") ;
      return NULL ;
   }

   if( ar->verbose ) fprintf(stderr,"--MREN: setup for rendering\n") ;

   /*-- if have new voxel array, must tell VolPack all about it --*/

   if( ar->newvox || (isrgb && ar->newcmap) ){
      int nvox = ar->nx * ar->ny * ar->nz ;
      int oct_range ;
      rgbvox vv , *rv = &vv ;

      /* 3D dimensions */

      if( ar->verbose ) fprintf(stderr,"  call vpSetVolumeSize\n") ;

      fred = vpSetVolumeSize( ar->vpc , ar->nx , ar->ny , ar->nz ) ;
      if( fred != VP_OK ){
         fprintf(stderr,"**MREN: vpSetVolumeSize failed: code=%d\n",(int)fred) ;
         return NULL ;
      }

      /* each voxel has 2 data fields; 1 for shading and 1 for opacity */

      if( ar->verbose ) fprintf(stderr,"  call vpSetVoxelSize\n") ;

      fred = vpSetVoxelSize( ar->vpc , sizeof(rgbvox) , 2 , 1 , 1 ) ;
      if( fred != VP_OK ){
         fprintf(stderr,"**MREN: vpSetVoxelSize failed: code=%d\n",(int)fred) ;
         return NULL ;
      }

      /* voxel field 1 (alpha) is an index into MREN_opatable */

      if( ar->verbose ) fprintf(stderr,"  call vpSetVoxelField(1)\n") ;

      fred = vpSetVoxelField( ar->vpc, 1, sizeof(short),
                              vpFieldOffset(rv,alpha), MREN_MAX_GRAYS-1 );
      if( fred != VP_OK ){
         fprintf(stderr,"**MREN: vpSetVoxelField(1) failed: code=%d\n",(int)fred) ;
         return NULL ;
      }

      /* tell VolPack where the voxels are */

      if( ar->verbose ) fprintf(stderr,"  call vpSetRawVoxels\n") ;

      fred = vpSetRawVoxels( ar->vpc , ar->vox ,
                             sizeof(rgbvox)*nvox ,
                             sizeof(rgbvox) ,
                             sizeof(rgbvox)*(ar->nx) ,
                             sizeof(rgbvox)*(ar->nx * ar->ny) ) ;
      if( fred != VP_OK ){
         fprintf(stderr,"**MREN: vpSetRawVoxels failed: code=%d\n",(int)fred) ;
         return NULL ;
      }

      /*--- gray scale input ---*/

      if( isgray ){

         /* voxel field 0 (rgb) is an index into MREN_graytable */

         if( ar->verbose ) fprintf(stderr,"  call vpSetVoxelField(grays)\n") ;

         fred = vpSetVoxelField( ar->vpc, 0, sizeof(short),
                                 vpFieldOffset(rv,rgb), MREN_MAX_GRAYS-1 );
         if( fred != VP_OK ){
            fprintf(stderr,"**MREN: vpSetVoxelField(0) failed: code=%d\n",(int)fred) ;
            return NULL ;
         }

         /* setup MREN_graytable to hold the colormap */

         if( ar->verbose ) fprintf(stderr,"  call vpSetLookupShader(graytable)\n") ;

         fred = vpSetLookupShader( ar->vpc , 1 , 1 , 0 ,
                                   MREN_graytable , sizeof(float)*MREN_MAX_GRAYS ,
                                   0 , NULL , 0 ) ;
         if( fred != VP_OK ){
            fprintf(stderr,"**MREN: vpSetLookupShader failed: code=%d\n",(int)fred) ;
            return NULL ;
         }

      /*--- color input ---*/

      } else if( isrgb ){

         /* There are 3 possible cases for the colormap:
              a) user supplied colormap
              b) standard 16 bit index colormap MREN_colorshorts
              c) standard 8 bit index colormap MREN_colorbytes    */

         if( ar->cmap != NULL && ar->ncmap > 1 ){   /* user supplied colormap */

            if( ar->verbose ) fprintf(stderr,"  call vpSetVoxelField(cmap)\n") ;

            fred = vpSetVoxelField( ar->vpc, 0, sizeof(short),
                                    vpFieldOffset(rv,rgb), ar->ncmap-1 );

            if( ar->verbose ) fprintf(stderr,"  call vpSetLookupShader(cmap)\n") ;

            fred = vpSetLookupShader( ar->vpc , 3 , 1 , 0 ,
                                      ar->cmap , sizeof(float)*3*ar->ncmap ,
                                      0 , NULL , 0 ) ;

         } else if( ar->rgbset == 2 ){          /* MREN_colorshorts */

            if( ar->verbose ) fprintf(stderr,"  call vpSetVoxelField(rgb shorts)\n") ;

            fred = vpSetVoxelField( ar->vpc, 0, sizeof(short),
                                    vpFieldOffset(rv,rgb), TOT_COLORS-1 );

            if( ar->verbose ) fprintf(stderr,"  call vpSetLookupShader(colorshorts)\n") ;

            fred = vpSetLookupShader( ar->vpc , 3 , 1 , 0 ,
                                      MREN_colorshorts , sizeof(float)*TOT_COLORS*3 ,
                                      0 , NULL , 0 ) ;

         } else {                               /* MREN_colorbytes */

            if( ar->verbose ) fprintf(stderr,"  call vpSetVoxelField(rgb bytes)\n") ;

            fred = vpSetVoxelField( ar->vpc, 0, sizeof(short),
                                    vpFieldOffset(rv,rgb), MREN_MAX_GRAYS-1 );

            if( ar->verbose ) fprintf(stderr,"  call vpSetLookupShader(colorbytes)\n") ;

            fred = vpSetLookupShader( ar->vpc , 3 , 1 , 0 ,
                                      MREN_colorbytes , sizeof(float)*MREN_MAX_GRAYS*3 ,
                                      0 , NULL , 0 ) ;
         }

         if( fred != VP_OK ){
            fprintf(stderr,"**MREN: vpSetLookupShader failed: code=%d\n",(int)fred) ;
            return NULL ;
         }
      }

      /*-- in all cases, voxel field 1 (alpha) is an index into MREN_opatable --*/

      if( ar->verbose ) fprintf(stderr,"  call vpSetClassifierTable\n") ;

      fred = vpSetClassifierTable( ar->vpc, 0, 1, MREN_opatable, sizeof(float)*MREN_MAX_GRAYS ) ;
      if( fred != VP_OK ){
         fprintf(stderr,"**MREN: vpSetClassifierTable failed: code=%d\n",(int)fred) ;
         return NULL ;
      }

      /* threshold for octree bins: 12 = 5% of possible opacity range */

      if( ar->verbose ) fprintf(stderr,"  call vpMinMaxOctreeThreshold\n") ;

      fred = vpMinMaxOctreeThreshold( ar->vpc , 0 , 12 ) ;
      if( fred != VP_OK ){
         fprintf(stderr,"**MREN: vpMinMaxOctreeThreshold failed: code=%d\n",(int)fred) ;
         return NULL ;
      }

      ar->newopac = 1 ;  /* for the precalculations below */
      ar->newvox  = 0 ;
      ar->newcmap = 0 ;
   }

   /*-- if have new data in the voxel array, must do precalculations --*/

   vpSetd( ar->vpc , VP_MAX_RAY_OPACITY   , 0.95 ) ;
   vpSetd( ar->vpc , VP_MIN_VOXEL_OPACITY , ar->min_opacity ) ;

   if( ar->newopac ){

      (void) vpDestroyMinMaxOctree( ar->vpc ) ;      /* toss previous work */
      (void) vpDestroyClassifiedVolume( ar->vpc ) ;

      if( ar->pmode == PMODE_MEDIUM ){

         if( ar->verbose ) fprintf(stderr,"--MREN: computing octree\n") ;

         /* make octree, down to 4x4x4 voxel bins */

         if( ar->verbose ) fprintf(stderr,"  call vpCreateMinMaxOctree\n") ;

         fred = vpCreateMinMaxOctree( ar->vpc , 0 , 4 ) ;
         if( fred != VP_OK ){
            fprintf(stderr,"**MREN: vpCreateMinMaxOctree failed: code=%d\n",(int)fred) ;
            return NULL ;
         }

      } else if( ar->pmode == PMODE_HIGH ){

         if( ar->verbose ) fprintf(stderr,"--MREN: computing classified volume\n") ;

         /* classify volume (slower than octree, but may do faster rendering) */

         if( ar->verbose ) fprintf(stderr,"  call vpClassifyVolume\n") ;

         fred = vpClassifyVolume( ar->vpc ) ;
         if( fred != VP_OK ){
            fprintf(stderr,"**MREN: vpClassifyVolume failed: code=%d\n",(int)fred) ;
            return NULL ;
         }
      }

      ar->newopac = 0 ;
   }

   /*-- create the output image --*/

#undef GET_ALPHA  /* for debugging: compute the opacity image */

   if( isgray ){
      im   = mri_new( npix , npix , MRI_byte ) ;
      imar = MRI_BYTE_PTR(im) ;
#ifndef GET_ALPHA
      if( ar->verbose ) fprintf(stderr,"  call vpSetImage(LUMINANCE)\n") ;
      vpSetImage( ar->vpc , imar , npix,npix,npix , VP_LUMINANCE ) ;
#else
      if( ar->verbose ) fprintf(stderr,"  call vpSetImage(ALPHA)\n") ;
      vpSetImage( ar->vpc , imar , npix,npix,npix , VP_ALPHA ) ;
#endif
   } else if( isrgb ){
#ifndef GET_ALPHA
      im   = mri_new( npix , npix , MRI_rgb ) ;
      imar = MRI_RGB_PTR(im) ;
      if( ar->verbose ) fprintf(stderr,"  call vpSetImage(RGB)\n") ;
      vpSetImage( ar->vpc , imar , npix,npix,3*npix , VP_RGB ) ;
#else
      im   = mri_new( npix , npix , MRI_byte ) ;
      imar = MRI_BYTE_PTR(im) ;
      if( ar->verbose ) fprintf(stderr,"  call vpSetImage(ALPHA)\n") ;
      vpSetImage( ar->vpc , imar , npix,npix,npix , VP_ALPHA ) ;
#endif
   }

   if( ar->verbose ) fprintf(stderr,"--MREN: rendering image\n") ;

   if( ar->pmode == PMODE_HIGH ){
      if( ar->verbose ) fprintf(stderr,"  call vpRenderClassifiedVolume\n") ;
      fred = vpRenderClassifiedVolume(ar->vpc) ;
      if( fred != VP_OK ){
         fprintf(stderr,"**MREN: vpRenderClassifiedVolume failed: code=%d\n",(int)fred) ;
         mri_free(im) ; return NULL ;
      }
   } else if( ar->pmode == PMODE_MEDIUM ){
      if( ar->verbose ) fprintf(stderr,"  call vpRenderRawVolume\n") ;
      fred = vpRenderRawVolume(ar->vpc) ;
      if( fred != VP_OK ){
         fprintf(stderr,"**MREN: vpRenderRawVolume failed: code=%d\n",(int)fred) ;
         mri_free(im) ; return NULL ;
      }
   } else {
      if( ar->verbose ) fprintf(stderr,"  call vpBruteForceRender\n") ;
      fred = vpBruteForceRender(ar->vpc) ;
      if( fred != VP_OK ){
         fprintf(stderr,"**MREN: vpBruteForceRender failed: code=%d\n",(int)fred) ;
         mri_free(im) ; return NULL ;
      }
   }

   return im ;
}


syntax highlighted by Code2HTML, v. 0.9.1