#ifndef RTSIMPLICIT3D_H
#define RTSIMPLICIT3D_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 rtsImplicit3D
{
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 rtsImplicit3D<T> source, rtsImplicit3D<T> &dest);
	inline point3D<double> getParameter(indextype i);

	inline float isosurface_distance(point3D<double> p0, point3D<double> p1, T isovalue);
	inline float manhattan_distance(rtsImplicit3D<float>* function, point3D<indextype> p, bool sdf = false);
	void compute_distance_function_boundary(T isovalue, rtsImplicit3D<float>* &result, rtsImplicit3D<bool>* &mask, bool sdf = false);


public:
	//construct an implicit function with a size of 1
	rtsImplicit3D();			///<Create an empty implicit function
	//construct an implicit function of the specified resolution
	rtsImplicit3D(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
	rtsImplicit3D(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
	rtsImplicit3D(vector3D<int> resolution, T boundary, point3D<double> min_domain, point3D<double> max_domain);
	//full copy constructor, defines all variables
	rtsImplicit3D(T* data, vector3D<int> resolution, T boundary, point3D<double> min_domain, point3D<double> max_domain);
	rtsImplicit3D(const rtsImplicit3D<T> &original);	//copy constructor
	~rtsImplicit3D();		//destructor

	//overloaded operators
	rtsImplicit3D<T>& operator=(const rtsImplicit3D<T>& original);		///<Overloaded operator creates a copy of an implicit function
	rtsImplicit3D<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.
	rtsImplicit3D<T>& operator*=(const T constant);						///<Multiplies the values at all sample points by a constant.
	rtsImplicit3D<T>& operator+=(const T constant);						///<Adds a constant to the values at all sample points.
	rtsImplicit3D<T>& operator-=(const T constant);						///<Subtracts a constant from the values at all sample points.
	rtsImplicit3D<T>& operator/=(const T constant);						///<Divides all values by a constant.
	const rtsImplicit3D<T> operator+(const T constant);					///<Adds a constant to an implicit function and returns a new function.
	const rtsImplicit3D<T> operator-(const T constant);					///<Subtracts a constant from an implicit function and returns a new function.
	const rtsImplicit3D<T> operator*(const T constant);					///<Multiplies an implicit function by a constant and returns a new function.
	const rtsImplicit3D<T> operator/(const T constant);					///<Divides an implicit function by a constant and returns a new function.

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

	//friend classes for overloading "backwards" operations (like 3*function)
	friend rtsImplicit3D<T> operator*(const T lhs, rtsImplicit3D<T> rhs){return rhs*lhs;}	///<Allows associative multiplication.
	friend rtsImplicit3D<T> operator+(const T lhs, rtsImplicit3D<T> rhs){return rhs+lhs;}	///<Allows associative addition.
	friend rtsImplicit3D<T> operator-(const T lhs, rtsImplicit3D<T> rhs)					///<Allows associative subtraction.
	{
		rtsImplicit3D<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 rtsImplicit3D<T> operator/(const T lhs, rtsImplicit3D<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.

	//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;}
	vector<point3D<indextype>> getEdgeNodes(T isovalue, bool protrusions = true); ///<Returns a vector nodes lying on the edge of an implicit surface.
	rtsImplicit3D<T> Project2D();			//<Projects the data along the z-axis using a maximum-intensity projection.
	unsigned int BackgroundComponents6(indextype x, indextype y, indextype z, T threshold, int n=18);	///<Returns the number of background components 6-connected to the specified node.
	unsigned int Neighbors6(indextype x, indextype y, indextype z, T threshold); //<Returns the number of 6-connected neighbors above threshold.

	//data input methods
	void Insert(rtsDTGrid3D<T>* dt_grid, double factor, indextype x, indextype y, indextype z);
	void Insert(rtsImplicit3D<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 Erode(T isovalue, T fill_value);			///<Erodes the image around the edges defined by an isovalue.
	unsigned int Thin(T isovalue);		///<Finds the skeleton of the given isosurface using erosion.
	bool TestTopology(T isovalue, unsigned int x, unsigned int y, unsigned int z);	///<Tests if the specified voxel is necessary to the topology of the specified isosurface.
	int Neighbors26(indextype x, indextype y, indextype z, T isovalue);		///<Returns the number of neighbors above the given isovalue.
	void FloodFill26(indextype x, indextype y, indextype z, T new_value);	///<Fills the connected component at (x, y, z) with the new value.
	void FloodFill6(indextype x, indextype y, indextype z, T new_value);	///<Performs a 6-connected flood-fill operation starting at the specified node.
	void MedianFilter(int dist_x, int dist_y, int dist_z, double factor = 0.5);		///<Performs a median filter on the data set.
	void ClampMax(T max);		///<Clamps the function to the given maximum value.
	void ClampMin(T min);
	void ClampMin(T min, T value);	///<Sets all function values below min to value.

	//create new data
	rtsImplicit3D<T>* Resample(indextype newres_x, indextype newres_y, indextype newres_z);
	rtsImplicit3D<float>* Isodistance_Manhattan(T isovalue, bool sdf = false);
	rtsImplicit3D<vector3D<T>>* Gradient();
	rtsImplicit3D<float>* EstimateAmbient(T threshold);
	rtsImplicit3D<float>* EstimateAttenuatedAmbient(T surface, T transparent, float attenuation);

	//implicit shapes
	//(these functions create some basic implicit shapes just for fun)
	void Sphere(double center_i, double center_j, double center_k, double radius, T in_value);

	//output functions
	void toConsole();

};

template <class T>
void rtsImplicit3D<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 rtsImplicit3D<T>::shallow_copy(const rtsImplicit3D<T> source, rtsImplicit3D<T> &dest)
{
	dest = rtsImplicit3D<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>
rtsImplicit3D<T>::rtsImplicit3D(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>
rtsImplicit3D<T>::rtsImplicit3D(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>
rtsImplicit3D<T>::rtsImplicit3D()
{
	m_resolution.x = 1;
	m_resolution.y = 1;
	m_resolution.z = 1;
	m_data = new T[1];
	//m_boundary = 0;				//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>
rtsImplicit3D<T>::rtsImplicit3D(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
	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>
rtsImplicit3D<T>::rtsImplicit3D(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];
	m_boundary = 0;							//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>
rtsImplicit3D<T>::rtsImplicit3D(const rtsImplicit3D<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>
rtsImplicit3D<T>::~rtsImplicit3D()
{
	delete m_data;
}

template <class T>
typename rtsImplicit3D<T>& rtsImplicit3D<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 rtsImplicit3D<T>& rtsImplicit3D<T>::operator=(const rtsImplicit3D<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& rtsImplicit3D<T>::operator ()(indextype x, indextype y, indextype z)
{
	return xyz(x, y, z);
}

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

template <class T>
rtsImplicit3D<T>& rtsImplicit3D<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>
rtsImplicit3D<T>& rtsImplicit3D<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>
rtsImplicit3D<T>& rtsImplicit3D<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>
rtsImplicit3D<T>& rtsImplicit3D<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 rtsImplicit3D<T> rtsImplicit3D<T>::operator *(const T constant)
{
	rtsImplicit3D<T> result = (*this);
	result *= constant;

	return result;
}

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

	return result;
}

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

	return result;
}

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

	return result;
}

template <class T>
template <class U>
rtsImplicit3D<T>::operator rtsImplicit3D<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];
	rtsImplicit3D<U> cast_result(new_data, m_resolution, m_boundary, m_domain_min, m_domain_max);

	return cast_result;
}

template <class T>
inline T& rtsImplicit3D<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> rtsImplicit3D<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 rtsImplicit3D<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;

	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 = new T[size];
	infile.read((char*)m_data, size*sizeof(T));

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

template <class T>
void rtsImplicit3D<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 rtsImplicit3D<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(char)*vol_size.x*vol_size.y*vol_size.z);
}

template <class T>
void rtsImplicit3D<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 rtsImplicit3D<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 rtsImplicit3D<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> rtsImplicit3D<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> rtsImplicit3D<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> rtsImplicit3D<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* rtsImplicit3D<T>::GetBits()
{
	/*Returns bit data in lexocographical order (possibly for 3D texture mapping)*/
	return m_data;
}

template <class T>
rtsImplicit3D<T>* rtsImplicit3D<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.
	*/

	rtsImplicit3D<T>* result = new rtsImplicit3D<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 rtsImplicit3D<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 rtsImplicit3D<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 rtsImplicit3D<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 rtsImplicit3D<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 rtsImplicit3D<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 rtsImplicit3D<T>::Insert(rtsDTGrid3D<T>* dt_grid, double factor, indextype pos_x, indextype pos_y, indextype pos_z)
{
/*	This function copies a 3D DT-Grid into the implicit function at the specified position.
*/
	rtsDTGrid3D<T>::iterator i;
	indextype x, y, z;
	for(i=dt_grid->begin(); i != dt_grid->end(); i.increment())	//for each node in the grid
	{
		x = pos_x + i.getX();
		y = pos_y + i.getY();
		z = pos_z + i.getZ();
		if(x >= 0 || x < m_resolution.x ||
			y >= 0 || y < m_resolution.y ||
			y >= 0 || y < m_resolution.z)
			xyz(x, y, z) = max(i.value* factor, (double)xyz(x, y, z));
	}
}

template <class T>
void rtsImplicit3D<T>::Insert(rtsImplicit3D<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>
inline float rtsImplicit3D<T>::manhattan_distance(rtsImplicit3D<float>* function, point3D<indextype> point, bool sdf)
{
	/*This function updates the manhattan distance from a surface using the manhattan
	distance of its neighboring points.
	*/
	indextype x, y, z;
	x=point.x; y=point.y, z=point.z;
	int sign = 1;
	float result = DIST_MAX;
	float near_value;				//the value of the neighbor being considered
	float possible_value;			
	if(x!=0)
	{ 
		near_value = (*function)(x-1, y, z);
		if(!sdf)
			result = min(result, near_value + (float)m_voxel_size.x);
		else
		{
			if(near_value<0) sign = -1; else sign = 1;	//determine if the value is inside or outside
			possible_value = sign*(fabs(near_value) + m_voxel_size.x);
			if(fabs(possible_value) < fabs(result))
				result = possible_value;
		}

	}
	if(x!=function->DimX()-1)
	{
		near_value = (*function)(x+1, y, z);
		if(!sdf)
			result = min(result, near_value + (float)m_voxel_size.x);
		else
		{
			if(near_value<0) sign = -1; else sign = 1;	//determine if the value is inside or outside
			possible_value = sign*(fabs(near_value) + m_voxel_size.x);
			if(fabs(possible_value) < fabs(result))
				result = possible_value;
		}
	}
	if(y!=0)
	{
		near_value = (*function)(x, y-1, z);
		if(!sdf)
			result = min(result, near_value + (float)m_voxel_size.y);
		else
		{
			if(near_value<0) sign = -1; else sign = 1;	//determine if the value is inside or outside
			possible_value = sign*(fabs(near_value) + m_voxel_size.y);
			if(fabs(possible_value) < fabs(result))
				result = possible_value;
		}
	}
	if(y!=function->DimY()-1)
	{
		near_value = (*function)(x, y+1, z);
		if(!sdf)
			result = min(result, near_value + (float)m_voxel_size.y);
		else
		{
			if(near_value<0) sign = -1; else sign = 1;	//determine if the value is inside or outside
			possible_value = sign*(fabs(near_value) + m_voxel_size.y);
			if(fabs(possible_value) < fabs(result))
				result = possible_value;
		}
	}
	if(z!=0)
	{
		near_value = (*function)(x, y, z-1);
		if(!sdf)
			result = min(result, near_value + (float)m_voxel_size.z);
		else
		{
			if(near_value<0) sign = -1; else sign = 1;	//determine if the value is inside or outside
			possible_value = sign*(fabs(near_value) + m_voxel_size.z);
			if(fabs(possible_value) < fabs(result))
				result = possible_value;
		}
	}
	if(z!=function->DimZ()-1)
	{
		near_value = (*function)(x, y, z+1);
		if(!sdf)
			result = min(result, near_value + (float)m_voxel_size.z);
		else
		{
			if(near_value<0) sign = -1; else sign = 1;	//determine if the value is inside or outside
			possible_value = sign*(fabs(near_value) + m_voxel_size.z);
			if(fabs(possible_value) < fabs(result))
				result = possible_value;
		}
	}
	return result;
}

template <class T>
inline float rtsImplicit3D<T>::isosurface_distance(point3D<double> p0, point3D<double> p1, T isovalue)
{
	/*This function computes the distance from p0 to the surface, given two points p0 and p1
	on either side of the surface.  isovalue specifies
	the value at the surface.  Right now, this function returns a float.  I'll have to think
	of something better to do in the future.
	*/

	//compute the normalized position of the surface between p0 and p1
	float val0 = ijk(p0.x, p0.y, p0.z);
	float val1 = ijk(p1.x, p1.y, p1.z);
	float isovalue_norm_pos = (isovalue - val0) / (val1 - val0);
	//compute the actual position of the surface
	point3D<double> s_pos = p0 + isovalue_norm_pos * (p1 - p0);
	//compute the distance from p0 to the surface
	float result = (s_pos - p0).Length();
	//cout<<"distance: "<<result<<endl;
	return result;
}

template <class T>
void rtsImplicit3D<T>::compute_distance_function_boundary(T isovalue, 
														  rtsImplicit3D<float>* &result, 
														  rtsImplicit3D<bool>* &mask, bool sdf)
{
	/*This function creates an initial signed distance function from a threshold image.
	All voxels adjacent to the surface specified by the threshold are initialized with a
	distance value.  Low values are inside, high values are outside.
	*/
	//current and neighboring voxel flags (false = inside, true = outside)
	bool c, x_p, x_n, y_p, y_n, z_p, z_n;
	float d_xp, d_xn, d_yp, d_yn, d_zp, d_zn;
	float in_out = 1;

	//boundary condition function and the mask
	result = new rtsImplicit3D<float>(m_resolution.x, m_resolution.y, m_resolution.z);
	//get the parameterization
	result->Parameterize(m_domain_min.x, m_domain_max.x, m_domain_min.y, m_domain_max.y, m_domain_min.z, m_domain_max.z);
	(*result) = DIST_MAX;
	result->setBoundary(DIST_MAX);	
	//create a mask
	mask = new rtsImplicit3D<bool>(m_resolution.x, m_resolution.y, m_resolution.z);
	(*mask) = false;

	cout<<"done making boundary condition function"<<endl;
	//for each voxel
	int 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++)
			{
				//reset flags
				c=x_p=x_n=y_p=y_n=z_p=z_n=true;
				in_out = 1.0;
				//look at the current voxel
				if(xyz(x, y, z) < isovalue)
					c=false;
				else c=true;
				//if the voxel is outside the domain, assume that it is equal to the current voxel
				if(x-1 < 0) x_n = c;	//X
				else if(xyz(x-1, y, z) < isovalue) x_n = false;
				if(x+1 >= m_resolution.x) x_p = c;
				else if(xyz(x+1, y, z) < isovalue) x_p = false;
				if(y-1 < 0) y_n = c;	//Y
				else if(xyz(x, y-1, z) < isovalue) y_n = false;
				if(y+1 >= m_resolution.y) y_p = c;
				else if(xyz(x, y+1, z) < isovalue) y_p = false;
				if(z-1 < 0) z_n = c;	//Z
				else if(xyz(x, y, z-1) < isovalue) z_n = false;
				if(z+1 >= m_resolution.z) z_p = c;
				else if(xyz(x, y, z+1) < isovalue) z_p = false;

				//set the distance from the isosurface
				if(c == false && sdf)
					in_out = -1.0;
				if(x_n != c)
					(*result)(x, y, z) = min((*result)(x,y,z),
										 isosurface_distance(getParameter(x, y, z), 
													  getParameter(x-1, y, z),
													  isovalue) * in_out);
				if(x_p != c)
					(*result)(x, y, z) = min((*result)(x,y,z),
											isosurface_distance(getParameter(x, y, z), 
													  getParameter(x+1, y, z),
													  isovalue) * in_out);
				if(y_n != c)
					(*result)(x, y, z) = min((*result)(x,y,z),
											isosurface_distance(getParameter(x, y, z), 
													  getParameter(x, y-1, z),
													  isovalue) * in_out);
				if(y_p != c)
					(*result)(x, y, z) = min((*result)(x,y,z),
											isosurface_distance(getParameter(x, y, z), 
													  getParameter(x, y+1, z),
													  isovalue) * in_out);
				if(z_n != c)
					(*result)(x, y, z) = min((*result)(x,y,z),
											isosurface_distance(getParameter(x, y, z-1), 
													  getParameter(x, y, z),
													  isovalue) * in_out);
				if(z_p != c)
					(*result)(x, y, z) = min((*result)(x,y,z),
											isosurface_distance(getParameter(x, y, z), 
													  getParameter(x, y, z+1),
													  isovalue) * in_out);

				//set the mask to 1 if the voxel is on an edge node
				if(x_n != c || x_p != c || y_n != c || y_p != c || z_n != c || z_p != c)
					(*mask)(x, y, z) = true;
			}
				

	//if a line between the two voxels crosses the surface
		//find the distance between the voxel center and the surface


			cout<<"done computing boundary conditions"<<endl;
}

template <class T>
rtsImplicit3D<float>* rtsImplicit3D<T>::EstimateAmbient(T threshold)
{
	rtsImplicit3D<float>* result = new rtsImplicit3D<float>(m_resolution.x, m_resolution.y, m_resolution.z);		//create a new implicit function
	(*result) = 0.0f;
	rtsImplicit3D<float>* temp = new rtsImplicit3D<float>(m_resolution.x, m_resolution.y, m_resolution.z);	//temp buffer for current lighting iteration
	(*temp) = 0.0f;
	temp->setBoundary(1.0);
	
	cout<<"first iteration..."<<endl;
	int x,y,z;
	float ambient;
	
	for(x=0; x<m_resolution.x; x++)
		for(y=0; y<m_resolution.y; y++)
			for(z=0; z<m_resolution.z; z++)
			{
				ambient = 0.0;
				if(xyz(x-1, y, z) < threshold)
					ambient += (*temp)(x-1, y, z);
				if(xyz(x, y-1, z) < threshold)
					ambient += (*temp)(x, y-1, z);
				if(xyz(x, y-1, z) < threshold)
					ambient += (*temp)(x, y, z-1);

				(*temp)(x, y, z) += ambient/3.0;
				(*result)(x, y, z) += ambient/3.0;
				if(ambient > 3.0)
					cout<<"error"<<endl;
			}
	(*temp) = 0.0f;
	temp->setBoundary(1.0);
	cout<<"done."<<endl;

	cout<<"second iteration..."<<endl;
	for(x=0; x<m_resolution.x; x++)
		for(y=0; y<m_resolution.y; y++)
			for(z=m_resolution.z-1; z>=0; z--)
			{
				ambient = 0.0;
				if(xyz(x-1, y, z) < threshold)
					ambient += (*temp)(x-1, y, z);
				if(xyz(x, y-1, z) < threshold)
					ambient += (*temp)(x, y-1, z);
				if(xyz(x, y, z+1) < threshold)
					ambient += (*temp)(x, y, z+1);

				(*temp)(x, y, z) += ambient/3.0;
				(*result)(x, y, z) += ambient/3.0;
				if(ambient > 3.0)
					cout<<"error"<<endl;
			}
	(*temp) = 0.0f;
	temp->setBoundary(1.0);
	cout<<"done."<<endl;

	cout<<"third iteration..."<<endl;
	for(x=0; x<m_resolution.x; x++)
		for(y=m_resolution.y-1; y>=0; y--)
			for(z=0; z<m_resolution.z; z++)
			{
				ambient = 0.0;
				if(xyz(x-1, y, z) < threshold)
					ambient += (*temp)(x-1, y, z);
				if(xyz(x, y+1, z) < threshold)
					ambient += (*temp)(x, y+1, z);
				if(xyz(x, y, z-1) < threshold)
					ambient += (*temp)(x, y, z-1);

				(*temp)(x, y, z) += ambient/3.0;
				(*result)(x, y, z) += ambient/3.0;
				if(ambient > 3.0)
					cout<<"error"<<endl;
			}
	(*temp) = 0.0f;
	temp->setBoundary(1.0);
	cout<<"done."<<endl;

	cout<<"fourth iteration..."<<endl;
	for(x=0; x<m_resolution.x; x++)
		for(y=m_resolution.y-1; y>=0; y--)
			for(z=m_resolution.z-1; z>=0; z--)
			{
				ambient = 0.0;
				if(xyz(x-1, y, z) < threshold)
					ambient += (*temp)(x-1, y, z);
				if(xyz(x, y+1, z) < threshold)
					ambient += (*temp)(x, y+1, z);
				if(xyz(x, y, z+1) < threshold)
					ambient += (*temp)(x, y, z+1);

				(*temp)(x, y, z) += ambient/3.0;
				(*result)(x, y, z) += ambient/3.0;
				if(ambient > 3.0)
					cout<<"error"<<endl;
			}
	(*temp) = 0.0f;
	temp->setBoundary(1.0);
	cout<<"done."<<endl;

	cout<<"fifth iteration..."<<endl;
	for(x=m_resolution.x-1; x>=0; x--)
		for(y=0; y<m_resolution.y; y++)
			for(z=0; z<m_resolution.z; z++)
			{
				ambient = 0.0;
				if(xyz(x+1, y, z) < threshold)
					ambient += (*temp)(x+1, y, z);
				if(xyz(x, y-1, z) < threshold)
					ambient += (*temp)(x, y-1, z);
				if(xyz(x, y, z-1) < threshold)
					ambient += (*temp)(x, y, z-1);

				(*temp)(x, y, z) += ambient/3.0;
				(*result)(x, y, z) += ambient/3.0;
				if(ambient > 3.0)
					cout<<"error"<<endl;
			}
	(*temp) = 0.0f;
	temp->setBoundary(1.0);
	cout<<"done."<<endl;

	cout<<"sixth iteration..."<<endl;
	for(x=m_resolution.x-1; x>=0; x--)
		for(y=0; y<m_resolution.y; y++)
			for(z=m_resolution.z-1; z>=0; z--)
			{
				ambient = 0.0;
				if(xyz(x+1, y, z) < threshold)
					ambient += (*temp)(x+1, y, z);
				if(xyz(x, y-1, z) < threshold)
					ambient += (*temp)(x, y-1, z);
				if(xyz(x, y, z+1) < threshold)
					ambient += (*temp)(x, y, z+1);

				(*temp)(x, y, z) += ambient/3.0;
				(*result)(x, y, z) += ambient/3.0;
				if(ambient > 3.0)
					cout<<"error"<<endl;
			}
	(*temp) = 0.0f;
	temp->setBoundary(1.0);
	cout<<"done."<<endl;

	cout<<"seventh iteration..."<<endl;
	for(x=m_resolution.x-1; x>=0; x--)
		for(y=m_resolution.y-1; y>=0; y--)
			for(z=0; z<m_resolution.z; z++)
			{
				ambient = 0.0;
				if(xyz(x+1, y, z) < threshold)
					ambient += (*temp)(x+1, y, z);
				if(xyz(x, y+1, z) < threshold)
					ambient += (*temp)(x, y+1, z);
				if(xyz(x, y, z-1) < threshold)
					ambient += (*temp)(x, y, z-1);

				(*temp)(x, y, z) += ambient/3.0;
				(*result)(x, y, z) += ambient/3.0;
				if(ambient > 3.0)
					cout<<"error"<<endl;
			}
	(*temp) = 0.0f;
	temp->setBoundary(1.0);
	cout<<"done."<<endl;

	cout<<"eighth iteration..."<<endl;
	for(x=m_resolution.x-1; x>=0; x--)
		for(y=m_resolution.y-1; y>=0; y--)
			for(z=m_resolution.z-1; z>=0; z--)
			{
				ambient = 0.0;
				if(xyz(x+1, y, z) < threshold)
					ambient += (*temp)(x+1, y, z);
				if(xyz(x, y+1, z) < threshold)
					ambient += (*temp)(x, y+1, z);
				if(xyz(x, y, z+1) < threshold)
					ambient += (*temp)(x, y, z+1);

				(*temp)(x, y, z) += ambient/3.0;
				(*result)(x, y, z) += ambient/3.0;
				if(ambient > 3.0)
					cout<<"error"<<endl;
			}
	(*temp) = 0.0f;
	temp->setBoundary(1.0);
	cout<<"done."<<endl;

	(*result)/=8.0;
	
	return result;
}


template <class T>
rtsImplicit3D<float>* rtsImplicit3D<T>::EstimateAttenuatedAmbient(T threshold, T transparent, float attenuation)
{
	rtsImplicit3D<float>* result = new rtsImplicit3D<float>(m_resolution.x, m_resolution.y, m_resolution.z);		//create a new implicit function
	(*result) = 0.0f;
	rtsImplicit3D<float>* temp = new rtsImplicit3D<float>(m_resolution.x, m_resolution.y, m_resolution.z);	//temp buffer for current lighting iteration
	(*temp) = 0.0f;
	temp->setBoundary(1.0);
	
	cout<<"first iteration..."<<endl;
	int x,y,z;
	float ambient;
	
	for(x=0; x<m_resolution.x; x++)
		for(y=0; y<m_resolution.y; y++)
			for(z=0; z<m_resolution.z; z++)
			{
				ambient = 0.0;
				if(xyz(x-1, y, z) < threshold)
					ambient += (*temp)(x-1, y, z)*(1.0 - (xyz(x-1, y, z)/255.0)*attenuation);
				if(xyz(x, y-1, z) < threshold)
					ambient += (*temp)(x, y-1, z)*(1.0 - (xyz(x, y-1, z)/255.0)*attenuation);
				if(xyz(x, y-1, z) < threshold)
					ambient += (*temp)(x, y, z-1)*(1.0 - (xyz(x, y, z-1)/255.0)*attenuation);

				(*temp)(x, y, z) += ambient/3.0;
				(*result)(x, y, z) += min(1.0, ambient/3.0);
				//if(ambient > 3.0)
				//	cout<<"error"<<endl;
			}
	(*temp) = 0.0f;
	temp->setBoundary(1.0);
	cout<<"done."<<endl;

	cout<<"second iteration..."<<endl;
	for(x=0; x<m_resolution.x; x++)
		for(y=0; y<m_resolution.y; y++)
			for(z=m_resolution.z-1; z>=0; z--)
			{
				ambient = 0.0;
				if(xyz(x-1, y, z) < threshold)
					ambient += (*temp)(x-1, y, z)*(1.0 - (xyz(x-1, y, z)/255.0)*attenuation);
				if(xyz(x, y-1, z) < threshold)
					ambient += (*temp)(x, y-1, z)*(1.0 - (xyz(x, y-1, z)/255.0)*attenuation);
				if(xyz(x, y, z+1) < threshold)
					ambient += (*temp)(x, y, z+1)*(1.0 - (xyz(x, y, z+1)/255.0)*attenuation);

				(*temp)(x, y, z) += ambient/3.0;
				(*result)(x, y, z) += ambient/3.0;
				if(ambient > 3.0)
					cout<<"error"<<endl;
			}
	(*temp) = 0.0f;
	temp->setBoundary(1.0);
	cout<<"done."<<endl;

	cout<<"third iteration..."<<endl;
	for(x=0; x<m_resolution.x; x++)
		for(y=m_resolution.y-1; y>=0; y--)
			for(z=0; z<m_resolution.z; z++)
			{
				ambient = 0.0;
				if(xyz(x-1, y, z) < threshold)
					ambient += (*temp)(x-1, y, z)*(1.0 - (xyz(x-1, y, z)/255.0)*attenuation);
				if(xyz(x, y+1, z) < threshold)
					ambient += (*temp)(x, y+1, z)*(1.0 - (xyz(x, y+1, z)/255.0)*attenuation);
				if(xyz(x, y, z-1) < threshold)
					ambient += (*temp)(x, y, z-1)*(1.0 - (xyz(x, y, z-1)/255.0)*attenuation);

				(*temp)(x, y, z) += ambient/3.0;
				(*result)(x, y, z) += ambient/3.0;
				if(ambient > 3.0)
					cout<<"error"<<endl;
			}
	(*temp) = 0.0f;
	temp->setBoundary(1.0);
	cout<<"done."<<endl;

	cout<<"fourth iteration..."<<endl;
	for(x=0; x<m_resolution.x; x++)
		for(y=m_resolution.y-1; y>=0; y--)
			for(z=m_resolution.z-1; z>=0; z--)
			{
				ambient = 0.0;
				if(xyz(x-1, y, z) < threshold)
					ambient += (*temp)(x-1, y, z)*(1.0 - (xyz(x-1, y, z)/255.0)*attenuation);
				if(xyz(x, y+1, z) < threshold)
					ambient += (*temp)(x, y+1, z)*(1.0 - (xyz(x, y+1, z)/255.0)*attenuation);
				if(xyz(x, y, z+1) < threshold)
					ambient += (*temp)(x, y, z+1)*(1.0 - (xyz(x, y, z+1)/255.0)*attenuation);

				(*temp)(x, y, z) += ambient/3.0;
				(*result)(x, y, z) += ambient/3.0;
				if(ambient > 3.0)
					cout<<"error"<<endl;
			}
	(*temp) = 0.0f;
	temp->setBoundary(1.0);
	cout<<"done."<<endl;

	cout<<"fifth iteration..."<<endl;
	for(x=m_resolution.x-1; x>=0; x--)
		for(y=0; y<m_resolution.y; y++)
			for(z=0; z<m_resolution.z; z++)
			{
				ambient = 0.0;
				if(xyz(x+1, y, z) < threshold)
					ambient += (*temp)(x+1, y, z)*(1.0 - (xyz(x+1, y, z)/255.0)*attenuation);
				if(xyz(x, y-1, z) < threshold)
					ambient += (*temp)(x, y-1, z)*(1.0 - (xyz(x, y-1, z)/255.0)*attenuation);
				if(xyz(x, y, z-1) < threshold)
					ambient += (*temp)(x, y, z-1)*(1.0 - (xyz(x, y, z-1)/255.0)*attenuation);

				(*temp)(x, y, z) += ambient/3.0;
				(*result)(x, y, z) += ambient/3.0;
				if(ambient > 3.0)
					cout<<"error"<<endl;
			}
	(*temp) = 0.0f;
	temp->setBoundary(1.0);
	cout<<"done."<<endl;

	cout<<"sixth iteration..."<<endl;
	for(x=m_resolution.x-1; x>=0; x--)
		for(y=0; y<m_resolution.y; y++)
			for(z=m_resolution.z-1; z>=0; z--)
			{
				ambient = 0.0;
				if(xyz(x+1, y, z) < threshold)
					ambient += (*temp)(x+1, y, z)*(1.0 - (xyz(x+1, y, z)/255.0)*attenuation);
				if(xyz(x, y-1, z) < threshold)
					ambient += (*temp)(x, y-1, z)*(1.0 - (xyz(x, y-1, z)/255.0)*attenuation);
				if(xyz(x, y, z+1) < threshold)
					ambient += (*temp)(x, y, z+1)*(1.0 - (xyz(x, y, z+1)/255.0)*attenuation);

				(*temp)(x, y, z) += ambient/3.0;
				(*result)(x, y, z) += ambient/3.0;
				if(ambient > 3.0)
					cout<<"error"<<endl;
			}
	(*temp) = 0.0f;
	temp->setBoundary(1.0);
	cout<<"done."<<endl;

	cout<<"seventh iteration..."<<endl;
	for(x=m_resolution.x-1; x>=0; x--)
		for(y=m_resolution.y-1; y>=0; y--)
			for(z=0; z<m_resolution.z; z++)
			{
				ambient = 0.0;
				if(xyz(x+1, y, z) < threshold)
					ambient += (*temp)(x+1, y, z)*(1.0 - (xyz(x+1, y, z)/255.0)*attenuation);
				if(xyz(x, y+1, z) < threshold)
					ambient += (*temp)(x, y+1, z)*(1.0 - (xyz(x, y+1, z)/255.0)*attenuation);
				if(xyz(x, y, z-1) < threshold)
					ambient += (*temp)(x, y, z-1)*(1.0 - (xyz(x, y, z-1)/255.0)*attenuation);

				(*temp)(x, y, z) += ambient/3.0;
				(*result)(x, y, z) += ambient/3.0;
				if(ambient > 3.0)
					cout<<"error"<<endl;
			}
	(*temp) = 0.0f;
	temp->setBoundary(1.0);
	cout<<"done."<<endl;

	cout<<"eighth iteration..."<<endl;
	for(x=m_resolution.x-1; x>=0; x--)
		for(y=m_resolution.y-1; y>=0; y--)
			for(z=m_resolution.z-1; z>=0; z--)
			{
				ambient = 0.0;
				if(xyz(x+1, y, z) < threshold)
					ambient += (*temp)(x+1, y, z)*(1.0 - (xyz(x+1, y, z)/255.0)*attenuation);
				if(xyz(x, y+1, z) < threshold)
					ambient += (*temp)(x, y+1, z)*(1.0 - (xyz(x, y+1, z)/255.0)*attenuation);
				if(xyz(x, y, z+1) < threshold)
					ambient += (*temp)(x, y, z+1)*(1.0 - (xyz(x, y, z+1)/255.0)*attenuation);

				(*temp)(x, y, z) += ambient/3.0;
				(*result)(x, y, z) += ambient/3.0;
				if(ambient > 3.0)
					cout<<"error"<<endl;
			}
	(*temp) = 0.0f;
	temp->setBoundary(1.0);
	cout<<"done."<<endl;

	(*result)/=8.0;
	
	return result;
}

template <class T>
rtsImplicit3D<float>* rtsImplicit3D<T>::Isodistance_Manhattan(T isovalue, bool sdf)
{
	rtsImplicit3D<float>* function;
	rtsImplicit3D<bool>* mask;
	compute_distance_function_boundary(isovalue, function, mask, sdf);

	//compute the manhattan distance for the entire function
	//use fast sweeping to compute the manhattan distance
	//0:X  0:Y  0:Z
	cout<<"first iteration..."<<endl;
	int 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 the current point is not a boundary value
				if(!(*mask)(x, y, z))
					(*function)(x,y,z) = manhattan_distance(function, point3D<indextype>(x, y, z), sdf);
	cout<<"done."<<endl;
	cout<<"second iteration..."<<endl;
	//0:X 0:Y Z:0
	for(x=0; x<m_resolution.x; x++)
		for(y=0; y<m_resolution.y; y++)
			for(z=m_resolution.z-1; z>=0; z--)
				//if the current point is not a boundary value
				if(!(*mask)(x, y, z))
					(*function)(x,y,z) = manhattan_distance(function, point3D<indextype>(x, y, z), sdf);
	cout<<"done."<<endl;
	cout<<"third iteration..."<<endl;
	//0:X Y:0 0:Z
	for(x=0; x<m_resolution.x; x++)
		for(y=m_resolution.y-1; y>=0; y--)
			for(z=0; z<m_resolution.z; z++)
				//if the current point is not a boundary value
				if(!(*mask)(x, y, z))
					(*function)(x,y,z) = manhattan_distance(function, point3D<indextype>(x, y, z), sdf);
	cout<<"done."<<endl;
	cout<<"fourth iteration..."<<endl;
	//0:X Y:0 Z:0
	for(x=0; x<m_resolution.x; x++)
		for(y=m_resolution.y-1; y>=0; y--)
			for(z=m_resolution.z-1; z>=0; z--)
				//if the current point is not a boundary value
				if(!(*mask)(x, y, z))
					(*function)(x,y,z) = manhattan_distance(function, point3D<indextype>(x, y, z), sdf);
	cout<<"done."<<endl;
	cout<<"fifth iteration..."<<endl;
	//X:0 0:Y 0:Z
	for(x=m_resolution.x-1; x>=0; x--)
		for(y=0; y<m_resolution.y; y++)
			for(z=0; z<m_resolution.z; z++)
				//if the current point is not a boundary value
				if(!(*mask)(x, y, z))
					(*function)(x,y,z) = manhattan_distance(function, point3D<indextype>(x, y, z), sdf);
	cout<<"done."<<endl;
	cout<<"sixth iteration..."<<endl;
	//X:0 0:Y Z:0
	for(x=m_resolution.x-1; x>=0; x--)
		for(y=0; y<m_resolution.y; y++)
			for(z=m_resolution.z-1; z>=0; z--)
				//if the current point is not a boundary value
				if(!(*mask)(x, y, z))
					(*function)(x,y,z) = manhattan_distance(function, point3D<indextype>(x, y, z), sdf);
	cout<<"done."<<endl;
	cout<<"seventh iteration..."<<endl;
	//X:0 Y:0 0:Z
	for(x=m_resolution.x-1; x>=0; x--)
		for(y=m_resolution.y-1; y>=0; y--)
			for(z=0; z<m_resolution.z; z++)
				//if the current point is not a boundary value
				if(!(*mask)(x, y, z))
					(*function)(x,y,z) = manhattan_distance(function, point3D<indextype>(x, y, z), sdf);
	cout<<"done."<<endl;
	cout<<"eighth iteration..."<<endl;
	//X:0 Y:0 Z:0
	for(x=m_resolution.x-1; x>=0; x--)
		for(y=m_resolution.y-1; y>=0; y--)
			for(z=m_resolution.z-1; z>=0; z--)
				//if the current point is not a boundary value
				if(!(*mask)(x, y, z))
					(*function)(x,y,z) = manhattan_distance(function, point3D<indextype>(x, y, z), sdf);
	cout<<"done."<<endl;
	

	return function;
}

//computes the gradient along all three dimensions and returns a vector field
template <class T>
rtsImplicit3D<vector3D<T>>* rtsImplicit3D<T>::Gradient()
{
	int x, y, z;
	rtsImplicit3D<vector3D<T>>* result = new rtsImplicit3D<vector3D<T>>(m_resolution.x, m_resolution.y, m_resolution.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++)
			{
				result->xyz(x, y, z).x = xyz(x-1, y, z) - xyz(x, y, z);
				result->xyz(x, y, z).y = xyz(x, y-1, z) - xyz(x, y, z);
				result->xyz(x, y, z).z = xyz(x, y, z-1) - xyz(x, y, z);
			}
	return result;
}


template <class T>
int rtsImplicit3D<T>::Neighbors26(indextype x, indextype y, indextype z, T isovalue)
{
	int neighbors = 0;
	int u,v,w;

	for(u=-1; u<=1; u++)
		for(v=-1; v<=1; v++)
			for(w=-1; w<=1; w++)
				if(xyz(x+u, y+v, z+w) >= isovalue)
					neighbors++;
	if(xyz(x, y, z) > isovalue)
		neighbors--;

	return neighbors;
}

template <class T>
unsigned int rtsImplicit3D<T>::Neighbors6(indextype x, indextype y, indextype z, T threshold)
{

	unsigned int neighbors = 0;
	if(xyz(x+1, y, z) >= threshold)
		neighbors++;
	if(xyz(x-1, y, z) >= threshold)
		neighbors++;
	if(xyz(x, y+1, z) >= threshold)
		neighbors++;
	if(xyz(x, y-1, z) >= threshold)
		neighbors++;
	if(xyz(x, y, z+1) >= threshold)
		neighbors++;
	if(xyz(x, y, z-1) >= threshold)
		neighbors++;

	return neighbors;
}

template <class T>
bool rtsImplicit3D<T>::TestTopology(T isovalue, unsigned int x, unsigned int y, unsigned int z)
{
	if(xyz(x,y,z) < isovalue)
		return false;
	//This function returns true if a voxel is necessary, otherwise it returns false
	unsigned int neighbors = Neighbors(x, y, z, isovalue);

	if(neighbors == 3)
		return false;
	if(neighbors == 0 || neighbors == 1 || neighbors == 4)
		return true;
	if(neighbors == 2)
	{
		if(xyz(x-1, y, z) >= isovalue && xyz(x+1, y, z) >= isovalue)
			return true;
		if(xyz(x, y-1, z) >= isovalue && xyz(x, y+1, z) >= isovalue)
			return true;
		return false;
	}
	
}

template <class T>
void rtsImplicit3D<T>::FloodFill6(indextype x, indextype y, indextype z, T target_value)
{
	T old_value = xyz(x, y, z);					//find the old value (the value being flood-filled)
	if(target_value == old_value)				//if the target value is the same as the old value, nothing to do
		return;

	queue<point3D<indextype>> Q;				//create a queue for neighboring points
	point3D<indextype> current(x, y, z);		//start with the current point
	xyz(current.x, current.y, current.z) = target_value;
	point3D<indextype> next;
	Q.push(current);
	indextype u, v, w;

	while(!Q.empty())							//continue until the queue is empty
	{
		current = Q.front();					//get the first element from the queue
		Q.pop();							
		
		if(current.x != m_resolution.x - 1)
			if(xyz(current.x + 1, current.y, current.z) == old_value)
			{
				xyz(current.x + 1, current.y, current.z) = target_value;
				Q.push(point3D<indextype>(current.x + 1, current.y, current.z));
			}
		if(current.x != 0)
			if(xyz(current.x - 1, current.y, current.z) == old_value)
			{
				xyz(current.x - 1, current.y, current.z) = target_value;
				Q.push(point3D<indextype>(current.x - 1, current.y, current.z));
			}
		if(current.y != m_resolution.y - 1)
			if(xyz(current.x, current.y +1, current.z) == old_value)
			{
				xyz(current.x, current.y+1, current.z) = target_value;
				Q.push(point3D<indextype>(current.x, current.y+1, current.z));
			}
		if(current.y != 0)
			if(xyz(current.x, current.y-1, current.z) == old_value)
			{
				xyz(current.x, current.y-1, current.z) = target_value;
				Q.push(point3D<indextype>(current.x, current.y-1, current.z));
			}
		if(current.z != m_resolution.z - 1)
			if(xyz(current.x, current.y, current.z+1) == old_value)
			{
				xyz(current.x, current.y, current.z+1) = target_value;
				Q.push(point3D<indextype>(current.x, current.y, current.z+1));
			}
		if(current.z != 0)
			if(xyz(current.x, current.y, current.z-1) == old_value)
			{
				xyz(current.x, current.y, current.z-1) = target_value;
				Q.push(point3D<indextype>(current.x, current.y, current.z-1));
			}

	}

}

template <class T>
void rtsImplicit3D<T>::FloodFill26(int x, int y, int z, T target_value)
{
	T old_value = xyz(x, y, z);
	if(target_value == old_value)
		return;

	queue<point3D<indextype>> Q;
	point3D<indextype> current(x, y, z);
	point3D<indextype> next;
	Q.push(current);
	indextype u, v, w;
	while(!Q.empty())
	{
		current = Q.front();
		if(xyz(current.x, current.y, current.z) == old_value)
			xyz(current.x, current.y, current.z) = target_value;
		Q.pop();
		for(u=-1; u<=1; u++)
			for(v=-1; v<=1; v++)
				for(w=-1; w<=1; w++)
				{
					next.x = current.x + u;
					next.y = current.y + v;
					next.z = current.z + w;

					if(next.x >= 0 && next.x < m_resolution.x &&
						next.y >= 0 && next.y < m_resolution.y &&
						next.z >= 0 && next.z < m_resolution.z &&
						xyz(next.x, next.y, next.z) == old_value)
						{
							xyz(next.x, next.y, next.z) = target_value;
							Q.push(next);
						}
				}
	}

				
}

template <class T>
void rtsImplicit3D<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>
vector<point3D<indextype>> rtsImplicit3D<T>::getEdgeNodes(T isovalue, bool protrusions = true)
{
	vector<point3D<indextype>> result;
	indextype x, y, z;
	int neighbors;
	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(xyz(x, y, z) >= isovalue)
				{
					neighbors = Neighbors26(x, y, z, isovalue);
					if(protrusions == false && neighbors < 1)
						continue;
					if(neighbors < 26)
						result.push_back(point3D<indextype>(x, y, z));
				}
			}

	return result;

}

template <class T>
unsigned int rtsImplicit3D<T>::BackgroundComponents6(indextype x, indextype y, indextype z, T threshold, int n = 18)
{
	/**
	This function computes the number of 6-connected background components in the local region of (x, y, z).
	This computation is performed by testing all 6 possible voxels that can connect to the specified node.  If
	a background node is found, the entire background component associated with that node is filled and the counter
	is incremented by 1.  The value n specifies the connectivity domain for the flood fill.
	The definition of background components is that specified by He, Kischell, Rioult and Holmes.
	**/

	//see if there is at least one BG component
	if(Neighbors6(x, y, z, threshold) == 6)
		return 0;

	
	//retrieve the local region of the function
	rtsImplicit3D<T> local(3, 3, 3);
	point3D<indextype> corner(x-1, y-1, z-1);
	indextype u, v, w;
	for(u=0; u<3; u++)
		for(v=0; v<3; v++)
			for(w=0; w<3; w++)
				local(u, v, w) = xyz(corner.x + u, corner.y + v, corner.z + w);

	//threshold the background to find inside/outside points
	local.Binary(threshold, 1);
	//fill points that are not in the connectivity domain
	if(n == 18)
	{
		local(0, 0, 0) = 1;
		local(0, 0, 2) = 1;
		local(0, 2, 0) = 1;
		local(0, 2, 2) = 1;
		local(2, 0, 0) = 1;
		local(2, 0, 2) = 1;
		local(2, 2, 0) = 1;
		local(2, 2, 2) = 1;
	}
	//local.toConsole();

	//search all 6 possible connected points.  If a background node is found, fill the component
	unsigned int components = 0;
	if(local(0, 1, 1) == 0)
	{
		components++;
		local.FloodFill6(0, 1, 1, 1);
	}
	if(local(2, 1, 1) == 0)
	{
		components++;
		local.FloodFill6(2, 1, 1, 1);
	}
	if(local(1, 0, 1) == 0)
	{
		components++;
		local.FloodFill6(1, 0, 1, 1);
	}
	if(local(1, 2, 1) == 0)
	{
		components++;
		local.FloodFill6(1, 2, 1, 1);
	}
	if(local(1, 1, 0) == 0)
	{
		components++;
		local.FloodFill6(1, 1, 0, 1);
	}
	if(local(1, 1, 2) == 0)
	{
		components++;
		local.FloodFill6(1, 1, 2, 1);
	}

	return components;
}

template <class T>
rtsImplicit3D<T> rtsImplicit3D<T>::Project2D()
{
	/**
	This function projects the entire 3D function onto a 2D function along the z-axis.
	**/
	rtsImplicit3D<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;
}

template <class T>
void rtsImplicit3D<T>::Erode(T isovalue, T fill_value)
{
	vector<point3D<indextype>> border_nodes;	//get the border nodes for the image
	indextype x, y, z;
	int condition;
	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(xyz(x, y, z) >= isovalue && BackgroundComponents6(x, y, z, isovalue) == 1)
				{
					condition = 0;
					//now find the border pairs.  A border point must meet two of the following conditions to be a border pair.
					//south border: s(p) is background
					if(xyz(x, y-1, z) < isovalue)
						condition++;
					//north border: n(p) is background, s(p) and s(s(p)) are foreground
					if(xyz(x, y+1, z) < isovalue && xyz(x, y-1, z) >= isovalue && xyz(x, y-2, z) >= isovalue)
						condition++;
					//west border: w(p) is background
					if(xyz(x-1, y, z) < isovalue)
						condition++;
					//east border: e(p) is background, w(p) and w(w(p)) are foreground
					if(xyz(x+1, y, z) < isovalue && xyz(x-1, y, z) >= isovalue && xyz(x-2, y, z) >= isovalue)
						condition++;
					//up border: u(p) is background
					if(xyz(x, y, z-1) < isovalue)
						condition++;
					//down border: d(p) is background, u(p) and u(u(p)) are foreground
					if(xyz(x, y, z+1) < isovalue && xyz(x, y, z-1) >= isovalue && xyz(x, y, z-2) >= isovalue)
						condition++;
					
					if(condition > 1)
						border_nodes.push_back(point3D<indextype>(x, y, z));
				}
	cout<<"Number of border nodes: "<<border_nodes.size()<<endl;
	
	vector<point3D<indextype>>::iterator i;
	for(i=border_nodes.begin(); i!= border_nodes.end(); i++)
		xyz((*i).x, (*i).y, (*i).z) = fill_value;


}

template <class T>
void rtsImplicit3D<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 rtsImplicit3D<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>
void rtsImplicit3D<T>::ClampMin(T min, T value)
{
	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] = value;
}

template <class T>
void rtsImplicit3D<T>::MedianFilter(int dist_x, int dist_y, int dist_z, double factor = 0.5)
{
	rtsImplicit3D<T> result = (*this);
	indextype x, y, z;
	indextype min_x, min_y, min_z;
	indextype max_x, max_y, max_z;
	indextype u, v, w;
	vector<T> region;
	for(x=0; x<m_resolution.x; x++)
		for(y=0; y<m_resolution.y; y++)
			for(z=0; z<m_resolution.z; z++)
			{
				region.clear();
				min_x = x - dist_x;
				min_y = y - dist_y;
				min_z = z - dist_z;
				max_x = x + dist_x;
				max_y = y + dist_y;
				max_z = z + dist_z;

				for(u=min_x; u<=max_x; u++)
					for(v=min_y; v<=max_y; v++)
						for(w=min_z; w<=max_z; w++)
						{
							region.push_back(xyz(u, v, w));
						}
				sort(region.begin(), region.end());
				result(x, y, z) = region[(int)(region.size()*factor)];
			}
	(*this) = result;
}

template <class T>
unsigned int rtsImplicit3D<T>::Thin(T isovalue)
{
	/**
	This function computes the skeleton of an isosurface embedded in the implicit function and
	described by the "isovalue" parameter.
	**/

	vector<point3D<indextype>> border_nodes;	//get the border nodes for the image
	indextype x, y, z;
	int condition;
	for(x=0; x<m_resolution.x; x++)
		for(y=0; y<m_resolution.y; y++)
			for(z=0; z<m_resolution.z; z++)
				//find the border nodes
				if(xyz(x, y, z) >= isovalue && BackgroundComponents6(x, y, z, isovalue) == 1 && Neighbors26(x, y, z, isovalue) != 1)
				{
					condition = 0;
					//now find the border pairs.  A border point must meet two of the following conditions to be a border pair.
					//south border: s(p) is background
					if(xyz(x, y-1, z) < isovalue)
						condition++;
					//north border: n(p) is background, s(p) and s(s(p)) are foreground
					if(xyz(x, y+1, z) < isovalue && xyz(x, y-1, z) >= isovalue && xyz(x, y-2, z) >= isovalue)
						condition++;
					//west border: w(p) is background
					if(xyz(x-1, y, z) < isovalue)
						condition++;
					//east border: e(p) is background, w(p) and w(w(p)) are foreground
					if(xyz(x+1, y, z) < isovalue && xyz(x-1, y, z) >= isovalue && xyz(x-2, y, z) >= isovalue)
						condition++;
					//up border: u(p) is background
					if(xyz(x, y, z-1) < isovalue)
						condition++;
					//down border: d(p) is background, u(p) and u(u(p)) are foreground
					if(xyz(x, y, z+1) < isovalue && xyz(x, y, z-1) >= isovalue && xyz(x, y, z-2) >= isovalue)
						condition++;
					
					if(condition > 1)
						border_nodes.push_back(point3D<indextype>(x, y, z));
				}
	cout<<"Number of border nodes: "<<border_nodes.size()<<endl;

	//determine if each edge node can be removed without changing the topology of the model
	//declare some initial variables
	rtsImplicit3D<T> local(3, 3, 3);	//store the region local to the current voxel
	int nodes_before, nodes_after;		//number of neighboring nodes before and after the filling operation
	point3D<indextype> fill_start;
	vector<point3D<indextype>>::iterator i;

	/*
	Here we determine if a point can be removed by looking at the number of foreground connected
	components in the local region.  If there is more than one connected component
	*/
	unsigned int removed = 0;
	for(i=border_nodes.begin(); i<border_nodes.end(); i++)
	{
		//get the local region around the current point
		for(x=-1; x<=1; x++)
			for(y=-1; y<=1; y++)
				for(z=-1; z<=1; z++)
					local(x+1, y+1, z+1) = xyz((*i).x + x, (*i).y + y, (*i).z + z);

		//deal with the degenerate case of all four sides being internal
		//if(local(0, 1, 0) >= isovalue && local(1, 0, 0) >= isovalue && local(1, 2, 0) >= isovalue && local(2, 1, 0) >= isovalue)
		//	continue;
		local(1, 1, 1) = 0;				//remove the center voxel
		local.Binary(isovalue, 1);
		nodes_before = local.Neighbors26(1, 1, 1, 1);
		//if(nodes_before == 1)			//prevent reducing ends
		//	continue;

		//find an interior voxel to fill
		for(x=0; x<3; x++)
			for(y=0; y<3; y++)
				for(z=0; z<3; z++)
					if(local(x, y, z) > 0)
						fill_start = point3D<indextype>(x, y, z);

		//fill the local region
		local.FloodFill26(fill_start.x, fill_start.y, fill_start.z, 2);
		//get the number of filled neighbors
		nodes_after = local.Neighbors26(1, 1, 1, 2);
		if(nodes_after == nodes_before)
		{
			xyz((*i).x, (*i).y, (*i).z) = 0;
			removed++;
			//cout<<"removed"<<endl;
		}
		//else
		//{
			/*for(x=-1; x<=1; x++)
			for(y=-1; y<=1; y++)
				for(z=-1; z<=1; z++)
					local(x+1, y+1, z+1) = xyz((*i).x + x, (*i).y + y, (*i).z + z);
			local.toConsole();
			cout<<"not removed:  "<<nodes_before<<"  "<<nodes_after<<endl;
			xyz((*i).x, (*i).y, (*i).z) = 100;
			char c;
			cin>>c;*/
		//}
	}
	return removed;
}

/*Shape functions*/
/*These functions create implicit shapes in the function.
Generally, the shape is based on the parameterization.*/
template <class T>
void rtsImplicit3D<T>::Sphere(double center_i, double center_j, double center_k, double radius, T in_value)
{
	point3D<double> center(center_i, center_j, center_k);
	//for each point in the function
	indextype x, y, z;
	double radius_squared = radius*radius;
	vector3D<double> point_to_point;
	point3D<double> result;
	double distance_squared;

	for(x=0; x<m_resolution.x; x++)
		for(y=0; y<m_resolution.y; y++)
			for(z=0; z<m_resolution.z; z++)
			{
				//get the parameterized value
				result = getParameter(x, y, z);
				//find the distance between the center of the sphere and the resulting point
				//double distance = (result - center).Length();
				point_to_point = result - center;
				distance_squared = point_to_point*point_to_point;
				//if the distance is less than the radius, the point is inside the sphere
				if(distance_squared < radius_squared)
					xyz(x, y, z) = in_value;
			}
			
}


template <class T>
void rtsImplicit3D<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;
	}

}


#endif