nrroGinsu.cpp

Go to the documentation of this file.

//------------------------------------------------------------------------
//
//   Joe Kniss
//     6-27-03
//                   ________    ____   ___ 
//                  |        \  /    | /  /
//                  +---+     \/     |/  /
//                  +--+|  |\    /|     < 
//                  |  ||  | \  / |  |\  \ 
//                  |      |  \/  |  | \  \ 
//                   \_____|      |__|  \__\
//                       Copyright  2003 
//                      Joe Michael Kniss
//                   <<< jmk@cs.utah.edu >>>
//               "All Your Base are Belong to Us"
//-------------------------------------------------------------------------

//nrroGinsu.cpp

// slice and dice

#include "nrro.h"
#include "nrroKernel.h"
#include "nrroDispatch.h"
#include <sstream>
#include <typeinfo>

using namespace gutz;
using namespace std;

/////////////////////////////////////////////////////////////////////////
// slice
/////////////////////////////////////////////////////////////////////////
NrroSP  Nrro::slice(int axis, int pos)
{
   if(!_nrrd)
      return NrroSP();

   if(!dim())
      return NrroSP();

   Nrrd *nnew = nrrdNew();

   if(_kind & PROXY)
   {
      nnew->dim = dim() - 1;
      int *axesC = new int[dim() -1];
      for(int i=0; i<axis; ++i)
         axesC[i] = i;
      for(int i=axis +1; i< int(dim()) -1; ++i)
         axesC[i-1] = i;
      

      nnew->data = 0;
      nnew->type = _nrrd->type;
      nrrdCommentCopy(nnew, _nrrd);
   }
   else if(nrrdSlice(nnew, _nrrd, axis, pos))
   {
      char *berr = biffGetDone(NRRD);
      NrroErr("slice(), failed", berr);
      nrrdNuke(nnew);
      return NrroSP();
   }

   NrroSP nret(new Nrro(nnew));
   nret->setKind(nret->guessKind());

   ostringstream ss;
   ss << _name << "_slc" << axis << "-" << pos;
   nret->_name = ss.str();

   nret->setModified();
   nret->update();
   return nret;
}

/////////////////////////////////////////////////////////////////////////
// crop
/////////////////////////////////////////////////////////////////////////
SmartPtr<Nrro>  Nrro::crop(int min[], int max[])
{
   if(!_nrrd)
      return NrroSP();

   ///TODO: could make this smart enough to pad
   /// where necessary, even though this violates
   /// the semantics of "crop"
   for(int i=0; i<dim(); ++i)
   {
      if(min[i] == -1)
      {
         min[i] = 0;
      }
      else if(min[i] < 0)
      {
         NrroErr("crop(), axis min less than 0, should pad this axis");
         return NrroSP();
      }
      if(max[i] == -1)
         max[i] = dim(i)-1;
      else if(max[i] > int(dim(i))-1)
      {
         NrroErr("crop(), axis max greater than axis dim, should pad this axis");
         return NrroSP();
      }
   }

   Nrrd *nnew = nrrdNew();
   if(nrrdCrop(nnew, _nrrd, min, max))
   {
      char *berr = biffGetDone(NRRD);
      NrroErr("crop(), crop failed", berr);
      nrrdNuke(nnew);
      return NrroSP();
   }
   NrroSP nret(new Nrro(nnew));

   nret->setKind(getKind());

   ///make sure the 3D position is still correct.
   if(dim() == 3)
   {
      float xp = float(AIR_AFFINE(0,min[0],dim(0)-1,_pos.x,_pos.x+axisSize(0)));
      float yp = float(AIR_AFFINE(0,min[1],dim(1)-1,_pos.y,_pos.y+axisSize(1)));
      float zp = float(AIR_AFFINE(0,min[2],dim(2)-1,_pos.z,_pos.z+axisSize(2)));
      nret->setPos(vec3f(xp,yp,zp));
   }
   else if(dim()==4)
   {
      float xp = float(AIR_AFFINE(0,min[1],dim(1)-1,_pos.x,_pos.x+axisSize(1)));
      float yp = float(AIR_AFFINE(0,min[2],dim(2)-1,_pos.y,_pos.y+axisSize(2)));
      float zp = float(AIR_AFFINE(0,min[3],dim(3)-1,_pos.z,_pos.z+axisSize(3)));
      nret->setPos(vec3f(xp,yp,zp));
   }

   /// update the pads so the axisSizeReal is correct
   for(int i=0; i<dim(); ++i)
   {
      if(min[i] < int(_pads[i].x)) // cropped through padded part
         nret->_pads[i].x = _pads[i].x - min[i];
      else // cropped through real part
         nret->_pads[i].x = 0;
      if(max[i] > int(_pads[i].y)) //cropped through padded part
         nret->_pads[i].y = nret->dim(i) - max[i] - _pads[i].y - 1;
      else // cropped through real part
         nret->_pads[i].y = nret->dim(i) - 1;
   }    

   ostringstream ss;
   ss << _name << "_c";
   for(int i=0; i<int(dim())-1; ++i)
      ss << min[i] << "-" << max[i] << "_";
   ss << min[int(dim())-1] << "-" << max[int(dim())-1];
   nret->_name = ss.str();
  
   nret->setModified(); /// forcing an update
   nret->update();
   return nret;
}

/////////////////////////////////////////////////////////////////////////
// Pad
/////////////////////////////////////////////////////////////////////////
SmartPtr<Nrro>  Nrro::pad(int min[], int max[], int bound, double padval)
{
   if(!_nrrd)
      return NrroSP();

   /// TODO: could make this smart enough to crop where needed!
   ///  even though this violates the sematics of "pad"
   for(int i=0; i<dim(); ++i)
   {
      if(min[i] > 0)
      {
         NrroErr("pad(), min setting greater than 0, should crop this axis");
         return NrroSP();
      }
      if(max[i] < int(dim(i))-1)
      {
         NrroErr("pad(), max setting less than dim of axis, should crop this axis");
         return NrroSP();
      }
   }

   Nrrd *nnew = nrrdNew();
   if(nrrdPad_nva(nnew, _nrrd, min, max, bound, padval))
   {
      char *berr = biffGetDone(NRRD);
      NrroErr("pad(), pad failed", berr);
      nrrdNuke(nnew);
      return NrroSP();
   }
   NrroSP nret(new Nrro(nnew));

   nret->setKind(getKind());

   /// TODO: shouldn't this handle the pos for back padding?
   /// set _pads so we know what part of the data is valid
   /// keep the real size correct
   for(int i=0; i<dim(); ++i)
   {
      nret->_pads[i].x = _pads[i].x - min[i];
      nret->_pads[i].y = nret->dim(i) - 1 - (max[i] - _pads[i].y);
   }    
   
   ostringstream ss;
   ss << _name << "_p_";
   for(int i=0; i<int(dim())-1; ++i)
      ss << min[i] << "-" << max[i] << "_";
   ss << min[int(dim())-1] << "-" << max[int(dim())-1];
   nret->_name = ss.str();

   nret->setModified(); /// forcing an update
   nret->update();
   return nret;
}

///////////////////////////////////////////////////
/// pad and or crop to exact size
NrroSP Nrro::cropPad(int min[], int max[], int bound, double padval)
{
   if(!_nrrd)
      return NrroSP();

   int *cmin = new int[dim()];
   int *cmax = new int[dim()];
   int *pmin = new int[dim()];
   int *pmax = new int[dim()];

   /// figure out sizes for crop and then padd
   for(int i=0; i< int(dim()); ++i)
   {
      /// just clamp the requested dims to the dims of this guy,
      /// ie. they must be inside
      cmin[i] = clamp(int(0), min[i], int(dim(i))-1);
      cmax[i] = clamp(int(0), max[i], int(dim(i))-1);
      
      /// we croped the lower part of the axis, can't pad
      if(cmin[i] > 0)
         pmin[i] = 0;
      else /// need to pad down
         pmin[i] = min[i];
   
      /// we cropped the upper part of the axis, can't pad
      if(cmax[i] < int(dim(i))-1)
         pmax[i] = cmax[i] - pmin[i];
      else // need to padd up, figure size of croped volume, add difference from original
         pmax[i] = cmax[i] - cmin[i] + max[i] - (dim(i)-1);
   }

   NrroSP ret = crop(cmin, cmax);
   return ret->pad(pmin, pmax, bound, padval);
}

/////////////////////////////////////////////////////////////////////////
// quantize
/////////////////////////////////////////////////////////////////////////
SmartPtr<Nrro> Nrro::quantize(int bits, double min, double max)
{
   if(!_nrrd)
      return NrroSP();

   if(min != AIR_NAN)
   {
      _nrrd->min = min;
   }
   if(max != AIR_NAN)
   {
      _nrrd->max = max;
   }

   Nrrd *nnew = nrrdNew();

   if(nrrdQuantize(nnew, _nrrd, bits))
   {
      char *berr = biffGetDone(NRRD);
      NrroErr("quantize(), quantize failed", berr);
      nrrdNuke(nnew);
      return NrroSP();
   }
   
   NrroSP nret(new Nrro(nnew));
   nret->_pads = _pads;
   nret->_pos = _pos;
   nret->_kind = _kind;

   ostringstream ss;
   ss << _name << "_q" << bits;
   nret->_name = ss.str();
 
   nret->setModified();  //< forcing an update
   nret->update();
   return nret;
}

/////////////////////////////////////////////////////////////////////////
// insert Axis
/////////////////////////////////////////////////////////////////////////
int  Nrro::insertAxis(unsigned int axis)
{
   if(!_nrrd)
      return 1;

   Nrrd *nnew = nrrdNew();
   if(nrrdAxesInsert(nnew, _nrrd, axis))
   {
      char *berr = biffGetDone(NRRD);
      NrroErr("insertAxis(), insert failed", berr);
      nrrdNuke(nnew);
      return 1;
   }
   nrrdNuke(_nrrd);
   _nrrd = nnew;
   for(int i=int(dim())-1; i>int(axis); --i)
   {
      _pads[i] = _pads[i-1];
   }
   _pads[axis] = vec2ui(0,0);

   setModified();
   update();
   return 0;
}

/////////////////////////////////////////////////////////////////////////
// min max
/////////////////////////////////////////////////////////////////////////

template<class T>
vec2d minmaxFunc(Nrro::NrroIter<T> ni)
{
   if(!ni.isValid()){
      cerr << "invalid volume " << typeid(ni).name() 
           << " pos = " << ni.pos() << endl;
      return vec2d(0,0);
   }

   T min = *ni;
   T max = *ni;
   ++ni;
   while(ni.isValid())
   {
      min = g_min(*ni,min);
      max = g_max(*ni,max);
      ++ni;
   }
   return vec2d(double(min),double(max));
}

void  Nrro::updateMinMax()
{
   /// we need a reference to ourselves to use
   /// an iterator internally!!!
   _incCount();
   retDispatchIter1(_minmax,minmaxFunc,this);
   _decCount();

   _nrrd->min = _minmax.x;
   _nrrd->max = _minmax.y;
}

/////////////////////////////////////////////////////////////////////////
// resample to specific sizes
/////////////////////////////////////////////////////////////////////////
template<class TO, class TI>
bool nearestResamp(Nrro::NrroIter<TO> nii, Nrro::NrroIter<TI> nio)
{
   if(nii->dim() != nio->dim()) return true;
}

gutz::SmartPtr<Nrro>  Nrro::resample(int dims[], const NrroKernel &kern,
                                     int bound, int typeOut)
{
   if(!_nrrd)
      return NrroSP();

   /// make sure we actually have something to do
   bool go = false;

   /// create resample info struct
   NrrdResampleInfo *rsmpInfo = nrrdResampleInfoNew();

   /// check statis of all axies compared to dims
   for(int i=0; i<int(dim()); ++i)
   {
      /// if we can do something to this axis, set it up
      if(-1 != dims[i])
      {
         rsmpInfo->kernel[i] = kern.getNrrdKernel();
         rsmpInfo->samples[i] = dims[i];
         memcpy(rsmpInfo->parm[i], kern.getParams(), 
                (kern.getNrrdKernel())->numParm*sizeof(double));

         if (!AIR_EXISTS(_nrrd->axis[i].min)) {
            _nrrd->axis[i].min = 0.0;
         }
         if (!AIR_EXISTS(_nrrd->axis[i].max)) {
            _nrrd->axis[i].max = dim(i) -1;
         }

         rsmpInfo->min[i] = _nrrd->axis[i].min;
         rsmpInfo->max[i] = _nrrd->axis[i].max;
         _nrrd->axis[i].center = nrrdCenterNode;
         go = true;
      }
      /// otherwise we can't change the axis, so mark it so
      else
      {
         rsmpInfo->kernel[i] = 0;
      }
   }

   rsmpInfo->boundary = bound;
   if(typeOut == UNKNOWN)
   {
      rsmpInfo->type = getType();
   }
   else
   {
      rsmpInfo->type = typeOut;
   }

   Nrrd *nnew = nrrdNew();
   
   if(go)
   {
      if (nrrdSpatialResample(nnew, _nrrd, rsmpInfo)) 
      {
         char *berr = biffGetDone(NRRD);
         NrroErr("resample(), resample failed", berr);
         nrrdNuke(nnew);
         return NrroSP();
      }
   }
   else
   {
      if (nrrdCopy(nnew, _nrrd)) 
      {
         char *berr = biffGetDone(NRRD);
         NrroErr("resample(), resample/copy failed", berr);
         nrrdNuke(nnew);
         return NrroSP();
      }
   }

   NrroSP nret = new Nrro(nnew);
   nret->setKind(nret->guessKind());

   return nret;
}

/////////////////////////////////////////////////////////////////////////
// resample by some percent
/////////////////////////////////////////////////////////////////////////
gutz::SmartPtr<Nrro>  Nrro::resample(float pct)
{
   cerr << " Nrro::resample(float pct) not implemented yet " << endl;
   return NrroSP();
}

---

Send questions, comments, and bug reports to: