/*  Volume.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 "Volume.h"

using namespace NEWMAT;
namespace FILM {
  
  void Volume::unthreshold()
     {  
       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 VolumeSeries::Dims& pdims,const ColumnVector& in)
  {
    dims = pdims;
    preThresholdPositions = in;
    unthreshold();
  }
  
  void Volume::threshold()
     {  
       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)
     {
       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::writeAsFloat(const string& fname)
     {    
       Tracer ts("Volume::writeAsFloat");
       
       float fmin, fmax;
       const ColumnVector& outputvol = *this;
       AVW* OP = AvwOpen(fname.c_str(), "wc");

       AvwSetDim(OP,dims.x, dims.y, dims.z, 1);
       AvwSetVoxDim(OP,dims.vx, dims.vy, dims.vz, 0);
       AvwSetDataType(OP, DT_FLOAT);

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

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

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

       delete [] qv;

       AvwClose(OP);
     }

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

       const ColumnVector& outputvol = *this;
       AVW* OP = AvwOpen(fname.c_str(), "wc");
 
       AvwSetDim(OP, dims.x, dims.y, dims.z, 1);  
       AvwSetVoxDim(OP, dims.vx, dims.vy, dims.vz, 0);  
       AvwSetDataType(OP, DT_SIGNED_SHORT);

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

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

       AvwWriteVolumes(OP, qv, 1);

       delete [] qv;

       AvwClose(OP);
     }
 
   void Volume::read(const string& fname)
     {
       Tracer ts("Volume::read");
      
       AVW* IP = AvwOpen(fname.c_str(), "r");
       ColumnVector& 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;
       AvwGetDataType(IP,&type);
  
       output.ReSize(x*y*z);
  
       switch(type)
	 {
	 case DT_SIGNED_SHORT:
	   {
	     short* sbuffer=new short[imagesize];
	     AvwReadVolumes(IP,sbuffer,v);
	     
	     for(size_t j = 1; j<=(size_t)x*y*z; j++)
	       {
		 output(j)=sbuffer[j-1];
	       }
	  
	    delete[] sbuffer;
	  }
	  break;
	case DT_FLOAT:
	  {
	    float* fbuffer=new float[imagesize];
	    AvwReadVolumes(IP,fbuffer,v);

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

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

      AvwClose(IP);
    
      return;
    }
}