volumeseries.cc

/*  volumeseries.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 "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 > 1e-10)
	     {
	       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());

      FSLIO* IP = FslOpen(fname.c_str(), "rb");
      Matrix& 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);
      FslGetIntent(IP, &(volinfo.intent_code), &(volinfo.intent_p1), &(volinfo.intent_p2), 
		   &(volinfo.intent_p3));
      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;

      volinfo.miscinfo = FslInit();
      FslCloneHeader(volinfo.miscinfo,IP);

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

      FslClose(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());

       FSLIO* OP = FslOpen(fname.c_str(), "wb");

       FslCloneHeader(OP,volinfo.miscinfo);
       
       FslSetCalMinMax(OP,Minimum(),Maximum());
       FslSetDim(OP,volinfo.x, volinfo.y, volinfo.z, volinfo.v);
       FslSetVoxDim(OP,volinfo.vx, volinfo.vy, volinfo.vz, volinfo.tr);
       FslSetDataType(OP, DT_FLOAT);
       FslSetIntent(OP, volinfo.intent_code, volinfo.intent_p1, volinfo.intent_p2, 
		    volinfo.intent_p3);

       int volSize = getNumSeries();
       int volNum = getNumVolumes();

       // This should assert() but I'm too scared to add that to production code!
       if ( (volinfo.x * volinfo.y * volinfo.z != volSize) || (volinfo.v != volNum) )
         cerr << "WARNING: Dimensions are inconsistent in VolumeSeries::writeAsFloat()\n"
              << "x, y, z, v = " << volinfo.x << ", " << volinfo.y << ", " << volinfo.z
              << ", " << volinfo.v << "\nvolSize, volNum = " << volSize << ", " 
              << volNum << "\nThis is probably NOT what you intended!!" << endl;

       FslWriteHeader(OP);

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

       delete [] qv;

       FslClose(OP);
       
    }

  void VolumeSeries::writeThresholdedSeriesAsFloat(const VolumeInfo& pvolinfo,const ColumnVector& in,const string& fname)
    {
       volinfo = pvolinfo;
       preThresholdPositions = in;

       Time_Tracer ts(string("VolumeSeries::writeThresholdedSeriesAsFloat" + fname).c_str());

       FSLIO* OP = FslOpen(fname.c_str(), "wb");

       FslCloneHeader(OP,volinfo.miscinfo);

       FslSetCalMinMax(OP,Minimum(),Maximum());
       FslSetDim(OP,volinfo.x, volinfo.y, volinfo.z, volinfo.v);
       FslSetVoxDim(OP,volinfo.vx, volinfo.vy, volinfo.vz, volinfo.tr);
       FslSetDataType(OP, DT_FLOAT);
       FslSetIntent(OP, volinfo.intent_code, volinfo.intent_p1, volinfo.intent_p2, volinfo.intent_p3);

       int volSize = getUnthresholdNumSeries();
       int numThresholdedSeries = getNumSeries();
       int volNum = getNumVolumes();

       FslWriteHeader(OP);

       float *qv = new float[volSize];
      
       for(int i = 1; i<= volNum; i++)
	 {
	   for(int j = 1; j <= volSize; j++)
	     qv[j-1]=0;
	   for(int j = 1; j <= numThresholdedSeries; j++)
	     qv[int(preThresholdPositions(j))-1] = (*this)(i,j);
	   FslWriteVolumes(OP, qv, 1);
	 }

       delete [] qv;

       FslClose(OP);
       
    }

}