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;
}
}