/**
 * @file  Fields.h
 * @brief Representation of fields of different types.
 *
 * Copyright (c) 2008, 2009, 2012 University of Pennsylvania.
 *
 * This file is part of DTI-DROID.
 *
 * DTI-DROID is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * DTI-DROID is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with DTI-DROID.  If not, see <http://www.gnu.org/licenses/>.
 *
 * Contact: SBIA Group <sbia-software at uphs.upenn.edu>
 */

#ifndef _SBIA_DTIDROID_FIELDS_H
#define _SBIA_DTIDROID_FIELDS_H


#include "FieldDRTG.h"
//#include "Ellipsoid.h" 
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_blas.h>
#include <math.h>

struct Pt3d{
  float x,y,z;
};


class ScalarField : public Field<float, 1>
{
 public:
  void ComputeSmoothingFilter(float);
  void ComputeSmoothingFilter(float sigma, float cen_x, float cen_y, float cen_z);
  void Smooth(ScalarField *smoothingFilter, ScalarField *out);
};

class VectorField : public Field<float, 3>
{
 public:
  void AddOffset();
  void SubtractOffset();
  void FlipComp(int);
  void Normalize();
  int getVectorAtAnyPosition(float px, float py, float pz, struct Pt3d & P);
  void reverseWarpField(VectorField *out);
  void reverseWarpFieldHAMMER(VectorField *out);
  void Smooth(ScalarField *smoothingFilter, VectorField *out);
};

class TensorField;

/*--
class RotationField : public Field<float, 9>
{
 public:
  void getTransMatrix(VectorField* wV, TensorField* TensorVol, VectorField* PD, VectorField *PD2,  TensorField* optDTI, int useOptimized);
  void getProcrusteanAB(int col, int row, int slice, VectorField* wV, VectorField* PD, int *records, int neib, float **A, float **B,int & nAB);
  int  weightedShapeAt_withOptE(int col, int row, int slice, VectorField *wV, 
		TensorField* TensorVol , int *records, int& neib,
 		float **F, float *eVal, float **eVec, Ellipsoid& optE); 
};
--*/

/* note that upon occasion, symmetric but not positive-definite matrices
   may also be stored in objects of the TensorField class */
class TensorField : public Field<float, 6>
{
 public:
  void ComputeTrace(ScalarField *traceField);
  void ComputeFA(ScalarField *faVol);
  void ComputePD(VectorField *PDVol, int d_index, int weight_flag);
  void ComputeEigD(ScalarField *e1Field,ScalarField *e2Field,ScalarField *e3Field,VectorField *PD1Field,VectorField *PD2Field,VectorField *PD3Field);
  void ComputeEigD(ScalarField *e1Field,ScalarField *e2Field,ScalarField *e3Field,VectorField *PD1Field,VectorField *PD2Field,VectorField *PD3Field, ScalarField *traceField,ScalarField *faField, int FAweighting);
  void logESmooth(ScalarField *smoothingFilter, TensorField *dti_out);
  void ESmooth(ScalarField *smoothingFilter, TensorField *dti_out);
  void ESmooth(ScalarField *smoothingFilter, TensorField *dti_out, unsigned char *mask) ;
  void logTField(TensorField *dti_out, unsigned char *mask);
  void expTField(TensorField *dti_out, unsigned char *mask);
  void computeNZmask(unsigned char *mask, float threshold_factor);
  float traceTensorThreshold;  

  void reOrientTensorFieldFS(VectorField *warpField, TensorField *dti_out);

  // this function is not yet fully implemented
  void reOrientTensorFieldPE(VectorField *warpField, TensorField *dti_out);
  
  // these are two slow, but accurate functions
  void displaceTensorFieldUsingFwdField(VectorField *warpField, TensorField *dti_out, float sigma, float filter_size_red_factor);
  void displaceTensorFieldUsingInvField(VectorField *warpField, TensorField *dti_out, float sigma, float filter_size_red_factor);

  // this is a faster function
  void displaceTensorFieldUsingInvFieldTL(VectorField *warpField, TensorField *dti_out);
  
  // don't use these four routines-- they cause interpolation artifacts
  void displaceTensorFieldUsingFwdFieldFixedKernel(VectorField *warpField, TensorField *dti_out, float sigma, float filter_size_red_factor);
  void displaceTensorFieldUsingInvFieldFixedKernel(VectorField *warpField, TensorField *dti_out, float sigma, float filter_size_red_factor);
  void displaceTensorFieldUsingFwdFieldFixedKernelHybrid(VectorField *warpField, TensorField *dti_out, float sigma, float filter_size_red_factor);
  void displaceTensorFieldUsingInvFieldFixedKernelHybrid(VectorField *warpField, TensorField *dti_out, float sigma, float filter_size_red_factor);


};


#endif // _SBIA_DTIDROID_FIELDS_H
