#ifndef RTSFUNCTION3D_H
#define RTSFUNCTION3D_H

#define DIST_MAX	255

#include "rtsLinearAlgebra.h"
//#include "rtsDTGrid3D.h"
#include <fstream>
//#include <iostream>
//#include <math.h>
//#include <queue>
//#include <algorithm>
using namespace std;

typedef int indextype;

///This class represents a 3D implicit function as a grid.  It provides methods for accessing values, interpolation, and several utilities.

template <class T> class rtsFunction3D
{
private:
	//pointer to store the data
	T* m_data;
	//resolution of the data (x, y, z) dimensional extents
	vector3D<indextype> m_resolution;
	T m_boundary;			//boundary condition
	point3D<double> m_domain_min;	//min and max range values (used for parametric access)
	point3D<double> m_domain_max;
	vector3D<double> m_voxel_size;
	
	//bit-blit function copies 3D data quickly from source to dest
	void blit3D(const T* source,
				   indextype s_px, indextype s_py, indextype s_pz,
				   indextype s_sx, indextype s_sy, indextype s_sz,
				   T* dest,
				   indextype d_px, indextype d_py, indextype d_pz,
				   indextype d_sx, indextype d_sy, indextype d_sz,
				   indextype blit_size_x, indextype blit_size_y, indextype blit_size_z);

	void shallow_copy(const rtsFunction3D<T> source, rtsFunction3D<T> &dest);
	inline point3D<double> getParameter(indextype i);
	void initialize_empty(indextype res_x, indextype res_y, indextype res_z);

public:
	//construct an implicit function with a size of 1
	rtsFunction3D();			///<Create an empty implicit function
	//construct an implicit function of the specified resolution
	rtsFunction3D(indextype res_x, indextype res_y, indextype res_z);	///<Create an implicit function with the specified resolution
	//construct an implicit function from sample data and a specified size
	rtsFunction3D(T* data, indextype res_x, indextype res_y, indextype res_z);	///<Create an implicit function from previous data at the specified resolution
	//shallow-copy constructor, defines all shallow variables
	rtsFunction3D(vector3D<int> resolution, T boundary, point3D<double> min_domain, point3D<double> max_domain);
	//full copy constructor, defines all variables
	rtsFunction3D(T* data, vector3D<int> resolution, T boundary, point3D<double> min_domain, point3D<double> max_domain);
	rtsFunction3D(const rtsFunction3D<T> &original);	//copy constructor
	~rtsFunction3D();		//destructor
	void Init(indextype res_x, indextype res_y, indextype res_z);

	//overloaded operators
	rtsFunction3D<T>& operator=(const rtsFunction3D<T>& original);		///<Overloaded operator creates a copy of an implicit function
	rtsFunction3D<T>& operator=(const T constant);						///<Overloaded operator sets all points in an implicit function to the given constant value
	inline T& operator()(indextype x, indextype y, indextype z);		///<Allows access to the sample point indexed by x, y, and z using the parenthesis operator
	inline T operator()(double i, double j, double k);					///<Allows access to the implicit function (based on the domain boundaries) at the position (i, j, k).  This class uses linear interpolation.
	rtsFunction3D<T>& operator*=(const T constant);						///<Multiplies the values at all sample points by a constant.
	rtsFunction3D<T>& operator+=(const T constant);						///<Adds a constant to the values at all sample points.
	rtsFunction3D<T>& operator-=(const T constant);						///<Subtracts a constant from the values at all sample points.
	rtsFunction3D<T>& operator/=(const T constant);						///<Divides all values by a constant.
	const rtsFunction3D<T> operator+(const T constant);					///<Adds a constant to an implicit function and returns a new function.
	const rtsFunction3D<T> operator-(const T constant);					///<Subtracts a constant from an implicit function and returns a new function.
	const rtsFunction3D<T> operator*(const T constant);					///<Multiplies an implicit function by a constant and returns a new function.
	const rtsFunction3D<T> operator/(const T constant);					///<Divides an implicit function by a constant and returns a new function.

	//casting operator
	//template <class U> friend class rtsFunction3D<U>;
	template <class U> operator rtsFunction3D<U>();						///<Casts between data types.

	//friend classes for overloading "backwards" operations (like 3*function)
	friend rtsFunction3D<T> operator*(const T lhs, rtsFunction3D<T> rhs){return rhs*lhs;}	///<Allows associative multiplication.
	friend rtsFunction3D<T> operator+(const T lhs, rtsFunction3D<T> rhs){return rhs+lhs;}	///<Allows associative addition.
	friend rtsFunction3D<T> operator-(const T lhs, rtsFunction3D<T> rhs)					///<Allows associative subtraction.
	{
		rtsFunction3D<T> result;
		rhs.shallow_copy(rhs, result);	//make a copy of all of the shallow variables and allocate memory
		indextype size = rhs.m_resolution.x * rhs.m_resolution.y * rhs.m_resolution.z;
		//iterate and subtract
		for(indextype i=0; i<size; i++)
			result.m_data[i] = lhs - rhs.m_data[i];

		return result;
	}
	//friend rtsFunction3D<T> operator/(const T lhs, rtsFunction3D<T> rhs);

	//loading/saving data to disk
	void LoadRAW(indextype header_size, indextype data_x, indextype data_y, indextype data_z, const char* filename);	///<Loads RAW data from a file with the specified header size and data size.
	void SaveRAW(const char* filename);							///<Save the data as RAW data to disk.
	void LoadVOL(const char* filename);							///<Load a VOL file from disk.
	void SaveVOL(const char* filename);							///<Save a VOL file to disk.

	rtsFunction3D<T> Project2D();			//<Projects the data along the z-axis using a maximum-intensity projection.

	//data access methods
	inline T& xyz(indextype x, indextype y, indextype z);
	inline T ijk(double i, double j, double k);
	void Parameterize(double x_min, double x_max, double y_min, double y_max, double z_min, double z_max);
	void setBoundary(T boundary){m_boundary = boundary;}
	T getBoundary(){return m_boundary;}
	T* GetBits();
	indextype DimX(){return m_resolution.x;}
	indextype DimY(){return m_resolution.y;}
	indextype DimZ(){return m_resolution.z;}
	inline point3D<double> getParameter(indextype x, indextype y, indextype z);
	inline point3D<indextype> getNearestIndex(double i, double j, double k);
	inline point3D<double> getFractionalIndex(double i, double j, double k);
	inline point3D<indextype> getNearestIndex(indextype i);
	point3D<double> getMinDomain(){return m_domain_min;}
	point3D<double> getMaxDomain(){return m_domain_max;}

	//data input methods
	void Insert(rtsFunction3D<T>* source, indextype x, indextype y, indextype z);

	//data massaging
	
	void Scale(T min, T max);
	void Crop(indextype x, indextype y, indextype z, indextype size_x, indextype size_y, indextype size_z);
	void Binary(T threshold, T true_value);		///<Turns the image into a binary image based on a threshold value T.  All values below T are set to 0, all values above are set to true_value.
	void Threshold(T min, T value);
	void Threshold(T min, T max, T value);
	void Threshold(T min, T max, T inside, T outside);
	void ClampMax(T max);		///<Clamps the function to the given maximum value.
	void ClampMin(T min);
	T getMin();
	T getMax();

	//create new data
	rtsFunction3D<T>* Resample(indextype newres_x, indextype newres_y, indextype newres_z);

	//output functions
	void toConsole();

};

template <class T>
void rtsFunction3D<T>::blit3D(const T* source,
				   indextype s_px, indextype s_py, indextype s_pz,
				   indextype s_sx, indextype s_sy, indextype s_sz,
				   T* dest,
				   indextype d_px, indextype d_py, indextype d_pz,
				   indextype d_sx, indextype d_sy, indextype d_sz,
				   indextype blit_size_x, indextype blit_size_y, indextype blit_size_z)
{
	indextype ps, pd;		//stores the mapping for the source point to the dest point
	//find the maximum points that can be blit to (in case source overlaps the edges of dest)
	blit_size_x = min(blit_size_x, min(s_sx - s_px, d_sx - d_px));
	blit_size_y = min(blit_size_y, min(s_sy - s_py, d_sy - d_py));
	blit_size_z = min(blit_size_z, min(s_sz - s_pz, d_sz - d_pz));

	indextype source_z_offset = s_sx * s_sy;
	indextype dest_z_offset = d_sx * d_sy;

	indextype z,y;
	for(z=0; z<blit_size_z; z++)
		for(y=0; y<blit_size_y; y++)
		{
			ps = (z + s_pz) * source_z_offset + (y + s_py) * s_sx + s_px;
			pd = (z + d_pz) * dest_z_offset + (y + d_py) * d_sx + d_px;
			memcpy((void*)(&dest[pd]), (void*)(&source[ps]), sizeof(T)*blit_size_x);
		}
}

template <class T>
void rtsFunction3D<T>::shallow_copy(const rtsFunction3D<T> source, rtsFunction3D<T> &dest)
{
	dest = rtsFunction3D<T>(source.m_resolution.x, source.m_resolution.y, source.m_resolution.z);
	dest.m_boundary = source.m_boundary;
	dest.m_domain_max = source.m_domain_max;
	dest.m_domain_min = source.m_domain_max;
	dest.m_voxel_size = source.m_voxel_size;
}

template <class T>
rtsFunction3D<T>::rtsFunction3D(vector3D<int> resolution, T boundary, point3D<double> domain_min, point3D<double> domain_max)
{
	//This function creates an implicit function based on all of the shallow variables
	m_resolution = resolution;
	m_boundary = boundary;
	m_domain_min = domain_min;
	m_domain_max = domain_max;
	m_voxel_size = domain_max - domain_min;
	m_voxel_size.x /= m_resolution.x;
	m_voxel_size.y /= m_resolution.y;
	m_voxel_size.z /= m_resolution.z;

	//allocate the data
	m_data = new T[m_resolution.x * m_resolution.y * m_resolution.z];
}

template <class T>
rtsFunction3D<T>::rtsFunction3D(T* data, vector3D<int> resolution, T boundary, point3D<double> domain_min, point3D<double> domain_max)
{
	//This function creates an implicit function based on ALL of the variables
	m_resolution = resolution;
	m_boundary = boundary;
	m_domain_min = domain_min;
	m_domain_max = domain_max;
	m_voxel_size = domain_max - domain_min;
	m_voxel_size.x /= m_resolution.x;
	m_voxel_size.y /= m_resolution.y;
	m_voxel_size.z /= m_resolution.z;

	//allocate the data
	indextype size = m_resolution.x * m_resolution.y * m_resolution.z;
	m_data = new T[size];
	memcpy(m_data, data, sizeof(T)*size);
	//for(int i=0; i<size; i++)
	//	m_data[i] = data[i];
}



template <class T>
rtsFunction3D<T>::rtsFunction3D()
{
	m_resolution.x = 1;
	m_resolution.y = 1;
	m_resolution.z = 1;
	m_data = new T[1];
	memset(&m_boundary, 0, sizeof(T));							//initialize boundary condition
	m_domain_min = point3D<double>(0.0, 0.0, 0.0);	//set range parameters
	m_domain_max = point3D<double>(1.0, 1.0, 1.0);
	m_voxel_size = vector3D<double>(1.0, 1.0, 1.0);
}

template <class T>
void rtsFunction3D<T>::initialize_empty(indextype res_x, indextype res_y, indextype res_z)
{
	m_resolution.x = res_x;					//set resolution vector
	m_resolution.y = res_y;
	m_resolution.z = res_z;
	m_data = (T*)calloc(res_x*res_y*res_z, sizeof(T));		//allocate data
}

template <class T>
rtsFunction3D<T>::rtsFunction3D(indextype res_x, indextype res_y, indextype res_z)
{
	initialize_empty(res_x, res_y, res_z);
	memset(&m_boundary, 0, sizeof(T));							//initialize boundary condition
	m_domain_min = point3D<double>(0.0, 0.0, 0.0);	//set range parameters
	m_domain_max = point3D<double>(1.0, 1.0, 1.0);

	m_voxel_size = m_domain_max - m_domain_min;
	m_voxel_size.x /= m_resolution.x;
	m_voxel_size.y /= m_resolution.y;
	m_voxel_size.z /= m_resolution.z;
}

template <class T>
void rtsFunction3D<T>::Init(indextype res_x, indextype res_y, indextype res_z)
{
	initialize_empty(res_x, res_y, res_z);
	m_domain_min = point3D<double>(0.0, 0.0, 0.0);	//set range parameters
	m_domain_max = point3D<double>(1.0, 1.0, 1.0);
	memset(&m_boundary, 0, sizeof(T));							//initialize boundary condition

	m_voxel_size = m_domain_max - m_domain_min;
	m_voxel_size.x /= m_resolution.x;
	m_voxel_size.y /= m_resolution.y;
	m_voxel_size.z /= m_resolution.z;
}

template <class T>
rtsFunction3D<T>::rtsFunction3D(T* data, indextype res_x, indextype res_y, indextype res_z)
{
	m_resolution.x = res_x;					//set resolution vector
	m_resolution.y = res_y;
	m_resolution.z = res_z;
	m_data = new T[res_x*res_y*res_z];		//allocate data
	//copy the sample data into the data array
	indextype size = res_x*res_y*res_z;
	for(indextype i=0; i<size; i++)
		m_data[i] = data[i];
	memset(&m_boundary, 0, sizeof(T));							//initialize boundary condition
	m_domain_min = point3D<double>(0.0, 0.0, 0.0);	//set range parameters
	m_domain_max = point3D<double>(1.0, 1.0, 1.0);

	m_voxel_size = domain_max - domain_min;
	m_voxel_size.x /= m_resolution.x;
	m_voxel_size.y /= m_resolution.y;
	m_voxel_size.z /= m_resolution.z;
}

template <class T>
rtsFunction3D<T>::rtsFunction3D(const rtsFunction3D<T>& original)
{
	//copy the shallow variables
	m_resolution = original.m_resolution;
	m_boundary = original.m_boundary;
	m_domain_min = original.m_domain_min;
	m_domain_max = original.m_domain_max;
	m_voxel_size = original.m_voxel_size;

	//allocate space for the data
	m_data = new T[m_resolution.x * m_resolution.y * m_resolution.z];
	//copy the data
	blit3D(original.m_data,
		   0, 0, 0,
		   m_resolution.x, m_resolution.y, m_resolution.z,
		   m_data,
		   0, 0, 0,
		   m_resolution.x, m_resolution.y, m_resolution.z,
		   m_resolution.x, m_resolution.y, m_resolution.z);
}

template <class T>
rtsFunction3D<T>::~rtsFunction3D()
{
	delete m_data;
}

template <class T>
typename rtsFunction3D<T>& rtsFunction3D<T>::operator=(const T rhs)
{
	indextype size = m_resolution.x*m_resolution.y*m_resolution.z;
	for(int i=0; i<size; i++)
		m_data[i] = rhs;

	return *this;
}

template <class T>
typename rtsFunction3D<T>& rtsFunction3D<T>::operator=(const rtsFunction3D<T>& rhs)
{
	//check for self-assignment
	if(this == &rhs)
		return *this;

	//deallocate memory
	if(m_data != NULL)
		delete m_data;

	//copy the shallow variables
	m_resolution = rhs.m_resolution;
	m_boundary = rhs.m_boundary;
	m_domain_min = rhs.m_domain_min;
	m_domain_max = rhs.m_domain_max;
	m_voxel_size = rhs.m_voxel_size;

	//allocate and copy memory
	m_data = new T[m_resolution.x * m_resolution.y * m_resolution.z];
	//copy the data
	blit3D(rhs.m_data,
		   0,0,0,
		   m_resolution.x, m_resolution.y, m_resolution.z,
		   m_data, 
		   0, 0, 0, 
		   m_resolution.x, m_resolution.y, m_resolution.z,
		   m_resolution.x, m_resolution.y, m_resolution.z);

	//return the left hand side
	return *this;
}

template <class T>
inline T& rtsFunction3D<T>::operator ()(indextype x, indextype y, indextype z)
{
	return xyz(x, y, z);
}

template <class T>
inline T rtsFunction3D<T>::operator()(double i, double j, double k)
{
	return ijk(i, j, k);
}

template <class T>
rtsFunction3D<T>& rtsFunction3D<T>::operator *=(const T constant)
{
	indextype size = m_resolution.x * m_resolution.y * m_resolution.z;
	for(indextype i = 0; i<size; i++)
		m_data[i] *= constant;

	return *this;
}

template <class T>
rtsFunction3D<T>& rtsFunction3D<T>::operator +=(const T constant)
{
	indextype size = m_resolution.x * m_resolution.y * m_resolution.z;
	for(indextype i = 0; i<size; i++)
		m_data[i] += constant;

	return *this;
}
template <class T>
rtsFunction3D<T>& rtsFunction3D<T>::operator -=(const T constant)
{
	indextype size = m_resolution.x * m_resolution.y * m_resolution.z;
	for(indextype i = 0; i<size; i++)
		m_data[i] -= constant;

	return *this;
}
template <class T>
rtsFunction3D<T>& rtsFunction3D<T>::operator /=(const T constant)
{
	indextype size = m_resolution.x * m_resolution.y * m_resolution.z;
	for(indextype i = 0; i<size; i++)
		m_data[i] /= constant;

	return *this;
}
template <class T>
const rtsFunction3D<T> rtsFunction3D<T>::operator *(const T constant)
{
	rtsFunction3D<T> result = (*this);
	result *= constant;

	return result;
}

template <class T>
const rtsFunction3D<T> rtsFunction3D<T>::operator +(const T constant)
{
	rtsFunction3D<T> result = (*this);
	result += constant;

	return result;
}

template <class T>
const rtsFunction3D<T> rtsFunction3D<T>::operator -(const T constant)
{
	rtsFunction3D<T> result = (*this);
	result -= constant;

	return result;
}

template <class T>
const rtsFunction3D<T> rtsFunction3D<T>::operator /(const T constant)
{
	rtsFunction3D<T> result = (*this);
	result /= constant;

	return result;
}

template <class T>
template <class U>
rtsFunction3D<T>::operator rtsFunction3D<U>()
{
	//cast one type to another
	//create the data pointer from the current function
	indextype size = m_resolution.x * m_resolution.y * m_resolution.z;
	U* new_data = new U[size];
	for(int i=0; i<size; i++)
		new_data[i] = m_data[i];
	rtsFunction3D<U> cast_result(new_data, m_resolution, m_boundary, m_domain_min, m_domain_max);

	return cast_result;
}

template <class T>
inline T& rtsFunction3D<T>::xyz(indextype x, indextype y, indextype z)
{
	if(x<0 || y<0 || z<0 || x>=m_resolution.x || y>=m_resolution.y || z>=m_resolution.z)
		return m_boundary;
	//return m_data[(z * m_resolution.x * m_resolution.y) + (y * m_resolution.x) + x];
	return m_data[x + m_resolution.x * (y + z * m_resolution.y)];

}

template <class T>
inline point3D<indextype> rtsFunction3D<T>::getNearestIndex(indextype i)
{
	point3D<indextype> result;
	result.z = i/(m_resolution.x*m_resolution.y);
	indextype mod = i%(m_resolution.x*m_resolution.y);
	result.y = mod/m_resolution.x;
	result.x = mod%m_resolution.x;

	return result;

}

template <class T>
void rtsFunction3D<T>::LoadRAW(indextype header_size, indextype size_x,
							   indextype size_y, indextype size_z, const char *filename)
{
	//set the data size
	m_resolution = vector3D<indextype>(size_x, size_y, size_z);
	//delete any previous data
	if(m_data != NULL)
	{
		delete m_data;
		m_data = NULL;
	}

	ifstream infile(filename, ios::in | ios::binary);

	//load the header
	unsigned char* header = new unsigned char[header_size];
	infile.read((char*)header, header_size);

	//load the actual data
	indextype size = m_resolution.x * m_resolution.y * m_resolution.z;
	//m_data = (T*)malloc(size*sizeof(T));
	initialize_empty(m_resolution.x, m_resolution.y, m_resolution.z);
	infile.read((char*)m_data, size*sizeof(T));

	//calculate min and maxes
	infile.close();
}

template <class T>
void rtsFunction3D<T>::LoadVOL(const char *filename)
{
	ifstream infile(filename, ios::in | ios::binary);	//create the files stream
	if(!infile)
		return;

	indextype size_x, size_y, size_z;				//create variables to store the size of the data set
	//load the dimensions of the data set
	infile.read((char*)&size_x, sizeof(int));			//load the file header
	infile.read((char*)&size_y, sizeof(int));
	infile.read((char*)&size_z, sizeof(int));

	//close the file
	infile.close();
	//load the raw data
	LoadRAW(12, size_x, size_y, size_z, filename);
}

template <class T>
void rtsFunction3D<T>::SaveVOL(const char *filename)
{
	ofstream outfile(filename, ios::out | ios::binary);	//create the binary file stream

	//write the volume size to the file
	vector3D<int> vol_size = m_resolution;
	outfile.write((char*)&vol_size.x, sizeof(int));
	outfile.write((char*)&vol_size.y, sizeof(int));
	outfile.write((char*)&vol_size.z, sizeof(int));

	outfile.write((char*)m_data, sizeof(T)*vol_size.x*vol_size.y*vol_size.z);
}

template <class T>
void rtsFunction3D<T>::SaveRAW(const char *filename)
{
	ofstream outfile(filename, ios::out | ios::binary);	//create the binary file stream

	//write the volume data
	outfile.write((char*)m_data, sizeof(T)*m_resolution.x*m_resolution.y*m_resolution.z);
}

template <class T>
inline T rtsFunction3D<T>::ijk(double i, double j, double k)
{
	/*This function determines the value at the specified parametric points
	defined by the m_domain_min and m_domain_max parameter values.*/

	//if the parameter is outside the range, return the boundary value
	if(i<m_domain_min.x || j<m_domain_min.y || k<m_domain_min.z ||
	   i>m_domain_max.x || j>m_domain_max.y || k>m_domain_max.z)
	   return m_boundary;
	
	point3D<double> index = getFractionalIndex(i, j, k);

	//cout<<index.x<<","<<index.y<<","<<index.z<<endl;

	//interpolate the values
	int f_x = (int)floor(index.x);				//calculate floor and ceiling values
	int f_y = (int)floor(index.y);
	int f_z = (int)floor(index.z);
	int c_x = (int)ceil(index.x);
	int c_y = (int)ceil(index.y);
	int c_z = (int)ceil(index.z);

	double x_d = index.x - f_x;			//find the point within the voxel
	double y_d = index.y - f_y;
	double z_d = index.z - f_z;

	T i_1 = xyz(f_x, f_y, f_z)*(1.0 - z_d) + xyz(f_x, f_y, c_z)*(z_d);	//interpolate along z
	T i_2 = xyz(f_x, c_y, f_z)*(1.0 - z_d) + xyz(f_x, c_y, c_z)*(z_d);
	T j_1 = xyz(c_x, f_y, f_z)*(1.0 - z_d) + xyz(c_x, f_y, c_z)*(z_d);
	T j_2 = xyz(c_x, c_y, f_z)*(1.0 - z_d) + xyz(c_x, c_y, c_z)*(z_d);

	T w_1 = i_1*(1.0 - y_d) + i_2*(y_d);
	T w_2 = j_1*(1.0 - y_d) + j_2*(y_d);

	return w_1*(1.0 - x_d) + w_2*(x_d);
}


template <class T>
void rtsFunction3D<T>::Parameterize(double x_min, double x_max, double y_min, double y_max, double z_min, double z_max)
{
	m_domain_min = point3D<double>(x_min, y_min, z_min);
	m_domain_max = point3D<double>(x_max, y_max, z_max);
	m_voxel_size = m_domain_max - m_domain_min;
	m_voxel_size.x /= m_resolution.x;
	m_voxel_size.y /= m_resolution.y;
	m_voxel_size.z /= m_resolution.z;
}

template <class T>
inline point3D<double> rtsFunction3D<T>::getParameter(indextype x, indextype y, indextype z)
{
	//get the value between 0 and 1
	point3D<double> normalized((double)x / (double)(m_resolution.x) + (1.0/(m_resolution.x*2.0)),
							   (double)y / (double)(m_resolution.y) + (1.0/(m_resolution.y*2.0)),
							   (double)z/(double)(m_resolution.z) + (1.0/(m_resolution.z*2.0)));

	point3D<double> result(normalized.x * (m_domain_max.x - m_domain_min.x) + m_domain_min.x,
						   normalized.y * (m_domain_max.y - m_domain_min.y) + m_domain_min.y,
						   normalized.z * (m_domain_max.z - m_domain_min.z) + m_domain_min.z);

	return result;
}

template <class T>
inline point3D<indextype> rtsFunction3D<T>::getNearestIndex(double i, double j, double k)
{
	//this function returns the index of the voxel containing the specified parameter point
	point3D<double> normalized((i - m_domain_min.x)/(m_domain_max.x-m_domain_min.x),
							   (j - m_domain_min.y)/(m_domain_max.y-m_domain_min.y),
							   (k - m_domain_min.z)/(m_domain_max.z-m_domain_min.z));
	
	point3D<indextype> result((normalized.x - (1.0/(m_resolution.x*2.0)))*(double)m_resolution.x+0.5,
							  (normalized.y - (1.0/(m_resolution.y*2.0)))*(double)m_resolution.y+0.5,
							  (normalized.z - (1.0/(m_resolution.z*2.0)))*(double)m_resolution.z+0.5);

	return result;
}

template <class T>
inline point3D<double> rtsFunction3D<T>::getFractionalIndex(double i, double j, double k)
{
	//this function returns the index of the voxel containing the specified parameter point
	point3D<double> normalized((i - m_domain_min.x)/(m_domain_max.x-m_domain_min.x),
							   (j - m_domain_min.y)/(m_domain_max.y-m_domain_min.y),
							   (k - m_domain_min.z)/(m_domain_max.z-m_domain_min.z));
	
	point3D<double> result((normalized.x - (1.0/(m_resolution.x*2.0)))*(double)m_resolution.x,
							  (normalized.y - (1.0/(m_resolution.y*2.0)))*(double)m_resolution.y,
							  (normalized.z - (1.0/(m_resolution.z*2.0)))*(double)m_resolution.z);
	return result;
}


template <class T>
T* rtsFunction3D<T>::GetBits()
{
	/*Returns bit data in lexocographical order (possibly for 3D texture mapping)*/
	return m_data;
}

template <class T>
rtsFunction3D<T>* rtsFunction3D<T>::Resample(indextype newres_x, indextype newres_y, indextype newres_z)
{
	/*This function resamples the current function at the specified resolution.
	No convolution is done for reducing he resolution.
	*/

	rtsFunction3D<T>* result = new rtsFunction3D<T>(vector3D<indextype>(newres_x, newres_y, newres_z),
						    m_boundary, m_domain_min, m_domain_max);

	//run through the entire resolution of the new function, sampling the current function
	int x, y, z;
	point3D<double> parametric;
	for(x = 0; x<newres_x; x++)
		for(y=0; y<newres_y; y++)
			for(z=0; z<newres_z; z++)
			{
				//compute the parametric point for the sample point
				parametric = result->getParameter(x, y, z);
				(*result)(x, y, z) = ijk(parametric.x, parametric.y, parametric.z);
			}

	return result;
}

template <class T>
void rtsFunction3D<T>::Scale(T new_min, T new_max)
{
	/*This function scales all values of the implicit function to within a specified range
	*/

	//find the minimum and maximum values in this function
	indextype data_size = m_resolution.x * m_resolution.y * m_resolution.z;
	T min = m_data[0];
	T max = m_data[0];
	for(indextype i=0; i<data_size; i++)
	{
		if(m_data[i] < min)
			min = m_data[i];
		if(m_data[i] > max)
			max = m_data[i];
	}

	//scale all values to the specified range
	T current_range = max - min;
	T new_range = new_max - new_min;
	for(indextype i=0; i<data_size; i++)
		m_data[i] = ((m_data[i] - min)/current_range)*(new_range) + new_min;
}

template <class T>
void rtsFunction3D<T>::Crop(indextype x, indextype y, indextype z, 
							indextype size_x, indextype size_y, indextype size_z)
{
	/*This function crops the implicit function at the specified nodes
	*/
	//create a pointer for the new data
	T* new_data = new T[size_x*size_y*size_z];

	//blit from the old data to the new data
	blit3D(m_data,
			x, y, z,
			m_resolution.x, m_resolution.y, m_resolution.z,
			new_data,
			0, 0, 0,
			size_x, size_y, size_z,
			size_x, size_y, size_z);

	//change the shallow variables
	vector3D<indextype> new_resolution = vector3D<indextype>(size_x, size_y, size_z);
	vector3D<double> voxel_size = getParameter(0,0,0) - getParameter(1,1,1);
	point3D<double> new_domain_min = getParameter(x, y, z) - 0.5*voxel_size;
	point3D<double> new_domain_max = getParameter(size_x-1, size_y - 1, size_z-1) + 0.5*voxel_size;
	//copy new shallow variables
	m_resolution = new_resolution;
	m_domain_min = new_domain_min;
	m_domain_max = new_domain_max;

	//copy data
	delete m_data;
	m_data = new_data;

}

template <class T>
void rtsFunction3D<T>::Threshold(T min, T value)
{
	/*This function sets all values between min and max to value.
	*/
	int x, y, z;
	T test_value;
	for(x=0; x<m_resolution.x; x++)
		for(y=0; y<m_resolution.y; y++)
			for(z=0; z<m_resolution.z; z++)
			{
				test_value = xyz(x, y, z);
				if(test_value >= min)
					xyz(x, y, z) = value;
			}
}

template <class T>
void rtsFunction3D<T>::Threshold(T min, T max, T value)
{
	/*This function sets all values between min and max to value.
	*/
	int x, y, z;
	T test_value;
	for(x=0; x<m_resolution.x; x++)
		for(y=0; y<m_resolution.y; y++)
			for(z=0; z<m_resolution.z; z++)
			{
				test_value = xyz(x, y, z);
				if(test_value >= min && test_value <= max)
					xyz(x, y, z) = value;
			}
}

template <class T>
void rtsFunction3D<T>::Threshold(T min, T max, T inside, T outside)
{
	/*This function sets all values between min and max to value.
	*/
	int x, y, z;
	T test_value;
	for(x=0; x<m_resolution.x; x++)
		for(y=0; y<m_resolution.y; y++)
			for(z=0; z<m_resolution.z; z++)
			{
				test_value = xyz(x, y, z);
				if(test_value >= min && test_value <= max)
					xyz(x, y, z) = inside;
				else
					xyz(x, y, z) = outside;
			}
}

template <class T>
void rtsFunction3D<T>::Insert(rtsFunction3D<T>* source, indextype x, indextype y, indextype z)
{
	blit3D(source->m_data, 0, 0, 0, source->m_resolution.x, source->m_resolution.y, source->m_resolution.z,
			m_data, x, y, z, m_resolution.x, m_resolution.y, m_resolution.z,
			source->m_resolution.x, source->m_resolution.y, source->m_resolution.z);
}




template <class T>
void rtsFunction3D<T>::Binary(T threshold, T true_value)
{
	/**
	This function converts an implicit function into a binary or characteristic function describing the solid represented by the level
	set at isovalue "threshold".  All values below threshold are set to zero while all values above threshold are set to the specified
	"true_value".  In order to use this function, the data type T must be able to be set to 0.
	**/
	int max_index = m_resolution.x * m_resolution.y * m_resolution.z;	//find the size of the data array
	int i;
	for(i=0; i<max_index; i++)
		if(m_data[i] >= threshold)
			m_data[i] = true_value;
		else
			m_data[i] = 0;
}



template <class T>
void rtsFunction3D<T>::ClampMax(T max)
{
	int i;
	int elements = m_resolution.x * m_resolution.y * m_resolution.z;
	for(i=0; i<elements; i++)
		if(m_data[i] > max)
			m_data[i] = max;
}

template <class T>
void rtsFunction3D<T>::ClampMin(T min)
{
	int i;
	int elements = m_resolution.x * m_resolution.y * m_resolution.z;
	for(i=0; i<elements; i++)
		if(m_data[i] < min)
			m_data[i] = min;
}

template <class T>
T rtsFunction3D<T>::getMin()
{
	int i;
	int elements = m_resolution.x * m_resolution.y * m_resolution.z;
	T current = m_data[0];
	for(i=1; i<elements; i++)
		if(m_data[i] < current)
			current = m_data[i];
	return current;
}

template <class T>
T rtsFunction3D<T>::getMax()
{
	int i;
	int elements = m_resolution.x * m_resolution.y * m_resolution.z;
	T current = m_data[0];
	for(i=1; i<elements; i++)
		if(m_data[i] > current)
			current = m_data[i];
	return current;
}


template <class T>
void rtsFunction3D<T>::toConsole()
{
	cout<<endl;
	int x, y, z;
	for(z=0; z<m_resolution.z; z++)
	{
		for(y=0; y<m_resolution.y; y++)
		{
			for(x=0; x<m_resolution.x; x++)
			{
				cout.width(7);
				cout.precision(3);
				cout<<(double)xyz(x, y, z);
			}
			cout<<endl;
		}
		cout<<"-----------------------------"<<endl;
	}

}

template <class T>
rtsFunction3D<T> rtsFunction3D<T>::Project2D()
{
	/**
	This function projects the entire 3D function onto a 2D function along the z-axis.
	**/
	rtsFunction3D<T> result(m_resolution.x, m_resolution.y, 1);
	result = 0;
	
	indextype x, y, z;
	for(x = 0; x<m_resolution.x; x++)
		for(y=0; y<m_resolution.y; y++)
			for(z=0; z<m_resolution.z; z++)
			{
				if(result(x, y, 0) < xyz(x, y, z))
					result(x, y, 0) = xyz(x, y, z);
			}
	return result;
}


#endif