codebase / stimlib

#ifndef RTS_DTGRID_3D_H
#define RTS_DTGRID_3D_H

#include <vector>
#include <iostream>
using namespace std;

#include "rtsDTGrid2D.h"

struct Coord3D
{
	int x1;
	int x2;
	int x3;

	bool operator==(Coord3D &rhs)
	{
		if( (x1 == rhs.x1) && (x2 == rhs.x2) && (x3 == rhs.x3) )
			return true;
		return false;
	}

	bool operator<(Coord3D &rhs)
	{
		if(x1 < rhs.x1)
			return true;
		else if( (x1 == rhs.x1) && (x2 < rhs.x2) )
			return true;
		else if( (x1 == rhs.x1) && (x2 == rhs.x2) && (x3 < rhs.x3) )
			return true;
		return false;

	}

};




/**************DT GRID**************************/
template<typename T>
class rtsDTGrid3D
{
private:
	//main arrays
	vector<T> value;
	vector<ConnectedComponent> conn;
	rtsDTGrid2D<IndexPair> proj2D;
	Coord2D current_column;
	bool randomIndex(rtsDTGrid2D<IndexPair>::iterator &iter2D, int &v_i, int &c_i, int x1, int x2, int x3);

	

	//variables to keep track of insertion
	int max_coord;
	bool grid_insertion_started;
	bool column_insertion_started;


public:
	rtsDTGrid3D<T>()
	{
		grid_insertion_started = false;
		column_insertion_started = false;
		proj2D.background.toConn = -1;
		max_coord = 0;
	}

	bool push(int x1, int x2, int x3, T v);
	T random(int x1, int x2, int x3);
	T background;

	
	T& back();
	
	void print();
	void operator=(T rhs);
	void getBounds(int &min_x1, int &min_x2, int &min_x3, int &max_x1, int &max_x2, int &max_x3);
	int getNumNodes(){return value.size();}
	void insert(rtsDTGrid3D<T> toInsert);

	//arithmetic
	rtsDTGrid3D<T> operator+(rtsDTGrid3D<T> &rhs);
	rtsDTGrid3D<T> operator-(rtsDTGrid3D<T> &rhs);


	//dilation
	friend class ColumnUnion;
	void dilate(int H);

	/*
	
	
	


	//other types of iteration
	iterator begin_pn();
	iterator end_pn();
	iterator begin_np();
	iterator end_np();
	*/

	//iterator
	class iterator;
	friend class iterator;
	class stencil_iterator;
	//iterator randomIterator(int x1, int x2);
	iterator begin();
	iterator before();
	iterator end();
	iterator after();

	//other types of iteration
	iterator begin_ppn();
	iterator begin_nnp();
	iterator begin_pnp();
	iterator begin_npn();
	iterator begin_pnn();
	iterator begin_npp();


	iterator end_ppn(){return begin_nnp();}
	iterator end_nnp(){return begin_ppn();}
	iterator end_pnp(){return begin_npn();}
	iterator end_npn(){return begin_pnp();}
	iterator end_pnn(){return begin_npp();}
	iterator end_npp(){return begin_pnn();}
};


/**********ITERATOR***********************/
template<typename T>
class rtsDTGrid3D<T>::iterator
{
	friend class rtsDTGrid3D;
	rtsDTGrid3D<T>* parent;
	rtsDTGrid2D<IndexPair>::iterator loc2D;
	int iv;
	int ic;
	int x3;

public:
	

	T Value(){return parent->value[iv];}
	int X1(){return loc2D.X1();}
	int X2(){return loc2D.X2();}
	int X3(){return x3;}
	Coord3D Coord()
	{
		Coord3D result;
		result.x1 = X1();
		result.x2 = X2();
		result.x3 = X3();
		return result;
	}
	iterator(){parent = NULL;}
	void SetValue(T value){parent->value[iv] = value;}
	

	void ppp()
	{
		//increment the value
		iv++;
		//if this exceeds the length of the value array, we are at the end of the grid
		if(iv == parent->value.size())
		{
			(*this) = parent->after();
			return;
		}

		//increment the current coordinate
		x3++;
		//if we are outside of the current connected component
		if(x3 > parent->conn[ic].coordMax)
		{
			//move to the next connected component
			ic++;
			//if this is the last connected component in the column
			if(ic == loc2D.Value().toConn + loc2D.Value().numConn)
				loc2D++;
			//if there are no more connected components
			if(ic == parent->conn.size())
			{
				//we're at the end, return end
				(*this) = parent->end();
				return;
			}
			x3 = parent->conn[ic].coordMin;
		}
	}
	void ppn()
	{
		//decrement the value
		iv--;

		//decrement the current coordinate
		x3--;
		//if we are outside of the current connected component
		if(x3 < parent->conn[ic].coordMin)
		{
			//move to the previous connected component
			ic--;
			//if this is before the first connected component in the column
			if(ic < loc2D.Value().toConn)
			{
				//increment to the next column
				loc2D++;
				//set the connected component to the last one in the column
				ic = loc2D.Value().toConn + loc2D.Value().numConn - 1;
				int num_values = parent->conn[ic].coordMax - parent->conn[ic].coordMin + 1;
				iv = parent->conn[ic].toValue + num_values - 1;
			}

			x3 = parent->conn[ic].coordMax;
		}


	}

	void nnp()
	{
		//increment the value
		iv++;

		//increment the current coordinate
		x3++;
		//if we are outside of the current connected component
		if(x3 > parent->conn[ic].coordMax)
		{
			//move to the next connected component
			ic++;
			//if this is after the last connected component in the column
			if(ic == loc2D.Value().toConn + loc2D.Value().numConn)
			{
				//decrement to the previous column
				loc2D--;
				//set the connected component to the first one in the column
				ic = loc2D.Value().toConn;
				iv = parent->conn[ic].toValue;
			}

			x3 = parent->conn[ic].coordMin;
		}


	}

	
	void nnn()
	{
		//decrement the value
		iv--;
		//if this is less than 0, we are at the beginning of the grid
		if(iv < 0)
		{
			(*this) = parent->before();
			return;
		}

		//decrement the current coordinate
		x3--;
		//if we are outside of the current connected component
		if(x3 < parent->conn[ic].coordMin)
		{
			//move to the previous connected component
			ic--;
			//if this is the first connected component in the column
			if(ic < loc2D.Value().toConn)
				loc2D--;
			//if there are no more connected components
			if(ic < 0)
			{
				//we're at the beginning, return begin
				(*this) = parent->before();
				return;
			}
			x3 = parent->conn[ic].coordMax;
		}

	}

	void pnp()
	{
		//increment the value
		iv++;

		//increment the current coordinate
		x3++;
		//if we are outside of the current connected component
		if(x3 > parent->conn[ic].coordMax)
		{
			//move to the next connected component
			ic++;
			//if this is after the last connected component in the column
			if(ic == loc2D.Value().toConn + loc2D.Value().numConn)
			{
				//decrement to the previous column
				loc2D.pn();
				//set the connected component to the first one in the column
				ic = loc2D.Value().toConn;
				iv = parent->conn[ic].toValue;
			}

			x3 = parent->conn[ic].coordMin;
		}


	}
	void npn()
	{
		//decrement the value
		iv--;

		//decrement the current coordinate
		x3--;
		//if we are outside of the current connected component
		if(x3 < parent->conn[ic].coordMin)
		{
			//move to the previous connected component
			ic--;
			//if this is before the first connected component in the column
			if(ic < loc2D.Value().toConn)
			{
				//increment to the next column
				loc2D.np();
				//set the connected component to the last one in the column
				ic = loc2D.Value().toConn + loc2D.Value().numConn - 1;
				int num_values = parent->conn[ic].coordMax - parent->conn[ic].coordMin + 1;
				iv = parent->conn[ic].toValue + num_values - 1;
			}

			x3 = parent->conn[ic].coordMax;
		}


	}
	void pnn()
	{
		//decrement the value
		iv--;

		//decrement the current coordinate
		x3--;
		//if we are outside of the current connected component
		if(x3 < parent->conn[ic].coordMin)
		{
			//move to the previous connected component
			ic--;
			//if this is before the first connected component in the column
			if(ic < loc2D.Value().toConn)
			{
				//increment to the next column
				loc2D.pn();
				//set the connected component to the last one in the column
				ic = loc2D.Value().toConn + loc2D.Value().numConn - 1;
				int num_values = parent->conn[ic].coordMax - parent->conn[ic].coordMin + 1;
				iv = parent->conn[ic].toValue + num_values - 1;
			}

			x3 = parent->conn[ic].coordMax;
		}


	}
	void npp()
	{
		//increment the value
		iv++;

		//increment the current coordinate
		x3++;
		//if we are outside of the current connected component
		if(x3 > parent->conn[ic].coordMax)
		{
			//move to the next connected component
			ic++;
			//if this is after the last connected component in the column
			if(ic == loc2D.Value().toConn + loc2D.Value().numConn)
			{
				//decrement to the previous column
				loc2D.np();
				//set the connected component to the first one in the column
				ic = loc2D.Value().toConn;
				iv = parent->conn[ic].toValue;
			}

			x3 = parent->conn[ic].coordMin;
		}


	}

	/*
	void pn()
	{
		//for the most part we will be going backwards through the array
		//so first decrement iv
		iv--;
		x2--;
		//if iv is less than the current connected component
		if(iv < parent->conn[ic].toValue)
		{
			//go to the previous connected component
			ic--;
			//if we are before the first connected component in the column
			if(ic < loc1D.Value().toConn)
			{
				//increment the 1D iterator
				loc1D++;
				//reset ic to the last component of the new column
				ic = loc1D.Value().toConn + loc1D.Value().numConn - 1;
				//find the new value identifier
				iv = parent->conn[ic].toValue + (parent->conn[ic].coordMax - parent->conn[ic].coordMin);
			}
			//compute the currect coordinate
			x2 = parent->conn[ic].coordMax;
		}
	}

	void np()
	{
		//for the most part we will be going forward through the grid
		//increment iv
		iv++;
		x2++;
		
		//if we are outside of the current connected component
		if(x2 > parent->conn[ic].coordMax)
		{
			//move to the next connected component
			ic++;
			//if this is the last connected component in the column
			if(ic == loc1D.Value().toConn + loc1D.Value().numConn)
			{
				loc1D--;
				ic = loc1D.Value().toConn;
				iv = parent->conn[ic].toValue;
			}

			x2 = parent->conn[ic].coordMin;
		}

	}
	

	

	friend bool operator<(iterator &left, iterator &right)
	{
		if(left.iv < right.iv)
			return true;
		return false;
	}
	friend bool operator<=(iterator &left, iterator &right)
	{
		if(left.iv <= right.iv)
			return true;
		return false;
	}
	
	*/
	
	void increment_until(int p1, int p2, int p3)
	{
		while((*this) != parent->end())
		{
			if(X1() > p1)
				return;
			if(X1() == p1 && X2() > p2)
				return;
			else if(X1() == p1 && X2() == p2 && X3() >= p3)
				return;

			ppp();
		}
	}
	void operator++(){ppp();}
	void operator--(){nnn();}

	//boolean operators for comparing iterators
	bool operator==(iterator &rhs)
	{
		if(parent == rhs.parent && iv == rhs.iv)
			return true;
	
		//if(loc2D == rhs.loc2D && x3 == rhs.x3)
		//	return true;
		return false;
	}
	bool operator!=(iterator &rhs){return !((*this) == rhs);}
	bool operator<(iterator rhs)
	{
		if(parent == rhs.parent && iv < rhs.iv)
			return true;
		//if(loc2D < rhs.loc2D && x3 < rhs.x3)
		//	return true;
		return false;

	}
};

/**************ITERATOR METHODS IN DT GRID*******************/
template<typename T>
typename rtsDTGrid3D<T>::iterator rtsDTGrid3D<T>::begin()
{
	//if the grid is empty, return an iterator to "after"
	if(value.size() == 0)
		return after();

	iterator result;
	result.parent = this;
	result.ic = 0;
	result.iv = 0;
	result.x3 = conn[0].coordMin;
	result.loc2D = proj2D.begin();

	return result;
}
template<typename T>
typename rtsDTGrid3D<T>::iterator rtsDTGrid3D<T>::end()
{
	//if the grid is empty, return after()
	if(value.size() == 0)
		return after();
	iterator result;
	result.parent = this;
	result.ic = conn.size() - 1;
	result.iv = value.size() - 1;
	result.x3 = conn[result.ic].coordMax;
	result.loc2D = proj2D.end();
	return result;
}
template<typename T>
typename rtsDTGrid3D<T>::iterator rtsDTGrid3D<T>::after()
{
	iterator result;
	result.parent = this;
	result.ic = conn.size() - 1;
	result.iv = value.size();
	//result.x2 = conn[result.ic].coordMax;
	result.loc2D = proj2D.after();
	return result;
}


template<typename T>
typename rtsDTGrid3D<T>::iterator rtsDTGrid3D<T>::before()
{
	if(value.size() == 0)
		return after();

	iterator result;
	result.parent = this;
	result.ic = 0;
	result.iv = -1;
	//result.x2 = conn[0].coordMin;
	result.loc2D = proj2D.before();

	return result;
}

template<typename T>
typename rtsDTGrid3D<T>::iterator rtsDTGrid3D<T>::begin_ppn()
{
	//if the grid is empty, return an iterator to "after"
	if(value.size() == 0)
		return after();

	iterator result;
	result.parent = this;
	result.loc2D = proj2D.begin();
	//find the index of the last connected component in the first column
	int last_conn = result.loc2D.Value().toConn + result.loc2D.Value().numConn - 1;
	result.ic = last_conn;
	result.iv = conn[last_conn].toValue + conn[last_conn].coordMax - conn[last_conn].coordMin;
	result.x3 = conn[last_conn].coordMax;
	

	return result;
}

template<typename T>
typename rtsDTGrid3D<T>::iterator rtsDTGrid3D<T>::begin_nnp()
{
	//if the grid is empty, return an iterator to "after"
	if(value.size() == 0)
		return after();

	iterator result;
	result.parent = this;
	result.loc2D = proj2D.end();
	//find the index of the first connected component in the last column
	int first_conn = result.loc2D.Value().toConn;
	result.ic = first_conn;
	result.iv = conn[first_conn].toValue;
	result.x3 = conn[first_conn].coordMin;
	

	return result;
}
template<typename T>
typename rtsDTGrid3D<T>::iterator rtsDTGrid3D<T>::begin_pnp()
{
	//if the grid is empty, return an iterator to "after"
	if(value.size() == 0)
		return after();

	iterator result;
	result.parent = this;
	result.loc2D = proj2D.begin_pn();
	//find the index of the first connected component in the column
	int first_conn = result.loc2D.Value().toConn;
	result.ic = first_conn;
	result.iv = conn[first_conn].toValue;
	result.x3 = conn[first_conn].coordMin;
	

	return result;
}
template<typename T>
typename rtsDTGrid3D<T>::iterator rtsDTGrid3D<T>::begin_npn()
{
	//if the grid is empty, return an iterator to "after"
	if(value.size() == 0)
		return after();

	iterator result;
	result.parent = this;
	result.loc2D = proj2D.begin_np();
	//find the index of the last connected component in the column
	int last_conn = result.loc2D.Value().toConn + result.loc2D.Value().numConn - 1;
	result.ic = last_conn;
	result.iv = conn[last_conn].toValue + conn[last_conn].coordMax - conn[last_conn].coordMin;
	result.x3 = conn[last_conn].coordMax;
	

	return result;
}
template<typename T>
typename rtsDTGrid3D<T>::iterator rtsDTGrid3D<T>::begin_pnn()
{
	//if the grid is empty, return an iterator to "after"
	if(value.size() == 0)
		return after();

	iterator result;
	result.parent = this;
	result.loc2D = proj2D.begin_pn();
	//find the index of the last connected component in the column
	int last_conn = result.loc2D.Value().toConn + result.loc2D.Value().numConn - 1;
	result.ic = last_conn;
	result.iv = conn[last_conn].toValue + conn[last_conn].coordMax - conn[last_conn].coordMin;
	result.x3 = conn[last_conn].coordMax;
	

	return result;
}
template<typename T>
typename rtsDTGrid3D<T>::iterator rtsDTGrid3D<T>::begin_npp()
{
	//if the grid is empty, return an iterator to "after"
	if(value.size() == 0)
		return after();

	iterator result;
	result.parent = this;
	result.loc2D = proj2D.begin_np();
	//find the index of the first connected component in the column
	int first_conn = result.loc2D.Value().toConn;
	result.ic = first_conn;
	result.iv = conn[first_conn].toValue;
	result.x3 = conn[first_conn].coordMin;
	

	return result;
}
/**************DT GRID METHODS**************************/
template<typename T>
bool rtsDTGrid3D<T>::push(int x1, int x2, int x3, T v)
{
	//run this code if we have to start a new column.  This happens when:
	//(a) we have just inserted the first value into the grid
	//(b) the insertion takes place in a different column
	if(grid_insertion_started == false || x1 != current_column.x1 || x2 != current_column.x2)
	{
		//assume that a new column is being created
		//create the column in proj1D
		IndexPair newPair;
		newPair.toConn = conn.size();
		//newPair.toValue = value.size();
		newPair.numConn = 0;

		//if this insertion throws an error, the value was inserted incorrectly
		//the error will get thrown by DTGrid2D, so just return
		if(!proj2D.push(x1, x2, newPair))
		{
			cout<<"Out-of-order insertion in D = 3: X1 = "<<x1<<" X2 = "<<x2<<" X3 = "<<x3<<endl;
			return false;
		}

		column_insertion_started = false;
		grid_insertion_started = true;
	}
	//If there is no new column, we still have to check to make sure the value is inserted
	//in the correct order in X3:
	else
	{
		if(column_insertion_started && x3 <= max_coord)
		{
			cout<<"Out-of-order insertion in D = 3: X1 = "<<x1<<" X2 = "<<x2<<" X3 = "<<x3<<endl;
			return false;
		}
	}
	

	//run this code if we have to start a new connected component.  This happens when:
	//(a) We insert the first value into the column
	//(b) There is empty space between the last insertion and this one
	if(column_insertion_started == false || x3 > (max_coord + 1))
	{
		//start a new connected component
		ConnectedComponent new_conn;
		new_conn.toValue = value.size();
		new_conn.coordMin = x3;
		new_conn.coordMax = x3;
		conn.push_back(new_conn);
		//increment the number of connected components for the column in X1
		proj2D.back().numConn++;
		column_insertion_started = true;
	}

	//insert the value into the grid:
	//(a) Insert the value at the end of the coord array
	//(b) Increment the coordMax value of the current connected component
	//(c) Change the maximum inserted coordinate to the new value
	value.push_back(v);
	conn[conn.size() - 1].coordMax = x3;
	//change max_coord to the new coordinate
	max_coord = x3;
	//set the current column
	current_column.x1 = x1;
	current_column.x2 = x2;
	return true;
}


template<typename T>
bool rtsDTGrid3D<T>::randomIndex(rtsDTGrid2D<IndexPair>::iterator &iter2D, int &v_i, int &c_i, int x1, int x2, int x3)
{
	//search along the 2nd dimension
	//get an iterator
	iter2D = proj2D.randomIterator(x1, x2);

	//if the after() iterator is returned, exit
	if(iter2D == proj2D.after())
	{
		c_i = conn.size();
		v_i = value.size();
		return false;
	}

	IndexPair i = iter2D.Value();

	//if the column does not exist
	if(iter2D.X1() != x1 || iter2D.X2() != x2)
	{
		//get the appropriate values from the 2D iterator
		//(which points to the next valid column)
		c_i = i.toConn;
		v_i = conn[c_i].toValue;
		return false;
	}
	
	//otherwise perform the binary search within the column
	int high, low, mid;
	low = i.toConn;
	high = i.toConn + i.numConn - 1;

	do
	{
		mid = low + (high - low)/2;
		if(x3 > conn[mid].coordMax)
			low = mid + 1;
		else if(x3 < conn[mid].coordMin)
			high = mid - 1;
		else break;
	}
	while(low <= high);

	//at this point, mid is either at the appropriate connected component,
	//or x2 is not in the grid
	if(x3 >= conn[mid].coordMin && x3 <= conn[mid].coordMax)
	{
		int offset = x3 - conn[mid].coordMin;
		c_i = mid;
		v_i = conn[mid].toValue + offset;
		return true;
	}
	else
	{
		c_i = mid;
		v_i = conn[c_i].toValue;
		return false;
	}
}

template<typename T>
T rtsDTGrid3D<T>::random(int x1, int x2, int x3)
{
	int v_i, c_i;
	rtsDTGrid2D<IndexPair>::iterator iter2D;
	if(randomIndex(iter2D, v_i, c_i, x1, x2, x3))
		return value[v_i];
	else
		return background;
}



template<typename T>
T& rtsDTGrid3D<T>::back()
{
	return value[value.size()-1];
}
template<typename T>
void rtsDTGrid3D<T>::print()
{
	rtsDTGrid3D<T>::iterator i;
	for(i = begin(); i != after(); i++)
	{
		cout<<i.X1()<<","<<i.X2()<<","<<i.X3()<<":"<<i.Value()<<endl;
	}
	

}
template<typename T>
void rtsDTGrid3D<T>::operator =(T rhs)
{
	for(int i=0; i<value.size(); i++)
		value[i] = rhs;
}

template<typename T>
void rtsDTGrid3D<T>::getBounds(int &min_x1, int &min_x2, int &min_x3, int &max_x1, int &max_x2, int &max_x3)
{
	//if the grid is empty, return an empty bounding volume
	if(value.size() == 0)
	{
		min_x1 = min_x2 = min_x3 = max_x1 = max_x2 = max_x3 = 0;
		return;
	}

	//get the min and max values for the 1D grid (x1 coordinate)
	proj2D.getBounds(min_x1, min_x2, max_x1, max_x2);

	//initialize the min and max values
	min_x3 = conn[0].coordMin;
	max_x3 = conn.back().coordMax;

	//iterate through all columns finding the smallest and largest coordinate values
	rtsDTGrid2D<IndexPair>::iterator i;
	IndexPair col;
	for(i=proj2D.begin(); i!=proj2D.after(); i++)
	{
		col = i.Value();
		if(conn[col.toConn].coordMin < min_x3)
			min_x3 = conn[col.toConn].coordMin;
		if(conn[col.toConn + col.numConn - 1].coordMax > max_x3)
			max_x3 = conn[col.toConn + col.numConn - 1].coordMax;
	}
		

}
template<typename T>
void rtsDTGrid3D<T>::insert(rtsDTGrid3D<T> toInsert)
{
	//create source and destination iterators
	rtsDTGrid3D<T>::iterator source = toInsert.begin();
	rtsDTGrid3D<T>::iterator dest = begin();

	for(source = toInsert.begin(); source != toInsert.after(); source++)
	{
		//move the destination iterator to the current source position
		dest.increment_until(source.X1(), source.X2(), source.X3());
		//cout<<"source: "<<source.X1()<<" "<<source.X2()<<": "<<source.Value()<<endl;
		//cout<<"dest: "<<dest.X1()<<" "<<dest.X2()<<": "<<dest.Value()<<endl;
		//if the position exists in dest
		if(dest.X1() == source.X1() && dest.X2() == source.X2() && dest.X3() == source.X3())
			dest.SetValue(source.Value());
		//cout<<"dest: "<<dest.X1()<<" "<<dest.X2()<<": "<<dest.Value()<<endl;
	}

}

template<typename T>
void rtsDTGrid3D<T>::dilate(int H)
{
	//if the grid is empty, return unchanged
	if(value.size() == 0)
		return;

	ColumnUnion CUqueue;
	CUqueue.ConnectToGrid(&conn);


	//dilate the N-1 DT-Grid constituent
	rtsDTGrid2D<IndexPair> dilated_proj2D = proj2D;
	
	//an empty pair has a number of connected components equal to zero
	//this should not happen in reality and can therefore be used to check for new columns
	IndexPair empty_pair;
	empty_pair.numConn = 0;

	//set the background node to the empty node and dilate the 2D projection
	dilated_proj2D.background = empty_pair;
	dilated_proj2D.dilate(H);

	//create a new DT Grid that will replace this one
	rtsDTGrid3D<T> new_grid;
	new_grid.proj2D = dilated_proj2D;
	new_grid.proj2D.background.toConn = -1;

	//create iteratorDilate that iterates along the dilated N-1 grid
	rtsDTGrid2D<IndexPair>::stencil_iterator iteratorDilate;
	
	//create the template entrance nodes
	int n;
	for(n=-H; n<=H; n++)
		iteratorDilate.addPosition(n, H);
	//create the template exit nodes
	for(n=-H; n<=H; n++)
		iteratorDilate.addPosition(n, -H);
	
	//create an iterator to set the new proj2D values
	rtsDTGrid2D<IndexPair>::iterator iteratorNew;


	//variables for each iteration
	IndexPair new_pair;
	vector<ConnectedComponent>::iterator ccIterator;
	unsigned int numValues = 0;
	for(iteratorNew = new_grid.proj2D.begin(),
		iteratorDilate = dilated_proj2D.begin();
		iteratorDilate != dilated_proj2D.after();
		iteratorNew++,
		iteratorDilate++)
		{
			//if a column is entering the stencil
			for(n=0; n<=(2*H); n++)
				if(iteratorDilate.getValue(n).numConn)
					CUqueue.InsertColumn(iteratorDilate.getValue(n));
			
			//compute the union of all columns in the queue
			vector<ConnectedComponent> result = CUqueue.ComputeUnion(H);

			//compute the new IndexPair representing the column
			new_pair.toConn = new_grid.conn.size();
			new_pair.numConn = result.size();
			//store the index pair
			iteratorNew.SetValue(new_pair);

			//insert each of the connected components
			for(ccIterator = result.begin(); ccIterator!=result.end(); ccIterator++)
			{
				new_grid.conn.push_back(*ccIterator);
				new_grid.conn.back().toValue = numValues;
				numValues += (*ccIterator).coordMax - (*ccIterator).coordMin + 1;
			}

			//if a column is leaving the stencil
			for(n=(2*H + 1); n<=(4*H + 1); n++)
				if(iteratorDilate.getValue(n).numConn)
					CUqueue.RemoveColumn();
		}

	//allocate space for the new value array
	new_grid.value.resize(numValues, background);

	//copy the data from this grid into the new grid
	new_grid.insert(*this);

	//replace this grid with the new grid
	conn = new_grid.conn;
	value = new_grid.value;
	proj2D = new_grid.proj2D;

	
}



/***************ARITHMETIC********************************/
template<typename T>
rtsDTGrid3D<T> rtsDTGrid3D<T>::operator+(rtsDTGrid3D<T> &rhs)
{
	rtsDTGrid3D<T> result;

	//create an iterator for each DT Grid
	rtsDTGrid3D<T>::iterator left = begin();
	rtsDTGrid3D<T>::iterator right = rhs.begin();

	//iterate both until one iterator has hit after()
	while(left != after() && right != rhs.after())
	{
		//if the iterators are at the same coordinate
		if(left.Coord() == right.Coord())
		{
			//insert their sum into the new grid
			result.push(left.X1(), left.X2(), left.X3(), left.Value() + right.Value());
			//increment both
			left++; right++;
		}
		//add the lowest (lexicographically) value to the background, insert the result, and increment
		else if( (left.Coord() < right.Coord()) )
		{
			result.push(left.X1(), left.X2(), left.X3(), left.Value() + rhs.background);
			left++;
		}
		else if( (right.Coord() < left.Coord()) )
		{
			result.push(right.X1(), right.X2(), right.X3(), right.Value() + background);
			right++;
		}
	}

	//if the left iterator hasn't finished, iterate to finish it off
	while(left != after())
	{
		result.push(left.X1(), left.X2(), left.X3(), left.Value() + rhs.background);
		left++;
	}

	while(right != rhs.after())
	{
		result.push(right.X1(), right.X2(), right.X3(), right.Value() + rhs.background);
		right++;
	}


	return result;


}

template<typename T>
rtsDTGrid3D<T> rtsDTGrid3D<T>::operator-(rtsDTGrid3D<T> &rhs)
{
	rtsDTGrid3D<T> result;

	//create an iterator for each DT Grid
	rtsDTGrid3D<T>::iterator left = begin();
	rtsDTGrid3D<T>::iterator right = rhs.begin();

	//iterate both until one iterator has hit after()
	while(left != after() && right != rhs.after())
	{
		//if the iterators are at the same coordinate
		if(left.Coord() == right.Coord())
		{
			//insert their sum into the new grid
			result.push(left.X1(), left.X2(), left.X3(), left.Value() - right.Value());
			//increment both
			left++; right++;
		}
		//add the lowest (lexicographically) value to the background, insert the result, and increment
		else if( (left.Coord() < right.Coord()) )
		{
			result.push(left.X1(), left.X2(), left.X3(), left.Value() - rhs.background);
			left++;
		}
		else if( (right.Coord() < left.Coord()) )
		{
			result.push(right.X1(), right.X2(), right.X3(), background - right.Value());
			right++;
		}
	}

	//if the left iterator hasn't finished, iterate to finish it off
	while(left != after())
	{
		result.push(left.X1(), left.X2(), left.X3(), left.Value() - rhs.background);
		left++;
	}

	while(right != rhs.after())
	{
		result.push(right.X1(), right.X2(), right.X3(), background - right.Value());
		right++;
	}


	return result;


}
#endif