/**
 * @file  mVolume.h
 * @brief Volume of warping matrix, i.e. a matrix at each voxel, defining a
 *        transformation matrix.
 *
 * 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_MVOLUME_H
#define _SBIA_DTIDROID_MVOLUME_H


#include "Global.h"
#include "Basics.h"
#include "nrutil.h"
#include "WarpVolume.h"
#include "eig.h"

//--------------------------------------------------------------
// class mVolume: Volume of warping matrix, i.e. a matrix at each
//		  voxel, defining a transformation matrix
//
//   warp volume dimensions:  X = 256, Y = 256, Z=?
// Members:
//    dimX,Y,Z: dimensions of the volume, file size would be
//              dimX*dimY*dimZ*(6*4): 6 float components at a voxel
//
//
// Memeber functions:
//   - getTransMatrixOld(WarpVolume* &wV): Procrustean Method
//   - getTransMatrix(WarpVolume* &wV): Weighted Procrustean Method
//	 initialize the tranformation matrix volume based on
//	 the warping volume (a vector field)
//
//   - getMatrix(int ind)
//   - getMatrix(int xx,int yy, int zz)
//	 retrieve matrix at (xx,yy,zz) or index in the whole array
//   -int  weightedShapeAt(int col, int row, int slice,VoxMap* &wV, 
//		DTVolume* TensorVol , Ellipsoid& E, int& neib);
//	   iterate at (col,row,slice) to get an Ellipsoid which 
//	     represents for the weighted Procrustean shape
//	   return 0 if the Ellipsoid is not defined
//	   return 1 if the Ellipsoid is properly defined
//
//   -void getProcrusteanAB(int col, int row, int slice, VoxMap* &wV, 
//	    VoxMap* &PD, int* records, int neib, float **A, float **B,
//	    int & nAB);
//	   based on the Procrustean weighted Ellipsoid shape
//	   to form the A & B of A = U*B
//	   and return 'nAB' to tell the actual dimension of **A & **B 
//	   which will be of (3 x nAB): 3 rows & nAB columns
//   -int  weightedShapeAt(int col, int row, int slice,VoxMap* &wV, 
//		DTVolume* TensorVol , int* records, int& neib);
//	   iterate at (col,row,slice) to get an Ellipsoid shape which 
//	     represents for the weighted Procrustean shape
//	   the result is recorded in an array "records", which marks
//	     in the neighborhood the grid points located inside in the
//	     Ellipsoid.
//	   return 0 if the Ellipsoid is not defined
//	   return 1 if the Ellipsoid is properly defined
//
//   -void getProcrusteanAB(int col, int row, int slice, VoxMap* &wV, 
//	    VoxMap* &PD, int* records, int neib, float **A, float **B,
//	    int & nAB);
//	   based on the Procrustean weighted Ellipsoid shape
//	     to form the A & B of A = U*B;
//	   the Ellipsoid shape is expressed via an int array "records"
//	     which is generated by  weightedShapeAt()
//	   and return 'nAB' to tell the actual dimension of **A & **B 
//	   which will be of (3 x nAB): 3 rows & nAB columns

///------------------------------------------------------------------------------

class mVolume{
  int dimX,dimY,dimZ;
  Matrix * mv;
  int CPUs; //how many CPUs to be used

public:
   mVolume();
   mVolume(int xx, int yy, int zz);
  ~mVolume();

   void getSpace();

   //int  weightedShapeAt(int col, int row, int slice,VoxMap* &wV, 
		//DTVolume* TensorVol , Ellipsoid& E, int& neib);
   int  weightedShapeAt(int col, int row, int slice,VoxMap* &wV, 
		DTVolume* TensorVol , int *records, int& neib,
 		float **F, float *eVal, float **eVec); //, float **T);
   int  weightedShapeAt_withOptE(int col, int row, int slice,VoxMap* &wV, 
		DTVolume* TensorVol , int *records, int& neib,
 		float **F, float *eVal, float **eVec,Ellipsoid& optE); 
		//optE is the estimated (optimized)orientation at (col,row,slice)

  //void getProcrusteanAB(int col, int row, int slice, VoxMap* &wV, VoxMap* &PD,
	    //Ellipsoid E, int neib, float **A, float **B,int & nAB);
   void getProcrusteanAB(int col, int row, int slice, VoxMap* &wV, VoxMap* &PD,
	    int *records, int neib, float **A, float **B,int & nAB);
   void getProcrusteanAB4PD2(int col,  int row, int slice, VoxMap* &wV, 
	    Vector warpedPD1, VoxMap* &PD2, int *records, int neib, 
	    float **A, float **B,int & nAB);

   void getTransMatrix(VoxMap* &wV, DTVolume* TensorVol, VoxMap* &PD,VoxMap* &PD2, 
			DTVolume* &optDTI,int useOptimized=1);
   // get transformation matrix: new version

   void getTransMatrixOld(VoxMap* &wV); //get transformation matrix: old version

   void loadIdentity(); //reset all matrices to Identity I
   Matrix & getMatrix(int ind);
   Matrix & getMatrix(int xx,int yy, int zz);

};


#endif // _SBIA_DTIDROID_MVOLUME_H
