From 5ae8d69d31a67b74c41a80c4a2e237bc8598ad06 Mon Sep 17 00:00:00 2001 From: dlegland Date: Sun, 22 Jul 2018 20:05:34 +0200 Subject: [PATCH] add max feret diameter computation classes --- src/main/java/inra/ijpb/geometry/Point3D.java | 8 + .../java/inra/ijpb/geometry/PointPair3D.java | 34 +++ .../measure/region2d/MaxFeretDiameter.java | 2 +- .../measure/region3d/MaxFeretDiameter3D.java | 205 ++++++++++++++++++ .../measure/region3d/RegionBoundaries3D.java | 157 ++++++++++++++ .../inra/ijpb/measure/region3d/AllTests.java | 10 +- .../region3d/MaxFeretDiameter3DTest.java | 42 ++++ 7 files changed, 453 insertions(+), 5 deletions(-) create mode 100644 src/main/java/inra/ijpb/geometry/PointPair3D.java create mode 100644 src/main/java/inra/ijpb/measure/region3d/MaxFeretDiameter3D.java create mode 100644 src/main/java/inra/ijpb/measure/region3d/RegionBoundaries3D.java create mode 100644 src/test/java/inra/ijpb/measure/region3d/MaxFeretDiameter3DTest.java diff --git a/src/main/java/inra/ijpb/geometry/Point3D.java b/src/main/java/inra/ijpb/geometry/Point3D.java index 8f94f8f1..b60f31d5 100644 --- a/src/main/java/inra/ijpb/geometry/Point3D.java +++ b/src/main/java/inra/ijpb/geometry/Point3D.java @@ -29,6 +29,14 @@ public Point3D(double x, double y, double z) this.z = z; } + public double distance(Point3D point) + { + double dx = (point.x - x); + double dy = (point.y - y); + double dz = (point.z - z); + return Math.hypot(Math.hypot(dx, dy), dz); + } + public double getX() { return x; diff --git a/src/main/java/inra/ijpb/geometry/PointPair3D.java b/src/main/java/inra/ijpb/geometry/PointPair3D.java new file mode 100644 index 00000000..3774b641 --- /dev/null +++ b/src/main/java/inra/ijpb/geometry/PointPair3D.java @@ -0,0 +1,34 @@ +/** + * + */ +package inra.ijpb.geometry; + + +/** + * A pair of points in the 3D space, useful for representing result of Max Feret + * Diameter computation or similar problems. + * + * Simply contains the reference to each extremity. + * + * @author dlegland + * + */ +public class PointPair3D +{ + public final Point3D p1; + public final Point3D p2; + + public PointPair3D(Point3D p1, Point3D p2) + { + this.p1 = p1; + this.p2 = p2; + } + + public double diameter() + { + double dx = p2.getX() - p1.getX(); + double dy = p2.getY() - p1.getY(); + double dz = p2.getZ() - p1.getZ(); + return Math.hypot(Math.hypot(dx, dy), dz); + } +} diff --git a/src/main/java/inra/ijpb/measure/region2d/MaxFeretDiameter.java b/src/main/java/inra/ijpb/measure/region2d/MaxFeretDiameter.java index d7b2fc5c..71a4c0f9 100644 --- a/src/main/java/inra/ijpb/measure/region2d/MaxFeretDiameter.java +++ b/src/main/java/inra/ijpb/measure/region2d/MaxFeretDiameter.java @@ -130,7 +130,7 @@ public PointPair2D[] analyzeRegions(ImageProcessor image, int[] labels, Calibrat // Extract spatial calibration of image double sx = 1, sy = 1; - double ox = 0, oy = 1; + double ox = 0, oy = 0; if (calib != null) { sx = calib.pixelWidth; diff --git a/src/main/java/inra/ijpb/measure/region3d/MaxFeretDiameter3D.java b/src/main/java/inra/ijpb/measure/region3d/MaxFeretDiameter3D.java new file mode 100644 index 00000000..550ba366 --- /dev/null +++ b/src/main/java/inra/ijpb/measure/region3d/MaxFeretDiameter3D.java @@ -0,0 +1,205 @@ +/** + * + */ +package inra.ijpb.measure.region3d; + +import java.util.ArrayList; +import java.util.Map; + +import ij.measure.Calibration; +import ij.measure.ResultsTable; +import ij.ImageStack; +import inra.ijpb.geometry.Point3D; +import inra.ijpb.geometry.PointPair3D; + +/** + * Computes maximum Feret Diameter for each region of a 3D binary or label + * image. + * + * @author dlegland + * + */ +public class MaxFeretDiameter3D extends RegionAnalyzer3D +{ + // ================================================== + // Static methods + + public final static PointPair3D[] maxFeretDiameters(ImageStack image, int[] labels, Calibration calib) + { + return new MaxFeretDiameter3D().analyzeRegions(image, labels, calib); + } + + /** + * Computes Maximum Feret diameter of a set of points. + * + * Note: it is often a good idea to compute convex hull before computing + * Feret diameter. + * + * @param points + * a collection of planar points + * @return the maximum Feret diameter of the point set + */ + public final static PointPair3D maxFeretDiameter(ArrayList points) + { + double distMax = Double.NEGATIVE_INFINITY; + PointPair3D maxDiam = null; + + for (Point3D p1 : points) + { + for (Point3D p2 : points) + { + double dist = p1.distance(p2); + if (dist > distMax) + { + maxDiam = new PointPair3D(p1, p2); + distMax = dist; + } + } + } + + return maxDiam; + } + + + // ================================================== + // Constructor + + /** + * Default constructor + */ + public MaxFeretDiameter3D() + { + } + + // ================================================== + // Implementation of RegionAnalyzer interface + + /** + * Converts the result of maximum Feret diameters computation to a + * ResultsTable that can be displayed within ImageJ. + * + * @param maxDiamsMap + * the map of PointPair3D for each label within a label image + * @return a ResultsTable instance + */ + public ResultsTable createTable(Map maxDiamsMap) + { + // Create data table + ResultsTable table = new ResultsTable(); + + // compute ellipse parameters for each region + for (int label : maxDiamsMap.keySet()) + { + table.incrementCounter(); + table.addLabel(Integer.toString(label)); + + // add coordinates of origin pixel (IJ coordinate system) + PointPair3D maxDiam = maxDiamsMap.get(label); + table.addValue("Diameter", maxDiam.diameter()); + table.addValue("P1.x", maxDiam.p1.getX()); + table.addValue("P1.y", maxDiam.p1.getY()); + table.addValue("P1.z", maxDiam.p1.getZ()); + table.addValue("P2.x", maxDiam.p2.getX()); + table.addValue("p2.y", maxDiam.p2.getY()); + table.addValue("p2.z", maxDiam.p2.getZ()); + } + + // return the created array + return table; + } + + /** + * Computes maximum Feret Diameter for each label of the input label image. + * + * Computes diameter between corners of image pixels, so the result is + * always greater than or equal to one. + * + * @param image + * a label image (8, 16 or 32 bits) + * @param labels + * the set of labels within the image + * @param calib + * the spatial calibration of the image + * @return an array of PointPair3D representing the coordinates of extreme + * points, in calibrated coordinates. + */ + public PointPair3D[] analyzeRegions(ImageStack image, int[] labels, Calibration calib) + { + // Check validity of parameters + if (image == null) + return null; + + // Extract spatial calibration of image + double sx = 1, sy = 1, sz = 1; + double ox = 0, oy = 0, oz = 0; + if (calib != null) + { + sx = calib.pixelWidth; + sy = calib.pixelHeight; + sz = calib.pixelDepth; + ox = calib.xOrigin; + oy = calib.yOrigin; + oz = calib.zOrigin; + } + + int nLabels = labels.length; + + // For each label, create a list of corner points + fireStatusChanged(this, "Find Label Corner Points"); + ArrayList[] labelCornerPointsArray = RegionBoundaries3D.regionsCornersArray(image, labels); + + // Compute the oriented box of each set of corner points + PointPair3D[] labelMaxDiams = new PointPair3D[nLabels]; + fireStatusChanged(this, "Compute feret Diameters"); + for (int i = 0; i < nLabels; i++) + { + this.fireProgressChanged(this, i, nLabels); + + ArrayList hull = labelCornerPointsArray[i]; + // calibrate coordinates of hull vertices + for (int iv = 0; iv < hull.size(); iv++) + { + Point3D vertex = hull.get(iv); + vertex = new Point3D(vertex.getX() * sx + ox, vertex.getY() * sy + oy, vertex.getZ() * sz + oz); + hull.set(iv, vertex); + } + + // compute Feret diameter of calibrated hull + labelMaxDiams[i] = maxFeretDiameter(hull); + } + + fireProgressChanged(this, 1, 1); + fireStatusChanged(this, ""); + return labelMaxDiams; + } + + +// /** +// * Computes Maximum Feret diameter from a single region in a binary image. +// * +// * Computes diameter between corners of image pixels, so the result is +// * always greater than or equal to one. +// * +// * @param image +// * a binary image representing the particle. +// * @param calib +// * the spatial calibration +// * @return the maximum Feret diameter of the binary region +// */ +// public PointPair3D analyzeBinary(ImageStack image, double[] calib) +// { +// ArrayList points = RegionBoundaries.binaryParticleCorners(image); +// ArrayList convHull = Polygons2D.convexHull_jarvis(points); +// +// // calibrate coordinates of convex hull vertices +// for (int i = 0; i < convHull.size(); i++) +// { +// Point2D vertex = convHull.get(i); +// vertex = new Point2D.Double(vertex.getX() * calib[0], vertex.getY() * calib[1]); +// convHull.set(i, vertex); +// } +// +// // compute Feret diameter of calibrated vertices +// return maxFeretDiameter(convHull); +// } +} diff --git a/src/main/java/inra/ijpb/measure/region3d/RegionBoundaries3D.java b/src/main/java/inra/ijpb/measure/region3d/RegionBoundaries3D.java new file mode 100644 index 00000000..e65a7041 --- /dev/null +++ b/src/main/java/inra/ijpb/measure/region3d/RegionBoundaries3D.java @@ -0,0 +1,157 @@ +/** + * + */ +package inra.ijpb.measure.region3d; + +import ij.ImageStack; +import inra.ijpb.geometry.Point3D; + +import java.util.ArrayList; +import java.util.Map; +import java.util.TreeMap; + +/** + * Utility functions for computing position of boundary points/corners of + * regions within binary or label images. + * + * @author dlegland + * + */ +public class RegionBoundaries3D +{ + /** + * private constructor to prevent instantiations. + */ + private RegionBoundaries3D() + { + } + + /** + * Returns a set of points located at the corners of a binary particle. + * Point coordinates are integer (ImageJ locates pixels in a [0 1]^d area. + * + * @param image + * a binary image representing the particle + * @param labels + * the list of labels to process + * @return a list of points that can be used for convex hull computation + */ + public final static Map> regionsCorners(ImageStack image, int[] labels) + { + int sizeX = image.getWidth(); + int sizeY = image.getHeight(); + int sizeZ = image.getSize(); + + // For each label, create an empty list of corner points + Map> labelCornerPoints = new TreeMap>(); + for (int label : labels) + { + labelCornerPoints.put(label, new ArrayList()); + } + + // for each row, add corner point for first and last pixel of each run-length + for (int z = 0; z < sizeZ; z++) + { + for (int y = 0; y < sizeY; y++) + { + // start from background + int currentLabel = 0; + + // Identify transition inside and outside the each label + for (int x = 0; x < sizeX; x++) + { + int pixel = (int) image.getVoxel(x, y, z); + if (pixel != currentLabel) + { + // creates the four corners corresponding to the + // transition between previous and current voxels + ArrayList newPoints = new ArrayList(4); + newPoints.add(new Point3D(x, y, z)); + newPoints.add(new Point3D(x, y + 1, z)); + newPoints.add(new Point3D(x, y, z + 1)); + newPoints.add(new Point3D(x, y + 1, z + 1)); + + // if leave a region, add a new corner points for the end of the region + if (currentLabel > 0) + { + ArrayList corners = labelCornerPoints.get(currentLabel); + for (Point3D p : newPoints) + { + if (!corners.contains(p)) + { + corners.add(p); + } + } + } + + // transition into a new region + if (pixel > 0) + { + ArrayList corners = labelCornerPoints.get(pixel); + for (Point3D p : newPoints) + { + if (!corners.contains(p)) + { + corners.add(p); + } + } + } + } + currentLabel = pixel; + } + + // if particle touches right border, add another point + if (currentLabel > 0) + { + // creates the four corners corresponding to the + // transition between current label and background + ArrayList newPoints = new ArrayList(4); + newPoints.add(new Point3D(sizeX, y, z)); + newPoints.add(new Point3D(sizeX, y + 1, z)); + newPoints.add(new Point3D(sizeX, y, z + 1)); + newPoints.add(new Point3D(sizeX, y + 1, z + 1)); + + ArrayList corners = labelCornerPoints.get(currentLabel); + for (Point3D p : newPoints) + { + if (!corners.contains(p)) + { + corners.add(p); + } + } + } + } + } + + return labelCornerPoints; + } + + /** + * Returns a set of points located at the corners of a binary particle. + * Point coordinates are integer (ImageJ locates pixels in a [0 1]^d area. + * + * @param image + * a binary image representing the particle + * @param labels + * the list of labels to process + * @return for each label, an array of points + */ + public final static ArrayList[] regionsCornersArray(ImageStack image, int[] labels) + { + // Compute corner points for each label + Map> cornerPointsMap = regionsCorners(image, labels); + + // allocate array + int nLabels = labels.length; + @SuppressWarnings("unchecked") + ArrayList[] labelCornerPoints = (ArrayList[]) new ArrayList[nLabels]; + + // convert map to array + for (int i = 0; i < nLabels; i++) + { + labelCornerPoints[i] = cornerPointsMap.get(labels[i]); + } + + return labelCornerPoints; + } +} diff --git a/src/test/java/inra/ijpb/measure/region3d/AllTests.java b/src/test/java/inra/ijpb/measure/region3d/AllTests.java index 32ee0b8e..ae0c2f3e 100644 --- a/src/test/java/inra/ijpb/measure/region3d/AllTests.java +++ b/src/test/java/inra/ijpb/measure/region3d/AllTests.java @@ -28,8 +28,10 @@ @RunWith(Suite.class) @Suite.SuiteClasses({ // generic classes - InertiaEllipsoidTest.class, - }) -public class AllTests { - //nothing + InertiaEllipsoidTest.class, + MaxFeretDiameter3DTest.class +}) +public class AllTests +{ + // nothing } diff --git a/src/test/java/inra/ijpb/measure/region3d/MaxFeretDiameter3DTest.java b/src/test/java/inra/ijpb/measure/region3d/MaxFeretDiameter3DTest.java new file mode 100644 index 00000000..412518e3 --- /dev/null +++ b/src/test/java/inra/ijpb/measure/region3d/MaxFeretDiameter3DTest.java @@ -0,0 +1,42 @@ +/** + * + */ +package inra.ijpb.measure.region3d; + +import static org.junit.Assert.assertEquals; +import ij.ImageStack; +import ij.measure.Calibration; +import inra.ijpb.geometry.PointPair3D; + +import org.junit.Test; + +/** + * @author dlegland + * + */ +public class MaxFeretDiameter3DTest +{ + + /** + * Test method for {@link inra.ijpb.measure.region3d.MaxFeretDiameter3D#analyzeRegions(ij.ImageStack, int[], ij.measure.Calibration)}. + */ + @Test + public void testAnalyzeRegionsImageStackIntArrayCalibration() + { + ImageStack stack = ImageStack.create(10, 10, 10, 8); + for (int i = 1; i < 9; i++) + { + stack.setVoxel(i, i, i, 255); + } + int[] labels = new int[]{255}; + Calibration calib = new Calibration(); + + MaxFeretDiameter3D algo = new MaxFeretDiameter3D(); + PointPair3D[] result = algo.analyzeRegions(stack, labels, calib); + + assertEquals(1, result.length); + PointPair3D pair = result[0]; + assertEquals(8*Math.sqrt(3), pair.diameter(), .1); + } + +}