VolumeSeries.cc

/*  VolumeSeries.cc

    Mark Woolrich, FMRIB Image Analysis Group

    Copyright (C) 1999-2000 University of Oxford  */

/*  CCOPYRIGHT  */

#include <iostream>
#include <cstdlib>
#include <avwio.h>

#include "newmatap.h"
#include "newmatio.h"
#include "VolumeSeries.h"
#include "miscmaths.h"

using namespace NEWMAT;
namespace TACO {

  void VolumeSeries::replaceMeans()
  {
    Tracer ts("VolumeSeries::replaceMeans");
    for(int i = 1; i <= getNumSeries(); i++)
      {
	Column(i) = Column(i) + means(i);
      }
  }

  void VolumeSeries::unthresholdSeries()
    {
       Tracer ts("VolumeSeries::unthresholdSeries");

       cerr << "Unthresholding Series... ";

       int numUnThresholdedSeries = getUnthresholdNumSeries();
       int numThresholdedSeries = getNumSeries();
       int numVolumes = getNumVolumes();
      
       // Increase this to required size and set to 0
       Release();        
       Matrix X=*this; 
       ReSize(numVolumes, numUnThresholdedSeries);
       Matrix::operator=(0);
       
       // Loop through restoring non thresholded series: 
       for(int i = 1; i <= numThresholdedSeries; i++)
	 {
	   Column(preThresholdPositions(i))=X.Column(i);
	 }

       cerr << "Completed" << endl;
    }

  void VolumeSeries::thresholdSeries()
    {
       Tracer ts("VolumeSeries::thresholdSeries");

       cerr << "Thresholding Series... ";
       int numSeries = preThresholdPositions.Nrows();
       int numVolumes = getNumVolumes();

       Matrix X(numVolumes, numSeries); 
      
       for(int i = 1; i <= numSeries; i++)
	 {
	   X.Column(i)=Column(preThresholdPositions(i));
	 }

       *this = X;
       cerr << "Completed" << endl;
    }
  
   void VolumeSeries::thresholdSeries(float thresh, bool removeMean)
     {
       Tracer ts("VolumeSeries::thresholdSeries");

       cerr << "Applying threshold to Series... ";

       int numSeries = getNumSeries();
       int j = 0;

       if(removeMean)
	 {
	   means.ReSize(numSeries);
	   means = 0;
	 }

       preThresholdPositions.ReSize(numSeries);
       
       float m = 0;
       for(int i = 1; i <= numSeries; i++)
	 {
	   m = MISCMATHS::mean(getSeries(i));
	   
	   if(m > thresh)
	     {
	       j++;	
	       preThresholdPositions(j) = i;
	     
	       if(removeMean)
		 {
		   Column(j) = getSeries(i) - m;
		   means(i) = m;
		 }
	       else
		 Column(j) = getSeries(i);
	     }
	 }

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

       preThresholdPositions = preThresholdPositions.Rows(1,j);

       cerr << "Completed" << endl;
     }

  void VolumeSeries::removeSeriesMeans()
    {
      int numSeries = getNumSeries();
    
      for(int i = 1; i <= numSeries; i++)
	{
	  float m = MISCMATHS::mean(getSeries(i));
	  Column(i) = getSeries(i) - m;
	}
    }

  void VolumeSeries::read(const string& fname)
    {
      Tracer ts("VolumeSeries::read");
      cerr << "Loading image...";

      AVW* IP = AvwOpen(fname.c_str(), "r");
      Matrix& output = *this;

      short x,y,z,v,type;
      float vx,vy,vz,tr;

      AvwGetDim(IP,&x,&y,&z,&v);
      AvwGetVoxDim(IP,&vx,&vy,&vz,&tr);

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

      size_t imagesize=x*y*z*v;
      AvwGetDataType(IP,&type);
  
      output.ReSize(v,x*y*z);
  
      switch(type)
	{
	case DT_SIGNED_SHORT:
	  {
	    short* sbuffer=new short[imagesize];
	
	    AvwReadVolumes(IP, sbuffer,v);
	
	    size_t volsize = dims.x*dims.y*dims.z;
	    size_t volstart = 1;

	    for(size_t i=1; i<=(size_t)v; i++)
	      {
		volstart = (i-1)*volsize;
		for(size_t j = 1; j<=(size_t)x*y*z; j++)
		  {
		    output(i,j)=sbuffer[volstart+j-1];
		  }
	      }
	  
	    delete[] sbuffer;
	  }
	  break;
	case DT_FLOAT:
	  {
	    float* fbuffer=new float[imagesize];
	    AvwReadVolumes(IP,fbuffer,v);

	    size_t volsize = dims.x*dims.y*dims.z;
	    size_t volstart = 1;
	    for(size_t i=1; i<=(size_t)v; i++)
	      {
		volstart = (i-1)*volsize;
		for(size_t j = 1; j<=(size_t)x*y*z; j++)
		  {
		    output(i,j)=fbuffer[volstart+j-1];
		  }
	      }
	    delete[] fbuffer;
	  }
	  break;
	case DT_UNSIGNED_CHAR:
	  {
	    unsigned char* cbuffer=new unsigned char[imagesize];
	    AvwReadVolumes(IP,cbuffer,v);

	    size_t volsize = dims.x*dims.y*dims.z;
	    size_t volstart = 1;
	    for(size_t i=1; i<=(size_t)v; i++)
	      {
		volstart = (i-1)*volsize;
		for(size_t j = 1; j<=(size_t)x*y*z; j++)
		  {
		    output(i,j)=cbuffer[volstart+j-1];
		  }
	      }
	    delete[] cbuffer;
	  }
	  break;
	default:
	  perror("AvwRead: DT not supported");
	}

      cerr << " Completed" << endl;

      AvwClose(IP);

      return;
    }

  void VolumeSeries::writeAsFloat(const string& fname)
    {
       Tracer ts("VolumeSeries::writeAsFloat");

       cerr << "Saving image...";
 
       AVW* OP = AvwOpen(fname.c_str(), "wc");

       AvwSetDim(OP,dims.x, dims.y, dims.z, dims.v);
       AvwSetVoxDim(OP,dims.vx, dims.vy, dims.vz, dims.tr);
       AvwSetDataType(OP, DT_FLOAT);
       
       int volStart = 1;
       int volSize = getNumSeries();
       int volNum = getNumVolumes();

       float *qv = new float[volSize*volNum];

       for(int i = 1; i<= volNum; i++)
	 {
	   volStart = (i-1)*volSize;
	   for(int j = 1; j <= volSize; j++)
	     {
	       qv[volStart+j-1] = (*this)(i,j);
	     }
	 }

       AvwWriteVolumes(OP, qv, volNum);

       delete [] qv;

       AvwClose(OP);
       
       cerr << " Completed" << endl;
    }

}