tira/envi/agilent_binary.h

//make sure that this header file is only loaded once
#ifndef STIM_AGILENT_BINARY_H
#define STIM_AGILENT_BINARY_H
#include <string>
#include <fstream>
#include <complex>
#include <cstring>
#include <chrono>
//CUDA
//#ifdef CUDA_FOUND
#include <cuda_runtime.h>
#include "cufft.h"
#include <stim/cuda/cudatools/error.h>
#include <stim/envi/envi_header.h>
//#endif
namespace stim{
template<typename T>
class agilent_binary{
protected:
	std::string fname;
	T* ptr;														//pointer to the image data
	size_t R[3];												//size of the binary image in X, Y, and Z
	static const size_t header = 1020;							//header size
	double Z[2];												//range of z values (position or wavelength)
public:
	size_t size(){
		return (size_t)R[0] * (size_t)R[1] * (size_t)R[2];
	}
	size_t bytes(){
		return size() * sizeof(T);
	}
	void alloc(){
		if (ptr != NULL) free(ptr);
		ptr = NULL;
		ptr = (T*) malloc(bytes());
	}
	void alloc(size_t x, size_t y, size_t z){		
		R[0] = x;
		R[1] = y;
		R[2] = z;
		alloc();
	}
	char* data() {
		return (char*)ptr;
	}
	size_t dim(size_t i){
		return R[i];
	}
	/// Create a deep copy of an agileng_binary object
	void deep_copy(agilent_binary<T>* dst, const agilent_binary<T>* src){
		dst->alloc(src->R[0], src->R[1], src->R[2]);			//allocate memory
		memcpy(dst->ptr, src->ptr, bytes());					//copy the data
		memcpy(dst->Z, src->Z, sizeof(double) * 2);				//copy the data z range
	}
	/// Default constructor, sets the resolution to zero and the data pointer to NULL
	agilent_binary(){
		ptr = NULL;
		memset(R, 0, sizeof(size_t) * 3);				//set the resolution to zero
		memset(Z, 0, sizeof(double) * 2);		
	}
	/// Constructor with resolution
	agilent_binary(size_t x, size_t y, size_t z){
		ptr = NULL;
		alloc(x, y, z);
		memset(Z, 0, sizeof(double) * 2);
	}
	/// Constructor with filename
	agilent_binary(std::string filename){
		ptr = NULL;
		memset(Z, 0, sizeof(double) * 2);
		load(filename);
	}
	
	/// Copy constructor
	agilent_binary(const agilent_binary<T> &obj){
		ptr = NULL;
		deep_copy(this, &obj);
	}
	agilent_binary<T>& operator=(const agilent_binary<T> rhs){
		if(this != &rhs){								//check for self-assignment
			deep_copy(this, &rhs);						//make a deep copy
		}
		return *this;									//return the result
	}
	operator bool() {
		if (R[0] == 0 || R[1] == 0 || R[2] == 0)	return false;
		else return true;
	}
	~agilent_binary(){
		if(ptr != NULL)
			free(ptr);
	}
	void load(std::string filename){
		if(ptr != NULL) free(ptr);									//if memory has been allocated, free it
		ptr = NULL;
		fname = filename;											//save the filename
		short x, y, z;
		std::ifstream infile(fname, std::ios::binary);				//open the input file
		if (infile) {
			infile.seekg(9, std::ios::beg);							//seek past 9 bytes from the beginning of the file
			infile.read((char*)(&z), 2);							//read two bytes of data (the number of samples is stored as a 16-bit integer)
			infile.seekg(13, std::ios::cur);						//skip another 13 bytes
			infile.read((char*)(&x), 2);							//read the X and Y dimensions
			infile.read((char*)(&y), 2);
			infile.seekg(header, std::ios::beg);					//seek to the start of the data
			alloc(x, y, z);											//allocate the data
			infile.read((char*)ptr, bytes());						//read the data		
			infile.close();											//close the file
			Z[0] = 1;
			Z[1] = (double)R[2];
		}
	}
	void save(std::string filename){
		std::ofstream outfile(filename, std::ios::binary);			//create an output file
		char zero = 0;
		for(size_t i = 0; i < 9; i++) outfile.write(&zero, 1);		//write 9 zeros
		outfile.write((char*)&R[2], 2);
		for(size_t i = 0; i < 13; i++) outfile.write(&zero, 1);		//write 13 zeros
		outfile.write((char*)&R[0], 2);
		outfile.write((char*)&R[1], 2);
		for(size_t i = 0; i < 992; i++) outfile.write(&zero, 1);		//write 992 zeros
		size_t b = bytes();
		outfile.write((char*)ptr, b);							//write the data to the output file
		outfile.close();
	}
	stim::envi_header create_header(){
		stim::envi_header header;
		header.samples = R[0];
		header.lines = R[1];
		header.bands = R[2];
		double z_delta = (double)(Z[1] - Z[0]) / (double)(R[2] - 1);
		header.wavelength.resize(R[2]);
		for(size_t i = 0; i < R[2]; i++)
			header.wavelength[i] = i * z_delta + Z[0];
		return header;
	}
	/// Subtract the mean from each pixel. Generally used for centering an interferogram.
	void meancenter(){
		size_t Z = R[2];											//store the number of bands
		size_t XY = R[0] * R[1];									//store the number of pixels in the image
		T sum = (T)0;
		T mean;
		for(size_t xy = 0; xy < XY; xy++){							//for each pixel
			sum = 0;
			for(size_t z = 0; z < Z; z++){							//for each band
				sum += ptr[ z * XY + xy ];							//add the band value to a running sum
			}
			mean = sum / (T)Z;										//calculate the pixel mean
			for(size_t z = 0; z < Z; z++){
				ptr[ z * XY + xy ] -= mean;							//subtract the mean from each band
			}
		}
	}
	/// adds n bands of zero padding to the end of the file
	void zeropad(size_t n){
		size_t newZ = R[2] + n;
		T* temp = (T*) calloc(R[0] * R[1] * newZ, sizeof(T));	//allocate space for the new image
		memcpy(temp, ptr, size() * sizeof(T));					//copy the old data to the new image
		
		free(ptr);												//free the old data
		ptr = temp;												//swap in the new data
		R[2] = newZ;											//set the z-dimension to the new zero value
	}
	//pads to the nearest power-of-two
	void zeropad(){
		size_t newZ = (size_t)pow(2, ceil(log(R[2])/log(2)));			//find the nearest power-of-two
		size_t n = newZ - R[2];									//calculate the number of bands to add
		zeropad(n);												//add the padding
	}
	/// Calculate the absorbance spectrum from the transmission spectrum given a background
	void absorbance(stim::agilent_binary<T>* background){
		size_t N = size();											//calculate the number of values to be ratioed
		if(N != background->size()){
			std::cerr<<"ERROR in stim::agilent_binary::absorbance() - transmission image size doesn't match background"<<std::endl;
			exit(1);
		}
		for(size_t i = 0; i < N; i++)
			ptr[i] = -log10(ptr[i] / background->ptr[i]);
	}
	//crops the image down to a set number of samples
	void crop(size_t n) {
		if (n < R[2]) {											//if the requested size is smaller than the image
			R[2] = n;											//update the number of bands
			T* old_ptr = ptr;									//store the old pointer
			alloc();											//allocate space for the new image
			memcpy(ptr, old_ptr, bytes());						//copy the old data to the new image
			free(old_ptr);										//free the old data
		}
	}
//#ifdef CUDA_FOUND
	/// Perform an FFT and return a binary file with bands in the specified range
	agilent_binary<T> fft(double band_min, double band_max, double ELWN = 15798, int UDR = 2, int device = 0){
		auto total_start = std::chrono::high_resolution_clock::now();
		auto start = std::chrono::high_resolution_clock::now();
		T* cpu_data = (T*) malloc( bytes() );										//allocate space for the transposed data
		for(size_t b = 0; b < R[2]; b++){
			for(size_t x = 0; x < R[0] * R[1]; x++){
				cpu_data[x * R[2] + b] = ptr[b * R[0] * R[1] + x];
			}
		}
		auto end = std::chrono::high_resolution_clock::now();
		std::chrono::duration<double> diff = end-start;
       // std::cout << "Transpose data: " << diff.count() << " s\n";
		start = std::chrono::high_resolution_clock::now();
		if (device >= 0) {														//if a CUDA device is specified
			int dev_count;
			HANDLE_ERROR(cudaGetDeviceCount(&dev_count));						//get the number of CUDA devices
			//std::cout << "Number of CUDA devices: " << dev_count << std::endl;		//output the number of CUDA devices
			cudaDeviceProp prop;
			//std::cout << "CUDA devices----" << std::endl;
			for (int d = 0; d < dev_count; d++) {									//for each CUDA device
				cudaGetDeviceProperties(&prop, d);									//get the property of the first device
																					//float cc = prop.major + prop.minor / 10.0f;						//calculate the compute capability
				//std::cout << d << ":  [" << prop.major << "." << prop.minor << "]      " << prop.name << std::endl;	//display the device information
																													//if(cc > best_device_cc){
																													//	best_device_cc = cc;										//if this is better than the previous device, use it
																													//	best_device_id = d;
																													//}
			}
			if (dev_count > 0 && dev_count > device) {							//if the first device is not an emulator
				cudaGetDeviceProperties(&prop, device);						//get the property of the requested CUDA device
				if (prop.major != 9999) {
					//std::cout << "Using device " << device << std::endl;
					HANDLE_ERROR(cudaSetDevice(device));
				}
			}
		}
		cufftHandle plan;															//allocate space for a cufft plan
		cufftReal* gpu_data;														//create a pointer to the data
		size_t batch = R[0] * R[1];													//calculate the batch size (X * Y)
		HANDLE_ERROR(cudaMalloc((void**)&gpu_data, bytes()));						//allocate space on the GPU
		HANDLE_ERROR(cudaMemcpy(gpu_data, cpu_data, bytes(), cudaMemcpyHostToDevice));	//copy the data to the GPU
		cufftComplex* gpu_fft;
		HANDLE_ERROR(cudaMalloc((void**)&gpu_fft, R[0] * R[1] * (R[2]/2 + 1) * sizeof(cufftComplex)));
		end = std::chrono::high_resolution_clock::now();
		diff = end-start;
		//std::cout << "Allocate/copy: " << diff.count() << " s\n";
		start = std::chrono::high_resolution_clock::now();
		int N[1];					//create an array with the interferogram size (required for cuFFT input)
		N[0] = (int)R[2];				//set the only array value to the interferogram size
		if(cufftPlanMany(&plan, 1, N, NULL, 1, (int)R[2], NULL, 1, (int)R[2], CUFFT_R2C, (int)batch) != CUFFT_SUCCESS){
			std::cout<<"cuFFT Error: unable to create 1D plan."<<std::endl;
			exit(1);
		}
		end = std::chrono::high_resolution_clock::now();
		diff = end-start;
		//std::cout << "Create a plan: " << diff.count() << " s\n";
		start = std::chrono::high_resolution_clock::now();
		if (cufftExecR2C(plan, gpu_data, gpu_fft) != CUFFT_SUCCESS){		//execute the (implicitly forward) transform
			std::cout<<"CUFFT error: ExecR2C Forward failed";
			exit(1);
		}
		end = std::chrono::high_resolution_clock::now();
		diff = end-start;
		//std::cout << "Perform FFT: " << diff.count() << " s\n";
		start = std::chrono::high_resolution_clock::now();
		std::complex<T>* cpu_fft = (std::complex<T>*) malloc( R[0] * R[1] * (R[2]/2+1) * sizeof(std::complex<T>) );
		HANDLE_ERROR(cudaMemcpy(cpu_fft, gpu_fft, R[0] * R[1] * (R[2]/2+1) * sizeof(cufftComplex), cudaMemcpyDeviceToHost));	//copy data from the host to the device
		//double int_delta = 0.00012656;									//interferogram sample spacing in centimeters
		double int_delta = (1.0 / ELWN) * ((double)UDR / 2.0);				//calculate the interferogram spacing
		double int_length = int_delta * R[2];								//interferogram length in centimeters
		double fft_delta = 1/int_length;									//spectrum spacing (in inverse centimeters, wavenumber)
		double fft_max = fft_delta * R[2]/2;								//get the maximum wavenumber value supported by the specified number of interferogram samples
		if(band_max > fft_max) band_max = fft_max;							//the user gave a band outside of the FFT range, reset the band to the maximum available
		if (band_min < 0) band_min = 0;
		size_t start_i = (size_t)std::ceil(band_min / fft_delta);				//calculate the first band to store
		size_t size_i = (size_t)std::floor(band_max / fft_delta) - start_i;		//calculate the number of bands to store
		size_t end_i = start_i + size_i;										//last band number
		agilent_binary<T> result(R[0], R[1], size_i);
		result.Z[0] = start_i * fft_delta;										//set the range for the FFT result
		result.Z[1] = end_i * fft_delta;
		for(size_t b = start_i; b < end_i; b++){
			for(size_t x = 0; x < R[0] * R[1]; x++){
				result.ptr[(b - start_i) * R[0] * R[1] + x] = abs(cpu_fft[x * (R[2]/2+1) + b]);
			}
		}
		end = std::chrono::high_resolution_clock::now();
		diff = end-start;
		//std::cout << "Transpose/Crop: " << diff.count() << " s\n";
		auto total_end = std::chrono::high_resolution_clock::now();
		diff = total_end-total_start;
		cufftDestroy(plan);
		HANDLE_ERROR(cudaFree(gpu_data));
		HANDLE_ERROR(cudaFree(gpu_fft));
		free(cpu_data);
		free(cpu_fft);
		return result;
	}
	//saves the binary as an ENVI file with a BIP interleave format
	int bip(T* bip_ptr){
		//std::ofstream out(outfile.c_str(), std::ios::binary);			//create a binary file stream for output
		size_t XY = R[0] * R[1];
		size_t B = R[2];
		size_t b;
		for(size_t xy = 0; xy < XY; xy++){
			for(b = 0; b < B; b++){
				bip_ptr[xy * B + b] = ptr[b * XY + xy];
			}
		}
		return 0;
	}
//#endif
};
}
#endif