codebase / stimlib

  #ifndef STIM_IMAGE_H
  #define STIM_IMAGE_H
  
  #ifdef _WIN32
  #undef max
  #endif
  
  #include <vector>
  #include <iostream>
  #include <limits>							//use limits and remove the MIN and MAX macros
  #include <typeinfo>
  #include <fstream>
  #include <cstring>
  #include <stim/image/CImg.h>
  
  
  #include <stim/parser/filename.h>
  
  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).
  
  template <class T>
  class image{
  
  	T* img;										//pointer to the image data (interleaved RGB for color)
  	size_t R[3];
  	bool interleaved = true;					//by default the data is interleaved
  
  	inline size_t X() const { return R[1]; }
  	inline size_t Y() const { return R[2]; }
  	inline size_t C() const { return R[0]; }
  
  	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
  	}
  
  
  	void clear(){								//clears all image data
  		unalloc();								//unallocate previous memory
  		init();									//re-initialize the variables
  	}
  
  	void allocate(){
  		unalloc();
  		img = (T*) malloc( sizeof(T) * R[0] * R[1] * R[2] );	//allocate memory
  		if (img == NULL) {
  			std::cout << "stim::image ERROR - failed to allocate memory for image" << std::endl;
  			exit(1);
  		}
  	}
  
  	void allocate(size_t x, size_t y, size_t c){	//allocate memory based on the resolution
  		R[0] = c; R[1] = x; R[2] = y;				//set the resolution
  		allocate();									//allocate memory
  	}
  
  	inline size_t idx(size_t x, size_t y, size_t c = 0) const {
  		return y * R[0] * R[1] + x * R[0] + c;
  	}
  
  	/// Returns the value for "white" based on the dynamic range (assumes white is 1.0 for floating point images)
  	T white(){
  		if (typeid(T) == typeid(double) || typeid(T) == typeid(float))
  			return (T)1.0;
  		else
  			return std::numeric_limits<T>::max();
  	}
  
  public:
  
  	size_t bytes() { return size() * sizeof(T); }
  
  	/// 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;
  		init();
  		allocate(I.X(), I.Y(), I.C());
  		memcpy(img, I.img, bytes());
  		return *this;
  	}
  
  	//determines if a filename represents a valid file format that can be loaded/saved
  	static bool test_filename(std::string f) {
  		stim::filename fname = f;
  		std::string ext = fname.extension();
  		if (ext == "bmp" ||
  			ext == "jpg" ||
  			ext == "png" ||
  			ext == "pbm" ||
  			ext == "tif" )
  			return true;
  		else
  			return false;
  	}
  
  	//save a Netpbm file
  	void load_netpbm(std::string filename) {
  		std::ifstream infile(filename.c_str(), std::ios::in | std::ios::binary);		//open an output file
  		if (!infile) {
  			std::cout << "Error opening input file in image::load_netpbm()" << std::endl;
  			exit(1);
  		}
  
  		size_t nc;													//allocate space for the number of channels
  		char format[2];												//allocate space to hold the image format tag
  		infile.read(format, 2);										//read the image format tag
  		infile.seekg(1, std::ios::cur);								//skip the newline character
  
  		if (format[0] != 'P') {
  			std::cout << "Error in image::load_netpbm() - file format tag is invalid: " << format[0] << format[1] << std::endl;
  			exit(1);
  		}
  		if (format[1] == '5') nc = 1;								//get the number of channels from the format flag
  		else if (format[1] == '6') nc = 3;
  		else {
  			std::cout << "Error in image::load_netpbm() - file format tag is invalid: " << format[0] << format[1] << std::endl;
  			exit(1);
  		}
  
  		unsigned char c;								//stores a character
  		while (infile.peek() == '#') {					//if the next character indicates the start of a comment
  			while (true) {
  				c = infile.get();
  				if (c == 0x0A) break;
  			}
  		}
  
  		std::string sw;									//create a string to store the width of the image
  		while(true){
  			c = infile.get();							//get a single character
  			if (c == ' ') break;						//exit if we've encountered a space
  			sw.push_back(c);							//push the character on to the string
  		}
  		size_t w = atoi(sw.c_str());					//convert the string into an integer
  
  		std::string sh;
  		while (true) {
  			c = infile.get();
  			if (c == 0x0A) break;
  			sh.push_back(c);
  		}
  
  		while (true) {									//skip the maximum value
  			c = infile.get();
  			if (c == 0x0A) break;
  		}
  		size_t h = atoi(sh.c_str());					//convert the string into an integer
  
  		allocate(w, h, nc);													//allocate space for the image
  		unsigned char* buffer = (unsigned char*)malloc(w * h * nc);			//create a buffer to store the read data
  		infile.read((char*)buffer, size());									//copy the binary data from the file to the image
  		infile.close();														//close the file
  		for (size_t n = 0; n < size(); n++) img[n] = (T)buffer[n];			//copy the buffer data into the image
  		free(buffer);														//free the buffer array
  	}
  	
  
  	//Copy N data points from source to dest, casting while doing so
  	template<typename S, typename D>
  	void type_copy(S* source, D* dest, size_t N) {
  		if (typeid(S) == typeid(D))						//if both types are the same
  			memcpy(dest, source, N * sizeof(S));		//just use a memcpy
  		for (size_t n = 0; n < N; n++)					//otherwise, iterate through each element
  			dest[n] = (D)source[n];							//copy and cast
  	}
  	/// Load an image from a file
  	void load(std::string filename){
  		//Use CImg to load the file
  		cimg_library::CImg<T> cimg(filename.c_str());	//create a CImg object for the image file
  
  		set_noninterleaved(cimg.data(), cimg.width(), cimg.height(), cimg.spectrum());
  	}
  
  
  
  	//save a Netpbm file
  	void save_netpbm(std::string filename) {
  		std::ofstream outfile(filename.c_str(), std::ios::out | std::ios::binary);		//open an output file
  		if(!outfile) {
  			std::cout << "Error generating output file in image::save_netpbm()" << std::endl;
  			exit(1);
  		}
  		if (sizeof(T) != 1) {
  			std::cout << "Error in image::save_netpbm() - data type must be 8-bit integer." << std::endl;
  			exit(1);
  		}
  		std::string format;
  		if (channels() == 1) outfile << "P5" << (char)0x0A;			//output P5 if the file is grayscale
  		else if (channels() == 3) outfile << "P6" << (char)0x0A;		//output P6 if the file is color
  		else {
  			std::cout << "Error in image::save_netpbm() - data must be grayscale or RGB." << std::endl;
  			exit(1);
  		}
  		size_t w = width();
  		size_t h = height();
  		outfile << w << " " << h << (char)0x0A;			//save the width and height
  		outfile << "255" << (char)0x0A;								//output the maximum value
  		outfile.write((const char*)img, size());			//write the binary data
  		outfile.close();
  	}
  
  	//save a file
  	void save(std::string filename){
  		stim::filename file(filename);
  		if (file.extension() == "raw" || file.extension() == "") {
  			std::ofstream outfile(filename.c_str(), std::ios::binary);
  			outfile.write((char*)img, sizeof(T) * R[0] * R[1] * R[2]);
  			outfile.close();
  		}
  		else {
  			cimg_library::CImg<T> cimg((unsigned int)R[1], (unsigned int)R[2], 1, (unsigned int)R[0]);
  			get_noninterleaved(cimg.data());
  			cimg.save(filename.c_str());
  		}
  	}
  
  	/// Returns an image cast to the specified format
  	template<typename U>
  	image<U> convert() {
  		
  		image<U> new_image(R[1], R[2], R[0]);					//create a new image with the destination data type
  
  		size_t ni = R[0] * R[1] * R[2];							//calculate the number of data points in the image
  
  		double inmax = (std::numeric_limits<T>::max)();				//get the maximum value for the input image
  		double outmax = (std::numeric_limits<U>::max)();				//get the maximum value for the output image
  		for (size_t i = 0; i < ni; i++) {							//for each pixel in the image
  			if (img[i] > outmax) new_image(i) = outmax;			//if the source pixel is greater than the maximum destination pixel, set the output to maximum
  			else new_image(i) = img[i];							//otherwise, copy the source value and cast it to the destination value		
  		}
  		return new_image;
  	}
  
  	void set_interleaved(T* buffer, size_t width, size_t height, size_t channels){
  		allocate(width, height, channels);
  		memcpy(img, buffer, bytes());
  	}
  
  	//create an image from an interleaved buffer
  	void set_interleaved_rgb(T* buffer, size_t width, size_t height){
  		set_interleaved(buffer, width, height, 3);
  	}
  
  	void set_interleaved_bgr(T* buffer, size_t width, size_t height){
  		allocate(width, height, 3);
  		T value;
  		size_t i;
  		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++){
  					i = y * X() * C() + x * C() + (2-c);
  					value = buffer[i];
  					img[idx(x, y, c)] = value;
  				}
  			}
  		}
  	}
  
  	void set_interleaved(T* buffer, size_t width, size_t height){
  		set_interleaved_rgb(buffer, width, height);
  	}
  
  	//copies data in the given channel order as a non-interleaved image
  	void set_noninterleaved(T* data, size_t width, size_t height, size_t chan) {
  		allocate(width, height, chan);
  
  		//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++) {
  					img[idx(x, y, c)] = data[c * Y() * X() + y * X() + x];
  				}
  			}
  		}
  	}
  
  	void get_interleaved_bgr(T* data){
  		//for each channel
  		T* source;
  		if (C() == 3) {
  			source = img;														//if the image has 3 channels, interleave all three
  		}
  		else if (C() == 4) {
  			source = img + X() * Y();											//if the image has 4 channels, skip the alpha channel
  		}
  		else {
  			throw std::runtime_error("ERROR: a BGR image must be 3 or 4 channels");	//throw an error if any other number of channels is provided
  		}
  		for(size_t y = 0; y < Y(); y++){
  			for(size_t x = 0; x < X(); x++){
  				for(size_t c = 0; c < 3; c++){
  					data[y * X() * 3 + x * 3 + (2-c)] = source[y*3*R[1] + x*3 + 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)];
  				}
  			}
  		}
  	}
  
  	/// Return an image representing a specified channel
  	/// @param c is the channel to be returned
  	image<T> channel(size_t c) const {		
  		image<T> r(X(), Y(), 1);				//create a new image
  		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;
  	}
  
  	/// Returns an std::vector containing each channel as a separate image
  	std::vector< image<T> > split() const {
  		std::vector< image<T> > r;			//create an image array
  		r.resize(C());						//create images for each channel
  
  		for (size_t c = 0; c < C(); c++) {	//for each channel
  			r[c] = channel(c);				//copy the channel image to the array
  		}
  		return r;
  	}
  
  	/// Merge a series of single-channel images into a multi-channel image
  	void merge(std::vector< image<T> >& list) {
  		size_t x = list[0].width();				//calculate the size of the image
  		size_t y = list[0].height();
  		allocate(x, y, list.size());			//re-allocate the image
  		for (size_t c = 0; c < list.size(); c++)		//for each channel
  			set_channel(list[c].channel(0).data(), c);	//insert the channel into the output image
  	}
  
  	T& operator()(size_t x, size_t y, size_t c = 0){
  		return img[idx(x, y, c)];
  	}
  
  	/// This function returns a pixel reference based on a 1D index into the image
  	T& operator()(size_t i) {
  		return img[i];
  	}
  
  	/// Set all elements in the image to a given scalar value
  
  	/// @param v is the value used to set all values in the image
  	void set_all(T v) {														//set all elements of the image to a given value v
  		size_t N = size();
  		for (size_t n = 0; n < N; n++) img[n] = v;
  	}
  	image<T> operator=(T v){
  		set_all(v);
  		return *this;
  	}
  
  	/// invert the image, given a specified maximum value (ex. maxval = 255, I' = 255 - I)
  	/*image<T> invert(T maxval) {
  		image<T> result(width(), height(), channels());		//create a new image
  		size_t N = size();									//get the number of elements in the image
  		for (size_t n = 0; n < N; n++)
  			result.data()[n] = maxval - img[n];				//perform the inversion and save the result to the new image
  		return result;
  	}*/
  
  	/// Stretch the contrast of the image such that the minimum and maximum intensity match the given values
  	image<T> stretch(T low, T high) {
  		T maxval = maxv();
  		T minval = minv();
  		image<T> result = *this;				//create a new image for output
  		if (maxval == minval) {					//if the minimum and maximum values are the same, return an image composed of low
  			result = low;
  			return result;
  		}	
  		
  		size_t N = size();						//get the number of values in the image
  		T range = maxval - minval;			//calculate the current range of the image
  		T desired_range = high - low;		//calculate the desired range of the image
  		for (size_t n = 0; n < N; n++) {		//for each element in the image
  			result.data()[n] = desired_range * (img[n] - minval) / range + low;
  		}
  		return result;
  	}
  
  	/// Add a border of width w with the given value around the image
  	/// @param w specifies the total size of the border
  	/// @param T is the pixel value (all channels will be the same)
  	image<T> border(size_t w, T value = 0) {
  		image<T> result(width() + w * 2, height() + w * 2, channels());						//create an output image
  		result = value;														//assign the border value to all pixels in the new image
  		for (size_t y = 0; y < height(); y++) {								//for each pixel in the original image
  			for (size_t x = 0; x < width(); x++) {
  				size_t n = (y + w) * (width() + w * 2) + x + w;				//calculate the index of the corresponding pixel in the result image
  				size_t n0 = idx(x,y);										//calculate the index for this pixel in the original image
  				result.data()[n] = img[n0];									// copy the original image to the result image afer the border area
  			}
  		}
  		return result;
  	}
  
  	/// Adds curcular padding for the specified number of pixels - in this case replicating the boundary pixels
  	image<T> pad_replicate(size_t p) {
  		image<T> result(width() + p * 2, height() + p * 2, channels());						//create an output image
  		result = 0;
  		//result = value;														//assign the border value to all pixels in the new image
  		for (size_t y = 0; y < height(); y++) {								//for each pixel in the original image
  			for (size_t x = 0; x < width(); x++) {
  				size_t n = (y + p) * (width() + p * 2) + x + p;				//calculate the index of the corresponding pixel in the result image
  				size_t n0 = idx(x, y);										//calculate the index for this pixel in the original image
  				result.data()[n] = img[n0];									// copy the original image to the result image afer the border area
  			}
  		}
  		size_t l = p;
  		size_t r = p + width() - 1;
  		size_t t = p;
  		size_t b = p + height() - 1;
  		for (size_t y = 0; y < p; y++) for (size_t x = l; x <= r; x++) result(x, y) = result(x, t);						//pad the top
  		for (size_t y = b + 1; y < result.height(); y++) for (size_t x = l; x <= r; x++) result(x, y) = result(x, b);	//pad the bottom
  		for (size_t y = t; y <= b; y++) for (size_t x = 0; x < l; x++) result(x, y) = result(l, y);						//pad the left
  		for (size_t y = t; y <= b; y++) for (size_t x = r+1; x < result.width(); x++) result(x, y) = result(r, y);		//pad the right
  		for (size_t y = 0; y < t; y++) for (size_t x = 0; x < l; x++) result(x, y) = result(l, t);						//pad the top left
  		for (size_t y = 0; y < t; y++) for (size_t x = r+1; x < result.width(); x++) result(x, y) = result(r, t);		//pad the top right
  		for (size_t y = b+1; y < result.height(); y++) for (size_t x = 0; x < l; x++) result(x, y) = result(l, b);		//pad the bottom left
  		for (size_t y = b+1; y < result.height(); y++) for (size_t x = r + 1; x < result.width(); x++) result(x, y) = result(r, b);		//pad the bottom right
  		return result;
  	}
  
  	/// 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[y * X() + x];
  			}
  		}
  	}
  
  	/// Set the specified channel to a constant value
  
  	/// @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 val, size_t c) {
  		size_t x, y;
  		for (y = 0; y < Y(); y++) {
  			for (x = 0; x < X(); x++) {
  				img[idx(x, y, c)] = val;
  			}
  		}
  	}
  
  	size_t channels() const{
  		return C();
  	}
  
  	size_t width() const{
  		return X();
  	}
  
  	size_t height() const{
  		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
  	}
  
  	image<T> crop(size_t x0, size_t y0, size_t w, size_t h){
  		image<T> result(w, h, C());								//create the output cropped image
  
  		size_t srci;
  		size_t dsti;
  		size_t line_bytes = w * C();							//calculate the number of bytes in a line
  		for (size_t yi = 0; yi < h; yi++) {						//for each row in the cropped image
  			srci = (y0 + yi) * X() * C() + x0 * C();			//calculate the source index
  			dsti = yi * w * C();								//calculate the destination index
  			memcpy(&result.img[dsti], &img[srci], line_bytes);	//copy the data
  		}
  		return result;
  	}
  
  	//crop regions given by an array of 1D index values
  	std::vector< image<T> > crop_idx(size_t w, size_t h, std::vector<size_t> idx) {
  		std::vector< image<T> > result(idx.size());										//create an array of image files to return
  		for (size_t i = 0; i < idx.size(); i++) {										//for each specified index point
  			size_t y = idx[i] / X();													//calculate the y coordinate from the 1D index (center of ROI)
  			size_t x = idx[i] - y * X();												//calculate the x coordinate (center of ROI)
  			y -= w / 2;																	//update x and y values to reflect the lower corner of the ROI
  			x -= h / 2;
  			result[i] = crop(x, y, w, h);												//get the cropped image and store it in the result array
  		}
  		return result;
  	}
  
  	//operator functions
  	image<T> operator+(image<T> rhs) {
  		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] = img[n] + rhs.img[n];		//perform the inversion
  		return r;								//return the inverted image
  	}
  
  	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);
  	}
  
  	/// Casting operator, casts every value in an image to a different data type V
  	template<typename V>
  	operator image<V>() {
  		image<V> r(X(), Y(), C());					//create a new image
  		std::copy(img, img + size(), r.data());		//copy and cast the data
  		return r;									//return the new image
  	}
  
  };
  
  };		//end namespace stim
  
  
  #endif