/**
 * @file  eig.h
 * @brief Eigen vectors and eigen values decomposition.
 *
 * Copyright (c) 2008, 2009, 2012, 2013 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 _DTIDROID_EIG_H
#define _DTIDROID_EIG_H


// #### this functions have been claimed in "nrutil.h" ####

void gaucof(int n, float a[], float b[], float amu0, float x[], float w[]);
void tred2(float **a, int n, float d[], float e[]);
int tqli(float d[], float e[], int n, float **z);
void eigsrt(float d[], float **v, int n);

// eigOfTensor is especially for Diffusion Tensor
// return 0 if result bad
int eigOfTensor(float **inMatrix,      //input squre symmetric matrix
	 float *eigenvalue,  		//allocated 1D space 4 eigenValues
	 float **eigenvector,  		//allocated 2D space 4 eigVector
	 int dim);          		//dimension

// calculate FA from 3 eigenvalues
double getFA(float *eigenvalue);

// eigOfData() will call gaucof() & tred2() & tqli() & eigsrt()
// This is copied from Songyang, revised from "C Nemuric Receipt"
void eigOfData(float **feature,   	    // matrix to be decomposed
	 float *average, 	    // return value
	 float *eigenvlaue,         // return value
	 float **eigenvector,       // return value
	 int number_of_features,    // X-dim
	 int number_of_samples);    // Y-dim

// function eig () decompose 2D matrix into eigen vectors and eigen values
// "feature" matrix:
//    +-----------------------> number_of_features, e.g. X,Y,Z coordinates
//    |				(X)
//    |
//    |
//    |
//    |
//    |
//    |
//    V
//   number_of_samples 
//   (Y)
//
//  NOTE:
//   a = matrix (1,3, 1, 5)
//  Then:
//	 a[1][1] a[1][2] a[1][3] a[1][4] a[1][5] 
//	 a[2][1] a[2][2] a[2][3] a[2][4] a[2][5] 
//	 a[3][1] a[3][2] a[3][3] a[3][4] a[3][5] 
//
//


#endif //_DTIDROID_EIG_H
