/**
 * @file  matrixSHENCLASSIC.h
 * @brief Linear algebra functions copied from CLASSIC.
 *
 * 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_MATRIXSHENCLASSIC_H
#define _SBIA_DTIDROID_MATRIXSHENCLASSIC_H


#include <stdio.h>
#include <stdlib.h>
#include <math.h>



/* Nov 2001, for Head brain image */
typedef struct 
{  
  unsigned char Edge;  
  unsigned char Tiss;  
  unsigned char Geom;   
  unsigned char BGvlm;  
  unsigned char CSFvlm;  
  unsigned char VNvlm;  
} HeadImgAttribute ;  



/* June 2001, for skull-stripped brain image */
typedef struct
{
  unsigned char Edge;
  unsigned char Tiss;
  unsigned char Geom; 
  unsigned char VNvlm;
  unsigned char CSFBG;
} ImgAttribute ;



typedef struct MatrixStruct 
{
	double **data;
	int height, width;
} Matrix;

typedef struct uc_MatrixStruct 
{
	unsigned char **data;
	int height, width;
} uc_Matrix;

typedef struct i_MatrixStruct 
{
	int **data;
	int height, width;
} i_Matrix;

#define TRUE     1
#define FREE_ARG char*
#define NR_END 1
#define PI 3.141592653589793115997963468544185161590576171875

static float tempr;
#define SWAP(a,b) {tempr=(a);(a)=(b);(b)=tempr;}

#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))

static float maxarg1,maxarg2;
#define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ? (maxarg1) : (maxarg2))

static float minarg1,minarg2;
#define FMIN(a,b) (minarg1=(a),minarg2=(b),(minarg1) < (minarg2) ?\
        (minarg1) : (minarg2))

static int iminarg1,iminarg2;
#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ?  (iminarg1) : (iminarg2))

static float sqrarg;
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)




void   nrerrorSHEN(const char *error_text);
double *dvectorSHEN(int nl, int nh);
void   free_dvectorSHEN(double *v, int nl, int nh);
float  *vectorSHEN(int nl, int nh);
void   free_vectorSHEN(float *v, int nl, int nh);
float  **matrixSHEN(int nrl,int nrh,int ncl,int nch);
void   free_matrixSHEN(float **m,int nrl,int nrh,int ncl,int nch);
double **dmatrixSHEN(int nrl, int nrh, int ncl, int nch);
void   free_dmatrixSHEN(double **m, int nrl, int nrh, int ncl, int nch);
double log2(double a);
void   sort(double *Y, int *I, double *A, int length);
void   minimun(double *Y, int *I, double *A, int length);
void   Mat_Abs(Matrix *A);
void   Mat_Mean(double *mean, Matrix *A);
void   Mat_Variance(double *variance, Matrix *A) ;
void   Mat_Vector(Matrix *A, float *a);
void   Mat_Shift(Matrix *A, Matrix *B, int side);
void   Mat_Zeros(Matrix *A);
void   Mat_Zeros_uc(uc_Matrix *A);
void   Mat_Zeros_i(i_Matrix *A);
void   CreateMatrix(Matrix **M, int hei, int wid);
void   new_matrix(Matrix *tmp, int row, int col);
void clear_matrix(Matrix *tmp);
void   FreeMatrix(Matrix *M);
void   Create_i_Matrix(i_Matrix **M, int hei, int wid);
void   Free_i_Matrix(i_Matrix *M);
void   Create_uc_Matrix(uc_Matrix **M, int hei, int wid);
void   Free_uc_Matrix(uc_Matrix *M);
void   Mat_FFT2(Matrix *Output_real, Matrix *Output_imag, Matrix *Input_real, Matrix *Input_imag);
void   Mat_IFFT2(Matrix *Output_real, Matrix *Output_imag, Matrix *Input_real, Matrix *Input_imag);
void   four2(double **fftr, double **ffti, double **rdata, double **idata, int rs, int cs, int isign);
void   four1(double *data, int nn, int isign);
void   Mat_Copy(Matrix *A, Matrix *B, int h_target, int w_target, int h_begin, int w_begin, int h_end,int w_end);
void   Mat_uc_Copy(uc_Matrix *A, uc_Matrix *B, int h_target, int w_target, int h_begin, int w_begin,int h_end, int w_end);
void   Mat_i_Copy(i_Matrix *A, i_Matrix *B, int h_target, int w_target, int h_begin, int w_begin, int h_end, int w_end);
void   Mat_Product(Matrix *A, Matrix *B, Matrix *C);
void   Mat_Sum(Matrix *A, Matrix *B, Matrix *C);
void   Mat_Substract(Matrix *A, Matrix *B, Matrix *C);
void   Mat_Fliplr(Matrix *A);
void   Mat_Flipud(Matrix *A);
void   Mat_uc_Fliplr(uc_Matrix *A);
void   Mat_uc_Flipud(uc_Matrix *A);

/*by SHEN in JHU*/
int    *ivectorSHEN(long nl, long nh) ;
void   free_ivectorSHEN(int *v, long nl, long nh) ;
void   Mat_Inverse(Matrix *A, Matrix *B) ; /*A in, B out*/
int    Mat_Inverse_Singular(Matrix *A, Matrix *B); /*A in, B out*/
void   Mat_times_Vector(float *Vout, Matrix *A, float *Vin) ;
void   Mat_A_equal_BxC(Matrix *A, Matrix *B, Matrix *C) ;
void   Mat_EqualCopy(Matrix *A, Matrix *B) ;
void   Mat_Print(Matrix *A) ;
void   Mat_Calculate_EigenVectors_EigenValues(Matrix *C, float *EigenValue, Matrix *EigenVector, int PRNorNOT) ;
void   svdcmp(float **a, int m, int n, float w[], float **v) ;
void   vector_Print(float *v, int size) ;
float  gasdev(long *idum) ;

Matrix * Mat_times_Mat(Matrix *A, Matrix *B);
void   Copy_Mat(Matrix *A, Matrix *B);

/* June 2001 */
ImgAttribute ***ImgAttributealloc3d(int i_size,int j_size,int k_size) ;
ImgAttribute *ImgAttributealloc1d(int k_size) ;

/* June 2001 */
HeadImgAttribute ***HeadImgAttributealloc3d(int i_size,int j_size,int k_size) ;
HeadImgAttribute *HeadImgAttributealloc1d(int k_size) ;


typedef unsigned char   byte;                   //!<  8 bit unsigned
typedef unsigned short  unshort;

#ifdef M_PI
#  undef M_PI
#endif
#define M_PI 3.14159265358979323846
#define M_2PI M_PI * 2

//#ifndef COMPLEX_DEFINED

typedef struct 
{
  float re;
  float im;
} complex;

//#define COMPLEX_DEFINED
//#endif

typedef struct
{
  int X;
  int Y;
  int Z;
  int D;
  int N;
  int *No_Of_Features;
  byte ***img;
  byte ***grd;
  byte ***tissue;
  byte ***edge;
  byte ***geom;
  byte ***flag;
  byte ****vbox;
} AVBox;

typedef struct
{
  int Size;
  short *sig;
} Short_Signal;

typedef struct
{
  int X;
  int Y;
  int Z;
  Short_Signal **vector;
} Vector_Series;

typedef struct
{
  int X;
  int Y;
  int Z;
  byte ***box;
} BBox;

typedef struct
{
  int X;
  int Y;
  int Z;
  float dx;
  float dy;
  float dz;
  unshort ***img;
} XZVolume;

typedef struct
{
  int X;
  int Y;
  int Z;
  float dx;
  float dy;
  float dz;
  float ***box;
} DBox;
typedef struct
{
  int X;
  int Y;
  float **img;
} Image;

typedef struct
{
  int Size;
  float *sig;
} Signal;

typedef struct
{
  int size;
  double *value;
} Vector;

typedef struct
{
  int Size;
  short *sig;
} SSignal;
typedef struct
{
  int X;
  int Y;
  int Z;
  int D;
  byte ***flag;
  byte ***img;
  float ****vbox;
} VectorBox;
typedef struct
{
  int X;
  int Y;
  int Z;
} Point;

typedef struct
{
  short X;
  short Y;
  short Z;
} ShortPoint;

typedef struct
{
  float X;
  float Y;
  float Z;
} FPoint;

typedef struct
{
  FPoint *point;
  int Num;
} FPoints;


typedef struct
{
  Point *point;
  int Num;
} Points;

typedef struct
{
  Point *point;
  int Num;
  float *degree;
} DegreePoints;
typedef struct
{
  int X;
  int Y;
  int Z;
  short ***box;
} SBox;

typedef struct
{
  int X;
  int Y;
  int Z;
  int BX;
  int BY;
  int BZ;
  int dx;
  int dy;
  int dz;
  Point ***point;
} PointArray3D;

typedef struct
{
  int X;
  int Y;
} Point_2D;

typedef struct
{
  int X;
  int Y;
  int BX;
  int BY;
  int dx;
  int dy;
  Point_2D **point;
} PointArray2D;

typedef struct Focus_Point
{
  struct Focus_Point *previous;
  struct Focus_Point *next;
  Point point;
  float consistency;
  float actual_dist;
} Focus_Point;

typedef struct
{
  int X;
  int Y;
  int Z;
  int D;
  int N;
  int *No_Of_Features;
  byte ***img; // original image intensity
  byte ***tissue; // tissue = 0 is background, tissue = 1 is image data;
  byte ***edge;  // edge type, 0: no edge, 1: between image and background, 2: in-image edge
  byte ***grd; // gradient magnitude of original image;
  float ***geom; // [0,1] percentage of image volume within a neighborhood, 1-geom is percentage of background (CSF,VN,BG) within the neighborhood;
  short ****dxdydz; // gradient directions, now no use.
  byte ***flag; // have wav data.
  float ****vbox; // wav
} FAVBox;

typedef struct
{
  int Num;
  int *neighbor_pts;
  float *distance_m;
  float *distance_s;
} Neighbor_PTS;

typedef struct
{
  int Num;
  Neighbor_PTS *n_pts;
  float neighborhood_size;
} All_Neighbor_PTS;

typedef struct Point_Link
{
  struct Point_Link *previous;
  struct Point_Link *next;
  int index;
  Point point;
} Point_Link;

typedef struct
{
  int X;
  int Y;
  int Z;
  ShortPoint ***point;
} Short_Deformation_Field;

typedef struct
{
  int X;
  int Y;
  int Z;
  FPoint ***point;
} Deformation_Field;
typedef struct
{
  int X;
  int Y;
  int Z;
  int T;
  byte ****data;
} Image_4D;

typedef struct
{
  Matrix *A;
  float  Tx[3];
  float  Ty[3];
} Transformation;

typedef struct
{
  int X;
  int Y;
  int Z;
  int T;
} Point_4D;

typedef struct
{
  Point_4D *point;
  int Num;
} Points_4D;

typedef struct
{
  int X;
  int Y;
  int Z;
  int D;
  int N;
  int *No_Of_Features;
  byte ***img; // original image intensity
  byte ***tissue; // tissue = 0 is background, tissue = 1 is image data;
  byte ***edge;  // edge type, 0: no edge, 1: between image and background, 2: in-image edge
  byte ***grd; // gradient magnitude of original image;
  float ***geom; // [0,1] percentage of image volume within a neighborhood, 1-geom is percentage of background (CSF,VN,BG) within the neighborhood;
  short ****dxdydz; // gradient directions, now no use.
  byte ***flag; // have wav data.
  float ****vbox; // wav
} FAVBox_Small;


int fileExists(char *filename);
float * new_image3d(int x, int y, int z);
void free_image3d(float *img);
void read_image3d(float *img, int x, int y, int z, char *filename, char *format);
void write_image3d(float *img, int x, int y, int z, char *filename, char *format);
double uniform_noise();
double gaussian_noise(float mean, float sigma);
void no_mem_exit(char *where);
float ***Falloc3d(int X,int Y,int Z);
void Ffree3d(float ***array,int Z,int Y);
unshort ***USalloc3d(int X,int Y,int Z);
void USfree3d(unshort ***array,int Z,int Y);
short ***Salloc3d(int X,int Y,int Z);
void Sfree3d(short ***array,int Z,int Y);
byte ***Balloc3d(int X,int Y,int Z);
void Bfree3d(byte ***array,int Z,int Y);
float **Falloc2d(int X,int Y);
void Ffree2d(float **array,int Y);
byte **Balloc2d(int X,int Y);
void Bfree2d(byte **array,int Y);
void display_dbox(DBox *image);
void threshold_dbox(DBox *image, float threshold);
DBox * new_dbox(int X, int Y, int Z);
void half_intensity_bbox(BBox *image);
void free_dbox(DBox *tmp);
void init_dbox(DBox * dbox, int X, int Y, int Z);
DBox * subsampling_dbox(DBox *oimage, int power, float sigma);
DBox * generate_3d_gaussian_mask(float sigma, float xres, float yres, float zres);
void generate_3d_gaussian_image(DBox *mask, int px, int py, int pz, float sigma, float xres, float yres, float zres);
DBox * filtering_dbox_with_mask(DBox *oimg, DBox *mask);
float interpolate_dbox(DBox *db, float xx, float yy, float zz);
void clear_dbox(DBox *dbox);
void clear_bbox(BBox *dbox);
void smooth_dbox(DBox *dbox);
void iwt_3d_base_same(DBox * image, int x_l, int x_u, int y_l, int y_u, int z_l, int z_u);
void wt_3d_base_same(DBox * image, int x_l, int x_u, int y_l, int y_u, int z_l, int z_u);
Signal * new_signal(int Size);
void copy_signal(Signal *a, Signal *b);
void free_signal(Signal *signal);
void display_signal(Signal *sig);
void clear_signal(Signal *signal);
Deformation_Field *new_deformationfield(int X,int Y,int Z);
void free_deformationfield(Deformation_Field *df);
Points * read_points(char *filename);
Image * new_image(int X, int Y);
void free_image(Image *tmp);
float distance_3d_f(FPoint p1, FPoint p2);
FPoints * read_fpoints(char *filename);
void write_points(Points *pts, char *filename);
void write_fpoints(FPoints *pts, char *filename);
void display_points(Points *pts);
void display_fpoints(FPoints *pts);
Points * new_points(int N);
void free_points(Points *pts);
void copy_points(Points *pts1, Points *pts2);
void copy_fpoints(FPoints *pts1, FPoints *pts2);
FPoints * new_fpoints(int N);
void free_fpoints(FPoints *pts);
Points * generate_hood_byIncreasingRadius(int HoodSize, float xres,float yres,float zres) ;
float Gaussian_1d(float x, float sigma);
float Gaussian(int x, int y, int z, float sigma)     ;
float Gaussian_Similarity(int x, int y, int z, float sigma)     ;
void Smooth_DeformationField(Deformation_Field *df, float LocalRatio);
void Smooth_DeformationField_Without_BoundaryConstraints(Deformation_Field *df, float LocalRatio);
void Smooth_DeformationField_With_Edge_Preserving(Deformation_Field *df, float LocalRatio, FAVBox *mimg, FAVBox *simg);
DBox * volume_transform(DBox *tag, Deformation_Field *df);
void volume_transform_nointerpolation(DBox *tag, Deformation_Field *df, DBox *rst);
void volume_transform_interpolation(DBox *tag, Deformation_Field *df, DBox *rst);
DBox * volume_transform_forward(DBox *tag, Deformation_Field *df);
DBox * volume_transform_forward_bilge(DBox *tag, Deformation_Field *df);
DBox * volume_transform_forward_bilge_simple(DBox *tag, Deformation_Field *df, int do_segmentation);
void volume_transformation(DBox *org, DBox *tag, Deformation_Field *df);
void write_deformationfield(Deformation_Field *df, char *filename);
void write_deformationfield_with_switch(Deformation_Field *df, char *filename, int sw);
void read_deformationfield(Deformation_Field *df, char *filename);
void read_deformationfield_with_switch(Deformation_Field *df, char *filename, int sw);
BBox * generate_rotation_invariant_mask(int width);
BBox * generate_rotation_invariant_mask_shift(int width, int dx, int dy, int dz);
void copy_bbox(BBox *img1, BBox *img2);
BBox * new_bbox(int X, int Y, int Z);
void calloc_bbox(BBox *temp, int X, int Y, int Z);
void minux_bbox(BBox *bbox1, BBox*bbox2, BBox *rst);
long difference_bbox(BBox *bbox1, BBox*bbox2);
void free_bbox(BBox * tmp);
void free_bbox_data(BBox * tmp);
void display_bbox(BBox *tmp);
void mask_bbox(BBox *iimg,BBox *msk,BBox *oimg, int threshold);
void read_bbox(BBox *dbox, char *filename);
void write_bbox(BBox *dbox, char *filename);
float interpolate_bbox(BBox *db, float xx, float yy, float zz);
float * calculate_histogram_bbox(BBox *dbox);
void calculate_histogram_bbox_void(BBox *dbox, float *histo);
float * calculate_histogram(DBox *dbox);
float * calculate_histogram_mask(DBox *dbox, BBox *mask);
float * calculate_histogram_mask128(DBox *dbox, BBox *mask);
void calculate_histogram_void(DBox *dbox, float *histo);
float vec_dist2(float *a, float *b, int N);
DBox * calculate_3d_gradient_dd(DBox *dbox);
Image * histogram_gi(DBox *dbox, DBox *flag);
void write_matrix(Matrix *matrix, char *filename);
Matrix * read_matrix(char *filename);
void copy_matrix(Matrix *a, Matrix *b);
Transformation new_transformation(void);
Transformation read_transformation_fsl(char *filename);
void adjust_fsl_transform_to_myformat(Transformation *T, int XX, int YY, int ZZ);
void perform_global_image_transformation(BBox *iimg, BBox *oimg, Transformation T, int is_forward, int do_interpolation);
Transformation read_transformation(char *filename);
void write_transformation(Transformation T, char *filename);
void display_transformation(Transformation T);
FPoint subtract_fpoint(FPoint A, FPoint B);
FPoint normalize_fpoint(FPoint A);
FPoint scale_fpoint(FPoint A, float s);
void add_deformationfield(Deformation_Field *df, Deformation_Field *pdf, Deformation_Field *rst);
void subtract_deformationfield(Deformation_Field *df, Deformation_Field *pdf, Deformation_Field *rst);
void subtract_deformationfield_smallsize(Deformation_Field *df, Deformation_Field *pdf, Deformation_Field *rst, int BX, int EX, int BY, int EY, int BZ, int EZ);
void minus_deformationfield(Deformation_Field *odf, Deformation_Field *sdf, Deformation_Field *rst);
FPoint get_deformation_point(Deformation_Field *df, Point p);
FPoint get_deformation_fpoint(Deformation_Field *df, FPoint p);
void copy_deformationfield(Deformation_Field *dfa, Deformation_Field *dfb);
void clear_deformationfield(Deformation_Field *df);
void copy_dbox(DBox *dbx1, DBox *dbx2);
void XZPCA(Matrix *Samples, Matrix * Phi, Vector * lambda, Vector *AX);
void XZPCA_small_memory(Matrix *Samples, Matrix * Phi, Vector * lambda, Vector *AX);
void XZPCA1(Matrix *Samples, Matrix * Phi, Vector * lambda, Vector *AX);
int clip(Vector * lambda, double rate);
int return_small(int x, int y);
void PrincipalComponentAnalysis(Matrix *Samples, Matrix *Phi, Vector *lambda, Vector *AX);
void PrincipalComponentAnalysis_NoNewMemory(Matrix *Samples, Matrix *Phi, Vector *lambda, Vector *AX);
void free_vector(Vector *vec);
void display_vector(Vector *vec);
Vector * new_vector(int N);
Deformation_Field * inverse_deformationfield(Deformation_Field *df);
void get_a_smaller_size(DBox *img, int *oBX, int *oEX, int *oBY, int *oEY, int *oBZ, int *oEZ);
void get_a_smaller_size_bbox(BBox *img, int *oBX, int *oEX, int *oBY, int *oEY, int *oBZ, int *oEZ);
void get_a_smaller_size_bbox_samexysize(BBox *img, int *oBX, int *oEX, int *oBY, int *oEY, int *oBZ, int *oEZ);
float image_similarity(DBox *mdl, DBox *warped_sub);
Short_Deformation_Field * new_short_deformationfield(int X,int Y,int Z);
void free_short_deformationfield(Short_Deformation_Field * df);
void read_short_deformationfield_with_switch(Short_Deformation_Field *df, char *filename, char *format, int swch);
void read_short_deformationfield_with_switch_and_scale(Short_Deformation_Field *df, char *filename, char *format, int swch, float scale);
void Calculate_Jakobian(Deformation_Field *df, DBox *jcb);
Signal * fcm(float *histo, int Mx);
Signal * fcm_center_no(float *histo, int Mx, int C);
Signal * fcm_center_no_prior(float *histo, int Mx, int C, float *prior, float lambda);
void read_small_bbox(BBox *smbox, char *filename, int BX, int EX, int BY, int EY, int BZ, int EZ, int XX, int YY, int ZZ);
void read_small_dbox(DBox *smbox, char *filename, int BX, int EX, int BY, int EY, int BZ, int EZ, int XX, int YY, int ZZ);
void read_small_deformationfield(Deformation_Field *smdf, char *filename, int BX, int EX, int BY, int EY, int BZ, int EZ, int XX, int YY, int ZZ);
void xz_randomize();
void file_open_error(char *filename);
void Markov_Smoothing(Deformation_Field *tdf, Deformation_Field *tndf, float lambda, int Nsize);
void Markov_Smoothing_Prior(Deformation_Field *tdf, Deformation_Field *tndf, float lambda, int Nsize, DBox *sigma);
void Markov_Smoothing_Weights(Deformation_Field *tdf, Deformation_Field *tndf, float lambda, int Nsize, DBox *sigma);
void vectorfcm_center_no(float **data, int Mx, int My, int C, float **center);
FPoint fpoint_affine_transform(FPoint M, Matrix *A, FPoint *T);
DBox * calculate_brainboundary(DBox *img);
void set_I_matrix(Matrix *A);


#endif // _SBIA_DTIDROID_MATRIXSHENCLASSIC_H
