gliftNrrd.cpp

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////
// 9/6/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
////////////////////////////////////////////////////////////////////////


#ifdef GLIFT_NRRD

#include <util/gliftNrrd.h>
#include <iostream>
#include <stdio.h>

using namespace std;
using namespace gutz;

using namespace glift;

/// Convert a nrrdType to an OpenGL type
GLenum typeN2GL( int nrrdType )
{
   GLenum glType;

   switch( nrrdType )
   {
   case nrrdTypeChar:   glType = GL_BYTE;                       break;  /*   signed 1-byte integer */
   case nrrdTypeUChar: glType = GL_UNSIGNED_BYTE;       break;  /* unsigned 1-byte integer */
   case nrrdTypeShort: glType = GL_SHORT;                       break;  /*   signed 2-byte integer */
   case nrrdTypeUShort:glType = GL_UNSIGNED_SHORT;      break;  /* unsigned 2-byte integer */
   case nrrdTypeInt:    glType = GL_INT;                        break;  /*   signed 4-byte integer */
   case nrrdTypeUInt:   glType = GL_UNSIGNED_INT;       break;  /* unsigned 4-byte integer */
   case nrrdTypeFloat: glType = GL_FLOAT;                       break;  /*          4-byte floating point */
   default:
      cerr << "typeN2GL(" << nrrdType << ") Error: OpenGL does not support that type\n";
      exit(1);
   }

   return glType;
}

/// Convert an OpenGL type to a nrrdType
GLenum typeGL2N( GLenum glType)
{
   int nrrdType;

   switch( glType )
   {
   case GL_BYTE:                        nrrdType = nrrdTypeChar;        break;  /*   signed 1-byte integer */
   case GL_UNSIGNED_BYTE:       nrrdType = nrrdTypeUChar;       break;  /* unsigned 1-byte integer */
   case GL_SHORT:                       nrrdType = nrrdTypeShort;       break;  /*   signed 2-byte integer */
   case GL_UNSIGNED_SHORT:      nrrdType = nrrdTypeUShort ;     break;  /* unsigned 2-byte integer */
   case GL_INT:                 nrrdType = nrrdTypeInt;         break;  /*   signed 4-byte integer */
   case GL_UNSIGNED_INT:        nrrdType = nrrdTypeUInt;        break;  /* unsigned 4-byte integer */
   case GL_FLOAT:                       nrrdType = nrrdTypeFloat;       break;  /*          4-byte floating point */
   default:
      cerr << "typeGL2N(" << glType << ") Error: Nrrd does not support that type\n";
      exit(1);
   }

   return nrrdType;
}

// Return GliftFormatPGM, GliftFormatPPM, or GliftFormatUnknown
GliftFormat gliftNrrdTypePNM( Nrrd* nIn )
{
   int format = nrrdFitsInFormat(nIn, nrrdEncodingRaw, nrrdFormatPNM, AIR_FALSE);
   GliftFormat type = GliftFormatUnknown; 

   if( format == 2 ) {
      type = GliftFormatPGM;
   }
   else if( format == 3 ) {
      type = GliftFormatPPM;
   }

   return type;
}

// Read just the Nrrd header of 'fileName'
int gliftNrrdLoadHeader( Nrrd* nIn, const char* fileName )
{
   char me[] = "nrrdLoadHeader_glift", err[AIR_STRLEN_MED];
   airArray *mop;

   mop = airMopNew();
   FILE *inFile = fopen( fileName, "rb" );
   if( !inFile ) {
      fprintf(stderr, "nrrdLoadHeader_glift(...) Error: Can't open input file '%s'.\nProgram aborted.\n\n", fileName);
      exit(1);
   }
   airMopAdd(mop, inFile, (airMopper)airFclose, airMopAlways);

   // Allocate a nrrdStruct
   NrrdIO* nIO = nrrdIONew();
   airMopAdd(mop, nIO, (airMopper)nrrdIONix, airMopAlways);
   nIO->skipData = true;        // Only read the header
   nIO->dir = airStrdup("This is a hack that should not need to exist!");

   // Read the header
   if( nrrdRead(nIn, inFile, nIO) ) {
      sprintf(err, "%s: Error in reading header of file \"%s\"", me, fileName);
      biffAdd(NRRD, err); airMopError(mop); return 1;
   }

   // Cleanup
   airMopOkay(mop);
   //airArrayNix(mop);
   return 0;
}

// Flip the y-axis for data in 'data'
void gliftNrrdFlipYub( arrayo2ub& data )
{       
   int e = 0;
   Nrrd* nIn = nrrdNew();
   Nrrd* nOut= nrrdNew();

   nerr( nrrdWrap( nIn, (void*)data.data(), nrrdTypeUChar, 2, data.dim(1), data.dim(0) ), e );
   nerr( nrrdFlip( nOut, nIn, 1 ), e );

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

   data.transfer( data.dim(0), data.dim(1), (ubyte*)nOut->data, false);

   nrrdNix(nIn);
   nrrdNix(nOut);
}


// Tile 3D scalar data only in X and Y
void gliftNrrdTileDataXYub( const arrayw3ub& srcData3, const vec2i& tileSize, arrayo5ub& tiledData )
{
   vec3i srcSize( srcData3.dim(2), srcData3.dim(1), srcData3.dim(0) );
   vec2i numTiles = vec2i(srcSize) / tileSize;
   assert( srcSize.x % tileSize.x == 0 );
   assert( srcSize.y % tileSize.y == 0 );

   int e = 0;
   Nrrd* nIn  = nrrdNew();
   Nrrd* nOut = nrrdNew();

   // Create nrrd of: sx, tx, sy, ty, z
   nerr( nrrdWrap( nIn, (void*)srcData3.data(), nrrdTypeUChar, 3, srcSize.x, srcSize.y, srcSize.z ), e );
   nerr( nrrdReshape( nIn, nIn, 5, tileSize.x, numTiles.x, tileSize.y, numTiles.y, srcSize.z ), e );

   // Change to: sx, sy, tx, ty, z
   int axes[5] = {0, 2, 1, 3, 4};
   nerr( nrrdAxesPermute( nOut, nIn, axes ), e );

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

   // Store tiled monster in arrayo5ub
   tiledData.transfer( srcSize.z, numTiles.y, numTiles.x, tileSize.y, tileSize.x, (ubyte*)nOut->data, false); 

   nrrdNix(nIn);
   nrrdNix(nOut);
}

// Untile XY-tiled (5D), 3D data into a standard 3D scalar format
// - Opposite of "gliftNrrdTileDataXYub"
void gliftNrrdUntileXYub( const arrayw5ub& tiledData, arrayo3ub& untiledData )
{
   vec2i numTiles( tiledData.dim(2), tiledData.dim(1) );
   vec2i tileSize( tiledData.dim(4), tiledData.dim(3) );
   vec3i untiledSize( numTiles.x * tileSize.x, numTiles.y * tileSize.y, tiledData.dim(0) );

   int e = 0;
   Nrrd* nIn = nrrdNew();
   Nrrd* nOut= nrrdNew();

   // Create Nrrd of: sx, sy, tx, ty, z
   nerr( nrrdWrap( nIn, (void*)tiledData.data(), nrrdTypeUChar, 5, 
      tiledData.dim(4), tiledData.dim(3), tiledData.dim(2), tiledData.dim(1), tiledData.dim(0) ), e );

   // Change to: sx, tx, sy, ty, z
   int axes[5] = {0, 2, 1, 3, 4};
   nerr( nrrdAxesPermute( nOut, nIn, axes ), e );

   // Reshape to: x, y, z where x = sx*tx and y = sy*ty
   nerr( nrrdReshape( nOut, nOut, 3, untiledSize.x, untiledSize.y, untiledSize.z ), e );

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

   // Transfer ownership of untiled data to an arrayo3ub
   untiledData.transfer( untiledSize.z, untiledSize.y, untiledSize.x, (ubyte*)nOut->data, false); 

   nrrdNix(nIn);
   nrrdNix(nOut);
}

// Helper function for "resampleVolBoxXY"
// - Private to this file
// - Set 'rescale' to 1.0 if you do NOT want the axis resampled
NrrdResampleInfo* f_buildResampleInfo( const vec3d& rescale, Nrrd* nIn )
{
   // Setup the sampleInfo struct
   NrrdResampleInfo *info = nrrdResampleInfoNew();

   // 1) Set the number of samples on each axis
   // 2) Set the min/max for each Axis
   for(int d=0; d < 3; d++) {

      int numSamples = nIn->axis[d].size  * rescale[d];
      info->samples[d] = numSamples > 0 ? numSamples : 1;
      nrrdAxisMinMaxSet(nIn, d, nrrdDefCenter);

      // Resample with box filter unless numSamples==1 or rescale==1
      info->kernel[d] = ((info->samples[d] > 1) && (fabs(rescale[d] - 1.0) > 1e-6)) ? nrrdKernelBox : NULL;
      info->min[d] = nIn->axis[d].min;
      info->max[d] = nIn->axis[d].max;
   }

   // Output a float Nrrd
   //info->type = nrrdTypeFloat;

   return info;
}

// Downsample 'scalarVol' to 'resampledVol' with a box filter
//              unu resample -k box -s x(rescale.x) x(rescale.y) x(rescale.z)
//              - Set 'rescale' to 1.0 if you do NOT want the axis resampled
void gliftNrrdResampleVolBoxf( const vec3d& rescale, const arrayw3f& scalarVol, arrayo3f& resampledVol)
{
   int e = 0;           // Nrrd error accumulator

   // Turn scalarVol into a Nrrd
   Nrrd* nIn  = nrrdNew();
   nerr( nrrdWrap( nIn, (float*)scalarVol.data(), nrrdTypeFloat, 3, scalarVol.dim(2), scalarVol.dim(1), scalarVol.dim(0) ), e );

   // Build resample info structure
   // - Use box kernel on X and Y axes
   NrrdResampleInfo* info = f_buildResampleInfo( rescale, nIn );

   // Do the resampling
   Nrrd* nOut = nrrdNew();
   nIn->axis[0].center = nrrdCenterCell;
   nIn->axis[1].center = nrrdCenterCell;
   nIn->axis[2].center = nrrdCenterCell;
   nerr( nrrdSpatialResample( nOut, nIn, info ), e );

   // Check for Nrrd Errors
   if( e ) {
      char* err = biffGetDone(NRRD);
      cerr << (err ? err : "") << endl
         << "\agliftNrrdResampleVolBox() Error:\n";
      exit(1);
   }

   // Nuke contents of "resampled" and transfer downsampled data to "resampled"
   resampledVol.transfer( nOut->axis[2].size, nOut->axis[1].size, nOut->axis[0].size, (float*)nOut->data, false);

   // Clean up
   nrrdNix(  nIn );
   nrrdNix( nOut );
   nrrdResampleInfoNix( info );
}

// Downsample 'scalarVol' to 'resampledVol' with a box filter
//              unu resample -k box -s x(rescale.x) x(rescale.y) x(rescale.z)
//              - Set 'rescale' to 1.0 if you do NOT want the axis resampled
void gliftNrrdResampleVolBoxub( const vec3d& rescale, const arrayw3ub& scalarVol, arrayo3ub& resampledVol)
{
   int e = 0;           // Nrrd error accumulator

   // Turn scalarVol into a Nrrd
   Nrrd* nIn  = nrrdNew();
   nerr( nrrdWrap( nIn, (uchar*)scalarVol.data(), nrrdTypeUChar, 3, scalarVol.dim(2), scalarVol.dim(1), scalarVol.dim(0) ), e );

   // Build resample info structure
   // - Use box kernel on X and Y axes
   NrrdResampleInfo* info = f_buildResampleInfo( rescale, nIn );

   // Do the resampling
   Nrrd* nOut = nrrdNew();
   nIn->axis[0].center = nrrdCenterCell;
   nIn->axis[1].center = nrrdCenterCell;
   nIn->axis[2].center = nrrdCenterCell;
   nerr( nrrdSpatialResample( nOut, nIn, info ), e );

   // Check for Nrrd Errors
   if( e ) {
      char* err = biffGetDone(NRRD);
      cerr << (err ? err : "") << endl
         << "\agliftNrrdResampleVolBox() Error:\n";
      exit(1);
   }

   // Nuke contents of "resampled" and transfer downsampled data to "resampled"
   resampledVol.transfer( nOut->axis[2].size, nOut->axis[1].size, nOut->axis[0].size, (uchar*)nOut->data, false);

   // Clean up
   nrrdNix(  nIn );
   nrrdNix( nOut );
   nrrdResampleInfoNix( info );
}


#endif

---

Send questions, comments, and bug reports to: