package gdsc.utils;
import java.awt.AWTEvent;
import java.awt.Color;
import java.awt.Point;
import java.awt.Rectangle;
import java.awt.image.IndexColorModel;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.HashMap;
import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
import org.apache.commons.math3.analysis.MultivariateVectorFunction;
import org.apache.commons.math3.fitting.leastsquares.LeastSquaresBuilder;
import org.apache.commons.math3.fitting.leastsquares.LeastSquaresOptimizer.Optimum;
import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem;
import org.apache.commons.math3.fitting.leastsquares.LevenbergMarquardtOptimizer;
import org.apache.commons.math3.linear.DiagonalMatrix;
import org.apache.commons.math3.optim.ConvergenceChecker;
import org.apache.commons.math3.optim.PointVectorValuePair;
import org.apache.commons.math3.optim.SimplePointChecker;
import org.apache.commons.math3.util.Precision;
import gdsc.UsageTracker;
/*-----------------------------------------------------------------------------
* GDSC Plugins for ImageJ
*
* Copyright (C) 2011 Alex Herbert
* Genome Damage and Stability Centre
* University of Sussex, UK
*
* 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 2 of the License, or
* (at your option) any later version.
*---------------------------------------------------------------------------*/
import ij.IJ;
import ij.ImagePlus;
import ij.ImageStack;
import ij.WindowManager;
import ij.gui.DialogListener;
import ij.gui.GenericDialog;
import ij.gui.Overlay;
import ij.gui.PointRoi;
import ij.gui.PolygonRoi;
import ij.gui.Roi;
import ij.plugin.LutLoader;
import ij.plugin.ZProjector;
import ij.plugin.filter.EDM;
import ij.plugin.filter.ExtendedPlugInFilter;
import ij.plugin.filter.PlugInFilterRunner;
import ij.process.Blitter;
import ij.process.ByteProcessor;
import ij.process.FloatBlitter;
import ij.process.FloatPolygon;
import ij.process.FloatProcessor;
import ij.process.ImageProcessor;
import ij.process.ShortProcessor;
/**
* Outlines a circular cell using the optimal path through a membrane scoring map.
* <p>
* The centre of the cell must be chosen and an approximate cell radius. Only pixels within a range of the cell radius
* are processed. The score map is created using 2D convolutions of a curve through a NxN square. The curve is
* constructed as a difference-of-Gaussians approximation to the Laplacian allowing dark/light edge detection. The curve
* arc uses the cell radius and is rotated to produce membrane edge projections around the centre point.
* <p>
* The maximum intensity membrane projection is weighted using the current elliptical fit of the cell (initialised as a
* circle using the cell radius). A polygon outline is constructed from the projection using the highest value in each
* 10 degree segment around the selected centre. An elliptical shape consisting of two ellipses back-to-back (allowing
* egg-like shapes) is then fitted to the polygon outline and the process iterated.
*/
public class Cell_Outliner implements ExtendedPlugInFilter, DialogListener
{
private int flags = DOES_8G + DOES_16 + DOES_32 + NO_CHANGES;
private static final String TITLE = "Cell Outliner";
private static int cellRadius = 25;
private static double tolerance = 0.8;
private static boolean darkEdge = true;
private static int kernelWidth = 13;
private static double kernelSmoothing = 1;
private static int polygonSmoothing = 1;
private static double weightingGamma = 3;
private static int iterations = 3;
private static boolean ellipticalFit = false;
private static int dilate = 0;
private boolean moreOptions = false;
private boolean buildMaskOutput = false;
private boolean processAllFrames = false;
private boolean processAllSlices = false;
private boolean debug = false;
private ImagePlus imp;
private Rectangle bounds;
private PointRoi roi;
private int[] xpoints, ypoints;
private ArrayList<Integer> rotationAngles = new ArrayList<Integer>();
HashMap<Integer, float[]> kernels = null;
HashMap<Integer, FloatProcessor> convolved = null;
private int halfWidth;
private double maxDistance2;
/*
* (non-Javadoc)
*
* @see ij.plugin.filter.PlugInFilter#setup(java.lang.String, ij.ImagePlus)
*/
public int setup(String arg, ImagePlus imp)
{
UsageTracker.recordPlugin(this.getClass(), arg);
if (imp == null)
{
return DONE;
}
if (imp.getRoi() == null || imp.getRoi().getType() != ij.gui.PolygonRoi.POINT)
{
IJ.error("Please select a centre point using the ROI tool");
return DONE;
}
this.imp = imp;
kernels = null;
roi = (PointRoi) imp.getRoi();
xpoints = roi.getPolygon().xpoints;
ypoints = roi.getPolygon().ypoints;
return flags;
}
/*
* (non-Javadoc)
*
* @see ij.plugin.filter.ExtendedPlugInFilter#showDialog(ij.ImagePlus, java.lang.String,
* ij.plugin.filter.PlugInFilterRunner)
*/
public int showDialog(ImagePlus imp, String command, PlugInFilterRunner pfr)
{
GenericDialog gd = new GenericDialog(TITLE);
moreOptions = IJ.shiftKeyDown();
gd.addSlider("Cell_radius", 10, 30, cellRadius);
gd.addSlider("Tolerance", 0.5, 1.5, tolerance);
gd.addSlider("Kernel_width", 7, 19, kernelWidth);
gd.addCheckbox("Dark_edge", darkEdge);
gd.addSlider("Kernel_smoothing", 1.0, 5.5, kernelSmoothing);
gd.addSlider("Polygon_smoothing", 0, 5, polygonSmoothing);
gd.addSlider("Weighting_gamma", 0.05, 5, weightingGamma);
gd.addSlider("Iterations", 1, 10, iterations);
gd.addCheckbox("Show_elliptical_fit", ellipticalFit);
gd.addSlider("Dilate", 0, 5, dilate);
if (moreOptions)
gd.addCheckbox("Debug", debug);
else
debug = false;
gd.addHelp(gdsc.help.URL.UTILITY);
gd.addPreviewCheckbox(pfr);
gd.addDialogListener(this);
gd.showDialog();
// Final run: stop debug mode
debug = false;
if (gd.wasCanceled() || !dialogItemChanged(gd, null))
{
imp.setOverlay(null); // Clear preview overlay
return DONE;
}
if (imp.getNFrames() > 1 || imp.getNSlices() > 1)
{
gd = new GenericDialog("Process Stack?");
gd.addMessage("Process multiple slices (" + (imp.getNFrames() * imp.getNSlices()) + ")");
if (imp.getNFrames() > 1)
gd.addCheckbox("All_frames (" + imp.getNFrames() + ")", false);
if (imp.getNSlices() > 1)
gd.addCheckbox("All_slices (" + imp.getNSlices() + ")", false);
gd.enableYesNoCancel("Yes", "No");
gd.showDialog();
if (gd.wasCanceled())
return DONE;
if (imp.getNFrames() > 1)
processAllFrames = gd.getNextBoolean();
if (imp.getNSlices() > 1)
processAllSlices = gd.getNextBoolean();
}
// Switch from the image overlay to mask output
buildMaskOutput = true;
return flags; //IJ.setupDialog(imp, flags);
}
/*
* (non-Javadoc)
*
* @see ij.gui.DialogListener#dialogItemChanged(ij.gui.GenericDialog, java.awt.AWTEvent)
*/
public boolean dialogItemChanged(GenericDialog gd, AWTEvent e)
{
int oldCellRadius = cellRadius;
double oldTolerance = tolerance;
boolean oldDarkEdge = darkEdge;
double oldKernelSmoothing = kernelSmoothing;
int oldKernelWidth = kernelWidth;
cellRadius = (int) gd.getNextNumber();
tolerance = gd.getNextNumber();
kernelWidth = (int) gd.getNextNumber();
darkEdge = gd.getNextBoolean();
kernelSmoothing = gd.getNextNumber();
polygonSmoothing = (int) gd.getNextNumber();
weightingGamma = gd.getNextNumber();
iterations = (int) gd.getNextNumber();
ellipticalFit = gd.getNextBoolean();
dilate = (int) gd.getNextNumber();
if (moreOptions)
debug = gd.getNextBoolean();
// Round down to nearest odd number and check minimum size
if (kernelWidth % 2 == 0)
kernelWidth--;
if (kernelWidth < 3)
kernelWidth = 3;
if (cellRadius < 1)
cellRadius = 1;
if (iterations < 1)
iterations = 1;
// Check if the convolutions need recalculating
if (oldCellRadius != cellRadius || oldTolerance != tolerance || oldKernelWidth != kernelWidth ||
oldDarkEdge != darkEdge || oldKernelSmoothing != kernelSmoothing)
kernels = null;
return true;
}
/*
* (non-Javadoc)
*
* @see ij.plugin.filter.ExtendedPlugInFilter#setNPasses(int)
*/
public void setNPasses(int nPasses)
{
// Do nothing
}
/*
* (non-Javadoc)
*
* @see ij.plugin.filter.PlugInFilter#run(ij.process.ImageProcessor)
*/
public void run(ImageProcessor inputProcessor)
{
// This can be called by the preview or as the final run
// Final: Create a new mask image
if (buildMaskOutput)
{
int w = inputProcessor.getWidth();
int h = inputProcessor.getHeight();
boolean useShort = (xpoints.length > 255);
// Build a list of the frames to process
int channel = imp.getChannel();
int[] frames = (processAllFrames) ? sequence(imp.getNFrames()) : new int[] { imp.getFrame() };
int[] slices = (processAllSlices) ? sequence(imp.getNSlices()) : new int[] { imp.getSlice() };
int nFrames = frames.length;
int nSlices = slices.length;
int size = nFrames * nSlices;
boolean resetConvolved = (size > 1);
ImageStack stack = new ImageStack(w, h, size);
ImageStack inputStack = this.imp.getImageStack();
int n = 1;
for (int frame : frames)
{
for (int slice : slices)
{
String label = "t" + frame + "z" + slice;
IJ.showStatus("Processing " + label);
IJ.showProgress(n - 1, size);
ImageProcessor maskIp = (useShort) ? new ShortProcessor(w, h) : new ByteProcessor(w, h);
int index = imp.getStackIndex(channel, slice, frame);
ImageProcessor ip = inputStack.getProcessor(index);
if (resetConvolved)
convolved = null;
PolygonRoi[] cells = findCells(ip);
for (int i = 0; i < cells.length; i++)
{
PolygonRoi cell = cells[i];
Rectangle b = cell.getBounds();
cell.setLocation(0, 0); // Remove the cell offset to allow it to be created
ByteProcessor cellIp = createFilledCell(b.width, b.height, cell, i + 1);
maskIp.copyBits(cellIp, b.x, b.y, Blitter.COPY_ZERO_TRANSPARENT);
}
stack.setPixels(maskIp.getPixels(), n);
stack.setSliceLabel(label, n);
n++;
}
}
// Create the output image
String title = imp.getTitle() + " Cell Outline";
ImagePlus maskImp = WindowManager.getImage(title);
if (maskImp == null)
{
maskImp = new ImagePlus(title, stack);
maskImp.setOpenAsHyperStack(true);
maskImp.setDimensions(1, nSlices, nFrames);
maskImp.show();
}
else
{
maskImp.setStack(stack, 1, nSlices, nFrames);
}
maskImp.setRoi(roi);
applyColorModel(maskImp);
maskImp.setDisplayRange(0, xpoints.length);
maskImp.updateAndDraw();
}
// Preview: Draw the polygon outline as an overlay
else
{
Overlay overlay = new Overlay();
overlay.setStrokeColor(Color.green);
overlay.setFillColor(null);
PolygonRoi[] cells = findCells(inputProcessor);
IJ.showStatus("");
IJ.showProgress(1);
if (cells == null)
return;
for (PolygonRoi cell : cells)
{
overlay.add(cell);
}
imp.setOverlay(overlay);
}
}
private int[] sequence(int n)
{
int[] s = new int[n];
for (int i = 0; i < n; i++)
s[i] = i + 1;
return s;
}
private PolygonRoi[] findCells(ImageProcessor inputProcessor)
{
// Limit processing to where it is needed
bounds = createBounds(inputProcessor);
ImageProcessor ip = inputProcessor.duplicate();
ip.setRoi(bounds);
ip = ip.crop();
if (kernels == null)
{
kernels = createKernels();
convolved = null;
}
if (convolved == null)
{
convolved = convolveImage(ip, kernels);
//showConvolvedImages(convolved);
}
if (Thread.currentThread().isInterrupted())
return null;
FloatProcessor combinedIp = null;
Blitter b = null;
ImagePlus combinedImp = null;
if (debug)
{
combinedIp = new FloatProcessor(ip.getWidth(), ip.getHeight());
b = new FloatBlitter(combinedIp);
combinedImp = displayImage(combinedIp, "Combined edge projection");
}
PolygonRoi[] cells = new PolygonRoi[xpoints.length];
if (!this.buildMaskOutput)
IJ.showStatus("Finding cells ...");
// Process each point
for (int n = 0; n < xpoints.length; n++)
{
if (!this.buildMaskOutput)
IJ.showProgress(n, xpoints.length);
int cx = xpoints[n] - bounds.x;
int cy = ypoints[n] - bounds.y;
// Restrict bounds using the cell radius and tolerance
int extra = (int) Math.ceil(cellRadius + cellRadius * tolerance + kernelWidth / 2);
int minx = Math.max(0, cx - extra);
int miny = Math.max(0, cy - extra);
int maxx = Math.min(ip.getWidth() - 1, cx + extra);
int maxy = Math.min(ip.getHeight() - 1, cy + extra);
Rectangle pointBounds = new Rectangle(minx, miny, maxx - minx + 1, maxy - miny + 1);
// Calculate the angle
FloatProcessor angle = createAngleProcessor(ip, cx, cy, pointBounds);
//if (debug) displayImage(angle, "Angle");
if (Thread.currentThread().isInterrupted())
return null;
FloatProcessor edgeProjection = computeEdgeProjection(convolved, cx, cy, pointBounds, angle);
// Initialise the edge as a circle.
PolygonRoi cell = null;
double[] params = { cx - minx, cy - miny, cellRadius, cellRadius, cellRadius, 0 };
double range = cellRadius * 0.9;
// Iterate to find the best cell outline
boolean returnEllipticalFit = ellipticalFit;
for (int iter = 0; iter < iterations; iter++)
{
// Use the current elliptical edge to define the weights for the edge projection
FloatProcessor weights = createWeightMap(pointBounds, params, range);
if (Thread.currentThread().isInterrupted())
return null;
if (debug)
displayImage(weights, "Weight map");
FloatProcessor weightedEdgeProjection = applyWeights(edgeProjection, weights);
if (debug)
{
b.copyBits(weightedEdgeProjection, pointBounds.x, pointBounds.y, Blitter.ADD);
combinedIp.resetMinAndMax();
combinedImp.updateAndDraw();
displayImage(weightedEdgeProjection, "Weighted edge projection");
}
cell = findPolygonalCell((int) Math.round(params[0]), (int) Math.round(params[1]),
weightedEdgeProjection, angle);
FloatProcessor weightMap = weightedEdgeProjection; // weights
double[] newParams = fitPolygonalCell(cell, params, weightMap, angle);
if (newParams == null)
{
returnEllipticalFit = false;
break;
}
// Set the parameters for the weight map
params = newParams;
// Update the range to become smaller for a tighter fit
//range = Math.max(3, range * 0.9);
}
// Return either the fitted elliptical cell or the last polygon outline
if (returnEllipticalFit)
{
EllipticalCell e = new EllipticalCell();
FloatPolygon ellipse = e.drawEllipse(params);
cell = new PolygonRoi(ellipse.xpoints, ellipse.ypoints, ellipse.npoints, PolygonRoi.POLYGON);
}
PolygonRoi finalCell = cell;
if (dilate > 0)
{
// Dilate the cell and then trace the new outline
ByteProcessor bp = new ByteProcessor(pointBounds.width, pointBounds.height);
bp.setColor(CELL & 0xff);
bp.draw(cell);
for (int i = 0; i < dilate; i++)
dilate(bp);
cell = traceOutline(bp, params[0], params[1]);
if (cell != null)
finalCell = cell;
}
Rectangle pos = finalCell.getBounds();
// Does not work in IJ 1.46+
//finalCell.setLocation(pos.x + bounds.x + pointBounds.x, pos.y + bounds.y + pointBounds.y);
// Create a new Polygon with the correct coordinates. This is required since IJ 1.46
// since setting the location is not passed through when drawing an overlay
int[] xCoords = finalCell.getXCoordinates();
int[] yCoords = finalCell.getYCoordinates();
int nPoints = finalCell.getNCoordinates();
for (int i = 0; i < nPoints; i++)
{
xCoords[i] += pos.x + bounds.x + pointBounds.x;
yCoords[i] += pos.y + bounds.y + pointBounds.y;
}
finalCell = new PolygonRoi(xCoords, yCoords, nPoints, Roi.POLYGON);
cells[n] = finalCell;
}
return cells;
}
private void applyColorModel(ImagePlus imp)
{
// Load the spectrum LUT
WindowManager.setTempCurrentImage(imp);
LutLoader lut = new LutLoader();
lut.run("spectrum");
WindowManager.setTempCurrentImage(null);
// Set zero to black
ImageProcessor ip = imp.getProcessor();
IndexColorModel cm = (IndexColorModel) ip.getColorModel();
byte[] r = new byte[256];
byte[] b = new byte[256];
byte[] g = new byte[256];
cm.getReds(r);
cm.getBlues(b);
cm.getGreens(g);
r[0] = 0;
b[0] = 0;
g[0] = 0;
cm = new IndexColorModel(8, 256, r, g, b);
ip.setColorModel(cm);
}
/**
* Draw the current elliptical cell. Then create an Euclidian Distance Map to use as the weights
* within the provided range. All other points a zeros.
*
* @param pointBounds
* @param params
* @param range
* @return
*/
private FloatProcessor createWeightMap(Rectangle pointBounds, double[] params, double range)
{
EllipticalCell cell = new EllipticalCell();
FloatPolygon ellipse = cell.drawEllipse(params);
ByteProcessor mask = new ByteProcessor(pointBounds.width, pointBounds.height);
mask.setValue(255);
mask.draw(new PolygonRoi(ellipse.xpoints, ellipse.ypoints, ellipse.npoints, PolygonRoi.POLYGON));
EDM edm = new EDM();
FloatProcessor map = edm.makeFloatEDM(mask, (byte) 255, false);
if (map == null)
{
// Preview thread has been interrupted
return null;
}
map.invert();
double max = map.getMax();
double min = Math.max(max - range, 0);
float[] data = (float[]) map.getPixels();
for (int i = 0; i < data.length; i++)
data[i] = (float) ((data[i] > min) ? (data[i] - min) / range : 0f);
// Apply a gamma function for a smoother roll-off
double g = weightingGamma;
if (g > 0)
{
g = 1.0 / g;
for (int i = 0; i < data.length; i++)
data[i] = (float) Math.pow(data[i], g);
}
if (debug)
{
map.resetMinAndMax();
ImagePlus mapImp = displayImage(map, "Current weight map");
mapImp.updateAndDraw();
}
return map;
}
private FloatProcessor applyWeights(FloatProcessor edgeProjection, FloatProcessor weights)
{
FloatProcessor ip = new FloatProcessor(edgeProjection.getWidth(), edgeProjection.getHeight());
float[] e = (float[]) edgeProjection.getPixels();
float[] w = (float[]) weights.getPixels();
float[] p = (float[]) ip.getPixels();
for (int i = 0; i < p.length; i++)
p[i] = e[i] * w[i];
return ip;
}
private ImagePlus displayImage(ImageProcessor maskIp, String title)
{
ImagePlus maskImp = WindowManager.getImage(title);
if (maskImp == null)
{
maskImp = new ImagePlus(title, maskIp);
maskImp.show();
}
else
{
maskImp.setProcessor(maskIp);
//maskImp.updateAndDraw();
}
return maskImp;
}
/**
* Find the bounding rectangle that contains all the pixels that will be required for processing
*
* @param ip
* @return
*/
private Rectangle createBounds(ImageProcessor ip)
{
// Get the range of clicked points
int minx = Integer.MAX_VALUE;
int maxx = 0;
for (int x : xpoints)
{
if (minx > x)
minx = x;
if (maxx < x)
maxx = x;
}
int miny = Integer.MAX_VALUE;
int maxy = 0;
for (int y : ypoints)
{
if (miny > y)
miny = y;
if (maxy < y)
maxy = y;
}
// Add the cell width, tolerance and kernel size to get the total required limits
int extra = (int) Math.ceil(cellRadius + cellRadius * tolerance + kernelWidth / 2);
minx -= extra;
miny -= extra;
maxx += extra;
maxy += extra;
minx = Math.max(0, minx);
miny = Math.max(0, miny);
maxx = Math.min(ip.getWidth() - 1, maxx);
maxy = Math.min(ip.getHeight() - 1, maxy);
return new Rectangle(minx, miny, maxx - minx + 1, maxy - miny + 1);
}
/**
* Build the convolution kernels
*
* @return
*/
private HashMap<Integer, float[]> createKernels()
{
HashMap<Integer, float[]> kernels = new HashMap<Integer, float[]>();
rotationAngles.clear();
// Used to weight the kernel from the distance to the centre
halfWidth = kernelWidth / 2;
maxDistance2 = halfWidth * halfWidth * 2;
// Build a set convolution kernels at 6 degree intervals
int cx = halfWidth, cy = halfWidth;
double degreesToRadians = Math.PI / 180.0;
for (int rotation = 0; rotation < 180; rotation += 12)
{
rotationAngles.add(rotation);
FloatProcessor fp = new FloatProcessor(kernelWidth, kernelWidth);
//createAliasedLines(halfWidth, cx, cy, degreesToRadians, rotation, fp);
// Build curves using the cell size.
double centreX = cx - Math.cos(rotation * degreesToRadians) * cellRadius;
double centreY = cy - Math.sin(rotation * degreesToRadians) * cellRadius;
createDoGKernel(fp, centreX, centreY, cellRadius, 1.0);
// Draw circles to use as a directional edge filter
//drawCircle(fp, centreX, centreY, cellRadius + 1, -1 * value);
//drawCircle(fp, centreX, centreY, cellRadius, 2 * value);
//drawCircle(fp, centreX, centreY, cellRadius - 1, -1 * value);
kernels.put(rotation, (float[]) fp.getPixels());
}
// Copy the kernels for the second half of the circle
for (int rotation : rotationAngles.toArray(new Integer[0]))
{
rotationAngles.add(rotation + 180);
FloatProcessor fp = new FloatProcessor(kernelWidth, kernelWidth, kernels.get(rotation), null);
fp = (FloatProcessor) fp.duplicate();
fp.flipHorizontal();
fp.flipVertical();
kernels.put(rotation + 180, (float[]) fp.getPixels());
}
// Show for debugging
//for (int rotation : rotationAngles)
//{
// FloatProcessor fp = new FloatProcessor(kernelWidth, kernelWidth, kernels.get(rotation), null);
// new ImagePlus("Rotation " + rotation, fp).show();
//}
return kernels;
}
/**
* Create a 1D Difference-of-Gaussians kernel using the distance from the centre minus the cell radius to determine
* the x distance.
*
* @param fp
* @param centreX
* @param centreY
* @param cellRadius
* @param width
*/
private void createDoGKernel(FloatProcessor fp, double centreX, double centreY, int cellRadius, double width)
{
// 1.6:1 is an approximation of the Laplacian of Gaussians
double sigma1 = 1.6 * width;
double sigma2 = width;
if (kernelSmoothing > 0)
{
sigma1 *= kernelSmoothing;
sigma2 *= kernelSmoothing;
}
float[] data = (float[]) fp.getPixels();
for (int y = 0, i = 0; y < fp.getHeight(); y++)
{
for (int x = 0; x < fp.getWidth(); x++, i++)
{
double dx = centreX - x;
double dy = centreY - y;
double x0 = cellRadius - Math.sqrt(dx * dx + dy * dy);
float value = (float) (gaussian(x0, sigma1) - gaussian(x0, sigma2));
if (!darkEdge)
value = -value;
// Weight the amount using the distance from the centre
double d = (x - halfWidth) * (x - halfWidth) + (y - halfWidth) * (y - halfWidth);
value *= 1 - (d / maxDistance2);
data[i] = value;
}
}
}
private double gaussian(double x0, double sigma)
{
return 1.0 / (2 * Math.PI * sigma * sigma) * Math.exp(-0.5 * x0 * x0 / (sigma * sigma));
}
/**
* Old method used to create straight line kernels
*
* @param halfWidth
* @param cx
* @param cy
* @param degreesToRadians
* @param rotation
* @param fp
*/
@SuppressWarnings("unused")
private void createAliasedLines(int halfWidth, int cx, int cy, double degreesToRadians, int rotation,
FloatProcessor fp)
{
// Calculate the direction of the line
double dx = Math.sin(rotation * degreesToRadians);
double dy = Math.cos(rotation * degreesToRadians);
// Normalise so that each increment moves one pixel in either X or Y
double norm = Math.max(Math.abs(dx), Math.abs(dy));
dx /= norm;
dy /= norm;
// Create aliased lines
for (int i = 0; i <= halfWidth; i++)
add(fp, cx + i * dx, cy + i * dy, 1);
for (int i = 1; i <= halfWidth; i++)
add(fp, cx - i * dx, cy - i * dy, 1);
}
/**
* Spread a value of 1 using bilinear weighting around the point
*
* @param fp
* @param x
* @param y
* @param value
*/
private void add(FloatProcessor fp, double x, double y, float value)
{
int x1 = (int) Math.floor(x);
int y1 = (int) Math.floor(y);
// Ignore out-of-bounds
if (x1 < -1 || x1 > fp.getWidth() || y1 < -1 || y1 > fp.getHeight())
return;
// Weight the amount using the distance from the centre
double d = (x - halfWidth) * (x - halfWidth) + (y - halfWidth) * (y - halfWidth);
value *= 1 - (d / maxDistance2);
double dx = x - x1;
double dy = y - y1;
add(fp, x1 + 1, y1 + 1, value * dx * dy);
add(fp, x1 + 1, y1, value * dx * (1 - dy));
add(fp, x1, y1 + 1, value * (1 - dx) * dy);
add(fp, x1, y1, value * (1 - dx) * (1 - dy));
}
private void add(FloatProcessor fp, int x, int y, double weight)
{
float value = fp.getPixelValue(x, y);
fp.putPixelValue(x, y, value + weight);
}
private HashMap<Integer, FloatProcessor> convolveImage(ImageProcessor ip, HashMap<Integer, float[]> kernels)
{
HashMap<Integer, FloatProcessor> convolved = new HashMap<Integer, FloatProcessor>();
// Convolve image with each
for (int rotation : rotationAngles)
{
if (!this.buildMaskOutput)
IJ.showStatus("Convolving " + rotation);
float[] kernel = kernels.get(rotation);
FloatProcessor fp = (ip instanceof FloatProcessor) ? (FloatProcessor) ip.duplicate() : ip.toFloat(0, null);
fp.convolve(kernel, kernelWidth, kernelWidth);
convolved.put(rotation, fp);
if (Thread.currentThread().isInterrupted())
{
convolved = null;
break;
}
}
if (!this.buildMaskOutput)
{
IJ.showProgress(1); // Convolver modifies the progress tracker
IJ.showStatus("");
}
return convolved;
}
/**
* Debugging method to show the results of convolution
*
* @param convolved
*/
@SuppressWarnings("unused")
private void showConvolvedImages(HashMap<Integer, FloatProcessor> convolved)
{
ImageProcessor ip = convolved.get(0);
ImageStack stack = new ImageStack(ip.getWidth(), ip.getHeight());
for (int rotation : rotationAngles)
stack.addSlice("Rotation " + rotation, convolved.get(rotation));
ImagePlus imp = new ImagePlus("Membrane filter", stack);
imp.show();
// Output different projections of the results
ZProjector projector = new ZProjector(imp);
for (int projectionMethod = 0; projectionMethod < ZProjector.METHODS.length; projectionMethod++)
{
IJ.showStatus("Projecting " + ZProjector.METHODS[projectionMethod]);
projector.setMethod(projectionMethod);
projector.doProjection();
ImagePlus projImp = projector.getProjection();
projImp.show();
}
}
/**
* For the given point, compute the angle between pixels and the centre.
*
* @param ip
* @param cx
* @param cy
* @param pointBounds
*/
private FloatProcessor createAngleProcessor(ImageProcessor ip, int cx, int cy, Rectangle pointBounds)
{
// Find the bounds to process
int minx = pointBounds.x;
int miny = pointBounds.y;
int maxx = minx + pointBounds.width;
int maxy = miny + pointBounds.height;
FloatProcessor angle = new FloatProcessor(pointBounds.width, pointBounds.height);
for (int y = miny, index = 0; y < maxy; y++)
{
for (int x = minx; x < maxx; x++, index++)
{
int dx = cx - x;
int dy = cy - y;
float a = (float) (Math.atan2(dy, dx) * 180.0 / Math.PI);
angle.setf(index, a + 180f); // Convert to 0-360 domain
}
}
return angle;
}
private FloatProcessor computeEdgeProjection(HashMap<Integer, FloatProcessor> convolved, int cx, int cy,
Rectangle pointBounds, FloatProcessor angle)
{
float[] a = (float[]) angle.getPixels();
// Do a projection of all membrane filters
float[][] stack = new float[rotationAngles.size()][];
for (int i = 0; i < rotationAngles.size(); i++)
{
int rotation = rotationAngles.get(i);
FloatProcessor fp = (FloatProcessor) convolved.get(rotation).duplicate();
fp.setRoi(pointBounds);
fp = (FloatProcessor) fp.crop();
if (darkEdge)
fp.invert();
float[] p = (float[]) fp.getPixels();
// Do a projection of membrane filters convolved with a filter roughly perpendicular to the edge
for (int j = 0; j < p.length; j++)
{
// Check if angle is close enough
float diff = (a[j] > rotation) ? a[j] - rotation : rotation - a[j];
if (diff >= 180)
diff -= 360;
if (Math.abs(diff) > 30)
p[j] = 0;
}
stack[i] = p;
}
// Perform Z-projection
FloatProcessor ip2 = new FloatProcessor(pointBounds.width, pointBounds.height);
float[] pixels = (float[]) ip2.getPixels();
for (float[] tmp : stack)
{
// Max intensity
for (int i = 0; i < tmp.length; i++)
{
if (pixels[i] < tmp[i])
pixels[i] = tmp[i];
}
}
return ip2;
}
private static final byte CELL = (byte) 255;
private FloatProcessor cropToValues(FloatProcessor fp, Rectangle cropBounds)
{
// Find the bounds
int minx = fp.getWidth();
int miny = fp.getHeight();
int maxx = 0;
int maxy = 0;
float[] data = (float[]) fp.getPixels();
for (int i = 0; i < data.length; i++)
{
if (data[i] != 0)
{
int x = i % fp.getWidth();
int y = i / fp.getWidth();
if (minx > x)
minx = x;
if (maxx < x)
maxx = x;
if (miny > y)
miny = y;
if (maxy < y)
maxy = y;
}
}
// Crop to the bounds
cropBounds.x = minx;
cropBounds.y = miny;
cropBounds.width = maxx - minx + 1;
cropBounds.height = maxy - miny + 1;
fp.setRoi(cropBounds);
FloatProcessor fp2 = (FloatProcessor) fp.crop();
fp.setRoi((Rectangle) null);
return fp2;
}
/**
* Provides methods for drawing an elliptical path composed of two ellipses placed back-to-back.
* <p>
* The ellipse is constructed using two values for the major axis. This allows the function to draw two ellipses
* back-to-back. The input parameters can be converted into the actual first and second major axis values using the
* helper functions.
*/
public class EllipticalCell
{
/**
* Convert the input function parameters for the major axis into the first actual major axis
*
* @param axis1
* @param axis2
* @return
*/
public double getMajor1(double axis1, double axis2)
{
return (2 * axis1 + axis2) / 3;
}
/**
* Convert the input function parameters for the major axis into the second actual major axis
*
* @param axis1
* @param axis2
* @return
*/
public double getMajor2(double axis1, double axis2)
{
return (2 * axis2 + axis1) / 3;
}
/**
* Draws the elliptical cell
*
* @param params
* @return
*/
public FloatPolygon drawEllipse(final double[] params)
{
final double centreX = params[0];
final double centreY = params[1];
final double axis1 = params[2];
final double axis2 = params[3];
final double minor = params[4];
final double phi = params[5];
return drawEllipse(centreX, centreY, axis1, axis2, minor, phi);
}
/**
* Draw an elliptical path of points. The ellipse is constructed as two half ellipses allowing a two different
* major axes lengths, one for each side (i.e. an egg shape).
*
* @param centreX
* @param centreY
* @param axis1
* @param axis2
* @param minor
* @param phi
* The angle from X-axis and the major axis of the ellipse
* @return
*/
public FloatPolygon drawEllipse(final double centreX, final double centreY, final double axis1,
final double axis2, final double minor, final double phi)
{
int nPoints = 90;
double arcAngle = 2 * Math.PI / nPoints;
float[] xPoints = new float[nPoints];
float[] yPoints = new float[nPoints];
int n = 0;
double major1 = getMajor1(axis1, axis2);
double major2 = getMajor2(axis1, axis2);
double angleLimit = Math.PI / 2;
double a1 = major1 * Math.cos(phi);
double a2 = major1 * Math.sin(phi);
double b1 = minor * Math.sin(phi);
double b2 = minor * Math.cos(phi);
// Create points around the start of the first half
for (double angle = -Math.PI; angle < -angleLimit; angle += arcAngle)
{
xPoints[n] = (float) (centreX + a1 * Math.cos(angle) - b1 * Math.sin(angle));
yPoints[n] = (float) (centreY + a2 * Math.cos(angle) + b2 * Math.sin(angle));
n++;
}
// Create points around the second half
a1 = major2 * Math.cos(phi);
a2 = major2 * Math.sin(phi);
for (double angle = -angleLimit; angle < angleLimit; angle += arcAngle)
{
xPoints[n] = (float) (centreX + a1 * Math.cos(angle) - b1 * Math.sin(angle));
yPoints[n] = (float) (centreY + a2 * Math.cos(angle) + b2 * Math.sin(angle));
n++;
}
// Create points around the rest of the first half
a1 = major1 * Math.cos(phi);
a2 = major1 * Math.sin(phi);
for (double angle = angleLimit; angle < Math.PI && n < nPoints; angle += arcAngle)
{
xPoints[n] = (float) (centreX + a1 * Math.cos(angle) - b1 * Math.sin(angle));
yPoints[n] = (float) (centreY + a2 * Math.cos(angle) + b2 * Math.sin(angle));
n++;
}
return new FloatPolygon(xPoints, yPoints, nPoints);
}
}
/**
* Try and find a polygon that traces a path around the detected edges.
*
* @param cx
* @param cy
* @param fp
* @param angle
* @return
*/
private PolygonRoi findPolygonalCell(int cx, int cy, FloatProcessor fp, FloatProcessor angle)
{
Rectangle cropBounds = new Rectangle();
FloatProcessor cropFp = cropToValues(fp, cropBounds);
angle.setRoi(cropBounds);
FloatProcessor cropAngles = (FloatProcessor) angle.crop();
if (debug)
{
ImagePlus fpImp = displayImage(cropFp, "Current edge outline");
fpImp.updateAndDraw();
}
//ImagePlus mapImp = displayImage(cropAngles, "Current angles");
//mapImp.updateAndDraw();
// Divide the region around the centre into segments and
// find the maxima in each segment
float[] a = (float[]) cropAngles.getPixels();
int segmentArcWidth = 5;
float[] max = new float[360 / segmentArcWidth + 1];
int[] index = new int[max.length];
for (int i = 0; i < a.length; i++)
{
int segment = (int) (a[i] / segmentArcWidth);
if (max[segment] < cropFp.getf(i))
{
max[segment] = cropFp.getf(i);
index[segment] = i;
}
}
// Count the number of points.
// Skip consecutive segments to ensure the sampled points are spread out.
int nPoints = 0;
for (int i = 0; i < 360 / segmentArcWidth; i += 2)
if (max[i] > 0)
nPoints++;
// Return the polygon
float[] xPoints = new float[nPoints];
float[] yPoints = new float[nPoints];
nPoints = 0;
for (int i = 0; i < 360 / segmentArcWidth; i += 2)
if (max[i] > 0)
{
xPoints[nPoints] = index[i] % cropBounds.width + cropBounds.x;
yPoints[nPoints] = index[i] / cropBounds.width + cropBounds.y;
nPoints++;
}
// Perform coordinate smoothing
int window = polygonSmoothing;
if (window > 0)
{
float[] newXPoints = new float[nPoints];
float[] newYPoints = new float[nPoints];
window = Math.min(nPoints / 4, window);
for (int i = 0; i < nPoints; i++)
{
double sumx = 0, sumy = 0, n = 0;
for (int j = -window; j <= window; j++)
{
double w = 1; //(window - Math.abs(j) + 1.0) / (window + 1);
int ii = (i + j + nPoints) % nPoints;
sumx += xPoints[ii] * w;
sumy += yPoints[ii] * w;
n += w;
}
newXPoints[i] = (float) (sumx / n);
newYPoints[i] = (float) (sumy / n);
}
xPoints = newXPoints;
yPoints = newYPoints;
}
return new PolygonRoi(xPoints, yPoints, nPoints, PolygonRoi.POLYGON);
}
/**
* Create a byte processor of the specified dimensions and fill the ROI
*
* @param width
* @param height
* @param roi
* @return
*/
private ByteProcessor createFilledCell(int width, int height, PolygonRoi roi)
{
return createFilledCell(width, height, roi, CELL & 0xff);
}
/**
* Create a byte processor of the specified dimensions and fill the ROI
*
* @param width
* @param height
* @param roi
* @param value
* @return
*/
private ByteProcessor createFilledCell(int width, int height, PolygonRoi roi, int value)
{
ByteProcessor bp = new ByteProcessor(width, height);
bp.setColor(value);
bp.fill(roi);
return bp;
}
/**
* Find an ellipse that optimises the fit to the polygon detected edges.
*
* @param roi
* @param params
* @param weightMap
* @param angle
* @return
*/
private double[] fitPolygonalCell(PolygonRoi roi, double[] params, FloatProcessor weightMap, FloatProcessor angle)
{
// Get an estimate of the starting parameters using the current polygon
double[] startPoint = params;
startPoint = estimateStartPoint(roi, weightMap.getWidth(), weightMap.getHeight());
int maxEval = 2000;
final DifferentiableEllipticalFitFunction func = new DifferentiableEllipticalFitFunction(roi, weightMap);
double relativeThreshold = 100 * Precision.EPSILON;
double absoluteThreshold = 100 * Precision.SAFE_MIN;
ConvergenceChecker<PointVectorValuePair> checker = new SimplePointChecker<PointVectorValuePair>(
relativeThreshold, absoluteThreshold);
double initialStepBoundFactor = 10;
double costRelativeTolerance = 1e-10;
double parRelativeTolerance = 1e-10;
double orthoTolerance = 1e-10;
double threshold = Precision.SAFE_MIN;
LevenbergMarquardtOptimizer optimiser = new LevenbergMarquardtOptimizer(initialStepBoundFactor,
costRelativeTolerance, parRelativeTolerance, orthoTolerance, threshold);
try
{
//@formatter:off
LeastSquaresProblem problem = new LeastSquaresBuilder()
.maxEvaluations(Integer.MAX_VALUE)
.maxIterations(maxEval)
.start(startPoint)
.target(func.calculateTarget())
.weight(new DiagonalMatrix(func.calculateWeights()))
.model(func, new MultivariateMatrixFunction() {
public double[][] value(double[] point) throws IllegalArgumentException
{
return func.jacobian(point);
}} )
.checkerPair(checker)
.build();
//@formatter:on
Optimum solution = optimiser.optimize(problem);
if (debug)
IJ.log(String.format("Eval = %d (Iter = %d), RMS = %f", solution.getEvaluations(),
solution.getIterations(), solution.getRMS()));
return solution.getPoint().toArray();
}
catch (Exception e)
{
IJ.log("Failed to find an elliptical solution, defaulting to polygon");
e.printStackTrace();
}
return null;
}
/**
* Estimate the starting ellipse using the eccentricity around the central moments
* <p>
*
* @see Burger & Burge, Digital Image Processing, An Algorithmic Introduction using Java (1st Edition), pp231
*
* @param roi
* @param height
* @param width
* @return
*/
private double[] estimateStartPoint(PolygonRoi roi, int width, int height)
{
ByteProcessor ip = createFilledCell(width, height, roi);
byte[] data = (byte[]) ip.getPixels();
//// Default processing
//double m00 = moment(ip, 0, 0); // region area
//double xCtr = moment(ip, 1, 0) / m00;
//double yCtr = moment(ip, 0, 1) / m00;
//
//double u00 = m00; //centralMoment(ip, 0, 0);
//double u11 = centralMoment(ip, 1, 1);
//double u20 = centralMoment(ip, 2, 0);
//double u02 = centralMoment(ip, 0, 2);
// Speed up processing of the moments
double u00 = 0;
for (byte b : data)
if (b != 0)
u00++;
double xCtr = 0, yCtr = 0;
for (int v = 0, i = 0; v < ip.getHeight(); v++)
{
for (int u = 0; u < ip.getWidth(); u++, i++)
{
if (data[i] != 0)
{
xCtr += u;
yCtr += v;
}
}
}
xCtr /= u00;
yCtr /= u00;
double u11 = 0;
double u20 = 0;
double u02 = 0;
for (int v = 0, i = 0; v < ip.getHeight(); v++)
{
for (int u = 0; u < ip.getWidth(); u++, i++)
{
if (data[i] != 0)
{
double dx = u - xCtr;
double dy = v - yCtr;
u11 += dx * dy;
u20 += dx * dx;
u02 += dy * dy;
}
}
}
// Calculate the ellipsoid
double A = 2 * u11;
double B = u20 - u02;
double angle = 0.5 * Math.atan2(A, B);
double a1 = u20 + u02 + Math.sqrt((u20 - u02) * (u20 - u02) + 4 * u11 * u11);
double a2 = u20 + u02 - Math.sqrt((u20 - u02) * (u20 - u02) + 4 * u11 * u11);
double ra = Math.sqrt(2 * a1 / u00);
double rb = Math.sqrt(2 * a2 / u00);
double[] params = new double[] { xCtr, yCtr, ra, ra, rb, angle };
if (debug)
{
EllipticalCell cell = new EllipticalCell();
FloatPolygon ellipse = cell.drawEllipse(params);
ip = createFilledCell(width, height,
new PolygonRoi(ellipse.xpoints, ellipse.ypoints, ellipse.npoints, PolygonRoi.POLYGON));
ip.setMinAndMax(0, CELL);
displayImage(ip, "Start estimate");
}
return params;
}
static final int BACKGROUND = 0;
public static double moment(ImageProcessor ip, int p, int q)
{
double Mpq = 0.0;
for (int v = 0, i = 0; v < ip.getHeight(); v++)
{
for (int u = 0; u < ip.getWidth(); u++, i++)
{
if (ip.get(i) != BACKGROUND)
{
Mpq += Math.pow(u, p) * Math.pow(v, q);
}
}
}
return Mpq;
}
public static double centralMoment(ImageProcessor ip, int p, int q)
{
double m00 = moment(ip, 0, 0); // region area
double xCtr = moment(ip, 1, 0) / m00;
double yCtr = moment(ip, 0, 1) / m00;
double cMpq = 0.0;
for (int v = 0, i = 0; v < ip.getHeight(); v++)
{
for (int u = 0; u < ip.getWidth(); u++, i++)
{
if (ip.get(i) != BACKGROUND)
{
cMpq += Math.pow(u - xCtr, p) * Math.pow(v - yCtr, q);
}
}
}
return cMpq;
}
/**
* Provides a function to score the ellipse for use in gradient based optimisation methods.
*
* @see http://commons.apache.org/math/userguide/optimization.html
*/
public class DifferentiableEllipticalFitFunction extends EllipticalCell implements MultivariateVectorFunction
{
FloatProcessor weightMap;
int nPoints;
int[] xPoints;
int[] yPoints;
// Debugging variables
int iter = 0;
/**
* @param roi
* The polygon to fit
* @param weightMap
*/
public DifferentiableEllipticalFitFunction(PolygonRoi roi, FloatProcessor weightMap)
{
// These methods try to minimise the difference between a target value and your model value.
// The target value is the polygon outline. The model is currently an elliptical path.
this.weightMap = weightMap;
nPoints = roi.getNCoordinates();
xPoints = Arrays.copyOf(roi.getXCoordinates(), nPoints);
yPoints = Arrays.copyOf(roi.getYCoordinates(), nPoints);
Rectangle bounds = roi.getBounds();
for (int i = 0; i < nPoints; i++)
{
xPoints[i] += bounds.x;
yPoints[i] += bounds.y;
}
}
/**
* @return Each point to be fitted
*/
public double[] calculateTarget()
{
ByteProcessor bp = new ByteProcessor(weightMap.getWidth(), weightMap.getHeight());
for (int i = 0; i < nPoints; i++)
bp.putPixel(xPoints[i], yPoints[i], nPoints);
if (debug)
{
bp.setMinAndMax(1, nPoints);
displayImage(bp, "Ellipse target");
}
// We want to minimise the distance to the polygon points.
// Our values will be the distances so the target can be zeros.
double[] target = new double[nPoints];
return target;
}
/**
* @return The weights for each point to be fitted
*/
public double[] calculateWeights()
{
double[] weights = new double[nPoints];
for (int i = 0; i < weights.length; i++)
{
weights[i] = weightMap.getPixelValue(yPoints[i], yPoints[i]);
//weights[i] = 1;
}
return weights;
}
/*
* (non-Javadoc)
*
* @see org.apache.commons.math3.analysis.MultivariateVectorFunction#value(double[])
*/
public double[] value(double[] point) throws IllegalArgumentException
{
if (debug)
System.out.printf("%f,%f %f,%f,%f %f\n", point[0], point[1], point[2], point[3], point[4],
point[5] * 180.0 / Math.PI);
return getValue(point);
}
private double[] getValue(double[] point) throws IllegalArgumentException
{
double[] values = new double[nPoints];
double cx = point[0];
double cy = point[1];
double axis1 = point[2];
double axis2 = point[3];
double minor = point[4];
double phi = point[5];
// Link the axes parameters so that each one will have a gradient for all points.
double major1 = getMajor1(axis1, axis2);
double major2 = getMajor2(axis1, axis2);
ByteProcessor bp = new ByteProcessor(weightMap.getWidth(), weightMap.getHeight());
// Calculate the distance of the ellipse to each point
for (int i = 0; i < values.length; i++)
{
// Get the angle between the point and the centre
double absAngle = Math.atan2(yPoints[i] - cy, xPoints[i] - cx);
// Adjust for the ellipse rotation
double angle = absAngle - phi;
// Check domain
if (angle < -Math.PI)
angle += 2 * Math.PI;
if (angle > Math.PI)
angle -= 2 * Math.PI;
double a1, a2;
double b1 = minor * Math.sin(phi);
double b2 = minor * Math.cos(phi);
// Create ellipse point. The shape is two ellipse halves so use the correct part.
if (angle < -Math.PI / 2 || angle > Math.PI / 2)
{
a1 = major1 * Math.cos(phi);
a2 = major1 * Math.sin(phi);
}
else
{
a1 = major2 * Math.cos(phi);
a2 = major2 * Math.sin(phi);
}
double x = cx + a1 * Math.cos(angle) - b1 * Math.sin(angle);
double y = cy + a2 * Math.cos(angle) + b2 * Math.sin(angle);
bp.putPixel((int) x, (int) y, nPoints);
// Get the distance
double dx = x - xPoints[i];
double dy = y - yPoints[i];
values[i] = Math.sqrt(dx * dx + dy * dy);
// Check if it is inside or outside the ellipse
dx = cx - xPoints[i];
dy = cy - yPoints[i];
double pointToCentre = dx * dx + dy * dy;
dx = cx - x;
dy = cy - y;
double ellipseToCentre = dx * dx + dy * dy;
if (pointToCentre < ellipseToCentre)
values[i] = -values[i];
}
if (debug)
{
bp.setMinAndMax(1, nPoints);
displayImage(bp, "Ellipse points " + iter);
}
return values;
}
private double[][] jacobian(double[] variables)
{
double[][] jacobian = new double[nPoints][variables.length];
// Compute numerical differentiation
// Param order:
// centreX, centreY, major1, major2, minor, phi
double[] delta = { 1, 1, 2, 2, 2, 5 * Math.PI / 180.0 };
for (int i = 0; i < variables.length - 1; i++)
{
double[] upper = getValue(updateVariables(variables, i, delta, 1));
double[] lower = getValue(updateVariables(variables, i, delta, -1));
for (int j = 0; j < jacobian.length; ++j)
{
jacobian[j][i] = (upper[j] - lower[j]) / (2 * delta[i]);
}
}
// Only compute the angle gradient if the ellipse is not a circle
if (variables[2] != variables[3] || variables[2] != variables[4])
{
int i = variables.length - 1;
double[] upper = getValue(updateVariables(variables, i, delta, 1));
double[] lower = getValue(updateVariables(variables, i, delta, -1));
for (int j = 0; j < jacobian.length; ++j)
{
jacobian[j][i] = (upper[j] - lower[j]) / (2 * delta[i]);
}
}
return jacobian;
}
private double[] updateVariables(double[] variables, int i, double[] delta, int sign)
{
double[] newVariables = Arrays.copyOf(variables, variables.length);
newVariables[i] += delta[i] * sign;
return newVariables;
}
/*
* (non-Javadoc)
*
* @see org.apache.commons.math3.analysis.DifferentiableMultivariateVectorFunction#jacobian()
*/
public MultivariateMatrixFunction jacobian()
{
return new MultivariateMatrixFunction()
{
public double[][] value(double[] point)
{
return jacobian(point);
}
};
}
}
private void dilate(ByteProcessor bp)
{
byte[] data = (byte[]) bp.getPixels();
byte[] newData = new byte[data.length];
initialise(bp);
for (int index = 0; index < data.length; index++)
{
if (data[index] != 0)
{
newData[index] = CELL;
int x = index % maxx;
int y = index / maxx;
boolean isInnerXY = (y != 0 && y != ylimit) && (x != 0 && x != xlimit);
// Use 4-connected cells
if (isInnerXY)
{
for (int d = 0; d < 8; d += 2)
{
newData[index + offset[d]] = CELL;
}
}
else
{
for (int d = 0; d < 8; d += 2)
{
if (isWithinXY(x, y, d))
newData[index + offset[d]] = CELL;
}
}
}
}
bp.setPixels(newData);
}
private PolygonRoi traceOutline(ByteProcessor bp, double cx, double cy)
{
byte[] data = (byte[]) bp.getPixels();
initialise(bp);
// Find first pixel
int startIndex = 0;
while (data[startIndex] == 0 && startIndex < data.length)
startIndex++;
if (startIndex == data.length)
return null;
ArrayList<Point> coords = new ArrayList<Point>(100);
addPoint(coords, startIndex);
// Set start direction for search
int searchDirection = 7;
int index = startIndex;
// Safety limit - The outline shouldn't be greater than the image perimeter
int limit = (bp.getWidth() + bp.getHeight()) * 2 - 2;
while (limit-- > 0)
{
int nextDirection = findNext(data, index, searchDirection);
if (nextDirection >= 0)
{
index += offset[nextDirection];
if (index == startIndex)
break; // End of the outline
addPoint(coords, index);
searchDirection = (nextDirection + 6) % 8;
}
else
break; // Single point with no move direction
}
if (limit <= 0)
return null;
// Return the outline
int nPoints = 0;
int[] xPoints = new int[coords.size()];
int[] yPoints = new int[coords.size()];
for (Point p : coords)
{
xPoints[nPoints] = p.x;
yPoints[nPoints] = p.y;
nPoints++;
}
return new PolygonRoi(xPoints, yPoints, nPoints, PolygonRoi.POLYGON);
}
private void addPoint(ArrayList<Point> coords, int index)
{
Point p = new Point(index % maxx, index / maxx);
//System.out.printf("%d, %d\n", p.x, p.y);
coords.add(p);
}
private int findNext(byte[] data, int index, int direction)
{
int x = index % maxx;
int y = index / maxx;
boolean isInnerXY = (y != 0 && y != ylimit) && (x != 0 && x != xlimit);
// Process the neighbours
for (int d = 0; d < 7; d++)
{
if (isInnerXY || isWithinXY(x, y, direction))
{
int index2 = index + offset[direction];
// Check if foreground
if (data[index2] != 0)
{
return direction;
}
}
direction = (direction + 1) % 8;
}
return -1;
}
private int maxx, maxy, xlimit, ylimit;
private int[] offset;
private final int[] DIR_X_OFFSET = new int[] { 0, 1, 1, 1, 0, -1, -1, -1 };
private final int[] DIR_Y_OFFSET = new int[] { -1, -1, 0, 1, 1, 1, 0, -1 };
/**
* Initialises the global width and height variables. Creates the direction offset tables.
*/
private void initialise(ImageProcessor ip)
{
maxx = ip.getWidth();
maxy = ip.getHeight();
xlimit = maxx - 1;
ylimit = maxy - 1;
// Create the offset table (for single array 3D neighbour comparisons)
offset = new int[DIR_X_OFFSET.length];
for (int d = offset.length; d-- > 0;)
{
offset[d] = DIR_X_OFFSET[d] + maxx * DIR_Y_OFFSET[d];
}
}
/**
* returns whether the neighbour in a given direction is within the image. NOTE: it is assumed that the pixel x,y
* itself is within the image! Uses class variables xlimit, ylimit: (dimensions of the image)-1
*
* @param x
* x-coordinate of the pixel that has a neighbour in the given direction
* @param y
* y-coordinate of the pixel that has a neighbour in the given direction
* @param direction
* the direction from the pixel towards the neighbour
* @return true if the neighbour is within the image (provided that x, y is within)
*/
private boolean isWithinXY(int x, int y, int direction)
{
switch (direction)
{
case 0:
return (y > 0);
case 1:
return (y > 0 && x < xlimit);
case 2:
return (x < xlimit);
case 3:
return (y < ylimit && x < xlimit);
case 4:
return (y < ylimit);
case 5:
return (y < ylimit && x > 0);
case 6:
return (x > 0);
case 7:
return (y > 0 && x > 0);
case 8:
return true;
}
return false;
}
}