diff --git a/src/main/java/net/imagej/ops/geom/GeomNamespace.java b/src/main/java/net/imagej/ops/geom/GeomNamespace.java index 6e0d449f0..029d092e8 100644 --- a/src/main/java/net/imagej/ops/geom/GeomNamespace.java +++ b/src/main/java/net/imagej/ops/geom/GeomNamespace.java @@ -278,13 +278,14 @@ public List convexHull(final Mesh in) { net.imagej.ops.Ops.Geometric.ConvexHull.class, in); return result; } - + + @Deprecated @OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class) public RandomAccessibleInterval voxelization(final Mesh in, final int width, final int height, final int depth ) { final RandomAccessibleInterval result = (RandomAccessibleInterval) ops().run( Voxelization.class, in, width, height, depth ); - return result; - } + return result; + } @OpMethod(op = net.imagej.ops.geom.geom2d.DefaultConvexityPolygon.class) public DoubleType convexity(final Polygon2D in) { @@ -675,27 +676,46 @@ public double[] vertexInterpolator(final int[] p1, final int[] p2, p1Value, p2Value); return result; } - - @OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class) + + @OpMethod(ops = { + net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class, + net.imagej.ops.geom.geom3d.EuclideanDistanceVoxelization3D.class + }) public RandomAccessibleInterval voxelization(final Mesh in) { @SuppressWarnings("unchecked") final RandomAccessibleInterval result = - (RandomAccessibleInterval) ops().run(net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class, in); + (RandomAccessibleInterval) ops().run(net.imagej.ops.Ops.Geometric.Voxelization.class, in); return result; } + @Deprecated @OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class) public RandomAccessibleInterval voxelization(final Mesh in, final int width) { @SuppressWarnings("unchecked") final RandomAccessibleInterval result = - (RandomAccessibleInterval) ops().run(net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class, in, width); + (RandomAccessibleInterval) ops().run(net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class, in, width); return result; } + @Deprecated @OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class) public RandomAccessibleInterval voxelization(final Mesh in, final int width, final int height) { final RandomAccessibleInterval result = - (RandomAccessibleInterval) ops().run(net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class, in, width, height); + (RandomAccessibleInterval) ops().run(net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class, in, width, height); + return result; + } + + @OpMethod(op = net.imagej.ops.geom.geom3d.EuclideanDistanceVoxelization3D.class) + public RandomAccessibleInterval voxelization(final RandomAccessibleInterval out, final Mesh in) { + final RandomAccessibleInterval result = + (RandomAccessibleInterval) ops().run(net.imagej.ops.geom.geom3d.EuclideanDistanceVoxelization3D.class, out, in); + return result; + } + + @OpMethod(op = net.imagej.ops.geom.geom3d.EuclideanDistanceVoxelization3D.class) + public RandomAccessibleInterval voxelization(final RandomAccessibleInterval out, final Mesh in, final double wallThickness) { + final RandomAccessibleInterval result = + (RandomAccessibleInterval) ops().run(net.imagej.ops.geom.geom3d.EuclideanDistanceVoxelization3D.class, out, in, wallThickness); return result; } diff --git a/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java b/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java index 8837875db..73fd13e8b 100644 --- a/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java +++ b/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java @@ -58,7 +58,9 @@ *

* * @author Kyle Harrington (University of Idaho) + * @deprecated Use {@link EuclideanDistanceVoxelization3D} instead. */ +@Deprecated @Plugin(type = Ops.Geometric.Voxelization.class) public class DefaultVoxelization3D extends AbstractUnaryFunctionOp> implements Ops.Geometric.Voxelization { diff --git a/src/main/java/net/imagej/ops/geom/geom3d/EuclideanDistanceVoxelization3D.java b/src/main/java/net/imagej/ops/geom/geom3d/EuclideanDistanceVoxelization3D.java new file mode 100644 index 000000000..7a9409dfa --- /dev/null +++ b/src/main/java/net/imagej/ops/geom/geom3d/EuclideanDistanceVoxelization3D.java @@ -0,0 +1,251 @@ +/* + * #%L + * ImageJ2 software for multidimensional image processing and analysis. + * %% + * Copyright (C) 2014 - 2024 ImageJ2 developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ +package net.imagej.ops.geom.geom3d; + +import net.imagej.mesh.Mesh; +import net.imagej.mesh.Meshes; +import net.imagej.mesh.Triangle; +import net.imagej.ops.OpService; +import net.imagej.ops.Ops; +import net.imagej.ops.special.hybrid.AbstractUnaryHybridCF; +import net.imglib2.FinalInterval; +import net.imglib2.Interval; +import net.imglib2.RandomAccess; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.iterator.LocalizingIntervalIterator; +import net.imglib2.type.logic.BitType; +import net.imglib2.util.Intervals; +import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; +import org.scijava.ItemIO; +import org.scijava.Priority; +import org.scijava.plugin.Parameter; +import org.scijava.plugin.Plugin; + +/** + *

+ * This is a voxelizer that produces a binary image with values set to true along the surface of the mesh. + * Calculates the Euclidean distance between each pixel in the output image and the nearest point in the mesh, + * and sets the pixel to true if it is within 'wallThickness'/2 pixel units of the mesh. + * Calling this op with the default wallThickness of 1.0, and subsequently performing a + * {@link net.imagej.ops.commands.morphology.FillHoles} command using + * {@link net.imglib2.algorithm.neighborhood.DiamondShape} of size 1 is a functional inverse of + * {@link DefaultMarchingCubes} for any object in the original data that is successfully converted + * to a mesh via marching cubes. + *

+ * + * @see DefaultVoxelization3D for the original Voxelization algorithm. This algorithm was deprecated as it generated results + * that could often differ greatly from the surface of the input mesh, particularly with larger triangular facets. + * @author Andrew McCall (University at Buffalo) + */ +@Plugin(type = Ops.Geometric.Voxelization.class, priority = Priority.HIGH) +public class EuclideanDistanceVoxelization3D> extends AbstractUnaryHybridCF + implements Ops.Geometric.Voxelization { + + @Parameter + private OpService ops; + + @Parameter(type = ItemIO.INPUT, required = false, description = "Average pixel thickness of the resulting voxelized" + + " image surface. Any pixel within 'wallThickness'/2 of the mesh will be set to true. Default value is 1.") + private double wallThickness = 1.0; + + @Override + public void compute(final Mesh input, final O output){ + RandomAccess ra = output.randomAccess(); + + input.triangles().forEach((Triangle t) -> { + Vector3D[] vector3DTriangle = triangleToVector3DTriangle(t); + Vector3D normal = getNormalizedNormal(vector3DTriangle); + //check if triangle is degenerate + if(normal.getX() == Double.NEGATIVE_INFINITY) + return; + Interval triangleBox = boundingBox(vector3DTriangle); + LocalizingIntervalIterator iterator = new LocalizingIntervalIterator(triangleBox); + while (iterator.hasNext()) { + //For each point within triangle bounding box, set to true based on Euclidean distance of point to triangle surface + iterator.fwd(); + if(Intervals.contains(output, iterator.positionAsPoint())) { + if (pointToTriangleDist(new Vector3D(iterator.getDoublePosition(0), iterator.getDoublePosition(1), iterator.getDoublePosition(2)), vector3DTriangle, normal) <= wallThickness/2) { + ra.setPositionAndGet(iterator.positionAsPoint()).set(true); + } + } + } + }); + } + + @Override + public O createOutput(Mesh input){ + + float[] bounds = Meshes.boundingBox(input); + long[] min = new long[3]; + long[] max = new long[3]; + for (int i = 0; i < 3; i++) { + min[i] = (long)Math.floor(bounds[i]-wallThickness/2); + max[i] = (long)Math.ceil(bounds[i+3]+wallThickness/2); + } + return (O) ops.create().img(new FinalInterval(min,max), new BitType()); + } + + private Vector3D[] triangleToVector3DTriangle(Triangle t){ + Vector3D[] o = new Vector3D[3]; + + o[0] = new Vector3D(t.v0x(), t.v0y(), t.v0z()); + o[1] = new Vector3D(t.v1x(), t.v1y(), t.v1z()); + o[2] = new Vector3D(t.v2x(), t.v2y(), t.v2z()); + + return o; + } + + private Vector3D getNormalizedNormal(Vector3D[] t){ + Vector3D ab = t[1].subtract(t[0]); + Vector3D ac = t[2].subtract(t[0]); + + // Find the normal to the plane: n = ab x ac + Vector3D n = ab.crossProduct(ac); + + // Normalize normal vector + try{ + n = n.normalize(); + } + catch(Exception e){ + return new Vector3D(Double.NEGATIVE_INFINITY,Double.NEGATIVE_INFINITY,Double.NEGATIVE_INFINITY); // Triangle is degenerate + } + return n; + } + + private Interval boundingBox(Vector3D[] t){ + long [] min = new long[3]; + long [] max = new long[3]; + + min[0] = (long) Math.floor(Math.min(t[0].getX(), Math.min(t[1].getX(), t[2].getX()))); + min[1] = (long) Math.floor(Math.min(t[0].getY(), Math.min(t[1].getY(), t[2].getY()))); + min[2] = (long) Math.floor(Math.min(t[0].getZ(), Math.min(t[1].getZ(), t[2].getZ()))); + + max[0] = (long) Math.ceil(Math.max(t[0].getX(), Math.max(t[1].getX(), t[2].getX()))); + max[1] = (long) Math.ceil(Math.max(t[0].getY(), Math.max(t[1].getY(), t[2].getY()))); + max[2] = (long) Math.ceil(Math.max(t[0].getZ(), Math.max(t[1].getZ(), t[2].getZ()))); + return new FinalInterval(min, max); + } + + private double pointToTriangleDist(Vector3D p, Vector3D[] t, Vector3D n){ + Vector3D tPoint = nearestPointInTriangle3D(p, t, n); + return p.distance(tPoint); + } + + private Vector3D nearestPointInTriangle3D(Vector3D p, Vector3D[] t, Vector3D n) { + /* + Need to project point 'p' onto the plane of triangle 't' as first step. This allows + use of the barycentric coordinate system to locate the nearest point in 't.' + */ + Vector3D ab = t[1].subtract(t[0]); + Vector3D ac = t[2].subtract(t[0]); + + //region Use normal to obtain projection of point p onto plane of triangle + // Project point p onto the plane spanned by a->b and a->c. + double dist = p.dotProduct(n) - t[0].dotProduct(n); + Vector3D projection = p.add(n.scalarMultiply(-dist)); + //endregion + + //Define projection of 'p' onto triangle plane as vector relative to point 'a' on triangle + Vector3D ap = projection.subtract(t[0]); + + /* + Vector ap can now be compared to vector sides of triangle 't' using dot products to + determine where it lies in relation to the triangle, as shown in this image: + https://i.sstatic.net/tPiEB.png + Regions from this image will be referenced below + as discussed here: + https://stackoverflow.com/questions/2924795/fastest-way-to-compute-point-to-triangle-distance-in-3d + */ + + //region nearest point in triangle t is corners + final double abDOTap = ab.dotProduct(ap); + final double acDOTap = ac.dotProduct(ap); + + //#1 in https://i.sstatic.net/tPiEB.png + if (abDOTap <= 0d && acDOTap <= 0d) return t[0]; + + final Vector3D bc = t[2].subtract(t[1]); + final Vector3D bp = projection.subtract(t[1]); + + final double baDOTbp = ab.negate().dotProduct(bp); + final double bcDOTbp = bc.dotProduct(bp); + //#2 in https://i.sstatic.net/tPiEB.png + if (baDOTbp <= 0d && bcDOTbp <= 0d) return t[1]; + + final Vector3D cp = projection.subtract(t[2]); + final double cbDOTcp = bc.negate().dotProduct(cp); + final double caDOTcp = ac.negate().dotProduct(cp); + //#3 in https://i.sstatic.net/tPiEB.png + if (cbDOTcp <= 0d && caDOTcp <= 0d) return t[2]; + //endregion + + // Compute barycentric coordinates (v, w) of projection point + double acDOTac = ac.dotProduct(ac); + double abDOTac = ab.dotProduct(ac); + double abDOTab = ab.dotProduct(ab); + + double denom = (acDOTac * abDOTab - abDOTac * abDOTac); + if (Math.abs(denom) < 1.0e-30) { + return new Vector3D(Double.NEGATIVE_INFINITY,Double.NEGATIVE_INFINITY,Double.NEGATIVE_INFINITY); // Triangle is degenerate + } + + double w = (acDOTac * abDOTap - abDOTac * acDOTap)/denom; //coordinate towards b from a + double v = (abDOTab * acDOTap - abDOTac * abDOTap)/denom; //coordinate towards c from a + + // Check barycentric coordinates + if ((v >= 0) && (w >= 0) && (v + w <= 1)) { + // Projection point is in triangle; #0 in https://i.sstatic.net/tPiEB.png + return projection; + } + + //region nearest point in triangle t is on side + if(w <= 0 && v > w){ + //#4 in https://i.sstatic.net/tPiEB.png + return t[0].add(ab.scalarMultiply(v)); + } + + if(v <= 0 && w > v){ + //#5 in https://i.sstatic.net/tPiEB.png + return t[0].add(ac.scalarMultiply(w)); + } + + if(v + w > 1){ + //#6 in https://i.sstatic.net/tPiEB.png + final double scalarValue = bcDOTbp/bc.getNormSq(); + return t[1].add(bc.scalarMultiply(scalarValue)); + } + //endregion + + if (v <=0 && w <= 0){ //this should be redundant, but for some reason isn't + //#1 in https://i.sstatic.net/tPiEB.png + return t[0]; + } + return new Vector3D(Double.NEGATIVE_INFINITY,Double.NEGATIVE_INFINITY,Double.NEGATIVE_INFINITY); + } +} diff --git a/src/test/java/net/imagej/ops/geom/MeshFeatureTests.java b/src/test/java/net/imagej/ops/geom/MeshFeatureTests.java index b1bf3df1a..7f4b3fb34 100644 --- a/src/test/java/net/imagej/ops/geom/MeshFeatureTests.java +++ b/src/test/java/net/imagej/ops/geom/MeshFeatureTests.java @@ -31,28 +31,22 @@ import static org.junit.Assert.assertEquals; import static org.junit.Assert.assertTrue; +import java.util.Arrays; import java.util.Iterator; import net.imagej.mesh.Mesh; import net.imagej.mesh.Triangle; import net.imagej.ops.Ops; import net.imagej.ops.features.AbstractFeatureTest; -import net.imagej.ops.geom.geom3d.DefaultBoxivityMesh; -import net.imagej.ops.geom.geom3d.DefaultCompactness; -import net.imagej.ops.geom.geom3d.DefaultConvexityMesh; -import net.imagej.ops.geom.geom3d.DefaultMainElongation; -import net.imagej.ops.geom.geom3d.DefaultMarchingCubes; -import net.imagej.ops.geom.geom3d.DefaultMedianElongation; -import net.imagej.ops.geom.geom3d.DefaultSolidityMesh; -import net.imagej.ops.geom.geom3d.DefaultSparenessMesh; -import net.imagej.ops.geom.geom3d.DefaultSphericity; -import net.imagej.ops.geom.geom3d.DefaultSurfaceArea; -import net.imagej.ops.geom.geom3d.DefaultSurfaceAreaConvexHullMesh; -import net.imagej.ops.geom.geom3d.DefaultVerticesCountConvexHullMesh; -import net.imagej.ops.geom.geom3d.DefaultVerticesCountMesh; -import net.imagej.ops.geom.geom3d.DefaultVolumeConvexHullMesh; -import net.imagej.ops.geom.geom3d.DefaultVolumeMesh; +import net.imagej.ops.geom.geom3d.*; +import net.imglib2.Cursor; +import net.imglib2.img.Img; +import net.imglib2.img.array.ArrayImgFactory; +import net.imglib2.roi.Regions; +import net.imglib2.roi.boundary.Boundary; import net.imglib2.roi.labeling.LabelRegion; +import net.imglib2.type.logic.BitType; +import net.imglib2.type.logic.BoolType; import net.imglib2.type.numeric.real.DoubleType; import org.junit.BeforeClass; @@ -200,7 +194,28 @@ public void verticesCountMesh() { } @Test - public void voxelization3D() { - // https://github.com/imagej/imagej-ops/issues/422 + public void voxelization3D(){ + /** + An ideal voxelization algorithm should be able to convert a {@link Mesh} generated from a binary image back into + a surface-pixel outline of the original binary image. This surface-pixel image should match the result of + processing the original binary image with {@link Boundary} using + {@link Boundary.StructuringElement.FOUR_CONNECTED}. When working with real images these can mismatch due to + isolated single-pixel objects in the original binary image not being incorporated into the {@link Mesh} + */ + final Img out = new ArrayImgFactory<>(new BitType()).create(getTestImage3D()); + ops.run(EuclideanDistanceVoxelization3D.class,out, mesh, 1.0); + final Boundary compareTo = new Boundary(ops.convert().bit(getTestImage3D()), Boundary.StructuringElement.FOUR_CONNECTED); + boolean matches = true; + Cursor voxelizedPositive = Regions.iterable(out).localizingCursor(); + Cursor boundaryPositive = compareTo.localizingCursor(); + while(voxelizedPositive.hasNext() && boundaryPositive.hasNext() && matches){ + voxelizedPositive.fwd(); + boundaryPositive.fwd(); + if(!Arrays.equals(voxelizedPositive.positionAsLongArray(),boundaryPositive.positionAsLongArray())) + matches = false; + } + if(voxelizedPositive.hasNext() || boundaryPositive.hasNext()) + matches = false; + assertTrue(matches); } }