package edu.stanford.rsl.apps.gui.blobdetection;

/**
 * 
 * This code was taken from the BoneJ GIT Hub.
 * <a href="https://github.com/mdoube/BoneJ/"></a>
 * The code was modified such that it no longer provides an ImageJ
 * plugin. It is now to be used inside CONRAD for detecting 3D
 * blobs.
 *
 * Please find the original copyright notice below.
 * 
 * @author Martin Berger
 * 
 *
 * 
 * ParticleCounter Copyright 2009 2010 2011 Michael Doube
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program. If not, see <http://www.gnu.org/licenses/>.
 */

import java.awt.AWTEvent;
import java.awt.Checkbox;
import java.awt.Choice;
import java.awt.TextField;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.HashSet;
import java.util.Iterator;
import java.util.List;
import java.util.ListIterator;
import java.util.Vector;
import java.util.concurrent.atomic.AtomicInteger;


import org.doube.geometry.FitEllipsoid;
import org.doube.jama.EigenvalueDecomposition;
import org.doube.jama.Matrix;
import org.doube.util.DialogModifier;
import org.doube.util.Multithreader;

import edu.stanford.rsl.conrad.geometry.shapes.simple.Point3D;

import ij.IJ;
import ij.ImagePlus;
import ij.ImageStack;
import ij.gui.GenericDialog;
import ij.measure.Calibration;
import ij.process.ImageProcessor;

/**
 * <p>
 * This class implements multithreaded and linear O(n) 3D particle
 * identification and shape analysis. Surface meshing and 3D visualisation are
 * provided by Bene Schmid's ImageJ 3D Viewer.
 * </p>
 * <p>
 * This plugin is based on Object_Counter3D by Fabrice P Cordelires and Jonathan
 * Jackson, but with significant speed increases through reduction of recursion
 * and multi-threading. Thanks to Robert Barbour for the suggestion to 'chunk'
 * the stack. Chunking works as follows:
 * </p>
 * <ol>
 * <li>Perform initial labelling on the whole stack in a single thread</li>
 * <li>for <i>n</i> discrete, contiguous chunks within the labelling array,
 * connectStructures()
 * <ol type="a">
 * <li>connectStructures() can run in a separate thread for each chunk</li>
 * <li>chunks are approximately equal-sized sets of slices</li>
 * </ol>
 * <li>stitchChunks() for the pixels on the first slice of each chunk, except
 * for the first chunk, restricting replaceLabels() to the current and all
 * previous chunks.
 * <ol type="a">
 * <li>stitchChunks() iterates through the slice being stitched in a single
 * thread</li>
 * </ol>
 * </li>
 *
 * </ol>
 * <p>
 * The performance improvement should be in the region of a factor of <i>n</i>
 * if run linearly, and if multithreaded over <i>c</i> processors, speed
 * increase should be in the region of <i>n</i> * <i>c</i>, minus overhead.
 * </p>
 *
 * @author Michael Doube
 * @author Jonathan Jackson
 * @author Fabrice Cordelires
 * @author Michal Klosowski
 * @see <p>
 * <a href="http://rsbweb.nih.gov/ij/plugins/track/objects.html">3D Object
 * Counter</a>
 * </p>
 *
 */
public class ConnectedComponent3D {

	/** Foreground value */
	public final static int FORE = -1;

	/** Background value */
	public final static int BACK = 0;

	/** Particle joining method */
	public final static int MULTI = 0, LINEAR = 1, MAPPED = 2;

	/** Surface colour style */
	//private final static int GRADIENT = 0, SPLIT = 1;

	private String sPhase = "";

	private String chunkString = "";

	private int labelMethod = MAPPED;
	

	public Object[][] getEllipsoids(ArrayList<List<Point3D>> surfacePoints) {
		Object[][] ellipsoids = new Object[surfacePoints.size()][];
		int p = 0;
		Iterator<List<Point3D>> partIter = surfacePoints.iterator();
		while (partIter.hasNext()) {
			List<Point3D> points = partIter.next();
			if (points == null) {
				p++;
				continue;
			}
			Iterator<Point3D> pointIter = points.iterator();
			double[][] coOrdinates = new double[points.size()][3];
			int i = 0;
			while (pointIter.hasNext()) {
				Point3D point = pointIter.next();
				coOrdinates[i][0] = point.getX();
				coOrdinates[i][1] = point.getY();
				coOrdinates[i][2] = point.getZ();
				i++;
			}
			try {
				ellipsoids[p] = FitEllipsoid.yuryPetrov(coOrdinates);
			} catch (RuntimeException re) {
				IJ.log("Could not fit ellipsoid to surface " + p);
				ellipsoids[p] = null;
			}
			p++;
		}
		return ellipsoids;
	}

	/**
	 * Get the mean and standard deviation of pixel values above a minimum value
	 * for each particle in a particle label work array
	 *
	 * @param imp
	 * Input image containing pixel values
	 * @param particleLabels
	 * workArray containing particle labels
	 * @param particleSizes
	 * array of particle sizes as pixel counts
	 * @param threshold
	 * restrict calculation to values > i
	 * @return array containing mean, std dev and max pixel values for each
	 * particle
	 */
	public double[][] getMeanStdDev(ImagePlus imp, int[][] particleLabels,
			long[] particleSizes, final int threshold) {
		final int nParticles = particleSizes.length;
		final int d = imp.getImageStackSize();
		final int wh = imp.getWidth() * imp.getHeight();
		ImageStack stack = imp.getImageStack();
		double[] sums = new double[nParticles];
		for (int z = 0; z < d; z++) {
			float[] pixels = (float[]) stack.getPixels(z + 1);
			int[] labelPixels = particleLabels[z];
			for (int i = 0; i < wh; i++) {
				final double value = pixels[i];
				if (value > threshold) {
					sums[labelPixels[i]] += value;
				}
			}
		}
		double[][] meanStdDev = new double[nParticles][3];
		for (int p = 1; p < nParticles; p++) {
			meanStdDev[p][0] = sums[p] / particleSizes[p];
		}

		double[] sumSquares = new double[nParticles];
		for (int z = 0; z < d; z++) {
			float[] pixels = (float[]) stack.getPixels(z + 1);
			int[] labelPixels = particleLabels[z];
			for (int i = 0; i < wh; i++) {
				final double value = pixels[i];
				if (value > threshold) {
					final int p = labelPixels[i];
					final double residual = value - meanStdDev[p][0];
					sumSquares[p] += residual * residual;
					meanStdDev[p][2] = Math.max(meanStdDev[p][2], value);
				}
			}
		}
		for (int p = 1; p < nParticles; p++) {
			meanStdDev[p][1] = Math.sqrt(sumSquares[p] / particleSizes[p]);
		}
		return meanStdDev;
	}

	public int getNCavities(ImagePlus imp) {
		Object[] result = getParticles(imp, 4, BACK);
		long[] particleSizes = (long[]) result[2];
		final int nParticles = particleSizes.length;
		final int nCavities = nParticles - 2; // 1 particle is the background
		return nCavities;
	}

	/**
	 * Get the minimum and maximum x, y and z coordinates of each particle
	 *
	 * @param imp
	 * ImagePlus (used for stack size)
	 * @param particleLabels
	 * work array containing labelled particles
	 * @param nParticles
	 * number of particles in the stack
	 * @return int[][] containing x, y and z minima and maxima.
	 */
	public int[][] getParticleLimits(ImagePlus imp, int[][] particleLabels,
			int nParticles) {
		final int w = imp.getWidth();
		final int h = imp.getHeight();
		final int d = imp.getImageStackSize();
		int[][] limits = new int[nParticles][6];
		for (int i = 0; i < nParticles; i++) {
			limits[i][0] = Integer.MAX_VALUE; // x min
			limits[i][1] = 0; // x max
			limits[i][2] = Integer.MAX_VALUE; // y min
			limits[i][3] = 0; // y max
			limits[i][4] = Integer.MAX_VALUE; // z min
			limits[i][5] = 0; // z max
		}
		for (int z = 0; z < d; z++) {
			for (int y = 0; y < h; y++) {
				final int index = y * w;
				for (int x = 0; x < w; x++) {
					final int i = particleLabels[z][index + x];
					limits[i][0] = Math.min(limits[i][0], x);
					limits[i][1] = Math.max(limits[i][1], x);
					limits[i][2] = Math.min(limits[i][2], y);
					limits[i][3] = Math.max(limits[i][3], y);
					limits[i][4] = Math.min(limits[i][4], z);
					limits[i][5] = Math.max(limits[i][5], z);
				}
			}
		}
		return limits;
	}

	public EigenvalueDecomposition[] getEigens(ImagePlus imp,
			int[][] particleLabels, double[][] centroids) {
		Calibration cal = imp.getCalibration();
		final double vW = cal.pixelWidth;
		final double vH = cal.pixelHeight;
		final double vD = cal.pixelDepth;
		final double voxVhVd = (vH * vH + vD * vD) / 12;
		final double voxVwVd = (vW * vW + vD * vD) / 12;
		final double voxVhVw = (vH * vH + vW * vW) / 12;
		final int w = imp.getWidth();
		final int h = imp.getHeight();
		final int d = imp.getImageStackSize();
		final int nParticles = centroids.length;
		EigenvalueDecomposition[] eigens = new EigenvalueDecomposition[nParticles];
		double[][] momentTensors = new double[nParticles][6];
		for (int z = 0; z < d; z++) {
			IJ.showStatus("Calculating particle moments...");
			IJ.showProgress(z, d);
			final double zVd = z * vD;
			for (int y = 0; y < h; y++) {
				final double yVh = y * vH;
				final int index = y * w;
				for (int x = 0; x < w; x++) {
					final int p = particleLabels[z][index + x];
					if (p > 0) {
						final double xVw = x * vW;
						final double dx = xVw - centroids[p][0];
						final double dy = yVh - centroids[p][1];
						final double dz = zVd - centroids[p][2];
						momentTensors[p][0] += dy * dy + dz * dz + voxVhVd; // Ixx
						momentTensors[p][1] += dx * dx + dz * dz + voxVwVd; // Iyy
						momentTensors[p][2] += dy * dy + dx * dx + voxVhVw; // Izz
						momentTensors[p][3] += dx * dy; // Ixy
						momentTensors[p][4] += dx * dz; // Ixz
						momentTensors[p][5] += dy * dz; // Iyz
					}
				}
			}
			for (int p = 1; p < nParticles; p++) {
				double[][] inertiaTensor = new double[3][3];
				inertiaTensor[0][0] = momentTensors[p][0];
				inertiaTensor[1][1] = momentTensors[p][1];
				inertiaTensor[2][2] = momentTensors[p][2];
				inertiaTensor[0][1] = -momentTensors[p][3];
				inertiaTensor[0][2] = -momentTensors[p][4];
				inertiaTensor[1][0] = -momentTensors[p][3];
				inertiaTensor[1][2] = -momentTensors[p][5];
				inertiaTensor[2][0] = -momentTensors[p][4];
				inertiaTensor[2][1] = -momentTensors[p][5];
				Matrix inertiaTensorMatrix = new Matrix(inertiaTensor);
				EigenvalueDecomposition E = new EigenvalueDecomposition(
						inertiaTensorMatrix);
				eigens[p] = E;
			}
		}
		return eigens;
	}

	/**
	 * Get the maximum distances from the centroid in x, y, and z axes, and
	 * transformed x, y and z axes
	 *
	 * @param imp
	 * @param particleLabels
	 * @param centroids
	 * @param E
	 * @return array containing two nPoints * 3 arrays with max and max
	 * transformed distances respectively
	 *
	 */
	public Object[] getMaxDistances(ImagePlus imp, int[][] particleLabels,
			double[][] centroids, EigenvalueDecomposition[] E) {
		Calibration cal = imp.getCalibration();
		final double vW = cal.pixelWidth;
		final double vH = cal.pixelHeight;
		final double vD = cal.pixelDepth;
		final int w = imp.getWidth();
		final int h = imp.getHeight();
		final int d = imp.getImageStackSize();
		final int nParticles = centroids.length;
		double[][] maxD = new double[nParticles][3];
		double[][] maxDt = new double[nParticles][3];
		for (int z = 0; z < d; z++) {
			for (int y = 0; y < h; y++) {
				final int index = y * w;
				for (int x = 0; x < w; x++) {
					final int p = particleLabels[z][index + x];
					if (p > 0) {
						final double dX = x * vW - centroids[p][0];
						final double dY = y * vH - centroids[p][1];
						final double dZ = z * vD - centroids[p][2];
						maxD[p][0] = Math.max(maxD[p][0], Math.abs(dX));
						maxD[p][1] = Math.max(maxD[p][1], Math.abs(dY));
						maxD[p][2] = Math.max(maxD[p][2], Math.abs(dZ));
						final double[][] eV = E[p].getV().getArray();
						final double dXt = dX * eV[0][0] + dY * eV[0][1] + dZ
								* eV[0][2];
						final double dYt = dX * eV[1][0] + dY * eV[1][1] + dZ
								* eV[1][2];
						final double dZt = dX * eV[2][0] + dY * eV[2][1] + dZ
								* eV[2][2];
						maxDt[p][0] = Math.max(maxDt[p][0], Math.abs(dXt));
						maxDt[p][1] = Math.max(maxDt[p][1], Math.abs(dYt));
						maxDt[p][2] = Math.max(maxDt[p][2], Math.abs(dZt));
					}
				}
			}
		}
		for (int p = 0; p < nParticles; p++) {
			Arrays.sort(maxDt[p]);
			double[] temp = new double[3];
			for (int i = 0; i < 3; i++) {
				temp[i] = maxDt[p][2 - i];
			}
			maxDt[p] = temp.clone();
		}
		final Object[] maxDistances = { maxD, maxDt };
		return maxDistances;
	}

	
	/**
	 * Get the Feret diameter of a surface. Uses an inefficient brute-force
	 * algorithm.
	 *
	 * @param particleSurfaces
	 * @return the array of ferets
	 */
	public double[] getFerets(ArrayList<List<Point3D>> particleSurfaces) {
		int nParticles = particleSurfaces.size();
		double[] ferets = new double[nParticles];
		ListIterator<List<Point3D>> it = particleSurfaces.listIterator();
		int i = 0;
		Point3D a;
		Point3D b;
		List<Point3D> surface;
		ListIterator<Point3D> ita;
		ListIterator<Point3D> itb;
		while (it.hasNext()) {
			IJ.showStatus("Finding Feret diameter...");
			IJ.showProgress(it.nextIndex(), nParticles);
			surface = it.next();
			if (surface == null) {
				ferets[i] = Double.NaN;
				i++;
				continue;
			}
			ita = surface.listIterator();
			while (ita.hasNext()) {
				a = ita.next();
				itb = surface.listIterator(ita.nextIndex());
				while (itb.hasNext()) {
					b = itb.next();
					ferets[i] = Math.max(ferets[i], a.euclideanDistance(b));
				}
			}
			i++;
		}
		return ferets;
	}

	/**
	 * create a binary ImagePlus containing a single particle and which 'just
	 * fits' the particle
	 *
	 * @param p
	 * The particle ID to get
	 * @param imp
	 * original image, used for calibration
	 * @param particleLabels
	 * work array of particle labels
	 * @param limits
	 * x,y and z limits of each particle
	 * @param padding
	 * amount of empty space to pad around each particle
	 * @return the imageplus
	 */
	public static ImagePlus getBinaryParticle(int p, ImagePlus imp,
			int[][] particleLabels, int[][] limits, int padding) {

		final int w = imp.getWidth();
		final int h = imp.getHeight();
		final int d = imp.getImageStackSize();
		final int xMin = Math.max(0, limits[p][0] - padding);
		final int xMax = Math.min(w - 1, limits[p][1] + padding);
		final int yMin = Math.max(0, limits[p][2] - padding);
		final int yMax = Math.min(h - 1, limits[p][3] + padding);
		final int zMin = Math.max(0, limits[p][4] - padding);
		final int zMax = Math.min(d - 1, limits[p][5] + padding);
		final int stackWidth = xMax - xMin + 1;
		final int stackHeight = yMax - yMin + 1;
		final int stackSize = stackWidth * stackHeight;
		ImageStack stack = new ImageStack(stackWidth, stackHeight);
		for (int z = zMin; z <= zMax; z++) {
			byte[] slice = new byte[stackSize];
			int i = 0;
			for (int y = yMin; y <= yMax; y++) {
				final int sourceIndex = y * w;
				for (int x = xMin; x <= xMax; x++) {
					if (particleLabels[z][sourceIndex + x] == p) {
						slice[i] = (byte) (255 & 0xFF);
					}
					i++;
				}
			}
			stack.addSlice(imp.getStack().getSliceLabel(z + 1), slice);
		}
		ImagePlus binaryImp = new ImagePlus("Particle_" + p, stack);
		Calibration cal = imp.getCalibration();
		binaryImp.setCalibration(cal);
		return binaryImp;
	}

	/**
	 * Create an image showing some particle measurement
	 *
	 * @param imp
	 * @param particleLabels
	 * @param values
	 * list of values whose array indices correspond to
	 * particlelabels
	 * @param title
	 * tag stating what we are displaying
	 * @return ImagePlus with particle labels substituted with some value
	 */
	public ImagePlus displayParticleValues(ImagePlus imp,
			int[][] particleLabels, double[] values, String title) {
		final int w = imp.getWidth();
		final int h = imp.getHeight();
		final int d = imp.getImageStackSize();
		final int wh = w * h;
		float[][] pL = new float[d][wh];
		values[0] = 0; // don't colour the background
		ImageStack stack = new ImageStack(w, h);
		for (int z = 0; z < d; z++) {
			for (int i = 0; i < wh; i++) {
				final int p = particleLabels[z][i];
				pL[z][i] = (float) values[p];
			}
			stack.addSlice(imp.getImageStack().getSliceLabel(z + 1), pL[z]);
		}
		final int nValues = values.length;
		double max = 0;
		for (int i = 0; i < nValues; i++) {
			max = Math.max(max, values[i]);
		}
		ImagePlus impOut = new ImagePlus(imp.getShortTitle() + "_" + title,
				stack);
		impOut.setCalibration(imp.getCalibration());
		impOut.getProcessor().setMinAndMax(0, max);
		return impOut;
	}

	/**
	 * Get the centroids of all the particles in real units
	 *
	 * @param imp
	 * @param particleLabels
	 * @param particleSizes
	 * @return double[][] containing all the particles' centroids
	 */
	public double[][] getCentroids(ImagePlus imp, int[][] particleLabels,
			long[] particleSizes) {
		final int nParticles = particleSizes.length;
		final int w = imp.getWidth();
		final int h = imp.getHeight();
		final int d = imp.getImageStackSize();
		double[][] sums = new double[nParticles][3];
		for (int z = 0; z < d; z++) {
			for (int y = 0; y < h; y++) {
				final int index = y * w;
				for (int x = 0; x < w; x++) {
					final int particle = particleLabels[z][index + x];
					sums[particle][0] += x;
					sums[particle][1] += y;
					sums[particle][2] += z;
				}
			}
		}
		Calibration cal = imp.getCalibration();
		double[][] centroids = new double[nParticles][3];
		for (int p = 0; p < nParticles; p++) {
			centroids[p][0] = cal.pixelWidth * sums[p][0] / particleSizes[p];
			centroids[p][1] = cal.pixelHeight * sums[p][1] / particleSizes[p];
			centroids[p][2] = cal.pixelDepth * sums[p][2] / particleSizes[p];
		}
		return centroids;
	}

	public double[] getVolumes(ImagePlus imp, long[] particleSizes) {
		Calibration cal = imp.getCalibration();
		final double voxelVolume = cal.pixelWidth * cal.pixelHeight
				* cal.pixelDepth;
		final int nLabels = particleSizes.length;
		double[] particleVolumes = new double[nLabels];
		for (int i = 0; i < nLabels; i++) {
			particleVolumes[i] = voxelVolume * particleSizes[i];
		}
		return particleVolumes;
	}

	/**
	 * Get particles, particle labels and particle sizes from a 3D ImagePlus
	 *
	 * @param imp
	 * Binary input image
	 * @param slicesPerChunk
	 * number of slices per chunk. 2 is generally good.
	 * @param minVol
	 * minimum volume particle to include
	 * @param maxVol
	 * maximum volume particle to include
	 * @param phase
	 * foreground or background (FORE or BACK)
	 * @param doExclude
	 * if true, remove particles touching sides of the stack
	 * @return Object[] {byte[][], int[][]} containing a binary workArray and
	 * particle labels.
	 */
	public Object[] getParticles(ImagePlus imp, int slicesPerChunk,
			double minVol, double maxVol, int phase, boolean doExclude) {
		byte[][] workArray = makeWorkArray(imp);
		return getParticles(imp, workArray, slicesPerChunk, minVol, maxVol,
				phase, doExclude);
	}

	public Object[] getParticles(ImagePlus imp, int slicesPerChunk,
			double minVol, double maxVol, int phase) {
		byte[][] workArray = makeWorkArray(imp);
		return getParticles(imp, workArray, slicesPerChunk, minVol, maxVol,
				phase, false);
	}

	public Object[] getParticles(ImagePlus imp, int slicesPerChunk, int phase) {
		byte[][] workArray = makeWorkArray(imp);
		double minVol = 0;
		double maxVol = Double.POSITIVE_INFINITY;
		return getParticles(imp, workArray, slicesPerChunk, minVol, maxVol,
				phase, false);
	}

	public Object[] getParticles(ImagePlus imp, byte[][] workArray,
			int slicesPerChunk, int phase, int method) {
		double minVol = 0;
		double maxVol = Double.POSITIVE_INFINITY;
		return getParticles(imp, workArray, slicesPerChunk, minVol, maxVol,
				phase, false);
	}

	public Object[] getParticles(ImagePlus imp, byte[][] workArray,
			int slicesPerChunk, double minVol, double maxVol, int phase) {
		return getParticles(imp, workArray, slicesPerChunk, minVol, maxVol,
				phase, false);
	}

	/**
	 * Get particles, particle labels and sizes from a workArray using an
	 * ImagePlus for scale information
	 *
	 * @param imp
	 * input binary image
	 * @param workArray
	 * work array
	 * @param slicesPerChunk
	 * number of slices to use for each chunk
	 * @param minVol
	 * minimum volume particle to include
	 * @param maxVol
	 * maximum volume particle to include
	 * @param phase
	 * FORE or BACK for foreground or background respectively
	 * @return Object[] array containing a binary workArray, particle labels and
	 * particle sizes
	 */
	public Object[] getParticles(ImagePlus imp, byte[][] workArray,
			int slicesPerChunk, double minVol, double maxVol, int phase,
			boolean doExclude) {
		if (phase == FORE) {
			this.sPhase = "foreground";
		} else if (phase == BACK) {
			this.sPhase = "background";
		} else {
			throw new IllegalArgumentException();
		}
		if (slicesPerChunk < 1) {
			throw new IllegalArgumentException();
		}
		// Set up the chunks
		final int nChunks = getNChunks(imp, slicesPerChunk);
		final int[][] chunkRanges = getChunkRanges(imp, nChunks, slicesPerChunk);
		final int[][] stitchRanges = getStitchRanges(imp, nChunks,
				slicesPerChunk);

		int[][] particleLabels = firstIDAttribution(imp, workArray, phase);
		final int nParticles = getParticleSizes(particleLabels).length;

		if (labelMethod == MULTI) {
			// connect particles within chunks
			final int nThreads = Runtime.getRuntime().availableProcessors();
			ConnectStructuresThread[] cptf = new ConnectStructuresThread[nThreads];
			for (int thread = 0; thread < nThreads; thread++) {
				cptf[thread] = new ConnectStructuresThread(thread, nThreads,
						imp, workArray, particleLabels, phase, nChunks,
						chunkRanges);
				cptf[thread].start();
			}
			try {
				for (int thread = 0; thread < nThreads; thread++) {
					cptf[thread].join();
				}
			} catch (InterruptedException ie) {
				IJ.error("A thread was interrupted.");
			}

			// connect particles between chunks
			if (nChunks > 1) {
				chunkString = ": stitching...";
				connectStructures(imp, workArray, particleLabels, phase,
						stitchRanges);
			}
		} else if (labelMethod == LINEAR) {
			joinStructures(imp, particleLabels, phase);
		} else if (labelMethod == MAPPED) {
			joinMappedStructures(imp, particleLabels, nParticles, phase);
		}
		filterParticles(imp, workArray, particleLabels, minVol, maxVol, phase);
		if (doExclude)
			excludeOnEdges(imp, particleLabels, workArray);
		minimiseLabels(particleLabels);
		long[] particleSizes = getParticleSizes(particleLabels);
		Object[] result = { workArray, particleLabels, particleSizes };
		return result;

	}

	/**
	 * Remove particles outside user-specified volume thresholds
	 *
	 * @param imp
	 * ImagePlus, used for calibration
	 * @param workArray
	 * binary foreground and background information
	 * @param particleLabels
	 * Packed 3D array of particle labels
	 * @param minVol
	 * minimum (inclusive) particle volume
	 * @param maxVol
	 * maximum (inclusive) particle volume
	 * @param phase
	 * phase we are interested in
	 */
	private void filterParticles(ImagePlus imp, byte[][] workArray,
			int[][] particleLabels, double minVol, double maxVol, int phase) {
		if (minVol == 0 && maxVol == Double.POSITIVE_INFINITY)
			return;
		final int d = imp.getImageStackSize();
		final int wh = workArray[0].length;
		long[] particleSizes = getParticleSizes(particleLabels);
		double[] particleVolumes = getVolumes(imp, particleSizes);
		byte flip = 0;
		if (phase == FORE) {
			flip = (byte) 0;
		} else {
			flip = (byte) 255;
		}
		for (int z = 0; z < d; z++) {
			for (int i = 0; i < wh; i++) {
				final int p = particleLabels[z][i];
				final double v = particleVolumes[p];
				if (v < minVol || v > maxVol) {
					workArray[z][i] = flip;
					particleLabels[z][i] = 0;
				}
			}
		}
	}

	/**
	 * Gets rid of redundant particle labels
	 *
	 * @param particleLabels
	 * @return
	 */
	private void minimiseLabels(int[][] particleLabels) {
		IJ.showStatus("Minimising labels...");
		final int d = particleLabels.length;
		long[] particleSizes = getParticleSizes(particleLabels);
		final int nLabels = particleSizes.length;
		int[] newLabel = new int[nLabels];
		int minLabel = 0;
		// find the minimised labels
		for (int i = 0; i < nLabels; i++) {
			if (particleSizes[i] > 0) {
				if (i == minLabel) {
					newLabel[i] = i;
					minLabel++;
					continue;
				} else {
					newLabel[i] = minLabel;
					particleSizes[minLabel] = particleSizes[i];
					particleSizes[i] = 0;
					minLabel++;
				}
			}
		}
		// now replace labels
		final int wh = particleLabels[0].length;
		for (int z = 0; z < d; z++) {
			IJ.showStatus("Replacing with minimised labels...");
			IJ.showProgress(z, d);
			int[] slice = particleLabels[z];
			for (int i = 0; i < wh; i++) {
				final int p = slice[i];
				if (p > 0) {
					slice[i] = newLabel[p];
				}
			}
		}
		return;
	}

	/**
	 * Scans edge voxels and set all touching particles to background
	 *
	 * @param particleLabels
	 * @param nLabels
	 * @param w
	 * @param h
	 * @param d
	 */
	private void excludeOnEdges(ImagePlus imp, int[][] particleLabels,
			byte[][] workArray) {
		final int w = imp.getWidth();
		final int h = imp.getHeight();
		final int d = imp.getImageStackSize();
		long[] particleSizes = getParticleSizes(particleLabels);
		final int nLabels = particleSizes.length;
		int[] newLabel = new int[nLabels];
		for (int i = 0; i < nLabels; i++)
			newLabel[i] = i;

		// scan faces
		// top and bottom faces
		for (int y = 0; y < h; y++) {
			final int index = y * w;
			for (int x = 0; x < w; x++) {
				final int pt = particleLabels[0][index + x];
				if (pt > 0)
					newLabel[pt] = 0;
				final int pb = particleLabels[d - 1][index + x];
				if (pb > 0)
					newLabel[pb] = 0;
			}
		}

		// west and east faces
		for (int z = 0; z < d; z++) {
			for (int y = 0; y < h; y++) {
				final int pw = particleLabels[z][y * w];
				final int pe = particleLabels[z][y * w + w - 1];
				if (pw > 0)
					newLabel[pw] = 0;
				if (pe > 0)
					newLabel[pe] = 0;
			}
		}

		// north and south faces
		final int lastRow = w * (h - 1);
		for (int z = 0; z < d; z++) {
			for (int x = 0; x < w; x++) {
				final int pn = particleLabels[z][x];
				final int ps = particleLabels[z][lastRow + x];
				if (pn > 0)
					newLabel[pn] = 0;
				if (ps > 0)
					newLabel[ps] = 0;
			}
		}

		// replace labels
		final int wh = w * h;
		for (int z = 0; z < d; z++) {
			for (int i = 0; i < wh; i++) {
				final int p = particleLabels[z][i];
				final int nL = newLabel[p];
				if (nL == 0) {
					particleLabels[z][i] = 0;
					workArray[z][i] = (byte) 0;
				}
			}
		}

		return;
	}

	/**
	 * Gets number of chunks needed to divide a stack into evenly-sized sets of
	 * slices.
	 *
	 * @param imp
	 * input image
	 * @param slicesPerChunk
	 * number of slices per chunk
	 * @return number of chunks
	 */
	public int getNChunks(ImagePlus imp, int slicesPerChunk) {
		final int d = imp.getImageStackSize();
		int nChunks = (int) Math.floor((double) d / (double) slicesPerChunk);

		int remainder = d % slicesPerChunk;

		if (remainder > 0) {
			nChunks++;
		}
		return nChunks;
	}

	/**
	 * Go through all pixels and assign initial particle label
	 *
	 * @param workArray
	 * byte[] array containing pixel values
	 * @param phase
	 * FORE or BACK for foreground of background respectively
	 * @return particleLabels int[] array containing label associating every
	 * pixel with a particle
	 */
	private int[][] firstIDAttribution(ImagePlus imp, final byte[][] workArray,
			final int phase) {
		final int w = imp.getWidth();
		final int h = imp.getHeight();
		final int d = imp.getImageStackSize();
		final int wh = w * h;
		IJ.showStatus("Finding " + sPhase + " structures");
		int[][] particleLabels = new int[d][wh];
		int ID = 1;

		if (phase == FORE) {
			for (int z = 0; z < d; z++) {
				for (int y = 0; y < h; y++) {
					final int rowIndex = y * w;
					for (int x = 0; x < w; x++) {
						final int arrayIndex = rowIndex + x;
						if (workArray[z][arrayIndex] == phase) {
							particleLabels[z][arrayIndex] = ID;
							int minTag = ID;
							// Find the minimum particleLabel in the
							// neighbouring pixels
							for (int vZ = z - 1; vZ <= z + 1; vZ++) {
								for (int vY = y - 1; vY <= y + 1; vY++) {
									for (int vX = x - 1; vX <= x + 1; vX++) {
										if (withinBounds(vX, vY, vZ, w, h, 0, d)) {
											final int offset = getOffset(vX,
													vY, w);
											if (workArray[vZ][offset] == phase) {
												final int tagv = particleLabels[vZ][offset];
												if (tagv != 0 && tagv < minTag) {
													minTag = tagv;
												}
											}
										}
									}
								}
							}
							// assign the smallest particle label from the
							// neighbours to the pixel
							particleLabels[z][arrayIndex] = minTag;
							// increment the particle label
							if (minTag == ID) {
								ID++;
							}
						}
					}
				}
				IJ.showProgress(z, d);
			}
			ID++;
		} else if (phase == BACK) {
			for (int z = 0; z < d; z++) {
				for (int y = 0; y < h; y++) {
					final int rowIndex = y * w;
					for (int x = 0; x < w; x++) {
						final int arrayIndex = rowIndex + x;
						if (workArray[z][arrayIndex] == phase) {
							particleLabels[z][arrayIndex] = ID;
							int minTag = ID;
							// Find the minimum particleLabel in the
							// neighbouring pixels
							int nX = x, nY = y, nZ = z;
							for (int n = 0; n < 7; n++) {
								switch (n) {
								case 0:
									break;
								case 1:
									nX = x - 1;
									break;
								case 2:
									nX = x + 1;
									break;
								case 3:
									nY = y - 1;
									nX = x;
									break;
								case 4:
									nY = y + 1;
									break;
								case 5:
									nZ = z - 1;
									nY = y;
									break;
								case 6:
									nZ = z + 1;
									break;
								}
								if (withinBounds(nX, nY, nZ, w, h, 0, d)) {
									final int offset = getOffset(nX, nY, w);
									if (workArray[nZ][offset] == phase) {
										final int tagv = particleLabels[nZ][offset];
										if (tagv != 0 && tagv < minTag) {
											minTag = tagv;
										}
									}
								}
							}
							// assign the smallest particle label from the
							// neighbours to the pixel
							particleLabels[z][arrayIndex] = minTag;
							// increment the particle label
							if (minTag == ID) {
								ID++;
							}
						}
					}
				}
				IJ.showProgress(z, d);
			}
			ID++;
		}
		return particleLabels;
	}

	/**
	 * Connect structures = minimisation of IDs
	 *
	 * @param workArray
	 * @param particleLabels
	 * @param phase
	 * foreground or background
	 * @param scanRanges
	 * int[][] listgetPixel(particleLabels, x, y, z, w, h, d);ing
	 * ranges to run connectStructures on
	 * @return particleLabels with all particles connected
	 */
	private void connectStructures(ImagePlus imp, final byte[][] workArray,
			int[][] particleLabels, final int phase, final int[][] scanRanges) {
		IJ.showStatus("Connecting " + sPhase + " structures" + chunkString);
		final int w = imp.getWidth();
		final int h = imp.getHeight();
		final int d = imp.getImageStackSize();
		for (int c = 0; c < scanRanges[0].length; c++) {
			final int sR0 = scanRanges[0][c];
			final int sR1 = scanRanges[1][c];
			final int sR2 = scanRanges[2][c];
			final int sR3 = scanRanges[3][c];
			if (phase == FORE) {
				for (int z = sR0; z < sR1; z++) {
					for (int y = 0; y < h; y++) {
						final int rowIndex = y * w;
						for (int x = 0; x < w; x++) {
							final int arrayIndex = rowIndex + x;
							if (workArray[z][arrayIndex] == phase
									&& particleLabels[z][arrayIndex] > 1) {
								int minTag = particleLabels[z][arrayIndex];
								// Find the minimum particleLabel in the
								// neighbours' pixels
								for (int vZ = z - 1; vZ <= z + 1; vZ++) {
									for (int vY = y - 1; vY <= y + 1; vY++) {
										for (int vX = x - 1; vX <= x + 1; vX++) {
											if (withinBounds(vX, vY, vZ, w, h,
													sR2, sR3)) {
												final int offset = getOffset(
														vX, vY, w);
												if (workArray[vZ][offset] == phase) {
													final int tagv = particleLabels[vZ][offset];
													if (tagv != 0
															&& tagv < minTag) {
														minTag = tagv;
													}
												}
											}
										}
									}
								}
								// Replacing particleLabel by the minimum
								// particleLabel found
								for (int vZ = z - 1; vZ <= z + 1; vZ++) {
									for (int vY = y - 1; vY <= y + 1; vY++) {
										for (int vX = x - 1; vX <= x + 1; vX++) {
											if (withinBounds(vX, vY, vZ, w, h,
													sR2, sR3)) {
												final int offset = getOffset(
														vX, vY, w);
												if (workArray[vZ][offset] == phase) {
													final int tagv = particleLabels[vZ][offset];
													if (tagv != 0
															&& tagv != minTag) {
														replaceLabel(
																particleLabels,
																tagv, minTag,
																sR2, sR3);
													}
												}
											}
										}
									}
								}
							}
						}
					}
					IJ.showStatus("Connecting foreground structures"
							+ chunkString);
					IJ.showProgress(z, d);
				}
			} else if (phase == BACK) {
				for (int z = sR0; z < sR1; z++) {
					for (int y = 0; y < h; y++) {
						final int rowIndex = y * w;
						for (int x = 0; x < w; x++) {
							final int arrayIndex = rowIndex + x;
							if (workArray[z][arrayIndex] == phase) {
								int minTag = particleLabels[z][arrayIndex];
								// Find the minimum particleLabel in the
								// neighbours' pixels
								int nX = x, nY = y, nZ = z;
								for (int n = 0; n < 7; n++) {
									switch (n) {
									case 0:
										break;
									case 1:
										nX = x - 1;
										break;
									case 2:
										nX = x + 1;
										break;
									case 3:
										nY = y - 1;
										nX = x;
										break;
									case 4:
										nY = y + 1;
										break;
									case 5:
										nZ = z - 1;
										nY = y;
										break;
									case 6:
										nZ = z + 1;
										break;
									}
									if (withinBounds(nX, nY, nZ, w, h, sR2, sR3)) {
										final int offset = getOffset(nX, nY, w);
										if (workArray[nZ][offset] == phase) {
											final int tagv = particleLabels[nZ][offset];
											if (tagv != 0 && tagv < minTag) {
												minTag = tagv;
											}
										}
									}
								}
								// Replacing particleLabel by the minimum
								// particleLabel found
								for (int n = 0; n < 7; n++) {
									switch (n) {
									case 0:
										nZ = z;
										break; // last switch block left nZ = z
										// + 1;
									case 1:
										nX = x - 1;
										break;
									case 2:
										nX = x + 1;
										break;
									case 3:
										nY = y - 1;
										nX = x;
										break;
									case 4:
										nY = y + 1;
										break;
									case 5:
										nZ = z - 1;
										nY = y;
										break;
									case 6:
										nZ = z + 1;
										break;
									}
									if (withinBounds(nX, nY, nZ, w, h, sR2, sR3)) {
										final int offset = getOffset(nX, nY, w);
										if (workArray[nZ][offset] == phase) {
											final int tagv = particleLabels[nZ][offset];
											if (tagv != 0 && tagv != minTag) {
												replaceLabel(particleLabels,
														tagv, minTag, sR2, sR3);
											}
										}
									}
								}
							}
						}
					}
					IJ.showStatus("Connecting background structures"
							+ chunkString);
					IJ.showProgress(z, d + 1);
				}
			}
		}
		return;
	}

	class ConnectStructuresThread extends Thread {
		final ImagePlus imp;

		final int thread, nThreads, nChunks, phase;

		final byte[][] workArray;

		final int[][] particleLabels;

		final int[][] chunkRanges;

		public ConnectStructuresThread(int thread, int nThreads, ImagePlus imp,
				byte[][] workArray, int[][] particleLabels, final int phase,
				int nChunks, int[][] chunkRanges) {
			this.imp = imp;
			this.thread = thread;
			this.nThreads = nThreads;
			this.workArray = workArray;
			this.particleLabels = particleLabels;
			this.phase = phase;
			this.nChunks = nChunks;
			this.chunkRanges = chunkRanges;
		}

		public void run() {
			for (int k = this.thread; k < this.nChunks; k += this.nThreads) {
				// assign singleChunkRange for chunk k from chunkRanges
				int[][] singleChunkRange = new int[4][1];
				for (int i = 0; i < 4; i++) {
					singleChunkRange[i][0] = this.chunkRanges[i][k];
				}
				chunkString = ": chunk " + (k + 1) + "/" + nChunks;
				connectStructures(this.imp, this.workArray,
						this.particleLabels, this.phase, singleChunkRange);
			}
		}
	}// ConnectStructuresThread

	/**
	 * Joins semi-labelled particles using a non-recursive algorithm
	 *
	 * @param imp
	 * @param particleLabels
	 */
	private void joinStructures(ImagePlus imp, int[][] particleLabels, int phase) {
		final int w = imp.getWidth();
		final int h = imp.getHeight();
		final int d = imp.getImageStackSize();
		long[] particleSizes = getParticleSizes(particleLabels);
		final int nBlobs = particleSizes.length;
		ArrayList<ArrayList<short[]>> particleLists = getParticleLists(
				particleLabels, nBlobs, w, h, d);
		switch (phase) {
		case FORE: {
			for (int b = 1; b < nBlobs; b++) {
				IJ.showStatus("Joining substructures...");
				IJ.showProgress(b, nBlobs);
				if (particleLists.get(b).isEmpty()) {
					continue;
				}

				for (int l = 0; l < particleLists.get(b).size(); l++) {
					final short[] voxel = particleLists.get(b).get(l);
					final int x = voxel[0];
					final int y = voxel[1];
					final int z = voxel[2];
					// find any neighbours with bigger labels
					for (int zN = z - 1; zN <= z + 1; zN++) {
						for (int yN = y - 1; yN <= y + 1; yN++) {
							final int index = yN * w;
							for (int xN = x - 1; xN <= x + 1; xN++) {
								if (!withinBounds(xN, yN, zN, w, h, d))
									continue;
								final int iN = index + xN;
								int p = particleLabels[zN][iN];
								if (p > b) {
									joinBlobs(b, p, particleLabels,
											particleLists, w);
								}
							}
						}
					}
				}
			}
		}
		case BACK: {
			for (int b = 1; b < nBlobs; b++) {
				IJ.showStatus("Joining substructures...");
				IJ.showProgress(b, nBlobs);
				if (particleLists.get(b).isEmpty()) {
					continue;
				}
				for (int l = 0; l < particleLists.get(b).size(); l++) {
					final short[] voxel = particleLists.get(b).get(l);
					final int x = voxel[0];
					final int y = voxel[1];
					final int z = voxel[2];
					// find any neighbours with bigger labels
					int xN = x, yN = y, zN = z;
					for (int n = 1; n < 7; n++) {
						switch (n) {
						case 1:
							xN = x - 1;
							break;
						case 2:
							xN = x + 1;
							break;
						case 3:
							yN = y - 1;
							xN = x;
							break;
						case 4:
							yN = y + 1;
							break;
						case 5:
							zN = z - 1;
							yN = y;
							break;
						case 6:
							zN = z + 1;
							break;
						}
						if (!withinBounds(xN, yN, zN, w, h, d))
							continue;
						final int iN = yN * w + xN;
						int p = particleLabels[zN][iN];
						if (p > b) {
							joinBlobs(b, p, particleLabels, particleLists, w);
						}
					}
				}
			}
		}
		}
		return;
	}

	public ArrayList<ArrayList<short[]>> getParticleLists(
			int[][] particleLabels, int nBlobs, int w, int h, int d) {
		ArrayList<ArrayList<short[]>> pL = new ArrayList<ArrayList<short[]>>(
				nBlobs);
		long[] particleSizes = getParticleSizes(particleLabels);
		ArrayList<short[]> background = new ArrayList<short[]>(0);
		pL.add(0, background);
		for (int b = 1; b < nBlobs; b++) {
			ArrayList<short[]> a = new ArrayList<short[]>(
					(int) particleSizes[b]);
			pL.add(b, a);
		}
		// add all the particle coordinates to the appropriate list
		for (short z = 0; z < d; z++) {
			IJ.showStatus("Listing substructures...");
			IJ.showProgress(z, d);
			final int[] sliceLabels = particleLabels[z];
			for (short y = 0; y < h; y++) {
				final int i = y * w;
				for (short x = 0; x < w; x++) {
					final int p = sliceLabels[i + x];
					if (p > 0) { // ignore background
						final short[] voxel = { x, y, z };
						pL.get(p).add(voxel);
					}
				}
			}
		}
		return pL;
	}

	/**
	 * Join particle p to particle b, relabelling p with b.
	 *
	 * @param b
	 * @param p
	 * @param particleLabels
	 * array of particle labels
	 * @param particleLists
	 * list of particle voxel coordinates
	 * @param w
	 * stack width
	 */
	public void joinBlobs(int b, int p, int[][] particleLabels,
			ArrayList<ArrayList<short[]>> particleLists, int w) {
		ListIterator<short[]> iterB = particleLists.get(p).listIterator();
		while (iterB.hasNext()) {
			short[] voxelB = iterB.next();
			particleLists.get(b).add(voxelB);
			final int iB = voxelB[1] * w + voxelB[0];
			particleLabels[voxelB[2]][iB] = b;
		}
		particleLists.get(p).clear();
	}

	private void joinMappedStructures(ImagePlus imp, int[][] particleLabels,
			int nParticles, int phase) {
		IJ.showStatus("Mapping structures and joining...");
		final int w = imp.getWidth();
		final int h = imp.getHeight();
		final int d = imp.getImageStackSize();

		ArrayList<HashSet<Integer>> map = new ArrayList<HashSet<Integer>>(
				nParticles + 1);

		int[] lut = new int[nParticles + 1];
		// set each label to be its own root
		final int initialCapacity = 1;
		for (int i = 0; i < nParticles + 1; i++) {
			lut[i] = i;
			Integer root = Integer.valueOf(i);
			HashSet<Integer> set = new HashSet<Integer>(initialCapacity);
			set.add(root);
			map.add(set);
		}

		// populate the first list with neighbourhoods
		int[] nbh = null;
		if (phase == FORE)
			nbh = new int[26];
		else if (phase == BACK)
			nbh = new int[6];
		for (int z = 0; z < d; z++) {
			IJ.showStatus("Building neighbourhood list");
			IJ.showProgress(z, d - 1);
			final int[] slice = particleLabels[z];
			for (int y = 0; y < h; y++) {
				final int yw = y * w;
				for (int x = 0; x < w; x++) {
					final int centre = slice[yw + x];
					// ignore background
					if (centre == 0)
						continue;
					if (phase == FORE)
						get26Neighborhood(nbh, particleLabels, x, y, z, w, h, d);
					else if (phase == BACK)
						get6Neighborhood(nbh, particleLabels, x, y, z, w, h, d);
					addNeighboursToMap(map, nbh, centre);
				}
			}
		}
		// map now contains for every value the set of first degree neighbours

		IJ.showStatus("Minimising list and generating LUT...");
		// place to store counts of each label
		int[] counter = new int[lut.length];

		// place to map lut values and targets
		//lutList lists the indexes which point to each transformed lutvalue
		//for quick updating
		ArrayList<HashSet<Integer>> lutList = new ArrayList<HashSet<Integer>>(
				nParticles);

		//initialise the lutList
		for (int i = 0; i <= nParticles; i++){
			HashSet<Integer> set = new HashSet<Integer>(2);
			lutList.add(set);
		}

		//set it up. ArrayList index is now the transformed value
		// list contains the lut indices that have the transformed value
		for (int i = 1; i < nParticles; i++){
			HashSet<Integer> list = lutList.get(lut[i]);
			list.add(Integer.valueOf(i));
		}

		// initialise LUT with minimal label in set
		updateLUTwithMinPosition(lut, map, lutList);

		// find the minimal position of each value
		findFirstAppearance(lut, map);

		// de-chain the lut array
		minimiseLutArray(lut);

		int duplicates = Integer.MAX_VALUE;
		boolean snowball = true;
		boolean merge = true;
		boolean update = true;
		boolean find = true;
		boolean minimise = true;
		boolean consistent = false;
		while ((duplicates > 0) && snowball && merge && update && find
				&& minimise && !consistent) {
			snowball = snowballLUT(lut, map, lutList);

			duplicates = countDuplicates(counter, map, lut);

			// merge duplicates
			merge = mergeDuplicates(map, counter, duplicates, lut, lutList);

			// update the LUT
			update = updateLUTwithMinPosition(lut, map, lutList);

			find = findFirstAppearance(lut, map);

			// minimise the LUT
			minimise = minimiseLutArray(lut);

			consistent = checkConsistence(lut, map);
		}

		// replace all labels with LUT values
		applyLUT(particleLabels, lut, w, h, d);
		IJ.showStatus("LUT applied");
	}

	private boolean checkConsistence(int[] lut, ArrayList<HashSet<Integer>> map) {
		final int l = lut.length;
		Integer val = null;
		for (int i = 1; i < l; i++) {
			val = Integer.valueOf(i);
			if (!map.get(lut[i]).contains(val))
				return false;
		}
		return true;
	}

	private boolean findFirstAppearance(int[] lut,
			ArrayList<HashSet<Integer>> map) {
		final int l = map.size();
		boolean changed = false;
		for (int i = 0; i < l; i++) {
			HashSet<Integer> set = map.get(i);
			for (Integer val : set) {
				// if the current lut value is greater
				// than the current position
				// update lut with current position
				final int v = val.intValue();
				if (lut[v] > i) {
					lut[v] = i;
					changed = true;
				}
			}
		}
		return changed;
	}

	private boolean updateLUTwithMinPosition(int[] lut,
			ArrayList<HashSet<Integer>> map, ArrayList<HashSet<Integer>> lutList) {
		final int l = lut.length;
		boolean changed = false;
		for (int i = 1; i < l; i++) {
			HashSet<Integer> set = map.get(i);
			if (set.isEmpty())
				continue;
			// find minimal value or lut value in the set
			int min = Integer.MAX_VALUE;
			int minLut = Integer.MAX_VALUE;
			for (Integer val : set) {
				int v = val.intValue();
				min = Math.min(min, v);
				minLut = Math.min(minLut, lut[v]);
			}
			// min now contains the smaller of the neighbours or their LUT
			// values
			min = Math.min(min, minLut);
			// add minimal value to lut
			HashSet<Integer> target = map.get(min);
			for (Integer val : set) {
				target.add(val);
				final int v = val.intValue();
				if (lut[v] > min)
					lut[v] = min;
			}
			set.clear();
			updateLUT(i, min, lut, lutList);
		}
		return changed;
	}

	private boolean mergeDuplicates(ArrayList<HashSet<Integer>> map,
			int[] counter, int duplicates, int[] lut,
			ArrayList<HashSet<Integer>> lutList) {
		boolean changed = false;
		// create a list of duplicate values to check for
		int[] dupList = new int[duplicates];
		final int l = counter.length;
		int dup = 0;
		for (int i = 1; i < l; i++) {
			if (counter[i] > 1)
				dupList[dup] = i;
			dup++;
		}

		// find duplicates hiding in sets which are greater than the lut
		// HashSet<Integer> set = null;
		// HashSet<Integer> target = null;
		// Iterator<Integer> iter = null;
		// Integer val = null;
		for (int i = 1; i < l; i++) {
			HashSet<Integer> set = map.get(i);
			if (set.isEmpty())
				continue;
			for (int d : dupList) {
				// if we are in the lut key of this value, continue
				final int lutValue = lut[d];
				if (lutValue == i)
					continue;
				// otherwise check to see if the non-lut set contains our dup
				if (set.contains(Integer.valueOf(d))) {
					// we found a dup, merge whole set back to lut
					changed = true;
					Iterator<Integer> iter = set.iterator();
					HashSet<Integer> target = map.get(lutValue);
					// if (target.isEmpty())
					// IJ.log("attempting to merge with empty target"
					// + lutValue);
					while (iter.hasNext()) {
						Integer val = iter.next();
						target.add(val);
						lut[val.intValue()] = lutValue;
					}
					// empty the set
					set.clear();
					updateLUT(i, lutValue, lut, lutList);
					// move to the next set
					break;

				}
			}
		}
		return changed;
	}

	/**
	 * Iterate backwards over map entries, moving set values to their new lut
	 * positions in the map. Updates LUT value of shifted values
	 *
	 * @param lut
	 * @param map
	 * @return false if nothing changed, true if something changed
	 */
	private boolean snowballLUT(final int[] lut,
			ArrayList<HashSet<Integer>> map, ArrayList<HashSet<Integer>> lutList) {
		// HashSet<Integer> set = null;
		// HashSet<Integer> target = null;
		boolean changed = false;
		for (int i = lut.length - 1; i > 0; i--) {
			IJ.showStatus("Snowballing labels...");
			IJ.showProgress(lut.length - i + 1, lut.length);
			final int lutValue = lut[i];
			if (lutValue < i) {
				changed = true;
				HashSet<Integer> set = map.get(i);
				HashSet<Integer> target = map.get(lutValue);
				// if (target.isEmpty())
				// IJ.log("merging with empty target " + lutValue);
				for (Integer n : set) {
					target.add(n);
					lut[n.intValue()] = lutValue;
				}
				// set is made empty
				// if later tests come across empty sets, then
				// must look up the lut to find the new location of the
				// neighbour network
				set.clear();
				// update lut so that anything pointing
				// to cleared set points to the new set
				updateLUT(i, lutValue, lut, lutList);
			}
		}
		return changed;
	}

	/**
	 * Replace old value with new value in LUT using map
	 *
	 * @param oldValue
	 * @param newValue
	 * @param lut
	 * @param lutlist
	 */
	private void updateLUT(int oldValue, int newValue, int[] lut,
			ArrayList<HashSet<Integer>> lutlist) {
		HashSet<Integer> list = lutlist.get(oldValue);
		HashSet<Integer> newList = lutlist.get(newValue);

		for (Integer in : list){
			lut[in.intValue()] = newValue;
			newList.add(in);
		}
		list.clear();
	}

	/**
	 * Find duplicated values and update the LUT
	 *
	 * @param counter
	 * @param map
	 * @param lut
	 * @return
	 */
	private int countDuplicates(int[] counter, ArrayList<HashSet<Integer>> map,
			int[] lut) {
		// reset to 0 the counter array
		final int l = counter.length;
		counter = new int[l];
		HashSet<Integer> set = null;
		for (int i = 1; i < map.size(); i++) {
			set = map.get(i);
			for (Integer val : set) {
				final int v = val.intValue();
				// every time a value is seen, log it
				counter[v]++;
				// update its LUT value if value was
				// found in a set with lower than current
				// lut value
				if (lut[v] > i)
					lut[v] = i;
			}
		}
		minimiseLutArray(lut);
		// all values should be seen only once,
		// count how many >1s there are.
		int count = 0;
		for (int i = 1; i < l; i++) {
			if (counter[i] > 1)
				count++;
		}
		return count;
	}

	/**
	 * Add all the neighbouring labels of a pixel to the map, except 0
	 * (background) and the pixel's own label, which is already in the map.
	 *
	 * The LUT gets updated with the minimum neighbour found, but this is only
	 * within the first neighbours and not the minimum label in the pixel's
	 * neighbour network
	 *
	 * @param map
	 * @param nbh
	 * @param centre
	 * current pixel's label
	 * @param lut
	 */
	private void addNeighboursToMap(ArrayList<HashSet<Integer>> map, int[] nbh,
			int centre) {
		HashSet<Integer> set = map.get(centre);
		final int l = nbh.length;
		for (int i = 0; i < l; i++) {
			final int val = nbh[i];
			// skip background and self-similar labels
			// adding them again is a redundant waste of time
			if (val == 0 || val == centre )
				continue;
			set.add(Integer.valueOf(val));
		}
	}

	private void applyLUT(int[][] particleLabels, final int[] lut, final int w,
			final int h, final int d) {
		for (int z = 0; z < d; z++) {
			IJ.showStatus("Applying LUT...");
			IJ.showProgress(z, d - 1);
			int[] slice = particleLabels[z];
			for (int y = 0; y < h; y++) {
				final int yw = y * w;
				for (int x = 0; x < w; x++) {
					final int i = yw + x;
					final int label = slice[i];
					if (label == 0)
						continue;
					slice[i] = lut[label];
				}
			}
		}
	}

	private boolean minimiseLutArray(int[] lutArray) {
		final int l = lutArray.length;
		boolean changed = false;
		for (int key = 1; key < l; key++) {
			final int value = lutArray[key];
			// this is a root
			if (value == key)
				continue;
			// otherwise update the current key with the
			// value from the referred key
			final int minValue = lutArray[value];
			lutArray[key] = minValue;
			changed = true;
		}
		return changed;
	}

	/**
	 * Get neighborhood of a pixel in a 3D image (0 border conditions)
	 *
	 * @param image
	 * 3D image (int[][])
	 * @param x
	 * x- coordinate
	 * @param y
	 * y- coordinate
	 * @param z
	 * z- coordinate (in image stacks the indexes start at 1)
	 * @return corresponding 26-pixels neighborhood (0 if out of image)
	 */
	private void get26Neighborhood(int[] neighborhood, final int[][] image,
			final int x, final int y, final int z, final int w, final int h,
			final int d) {
		// if (phase == FORE) {
		// int[] neighborhood = new int[26];

		neighborhood[0] = getPixel(image, x - 1, y - 1, z - 1, w, h, d);
		neighborhood[1] = getPixel(image, x, y - 1, z - 1, w, h, d);
		neighborhood[2] = getPixel(image, x + 1, y - 1, z - 1, w, h, d);

		neighborhood[3] = getPixel(image, x - 1, y, z - 1, w, h, d);
		neighborhood[4] = getPixel(image, x, y, z - 1, w, h, d);
		neighborhood[5] = getPixel(image, x + 1, y, z - 1, w, h, d);

		neighborhood[6] = getPixel(image, x - 1, y + 1, z - 1, w, h, d);
		neighborhood[7] = getPixel(image, x, y + 1, z - 1, w, h, d);
		neighborhood[8] = getPixel(image, x + 1, y + 1, z - 1, w, h, d);

		neighborhood[9] = getPixel(image, x - 1, y - 1, z, w, h, d);
		neighborhood[10] = getPixel(image, x, y - 1, z, w, h, d);
		neighborhood[11] = getPixel(image, x + 1, y - 1, z, w, h, d);

		neighborhood[12] = getPixel(image, x - 1, y, z, w, h, d);
		// neighborhood[13] = getPixel(image, x, y, z, w, h, d);
		neighborhood[13] = getPixel(image, x + 1, y, z, w, h, d);

		neighborhood[14] = getPixel(image, x - 1, y + 1, z, w, h, d);
		neighborhood[15] = getPixel(image, x, y + 1, z, w, h, d);
		neighborhood[16] = getPixel(image, x + 1, y + 1, z, w, h, d);

		neighborhood[17] = getPixel(image, x - 1, y - 1, z + 1, w, h, d);
		neighborhood[18] = getPixel(image, x, y - 1, z + 1, w, h, d);
		neighborhood[19] = getPixel(image, x + 1, y - 1, z + 1, w, h, d);

		neighborhood[20] = getPixel(image, x - 1, y, z + 1, w, h, d);
		neighborhood[21] = getPixel(image, x, y, z + 1, w, h, d);
		neighborhood[22] = getPixel(image, x + 1, y, z + 1, w, h, d);

		neighborhood[23] = getPixel(image, x - 1, y + 1, z + 1, w, h, d);
		neighborhood[24] = getPixel(image, x, y + 1, z + 1, w, h, d);
		neighborhood[25] = getPixel(image, x + 1, y + 1, z + 1, w, h, d);

		// return neighborhood;
	}

	private void get6Neighborhood(int[] neighborhood, final int[][] image,
			final int x, final int y, final int z, final int w, final int h,
			final int d) {
		// int[] neighborhood = new int[6];
		neighborhood[0] = getPixel(image, x - 1, y, z, w, h, d);
		neighborhood[1] = getPixel(image, x, y - 1, z, w, h, d);
		neighborhood[2] = getPixel(image, x, y, z - 1, w, h, d);

		neighborhood[3] = getPixel(image, x + 1, y, z, w, h, d);
		neighborhood[4] = getPixel(image, x, y + 1, z, w, h, d);
		neighborhood[5] = getPixel(image, x, y, z + 1, w, h, d);
		// return neighborhood;
	}

	/* ----------------------------------------------------------------------- */
	/**
	 * Get pixel in 3D image (0 border conditions)
	 *
	 * @param image
	 * 3D image
	 * @param x
	 * x- coordinate
	 * @param y
	 * y- coordinate
	 * @param z
	 * z- coordinate (in image stacks the indexes start at 1)
	 * @return corresponding pixel (0 if out of image)
	 */
	private int getPixel(final int[][] image, final int x, final int y,
			final int z, final int w, final int h, final int d) {
		if (withinBounds(x, y, z, w, h, d))
			return image[z][x + y * w];
		else
			return 0;
	} /* end getPixel */

	/**
	 * Create a work array
	 *
	 * @return byte[] work array
	 */
	private byte[][] makeWorkArray(ImagePlus imp) {
		final int s = imp.getStackSize();
		final int p = imp.getWidth() * imp.getHeight();
		byte[][] workArray = new byte[s][p];
		ImageStack stack = imp.getStack();
		for (int z = 0; z < s; z++) {
			ImageProcessor ip = stack.getProcessor(z + 1);
			for (int i = 0; i < p; i++) {
				workArray[z][i] = (byte) ip.get(i);
			}
		}
		return workArray;
	}

	/**
	 * Get a 2 d array that defines the z-slices to scan within while connecting
	 * particles within chunkified stacks.
	 *
	 * @param nC
	 * number of chunks
	 * @return scanRanges int[][] containing 4 limits: int[0][] - start of outer
	 * for; int[1][] end of outer for; int[3][] start of inner for;
	 * int[4] end of inner 4. Second dimension is chunk number.
	 */
	public int[][] getChunkRanges(ImagePlus imp, int nC, int slicesPerChunk) {
		final int nSlices = imp.getImageStackSize();
		int[][] scanRanges = new int[4][nC];
		scanRanges[0][0] = 0; // the first chunk starts at the first (zeroth)
		// slice
		scanRanges[2][0] = 0; // and that is what replaceLabel() will work on
		// first

		if (nC == 1) {
			scanRanges[1][0] = nSlices;
			scanRanges[3][0] = nSlices;
		} else if (nC > 1) {
			scanRanges[1][0] = slicesPerChunk;
			scanRanges[3][0] = slicesPerChunk;

			for (int c = 1; c < nC; c++) {
				for (int i = 0; i < 4; i++) {
					scanRanges[i][c] = scanRanges[i][c - 1] + slicesPerChunk;
				}
			}
			// reduce the last chunk to nSlices
			scanRanges[1][nC - 1] = nSlices;
			scanRanges[3][nC - 1] = nSlices;
		}
		return scanRanges;
	}

	/**
	 * Return scan ranges for stitching. The first 2 values for each chunk are
	 * the first slice of the next chunk and the last 2 values are the range
	 * through which to replaceLabels()
	 *
	 * Running replace labels over incrementally increasing volumes as chunks
	 * are added is OK (for 1st interface connect chunks 0 & 1, for 2nd connect
	 * chunks 0, 1, 2, etc.)
	 *
	 * @param nC
	 * number of chunks
	 * @return scanRanges list of scan limits for connectStructures() to stitch
	 * chunks back together
	 */
	private int[][] getStitchRanges(ImagePlus imp, int nC, int slicesPerChunk) {
		final int nSlices = imp.getImageStackSize();
		if (nC < 2) {
			return null;
		}
		int[][] scanRanges = new int[4][3 * (nC - 1)]; // there are nC - 1
		// interfaces

		for (int c = 0; c < nC - 1; c++) {
			scanRanges[0][c] = (c + 1) * slicesPerChunk;
			scanRanges[1][c] = (c + 1) * slicesPerChunk + 1;
			scanRanges[2][c] = c * slicesPerChunk; // forward and reverse
			// algorithm
			// scanRanges[2][c] = 0; //cumulative algorithm - reliable but O�
			// hard
			scanRanges[3][c] = (c + 2) * slicesPerChunk;
		}
		// stitch back
		for (int c = nC - 1; c < 2 * (nC - 1); c++) {
			scanRanges[0][c] = (2 * nC - c - 2) * slicesPerChunk - 1;
			scanRanges[1][c] = (2 * nC - c - 2) * slicesPerChunk;
			scanRanges[2][c] = (2 * nC - c - 3) * slicesPerChunk;
			scanRanges[3][c] = (2 * nC - c - 1) * slicesPerChunk;
		}
		// stitch forwards (paranoid third pass)
		for (int c = 2 * (nC - 1); c < 3 * (nC - 1); c++) {
			scanRanges[0][c] = (-2 * nC + c + 3) * slicesPerChunk;
			scanRanges[1][c] = (-2 * nC + c + 3) * slicesPerChunk + 1;
			scanRanges[2][c] = (-2 * nC + c + 2) * slicesPerChunk;
			scanRanges[3][c] = (-2 * nC + c + 4) * slicesPerChunk;
		}
		for (int i = 0; i < scanRanges.length; i++) {
			for (int c = 0; c < scanRanges[i].length; c++) {
				if (scanRanges[i][c] > nSlices) {
					scanRanges[i][c] = nSlices;
				}
			}
		}
		scanRanges[3][nC - 2] = nSlices;
		return scanRanges;
	}

	/**
	 * Check to see if the pixel at (m,n,o) is within the bounds of the current
	 * stack
	 *
	 * @param m
	 * x co-ordinate
	 * @param n
	 * y co-ordinate
	 * @param o
	 * z co-ordinate
	 * @param startZ
	 * first Z coordinate to use
	 *
	 * @param endZ
	 * last Z coordinate to use
	 *
	 * @return True if the pixel is within the bounds of the current stack
	 */
	private boolean withinBounds(int m, int n, int o, int w, int h, int startZ,
			int endZ) {
		return (m >= 0 && m < w && n >= 0 && n < h && o >= startZ && o < endZ);
	}

	private boolean withinBounds(final int m, final int n, final int o,
			final int w, final int h, final int d) {
		return (m >= 0 && m < w && n >= 0 && n < h && o >= 0 && o < d);
	}

	/**
	 * Find the offset within a 1D array given 2 (x, y) offset values
	 *
	 * @param m
	 * x difference
	 * @param n
	 * y difference
	 *
	 * @return Integer offset for looking up pixel in work array
	 */
	private int getOffset(int m, int n, int w) {
		return m + n * w;
	}

	/**
	 * Check whole array replacing m with n
	 *
	 * @param m
	 * value to be replaced
	 * @param n
	 * new value
	 * @param startZ
	 * first z coordinate to check
	 * @param endZ
	 * last+1 z coordinate to check
	 */
	public void replaceLabel(int[][] particleLabels, final int m, int n,
			int startZ, final int endZ) {
		final int s = particleLabels[0].length;
		for (int z = startZ; z < endZ; z++) {
			for (int i = 0; i < s; i++)
				if (particleLabels[z][i] == m) {
					particleLabels[z][i] = n;
				}
		}
	}

	/**
	 * Check whole array replacing m with n
	 *
	 * @param m
	 * value to be replaced
	 * @param n
	 * new value
	 * @param startZ
	 * first z coordinate to check
	 * @param endZ
	 * last+1 z coordinate to check
	 * @param multithreaded
	 * true if label replacement should happen in multiple threads
	 */
	public void replaceLabel(final int[][] particleLabels, final int m,
			final int n, int startZ, final int endZ, final boolean multithreaded) {
		if (!multithreaded) {
			replaceLabel(particleLabels, m, n, startZ, endZ);
			return;
		}
		final int s = particleLabels[0].length;
		final AtomicInteger ai = new AtomicInteger(startZ);
		Thread[] threads = Multithreader.newThreads();
		for (int thread = 0; thread < threads.length; thread++) {
			threads[thread] = new Thread(new Runnable() {
				public void run() {
					for (int z = ai.getAndIncrement(); z < endZ; z = ai
							.getAndIncrement()) {
						for (int i = 0; i < s; i++)
							if (particleLabels[z][i] == m) {
								particleLabels[z][i] = n;
							}
					}
				}
			});
		}
		Multithreader.startAndJoin(threads);
	}

	/**
	 * Get the sizes of all the particles as a voxel count
	 *
	 * @param particleLabels
	 * @return particleSizes
	 */
	public long[] getParticleSizes(final int[][] particleLabels) {
		IJ.showStatus("Getting " + sPhase + " particle sizes");
		final int d = particleLabels.length;
		final int wh = particleLabels[0].length;
		// find the highest value particleLabel
		int maxParticle = 0;
		for (int z = 0; z < d; z++) {
			final int[] slice = particleLabels[z];
			for (int i = 0; i < wh; i++) {
				maxParticle = Math.max(maxParticle, slice[i]);
			}
		}

		long[] particleSizes = new long[maxParticle + 1];
		for (int z = 0; z < d; z++) {
			final int[] slice = particleLabels[z];
			for (int i = 0; i < wh; i++) {
				particleSizes[slice[i]]++;
			}
			IJ.showProgress(z, d);
		}
		return particleSizes;
	}


	/**
	 * Return the value of this instance's labelMethod field
	 *
	 * @return the label method as int
	 */
	public int getLabelMethod() {
		return labelMethod;
	}

	/**
	 * Set the value of this instance's labelMethod field
	 *
	 * @param label
	 * one of ParticleCounter.MULTI or .LINEAR
	 */
	public void setLabelMethod(int label) {
		if (label != MULTI && label != LINEAR && label != MAPPED) {
			throw new IllegalArgumentException();
		}
		labelMethod = label;
		return;
	}

	public boolean dialogItemChanged(GenericDialog gd, AWTEvent e) {
		if (!DialogModifier.allNumbersValid(gd.getNumericFields()))
			return false;
		Vector<?> choices = gd.getChoices();
		Vector<?> checkboxes = gd.getCheckboxes();
		Vector<?> numbers = gd.getNumericFields();
		// link algorithm choice to chunk size field
		Choice choice = (Choice) choices.get(1);
		TextField num = (TextField) numbers.get(5);
		if (choice.getSelectedItem().contentEquals("Multithreaded")) {
			num.setEnabled(true);
		} else {
			num.setEnabled(false);
		}
		// link show stack 3d to volume resampling
		Checkbox box = (Checkbox) checkboxes.get(15);
		TextField numb = (TextField) numbers.get(4);
		if (box.getState()) {
			numb.setEnabled(true);
		} else {
			numb.setEnabled(false);
		}
		// link show surfaces, gradient choice and split value
		Checkbox surfbox = (Checkbox) checkboxes.get(11);
		Choice col = (Choice) choices.get(0);
		TextField split = (TextField) numbers.get(3);
		if (!surfbox.getState()) {
			col.setEnabled(false);
			split.setEnabled(false);
		} else {
			col.setEnabled(true);
			if (col.getSelectedIndex() == 1) {
				split.setEnabled(true);
			} else {
				split.setEnabled(false);
			}
		}
		DialogModifier.registerMacroValues(gd, gd.getComponents());
		return true;
	}
}