/*  volumeseries.h

    Mark Woolrich - FMRIB Image Analysis Group

    Copyright (C) 2002 University of Oxford  */

/*  CCOPYRIGHT  */

#include <iostream>
#include <fstream>
#define WANT_STREAM
#define WANT_MATH

#include "newmatap.h"
#include "newmatio.h"
#include "volumeinfo.h"

#include <string>

using namespace NEWMAT;

namespace MISCMATHS {
  
#if !defined(__volumeseries_h)
#define __volumeseries_h

  // Rows are volumes
  // Columns are (time) series
  // num Rows is size of (time) series
  // num Cols is size of volumes
  class VolumeSeries : public Matrix
    {
    public:

      VolumeSeries() : Matrix() {}
      
      VolumeSeries(const VolumeInfo& pvolinfo, const ColumnVector& in) : 
	Matrix(),
	volinfo(pvolinfo),
	preThresholdPositions(in)
	{}

      VolumeSeries(int pnumVols, int pnumSeries) : 
	Matrix(pnumVols, pnumSeries), 
	means(pnumSeries){}
      
      VolumeSeries(int pnumVols, int pnumSeries, const VolumeInfo& pvolinfo, const ColumnVector& in) : 
	Matrix(pnumVols, pnumSeries), 
	volinfo(pvolinfo),
	preThresholdPositions(in),
	means(pnumSeries){}

      VolumeSeries& operator=(const VolumeSeries& vol) { 
	Matrix::operator=(vol); 
	preThresholdPositions = vol.preThresholdPositions;
	volinfo = vol.volinfo;
	means = vol.means;
	return *this;
      }
      VolumeSeries& operator=(const Matrix& mat) {
	Matrix::operator=(mat);
	return *this;
      }

      VolumeSeries& operator=(float pin) { 
	Matrix::operator=(pin); 
	return *this;
      }
      VolumeSeries(const VolumeSeries& vol){operator=(vol);}
      VolumeSeries(const Matrix& pmat) : Matrix(pmat) {}
      
      void thresholdSeries(float thresh, bool removeMean);
      void thresholdSeries();
      void unthresholdSeries();
      void unthresholdSeries(const VolumeInfo& pvolinfo,const ColumnVector& in);
      void removeSeriesMeans();

      const ColumnVector& getPreThresholdPositions() const { return preThresholdPositions; } 
      void setPreThresholdPositions(const ColumnVector& in) { preThresholdPositions = in; }

      int getNumVolumes() const { return Nrows(); }
      int getNumSeries() const { return Ncols(); }
      int nvoxels() const { return Ncols(); }
      int tsize() const { return Nrows(); }

      const ReturnMatrix getSeries(int i) const { ColumnVector tmp = Column(i);tmp.Release();return tmp; }

      ReturnMatrix getSeries(int i) { ColumnVector tmp = Column(i);tmp.Release();return tmp; }

      void setSeries(const ColumnVector& in, int i) { Column(i) = in; }

      const ReturnMatrix getVolume(int i) const { RowVector tmp = Row(i);tmp.Release();return tmp; }
      ReturnMatrix getVolume(int i) { RowVector tmp = Row(i);tmp.Release();return tmp; }

      void setVolume(const ColumnVector& in, int i) { Row(i) = in.AsRow(); }
      void setVolume(const RowVector& in, int i) { Row(i) = in; }

      void read(const std::string& fname);
      void writeAsInt(const std::string& fname);
      void writeAsFloat(const std::string& fname);
      void writeThresholdedSeriesAsFloat(const VolumeInfo& pvolinfo,const ColumnVector& in,const std::string& fname);
      
      void replaceMeans();

      const VolumeInfo& getInfo() const { return volinfo; }
      void setInfo(const VolumeInfo& pvolinfo) { volinfo = pvolinfo; }

      int getUnthresholdNumSeries() const { return volinfo.x*volinfo.y*volinfo.z; }
      
    protected:
      VolumeInfo volinfo;
      ColumnVector preThresholdPositions; 
      ColumnVector means;
    };

#endif

}