stim/image/image.h

#ifndef STIM_IMAGE_H
#define STIM_IMAGE_H
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <vector>
#include <iostream>
#include <limits>
namespace stim{
/// This static class provides the STIM interface for loading, saving, and storing 2D images.
/// Data is stored in an interleaved (BIP) format (default for saving and loading is RGB).
//currently this interface uses CImg
//	T = data type (usually unsigned char)
template <class T>
class image{
	//cimg_library::CImg<T> img;
	T* img;										//pointer to the image data (assumes RGB for loading/saving)
	size_t R[3];
	size_t X() const { return R[1]; }
	size_t Y() const { return R[2]; }
	size_t C() const { return R[0]; }
	size_t bytes(){ return size() * sizeof(T); }
	void init(){								//initializes all variables, assumes no memory is allocated
		memset(R, 0, sizeof(size_t) * 3);		//set the resolution and number of channels to zero
		img = NULL;
	}
	void unalloc(){								//frees any resources associated with the image
		if(img)	free(img);						//if memory has been allocated, free it
		img=NULL;		
	}
	void clear(){								//clears all image data
		unalloc();								//unallocate previous memory
		init();									//re-initialize the variables
	}
	void allocate(){
		unalloc();
		img = (T*) malloc( bytes() );	//allocate memory
		memset(img, 0, bytes());
	}
	void allocate(size_t x, size_t y, size_t c){	//allocate memory based on the resolution
		unalloc();
		R[0] = c; R[1] = x; R[2] = y;				//set the resolution
		allocate();									//allocate memory
	}
	size_t idx(size_t x, size_t y, size_t c = 0){
		return y * C() * X() + x * C() + c;
	}
	int cv_type(){
		// The following is C++ 11 code, but causes problems on some compilers (ex. nvcc). Below is my best approximation to a solution
		//if(std::is_same<T, unsigned char>::value)	return CV_MAKETYPE(CV_8U, (int)C());
		//if(std::is_same<T, char>::value)			return CV_MAKETYPE(CV_8S, (int)C());
		//if(std::is_same<T, unsigned short>::value)	return CV_MAKETYPE(CV_16U, (int)C());
		//if(std::is_same<T, short>::value)			return CV_MAKETYPE(CV_16S, (int)C());
		//if(std::is_same<T, int>::value)				return CV_MAKETYPE(CV_32S, (int)C());
		//if(std::is_same<T, float>::value)			return CV_MAKETYPE(CV_32F, (int)C());
		//if(std::is_same<T, double>::value)			return CV_MAKETYPE(CV_64F, (int)C());
		if(typeid(T) == typeid(unsigned char))		return CV_MAKETYPE(CV_8U, (int)C());
		if(typeid(T) == typeid(char))				return CV_MAKETYPE(CV_8S, (int)C());
		if(typeid(T) == typeid(unsigned short))		return CV_MAKETYPE(CV_16U, (int)C());
		if(typeid(T) == typeid(short))				return CV_MAKETYPE(CV_16S, (int)C());
		if(typeid(T) == typeid(int))				return CV_MAKETYPE(CV_32S, (int)C());
		if(typeid(T) == typeid(float))				return CV_MAKETYPE(CV_32F, (int)C());
		if(typeid(T) == typeid(double))				return CV_MAKETYPE(CV_64F, (int)C());
		std::cout<<"ERROR in stim::image::cv_type - no valid data type found"<<std::endl;
		exit(1);
	}
	/// Returns the value for "white" based on the dynamic range (assumes white is 1.0 for floating point images)
	T white(){
		// The following is C++ 11 code, but causes problems on some compilers (ex. nvcc). Below is my best approximation to a solution
		//if(std::is_same<T, unsigned char>::value)		return UCHAR_MAX;
		//if(std::is_same<T, unsigned short>::value)		return SHRT_MAX;
		//if(std::is_same<T, unsigned>::value)			return UINT_MAX;
		//if(std::is_same<T, unsigned long>::value)		return ULONG_MAX;
		//if(std::is_same<T, unsigned long long>::value)	return ULLONG_MAX;
		//if(std::is_same<T, float>::value)				return 1.0f;
		//if(std::is_same<T, double>::value)				return 1.0;
		if(typeid(T) == typeid(unsigned char))		return UCHAR_MAX;
		if(typeid(T) == typeid(unsigned short))		return SHRT_MAX;
		if(typeid(T) == typeid(unsigned))			return UINT_MAX;
		if(typeid(T) == typeid(unsigned long))		return ULONG_MAX;
		if(typeid(T) == typeid(unsigned long long))	return ULLONG_MAX;
		if(typeid(T) == typeid(float))				return 1.0f;
		if(typeid(T) == typeid(double))				return 1.0;
		std::cout<<"ERROR in stim::image::white - no white value known for this data type"<<std::endl;
		exit(1);
	}
public:
	/// Default constructor - creates an empty image object
	image(){
		init();						//initialize all variables to zero, don't allocate any memory
	}							
	/// Constructor with a filename - loads the specified file
	image(std::string filename){				//constructor initialize the image with an image file
		init();
		load(filename);
	}
	/// Create a new image from scratch given a number of samples and channels
	image(size_t x, size_t y = 1, size_t c = 1){
		init();
		allocate(x, y, c);
	}
	/// Create a new image with the data given in 'data'
	image(T* data, size_t x, size_t y, size_t c = 1){
		init();
		allocate(x, y, c);
		memcpy(img, data, bytes());
	}
	/// Copy constructor - duplicates an image object
	image(const stim::image<T> &I){
		init();
		allocate(I.X(), I.Y(), I.C());
		memcpy(img, I.img, bytes());
	}
	/// Destructor - clear memory
	~image(){
		free(img);
	}
	stim::image<T>& operator=(const stim::image<T>& I){
		if(&I == this)									//handle self-assignment
			return *this;
		allocate(I.X(), I.Y(), I.C());
		memcpy(img, I.img, bytes());
		return *this;
	}
	/// Load an image from a file
	void load(std::string filename){
		cv::Mat cvImage = cv::imread(filename, CV_LOAD_IMAGE_UNCHANGED);	//use OpenCV to open the image file
		if(!cvImage.data){
			std::cout<<"ERROR stim::image::load() - unable to find image "<<filename<<" ["<<__FILE__<<" (line "<<__LINE__<<")]"<<std::endl;
			exit(1);
		}
		allocate(cvImage.cols, cvImage.rows, cvImage.channels());			//allocate space for the image
		T* cv_ptr = (T*) cvImage.data;
		if(C() == 1)
		{
	//if this is a single-color image, just copy the data
			memcpy(img, cv_ptr, bytes());     
		}
		if(C() == 3)
		{														//if this is a 3-color image, OpenCV uses BGR interleaving
			set_interleaved_bgr(cv_ptr, X(), Y());
		}
		cvImage.release();	
	}
	//save a file
	void save(std::string filename){
		//OpenCV uses an interleaved format, so convert first and then output
		T* buffer = (T*) malloc(bytes());
		if(C() == 1)
			memcpy(buffer, img, bytes());
		else if(C() == 3)
			get_interleaved_bgr(buffer);
		cv::Mat cvImage((int)Y(), (int)X(), cv_type(), buffer);
		cv::imwrite(filename, cvImage);
		cvImage.release();
		free(buffer);
	}
	//create an image from an interleaved buffer
	void set_interleaved_rgb(T* buffer, size_t width, size_t height, size_t channels = 3){
		allocate(width, height, channels);
		memcpy(img, buffer, bytes());
	}
	void set_interleaved_bgr(T* buffer, size_t width, size_t height, size_t channels = 3){
		allocate(width, height, channels);
		for(size_t c = 0; c < C(); c++){								//copy directly
			for(size_t y = 0; y < Y(); y++){
				for(size_t x = 0; x < X(); x++){
					img[idx(x, y, c)] = buffer[y * X() * C() + x * C() + (2-c)];
				}
			}
		}
	}
	void get_interleaved_bgr(T* data){
		//for each channel
		for(size_t y = 0; y < Y(); y++){
			for(size_t x = 0; x < X(); x++){
				for(size_t c = 0; c < C(); c++){
					data[y * X() * C() + x * C() + (2-c)] = img[idx(x, y, c)];
				}
			}
		}
	}
	void get_interleaved_rgb(T* data){
		memcpy(data, img, bytes());
	}
	//copies data in the given channel order as a non-interleaved image
	void get_noninterleaved(T* data){
		//for each channel
		for(size_t y = 0; y < Y(); y++){
			for(size_t x = 0; x < X(); x++){
				for(size_t c = 0; c < C(); c++){
					data[c * Y() * X() + y * X() + x] = img[idx(x, y, c)];
				}
			}
		}
	}
	image<T> channel(size_t c){
		//create a new image
		image<T> r(X(), Y(), 1);
		for(size_t x = 0; x < X(); x++){
			for(size_t y = 0; y < Y(); y++){
				r.img[r.idx(x, y, 0)] = img[idx(x, y, c)];
			}
		}
		return r;
	}
	T& operator()(size_t x, size_t y, size_t c = 0){
		return img[idx(x, y, c)];
	}
	/// Set all elements in the image to a given scalar value
	/// @param v is the value used to set all values in the image
	image<T> operator=(T v){
		size_t N = size();
		for(size_t n = 0; n < N; n++)
			img[n] = v;
		return *this;
	}
	/// Copy the given data to the specified channel
	/// @param c is the channel number that the data will be copied to
	/// @param buffer is a pointer to the image to be copied to channel c
	void set_channel(T* buffer, size_t c){
		size_t x, y;
		for(y = 0; y < Y(); y++){
			for(x = 0; x < X(); x++){
				img[idx(x, y, c)] = buffer[c];
			}
		}
	}
	size_t channels(){
		return C();
	}
	size_t width(){
		return X();
	}
	size_t height(){
		return Y();
	}
	T* data(){
		return img;
	}
	//returns the size (number of values) of the image
	size_t size(){ return C() * X() * Y(); }
	/// Returns the number of nonzero values
	size_t nnz(){
		size_t N = X() * Y() * C();
		size_t nz = 0;
		for(size_t n = 0; n < N; n++)
			if(img[n] != 0) nz++;
		return nz;	//return the number of nonzero pixels
	}
	//this function returns indices of pixels that have nonzero values
	std::vector<size_t> sparse_idx(){
		std::vector<size_t> s;				//allocate an array
		s.resize(nnz());					//allocate space in the array
		size_t N = size();
		//size_t C = channels();
		//T* ptr = img.data();				//get a pointer to the image data
		size_t i = 0;
		for(size_t n = 0; n < N; n++){
			if(img[n] != 0){
				s[i] = n;
				i++;
			}
		}
		return s;			//return the index list
	}
	/// Returns the maximum pixel value in the image
	T maxv(){
		T max_val = img[0];				//initialize the maximum value to the first one
		size_t N = size();	//get the number of pixels
		for (size_t n=0; n<N; n++){		//for every value
			if (img[n] > max_val){			//if the value is higher than the current max
				max_val = img[n];
			}
		}
		return max_val;
	}
	/// Returns the maximum pixel value in the image
	T minv(){
		T min_val = img[0];				//initialize the maximum value to the first one
		size_t N = size();	//get the number of pixels
		for (size_t n=0; n<N; n++){		//for every value
			if (img[n] < min_val){			//if the value is higher than the current max
				min_val = img[n];
			}
		}
		return min_val;
	}
	/// Invert an image by calculating I1 = alpha - I0, where alpha is the maximum image value
	image<T> invert(T white_val){
		size_t N = size();						//calculate the total number of values in the image
		image<T> r(X(), Y(), C());				//allocate space for the resulting image
		for(size_t n = 0; n < N; n++)
			r.img[n] = white_val - img[n];		//perform the inversion
		return r;								//return the inverted image
	}
	
	/// Invert an image by calculating I1 = alpha - I0, where alpha is the maximum image value
	image<T> invert(){
		size_t N = size();						//calculate the total number of values in the image
		image<T> r(X(), Y(), C());				//allocate space for the resulting image
		T white_val = maxv();
		for(size_t n = 0; n < N; n++)
			r.img[n] = white_val - img[n];		//perform the inversion
		return r;								//return the inverted image
	}
	///crops the image from x1 to x0 and y1 to y0 and returns a new (smaller) image.
	image<T> crop(int x0, int x1, int y0, int y1)
	{
		
		image<T> ret(x1-x0, y1-y0, C());
		int newWidth = x1-x0;
		int destidx, srcidx;
		///for each row, cut what amount of data from the original and put it into the new copy.
		for(int i = 0; i < (y1-y0); i++)
		{
			destidx = i*newWidth*C(); ///destination index one per each row
			srcidx = ((i+(y0))*X()+x0)*C(); ///source index, one per each row.
			memcpy(&ret.img[destidx], &img[srcidx], sizeof(T)*newWidth*C());
		}		
		return ret;
	}
	image<T> srgb2lab(){
		std::cout<<"ERROR stim::image::srgb2lab - function has been broken, re-implement."<<std::endl;
		exit(1);
	}
	image<T> convolve2(image<T> mask){
		std::cout<<"ERROR stim::image::convolve2 - function has been broken, and shouldn't really be in here."<<std::endl;
		exit(1);
	}
	image<T> rotate(float angle, float cx, float cy){
		std::cout<<"ERROR stim::image::rotate - function has been broken, and shouldn't really be in here."<<std::endl;
		exit(1);
	}
	// leila's code for non_interleaving data in 3D
	//create an data set from an interleaved buffer
	void set_interleaved3(T* buffer, size_t width, size_t height, size_t depth, size_t channels = 3){
		std::cout<<"ERROR stim::image::set_interleaved3 - stim::image no longer supports 3D images."<<std::endl;
		exit(1);
	}
};
};		//end namespace stim
#endif