gliftUtil.cpp

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////
// 8/11/02      Aaron Lefohn    Scientific Computing and Imaging Institute
// School of Computing          University of Utah
//
//  This library is free software; you can redistribute it and/or
//  modify it under the terms of the GNU Lesser General Public
//  License as published by the Free Software Foundation; either
//  version 2.1 of the License, or (at your option) any later version.
//
//  This library is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//  Lesser General Public License for more details.
//
//  You should have received a copy of the GNU Lesser General Public
//  License along with this library; if not, write to the Free Software
//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
////////////////////////////////////////////////////////////////////////

#include <util/gliftUtil.h>
#include <texture/subTex.h>
#include <texture/texNd.h>
#include <drawable/planarQuad.h>
#include <drawable/shadedPrim.h>
#include <texCoordGen/winTexGen2D.h>

#ifdef GLIFT_NRRD
#include <nrrd.h>
#endif
#include <string>
#include <vector>

using namespace std;
using namespace gutz;

using namespace glift;

/// The OpenGL texture unit names, indexable by texUnitNumber
/// - glim.max_tex_units should be checked for number of texture
///   units available on GPU.
const GLenum g_texUnitName[NUM_TEX_UNIT_NAMES] = 
{       
   GL_TEXTURE0, 
      GL_TEXTURE1,  
      GL_TEXTURE2, 
      GL_TEXTURE3, 
      GL_TEXTURE4, 
      GL_TEXTURE5,
      GL_TEXTURE6,
      GL_TEXTURE7,
      GL_TEXTURE8,
      GL_TEXTURE9,
      GL_TEXTURE10,
      GL_TEXTURE11,
      GL_TEXTURE12,
      GL_TEXTURE13,
      GL_TEXTURE14,
      GL_TEXTURE15,
      GL_TEXTURE16, 
      GL_TEXTURE17, 
      GL_TEXTURE18, 
      GL_TEXTURE19, 
      GL_TEXTURE20, 
      GL_TEXTURE21,
      GL_TEXTURE22,
      GL_TEXTURE23,
      GL_TEXTURE24,
      GL_TEXTURE25,
      GL_TEXTURE26,
      GL_TEXTURE27,
      GL_TEXTURE28,
      GL_TEXTURE29,
      GL_TEXTURE30,
      GL_TEXTURE31,
};

// Draw pixels in 'data' into the texture, 'tex'.
// Data is sent to texture via render-to-texture or by glCopySubTexImage2D(...). 
// Dimensions of pixel rectangle are 'data.dim(1)' x 'data.dim(0)'
// Number of color channels is determined by 'data.dim(2)'
void drawPixels( const arrayw3ub& data, SingleTex* tex, GLenum drawBuffer )
{
   drawSubPixels( vec2i(0), vec2i(0), vec2i(data.dim(1), data.dim(0)), data, tex, drawBuffer );
}

// Draw a 2D tile of pixels from 'data' into 'tex'
// - 'dataOrig' and 'dimen' denote the data tile
// - 'texOrig' and 'dimen' denote where in the texture the data is written
void drawSubPixels( const vec2i& dataOrig, const vec2i& texOrig, const vec2i& dimen,
                   const arrayw3ub& data, SingleTex* tex, GLenum drawBuffer )
{
   // Set the data format
   int numChannels = data.dim(2);
   GLenum glType;
   if(   numChannels == 1 ) { glType = GL_LUMINANCE; }
   else if( numChannels == 3 ) { glType = GL_RGB; }
   else if( numChannels == 4 ) { glType = GL_RGBA; }
   else {
      cerr << "glift::drawSubPixels(...) Error:\n\tCan only draw data with 1, 3, or 4 color channels.\n";
      exit(1);
   }

   // Crop data if necessary
   ubyte* dataPtr  = (uchar*)data.data();
   bool   nukeData = false;

#ifdef GLIFT_NRRD
   if( dataOrig != 0 ) {        
      vec3i minV( 0, dataOrig.x, dataOrig.y);
      vec3i maxV( numChannels - 1, dimen.x + dataOrig.x - 1, dimen.y + dataOrig.y - 1 );

      int e = 0; 
      Nrrd* nIn = nrrdNew();
      Nrrd* nOut = nrrdNew();
      nerr( nrrdWrap(nIn, dataPtr, nrrdTypeUChar, 3, data.dim(2), data.dim(1), data.dim(0)), e );
      nerr( nrrdCrop(nOut, nIn, minV.v(), maxV.v()), e );

      dataPtr = (ubyte*)nOut->data;
      nrrdNix(nIn);
      nrrdNix(nOut);
      nukeData = true;

      if( e ) {
         cerr << "drawSubPixels(...) Error\n"
            << biffGetDone(NRRD) << endl;
         exit(1);
      }
   }
#else
   if( dataOrig != 0 ) {
      cerr << "drawSubPixels(...) Error:\n\t"
         << "'dataOrig' != 0.\n"
         << "\tGlift must be compiled with GLIFT_NRRD defined for this to work.\n";
      exit(1);
   }
#endif

   // Enable pbuffer if possible
   PBuffGlift* pbuff = tex->pbuff();
   if( pbuff ) {
      pbuff->enable();//true);
   }

   // DEBUG:
   /*vec4i curView;
   glGetIntegerv(GL_VIEWPORT, curView.v());
   cerr << "drawSingleTexQuad(...) curViewport = " << curView.x << ", " << curView.y << ", " << curView.z << ", " << curView.w << endl;
   */

   // Save/set drawBuffer
   int saveDrawBuff = 0;
   glGetIntegerv( GL_DRAW_BUFFER, &saveDrawBuff ); 
   glDrawBuffer( drawBuffer ); 

   // Draw pixels
   if( SubTex* vtex = dynamic_cast<SubTex*>(tex) ) {
      vec2i origin = vtex->origin() + texOrig;
      setColorMask( vtex->channel() ); 
      glRasterPos2i( origin.x, origin.y );
      glDrawPixels( dimen.x, dimen.y, glType, GL_UNSIGNED_BYTE, dataPtr );

      glRasterPos2i(0,0);
      setColorMask( GL_RGBA );
   }
   else {
      glRasterPos2i(texOrig.x, texOrig.y);
      glDrawPixels( dimen.x, dimen.y, glType, GL_UNSIGNED_BYTE, dataPtr );
      glRasterPos2i(0,0);
   }

   // Clean up...
   glDrawBuffer( saveDrawBuff );
   if( pbuff) {
      //pbuff->disable();
   }
   else {
      if( Tex2D* tex2D = dynamic_cast<Tex2D*>(tex) ) {
         tex2D->copyToTex2D();
      }
   }

   if( nukeData ) {free(dataPtr); }
}

// - Draw pixels in 'data' into the texture, 'tex'
// - Dimensions of pixel rectangle are 'data.dim(1)' x 'data.dim(0)'
// - Number of color channels is determined by 'data.dim(2)'
// - Draws to the currenly set surface UNLESS 'tex->pbuff()' is a SubPBuff.
//   If that is the case, the surface defined in the SubPBuff is used instead
//
// WARNINGS:
// 1) Assumes an orthographic projection is already set up in the pbuffer context
// 2) Leaves pbuffer context enabled after drawing
void drawPixels( const arrayw3f& data, SingleTex* tex )
{
   vec2i size( data.dim(1), data.dim(0) );

   // Set the data format
   int numChannels = data.dim(2);
   GLenum glType;
   if(   numChannels == 1 ) { glType = GL_LUMINANCE; }
   else if( numChannels == 3 ) { glType = GL_RGB; }
   else if( numChannels == 4 ) { glType = GL_RGBA; }
   else {
      cerr << "glift::drawSubPixels(...) Error:\n\tCan only draw data with 1, 3, or 4 color channels.\n";
      exit(1);
   }

   // Enable pbuffer if possible
   PBuffGlift* pbuff = tex->pbuff();
   if( pbuff ) {
      pbuff->enable();//true);
   }

   // Save/set drawBuffer
   int saveDrawBuff = 0;
   glGetIntegerv( GL_DRAW_BUFFER, &saveDrawBuff ); 

   // Set surface to subTex surface if possible
   if( SubPBuff* subPB = dynamic_cast<SubPBuff*>(tex->pbuff()) ) {
      glDrawBuffer( subPB->glSurface() ); 
   }

   // Draw pixels
   glRasterPos2i(0,0);
   glDrawPixels( size.x, size.y, glType, GL_FLOAT, (float*)data.data() );
   glRasterPos2i(0,0);

   // Reset draw buffer
   glDrawBuffer( saveDrawBuff );
}

void drawSingleTexQuad( SingleTex* tex, bool drawToPbuff )
{
   drawSingleTexQuad(tex, tex->dimen(), drawToPbuff );
}

// Draw a quad to the screen that is of size 'tex->dimen()'
// Will draw to 'tex->pbuff()' if 'drawToPbuff' is true
// Assumes that render mode is in appropriate state ( glOrtho2D(...)? )
void drawSingleTexQuad( SingleTex* tex, const vec2i& dimen, bool drawToPbuff )
{
   // Don't try to draw if the texture is NULL!
   if( tex == NULL ) {
      return;
   }

   // Create a texture-mapped quad with the correct texture coordinates
   PlanarQuad* quadP = NULL;

   // Generate geom with either [0, winSize] texCoords or [0,1] coords
   if( dynamic_cast<TexRect*>(tex) ) {
      // Generate geometry with texCoords[0, dimen.x] x [0, dimen.y]
      WinTexGen2D* texGen = new WinTexGen2D();
      quadP = new PlanarQuad( vec2f(), vec2f((float)dimen.x, (float)dimen.y), 0.0f, true,
         false, texGen );
      delete texGen;
   }
   else {
      quadP = new PlanarQuad( vec2f(0.0f), vec2f( (float)dimen.x, (float)dimen.y), 0.0f, true );
   }
   Shader          shader = Shader( tex );
   ShadedPrim shadedQuad( quadP, &shader );

   PBuffGlift* pbuff = tex->pbuff();
   if( drawToPbuff ) {
      if( pbuff ) {
         pbuff->enable(true);
      }

      tex->tryToBindPbuff(false);
   }

   /*   glViewport(0, 0, dimen.x, dimen.y);
   glMatrixMode(GL_PROJECTION);
   glPushMatrix();
   glLoadIdentity();
   gluOrtho2D(0.0f, (float)dimen.x, 0.0f, (float)dimen.y);
   glMatrixMode(GL_MODELVIEW);
   */
   // Set ColorMask for this subtexture
   if( SubTex* subtex = dynamic_cast<SubTex*>(tex) ) {
      if(       subtex->channel() != GL_RGBA ) {        
         glClear(GL_COLOR_BUFFER_BIT);
         setColorMask( subtex->channel() );
      }
   }

   shadedQuad.draw(); //Draw textured-quad
   setColorMask( GL_RGBA );

   /*   glMatrixMode(GL_PROJECTION);
   glPopMatrix();
   glMatrixMode(GL_MODELVIEW);
   */
   if( drawToPbuff ) {
      tex->tryToBindPbuff(true);

      if( pbuff ) {
         pbuff->disable();
      }
   }
}

// Draw a quad of 'color' to lowerLeft=origin, upperRight=origin + dimen to 
// 'pbuff' or current context (if pbuff==NULL)
void drawColoredTile( const vec2i& origin, const vec2i& dimen, const vec2i& fullQuadDim, 
                     const vec3f& color, SingleTex* tex )
{
   assert(tex);
   PBuffGlift* pbuff = tex->pbuff();

   vec2i lowerLeft  = origin;
   vec2i upperRight = origin + dimen;

   vec2f lowerLeftF( 0.0f, 0.0f );
   vec2f upperRightF( (float)fullQuadDim.x, (float)fullQuadDim.y);

   static PlanarQuad quad( lowerLeftF, upperRightF, 0.0f );

   if( pbuff ) {
      pbuff->enable();//true);
   }
   vec4i oldView;                                                       // Save old viewport
   glGetIntegerv(GL_VIEWPORT, oldView.v());

   glColor3f(color.x, color.y, color.z);        // Draw quad
   glViewport(lowerLeft.x, lowerLeft.y, dimen.x, dimen.y);
   quad.draw();

   // Restore state
   glViewport(oldView[0], oldView[1], oldView[2], oldView[3]);
   if( pbuff ) {
      //pbuff->disable();
   }
   else {
      if( Tex2D* tex2D = dynamic_cast<Tex2D*>(tex) ) {
         tex2D->copyToTex2D();
      }
   }
}

/////////////////////////////////////////////////////////////////////////////////////
// Nrrd Utilities, only valid of complile define "GLIFT_NRRD" is set
/////////////////////////////////////////////////////////////////////////////////////
#ifdef GLIFT_NRRD

/* /// Get header info from PNM
int readDimenPNM( const char* fileName, vec3i& dimen )
{
// Allocate the destination nrrd
Nrrd* nIn = nrrdNew();

int e = nrrdLoadHeader_glift( nIn, fileName );

if( !e ) {
/// Is file PGM or PPM?
GliftFormat format = nrrdTypePNM_glift(nIn);
if( format == GliftFormatPGM ) { // PGM
dimen.x = nIn->axis[0].size;
dimen.y = nIn->axis[1].size;
}
else if( format == GliftFormatPPM ) { // PPM
dimen.x = nIn->axis[1].size;
dimen.y = nIn->axis[2].size;
}
else {
cerr << "readDimenPNM(...) Error:\n"
<< "\tUnknown data format. Only PGM and PPM files are recognized.\n";
exit(1);
}
}
else {
cerr << "readDimenPNM(...) Error:\n"
<< "\tCannot read " << fileName << endl;
}

nrrdNuke(  nIn );
return e;
}*/

/// Get header info from 'fileName's header
/// - 'fileName' can be a PGM, PPM, or any Nrrd-compatible file
/// - 'dimen' will hold the size of each dimension of the dataset
/// - 'dimen.size()' is the dimensionality of the dataset
/// - Note 1: Ordering of the dimensions is Glift-like. 
///   - i.e. slowest index first, fastest index last (opposite of Nrrd convention)
///   - Example: y, x, color
/// - Note 2: Scalar data is a dimension of size 1 (counter to Nrrd convention!) 
/// - Note 3: 'isPNM' is true if file is PGM or PPM
int readDimenNrrd( const char* fileName, arrayo1i& dimen, bool& isPNM)
{       
   // Allocate the destination nrrd
   Nrrd* nIn = nrrdNew();

   int e = gliftNrrdLoadHeader( nIn, fileName );

   if( !e ) {
      /// Is file PGM or PPM?
      GliftFormat format = gliftNrrdTypePNM(nIn);

      if(       format == GliftFormatPGM ) {     // PGM
         dimen = arrayo1i( 3, (int)0);
         dimen[0] = nIn->axis[1].size;   // y
         dimen[1] = nIn->axis[0].size;   // x
         dimen[2] = 1;                                   // greyscale value
         isPNM = true;
      }
      else if( format == GliftFormatPPM ) {// PPM
         dimen = arrayo1i( 3, (int)0 );
         dimen[0] = nIn->axis[2].size;
         dimen[1] = nIn->axis[1].size;
         dimen[2] = 3;
         isPNM = true;
      }
      else {    // Any dimension
         dimen = arrayo1i( nIn->dim, (int)0 );

         // Set each dimension...in slowest to fastest order
         for(int i=0; i < dimen.size(); i++) {
            int nrrdI = i;//dimen.size() - i - 1;// Nrrd is fastest to slowest order
            dimen[i] = nIn->axis[ nrrdI ].size;
         }
         isPNM = false;
      }
   }
   else { 
      cerr << "readDimenNrrd(...) Error:\n"
         << "\tCannot read " << fileName
         << biffGetDone(NRRD) << endl;
   }

   nrrdNuke(  nIn );
   return e;
}

/// Read the PNM file, 'inFileName', sets 'width' and 'height' based on input file,
/// and returns image in 'rgbaImage'
/// - If file is PGM, all 'channels' have the replicated, scalar value
/// - If file is PPM, up to 'channels' are output.
/// - If format is RGBA and file is PPM, alpha channel is set to 255
/// - If 'makeScalar' is true, PPMs are averaged down to scalar value (replicated as necessary)
///   If false, RGB data is retained and alpha channel is padded with 255 if necessary
/// - Returns nrrd error code (0 is alles klar)
int loadPNM( GLenum channels, bool makeScalar, const char* inFileName, gutz::arrayo3ub& data )
{       
   int e = 0;//Nrrd error accumulator
   int numChannels = getNumChannels(channels);
   vec3i newDim;
   vec3i minV(0,0,0);
   vec3i maxV;

   ///Load the file
   Nrrd* nIn = nrrdNew();
   nerr( nrrdLoad(nIn, inFileName), e );

   if( e ) { //Die gracefully if there was an error
      cerr << "loadPNM(...) Error:\n"
         << biffGetDone(NRRD) << endl;
      return e;
   }

   // Figure out what kind of file it was
   GliftFormat format = gliftNrrdTypePNM(nIn);

   // Correctly format data and put result in 'nOut'
   Nrrd* nOut = nrrdNew();      //The correctly-formatted data
   if( format == GliftFormatPGM ) {
      newDim.x = nIn->axis[0].size;
      newDim.y = nIn->axis[1].size;
      newDim.z = numChannels;
      maxV = vec3i( numChannels-1, newDim.x-1, newDim.y-1 );

      // Replicate scalar value out to 'numChannels' of color
      Nrrd* nT1 = nrrdNew();
      nerr( nrrdFlip(nT1, nIn, 1), e );/// Flip the y-axis for OpenGL
      nerr( nrrdPad( nOut, nT1, minV.v(), maxV.v(), nrrdBoundaryBleed ), e );
      nrrdNuke( nT1 );
   }
   else if( format == GliftFormatPPM ) {
      newDim.x = nIn->axis[1].size;
      newDim.y = nIn->axis[2].size;
      newDim.z = numChannels;
      maxV = vec3i( numChannels-1, newDim.x-1, newDim.y-1 );

      Nrrd* nT1 = nrrdNew();
      nerr( nrrdFlip(nT1, nIn, 2), e );/// Flip the y-axis for OpenGL

      if( makeScalar ) {
         ///TODO: This should be luminance calc instead of average!!!
         int axis = 0;//Color axis
         nerr( nrrdProject( nOut, nT1, axis, nrrdMeasureMean), e );
         nerr( nrrdReshape( nOut, nOut, 3, 1, newDim.x, newDim.y), e );//Reshape as 3D

         if( numChannels > 1 ) { /// Replicate scalar in all channels
            nerr( nrrdPad( nT1, nOut, minV.v(), maxV.v(), nrrdBoundaryBleed ), e );
            nerr( nrrdCopy( nOut, nT1 ), e );
         }
      }
      else {
         if( numChannels == 1 || numChannels == 2 ) {//Data reductions
            nrrdCrop(nOut, nT1, minV.v(), maxV.v() );
         }
         else if( numChannels == 3 )  {
            nerr( nrrdCopy( nOut, nT1 ), e);
         }
         else if( numChannels == 4 ) { /// Pad alpha channel with 255
            nerr( nrrdPad( nOut, nT1, minV.v(), maxV.v(), nrrdBoundaryPad, 255 ), e );
         }
      }

      nrrdNuke( nT1 );
   }

   if( !e ) {//If alles klar...
      //dimen = newDim;

      /// The image.
      // TODO: Use a "transfer" here!!
      data.transfer( newDim.y, newDim.x, newDim.z, (uchar*)nOut->data, false );
   }
   else {
      cerr << "loadPNM(...) Error\n"
         << biffGetDone(NRRD) << endl
         << "\tError processing input. No data returned.\n";
   }

   /// Clean up
   nrrdNix( nOut );
   nrrdNix( nIn );

   return e;
}

/// Read the Nrrd file, 'inFileName'.
/// - WARNING: 'inFileName' MUST be a scalar volume or this fcn will return an error
///             and 'dimen' and 'data' will be unaltered.
/// - 'channels' is GL_RED, GL_RGB, GL_RGBA, etc and sets how many times the scalar
///    voxel value is copied.
/// - 'dimen' contains the size of each dimension of the volume.
/// - 'data'  contains the scalar volume data. Accessors are: z, y, x, channel
int loadNrrd3D( GLenum channels, const char* inFileName, arrayo4ub& data )
{
   int e = 0;//Nrrd error accumulator
   int numChannels = getNumChannels(channels);

   // Load the file
   Nrrd* nIn = nrrdNew();
   nerr( nrrdLoad(nIn, inFileName), e );

   if( e ) {                      // Die gracefully if there was an error
      cerr << "loadNrrd3D(...) Error:\n"
         << biffGetDone(NRRD) << endl;
      return e;
   }

   if( nIn->dim != 3 &&
      !(nIn->dim == 4 && nIn->axis[3].size == 1) ) { // Die gracefully if dataset is not scalar volume
         cerr << "loadNrrd3D(...) Error:\n"
            << "\t" << inFileName << " is not a scalar volume file.\n";
         return e;
      }

      if( nIn->type != nrrdTypeUChar ) { // Die gracefully if dataset is not of uchar type
         cerr << "loadNrrd3D(...) Error:\n"
            << "\t" << inFileName << " is not in unsigned byte format.\n";
         return e;
      }


      // Set the dimensions and min/max vectors
      vec4i newDim;
      newDim.x = nIn->axis[0].size;
      newDim.y = nIn->axis[1].size;
      newDim.z = nIn->axis[2].size;
      newDim.w = numChannels;

      vec4i minV( 0, 0, 0, 0);
      vec4i maxV( numChannels-1, newDim.x-1, newDim.y-1, newDim.z-1 );

      /// Correctly format data and put result in 'nOut'
      Nrrd* nOut = nrrdNew();   // The correctly-formatted data
      Nrrd* nT1  = nrrdNew();   // A temporary Nrrd

      nerr( nrrdFlip(nT1, nIn, 1), e );// Flip the y-axis for OpenGL (TODO: Is this necessary?)

      // Replicate scalar value out to 'numChannels' of color
      int newDimen = 4;
      nerr( nrrdReshape( nIn, nT1, newDimen, 1, newDim.x, newDim.y, newDim.z ), e ); // Reshape to 4D   
      nerr( nrrdPad( nOut, nIn, minV.v(), maxV.v(), nrrdBoundaryBleed ), e );

      nrrdNuke( nT1 );

      if( !e ) {//If alles klar...
         //dimen = newDim;

         // The volume. Transfer the data to the data array 
         // - Nuke whatever is in 'data' and shallow copy in the new data.
         data.transfer( newDim.z, newDim.y, newDim.x, newDim.w, (uchar*)nOut->data, false );
      }
      else {
         cerr << "loadNrrd3D(...) Error\n"
            << biffGetDone(NRRD) << endl
            << "\tError processing input. No data returned.\n";
      }

      /// Clean up
      nrrdNix( nOut );
      nrrdNuke( nIn );

      return e;
}

/// Read framebuffer to CPU memory and write to file.
void writeFBtoPPM( GLenum channels, const vec2i& origin, const vec2i& dimen, 
                  char* baseFileName, bool sepChannels )
{
   int numChannels = getNumChannels(channels);
   long numElts = dimen.x * dimen.y * numChannels;
   GLubyte* tempImage = new GLubyte[numElts];
   glReadPixels(0, 0, dimen.x, dimen.y, channels, GL_UNSIGNED_BYTE, tempImage );

   arrayw3ub data( dimen.y, dimen.x, numChannels, tempImage );
   writePNM( baseFileName, data, sepChannels );
   delete tempImage;
}

/// Write to a pgm or ppm file called 'baseFileName.pgm' or 'baseFileName.ppm'
/// - 'data' has dimens: height, width, numChannels
/// - PGM is written if 'numChannels' == 1, otherwise a PPM is written
/// - If 'numChannels' == 4, the last channel (alpha) is written to a separate pgm.
void writePNM( char* baseFileName, const arrayw3ub& data, bool sepChannels )
{
   ubyte* imageData = (ubyte*)data.data();
   vec2i  dimen( data.dim(1), data.dim(0) );
   int    numChannels = data.dim(2);
   int     numFiles = 1;

   string baseName( baseFileName );
   vector<string> fileName(numChannels);

   vector<Nrrd*> nOut(numChannels,(Nrrd*)NULL);
   Nrrd* nIn = nrrdNew();
   int e = 0;                     //Nrrd error code

   /// Initialize the output nrrds
   int i;
   for(i=0; i < numChannels; i++) {
      nOut[i] = nrrdNew();
   }

   /// Partion the data appropriately
   if( numChannels == 1 ) {
      fileName[0] = baseName  + ".pgm";

      /// Wrap the data in a nrrd
      nerr( nrrdWrap( nIn, imageData, nrrdTypeUChar, 2, dimen.x, dimen.y), e );// "2" is for "2D"
      nerr( nrrdFlip( nOut[0], nIn, 1 ), e );   //Flip y-axis
      nrrdNix( nIn );
   }
   else {
      Nrrd* nT1 = nrrdNew();
      nerr( nrrdWrap( nT1, imageData, nrrdTypeUChar, 3, numChannels, dimen.x, dimen.y ), e);// "3" is for "3D"
      nerr( nrrdFlip( nIn, nT1, 2 ), e );               //Flip y-axis
      nrrdNix( nT1 ); // WARNING: Do NOT nuke here to protect "imageData" memory!

      if( numChannels == 2 ) {
         /// Wrap the data in a nrrd

         if( sepChannels ) {
            numFiles = 2;
            fileName[0] = baseName + "R.pgm";
            fileName[1] = baseName + "G.pgm";

            int axis = 0;//Color axis
            nerr( nrrdSlice( nOut[0], nIn, axis, 0 ), e );//Red channel
            nerr( nrrdSlice( nOut[1], nIn, axis, 1 ), e );//Green channel
         }
         else {
            numFiles = 1;
            fileName[0] = baseName + ".ppm";

            /// Pad the BLUE channel with 0's and write a PPM
            vec3i minV(0,0,0);
            vec3i maxV( 2, dimen.x-1, dimen.y-1 ); 
            int padVal = 0;

            nerr( nrrdPad( nOut[0], nIn, minV.v(), maxV.v(), nrrdBoundaryPad, padVal), e );
         }
      }
      else if( numChannels == 3 ) {
         if( sepChannels ) {
            numFiles = 3;
            fileName[0] = baseName + "R.pgm";
            fileName[1] = baseName + "G.pgm";
            fileName[2] = baseName + "B.pgm";

            int axis = 0;//Color axis
            nerr( nrrdSlice( nOut[0], nIn, axis, 0 ), e );//Red channel
            nerr( nrrdSlice( nOut[1], nIn, axis, 1 ), e );//Green channel
            nerr( nrrdSlice( nOut[2], nIn, axis, 2 ), e );//Blue channel
         }
         else {
            numFiles = 1;
            fileName[0] = baseName + ".ppm";
            nOut[0] = nIn;
         }
      }
      else if( numChannels == 4 ) {
         if( sepChannels ) {
            numFiles = 4;
            fileName[0] = baseName + "R.pgm";
            fileName[1] = baseName + "G.pgm";
            fileName[2] = baseName + "B.pgm";
            fileName[3] = baseName + "A.pgm";

            int axis = 0;//Color axis
            nerr( nrrdSlice( nOut[0], nIn, axis, 0 ), e );//Red channel
            nerr( nrrdSlice( nOut[1], nIn, axis, 1 ), e );//Green channel
            nerr( nrrdSlice( nOut[2], nIn, axis, 2 ), e );//Blue channel
            nerr( nrrdSlice( nOut[3], nIn, axis, 3 ), e );//Alpha channel
         }
         else {
            numFiles = 2;
            fileName[0] = baseName + "RGB.ppm";
            fileName[1] = baseName + "A.pgm";

            /// Crop out the alpha channel to make an RGB ppm file
            vec3i minV(0,0,0);
            vec3i maxV( 2, dimen.x-1, dimen.y-1 ); 

            nerr( nrrdCrop( nOut[0], nIn, minV.v(), maxV.v() ), e );

            /// Slice out alpha channel for pgm
            int axis = 0;//Color axis
            nerr( nrrdSlice( nOut[1], nIn, axis, 3 ), e );
         }
      }
      else {
         cerr << "writePNM(...) Error\n"
            << "\tChannels cannot equal " << numChannels << endl
            << "\tNo file written.\n";
      }

      nrrdNuke( nIn );//Okay to nuke here because it is new data
   }

   /// Process Accumulated errors
   if(e) { 
      cerr << "writePNM(...) Error\n" 
         << biffGetDone(NRRD) << endl; 
   }

   /// Write the output
   for(i=0; i < numFiles; i++) {
      nrrdSave( fileName[i].c_str(), nOut[i], NULL );
   }

   /// Release nrrd memory 
   for(i=0; i < numChannels; i++) {
      if( nOut[i] ) {
         nrrdNuke( nOut[i] ); //Duke nuke 'em (everything)
      }
   }
}

/// Write to a nrrd file called 'baseFileName.nrrd'
/// - 'data' has dimens: numSlabs, height, width, numChannels
void writeNrrd3D( const char* baseFileName, const arrayw4ub& data )
{
   vec3i dimen( data.dim(2), data.dim(1), data.dim(0) );
   int   numChannels = data.dim(3);

   Nrrd* nWrap = nrrdNew();
   int e = nrrdWrap( nWrap, (ubyte*)data.data(), nrrdTypeUChar, 
      4, dimen.x, dimen.y, dimen.z, numChannels );

   Nrrd* nOut = nrrdNew();
   nerr( nrrdFlip( nOut, nWrap, 1 ), e );

   /// Process Accumulated errors
   if(e) { 
      cerr << "writeNrrd3D(...) Error\n" 
         << biffGetDone(NRRD) << endl; 
   }

   /// Write the output
   string fileName = baseFileName + string(".nrrd");
   nerr( nrrdSave( fileName.c_str(), nOut, NULL ), e );

   if( e ) {
      cerr << "gliftUtil::writeNrrd3D(...) Error : \n\t"
         << "Cannot write file \"" << fileName << "\" because\n\t"
         << biffGetDone(NRRD) << endl;
   }

   // Clean up
   nrrdNix( nWrap );
   nrrdNuke(nOut  );
}

// Take a scalar volume and pack into an RGB volume
// - sizeZ is the number of slabs in 'data' to process 
// - sizeZ MUST be a multiple of 3!!!
// - Returned data MUST be deleted to avoid memory leaks :(
arrayw4ub packScalarToRGB( const arrayw4ub& data, int sizeZ )
{
   assert( sizeZ <= data.dim(0) );
   assert( sizeZ % 3 == 0 );
   assert( data.dim(3) == 1 );

   int   e = 0;                 // Nrrd error code
   int   numChannels = 1;
   vec3i dim( data.dim(2), data.dim(1), sizeZ );
   int   zPacked = sizeZ/3;

   // Pack scalar slabs into RGB
   // 1) Wrap phiData in nIn
   Nrrd* nIn = nrrdNew();
   nerr( nrrdWrap( nIn, (ubyte*)data.data(), nrrdTypeUChar,
      4, numChannels, dim.x, dim.y, dim.z ), e );

   // 2) Crop out every 3 slices into a single, a 3-slice volume with the color and z axes swapped
   Nrrd** nSlab = new Nrrd*[zPacked];
   Nrrd*  nOut  = nrrdNew();

   vec4i minV( 0, 0, 0, 0 );
   vec4i maxV( 0, dim.x-1, dim.y-1, 2 );
   int z;
   for(z=0; z < zPacked; z++) {
      nSlab[z] = nrrdNew();
      minV[3] = z*3;
      maxV[3] = (z+1)*3 - 1;
      nerr( nrrdCrop( nOut, nIn, minV.v(), maxV.v() ), e );
      nerr( nrrdAxesSwap( nSlab[z], nOut, 0, 3 ), e );
   }

   // 3) Rebuild volume
   int joinAxis = 3;
   nerr( nrrdJoin( nOut, nSlab, zPacked, joinAxis, false ), e );

   // 4) Create arrayWrap
   arrayw4ub outData( zPacked, dim.y, dim.x, 3, (ubyte*)nOut->data );

   // Clean up
   for(z=0; z < zPacked; z++) {
      nrrdNuke(nSlab[z]);
   }
   delete nSlab;
   nrrdNix(nIn);        // The input data
   nrrdNix(nOut);       // The output data

   if(e) {
      cerr << "SolverRLS::initPhiDenseMem(...) Nrrd Error:\n\t"
         << biffGetDone(NRRD) << endl;
      exit(1);
   }

   return outData;
}

// Take an RGB-packed, scalar volume and unpack it into scalar slabs
// - Returned data MUST be deleted to avoid memory leaks :(
//       TODO: This could be fixed with a transfer operation?
// - 'unpackZ' is the maxZ value for the unpacked data
gutz::arrayw4ub unpackRGBToScalar( const gutz::arrayw4ub& data, int unpackZ )
{
   assert( unpackZ <= data.dim(0)*3 );
   assert( data.dim(3) == 3 );

   int   e = 0;                 // Nrrd error code
   int   numPackChannels = 3;
   vec3i packDim( data.dim(2), data.dim(1), data.dim(0) );
   vec3i upackDim(data.dim(2), data.dim(2), unpackZ );

   // Pack scalar slabs into RGB
   // 1) Wrap phiData in nIn
   Nrrd* nIn = nrrdNew();
   nerr( nrrdWrap( nIn, (ubyte*)data.data(), nrrdTypeUChar,
      4, numPackChannels, packDim.x, packDim.y, packDim.z ), e );

   // 2) Crop out every RGB slice into a single, a 3-slice volume with the color and z axes swapped
   Nrrd** nSlab = new Nrrd*[packDim.z];
   Nrrd*  nOut  = nrrdNew();

   vec4i minV( 0, 0, 0, 0 );
   vec4i maxV( 2, packDim.x-1, packDim.y-1, 1 );
   int z;
   for(z=0; z < packDim.z; z++) {
      nSlab[z] = nrrdNew();
      minV[3] = z;
      maxV[3] = z;
      nerr( nrrdCrop( nOut, nIn, minV.v(), maxV.v() ), e );
      nerr( nrrdAxesSwap( nSlab[z], nOut, 0, 3 ), e );
   }

   // 3) Rebuild volume
   int joinAxis = 3;
   nerr( nrrdJoin( nOut, nSlab, packDim.z, joinAxis, false ), e );
   int newZ = nOut->axis[3].size;

   // 4) Crop top slabs if necessary
   if( upackDim.z < newZ ) {
      minV = vec4i(0,0,0,0);
      maxV = vec4i(0, upackDim.x-1, upackDim.y-1, upackDim.z-1);
      nerr( nrrdCrop( nSlab[0], nOut, minV.v(), maxV.v() ), e );
      nerr( nrrdCopy( nOut, nSlab[0]), e);
   }

   // 4) Create arrayWrap
   arrayw4ub outData( upackDim.z, upackDim.y, upackDim.x, 1, (ubyte*)nOut->data );

   // Clean up
   for(z=0; z < packDim.z; z++) {
      nrrdNuke(nSlab[z]);
   }
   delete nSlab;
   nrrdNix(nIn);        // The input data
   nrrdNix(nOut);       // The output data

   if(e) {
      cerr << "SolverRLS::initPhiDenseMem(...) Nrrd Error:\n\t"
         << biffGetDone(NRRD) << endl;
      exit(1);
   }

   return outData;
}


#endif //GLIFT_NRRD

---

Send questions, comments, and bug reports to: