volume.cc

/*  volume.cc

    Mark Woolrich, FMRIB Image Analysis Group

    Copyright (C) 2002 University of Oxford  */

/*  CCOPYRIGHT  */

#include <iostream>
#include <cstdlib>
#include "fslio/fslio.h"

#include "newmatap.h"
#include "newmatio.h"
#include "volume.h"
#include "utils/time_tracer.h"

using namespace NEWMAT;
using namespace Utilities;

namespace MISCMATHS {
  
  void Volume::unthreshold()
     {  
       Time_Tracer ts("Volume::unthreshold");
       
       int sizeVolUnthresholded = getUnthresholdSize();
       int sizeVolThresholded = getVolumeSize();

       // Increase this to required size and set to 0
       Release(); 
       ColumnVector X=*this; 
       ReSize(sizeVolUnthresholded);
       ColumnVector::operator=(0);
       
       // Loop through volume restoring non thresholded values:
       for(int i = 1; i <= sizeVolThresholded; i++)
	 {
	   (*this)(int(preThresholdPositions(i)))=X(i);
	 }
     }
  
  void Volume::unthreshold(const VolumeInfo& pvolinfo,const ColumnVector& in)
  {
    volinfo = pvolinfo;
    preThresholdPositions = in;
    unthreshold();
  }
  
  void Volume::threshold()
  {  
    Time_Tracer ts("Volume::threshold");
    
    int sizeVol = preThresholdPositions.Nrows();
    
    ColumnVector X(sizeVol); 
    
    // Loop through pulling out non thresholded values
    for(int i = 1; i <= sizeVol; i++)
      {
	X(i)=(*this)(int(preThresholdPositions(i)));
      }
    
    (*this) = X;
  }
  
  void Volume::threshold(float thresh)
  {
    Time_Tracer ts("Volume::threshold");
    
    int size = getVolumeSize();
    int j = 0;
    
    preThresholdPositions.ReSize(size);
       
    float m = 0;
    for(int i = 1; i <= size; i++)
      {
	m = (*this)(i);
	   
	if(m > thresh)
	  {
	    j++;	
	    preThresholdPositions(j) = i;
	     
	    (*this)(j) = (*this)(i);
	  }
      }

    // Discard rest:
    *this = Rows(1,j);

    preThresholdPositions = preThresholdPositions.Rows(1,j);
  }

//   void Volume::threshold(const volume<int>& pmask)
//   {
//     Time_Tracer ts("Volume::threshold");

//     int size = getVolumeSize();
//     int j = 0;
       
//     preThresholdPositions.ReSize(size);

//     int i = 0;
//     for(int x = 0; x < pmask.xsize(); x++)
//       for(int y = 0; y < pmask.ysize(); y++)
// 	for(int z = 0; z < pmask.zsize(); z++)
// 	  {
// 	    i++;
// 	    if(pmask(x,y,z))
// 	      {
// 		j++;	
// 		preThresholdPositions(j) = i;
		   
// 		(*this)(j) = (*this)(i);
// 	      }
// 	  }

//     // Discard rest:
//     *this = Rows(1,j);

//     preThresholdPositions = preThresholdPositions.Rows(1,j);
//   }

  void Volume::writeAsFloat(const string& fname)
  {    
    Time_Tracer ts(string("Volume::writeAsFloat" + fname).c_str());
       
    float fmin, fmax;
    const ColumnVector& outputvol = *this;
    FSLIO* OP = FslOpen(fname.c_str(), "wb");

    FslSetDim(OP,volinfo.x, volinfo.y, volinfo.z, 1);
    FslSetVoxDim(OP,volinfo.vx, volinfo.vy, volinfo.vz, 0);
    FslSetDataType(OP, DT_FLOAT);
    FslSetOriginator(OP, volinfo.originator);

    int sizeVol = outputvol.Nrows();

    float *qv = new float[sizeVol];

    fmin = outputvol(1);
    fmax = outputvol(1);  
    for (int i=1; i<=sizeVol; i++) 
      { 
	if(outputvol(i) < fmin)
	  fmin = outputvol(i);
	else if(outputvol(i) > fmax)
	  fmax = outputvol(i);
	   
	qv[i-1] = outputvol(i);
      }

    FslSetMinMax(OP, (short)fmin, (short)fmax);

    //       fwrite(qv,sizeVol*sizeof(float),1,OP->imgfp);
    //fwrite(&OP->header,sizeof(OP->header),1,OP->hdrfp);
    FslWriteHeader(OP);
    FslWriteVolumes(OP, qv, 1);

    delete [] qv;

    FslClose(OP);
  }

  void Volume::writeAsInt(const string& fname)
  {    
    Time_Tracer ts("Volume::writeAsInt");
    int fmin, fmax;

    const ColumnVector& outputvol = *this;
    FSLIO* OP = FslOpen(fname.c_str(), "wb");
 
    FslSetDim(OP, volinfo.x, volinfo.y, volinfo.z, 1);  
    FslSetVoxDim(OP, volinfo.vx, volinfo.vy, volinfo.vz, 0);  
    FslSetDataType(OP, DT_SIGNED_SHORT);
    FslSetOriginator(OP, volinfo.originator);

    int sizeVol = outputvol.Nrows();

    short *qv = new short[sizeVol];
  
    fmin = (int)outputvol(1);
    fmax = (int)outputvol(1);
    for (int i=1; i<=sizeVol; i++) 
      { 
	if(outputvol(i) < fmin)
	  fmin = (int)outputvol(i);
	else if(outputvol(i) > fmax)
	  fmax = (int)outputvol(i);
	   
	qv[i-1] = (short)outputvol(i);
      }

    FslSetMinMax(OP, (short)fmin, (short)(fmax+0.9999));

    FslWriteHeader(OP);
    FslWriteVolumes(OP, qv, 1);

    delete [] qv;

    FslClose(OP);
  }
 
  void Volume::read(const string& fname)
  {
    Time_Tracer ts(string("Volume::read" + fname).c_str());
      
    FSLIO* IP = FslOpen(fname.c_str(), "rb");
    ColumnVector& output = *this;
       
    short x,y,z,v,type;
    float vx,vy,vz,tr;
    float slope, intercept;
    int doscaling;


    FslGetDim(IP,&x,&y,&z,&v);
    FslGetVoxDim(IP,&vx,&vy,&vz,&tr);
    FslGetOriginator(IP,volinfo.originator);
    doscaling = FslGetIntensityScaling(IP,&slope,&intercept);

    volinfo.x = x; volinfo.y = y; volinfo.z = z; volinfo.v = v;
    volinfo.vx = vx; volinfo.vy = vy; volinfo.vz = vz; volinfo.tr = tr;

    size_t imagesize=x*y*z;
    FslGetDataType(IP,&type);
  
    output.ReSize(x*y*z);
  
    switch(type)
      {
      case DT_SIGNED_SHORT:
	{
	  short* sbuffer=new short[imagesize];
	  FslReadVolumes(IP,sbuffer,v);
	     
	  for(size_t j = 1; j<=(size_t)x*y*z; j++)
	    {
	      if (doscaling==0) { output(j)=sbuffer[j-1]; }
	      else { output(j)=(slope * sbuffer[j-1]) + intercept; }
	    }
	  
	  delete[] sbuffer;
	}
	break;
      case DT_FLOAT:
	{
	  float* fbuffer=new float[imagesize];
	  FslReadVolumes(IP,fbuffer,v);

	  for(size_t j = 1; j<=(size_t)x*y*z; j++)
	    {
	      if (doscaling==0) { output(j)=fbuffer[j-1]; }
	      else { output(j)=(slope * fbuffer[j-1]) + intercept; }
	    }
	      
	  delete[] fbuffer;
	}
	break;
      case DT_UNSIGNED_CHAR:
	{
	  unsigned char* cbuffer=new unsigned char[imagesize];
	  FslReadVolumes(IP,cbuffer,v);

	  for(size_t j = 1; j<=(size_t)x*y*z; j++)
	    {
	      if (doscaling==0) { output(j)=cbuffer[j-1]; }
	      else { output(j)=(slope * cbuffer[j-1]) + intercept; }
	    }
	      
	  delete[] cbuffer;
	}
	break;
      default:
	perror("FslRead: DT not supported");
      }

    FslClose(IP);
    
    return;
  }
}