Ziqi added functions to attain area

heziqi
1 parent 4a6f666c
Showing 4 changed files with 390 additions and 1 deletions Show diff stats
envi/bil.h
envi/bip.h
envi/bsq.h
envi/envi.h
@@ -387,6 +387,124 @@ public:
 		return true;
 	}
  
+
+	//providing the left and the right bound data, return baseline-corrected band height
+	bool baseline_band(double lb, double rb, T* lp, T* rp, double wavelength, T* result){
+
+		unsigned XY = R[0] * R[1];
+		band(result, wavelength);		//get band
+
+		//perform the baseline correction
+		double r = (double) (wavelength - lb) / (double) (rb - lb);
+		for(unsigned i=0; i < XY; i++){
+			result[i] =(T) (result[i] - (rp[i] - lp[i]) * r - lp[i] );
+		}
+		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* lp;
+		T* rp;
+		T* result;
+
+		lp = (T*) malloc(sizeof(R[0] * R[1] * sizeof(T)));			//memory allocation
+		rp = (T*) malloc(sizeof(R[0] * R[1] * sizeof(T)));
+		heigth = (T*) malloc(sizeof(R[0] * R[1] * sizeof(T)));
+
+		band(lp, lb);			//get the data of the left bounbd and the right bound
+		band(rp, rb);		
+
+		baseline_band(lb, rb, lp, rp, bandwavelength, result);
+
+		return true;
+	}
+
+	//untested
+	//calculate the area between two bound point(including baseline correction)
+	bool area(double lb, double rb){
+
+		T* lp;	//left band pointer
+		T* rp;	//right band pointer
+		T* cur;		//current band 1
+		T* cur2;	//current band 2
+		T* result;		//area result
+
+		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);
+		result = (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 (lb >= w[ai]){
+			ai++;
+		}
+		while (rb <= w[bi]){
+			bi--;
+		}
+
+		band(lp, lb);
+		band(rp, rb);
+
+		//calculate the beginning and the ending part
+		baseline_band(lb, rb, lp, rp, w[bi], cur);		//ending part
+		for(unsigned j = 0; j < XY; j++){
+			result[j] += (rb - w[bi]) * cur[j] / 2.0;
+		}
+		baseline_band(lb, rb, lp, rp, w[ai], cur);		//beginnning part
+		for(unsigned j = 0; j < XY; j++){	
+			result[j] += (w[ai] - lb) * cur[j] / 2.0;
+		}
+		//calculate the area
+		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]) * (cur[j] + cur2[j]) / 2.0; 
+			}
+			std::swap(cur,cur2);		//swap the band pointers
+		}
+		///////////////for testing use
+		std::ofstream text("area_result.txt",std::ios::out);
+		for (unsigned i = 0; i < 640; i++)	
+		{
+			for (unsigned j = 0; j < 384; j++)
+			{
+				text<<result[i*384+j]<<", ";
+			}
+			text<<std::endl;
+		}
+		///////////////
+		return true;
+	}
+
 	//create mask file
 	bool mask(unsigned char* p, double mask_band, double threshold){
  
@@ -477,6 +477,123 @@ public:
 		return true;
 	}
  
+	//providing the left and the right bound data, return baseline-corrected band height
+	bool baseline_band(double lb, double rb, T* lp, T* rp, double wavelength, T* result){
+
+		unsigned XY = R[0] * R[1];
+		band(result, wavelength);		//get band
+
+		//perform the baseline correction
+		double r = (double) (wavelength - lb) / (double) (rb - lb);
+		for(unsigned i=0; i < XY; i++){
+			result[i] =(T) (result[i] - (rp[i] - lp[i]) * r - lp[i] );
+		}
+		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* lp;
+		T* rp;
+		T* result;
+
+		lp = (T*) malloc(sizeof(R[0] * R[1] * sizeof(T)));			//memory allocation
+		rp = (T*) malloc(sizeof(R[0] * R[1] * sizeof(T)));
+		heigth = (T*) malloc(sizeof(R[0] * R[1] * sizeof(T)));
+
+		band(lp, lb);			//get the data of the left bounbd and the right bound
+		band(rp, rb);		
+
+		baseline_band(lb, rb, lp, rp, bandwavelength, result);
+
+		return true;
+	}
+
+	//untested
+	//calculate the area between two bound point(including baseline correction)
+	bool area(double lb, double rb){
+
+		T* lp;	//left band pointer
+		T* rp;	//right band pointer
+		T* cur;		//current band 1
+		T* cur2;	//current band 2
+		T* result;		//area result
+
+		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);
+		result = (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 (lb >= w[ai]){
+			ai++;
+		}
+		while (rb <= w[bi]){
+			bi--;
+		}
+
+		band(lp, lb);
+		band(rp, rb);
+
+		//calculate the beginning and the ending part
+		baseline_band(lb, rb, lp, rp, w[bi], cur);		//ending part
+		for(unsigned j = 0; j < XY; j++){
+			result[j] += (rb - w[bi]) * cur[j] / 2.0;
+		}
+		baseline_band(lb, rb, lp, rp, w[ai], cur);		//beginnning part
+		for(unsigned j = 0; j < XY; j++){	
+			result[j] += (w[ai] - lb) * cur[j] / 2.0;
+		}
+		//calculate the area
+		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]) * (cur[j] + cur2[j]) / 2.0; 
+			}
+			std::swap(cur,cur2);		//swap the band pointers
+		}
+		///////////////for testing use
+		std::ofstream text("area_result.txt",std::ios::out);
+		for (unsigned i = 0; i < 640; i++)	
+		{
+			for (unsigned j = 0; j < 384; j++)
+			{
+				text<<result[i*384+j]<<", ";
+			}
+			text<<std::endl;
+		}
+		///////////////
+		return true;
+	}
+
 	//create mask file
 	bool mask(unsigned char* p, double mask_band, double threshold){
  
@@ -309,7 +309,125 @@ public:
 		target.close();
 		return true;
 	}
-	
+
+	//providing the left and the right bound data, return baseline-corrected band height
+	bool baseline_band(double lb, double rb, T* lp, T* rp, double wavelength, T* result){
+
+		unsigned XY = R[0] * R[1];
+		band(result, wavelength);		//get band
+
+		//perform the baseline correction
+		double r = (double) (wavelength - lb) / (double) (rb - lb);
+		for(unsigned i=0; i < XY; i++){
+			result[i] =(T) (result[i] - (rp[i] - lp[i]) * r - lp[i] );
+		}
+		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* lp;
+		T* rp;
+		T* result;
+
+		lp = (T*) malloc(sizeof(R[0] * R[1] * sizeof(T)));			//memory allocation
+		rp = (T*) malloc(sizeof(R[0] * R[1] * sizeof(T)));
+		heigth = (T*) malloc(sizeof(R[0] * R[1] * sizeof(T)));
+
+		band(lp, lb);			//get the data of the left bounbd and the right bound
+		band(rp, rb);		
+
+		baseline_band(lb, rb, lp, rp, bandwavelength, result);
+
+		return true;
+	}
+
+	//untested
+	//calculate the area between two bound point(including baseline correction)
+	bool area(double lb, double rb){
+
+		T* lp;	//left band pointer
+		T* rp;	//right band pointer
+		T* cur;		//current band 1
+		T* cur2;	//current band 2
+		T* result;		//area result
+
+		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);
+		result = (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 (lb >= w[ai]){
+			ai++;
+		}
+		while (rb <= w[bi]){
+			bi--;
+		}
+
+		band(lp, lb);
+		band(rp, rb);
+
+		//calculate the beginning and the ending part
+		baseline_band(lb, rb, lp, rp, w[bi], cur);		//ending part
+		for(unsigned j = 0; j < XY; j++){
+			result[j] += (rb - w[bi]) * cur[j] / 2.0;
+		}
+		baseline_band(lb, rb, lp, rp, w[ai], cur);		//beginnning part
+		for(unsigned j = 0; j < XY; j++){	
+			result[j] += (w[ai] - lb) * cur[j] / 2.0;
+		}
+		//calculate the area
+		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]) * (cur[j] + cur2[j]) / 2.0; 
+			}
+			std::swap(cur,cur2);		//swap the band pointers
+		}
+		///////////////for testing use
+		std::ofstream text("area_result.txt",std::ios::out);
+		for (unsigned i = 0; i < 640; i++)	
+		{
+			for (unsigned j = 0; j < 384; j++)
+			{
+				text<<result[i*384+j]<<", ";
+			}
+			text<<std::endl;
+		}
+		///////////////
+		return true;
+	}
+
+
 	//create mask file
 	bool mask(unsigned char* p, double mask_band, double threshold){
  
@@ -310,6 +310,42 @@ public:
 		return false;
 	}
  
+	//do baseline correction and caculate the area between two bounds
+	bool area(double lb, double rb){
+		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
+			if(header.data_type ==envi_header::float32)
+				return ((bsq<float>*)file)->area(lb, rb);
+			else if(header.data_type == envi_header::float64)
+				return ((bsq<double>*)file)->area(lb, rb);
+			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)->area(lb, rb);
+			else if(header.data_type == envi_header::float64)
+				return ((bil<double>*)file)->area(lb, rb);
+			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)->area(lb, rb);
+			else if(header.data_type == envi_header::float64)
+				return ((bip<double>*)file)->area(lb, rb);
+			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)