package otherSupport;

import oneLiners.OneLiners;
import binaryMatrixFile.BinaryMatrixFile;

/** 3D Phase unwrapping ported from C
//Taken from http://www.ljmu.ac.uk/GERI/90208.htm
// License is special, but is basically non-profit free to share.
//
// oliford

	//This program was written by Hussein Abdul-Rahman and Munther Gdeisat to program the three-dimensional phase unwrapper
	//entitled "Fast three-dimensional phase-unwrapping algorithm based on sorting by 
	//reliability following a noncontinuous path"
	//by  Hussein Abdul-Rahman, Munther A. Gdeisat, David R. Burton, and Michael J. Lalor, 
	//published in the Proceedings of SPIE - 
	//The International Society for Optical Engineering, Vol. 5856, No. 1, 2005, pp. 32-40
	//This program was written by Munther Gdeisat, Liverpool John Moores University, United Kingdom.
	//Date 31st August 2007
	//The wrapped phase volume is assumed to be of floating point data type. The resultant unwrapped phase volume is also of floating point type.
	//Read the data from the file frame by frame
	//The mask is of byte data type. 
	//When the mask is 255 this means that the voxel is valid 
	//When the mask is 0 this means that the voxel is invalid (noisy or corrupted voxel)
	//This program takes into consideration the image wrap around problem encountered in MRI imaging.

*/
public class PhaseUnwrap3D {
	
	boolean x_connectivity;
	boolean y_connectivity;
	boolean z_connectivity;
	int No_of_edges = 0;

	//VoxelInfo information
	private class VoxelInfo {
		//int x;					//x coordinate of the voxel
	    //int y;					//y coordinate
		//int z;					//z coordinate
	    int increment;				//No. of 2*pi to add to the voxel to unwrap it
	    int number_of_voxels_in_group;	        //No. of voxel in the voxel group
	    double value;				//value of the voxel
	    double reliability;
	    boolean input_mask;		//false voxel is masked. true is not masked
	    boolean extended_mask;		//false voxel is masked. true is not masked
	    int group;					//group No.
	    int new_group;
	    VoxelInfo head;		//pointer to the first voxel in the group in the linked list
	    VoxelInfo last;		//pointer to the last voxel in the group
	    VoxelInfo next;		//pointer to the next voxel in the group
	};

	//the EdgeInfo is the line that connects two voxels.
	//if we have S voxels, then we have S horizontal edges and S vertical edges
	private class EdgeInfo 	{    
		double reliab;			//reliabilty of the edge and it depends on the two voxels
		VoxelInfo pointer_1;		//pointer to the first voxel
		VoxelInfo pointer_2;		//pointer to the second voxel
	    int increment;			//No. of 2*pi to add to one of the voxels to unwrap it with respect to the second 
	}; 

	
	//---------------start quicker_sort algorithm --------------------------------
	//#define swap(x,y) {EdgeInfo t; t=x; x=y; y=t;}
	//#define order(x,y) if (x.reliab > y.reliab) swap(x,y)
	//#define o2(x,y) order(x,y)
	//#define o3(x,y,z) o2(x,y); o2(x,z); o2(y,z)
	
	
	private double find_pivot(EdgeInfo edges[], int left, int right)	{
		int a, b, c;

		a = left;
		b = left + (right - left) /2;
		c = right;		
		if (edges[a].reliab > edges[b].reliab){ int t = a; a = b; b = t; }
		if (edges[a].reliab > edges[c].reliab){ int t = a; a = c; c = t; }
		if (edges[b].reliab > edges[c].reliab){ int t = b; b = c; c = t; }
		
		if (edges[a].reliab < edges[b].reliab) {
			return edges[b].reliab;
		}

		if (edges[b].reliab < edges[c].reliab){
			return edges[c].reliab;
		}

		for (int p = left + 1; p <= right; ++p) {
			if (edges[p].reliab != edges[left].reliab) {				
				return (edges[p].reliab < edges[left].reliab) ? edges[left].reliab : edges[p].reliab;
			}
			return Double.NaN;
		}
		
		return Double.NaN; //oliford: This wasn't here, which gave a warning. Not sure what to do with it.
	}

	private final int partition(int left, int right, double pivot) {
		while (left <= right)
		{
			while (edges[left].reliab < pivot)
				++left;
			while (edges[right].reliab >= pivot)
				--right;
			
			if (left < right) {
				EdgeInfo tmp = edges[left];
				edges[left] = edges[right];
				edges[right] = tmp;
				++left;
				--right;
			}
		}
		return left;
	}

	private final void quicker_sort(int left, int right) {
		double pivot = find_pivot(edges, left, right);

		if (!Double.isNaN(pivot)) {
			int p = partition(left, right, pivot);
			quicker_sort(left, p - 1);
			quicker_sort(p, right);
		}
	}
	//--------------end quicker_sort algorithm -----------------------------------

	//--------------------start initialse voxels ----------------------------------
	//initialse voxels. See the explanation of the voxel class above.
	//initially every voxel is assumed to belong to a group consisting of only itself
	private void initialiseVOXELs()
	{
		int idx = 0;
		
	    for (int n=0; n < volume_depth; n++) {
			for (int i=0; i < volume_height; i++) {
				for (int j=0; j < volume_width; j++) {
					voxels[idx] = new VoxelInfo();
					//voxels[idx].x = j;
	  				//voxels[idx].y = i;
					//voxels[idx].z = n;
					voxels[idx].increment = 0;
					voxels[idx].number_of_voxels_in_group = 1;		
					voxels[idx].value = wrappedVolume[idx];
					voxels[idx].reliability = 9999999 + idx;
					voxels[idx].input_mask = input_mask[idx];
					voxels[idx].extended_mask = extended_mask[idx];
					voxels[idx].head = voxels[idx];
					voxels[idx].last = voxels[idx];
					voxels[idx].next = null;			
					voxels[idx].new_group = 0;
					voxels[idx].group = -1;
					idx++;
				}
	        }
		}
	}
	//-------------------end initialise voxels -----------
	
	/** gamma function in the paper */
	private final double wrap(double voxel_value) {
		double wrapped_voxel_value;
		
		if (voxel_value > Math.PI)	
			wrapped_voxel_value = voxel_value - 2*Math.PI;
		
		else if (voxel_value < -Math.PI)	
			wrapped_voxel_value = voxel_value + 2*Math.PI;
		
		else 
			wrapped_voxel_value = voxel_value;
		
		return wrapped_voxel_value;
	}

	// pixelL_value is the left pixel,	pixelR_value is the right pixel
	private final int find_wrap(double voxelL_value, double voxelR_value) {
		int wrap_value;
		double difference = voxelL_value - voxelR_value;

		if (difference > Math.PI)	
			wrap_value = -1;
		
		else if (difference < -Math.PI)	
			wrap_value = 1;
		else 
			wrap_value = 0;

		return wrap_value;
	}

	private void extend_mask()	{
		//short hand
		final int vw = volume_width;
		final int vh = volume_height;
		final int vd = volume_depth;
		final int fs = frame_size;
		final int vs = volume_size;
		
		int n, i, j;
		int imp = fs + vw + 1;	//input mask pointer
		int emp = fs + vw + 1;	//extended mask pointer

		//extend the mask for the volume except borders
		for (n=1; n < vd - 1; n++) {
			for (i=1; i < vh - 1; i++) {
				for (j=1; j < vw - 1; j++) {
					
					if( input_mask[imp]				&& input_mask[imp - 1]				&& input_mask[imp + 1] && 
						input_mask[imp + vw]		&& input_mask[imp + vw - 1]			&& input_mask[imp + vw + 1] &&
						input_mask[imp - vw]		&& input_mask[imp - vw - 1]			&& input_mask[imp - vw + 1] &&
						input_mask[imp + fs]		&& input_mask[imp + fs - 1]			&& input_mask[imp + fs + 1] &&
						input_mask[imp + fs - vw]	&& input_mask[imp + fs - vw - 1]	&& input_mask[imp + fs - vw + 1] &&
						input_mask[imp + fs + vw]	&& input_mask[imp + fs + vw - 1]	&& input_mask[imp + fs + vw + 1] &&
						input_mask[imp - fs]		&& input_mask[imp - fs - 1]			&& input_mask[imp - fs + 1] &&
						input_mask[imp - fs - vw]	&& input_mask[imp - fs - vw - 1]	&& input_mask[imp - fs - vw + 1] &&
						input_mask[imp - fs + vw]	&& input_mask[imp - fs + vw - 1]	&& input_mask[imp - fs + vw + 1])
					{		
						extended_mask[emp] = true;
					}
					++emp;
					++imp;
				}
				emp += 2;
				imp += 2;
			}
			emp += 2 * vw;
			imp += 2 * vw;
		}		

		if (x_connectivity) {
			
			//extend the mask to the front side of the phase volume
			imp = fs + vw;	//input mask pointer
			emp = fs + vw;	//extended mask pointer
			for (n=1; n < vd - 1; n++) {
				for (i=1; i < vh - 1; i++) {
					
					if( input_mask[imp]						&& input_mask[imp + vw - 1]			&& input_mask[imp + 1] &&
						input_mask[imp - vw]				&& input_mask[imp + vw]	&& 
						input_mask[imp - fs]				&& input_mask[imp + fs] &&
						input_mask[imp - 1]					&& input_mask[imp + vw + 1] &&
						input_mask[imp - vw + 1]			&& input_mask[imp + 2 * vw - 1] && 
						input_mask[imp - fs - 1]			&& input_mask[imp + fs + vw + 1] && 
						input_mask[imp - fs - vw]			&& input_mask[imp + fs + vw] &&
						input_mask[imp - fs - vw + 1]		&& input_mask[imp + fs + 2 * vw - 1] && 
						input_mask[imp - fs + vw - 1]		&& input_mask[imp + fs + 1]	&& 
						input_mask[imp - fs + 1]			&& input_mask[imp + fs + vw - 1] &&
						input_mask[imp - fs + 2 * vw - 1]	&& input_mask[imp + fs - vw + 1] &&
						input_mask[imp - fs + vw]			&& input_mask[imp + fs - vw] && 
						input_mask[imp - fs + vw + 1]		&& input_mask[imp + fs - 1] )
					{		
						extended_mask[emp] = true;
					}
					emp += vw;
					imp += vw;
				}
				emp += 2 * vw;
				imp += 2 *vw;
			}

			//extend the mask to the rear side of the phase volume
			imp = fs + 2 * vw - 1;	//input mask pointer
			emp = fs + 2 * vw - 1;	//extended mask pointer
			for (n=1; n < vd - 1; n++)
			{
				for (i=1; i < vh - 1; i++)
				{
					if( input_mask[imp]						&& input_mask[imp - vw + 1]			&& input_mask[imp - 1] &&
						input_mask[imp - vw]				&& input_mask[imp + vw]	&& 
						input_mask[imp - fs]				&& input_mask[imp + fs] &&
						input_mask[imp - vw - 1]			&& input_mask[imp + 1] &&
						input_mask[imp + vw - 1]			&& input_mask[imp - 2 * vw + 1] && 
						input_mask[imp - fs - vw - 1]		&& input_mask[imp + fs + 1]	&& 
						input_mask[imp - fs - 2 * vw + 1]	&& input_mask[imp + fs + vw - 1] &&
						input_mask[imp - fs - 1]			&& input_mask[imp + fs - vw + 1] && 
						input_mask[imp - fs - vw + 1]		&& input_mask[imp + fs - 1]	&& 
						input_mask[imp - fs - vw]			&& input_mask[imp + fs + vw] &&
						input_mask[imp - fs + vw - 1]		&& input_mask[imp + fs - 2 * vw + 1] &&
						input_mask[imp - fs + vw]			&& input_mask[imp + fs - vw] && 
						input_mask[imp - fs + 1]			&& input_mask[imp + fs - vw - 1] )					
					{		
						extended_mask[emp] = true;
					}
					emp += vw;
					imp += vw;
				}
				emp += 2 * vw;
				imp += 2 *vw;
			}
		}		

		if (y_connectivity)
		{
			//extend the mask to the left side of the phase volume
			imp = fs + 1;	
			emp = fs + 1;
			for (n=1; n < vd - 1; n++)
			{
				for (j=1; j < vw - 1; j++)
				{
					if( input_mask[imp]						&& input_mask[imp - 1]			&& input_mask[imp + 1] &&
						input_mask[imp + fs - vw]			&& input_mask[imp + vw]	&& 
						input_mask[imp - fs]				&& input_mask[imp + fs] &&
						input_mask[imp + fs - vw - 1]		&& input_mask[imp + vw + 1] &&
						input_mask[imp + fs - vw + 1]		&& input_mask[imp + vw - 1] && 
						input_mask[imp - vw - 1]			&& input_mask[imp + fs + vw + 1] && 
						input_mask[imp - vw]				&& input_mask[imp + fs + vw] &&
						input_mask[imp - vw + 1]			&& input_mask[imp + fs + vw - 1] && 
						input_mask[imp - fs - 1]			&& input_mask[imp + fs + 1]	&& 
						input_mask[imp - fs + 1]			&& input_mask[imp + fs - 1] &&
						input_mask[imp - fs + vw - 1]		&& input_mask[imp + 2 * fs - vw + 1] &&
						input_mask[imp - fs + vw]			&& input_mask[imp + 2 * fs - vw] && 
						input_mask[imp - fs + vw + 1]		&& input_mask[imp + 2 * fs - vw - 1] )
					{		
						extended_mask[emp] = true;
					}
					emp++;
					imp++;
				}
				emp += fs - vw + 2;
				imp += fs - vw + 2;
			}

			//extend the mask to the right side of the phase volume
			imp = 2 * fs - vw + 1;
			emp = 2 * fs - vw + 1;
			for (n=1; n < vd - 1; n++)
			{
				for (j=1; j < vw - 1; j++)
				{
					if( input_mask[imp]						&& input_mask[imp + 1]			&& input_mask[imp - 1] &&
						input_mask[imp - vw]				&& input_mask[imp - fs + vw]	&& 
						input_mask[imp - fs]				&& input_mask[imp + fs] &&
						input_mask[imp - vw - 1]			&& input_mask[imp - fs + vw + 1] &&
						input_mask[imp - vw + 1]			&& input_mask[imp - fs + vw - 1] && 
						input_mask[imp - fs - vw - 1]		&& input_mask[imp + vw + 1]	&& 
						input_mask[imp - fs - vw + 1]		&& input_mask[imp + vw - 1] &&
						input_mask[imp - fs - vw]			&& input_mask[imp + vw] && 
						input_mask[imp - fs - 1]			&& input_mask[imp + fs + 1]	&& 
						input_mask[imp - fs + 1]			&& input_mask[imp + fs - 1] &&
						input_mask[imp - 2 * fs + vw - 1]	&& input_mask[imp + fs - vw + 1] &&
						input_mask[imp - 2 * fs + vw]		&& input_mask[imp + fs - vw] && 
						input_mask[imp - 2 * fs + vw + 1]	&& input_mask[imp + fs - vw - 1] )							
					{		
						extended_mask[emp] = true;
					}
					emp++;
					imp++;
				}
				emp += fs - vw + 2;
				imp += fs - vw + 2;
			}
		}

		if (z_connectivity)
		{
			//extend the mask to the bottom side of the phase volume
			imp = vw + 1;	
			emp = vw + 1;
			for (i=1; i < vh - 1; ++i)
			{
				for (j=1; j < vw - 1; ++j)
				{
					if( input_mask[imp]						&& input_mask[imp - 1]			&& input_mask[imp + 1] &&
						input_mask[imp - vw]				&& input_mask[imp + vw]	&& 
						input_mask[imp + fs]				&& input_mask[imp + vs  -  fs] &&
						input_mask[imp - vw - 1]			&& input_mask[imp + vw + 1] &&
						input_mask[imp - vw + 1]			&& input_mask[imp + vw - 1] && 
						input_mask[imp + vs - fs - vw - 1]	&& input_mask[imp + fs + vw + 1] && 
						input_mask[imp + vs - fs - vw]		&& input_mask[imp + fs + vw] &&
						input_mask[imp + vs - fs - vw + 1]	&& input_mask[imp + fs + vw - 1] && 
						input_mask[imp + vs - fs - 1]		&& input_mask[imp + fs + 1]	&& 
						input_mask[imp + vs - fs + 1]		&& input_mask[imp + fs - 1] &&
						input_mask[imp + vs - fs + vw - 1]	&& input_mask[imp + fs - vw + 1] &&
						input_mask[imp + vs - fs + vw]		&& input_mask[imp + fs - vw] && 
						input_mask[imp + vs - fs + vw + 1]	&& input_mask[imp + fs - vw - 1] )
					{		
						extended_mask[emp] = true;
					}
					emp++;
					imp++;
				}
				emp += 2;
				imp += 2;
			}

			//extend the mask to the top side of the phase volume
			imp = vs - fs + vw + 1;	
			emp = vs - fs + vw + 1;
			for (i=1; i < vh - 1; ++i)
			{
				for (j=1; j < vw - 1; ++j)
				{
					if( input_mask[imp]						&& input_mask[imp + 1]			&& input_mask[imp - 1] &&
						input_mask[imp - vw]				&& input_mask[imp - fs + vw]	&& 
						input_mask[imp - fs]				&& input_mask[imp - vs + fs] &&
						input_mask[imp - vw - 1]			&& input_mask[imp + vw + 1] &&
						input_mask[imp - vw + 1]			&& input_mask[imp + vw - 1] && 
						input_mask[imp - fs - vw - 1]		&& input_mask[imp - vs + fs + vw + 1]	&& 
						input_mask[imp - fs - vw + 1]		&& input_mask[imp - vs + fs + vw - 1] &&
						input_mask[imp - fs - vw]			&& input_mask[imp - vs + fs + vw] && 
						input_mask[imp - fs - 1]			&& input_mask[imp - vs + fs + 1]	&& 
						input_mask[imp - fs + 1]			&& input_mask[imp - vs + fs - 1] &&
						input_mask[imp - fs + vw - 1]		&& input_mask[imp - vs + fs - vw + 1] &&
						input_mask[imp - fs + vw]			&& input_mask[imp - vs + fs - vw] && 
						input_mask[imp - fs + vw + 1]		&& input_mask[imp - vs + fs - vw - 1] )		
					{		
						extended_mask[emp] = true;
					}
					emp++;
					imp++;
				}
				emp += 2;
				imp += 2;
			}
		}
	}

	private void calculate_reliability(){ 
		int voxel_index, wvp;
		double H, V, N, D1, D2, D3, D4, D5, D6, D7, D8, D9, D10;
		
		wvp = frame_size + volume_width + 1;
		voxel_index = frame_size + volume_width + 1;
		for (int n=1; n < volume_depth - 1; n++) {
			for (int i=1; i < volume_height - 1; i++) {
				for (int j=1; j < volume_width - 1; j++) {
					if (voxels[voxel_index].extended_mask) { 
						H  = wrap(wrappedVolume[wvp - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + 1]);
						V  = wrap(wrappedVolume[wvp - volume_width] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + volume_width]);
						N  = wrap(wrappedVolume[wvp - frame_size] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size]);
						D1 = wrap(wrappedVolume[wvp - volume_width - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + volume_width + 1]);
						D2 = wrap(wrappedVolume[wvp - volume_width + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + volume_width - 1]);
						D3 = wrap(wrappedVolume[wvp - frame_size - volume_width - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size + volume_width + 1]);
						D4 = wrap(wrappedVolume[wvp - frame_size - volume_width] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size + volume_width]);
						D5 = wrap(wrappedVolume[wvp - frame_size - volume_width + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size + volume_width - 1]);
						D6 = wrap(wrappedVolume[wvp - frame_size - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size + 1]);
						D7 = wrap(wrappedVolume[wvp - frame_size + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size - 1]);
						D8 = wrap(wrappedVolume[wvp - frame_size + volume_width - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size - volume_width + 1]);
						D9 = wrap(wrappedVolume[wvp - frame_size + volume_width] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size - volume_width]);
						D10 = wrap(wrappedVolume[wvp - frame_size + volume_width + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size - volume_width - 1]);
						voxels[voxel_index].reliability = H*H + V*V + N*N + D1*D1 + D2*D2  + D3*D3 + D4*D4  + D5*D5 + D6*D6  
							+ D7*D7 + D8*D8 + D9*D9 + D10*D10;
					}
					voxel_index++;
					wvp++;
				}
				voxel_index += 2;
				wvp += 2;
			}
			voxel_index += 2 * volume_width;
			wvp += 2 * volume_width;
		}

		if (x_connectivity)
		{
			//calculating reliability for the front side of the phase volume...add volume_width
			wvp =frame_size + volume_width;
			voxel_index = frame_size + volume_width;
			for (int n=1; n < volume_depth - 1; ++n)
			{
				for (int i=1; i < volume_height - 1; ++i)
				{
					if (voxels[voxel_index].extended_mask)
					{ 
						H  = wrap(wrappedVolume[wvp + volume_width - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + 1]);
						V  = wrap(wrappedVolume[wvp - volume_width] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + volume_width]);
						N  = wrap(wrappedVolume[wvp - frame_size] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size]);
						D1 = wrap(wrappedVolume[wvp - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + volume_width + 1]);
						D2 = wrap(wrappedVolume[wvp - volume_width + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + 2 * volume_width - 1]);
						D3 = wrap(wrappedVolume[wvp - frame_size - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size + volume_width + 1]);
						D4 = wrap(wrappedVolume[wvp - frame_size - volume_width] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size + volume_width]);
						D5 = wrap(wrappedVolume[wvp - frame_size - volume_width + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size + 2 * volume_width - 1]);
						D6 = wrap(wrappedVolume[wvp - frame_size + volume_width - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size + 1]);
						D7 = wrap(wrappedVolume[wvp - frame_size + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size + volume_width - 1]);
						D8 = wrap(wrappedVolume[wvp - frame_size + 2 * volume_width - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size - volume_width + 1]);
						D9 = wrap(wrappedVolume[wvp - frame_size + volume_width] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size - volume_width]);
						D10 = wrap(wrappedVolume[wvp - frame_size + volume_width + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size - 1]);
						voxels[voxel_index].reliability = H*H + V*V + N*N + D1*D1 + D2*D2  + D3*D3 + D4*D4  + D5*D5 + D6*D6  
							+ D7*D7 + D8*D8 + D9*D9 + D10*D10;
					}
					voxel_index += volume_width;
					wvp += volume_width;
				}
				voxel_index += 2 * volume_width;
				wvp += 2 * volume_width;
			}
			
			//calculating reliability for the rear side of the phase volume..... subtract volume_width
			wvp =frame_size + 2 * volume_width - 1;
			voxel_index = frame_size + 2 * volume_width - 1;
			for (int n=1; n < volume_depth - 1; ++n)
			{
				for (int i=1; i < volume_height - 1; ++i)
				{
					if (voxels[voxel_index].extended_mask)
					{ 
						H  = wrap(wrappedVolume[wvp - volume_width + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp - 1]);
						V  = wrap(wrappedVolume[wvp - volume_width] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + volume_width]);
						N  = wrap(wrappedVolume[wvp - frame_size] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size]);
						D1 = wrap(wrappedVolume[wvp - volume_width - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + 1]);
						D2 = wrap(wrappedVolume[wvp + volume_width - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp - 2 * volume_width + 1]);
						D3 = wrap(wrappedVolume[wvp - frame_size - volume_width - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size + 1]);
						D4 = wrap(wrappedVolume[wvp - frame_size - 2 * volume_width + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size + volume_width - 1]);
						D5 = wrap(wrappedVolume[wvp - frame_size - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size - volume_width + 1]);
						D6 = wrap(wrappedVolume[wvp - frame_size - volume_width + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size - 1]);
						D7 = wrap(wrappedVolume[wvp - frame_size - volume_width] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size + volume_width]);
						D8 = wrap(wrappedVolume[wvp - frame_size + volume_width - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size - 2 * volume_width + 1]);
						D9 = wrap(wrappedVolume[wvp - frame_size + volume_width] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size - volume_width]);
						D10 = wrap(wrappedVolume[wvp - frame_size + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size - volume_width - 1]);
						voxels[voxel_index].reliability = H*H + V*V + N*N + D1*D1 + D2*D2  + D3*D3 + D4*D4  + D5*D5 + D6*D6  
							+ D7*D7 + D8*D8 + D9*D9 + D10*D10;
					}
					voxel_index += volume_width;
					wvp += volume_width;
				}
				voxel_index += 2 * volume_width;
				wvp += 2 * volume_width;
			}		
		}

		if (y_connectivity)
		{
			//calculating reliability for the left side of the phase volume...add frame_size
			wvp =frame_size + 1;
			voxel_index = frame_size + 1;
			for (int n=1; n < volume_depth - 1; ++n)
			{
				for (int j=1; j < volume_width - 1; ++j)
				{
					if (voxels[voxel_index].extended_mask)
					{ 
						H  = wrap(wrappedVolume[wvp - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + 1]);
						V  = wrap(wrappedVolume[wvp + frame_size - volume_width] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + volume_width]);
						N  = wrap(wrappedVolume[wvp - frame_size] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size]);
						D1  = wrap(wrappedVolume[wvp + frame_size - volume_width - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + volume_width + 1]);
						D2  = wrap(wrappedVolume[wvp + frame_size - volume_width + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + volume_width - 1]);
						D3  = wrap(wrappedVolume[wvp - volume_width - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size + volume_width + 1]);
						D4  = wrap(wrappedVolume[wvp - volume_width] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size + volume_width]);
						D5  = wrap(wrappedVolume[wvp - volume_width + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size + volume_width - 1]);
						D6  = wrap(wrappedVolume[wvp - frame_size - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size + 1]);
						D7  = wrap(wrappedVolume[wvp - frame_size + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size - 1]);
						D8  = wrap(wrappedVolume[wvp - frame_size + volume_width - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + 2 * frame_size - volume_width + 1]);
						D9  = wrap(wrappedVolume[wvp - frame_size + volume_width] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + 2 * frame_size - volume_width]);
						D10  = wrap(wrappedVolume[wvp - frame_size + volume_width + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + 2 * frame_size - volume_width - 1]);
						voxels[voxel_index].reliability = H*H + V*V + N*N + D1*D1 + D2*D2  + D3*D3 + D4*D4  + D5*D5 + D6*D6  
							+ D7*D7 + D8*D8 + D9*D9 + D10*D10;
					}
					voxel_index++;
					wvp++;
				}
				voxel_index += frame_size - volume_width + 2;
				wvp += frame_size - volume_width + 2;
			}	

			//calculating reliability for the right side of the phase volume...subtract frame_size
			wvp =2 * frame_size - volume_width + 1;
			voxel_index = 2 * frame_size - volume_width + 1;
			for (int n=1; n < volume_depth - 1; ++n)
			{
				for (int j=1; j < volume_width - 1; ++j)
				{
					if (voxels[voxel_index].extended_mask)
					{ 
						H  = wrap(wrappedVolume[wvp + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp - 1]);
						V  = wrap(wrappedVolume[wvp - volume_width] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp - frame_size + volume_width]);
						N  = wrap(wrappedVolume[wvp - frame_size] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size]);
						D1  = wrap(wrappedVolume[wvp - volume_width - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp - frame_size + volume_width + 1]);
						D2  = wrap(wrappedVolume[wvp - volume_width + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp - frame_size + volume_width - 1]);
						D3  = wrap(wrappedVolume[wvp - frame_size - volume_width - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + volume_width + 1]);
						D4  = wrap(wrappedVolume[wvp - frame_size - volume_width + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + volume_width -	1]);
						D5  = wrap(wrappedVolume[wvp - frame_size - volume_width] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + volume_width]);
						D6  = wrap(wrappedVolume[wvp - frame_size - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size + 1]);
						D7  = wrap(wrappedVolume[wvp - frame_size + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size - 1]);
						D8  = wrap(wrappedVolume[wvp - 2 * frame_size + volume_width - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size - volume_width + 1]);
						D9  = wrap(wrappedVolume[wvp - 2 * frame_size + volume_width] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size - volume_width]);
						D10  = wrap(wrappedVolume[wvp - 2 * frame_size + volume_width + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size - volume_width - 1]);
						voxels[voxel_index].reliability = H*H + V*V + N*N + D1*D1 + D2*D2  + D3*D3 + D4*D4  + D5*D5 + D6*D6  
							+ D7*D7 + D8*D8 + D9*D9 + D10*D10;
					}
					voxel_index++;
					wvp++;
				}
				voxel_index += frame_size - volume_width + 2;
				wvp += frame_size - volume_width + 2;
			}		
		}

		if (z_connectivity)
		{
			//calculating reliability for the bottom side of the phase volume...add volume_size
			wvp =volume_width + 1;
			voxel_index = volume_width + 1;
			for (int i=1; i < volume_height - 1; ++i)
			{
				for (int j=1; j < volume_width - 1; ++j)
				{
					if (voxels[voxel_index].extended_mask)
					{ 
						H  = wrap(wrappedVolume[wvp - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + 1]);
						V  = wrap(wrappedVolume[wvp - volume_width] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + volume_width]);
						N  = wrap(wrappedVolume[wvp + frame_size] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + volume_size - frame_size]);
						D1  = wrap(wrappedVolume[wvp - volume_width - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + volume_width + 1]);
						D2  = wrap(wrappedVolume[wvp - volume_width + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + volume_width - 1]);
						D3  = wrap(wrappedVolume[wvp + volume_size - frame_size - volume_width - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size + volume_width + 1]);
						D4  = wrap(wrappedVolume[wvp + volume_size - frame_size - volume_width] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size + volume_width]);
						D5  = wrap(wrappedVolume[wvp + volume_size - frame_size - volume_width + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size + volume_width - 1]);
						D6  = wrap(wrappedVolume[wvp + volume_size - frame_size - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size + 1]);
						D7  = wrap(wrappedVolume[wvp + volume_size - frame_size + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size - 1]);
						D8  = wrap(wrappedVolume[wvp + volume_size - frame_size + volume_width - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size - volume_width + 1]);
						D9  = wrap(wrappedVolume[wvp + volume_size - frame_size + volume_width] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size - volume_width]);
						D10 = wrap(wrappedVolume[wvp + volume_size - frame_size + volume_width + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + frame_size - volume_width - 1]);
						voxels[voxel_index].reliability = H*H + V*V + N*N + D1*D1 + D2*D2  + D3*D3 + D4*D4  + D5*D5 + D6*D6  
							+ D7*D7 + D8*D8 + D9*D9 + D10*D10;
					}
					voxel_index++;
					wvp++;
				}
				voxel_index += 2;
				wvp += 2;
			}	

			//calculating reliability for the top side of the phase volume...subtract volume_size
			wvp =volume_size - frame_size + volume_width + 1;
			voxel_index = volume_size - frame_size + volume_width + 1;
			for (int i=1; i < volume_height - 1; ++i) {
				for (int j=1; j < volume_width - 1; ++j) {
					if (voxels[voxel_index].extended_mask) { 
						H  = wrap(wrappedVolume[wvp + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp - 1]);
						V  = wrap(wrappedVolume[wvp - volume_width] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + volume_width]);
						N  = wrap(wrappedVolume[wvp - frame_size] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp - volume_size + frame_size]);
						D1  = wrap(wrappedVolume[wvp - volume_width - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + volume_width + 1]);
						D2  = wrap(wrappedVolume[wvp - volume_width + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp + volume_width - 1]);
						D3  = wrap(wrappedVolume[wvp - frame_size - volume_width - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp - volume_size + frame_size + volume_width + 1]);
						D4  = wrap(wrappedVolume[wvp - frame_size - volume_width + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp - volume_size + frame_size + volume_width - 1]);
						D5  = wrap(wrappedVolume[wvp - frame_size - volume_width] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp - volume_size + frame_size + volume_width]);	
						D6  = wrap(wrappedVolume[wvp - frame_size - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp - volume_size + frame_size + 1]);
						D7  = wrap(wrappedVolume[wvp - frame_size + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp - volume_size + frame_size - 1]);
						D8  = wrap(wrappedVolume[wvp - frame_size + volume_width - 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp - volume_size + frame_size - volume_width + 1]);
						D9  = wrap(wrappedVolume[wvp - frame_size + volume_width] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp - volume_size + frame_size - volume_width]);
						D10 = wrap(wrappedVolume[wvp - frame_size + volume_width + 1] - wrappedVolume[wvp]) - wrap(wrappedVolume[wvp] - wrappedVolume[wvp - volume_size + frame_size - volume_width - 1]);
						voxels[voxel_index].reliability = H*H + V*V + N*N + D1*D1 + D2*D2  + D3*D3 + D4*D4  + D5*D5 + D6*D6  
							+ D7*D7 + D8*D8 + D9*D9 + D10*D10;
					}
					voxel_index++;
					wvp++;
				}
				voxel_index += 2;
				wvp += 2;
			}	
		}
	}

	//calculate the reliability of the horizontal edges of the volume
	//it is calculated by adding the reliability of voxel and the relibility of 
	//its right neighbour
	//edge is calculated between a voxel and its next neighbour
	private void horizontalEdgeInfos() {
		int edge_index = 0;
		int voxel_index = 0;
		
		for (int n=0; n < volume_depth; n++) {
			for (int i = 0; i < volume_height; i++) {
				for (int j = 0; j < volume_width - 1; j++) {
					//System.out.println("H\t" + edge_index + "\t" + n +"\t" + i + "\t" + voxel_index);
					if (voxels[voxel_index].input_mask && voxels[voxel_index+1].input_mask) {
						edges[edge_index] = new EdgeInfo();
						edges[edge_index].pointer_1 = voxels[voxel_index];
						edges[edge_index].pointer_2 = voxels[voxel_index + 1];
						edges[edge_index].reliab = voxels[voxel_index].reliability + voxels[voxel_index + 1].reliability;
						edges[edge_index].increment = find_wrap(voxels[voxel_index].value, voxels[voxel_index + 1].value);
						edge_index++;
						No_of_edges++;
					}
					voxel_index++;
				}
				voxel_index++;
			}
		}
		if (x_connectivity)
		{
			voxel_index = volume_width - 1;
			for (int n=0; n < volume_depth; n++)
			{
				for (int i = 0; i < volume_height; i++)
				{
					//System.out.println("Hx\t" + edge_index + "\t" + n +"\t" + i + "\t" + voxel_index);
					if (voxels[voxel_index].input_mask && voxels[voxel_index - volume_width + 1].input_mask )
					{
						edges[edge_index] = new EdgeInfo();
						edges[edge_index].pointer_1 = voxels[voxel_index];
						edges[edge_index].pointer_2 = voxels[voxel_index - volume_width + 1];
						edges[edge_index].reliab = voxels[voxel_index].reliability + voxels[voxel_index - volume_width + 1].reliability;
						edges[edge_index].increment = find_wrap(voxels[voxel_index].value, voxels[voxel_index - volume_width + 1].value);
						edge_index++;
						No_of_edges++;
					}
					voxel_index += volume_width;
				}
			}
		}
	}

	private void verticalEdgeInfos() {
		int voxel_index = 0;
		int edge_index = No_of_edges; 

		for (int n=0; n < volume_depth; n++) {
			for (int i=0; i<volume_height - 1; i++) {
				for (int j=0; j < volume_width; j++)  {
					//System.out.println("V\t" + edge_index + "\t" + n +"\t" + i + "\t" + voxel_index);
					if (voxels[voxel_index].input_mask && voxels[voxel_index + volume_width].input_mask )
					{
						edges[edge_index] = new EdgeInfo();
						edges[edge_index].pointer_1 = voxels[voxel_index];
						edges[edge_index].pointer_2 = voxels[voxel_index + volume_width];
						edges[edge_index].reliab = voxels[voxel_index].reliability + voxels[voxel_index + volume_width].reliability;
						edges[edge_index].increment = find_wrap(voxels[voxel_index].value, voxels[voxel_index + volume_width].value);
						edge_index++;
						No_of_edges++;
					}
					voxel_index++;
				}
			}
			voxel_index += volume_width;
		} 

		int frame_size = volume_width * volume_height;
		int next_voxel = frame_size - volume_width;
		if (y_connectivity)
		{
			voxel_index = frame_size - volume_width;
			for (int n=0; n < volume_depth; n++)
			{
				for (int i = 0; i < volume_width; i++)
				{
					//System.out.println("Vy\t" + edge_index + "\t" + n +"\t" + i + "\t" + voxel_index);
					if (voxels[voxel_index].input_mask && voxels[voxel_index - next_voxel].input_mask )
					{
						edges[edge_index] = new EdgeInfo();
						edges[edge_index].pointer_1 = voxels[voxel_index];
						edges[edge_index].pointer_2 = voxels[voxel_index - next_voxel];
						edges[edge_index].reliab = voxels[voxel_index].reliability + voxels[voxel_index - next_voxel].reliability;
						edges[edge_index].increment = find_wrap(voxels[voxel_index].value, voxels[voxel_index - next_voxel].value);
						edge_index++;
						No_of_edges++;
					}
					voxel_index++;
				}
				voxel_index += next_voxel + 1;
			}
		}
	}

	private void normalEdgeInfos() {
		int frame_size = volume_width * volume_height;
		int volume_size = volume_width * volume_height * volume_depth;
		int voxel_index = 0;
		int edge_index = No_of_edges; 

		for (int n=0; n < volume_depth - 1; n++) {
			for (int i=0; i<volume_height; i++) {
				for (int j=0; j < volume_width; j++) {
					//System.out.println("N\t" + edge_index + "\t" + n +"\t" + i + "\t" + voxel_index);
					if (voxels[voxel_index].input_mask && voxels[voxel_index + frame_size].input_mask )
					{
						edges[edge_index] = new EdgeInfo();
						edges[edge_index].pointer_1 = voxels[voxel_index];
						edges[edge_index].pointer_2 = voxels[voxel_index + frame_size];
						edges[edge_index].reliab = voxels[voxel_index].reliability + voxels[voxel_index + frame_size].reliability;
						edges[edge_index].increment = find_wrap(voxels[voxel_index].value, voxels[voxel_index + frame_size].value);
						edge_index++;
						No_of_edges++;
					}
					voxel_index++;
				}
			}
		}
		
		int next_voxel = volume_size - frame_size;
		if (z_connectivity)
		{
			voxel_index = next_voxel;
			for (int i=0; i < volume_height; i++)
			{
				for (int j = 0; j < volume_width; j++)
				{
					//System.out.println("Nz\t" + edge_index + "\t" + i +"\t" + j + "\t" + voxel_index);
					if (voxels[voxel_index].input_mask && voxels[voxel_index - next_voxel].input_mask)
					{
						edges[edge_index] = new EdgeInfo();
						edges[edge_index].pointer_1 = voxels[voxel_index];
						edges[edge_index].pointer_2 = voxels[voxel_index - next_voxel];
						edges[edge_index].reliab = voxels[voxel_index].reliability + voxels[voxel_index - next_voxel].reliability;
						edges[edge_index].increment = find_wrap(voxels[voxel_index].value, voxels[voxel_index - next_voxel].value);
						edge_index++;
						No_of_edges++;
					}
					voxel_index++;
				}
			}
		}
	}

	//gather the voxels of the volume into groups 
	private void gatherVOXELs() {		
		VoxelInfo VOXEL1;   
		VoxelInfo VOXEL2;
		VoxelInfo group1;
		VoxelInfo group2;
		int edge_index = 0;
		int incremento;

		for (int k = 0; k < No_of_edges; k++) {
			
			VOXEL1 = edges[edge_index].pointer_1;
			VOXEL2 = edges[edge_index].pointer_2;

			//VoxelInfo 1 and VoxelInfo 2 belong to different groups
			//initially each voxel is in a group by itself and one voxel can construct a group
			//no else or else if to this if
			if (VOXEL2.head != VOXEL1.head)
			{
				//VoxelInfo 2 is alone in its group
				//merge this voxel with VoxelInfo 1 group and find the number of 2 pi to add 
				//to or subtract to unwrap it
				if ((VOXEL2.next == null) && (VOXEL2.head == VOXEL2))
				{
					VOXEL1.head.last.next = VOXEL2;
					VOXEL1.head.last = VOXEL2;
					(VOXEL1.head.number_of_voxels_in_group)++;
					VOXEL2.head=VOXEL1.head;
					VOXEL2.increment = VOXEL1.increment - edges[edge_index].increment;
				}

				//VoxelInfo 1 is alone in its group
				//merge this voxel with VoxelInfo 2 group and find the number of 2 pi to add 
				//to or subtract to unwrap it
				else if ((VOXEL1.next == null) && (VOXEL1.head == VOXEL1))
				{
					VOXEL2.head.last.next = VOXEL1;
					VOXEL2.head.last = VOXEL1;
					(VOXEL2.head.number_of_voxels_in_group)++;
					VOXEL1.head = VOXEL2.head;
					VOXEL1.increment = VOXEL2.increment + edges[edge_index].increment;
				} 

				//VoxelInfo 1 and VoxelInfo 2 both have groups
				else
	            {
					group1 = VOXEL1.head;
	                group2 = VOXEL2.head;
					//if the no. of voxels in VoxelInfo 1 group is larger than the no. of voxels
					//in VoxelInfo 2 group.   Merge VoxelInfo 2 group to VoxelInfo 1 group
					//and find the number of wraps between VoxelInfo 2 group and VoxelInfo 1 group
					//to unwrap VoxelInfo 2 group with respect to VoxelInfo 1 group.
					//the no. of wraps will be added to VoxelInfo 2 grop in the future
					if (group1.number_of_voxels_in_group > group2.number_of_voxels_in_group)
					{
						//merge VoxelInfo 2 with VoxelInfo 1 group
						group1.last.next = group2;
						group1.last = group2.last;
						group1.number_of_voxels_in_group = group1.number_of_voxels_in_group + group2.number_of_voxels_in_group;
						incremento = VOXEL1.increment - edges[edge_index].increment - VOXEL2.increment;
						//merge the other voxels in VoxelInfo 2 group to VoxelInfo 1 group
						while (group2 != null)
						{
							group2.head = group1;
							group2.increment += incremento;
							group2 = group2.next;
						}
					} 

					//if the no. of voxels in VoxelInfo 2 group is larger than the no. of voxels
					//in VoxelInfo 1 group.   Merge VoxelInfo 1 group to VoxelInfo 2 group
					//and find the number of wraps between VoxelInfo 2 group and VoxelInfo 1 group
					//to unwrap VoxelInfo 1 group with respect to VoxelInfo 2 group.
					//the no. of wraps will be added to VoxelInfo 1 grop in the future
					else
	                {
						//merge VoxelInfo 1 with VoxelInfo 2 group
						group2.last.next = group1;
						group2.last = group1.last;
						group2.number_of_voxels_in_group = group2.number_of_voxels_in_group + group1.number_of_voxels_in_group;
						incremento = VOXEL2.increment + edges[edge_index].increment - VOXEL1.increment;
						//merge the other voxels in VoxelInfo 2 group to VoxelInfo 1 group
						while (group1 != null)
						{
							group1.head = group2;
							group1.increment += incremento;
							group1 = group1.next;
						} // while

	                } // else
	            } //else
	        } //if
	        edge_index++;
		}
	}

	//unwrap the volume 
	private void unwrapVolume() {
		int i;
		int volume_size = volume_width * volume_height * volume_depth;
		int voxel_index = 0;

		for (i = 0; i < volume_size; i++) {
			voxels[voxel_index].value += 2*Math.PI * (float)(voxels[voxel_index].increment);
			voxel_index++;
	    }
	}

	//set the masked voxels (mask = 0) to the minimum of the unwrapper phase
	private void maskVolume() {
		int volume_width_plus_one  = volume_width + 1;
		int volume_height_plus_one  = volume_height + 1;
		int volume_width_minus_one = volume_width - 1;
		int volume_height_minus_one = volume_height - 1;

		int voxel_index = 0;
		int imp = 0;	//input mask pointer
		double min = Double.POSITIVE_INFINITY;
		int i, j;
		int volume_size = volume_width * volume_height * volume_depth;

		//find the minimum of the unwrapped phase
		for (i = 0; i < volume_size; i++)
		{
			if ((voxels[voxel_index].value < min) && input_mask[imp]) 
				min = voxels[voxel_index].value;

				voxel_index++;
			imp++;
		}

		voxel_index = 0;
		imp = 0;

		//set the masked voxels to minimum
		for (i = 0; i < volume_size; i++){
			if(!input_mask[imp]) {
				voxels[voxel_index].value = min;
			}
			voxel_index++;
			imp++;
		}
	}

	//the input to this unwrapper is an array that contains the wrapped phase map. 
	//copy the volume on the buffer passed to this unwrapper to over-write the unwrapped 
	//phase map on the buffer of the wrapped phase map.
	private void returnVolume()
	{
		int i;
		int volume_size = volume_width * volume_height * volume_depth;
	    int uwv = 0;
	    int voxel_index = 0;

	    for (i=0; i < volume_size; i++) {
	        unwrappedVolume[uwv] = voxels[voxel_index].value;
	        voxel_index++;
	        uwv++;
		}
	}

	private int volume_width;
	private int volume_height;
	private int volume_depth;
	private int volume_size;
	private int frame_size;
	private double wrappedVolume[];
	private boolean input_mask[];
	private boolean extended_mask[];
	private double unwrappedVolume[];
	private VoxelInfo voxels[];
	private EdgeInfo edges[];
	
	public static double[] unwrap(double wrappedVolumeFlat[], boolean input_mask[], 
			int width, int height, int depth, 
			boolean x_connectivity, boolean y_connectivity, boolean z_connectivity){
		
		return (new PhaseUnwrap3D(wrappedVolumeFlat, input_mask, 
				width, height, depth, 
				x_connectivity, y_connectivity, z_connectivity)).getFlatUnwrappedVolume();
	}

	public PhaseUnwrap3D(double wrappedVolume[], boolean input_mask[], 
			int width, int height, int depth, 
			boolean x_connectivity, boolean y_connectivity, boolean z_connectivity) {

		this.x_connectivity = x_connectivity;
		this.y_connectivity = y_connectivity;
		this.z_connectivity = z_connectivity;
		this.wrappedVolume = wrappedVolume;
		this.input_mask = input_mask;
		this.volume_width = width;
		this.volume_height = height;
		this.volume_depth = depth;		
		frame_size  = volume_width * volume_height;		
		volume_size = frame_size * volume_depth;
		int No_of_Edges_initially = 3 * volume_width * volume_height * volume_depth;

		if(input_mask == null){
			this.input_mask = new boolean[volume_size];
			for(int i=0; i < volume_size; i++)
				this.input_mask[i] = true;
		}
		
		
		unwrappedVolume = new double[volume_size];
		
		extended_mask = new boolean[volume_size];
		extend_mask();

		voxels = new VoxelInfo[volume_size];
		edges = new EdgeInfo[No_of_Edges_initially];
		
		initialiseVOXELs();
		calculate_reliability();
		horizontalEdgeInfos();
		verticalEdgeInfos();
		normalEdgeInfos();

		//sort the EdgeInfos depending on their reiability. The VOXELs with higher relibility (small value) first
		quicker_sort(0, No_of_edges - 1);

		//gather VOXELs into groups
		gatherVOXELs();

		unwrapVolume();

		maskVolume();

		//copy the volume from VoxelInfo structure to the unwrapped phase array passed to this function
		returnVolume();
	}
	
	private double[] getFlatUnwrappedVolume() { return unwrappedVolume;  }

	public static void main(String[] args) {
		//csrc/phaseUnwrap3D/target-100x150x60.bin
		//double wrappedVolumeFlat[] = BinaryMatrixFile.mustLoad("csrc/phaseUnwrap3D/wrapped-100x150x60.bin", true)[0];
		double wrappedVolumeFlat[] = BinaryMatrixFile.mustLoad("/tmp/df.bin", true)[0];
		
		/*double maskData[] = BinaryMatrixFile.mustLoad("csrc/phaseUnwrap3D/mask3D-flatBMF.bin", true)[0];
		 * boolean input_mask =  new boolean[volume_size];
			for(int i=0; i < volume_size; i++)
				input_mask[i] = (maskData[i] != 0);		
		 */
		
		//int w = 60, h = 150, d = 100;
		int w = 320, h = 240, d = 242;
		
		double unwrappedVolumeFlat[] = PhaseUnwrap3D.unwrap(wrappedVolumeFlat, null, w, h, d, false, false, false);
		
		//BinaryMatrixFile.mustWrite("csrc/phaseUnwrap3D/unwrapped-flatBMF-java.bin", unwrappedVolumeFlat);
		BinaryMatrixFile.mustWrite("/tmp/u.bin", unwrappedVolumeFlat);
	}
}