From 517876d61f3637a07b5e36246be172cc1608ba0f Mon Sep 17 00:00:00 2001
From: heziqi <heziqi36@gmail.com>
Date: Tue, 21 Oct 2014 04:25:47 +0800
Subject: [PATCH] metrics finished and memory allocation error fixed

---
 envi/bil.h  | 122 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----
 envi/bip.h  | 131 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------
 envi/bsq.h  | 122 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---
 envi/envi.h |  36 ++++++++++++++++++++++++++++++++++++
 4 files changed, 398 insertions(+), 13 deletions(-)

diff --git a/envi/bil.h b/envi/bil.h
index 733c7f5..c1ae161 100644
--- a/envi/bil.h
+++ b/envi/bil.h
@@ -60,8 +60,7 @@ public:
 		unsigned int S = XY * sizeof(T);		//calculate the number of bytes of a band
 
 		unsigned page=0;                      //bands around the wavelength
-		T * p1;
-		T * p2;
+		
 
 		//get the bands numbers around the wavelength
 
@@ -82,6 +81,8 @@ public:
 		}
 		if ( wavelength < w[page] )
 		{
+			T * p1;
+			T * p2;
 			p1=(T*)malloc(S);                     //memory allocation
 			p2=(T*)malloc(S);
 			band_index(p1, page - 1);
@@ -90,6 +91,8 @@ public:
 				double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);
 				p[i] = (p2[i] - p1[i]) * r + p1[i];
 			}
+			free(p1);
+			free(p2);
 		}
 		else                           //if the wavelength is equal to a wavelength in header file
 		{
@@ -401,8 +404,6 @@ public:
 		}
 		return true;
 	}
-	//untested
-	//providing the left and the right bound wavelength, return baseline-corrected band height
 	bool height(double lb, double rb, double bandwavelength, T* result){
 
 		T* lp;
@@ -417,6 +418,8 @@ public:
 
 		baseline_band(lb, rb, lp, rp, bandwavelength, result);
 
+		free(lp);
+		free(rp);
 		return true;
 	}
 
@@ -479,6 +482,7 @@ public:
 		for(unsigned j = 0; j < XY; j++){	
 			result[j] += (w[ai] - lab) * (cur[j] + cur2[j]) / 2.0;
 		}
+
 		//calculate the area
 		ai++;
 		for(unsigned i = ai; i <= bi ;i++)
@@ -491,6 +495,10 @@ public:
 			std::swap(cur,cur2);		//swap the band pointers
 		}
 
+		free(lp);
+		free(rp);
+		free(cur);
+		free(cur2);
 		return true;
 	}
 
@@ -510,6 +518,9 @@ public:
 			else
 				result[i] = p1[i] / p2[i];
 		}
+
+		free(p1);
+		free(p2);
 		return true;	
 	}
 	
@@ -530,6 +541,9 @@ public:
 			else
 				result[i] = p1[i] / p2[i];
 		}
+
+		free(p1);
+		free(p2);
 		return true;	
 	}		
 	
@@ -550,9 +564,108 @@ public:
 			else
 				result[i] = p1[i] / p2[i];
 		}
+
+		free(p1);
+		free(p2);
 		return true;	
 	}		
 
+	//x * f(x)
+	bool x_area(double lb, double rb, double lab, double rab, T* result){
+		T* lp;	//left band pointer
+		T* rp;	//right band pointer
+		T* cur;		//current band 1
+		T* cur2;	//current band 2
+
+		unsigned XY = R[0] * R[1];
+		unsigned S = XY * sizeof(T);
+
+		lp = (T*) malloc(S);			//memory allocation
+		rp = (T*) malloc(S);
+		cur = (T*) malloc(S);
+		cur2 = (T*) malloc(S);
+
+		memset(result, (char)0, S);
+
+		//find the wavelenght position in the whole band
+		unsigned int n = w.size();
+		unsigned int ai = 0;		//left bound position
+		unsigned int bi = n - 1;		//right bound position
+
+		//to make sure the left and the right bound are in the bandwidth
+		if (lb < w[0] || rb < w[0] || lb > w[n-1] || rb >w[n-1]){
+			std::cout<<"ERROR: left bound or right bound out of bandwidth"<<std::endl;
+			exit(1);
+		}
+		//to make sure rigth bound is bigger than left bound
+		else if(lb > rb){
+			std::cout<<"ERROR: right bound should be bigger than left bound"<<std::endl;
+			exit(1);
+		}
+
+		//get the position of lb and rb
+		while (lab >= w[ai]){
+			ai++;
+		}
+		while (rab <= w[bi]){
+			bi--;
+		}
+
+		band(lp, lb);
+		band(rp, rb);
+
+		//calculate the beginning and the ending part
+		baseline_band(lb, rb, lp, rp, rab, cur2);		//ending part
+		baseline_band(lb, rb, lp, rp, w[bi], cur);
+		for(unsigned j = 0; j < XY; j++){
+			result[j] += (rab - w[bi]) * (rab + w[bi]) * (cur[j] + cur2[j]) / 4.0;
+		}
+		baseline_band(lb, rb, lp, rp, lab, cur2);		//beginnning part
+		baseline_band(lb, rb, lp, rp, w[ai], cur);
+		for(unsigned j = 0; j < XY; j++){	
+			result[j] += (w[ai] - lab) * (w[ai] + lab) * (cur[j] + cur2[j]) / 4.0;
+		}
+
+		//calculate f(x) times x
+		ai++;
+		for(unsigned i = ai; i <= bi ;i++)
+		{
+			baseline_band(lb, rb, lp, rp, w[ai], cur2);
+			for(unsigned j = 0; j < XY; j++)
+			{
+				result[j] += (w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * (cur[j] + cur2[j]) / 4.0; 
+			}
+			std::swap(cur,cur2);		//swap the band pointers
+		}
+
+		free(lp);
+		free(rp);
+		free(cur);
+		free(cur2);
+		return true;
+	}
+
+	//centroid point
+	bool cpoint(double lb, double rb, double lab, double rab, T* result){
+		T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T));	
+		T* p2 = (T*)malloc(R[0] * R[1] * sizeof(T));	
+		
+		//get the area and the peak band
+		x_area(lb, rb, lab, rab, p1);
+		area(lb, rb, lab, rab, p2);
+		//calculate the ratio in result
+		for(unsigned i = 0; i < R[0] * R[1]; i++){
+			if(p1[i] == 0 && p2[i] ==0)
+				result[i] = 1;
+			else
+				result[i] = p1[i] / p2[i];
+		}
+
+		free(p1);
+		free(p2);
+		return true;	
+	}
+
 	//create mask file
 	bool mask(unsigned char* p, double mask_band, double threshold){
 
@@ -566,6 +679,7 @@ public:
 					p[i] = 255;
 		}
 
+		free(temp);
 		return true;
 	}
 
diff --git a/envi/bip.h b/envi/bip.h
index 4198751..06b59f0 100644
--- a/envi/bip.h
+++ b/envi/bip.h
@@ -63,8 +63,7 @@ public:
 		unsigned int XY = R[0] * R[1];	//calculate the number of pixels in a band
 
 		unsigned page=0;                      //bands around the wavelength
-		T * p1;
-		T * p2;
+
 
 		//get the bands numbers around the wavelength
 
@@ -85,6 +84,8 @@ public:
 		}
 		if ( wavelength < w[page] )
 		{
+			T * p1;
+			T * p2;
 			p1=(T*)malloc( XY * sizeof(T));                     //memory allocation
 			p2=(T*)malloc( XY * sizeof(T));
 			band_index(p1, page - 1);
@@ -93,6 +94,8 @@ public:
 				double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);
 				p[i] = (p2[i] - p1[i]) * r + p1[i];
 			}
+			free(p1);
+			free(p2);
 		}
 		else                           //if the wavelength is equal to a wavelength in header file
 		{
@@ -108,8 +111,7 @@ public:
 		unsigned int B = R[2];
 
 		unsigned page=0;                      //samples around the wavelength
-		T * p1;
-		T * p2;
+
 
 		//get the bands numbers around the wavelength
 
@@ -136,6 +138,8 @@ public:
 		}
 		if ( wavelength < w[page] )
 		{
+			T * p1;
+			T * p2;
 			p1=(T*)malloc( X * sizeof(T));                     //memory allocation
 			p2=(T*)malloc( X * sizeof(T));
 			//band_index(p1, page - 1);
@@ -153,6 +157,8 @@ public:
 				double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);
 				p[i] = (p2[i] - p1[i]) * r + p1[i];
 			}
+			free(p1);
+			free(p2);
 		}
 		else                           //if the wavelength is equal to a wavelength in header file
 		{
@@ -488,8 +494,7 @@ public:
 		}
 		return true;
 	}
-	//untested
-	//providing the left and the right bound wavelength, return baseline-corrected band height
+	
 	bool height(double lb, double rb, double bandwavelength, T* result){
 
 		T* lp;
@@ -504,6 +509,8 @@ public:
 
 		baseline_band(lb, rb, lp, rp, bandwavelength, result);
 
+		free(lp);
+		free(rp);
 		return true;
 	}
 
@@ -566,6 +573,7 @@ public:
 		for(unsigned j = 0; j < XY; j++){	
 			result[j] += (w[ai] - lab) * (cur[j] + cur2[j]) / 2.0;
 		}
+
 		//calculate the area
 		ai++;
 		for(unsigned i = ai; i <= bi ;i++)
@@ -578,6 +586,10 @@ public:
 			std::swap(cur,cur2);		//swap the band pointers
 		}
 
+		free(lp);
+		free(rp);
+		free(cur);
+		free(cur2);
 		return true;
 	}
 
@@ -597,6 +609,9 @@ public:
 			else
 				result[i] = p1[i] / p2[i];
 		}
+
+		free(p1);
+		free(p2);
 		return true;	
 	}
 	
@@ -617,6 +632,9 @@ public:
 			else
 				result[i] = p1[i] / p2[i];
 		}
+
+		free(p1);
+		free(p2);
 		return true;	
 	}		
 	
@@ -637,8 +655,108 @@ public:
 			else
 				result[i] = p1[i] / p2[i];
 		}
+
+		free(p1);
+		free(p2);
 		return true;	
 	}		
+
+	//x * f(x)
+	bool x_area(double lb, double rb, double lab, double rab, T* result){
+		T* lp;	//left band pointer
+		T* rp;	//right band pointer
+		T* cur;		//current band 1
+		T* cur2;	//current band 2
+
+		unsigned XY = R[0] * R[1];
+		unsigned S = XY * sizeof(T);
+
+		lp = (T*) malloc(S);			//memory allocation
+		rp = (T*) malloc(S);
+		cur = (T*) malloc(S);
+		cur2 = (T*) malloc(S);
+
+		memset(result, (char)0, S);
+
+		//find the wavelenght position in the whole band
+		unsigned int n = w.size();
+		unsigned int ai = 0;		//left bound position
+		unsigned int bi = n - 1;		//right bound position
+
+		//to make sure the left and the right bound are in the bandwidth
+		if (lb < w[0] || rb < w[0] || lb > w[n-1] || rb >w[n-1]){
+			std::cout<<"ERROR: left bound or right bound out of bandwidth"<<std::endl;
+			exit(1);
+		}
+		//to make sure rigth bound is bigger than left bound
+		else if(lb > rb){
+			std::cout<<"ERROR: right bound should be bigger than left bound"<<std::endl;
+			exit(1);
+		}
+
+		//get the position of lb and rb
+		while (lab >= w[ai]){
+			ai++;
+		}
+		while (rab <= w[bi]){
+			bi--;
+		}
+
+		band(lp, lb);
+		band(rp, rb);
+
+		//calculate the beginning and the ending part
+		baseline_band(lb, rb, lp, rp, rab, cur2);		//ending part
+		baseline_band(lb, rb, lp, rp, w[bi], cur);
+		for(unsigned j = 0; j < XY; j++){
+			result[j] += (rab - w[bi]) * (rab + w[bi]) * (cur[j] + cur2[j]) / 4.0;
+		}
+		baseline_band(lb, rb, lp, rp, lab, cur2);		//beginnning part
+		baseline_band(lb, rb, lp, rp, w[ai], cur);
+		for(unsigned j = 0; j < XY; j++){	
+			result[j] += (w[ai] - lab) * (w[ai] + lab) * (cur[j] + cur2[j]) / 4.0;
+		}
+
+		//calculate f(x) times x
+		ai++;
+		for(unsigned i = ai; i <= bi ;i++)
+		{
+			baseline_band(lb, rb, lp, rp, w[ai], cur2);
+			for(unsigned j = 0; j < XY; j++)
+			{
+				result[j] += (w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * (cur[j] + cur2[j]) / 4.0; 
+			}
+			std::swap(cur,cur2);		//swap the band pointers
+		}
+
+		free(lp);
+		free(rp);
+		free(cur);
+		free(cur2);
+		return true;
+	}
+
+	//centroid point
+	bool cpoint(double lb, double rb, double lab, double rab, T* result){
+		T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T));	
+		T* p2 = (T*)malloc(R[0] * R[1] * sizeof(T));	
+		
+		//get the area and the peak band
+		x_area(lb, rb, lab, rab, p1);
+		area(lb, rb, lab, rab, p2);
+		//calculate the ratio in result
+		for(unsigned i = 0; i < R[0] * R[1]; i++){
+			if(p1[i] == 0 && p2[i] ==0)
+				result[i] = 1;
+			else
+				result[i] = p1[i] / p2[i];
+		}
+
+		free(p1);
+		free(p2);
+		return true;	
+	}
+
 	//create mask file
 	bool mask(unsigned char* p, double mask_band, double threshold){
 
@@ -652,6 +770,7 @@ public:
 					p[i] = 255;
 		}
 
+		free(temp);
 		return true;
 
 	}
diff --git a/envi/bsq.h b/envi/bsq.h
index 54bac04..089b420 100644
--- a/envi/bsq.h
+++ b/envi/bsq.h
@@ -64,8 +64,7 @@ public:
 		unsigned int XY = R[0] * R[1];	//calculate the number of pixels in a band
 
 		unsigned page=0;                      //bands around the wavelength
-		T * p1;
-		T * p2;
+
 
 		//get the bands numbers around the wavelength
 
@@ -86,6 +85,8 @@ public:
 		}
 		if ( wavelength < w[page] )
 		{
+			T * p1;
+			T * p2;
 			p1=(T*)malloc( XY * sizeof(T));                     //memory allocation
 			p2=(T*)malloc( XY * sizeof(T));
 			band_index(p1, page - 1);
@@ -94,6 +95,8 @@ public:
 				double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);
 				p[i] = (p2[i] - p1[i]) * r + p1[i];
 			}
+			free(p1);
+			free(p2);
 		}
 		else                           //if the wavelength is equal to a wavelength in header file
 		{
@@ -327,7 +330,7 @@ public:
 		}
 		return true;
 	}
-	//untested
+
 	//providing the left and the right bound wavelength, return baseline-corrected band height
 	bool height(double lb, double rb, double bandwavelength, T* result){
 
@@ -343,6 +346,8 @@ public:
 
 		baseline_band(lb, rb, lp, rp, bandwavelength, result);
 
+		free(lp);
+		free(rp);
 		return true;
 	}
 
@@ -418,6 +423,10 @@ public:
 			std::swap(cur,cur2);		//swap the band pointers
 		}
 
+		free(lp);
+		free(rp);
+		free(cur);
+		free(cur2);
 		return true;
 	}
 
@@ -437,6 +446,9 @@ public:
 			else
 				result[i] = p1[i] / p2[i];
 		}
+
+		free(p1);
+		free(p2);
 		return true;	
 	}
 	
@@ -457,6 +469,9 @@ public:
 			else
 				result[i] = p1[i] / p2[i];
 		}
+
+		free(p1);
+		free(p2);
 		return true;	
 	}		
 	
@@ -477,8 +492,108 @@ public:
 			else
 				result[i] = p1[i] / p2[i];
 		}
+
+		free(p1);
+		free(p2);
 		return true;	
 	}		
+
+	//x * f(x)
+	bool x_area(double lb, double rb, double lab, double rab, T* result){
+		T* lp;	//left band pointer
+		T* rp;	//right band pointer
+		T* cur;		//current band 1
+		T* cur2;	//current band 2
+
+		unsigned XY = R[0] * R[1];
+		unsigned S = XY * sizeof(T);
+
+		lp = (T*) malloc(S);			//memory allocation
+		rp = (T*) malloc(S);
+		cur = (T*) malloc(S);
+		cur2 = (T*) malloc(S);
+
+		memset(result, (char)0, S);
+
+		//find the wavelenght position in the whole band
+		unsigned int n = w.size();
+		unsigned int ai = 0;		//left bound position
+		unsigned int bi = n - 1;		//right bound position
+
+		//to make sure the left and the right bound are in the bandwidth
+		if (lb < w[0] || rb < w[0] || lb > w[n-1] || rb >w[n-1]){
+			std::cout<<"ERROR: left bound or right bound out of bandwidth"<<std::endl;
+			exit(1);
+		}
+		//to make sure rigth bound is bigger than left bound
+		else if(lb > rb){
+			std::cout<<"ERROR: right bound should be bigger than left bound"<<std::endl;
+			exit(1);
+		}
+
+		//get the position of lb and rb
+		while (lab >= w[ai]){
+			ai++;
+		}
+		while (rab <= w[bi]){
+			bi--;
+		}
+
+		band(lp, lb);
+		band(rp, rb);
+
+		//calculate the beginning and the ending part
+		baseline_band(lb, rb, lp, rp, rab, cur2);		//ending part
+		baseline_band(lb, rb, lp, rp, w[bi], cur);
+		for(unsigned j = 0; j < XY; j++){
+			result[j] += (rab - w[bi]) * (rab + w[bi]) * (cur[j] + cur2[j]) / 4.0;
+		}
+		baseline_band(lb, rb, lp, rp, lab, cur2);		//beginnning part
+		baseline_band(lb, rb, lp, rp, w[ai], cur);
+		for(unsigned j = 0; j < XY; j++){	
+			result[j] += (w[ai] - lab) * (w[ai] + lab) * (cur[j] + cur2[j]) / 4.0;
+		}
+
+		//calculate f(x) times x
+		ai++;
+		for(unsigned i = ai; i <= bi ;i++)
+		{
+			baseline_band(lb, rb, lp, rp, w[ai], cur2);
+			for(unsigned j = 0; j < XY; j++)
+			{
+				result[j] += (w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * (cur[j] + cur2[j]) / 4.0; 
+			}
+			std::swap(cur,cur2);		//swap the band pointers
+		}
+
+		free(lp);
+		free(rp);
+		free(cur);
+		free(cur2);
+		return true;
+	}
+
+	//centroid point
+	bool cpoint(double lb, double rb, double lab, double rab, T* result){
+		T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T));	
+		T* p2 = (T*)malloc(R[0] * R[1] * sizeof(T));	
+		
+		//get the area and the peak band
+		x_area(lb, rb, lab, rab, p1);
+		area(lb, rb, lab, rab, p2);
+		//calculate the ratio in result
+		for(unsigned i = 0; i < R[0] * R[1]; i++){
+			if(p1[i] == 0 && p2[i] ==0)
+				result[i] = 1;
+			else
+				result[i] = p1[i] / p2[i];
+		}
+
+		free(p1);
+		free(p2);
+		return true;	
+	}
+
 	//create mask file
 	bool mask(unsigned char* p, double mask_band, double threshold){
 
@@ -492,6 +607,7 @@ public:
 					p[i] = 255;
 		}
 
+		free(temp);
 		return true;
 
 	}
diff --git a/envi/envi.h b/envi/envi.h
index 6914b72..de7603f 100644
--- a/envi/envi.h
+++ b/envi/envi.h
@@ -419,6 +419,42 @@ public:
 		return false;
 	}
 
+		//center of gravity
+	bool cpoint(double lb1, double rb1, double lab1, double rab1, void* result){
+		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
+			if(header.data_type ==envi_header::float32)
+				return ((bsq<float>*)file)->cpoint(lb1, rb1, lab1, rab1, (float*)result);
+			else if(header.data_type == envi_header::float64)
+				return ((bsq<double>*)file)->cpoint(lb1, rb1, lab1, rab1, (double*)result);
+			else
+				std::cout<<"ERROR: unidentified data type"<<std::endl;
+		}
+
+		else if(header.interleave == envi_header::BIL){		//if the infile is bil file
+			if(header.data_type ==envi_header::float32)
+				return ((bil<float>*)file)->cpoint(lb1, rb1, lab1, rab1, (float*)result);
+			else if(header.data_type == envi_header::float64)
+				return ((bil<double>*)file)->cpoint(lb1, rb1, lab1, rab1, (double*)result);
+			else
+				std::cout<<"ERROR: unidentified data type"<<std::endl;
+		}
+
+		else if(header.interleave == envi_header::BIP){		//if the infile is bip file
+			if(header.data_type ==envi_header::float32)
+				return ((bip<float>*)file)->cpoint(lb1, rb1, lab1, rab1,(float*)result);
+			else if(header.data_type == envi_header::float64)
+				return ((bip<double>*)file)->cpoint(lb1, rb1, lab1, rab1, (double*)result);
+			else
+				std::cout<<"ERROR: unidentified data type"<<std::endl;
+		}
+
+		else{
+			std::cout<<"ERROR: unidentified file type"<<std::endl;
+			exit(1);
+		}
+		return false;
+	}
+
 	bool close(){
 		if(header.interleave == envi_header::BSQ){
 			if(header.data_type ==envi_header::float32)
--
libgit2 0.21.4