/**
 * @file  xzdgshenCLASSIC.h
 * @brief Image Processing Toolbox initally copied from CLASSIC.
 *
 * 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 _SBIA_DTIDROID_XZDGSHENCLASSIC_H
#define _SBIA_DTIDROID_XZDGSHENCLASSIC_H


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

#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;

double uniform_noise();
double gaussian_noise(float mean, float sigma);
void xz_randomize();
void file_open_error(char *filename);

/***************************

Part 1: Handling the memory

****************************/

int return_max(int x, int y);
int return_small(int x, int y);

void no_mem_exit(const char *where);

int get_mem2D(byte ***array2D, int rows, int columns);
void free_mem2D(byte **array2D, int rows);

int get_mem2Dint(int ***array2D, int rows, int columns);
void free_mem2Dint(int **array2D, int rows);

int get_mem2Dunshort(unshort ***array2D, int rows, int columns);
void free_mem2Dunshort(unshort **array2D, int rows);

int get_mem2Dfloat(float ***array2D, int rows, int columns);
void free_mem2Dfloat(float **array2D, int rows);

int get_mem3D(byte ****array3D, int frames, int rows, int columns);
void free_mem3D(byte ***array3D, int frames, int rows);

int get_mem3Dint(int ****array3D, int frames, int rows, int columns);
void free_mem3Dint(int ***array3D, int frames, int rows);

int get_mem3Dunshort(unshort ****array3D, int frames, int rows, int columns);
void free_mem3Dunshort(unshort ***array3D, int frames, int rows);

int get_mem3Dfloat(float ****array3D, int frames, int rows, int columns);
void free_mem3Dfloat(float ***array3D, int frames, int rows);

byte **Balloc2d(int X,int Y);
void Bfree2d(byte **array,int Y);

byte ***Balloc3d(int X,int Y,int Z);
void Bfree3d(byte ***array,int Z,int Y);

short ***Salloc3d(int X,int Y,int Z);
void Sfree3d(short ***array,int Z,int Y);

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

float **Falloc2d(int X,int Y);
void Ffree2d(float **array,int Y);


/*************************************

// XZVolume

// 3D volume image, unshort format.

**************************************/

XZVolume * new_volume(int, int, int);
void free_volume(XZVolume *);

void init_volume(XZVolume *, int, int, int);
void read_volume(XZVolume *, char *);
void write_volume(XZVolume *, char *);
XZVolume * subsampling_volume(XZVolume *oimage, int power);
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);

/************************************

DBox: Data Box

*************************************/

DBox * new_dbox(int, int, int);
void free_dbox(DBox *);
float dot_product_dbox(DBox *img1, DBox *img2);
void init_dbox(DBox *, int, int, int);
void get_dbox(XZVolume *, int, int, int, DBox *);
void clear_dbox(DBox *dbox);
void clear_bbox(BBox *dbox);
void read_dbox(DBox *dbox, const char *filename, const char *);
void write_dbox(DBox *dbox, const char *filename, const char *, int);
void display_dbox(DBox *image);
void set_dbox(DBox *image);
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);
DBox * subsampling_dbox(DBox *oimage, int power, float sigma);
DBox * filtering_dbox(DBox *oimage, float sigma);
DBox * filtering_dbox_with_mask(DBox *oimg, DBox *mask);
void filtering_dbox_2(DBox *oimage, float sigma, DBox *dimage);
void filtering_bbox_2(BBox *oimage, float sigma, DBox *dimage);

DBox * generatenewdbox(int xx, int yy, int zz);
void threshold_dbox(DBox *image, float threshold);

float interpolate_dbox(DBox *db, float jj, float ii, float kk);
void interpolate(DBox *low, int ratio, DBox *high);

/*************************************

XZVolume vs DBox

**************************************/

void volume2dbox(XZVolume *, DBox *);
void dbox2volume(DBox *, XZVolume *);
DBox *volume2newdbox(XZVolume *);
XZVolume *dbox2newvolume(DBox *, int);

/*************************************

 2D Image and 1D Signal

**************************************/

Image * new_image(int X, int Y);
void free_image(Image *);
void read_image(Image *image, char *filename, char *);
void write_image(Image *image, char *filename, char *);

Signal * new_signal(int Size);
void free_signal(Signal *);
void display_signal(Signal *sig);
void clear_signal(Signal *signal);
void copy_signal(Signal *, Signal *);

SSignal * new_ssignal(int Size);
void free_ssignal(SSignal *);

/***********************************************

VectorBox: 4D image.

 ***********************************************/

VectorBox * new_vectorbox(int X, int Y, int Z, int D);
void free_vectorbox(VectorBox *);
void display_vector(Vector *v);

AVBox * new_avbox(int X, int Y, int Z, int D, int N);
void free_avbox(AVBox *);
void write_avbox(AVBox *vb, char * filename);
AVBox * read_avbox(char * filename);
void get_img_from_avbox(AVBox *avbox, DBox *dbox);

void read_attribute_vector(VectorBox *vb, char *);
void write_attribute_vector(VectorBox *vb, char *);

void write_vectorbox(VectorBox *vb, char * filename);
void read_vectorbox(VectorBox *vb, char * filename);
VectorBox * read_new_vectorbox(char * filename);

int * get_vectorboxsize(char *filename);

void read_AV_from_vectorboxfile(char * filename, int XX, int YY, int ZZ, int DD, int xx, int yy, int zz, int ratio, byte *flag, float *vec);
//void read_AV_from_svectorboxfile(char * filename, int XX, int YY, int ZZ, int DD, int xx, int yy, int zz, short *vec);
void read_AV_from_fvectorboxfile(char * filename, int XX, int YY, int ZZ, int DD, int xx, int yy, int zz, byte *flag, float *vec);

void interpolate_vb(VectorBox *low, int ratio, VectorBox *high);
void interpolate_vbox(VectorBox *db, float xx, float yy, float zz, float *vector);

VectorBox * init_attribute_vector(int X, int Y, int Z, int level);
void set_attribute_vector(VectorBox *vb, DBox *wt, int level);
void normalize_attribute_vector(VectorBox *vb);
void normalize_attribute_vector_new(VectorBox *vb);
void normalize_attribute_vector_xz(VectorBox *vb);
//void normalize_AV(SVectorBox *vb);

//DBox * get_matching_map(VectorBox *vb, float *vector);
DBox * get_matching_map(VectorBox *vb, float *vector, float *co);
//DBox * get_matching_map_short(SVectorBox *vb, short *vector);
Image * get_matching_section(VectorBox *vb, char *axis, int xx, int yy, int zz, int skip_dc_flag);
Image * get_matching_section2(VectorBox *vb, char *axis, int xx, int yy, int zz, float *vector);
void search_matching_point(VectorBox *vb, float *vector, DBox *origimg);
void search_matching_point2(VectorBox *vb, float *vector, int xx, int yy, int zz, int r, DBox *origimg);
void mark_point(DBox *origimg, int xx, int yy, int zz);
float matching_degree_of_AVs(float *vector1, float *vector2, int D);

Image * get_section(DBox *vb, char *axis, int xx, int yy, int zz);

void set_AV(Signal *AV, DBox *wt, int level, int width, int dcflag);
void set_AV_new(Signal *AV, DBox *wt, int level);
int return_AV_size(int size, int level);
void set_AV2(Signal *AV, DBox *wt, int level, int width, int dcflag);
float norm_AV(Signal *av);
int return_AV_size_simple(int size, int level);
void set_AV_new_simple(Signal *AV, DBox *wt, int level);

/*************************************************************

FFT, IFFT

 *************************************************************/

complex Wn[512];

void FFT(complex *S,int lengh,char inv);
void FFT_3D(DBox *data_r, DBox *data_i, int inv);

/*************************************************************/
/*
/* Wavelet Transformation
/*
/*************************************************************/

// global variables

float *dwt_h, *dwt_g, *dwt_hi, *dwt_gi;
int dwt_filter_length;
int peak_pos_h;
int peak_pos_g;

// 1D/2D/3D DWT

Signal * wt_1d(int x_l, int x_u, Signal * signal);
void wt_1d_same(Signal * signal, int x_l, int x_u);
void wt_1d_same_no_downsampling(Signal * signal, Signal * high, int);

void wt_3d_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);
void wt_3d_same_no_downsampling(DBox * image, DBox *high, int);

Image * wt_2d(int x_l, int x_u, int y_l, int y_u, Image * image);
DBox * wt_3d(int x_l, int x_u, int y_l, int y_u, int z_l, int z_u, DBox * image);
DBox * wt_3d_volume(int x_l, int x_u, int y_l, int y_u, int z_l, int z_u, XZVolume * image);
void wt_3d_base(DBox * image, int, int, int);
void rotation_invariant_wt_3d_base(DBox * image, int dx, int dy, int dz);



//Signal * sliding_wt_3d(DBox * image, int xx, int yy, int zz, int width, int level);

void sliding_wt_3d(DBox * image, int xx, int yy, int zz, int width, int level, Signal *AV, byte *flag);
void sliding_wt_3d_new(DBox * image, int xx, int yy, int zz, int width, int level, Signal *AV, byte *flag);

void sliding_wt_3d_notsimple(DBox * image, int xx, int yy, int zz, int width, int level, Signal *AV, byte *flag);
void sliding_rotation_invariant(DBox * image, int xx, int yy, int zz, int width, SSignal *AV, BBox *mask);
BBox * generate_rotation_invariant_mask(int width);
BBox * generate_rotation_invariant_mask_shift(int width, int dx, int dy, int dz);

// 1D/2D/3D IDWT
void iwt_1d(Signal * wt, int x_l, int x_u, Signal * signal);

void iwt_1d_same(Signal *, int, int);

void iwt_2d(Image * wt, int x_l, int x_u, int y_l, int y_u, Image  * image);
void iwt_3d(DBox *, int x_l, int x_u, int y_l, int y_u, int z_l, int z_u, DBox  * image);

void iwt_3d_same(DBox *, int x_l, int x_u, int y_l, int y_u, int z_l, int z_u);
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);

// select the filter base to be used, input from 1,2,3,4,5,6,7,8:DD=6,8,10,12,14,16,18,20.

void select_filter(int);

/******************************************

Lump

 ******************************************/

void select_my_filter();
Signal * conv_1d(Signal *a, Signal *b);
Signal * xz_conv(Signal *m, Signal *s);
Signal * down_sampling(Signal *s);
Signal * up_sampling(Signal *s);
int wt_a(int x,int xsize);
void coef();

Vector_Series * new_vectorseries(int, int, int);
void write_vectorseries(Vector_Series *vs, char *filename);

Short_Signal * new_shortsignal(int Size);
void free_shortsignal(Short_Signal *);

DegreePoints * new_degreepoints(int);
void queque_degreepoints(DegreePoints *dpts);
void free_degreepoints(DegreePoints *);
void write_degreepoints(DegreePoints *pts, char *filename);
DegreePoints * read_degreepoints(char *filename);
void display_degreepoints(DegreePoints *pts);
void display_degreepoints_nozero(DegreePoints *pts);

Points * new_points(int N);
void copy_points(Points *pts1, Points *pts2);
void copy_fpoints(FPoints *pts1, FPoints *pts2);

Points * read_points(char *filename);
void write_points(Points *pts, char *filename);
FPoints * read_fpoints(char *filename);
int reorder_fpoints(FPoints *pts);
int reorder_fpoints_twosets(FPoints *pts, FPoints *spts);
void write_fpoints(FPoints *pts, char *filename);
void display_points(Points *pts);
void display_fpoints(FPoints *pts);
Points * subsampling_points(Points *opts, int scale);
Point upsampling_point(Point pt, int scale);
void free_points(Points *pts);

//Point get_matching_point(VectorBox *vb, float *vector);
//Point get_matching_point(VectorBox *vb, float *vector, int xx, int yy, int zz, int width);
Point get_matching_point(VectorBox *vb, VectorBox *vbo, Point P, Point C, Points *Hood, float * co, int searching_width);
Points * get_matching_points(VectorBox *vb, float *vector, int xx, int yy, int zz, int width, int number);

BBox * new_bbox(int, int, int);
void calloc_bbox(BBox *temp, int X, int Y, int Z);
void copy_bbox(BBox *img1, BBox *img2);
void display_bbox(BBox *tmp);
void free_bbox(BBox *);
void free_bbox_data(BBox * tmp);
void read_bbox(BBox *dbox, char *filename);
void mask_bbox(BBox *iimg,BBox *msk,BBox *oimg, int threshold);
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 write_bbox(BBox *dbox, char *filename);
float interpolate_bbox(BBox *db, float xx, float yy, float zz);
void half_intensity_bbox(BBox *image);
long difference_bbox(BBox *bbox1, BBox*bbox2);
void minux_bbox(BBox *bbox1, BBox*bbox2, BBox *rst);

//void MomentsOf3DImage(DBox *histo, int *center, int Range, Image *momentsMatrix3x3);

float Gaussian(int x, int y, int z, 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);
float Gaussian_Similarity(int x, int y, int z, float sigma);
float weight(Point p, float);
float weight_sum2one(Point p, float sigma);

void XZPCA(Matrix *Samples, Matrix * Phi, Vector * Lambda, Vector *AX);
void XZPCA_small_memory(Matrix *Samples, Matrix * Phi, Vector * lambda, Vector *AX);
int XZPCA_ondisk(char * samplebasefilename, int number_samples, int number_features, float percentage, char *PCAResultName);
int XZPCA_ondisk_splitfile(char * samplebasefilename, int number_samples, int number_features, float percentage, char *PCAResultName);
void UsePCA_ondisk(char *basefilename, Vector *bX, Vector *nX);
void UsePCA_ondisk_splitfile(char *basefilename, Vector *bX, Vector *nX, int number_of_pc, float);
void UsePCA_ondisk_calculatebX(char *basefilename, Vector *bX, Vector *nX);
void UsePCA_ondisk_calculatebX_splitfile(char *basefilename, Vector *bX, Vector *nX);
void XZPCA1(Matrix *Samples, Matrix * Phi, Vector * lambda, Vector *AX);
int clip(Vector * lambda, double rate);
//void free_vector(Vector *vec);
Vector * new_vector(int N);
void display_vector(Vector *vec);

//DegreePoints * get_matching_points_statistics(VectorBox *vb, Matrix *Phi, Vector *lambda, Vector *AX, Point cpt, Points *Hood, float, int width, int number);

//void XZPCA_FixedAX(Matrix *Samples, Matrix * Phi, Vector * lambda, Vector *AX);
void PrincipalComponentAnalisis_FixedAX(Matrix *Samples, Matrix *Phi, Vector *lambda, Vector *AX);

void PrincipalComponentAnalysis(Matrix *Samples, Matrix *Phi, Vector *lambda, Vector *AX);
void PrincipalComponentAnalysis_NoNewMemory(Matrix *Samples, Matrix *Phi, Vector *lambda, Vector *AX);

Points * generate_hood_byIncreasingRadius(int HoodSize, float xres,float yres,float zres) ;

PointArray3D *new_pointarray3d(int X,int Y,int Z,int BX,int BY,int BZ,int dx,int dy,int dz);
void free_pointarray3d(PointArray3D *pa);
void interpolate_pointarray(PointArray3D *pa);
DBox * get_warping_result(DBox *dbox1, PointArray3D *pa);
void init_pointarray3d(PointArray3D *pa, Points *mpts, Points *spts);

PointArray2D *new_pointarray2d(int X, int Y, int BX, int BY, int dx, int dy);
void free_pointarray2d(PointArray2D *pa);
void interpolate_pointarray2d(PointArray2D *pa);
Image * get_warping_result_2d(Image *img, PointArray2D *pa, PointArray2D *pad);
void init_pointarray2d(PointArray2D *pa, Points *pts, Points *grids);

//calculate the 3d gradient of a DBox

DBox * calculate_3d_gradient(DBox *dbox);
DBox * calculate_3d_gradient_bw(DBox *dbox);
DBox * calculate_3d_gradient_dd(DBox *dbox);

VectorBox * calculate_dxdydz(DBox *dbox);
Point get_avbox_matching_point(AVBox *vb, DBox *db, AVBox *vbo, Point P, Point C, Points *Hood, float *, float * co, int width, int *result_valide, int useedge);


Focus_Point * add_one_focus_point(Focus_Point *fp, Point p, float distance);
void display_all_focus_point(Focus_Point *fp);
float distance_3d(Point p1, Point p2);
float distance_3d_voxelsize(Point p1, Point p2, float vx, float vy, float vz);
void free_next(Focus_Point *focusp);
//int focus_point_distance_is_ok(Focus_Point *firstp, Point p, float distance);
int focus_point_distance_is_ok(Focus_Point *firstp, Point p, float distance, float *actual_dist);
int save_all_focus_point(Focus_Point *fp, char *filename);
void calcualte_minimal_distance(Focus_Point *currentp);
Focus_Point * read_all_focus_point(char *filename, int *pointnumber);
int transfer_focus_point_to_pts(Focus_Point *fp, Points *pts, int Number, float *minimal_distance);

float distance_3d_f(FPoint p1, FPoint p2);

FPoints * new_fpoints(int N);
void free_fpoints(FPoints *pts);

FPoint get_deformation_fpoint(Deformation_Field *df, FPoint p);
FPoint get_deformation_point(Deformation_Field *df, Point p);
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);
Point get_newmodelpoint(Deformation_Field *df, Point pt);
void display_deformationfield(Deformation_Field *df);
Deformation_Field *new_deformationfield(int X,int Y,int Z);
void free_deformationfield(Deformation_Field *df);
//void calculate_delta_deformationfield(Deformation_Field *delta_df, Points *originalpoints, Points *newpoints, float minimal_distance);

void update_deformationfield_float(Deformation_Field *df, Points *originalpoints, FPoints * newpoints, float distance_sigma, float rate, float lambda);
float update_deformationfield_combined(Deformation_Field *df, Points *originalpoints, FPoints * newpoints, FPoints * prevsubpts, FPoints * newsubpts, float distance_sigma, float rate, float lambda, float lambda1);

void read_deformationfield(Deformation_Field *df, char *filename);
void read_small_deformationfield(Deformation_Field *smfd, char *filename, int BX, int EX, int BY, int EY, int BZ, int EZ, int XX, int YY, int ZZ);
void read_deformationfield_with_switch(Deformation_Field *df, char *filename, int sw);
void write_deformationfield(Deformation_Field *df, char *filename);
void write_deformationfield_with_switch(Deformation_Field *df, char *filename, int sw);

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

// Focus Point Neighborhood


void update_deformationfield_float_allneighbor(Deformation_Field *df, Points *originalpoints, All_Neighbor_PTS all_neighbor_pts, FPoints * newpoints, float *, float distance_sigma, float rate, float lambda);
float update_deformationfield_combined_allneighbor(Deformation_Field *df, Points *originalpoints, All_Neighbor_PTS all_neighbor_pts, FPoints * newpoints, float *, FPoints * prevsubpts, FPoints * newsubpts, float distance_sigma, float rate, float lambda, float lambda1);

float update_deformationfield_allpoint_combined(Deformation_Field *df, Points *originalpoints, FPoints * newpoints, FPoints * prevsubpts, FPoints * newsubpts, float distance_sigma, float real_rate, float lambda, float lambda1, Deformation_Field *delta_df);

All_Neighbor_PTS new_all_neighbor_pts(Points *pts, int neighborhood_size);
void free_all_neighbor_pts(All_Neighbor_PTS tmp);
void display_all_neighbor_pts(All_Neighbor_PTS all_npts);

Point get_avbox_matching_point_newneighborhood(AVBox *vb, DBox *db, AVBox *vbo, Point P, Point C, Points *Hood, float *weights, float * co, int width, int *result_valide, int useedge, Deformation_Field *df, float *reliability);
Point get_avbox_matching_point_newneighborhood_fast(AVBox *vb, DBox *db, AVBox *vbo, Point P, Point C, Points *Hood, int hood_step, float *weights, float * co, int width, int grid_no, int useedge, Deformation_Field *df, float *reliability, float *similarity, Neighbor_PTS npts, Points * cur_sub_pts);

void affine_alignment(FPoints *modelpts, FPoints *subjectpts, Matrix *A, FPoint *T);

float degree_between_two_3D_vectors(FPoint one, FPoint two);

Image * get_matching_map_slice(VectorBox *, VectorBox *, Point , Points *Hood, float * co, int slice);

Point fast_matching(AVBox *vb, AVBox *vbo, Points *mod_pts, Points *sub_pts, int current_point, int width, Neighbor_PTS npts, float *co, Deformation_Field *df, float eta, float delta, float lambda_E, float lambda_D);

FAVBox * new_favbox(int X, int Y, int Z, int D);
void set_favbox_N(FAVBox *favbox, Signal *);
void free_favbox(FAVBox *);
void write_favbox(FAVBox *vb, char * filename);
FAVBox * read_favbox(char * filename);
void get_img_from_favbox(FAVBox *favbox, DBox *dbox);

void interpolate_favbox(FAVBox *low, int ratio, FAVBox *high);

void fast_matching_FAVBox(float *delta_E, FAVBox *vb, FAVBox *vbo, Points *mod_pts, Points *sub_pts, int width, All_Neighbor_PTS *all_npts, float *co, Deformation_Field *df, float eta, float delta, float lambda);
void fast_matching_FAVBox_new(float *delta_E, FAVBox *vb, FAVBox *vbo, Points *mod_pts, All_Neighbor_PTS *all_npts, float *co, FPoints *ddpp, float eta, float delta, float lambda, int hood_size);
void fast_matching_FAVBox_df(float *delta_E, FAVBox *vb, FAVBox *vbo, Points *mod_pts, All_Neighbor_PTS *all_npts, float *co, Deformation_Field *df, float eta, float delta, float lambda, int hood_size);
float AV_similarity_FAVBox_inter(FAVBox *vbo, Point p, float *co, FPoint np, FAVBox *vb, int *valide);
float AV_similarity_FAVBox(FAVBox *vbo, Point p, FAVBox *vb, Point np, float *co);

//Point get_avbox_matching_point_newneighborhood_fast_favbox(FAVBox *vb, DBox *db, FAVBox *vbo, Point P, Point C, Points *Hood, int hood_step, float *weights, float * co, int width, int grid_step, int useedge, Deformation_Field *df, float *reliability, float *similarity, Neighbor_PTS npts, Points * cur_sub_pts);

Point get_avbox_matching_point_newneighborhood_fast_favbox(FAVBox *vb, DBox *db, FAVBox *vbo, Point P, Point C, Points *Hood, int hood_step, float *weights, float * co, int width, int grid_step, int useedge, Deformation_Field *df, float *reliability, float *similarity, Neighbor_PTS npts, Points * cur_sub_pts, float coeff_distance, float coeff_similarity);

Point get_avbox_matching_point_newneighborhood_fast_favbox_fast(FAVBox *vb, DBox *db, FAVBox *vbo, Point P, Point C, Points *Hood, int hood_step, float *weights, float * co, int width, int grid_step, int useedge, Deformation_Field *df, float *reliability, float *similarity, Neighbor_PTS npts, Points * cur_sub_pts, float coeff_distance, float coeff_similarity);

Point get_avbox_matching_point_newneighborhood_fast_favbox_fast_with_costfunction(FAVBox *vb, DBox *db, FAVBox *vbo, Point P, Point C, Points *Hood, int hood_step, float *weights, float * co, int width, int grid_step, int useedge, Deformation_Field *df, float *reliability, float *similarity, Neighbor_PTS npts, Points * cur_sub_pts, float coeff_distance, float coeff_similarity);

void smooth_dbox(DBox *dbox);
void strong_smooth_dbox_mask(DBox *dbox, BBox *mask);
void adaptive_smooth_dbox(DBox *dbox, DBox *weighting);
void update_img_in_favbox(FAVBox *favbox, DBox *dbox);
void update_grd_in_favbox(FAVBox *favbox, DBox *dbox);
void update_edge_in_favbox(FAVBox *favbox, DBox *dbox);
void update_edge_tissue_in_favbox(FAVBox *favbox, DBox *edge, DBox *tissue);
void update_edge_tissue_in_favbox_separately(FAVBox *favbox, DBox *edge, DBox *tissue);
void update_dxdydzin_favbox(FAVBox *favbox, VectorBox *vbox);

Point get_favbox_matching_point(FAVBox *vb, DBox *db, FAVBox *vbo, Point P, Point C, Points *Hood, float *weights, float * co, int width, int *result_valide, int useedge, float *similarity);
void calculate_cost_function(FAVBox *vbmod, FAVBox *vbsub, Deformation_Field *df, float *co);
float calculate_cost_function_new(FAVBox *vb, FAVBox *vbo, Points *mod_pts, float *co, FPoints *ddpp,  All_Neighbor_PTS *all_npts, float lambda, int hood_size);
float calculate_cost_function_df(FAVBox *vb, FAVBox *vbo, Points *mod_pts, float *co, Deformation_Field *df, FPoints *ddpp,  All_Neighbor_PTS *all_npts, float lambda, int hood_size);
//DBox * get_matching_map_favbox(FAVBox *vb, float *vector, float *co);
DBox * get_matching_map_favbox(FAVBox *vb, FAVBox *vbo, Point P, float *co, Points *Hood, float *weights);
Image * get_matching_map_favbox_z(FAVBox *vb, FAVBox *vbo, Point P, Signal *co, Points *Hood, float *weights, int Z);
Vector * histogram_xz(DBox *dbox, DBox *flag);
Image * histogram_gi(DBox *dbox, DBox *flag);
void normalize_dbox_histogram(DBox *dbox, float *histo);
void normalize_dbox_mv(DBox *dbox, float *mean, float *variance);

void update_deformationfield_normalized(FAVBox *modav, Deformation_Field *df, Points *modelpoints, Points *newpoints, All_Neighbor_PTS *all_npts, float sigma);

void recalculate_favbox_flag(FAVBox *vb);
void recalculate_favbox_flag_D0(FAVBox *vb);

Point get_favbox_matching_point_use_edgeimage_x(FAVBox *vb, DBox *sedge, FAVBox *vbo, DBox *oedge, Point P, Point C, Points *Hood,float *weights, int width, Deformation_Field *df, float *reliability, float *similarity, Neighbor_PTS npts, Points * cur_sub_pts, float coeff_distance, float coeff_similarity, int, float);
Point get_favbox_matching_point_use_edgeimage_y(FAVBox *vb, DBox *sedge, FAVBox *vbo, DBox *oedge, Point P, Point C, Points *Hood,float *weights, int width, Deformation_Field *df, float *reliability, float *similarity, Neighbor_PTS npts, Points * cur_sub_pts, float coeff_distance, float coeff_similarity, int, float);
Point get_favbox_matching_point_use_edgeimage_z(FAVBox *vb, DBox *sedge, FAVBox *vbo, DBox *oedge, Point P, Point C, Points *Hood,float *weights, int width, Deformation_Field *df, float *reliability, float *similarity, Neighbor_PTS npts, Points * cur_sub_pts, float coeff_distance, float coeff_similarity, int, float);

Point get_favbox_matching_point_use_edgeimage(FAVBox *vb, DBox *sedge, DBox *bsimg, FAVBox *vbo, DBox *boimg, Point P, Point C, Points *Hood,float *weights, int width, Deformation_Field *df, float *reliability, float *similarity, Neighbor_PTS npts, Points * cur_sub_pts, float coeff_distance, float coeff_similarity);

int * determine_direction_2d(DBox *mimg, Points *mpts, char ch);
int search_here(int x0, int y0, int direction, int x, int y, float threshold);
int determine_direction(DBox *mimg, Point mpts, char ch);

void read_one_df_vector_with_switch(FPoint *FP, int position, FILE *fp, int sw);


Matrix * read_matrix(char *filename);
void write_matrix(Matrix *, char *);
void copy_matrix(Matrix *, Matrix *);

Transformation new_transformation(void);
Transformation read_transformation(char *filename);
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 *, BBox *, Transformation T, int, int);
void write_transformation(Transformation T, char *filename);
//void write_transformation(Transformation, char * filename);
void display_transformation(Transformation T);

void sliding_pdf(DBox * image, VectorBox *grad, int xx, int yy, int zz, int width_radium, int delta_width, int level, Signal *AV, int *valide);

VectorBox * calculate_3d_gradient_vector(DBox *dbox);

float Gaussian_1d(float x, float sigma);

Deformation_Field * sphere_shrunk(float R, float r, float b);
Deformation_Field * sphere_shrunk_directional(FPoint fpt, float R, float r, float b);
Deformation_Field * sphere_expand(float R, float r, float b);
FPoint subtract_fpoint(FPoint A, FPoint B);
FPoint normalize_fpoint(FPoint A);
FPoint scale_fpoint(FPoint A, float s);
float calculate_dt(float R, float r, float t);
float calculate_inverse_dt(float R, float r, float t);
void combine_deformationfield(Deformation_Field *df, Deformation_Field *sdf, Point CP);
void copy_dbox(DBox *dbx1, DBox *dbx2);
void clear_deformationfield(Deformation_Field *df);
void copy_deformationfield(Deformation_Field *df, Deformation_Field *);


Signal * get_AV_using_wtlh(DBox * wt_l, DBox * wt_h,  BBox *r_mask, int no_of_features);
void get_wt_low_high_featuresubimages(DBox *wt, int S, DBox *l, DBox *h);

void sliding_3ddwt(DBox * image, int xx, int yy, int zz, int width, int level, Signal *AV, byte *flag, BBox **mask, Signal *AVS);

Signal * get_AV_size(int width, int level, int);

float av_similarity(FAVBox *vb1, Point p1, FAVBox *vb2, Point p2, Signal *co);
float av_similarity_switch(FAVBox *vb1, Point p1, FAVBox *vb2, Point p2, Signal *co, int switchc);
float av_similarity_showmap(FAVBox *vb1, Point p1, FAVBox *vb2, Point p2, Signal *co);
float vector_similarity(float *v1, float *v2, int bgn, int edn);
void calculating_hoodfield(Point P, Points *hood, Deformation_Field *df, FPoints *hoodfield);

Point get_corresponding_point_favbox(Deformation_Field *df, FAVBox *vb, FAVBox *vbo, Point P, Point C, Points *Hood, float *weights, Signal * cosig, int width, int *result_valide, float *out_similarity);

DegreePoints * calculate_similarity_hari(Deformation_Field *df, FAVBox *vbo, FAVBox *vb, Points *cur_mod_pts, int cur_mod_idx, Points *Hood, float *weights, Signal * cosig, Points *search_hood);

Point get_corresponding_point_favbox_fast(FAVBox *vb, FAVBox *vbo, Point P, Point C, Points *Hood, FPoints *Hood_Field, float *weights, Signal * cosig, int width, int *result_valide, float *out_similarity);

Point get_corresponding_point_heart(Deformation_Field *df, FAVBox *vbo, FAVBox *vb, Points *cur_mod_pts, int cur_mod_idx, Points *Hood, float *weights, Signal * cosig, Points *search_hood, int *result_valide, float *out_similarity, float, int);

void get_corresponding_points_from_subject_to_model(Deformation_Field *df, Deformation_Field *idf, FAVBox *vbm, FAVBox *vbs, Points *cur_mod_pts, Points *all_sub_pts, Points *Hood, float hood_width_sigma, Signal * cosig, Points *search_hood, float search, FPoints *mdf, float SMR_THR, int swc);

void update_deformationfield(Deformation_Field *df, Points *cur_mod_pts, int cur_mod_idx, Point NC, float real_rate);

void update_deformationfield_modelptsdfm(Deformation_Field *df, Points *cur_mod_pts, FPoints *dfm, float real_rate, Deformation_Field *);

void constrain_deformationfield(Deformation_Field *df, Deformation_Field *pdf, Points *cur_mod_pts, Points *prev_mod_pts,  float lambda, float lambda1);

Point get_corresponding_point(Deformation_Field *df, FAVBox *vbo, FAVBox *vb, Points *cur_mod_pts, int cur_mod_idx,Points *Hood, float *weights, Signal * cosig, Points *search_hood, int *result_valide, float *out_similarity);

void normalize_favbox(FAVBox *wav);
float norm_similarity(float ax, float vx, float ay, float vy);
Signal * get_AV_moment_using_wtlh(DBox * wt_l, DBox * wt_h, BBox *r_mask);
void calculate_3_moments(BBox *r_mask, DBox *data, float *I0, float *I1, float *I2, float *I3);
void sliding_3ddwt_moments(DBox * image, int xx, int yy, int zz, int width, int level, Signal *AV, byte *flag, BBox **mask, Signal *AVS);
void normalize_favbox_moment(FAVBox *wav);
float av_similarity_moment(FAVBox *vb1, Point p1, FAVBox *vb2, Point p2, Signal *co);
Point get_corresponding_point_moment(Deformation_Field *df, FAVBox *vbo, FAVBox *vb, Points *cur_mod_pts, int cur_mod_idx,Points *Hood, float *weights, Signal * cosig, Points *search_hood, int *result_valide, float *out_similarity);

Deformation_Field * inverse_deformationfield(Deformation_Field *df);


void sliding_3d_moments(DBox * image, int xx, int yy, int zz, int width, Signal *AV, byte *flag, BBox *mask);
Image * get_matching_map_favbox_z_moment(FAVBox *vb, FAVBox *vbo, Point P, Signal *co, Points *Hood, float *weights, int Z);
Signal * get_AV_moment(DBox * temp, BBox *mask);

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 * calculate_histogram(DBox *dbox);
float * calculate_histogram_mask(DBox *dbox, BBox *mask);
float * calculate_histogram_mask128(DBox *dbox, BBox *mask);
float * calculate_histogram_bbox(BBox *dbox);
void calculate_histogram_bbox_void(BBox *dbox, float *histo);
BBox * calculate_3d_gradient_bbox(BBox *dbox);

void calculate_histogram_void(DBox *dbox, float *histo);

float image_similarity(DBox *mdl, DBox *warped_sub);

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

Point find_maxposition_dbox(DBox *db, float *maxvalue);



Image_4D * new_4dimage(int, int, int, int);
void read_4dimage(Image_4D *, char *);
void write_4dimage(Image_4D *, char *);
void free_4dimage(Image_4D *);
DBox * get_new_3dimage_from_4dimage(Image_4D *img4d, int t);
BBox * get_new_bbox_from_4dimage(Image_4D *img4d, int t);
void get_3dimage_from_4dimage(Image_4D *img4d, int t, DBox *img);
float * calculate_histogram_4dimage(Image_4D *img);


Points_4D * new_4dpoints(int N);
void free_4dpoints(Points_4D *tmp);


FAVBox * read_smallfavbox(char * filename);
void free_smallfavbox(FAVBox *tmp);

Short_Deformation_Field * new_short_deformationfield(int,int,int);
void free_short_deformationfield(Short_Deformation_Field *);
void read_short_deformationfield_with_switch(Short_Deformation_Field *df, char *filename, char *format, int swch);

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

DBox * calculate_jacobian_C(Deformation_Field *df);

void FFT(complex *S, int lengh, char inv);
int fileExists(char *filename);

typedef struct 
{
  int height;
  int width;
  int N;
  int *x;
  int *y;
  float *value;
} SparseMatrix;

typedef struct
{
  int Num;
  int * N;
  int ** index;
} sMCol;

SparseMatrix * determine_matrix_pts(Points *pts, int XX, int YY, int ZZ);
SparseMatrix * determine_matrix(int XX, int YY, int ZZ);
SparseMatrix * combine_smatrix(SparseMatrix *sm);
void free_smatrix(SparseMatrix *sM);
void display_smatrix(SparseMatrix *sMatrix);
SparseMatrix * smatrix_time_smatrix(SparseMatrix *sm1, SparseMatrix *sm2, sMCol *smc2);
void write_smatrix(SparseMatrix *sm, char * filename, int XX, int YY, int ZZ);
sMCol * make_index(SparseMatrix * sm);
sMCol * make_row_index(SparseMatrix * sm);
void display_sMCol(sMCol *smc, SparseMatrix *sm);
int is_it_ok(int *row, int *col);
int is_it_ok_cal(int *row, int *col, int *NN, int *common1, int *common2);
SparseMatrix * read_smatrix(char *filename, int *XX, int *YY, int *ZZ);
Matrix * smatrix_times_matrix(SparseMatrix *sm, Matrix *m);
Matrix * smatrix_times_matrix_scale3(SparseMatrix *sm, Matrix *m);
void return_xyz(int index, int X, int Y, int Z, int *x, int *y, int *z);
void free_sMCol(sMCol *smc);



float vec_dist(float *a, float *b, int N);
float vec_dist2(float *a, float *b, int N);
void vectorfcm_center_no(float **data, int Mx, int My, int C, float **center);

void determine_affine_transformation(DBox *imgm, DBox *imgs, Matrix *A, FPoint *T);
FPoint fpoint_affine_transform(FPoint M, Matrix *A, FPoint *T);
DBox * calculate_3d_gradient_eachaxis(DBox *dbox, int i);
void set_I_matrix(Matrix *A);
DBox * affine_transform_inverse(DBox *simg, Matrix *A, FPoint T, int);
void determine_affine_transformation_intensity_plus_gradient(DBox *imgm, DBox *imgs, Matrix *A, FPoint *T);
DBox * calculate_brainboundary(DBox *img);
void determine_affine_transformation_intensity_plus_brainboundary(DBox *imgm, DBox *imgs, Matrix *A, FPoint *T, float epsilon, float w_int, float w_grd, int dof);


/***
image3d
 ***/

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

BBox * Generate_Eppipsoid_Neighborhood(float a, float b, float c, float alpha, float beta, float gama);
void Generate_Eppipsoid_Neighborhood_GiveMask(BBox *mask, float a, float b, float c, float alpha, float beta, float gama);
float T_Value(float *va, int na, float *vb, int nb);
float Calculate_TValues_Of_Images(DBox **imga, int Na, DBox **imgb, int Nb, BBox *mask, int px, int py, int pz);
float Determine_EllipsoidParameters_Using_GradientDescent(DBox **imgsa, int Na, DBox **imgsb, int Nb, int px, int py, int pz, float inita, float initb, float initc, float initalpha, float initbeta, float initgama, float * outa, float * outb, float * outc, float * outalpha, float * outbeta, float * outgama);
BBox ** PreCalculationEppipsoid(int N_gama, int N_beta, float ini_r);
float optimize_filter_tscore(float radius, DBox **imgsa, int Na, DBox **imgsb, int Nb, int px, int py, int pz);
float optimize_filter_tscore_returnmask(float radius, DBox **imgsa, int Na, DBox **imgsb, int Nb, int px, int py, int pz, BBox *mask);


#endif // _SBIA_DTIDROID_XZDGSHENCLASSIC_H
