volumeseries.cc

/*  volumeseries.cc

    Mark Woolrich - FMRIB Image Analysis Group

    Copyright (C) 2002 University of Oxford  */

/*  CCOPYRIGHT  */

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

#include "newmatap.h"
#include "newmatio.h"
#include "volumeseries.h"
#include "miscmaths.h"
#include "utils/time_tracer.h"

using namespace NEWMAT;
using namespace Utilities;

namespace MISCMATHS {

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

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

       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(int(preThresholdPositions(i)))=X.Column(i);
	 }
    }

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

       int numSeries = preThresholdPositions.Nrows();
       int numVolumes = getNumVolumes();

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

       *this = X;
    }
  
   void VolumeSeries::thresholdSeries(float thresh, bool removeMean)
     {
       Time_Tracer ts("VolumeSeries::thresholdSeries");

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

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

       preThresholdPositions.ReSize(numSeries);
       
       float m = 0,s = 0;
       for(int i = 1; i <= numSeries; i++)
	 {
	   m = MISCMATHS::mean(ColumnVector(getSeries(i))).AsScalar();
	   s = MISCMATHS::var(ColumnVector(getSeries(i))).AsScalar();
	   /*
	   if(m > thresh && s == 0.0)
	     {
	       cerr << "m = " << m << endl;
	       cerr << "s = " << s << endl;
	       cerr << "i = " << i << endl;
	       cerr << "j = " << j+1 << endl;
	     }
	   */
	   if(m > thresh && s > 0.0)
	     {
	       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);
     }

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

  void VolumeSeries::read(const string& fname)
    {
      Time_Tracer ts(string("VolumeSeries::read-" + fname).c_str());

      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);
      AvwGetOriginator(IP,volinfo.originator);

      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*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 = volinfo.x*volinfo.y*volinfo.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 = volinfo.x*volinfo.y*volinfo.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 = volinfo.x*volinfo.y*volinfo.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");
	}

      AvwClose(IP);

      return;
    }
  
  void VolumeSeries::unthresholdSeries(const VolumeInfo& pvolinfo,const ColumnVector& in)
  {
    volinfo = pvolinfo;
    preThresholdPositions = in;
    unthresholdSeries();
  }

  void VolumeSeries::writeAsFloat(const string& fname)
    {
       Time_Tracer ts(string("VolumeSeries::writeAsFloat" + fname).c_str());

       AVW* OP = AvwOpen(fname.c_str(), "wc");

       AvwSetDim(OP,volinfo.x, volinfo.y, volinfo.z, volinfo.v);
       AvwSetVoxDim(OP,volinfo.vx, volinfo.vy, volinfo.vz, volinfo.tr);
       AvwSetDataType(OP, DT_FLOAT);
       AvwSetOriginator(OP, volinfo.originator);

       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);
       
    }

}