From 83b5c4c2eff32586968aafe6c0d8abc5fa26061f Mon Sep 17 00:00:00 2001 From: andmccall Date: Fri, 13 Feb 2026 14:06:29 -0500 Subject: [PATCH 01/13] Overhaul of Voxelization ops with new algorithm. Also included test --- .../net/imagej/ops/geom/GeomNamespace.java | 19 +- .../geom/geom3d/DefaultVoxelization3D.java | 426 ++++++------------ .../net/imagej/ops/geom/MeshFeatureTests.java | 29 +- 3 files changed, 154 insertions(+), 320 deletions(-) diff --git a/src/main/java/net/imagej/ops/geom/GeomNamespace.java b/src/main/java/net/imagej/ops/geom/GeomNamespace.java index 6e0d449f0..88d33d6ec 100644 --- a/src/main/java/net/imagej/ops/geom/GeomNamespace.java +++ b/src/main/java/net/imagej/ops/geom/GeomNamespace.java @@ -37,6 +37,7 @@ import net.imagej.ops.OpMethod; import net.imagej.ops.Ops.Geometric.Voxelization; import net.imagej.ops.geom.geom3d.mesh.VertexInterpolator; +import net.imglib2.Interval; import net.imglib2.IterableInterval; import net.imglib2.RandomAccessibleInterval; import net.imglib2.RealLocalizable; @@ -278,13 +279,6 @@ public List convexHull(final Mesh in) { net.imagej.ops.Ops.Geometric.ConvexHull.class, in); return result; } - - @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; - } @OpMethod(op = net.imagej.ops.geom.geom2d.DefaultConvexityPolygon.class) public DoubleType convexity(final Polygon2D in) { @@ -678,26 +672,25 @@ public double[] vertexInterpolator(final int[] p1, final int[] p2, @OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class) public RandomAccessibleInterval voxelization(final Mesh in) { - @SuppressWarnings("unchecked") final RandomAccessibleInterval result = (RandomAccessibleInterval) ops().run(net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class, in); return result; } @OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class) - public RandomAccessibleInterval voxelization(final Mesh in, final int width) { - @SuppressWarnings("unchecked") + public RandomAccessibleInterval voxelization(final Mesh in, final Interval dimensions) { 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, dimensions); return result; } @OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class) - public RandomAccessibleInterval voxelization(final Mesh in, final int width, final int height) { + public RandomAccessibleInterval voxelization(final Mesh in, final Interval dimensions, boolean scaleMeshToDimesions) { 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, dimensions, scaleMeshToDimesions); return result; } + @OpMethod(op = net.imagej.ops.geom.geom2d.DefaultVerticesCountPolygon.class) public DoubleType verticesCount(final Polygon2D in) { 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..884a8cec2 100644 --- a/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java +++ b/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java @@ -29,21 +29,20 @@ package net.imagej.ops.geom.geom3d; import net.imagej.mesh.Mesh; +import net.imagej.mesh.Meshes; import net.imagej.mesh.Triangle; -import net.imagej.mesh.Vertices; import net.imagej.ops.OpService; import net.imagej.ops.Ops; import net.imagej.ops.special.function.AbstractUnaryFunctionOp; -import net.imglib2.FinalInterval; -import net.imglib2.RandomAccess; -import net.imglib2.RandomAccessibleInterval; -import net.imglib2.RealLocalizable; -import net.imglib2.RealPoint; +import net.imglib2.*; import net.imglib2.img.Img; +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.log.LogService; import org.scijava.plugin.Parameter; import org.scijava.plugin.Plugin; @@ -52,353 +51,198 @@ * This is a voxelizer that produces a binary image with values filled in along * the surface of the mesh. *

- *

- * Thanks to Tomas Möller for sharing his public domain code: - * http://fileadmin.cs.lth.se/cs/personal/tomas_akenine-moller/code/tribox.txt - *

* - * @author Kyle Harrington (University of Idaho) + * @author Andrew McCall (University at Buffalo) */ @Plugin(type = Ops.Geometric.Voxelization.class) public class DefaultVoxelization3D extends AbstractUnaryFunctionOp> implements Ops.Geometric.Voxelization { @Parameter(type = ItemIO.INPUT, required = false) - private int width = 10; - - @Parameter(type = ItemIO.INPUT, required = false) - private int height = 10; + private Interval dimensions; @Parameter(type = ItemIO.INPUT, required = false) - private int depth = 10; + private boolean scaleMeshToDimesions = false; @Parameter private OpService ops; + @Parameter + private LogService logService; + + private final double wallThickness = 1.0; //This could be made into a parameter if needed + private final long[] offset = new long[] {0,0,0}; + private double scale = 1.0; + @Override public RandomAccessibleInterval calculate(Mesh input) { - Img outImg = ops.create().img(new FinalInterval(width, height, depth), new BitType()); - - Vertices verts = input.vertices(); + if(dimensions == null) { + float[] bounds = Meshes.boundingBox(input); + long[] outputInterval = new long[3]; + for (int i = 0; i < 3; i++) { + outputInterval[i] = (long)Math.ceil(bounds[i+3]-bounds[i]); + } + dimensions = new FinalInterval(outputInterval); + setScale(input, bounds); + scaleMeshToDimesions = false; + } - RealPoint minPoint = new RealPoint(verts.iterator().next()); - RealPoint maxPoint = new RealPoint(verts.iterator().next()); + if(scaleMeshToDimesions) + setScale(input); - for (RealLocalizable v : verts) { - if (v.getDoublePosition(0) < minPoint.getDoublePosition(0)) - minPoint.setPosition(v.getDoublePosition(0), 0); - if (v.getDoublePosition(1) < minPoint.getDoublePosition(1)) - minPoint.setPosition(v.getDoublePosition(1), 1); - if (v.getDoublePosition(2) < minPoint.getDoublePosition(2)) - minPoint.setPosition(v.getDoublePosition(2), 2); + Img outImg = ops.create().img(dimensions, new BitType()); - if (v.getDoublePosition(0) > maxPoint.getDoublePosition(0)) - maxPoint.setPosition(v.getDoublePosition(0), 0); - if (v.getDoublePosition(1) > maxPoint.getDoublePosition(1)) - maxPoint.setPosition(v.getDoublePosition(1), 1); - if (v.getDoublePosition(2) > maxPoint.getDoublePosition(2)) - maxPoint.setPosition(v.getDoublePosition(2), 2); - } + RandomAccess ra = outImg.randomAccess(); - RealPoint dimPoint = new RealPoint((maxPoint.getDoublePosition(0) - minPoint.getDoublePosition(0)), - (maxPoint.getDoublePosition(1) - minPoint.getDoublePosition(1)), - (maxPoint.getDoublePosition(2) - minPoint.getDoublePosition(2))); - - double[] stepSizes = new double[3]; - stepSizes[0] = dimPoint.getDoublePosition(0) / width; - stepSizes[1] = dimPoint.getDoublePosition(1) / height; - stepSizes[2] = dimPoint.getDoublePosition(2) / depth; - - double[] voxelHalfsize = new double[3]; - for (int k = 0; k < stepSizes.length; k++) - voxelHalfsize[k] = stepSizes[k] / 2.0; - - for (final Triangle tri : input.triangles()) { - final Vector3D v1 = new Vector3D(tri.v0x(), tri.v0y(), tri.v0z()); - final Vector3D v2 = new Vector3D(tri.v1x(), tri.v1y(), tri.v1z()); - final Vector3D v3 = new Vector3D(tri.v2x(), tri.v2y(), tri.v2z()); - - double[] minSubBoundary = new double[] { - Math.min(Math.min(v1.getX(), v2.getX()), v3.getX()) - minPoint.getDoublePosition(0), - Math.min(Math.min(v1.getY(), v2.getY()), v3.getY()) - minPoint.getDoublePosition(1), - Math.min(Math.min(v1.getZ(), v2.getZ()), v3.getZ()) - minPoint.getDoublePosition(2) }; - double[] maxSubBoundary = new double[] { - Math.max(Math.max(v1.getX(), v2.getX()), v3.getX()) - minPoint.getDoublePosition(0), - Math.max(Math.max(v1.getY(), v2.getY()), v3.getY()) - minPoint.getDoublePosition(1), - Math.max(Math.max(v1.getZ(), v2.getZ()), v3.getZ()) - minPoint.getDoublePosition(2) }; - - RandomAccess ra = outImg.randomAccess();// Should use the - // interval - // implementation - // for speed - - long[] indices = new long[3]; - for (indices[0] = (long) Math.floor(minSubBoundary[0] / stepSizes[0]); indices[0] < Math - .floor(maxSubBoundary[0] / stepSizes[0]); indices[0]++) { - for (indices[1] = (long) Math.floor(minSubBoundary[1] / stepSizes[1]); indices[1] < Math - .floor(maxSubBoundary[1] / stepSizes[1]); indices[1]++) { - for (indices[2] = (long) Math.floor(minSubBoundary[2] / stepSizes[2]); indices[2] < Math - .floor(maxSubBoundary[2] / stepSizes[2]); indices[2]++) { - ra.setPosition(indices); - if (!ra.get().get())// Don't check if voxel is already - // filled - { - double[] voxelCenter = new double[3]; - - for (int k = 0; k < 3; k++) - voxelCenter[k] = indices[k] * stepSizes[k] + voxelHalfsize[k]; - - if (triBoxOverlap(voxelCenter, voxelHalfsize, v1, v2, v3) == 1) { - ra.get().set(true); - } + input.triangles().forEach((Triangle t) -> { + Vector3D[] scaledT = scaleTriangleToOutput(t); + Interval triangleBox = boundingBox(scaledT); + LocalizingIntervalIterator it = new LocalizingIntervalIterator(triangleBox); + while (it.hasNext()) { + it.fwd(); + if(Intervals.contains(dimensions, it.positionAsPoint())) { + if (pointToTriangleDist(new Vector3D(it.getDoublePosition(0), it.getDoublePosition(1), it.getDoublePosition(2)), scaledT) < wallThickness/2) { + synchronized (ra) { + ra.setPositionAndGet(it.positionAsPoint()).set(true); } } } } - } + }); return outImg; } - private double findMin(double x0, double x1, double x2) { - return Math.min(Math.min(x0, x1), x2); - } + private Vector3D[] scaleTriangleToOutput(Triangle t){ + Vector3D[] o = new Vector3D[3]; - private double findMax(double x0, double x1, double x2) { - return Math.max(Math.max(x0, x1), x2); + o[0] = new Vector3D((t.v0x()-offset[0])*scale, (t.v0y()-offset[1])*scale, (t.v0z()-offset[2])*scale); + o[1] = new Vector3D((t.v1x()-offset[0])*scale, (t.v1y()-offset[1])*scale, (t.v1z()-offset[2])*scale); + o[2] = new Vector3D((t.v2x()-offset[0])*scale, (t.v2y()-offset[1])*scale, (t.v2z()-offset[2])*scale); + return o; } - private double dotArray(double[] v1, double[] v2) { - return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]; + private void setScale(Mesh input){ + setScale(input, Meshes.boundingBox(input)); } - private int planeBoxOverlap(double[] normalArray, double[] vertArray, double[] maxboxArray) { - double[] vminArray = new double[3]; - double[] vmaxArray = new double[3]; - for (int q = 0; q <= 2; q++) { - double v = vertArray[q]; - if (normalArray[q] > 0.0F) { - vminArray[q] = (-maxboxArray[q] - v); - maxboxArray[q] -= v; - } else { - maxboxArray[q] -= v; - vmaxArray[q] = (-maxboxArray[q] - v); - } - } - if (dotArray(normalArray, vminArray) > 0.0F) { - return 0; + private void setScale(Mesh input, float[] bounds) { + double[] axisScaling = new double[3]; + for (int i = 0; i < 3; i++) { + offset[i] = (long) Math.floor(bounds[i]) - dimensions.min(i); + axisScaling[i] = (dimensions.max(i) - dimensions.min(i)) / (bounds[i + 3] - bounds[i]);//dimensionRange/boundingboxRange } - if (dotArray(normalArray, vmaxArray) >= 0.0F) { - return 1; - } - return 0; + scale = Math.min(axisScaling[0], Math.min(axisScaling[1], axisScaling[2])); } - private int axisTest_x01(double e0, double e02, double fez, double fey, double[] v0, double[] v1, double[] v2, - double[] boxhalfsize) { - double p0 = e0 * v0[1] - e02 * v0[2]; - double p2 = e0 * v2[1] - e02 * v2[2]; - double max; - double min; - - if (p0 < p2) { - min = p0; - max = p2; - } else { - min = p2; - max = p0; - } - double rad = fez * boxhalfsize[1] + fey * boxhalfsize[2]; - if ((min > rad) || (max < -rad)) { - return 0; - } - return 1; - } - private int axisTest_x2(double a, double b, double fa, double fb, double[] v0, double[] v1, double[] v2, - double[] boxhalfsize) { - double p0 = a * v0[1] - b * v0[2]; - double p1 = a * v1[1] - b * v1[2]; - double max; - double min; - - if (p0 < p1) { - min = p0; - max = p1; - } else { - min = p1; - max = p0; - } - double rad = fa * boxhalfsize[1] + fb * boxhalfsize[2]; - if ((min > rad) || (max < -rad)) { - return 0; - } - return 1; - } + private Interval boundingBox(Vector3D[] t){ + long [] min = new long[3]; + long [] max = new long[3]; - private int axisTest_y02(double a, double b, double fa, double fb, double[] v0, double[] v1, double[] v2, - double[] boxhalfsize) { - double p0 = -a * v0[0] + b * v0[2]; - double p2 = -a * v2[0] + b * v2[2]; - double max; - double min; - - if (p0 < p2) { - min = p0; - max = p2; - } else { - min = p2; - max = p0; - } - double rad = fa * boxhalfsize[0] + fb * boxhalfsize[2]; - if ((min > rad) || (max < -rad)) { - return 0; - } - return 1; - } + 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 int axisTest_y1(double a, double b, double fa, double fb, double[] v0, double[] v1, double[] v2, - double[] boxhalfsize) { - double p0 = -a * v0[0] + b * v0[2]; - double p1 = -a * v1[0] + b * v1[2]; - double max; - double min; - - if (p0 < p1) { - min = p0; - max = p1; - } else { - min = p1; - max = p0; - } - double rad = fa * boxhalfsize[0] + fb * boxhalfsize[2]; - if ((min > rad) || (max < -rad)) { - return 0; - } - return 1; } - private int axisTest_z12(double a, double b, double fa, double fb, double[] v0, double[] v1, double[] v2, - double[] boxhalfsize) { - double p1 = a * v1[0] - b * v1[1]; - double p2 = a * v2[0] - b * v2[1]; - double max; - double min; - - if (p2 < p1) { - min = p2; - max = p1; - } else { - min = p1; - max = p2; - } - double rad = fa * boxhalfsize[0] + fb * boxhalfsize[1]; - if ((min > rad) || (max < -rad)) { - return 0; - } - return 1; + private double pointToTriangleDist(Vector3D p, Vector3D[] t){ + Vector3D tPoint = nearestPointInTriangle3D(p, t); + return p.distance(tPoint); } - private int axisTest_z0(double a, double b, double fa, double fb, double[] v0, double[] v1, double[] v2, - double[] boxhalfsize) { - double p0 = a * v0[0] - b * v0[1]; - double p1 = a * v1[0] - b * v1[1]; - double max; - double min; - - if (p0 < p1) { - min = p0; - max = p1; - } else { - min = p1; - max = p0; + + + private Vector3D nearestPointInTriangle3D(Vector3D p, Vector3D[] t) { + + Vector3D ab = t[1].subtract(t[0]); + Vector3D ac = t[2].subtract(t[0]); + + + //region Obtain projection (p) of origP onto plane of triangle + // Find the normal to the plane: n = (b - a) x (c - a) + Vector3D n = ab.crossProduct(ac); + + // Normalize normal vector + try{ + n = n.normalize(); } - double rad = fa * boxhalfsize[0] + fb * boxhalfsize[1]; - if ((min > rad) || (max < -rad)) { - return 0; + catch(Exception e){ + return new Vector3D(-100,-100,-100); // Triangle is degenerate } - return 1; - } - private void sub(double[] v0, double[] vert1, double[] boxcenter) { - vert1[0] -= boxcenter[0]; - vert1[1] -= boxcenter[1]; - vert1[2] -= boxcenter[2]; - } + // Project point origP 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 - private void cross(double[] dest, double[] v1, double[] v2) { - dest[0] = (v1[1] * v2[2] - v1[2] * v2[1]); - dest[1] = (v1[2] * v2[0] - v1[0] * v2[2]); - dest[2] = (v1[0] * v2[1] - v1[1] * v2[0]); - } + Vector3D ap = projection.subtract(t[0]); - private int triBoxOverlap(double[] boxcenter, double[] boxhalfsize, Vector3D pf1, Vector3D pf2, Vector3D pf3) { - double[] vert1 = pf1.toArray(); - double[] vert2 = pf2.toArray(); - double[] vert3 = pf3.toArray(); + //region nearest point is corners + final double abDOTap = ab.dotProduct(ap); + final double acDOTap = ac.dotProduct(ap); - double[] v0 = new double[3]; - double[] v1 = new double[3]; - double[] v2 = new double[3]; + if (abDOTap <= 0d && acDOTap <= 0d) return t[0]; - double[] normal = new double[3]; - double[] e0 = new double[3]; - double[] e1 = new double[3]; - double[] e2 = new double[3]; + final Vector3D bc = t[2].subtract(t[1]); + final Vector3D bp = projection.subtract(t[1]); - sub(v0, vert1, boxcenter); - sub(v1, vert2, boxcenter); - sub(v2, vert3, boxcenter); + final double baDOTbp = ab.negate().dotProduct(bp); + final double bcDOTbp = bc.dotProduct(bp); + if (baDOTbp <= 0d && bcDOTbp <= 0d) return t[1]; - sub(e0, v1, v0); - sub(e1, v2, v1); - sub(e2, v0, v2); - double fex = Math.abs(e0[0]); - double fey = Math.abs(e0[1]); - double fez = Math.abs(e0[2]); + final Vector3D cp = projection.subtract(t[2]); + final double cbDOTcp = bc.negate().dotProduct(cp); + final double caDOTcp = ac.negate().dotProduct(cp); + if (cbDOTcp <= 0d && caDOTcp <= 0d) return t[2]; + //endregion - axisTest_x01(e0[2], e0[1], fez, fey, v0, v1, v2, boxhalfsize); - axisTest_y02(e0[2], e0[0], fez, fex, v0, v1, v2, boxhalfsize); - axisTest_z12(e0[1], e0[0], fey, fex, v0, v1, v2, boxhalfsize); + double acDOTac = ac.dotProduct(ac); + double abDOTac = ab.dotProduct(ac); + double abDOTab = ab.dotProduct(ab); - fex = Math.abs(e1[0]); - fey = Math.abs(e1[1]); - fez = Math.abs(e1[2]); + // Compute barycentric coordinates (v, w) of projection point + double denom = (acDOTac * abDOTab - abDOTac *abDOTac); + if (Math.abs(denom) < 1.0e-30) { + return new Vector3D(-100,-100,-100); // Triangle is degenerate + } - axisTest_x01(e1[2], e1[1], fez, fey, v0, v1, v2, boxhalfsize); - axisTest_y02(e1[2], e1[0], fez, fex, v0, v1, v2, boxhalfsize); - axisTest_z0(e1[1], e1[0], fey, fex, v0, v1, v2, boxhalfsize); + double w = (acDOTac * abDOTap - abDOTac * acDOTap)/denom; //coordinate towards b from a + double v = (abDOTab * acDOTap - abDOTac * abDOTap)/denom; //coordinate towards c from a - fex = Math.abs(e2[0]); - fey = Math.abs(e2[1]); - fez = Math.abs(e2[2]); - axisTest_x2(e2[2], e2[1], fez, fey, v0, v1, v2, boxhalfsize); - axisTest_y1(e2[2], e2[0], fez, fex, v0, v1, v2, boxhalfsize); - axisTest_z12(e2[1], e2[0], fey, fex, v0, v1, v2, boxhalfsize); + // Check barycentric coordinates + if ((v >= 0) && (w >= 0) && (v + w <= 1)) { + // Nearest orthogonal projection point is in triangle + return projection; + } - double min = findMin(v0[0], v1[0], v2[0]); - double max = findMax(v0[0], v1[0], v2[0]); - if ((min > boxhalfsize[0]) || (max < -boxhalfsize[0])) { - return 0; + //region nearest point is on edge + if(w <= 0 && v > w){ + return t[0].add(ab.scalarMultiply(v)); } - min = findMin(v0[1], v1[1], v2[1]); - max = findMax(v0[1], v1[1], v2[1]); - if ((min > boxhalfsize[1]) || (max < -boxhalfsize[1])) { - return 0; + + if(v <= 0 && w > v){ + return t[0].add(ac.scalarMultiply(w)); } - min = findMin(v0[2], v1[2], v2[2]); - max = findMax(v0[2], v1[2], v2[2]); - if ((min > boxhalfsize[2]) || (max < -boxhalfsize[2])) { - return 0; + + if(v + w > 1){ + final double scalarValue = bcDOTbp/bc.getNormSq(); + return t[1].add(bc.scalarMultiply(scalarValue)); } - cross(normal, e0, e1); - if (planeBoxOverlap(normal, v0, boxhalfsize) != 1) { - return 0; + //endregion + + if (v <=0 && w <= 0){ //this should be redundant, but for some reason isn't + return t[0]; } - return 1; - } + return new Vector3D(-100,-100,-100); + } } diff --git a/src/test/java/net/imagej/ops/geom/MeshFeatureTests.java b/src/test/java/net/imagej/ops/geom/MeshFeatureTests.java index b1bf3df1a..5d5b98875 100644 --- a/src/test/java/net/imagej/ops/geom/MeshFeatureTests.java +++ b/src/test/java/net/imagej/ops/geom/MeshFeatureTests.java @@ -37,22 +37,12 @@ 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.RandomAccessibleInterval; +import net.imglib2.algorithm.neighborhood.DiamondShape; +import net.imglib2.roi.Regions; import net.imglib2.roi.labeling.LabelRegion; +import net.imglib2.type.logic.BitType; import net.imglib2.type.numeric.real.DoubleType; import org.junit.BeforeClass; @@ -201,6 +191,13 @@ public void verticesCountMesh() { @Test public void voxelization3D() { - // https://github.com/imagej/imagej-ops/issues/422 + /*Value of 184 here corresponds with: + RandomAccessibleInterval result = (RandomAccessibleInterval) ops.run(DefaultVoxelization3D.class, mesh, ROI); + ops.morphology().fillHoles(result, result, new DiamondShape(1)); + assertEquals(Ops.Geometric.Voxelization.NAME, ROI.size(), Regions.countTrue(result)); + */ + assertEquals(Ops.Geometric.Voxelization.NAME, 184, + Regions.countTrue((RandomAccessibleInterval) ops.run(DefaultVoxelization3D.class, mesh, ROI))); + } } From ff5da48e57251c645ccbc8624e9447219d726db3 Mon Sep 17 00:00:00 2001 From: andmccall <56968719+andmccall@users.noreply.github.com> Date: Fri, 13 Feb 2026 14:25:03 -0500 Subject: [PATCH 02/13] Remove unused LogService parameter Removed unused LogService parameter from DefaultVoxelization3D. --- .../java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java | 3 --- 1 file changed, 3 deletions(-) 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 884a8cec2..2b43a8458 100644 --- a/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java +++ b/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java @@ -67,9 +67,6 @@ public class DefaultVoxelization3D extends AbstractUnaryFunctionOp Date: Mon, 16 Feb 2026 11:21:46 -0500 Subject: [PATCH 03/13] Minor bug fix --- .../net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) 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 884a8cec2..209355189 100644 --- a/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java +++ b/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java @@ -67,9 +67,6 @@ public class DefaultVoxelization3D extends AbstractUnaryFunctionOp calculate(Mesh input) { float[] bounds = Meshes.boundingBox(input); long[] outputInterval = new long[3]; for (int i = 0; i < 3; i++) { - outputInterval[i] = (long)Math.ceil(bounds[i+3]-bounds[i]); + outputInterval[i] = (long)Math.ceil(bounds[i+3]-bounds[i])+1; } dimensions = new FinalInterval(outputInterval); setScale(input, bounds); From ff45913722d1deb793a79296d06d8dcdf4d8890a Mon Sep 17 00:00:00 2001 From: andmccall <56968719+andmccall@users.noreply.github.com> Date: Mon, 16 Feb 2026 11:25:48 -0500 Subject: [PATCH 04/13] Fix output interval calculation in voxelization Giving extra space on all sides when calculating the dimensions and offset seems to allow fillHoles to work more consistently after voxelization. --- .../java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 2b43a8458..2fd18da6a 100644 --- a/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java +++ b/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java @@ -78,7 +78,7 @@ public RandomAccessibleInterval calculate(Mesh input) { float[] bounds = Meshes.boundingBox(input); long[] outputInterval = new long[3]; for (int i = 0; i < 3; i++) { - outputInterval[i] = (long)Math.ceil(bounds[i+3]-bounds[i]); + outputInterval[i] = (long)Math.ceil(bounds[i+3]-bounds[i])+2; } dimensions = new FinalInterval(outputInterval); setScale(input, bounds); From b004747ed6618a669d82e6ff48ddbecba4937a16 Mon Sep 17 00:00:00 2001 From: andmccall Date: Mon, 16 Feb 2026 12:53:38 -0500 Subject: [PATCH 05/13] Minor bug fix --- .../java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 209355189..2fd18da6a 100644 --- a/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java +++ b/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java @@ -78,7 +78,7 @@ public RandomAccessibleInterval calculate(Mesh input) { float[] bounds = Meshes.boundingBox(input); long[] outputInterval = new long[3]; for (int i = 0; i < 3; i++) { - outputInterval[i] = (long)Math.ceil(bounds[i+3]-bounds[i])+1; + outputInterval[i] = (long)Math.ceil(bounds[i+3]-bounds[i])+2; } dimensions = new FinalInterval(outputInterval); setScale(input, bounds); From 3b49dafd16a69e43be430380187ab4a2054a7c90 Mon Sep 17 00:00:00 2001 From: andmccall Date: Mon, 16 Feb 2026 14:07:23 -0500 Subject: [PATCH 06/13] Minor bug fix --- .../net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 2fd18da6a..b45a6894e 100644 --- a/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java +++ b/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java @@ -78,7 +78,7 @@ public RandomAccessibleInterval calculate(Mesh input) { float[] bounds = Meshes.boundingBox(input); long[] outputInterval = new long[3]; for (int i = 0; i < 3; i++) { - outputInterval[i] = (long)Math.ceil(bounds[i+3]-bounds[i])+2; + outputInterval[i] = (long)Math.ceil(bounds[i+3]-bounds[i]); } dimensions = new FinalInterval(outputInterval); setScale(input, bounds); @@ -128,7 +128,7 @@ private void setScale(Mesh input, float[] bounds) { double[] axisScaling = new double[3]; for (int i = 0; i < 3; i++) { offset[i] = (long) Math.floor(bounds[i]) - dimensions.min(i); - axisScaling[i] = (dimensions.max(i) - dimensions.min(i)) / (bounds[i + 3] - bounds[i]);//dimensionRange/boundingboxRange + axisScaling[i] = (dimensions.max(i) - dimensions.min(i)) / ((bounds[i + 3]+2) - bounds[i]); } scale = Math.min(axisScaling[0], Math.min(axisScaling[1], axisScaling[2])); } From c3f0c41031acebb8cd30cf38deb01c0cb06e8f77 Mon Sep 17 00:00:00 2001 From: andmccall Date: Tue, 17 Feb 2026 13:52:25 -0500 Subject: [PATCH 07/13] Added wallThickness as a parameter, and changed automatically generated dimensions to make room for thicker walls --- .../net/imagej/ops/geom/GeomNamespace.java | 22 ++++++++++++++++ .../geom/geom3d/DefaultVoxelization3D.java | 25 ++++++++----------- 2 files changed, 33 insertions(+), 14 deletions(-) diff --git a/src/main/java/net/imagej/ops/geom/GeomNamespace.java b/src/main/java/net/imagej/ops/geom/GeomNamespace.java index 88d33d6ec..ac4f1a6dc 100644 --- a/src/main/java/net/imagej/ops/geom/GeomNamespace.java +++ b/src/main/java/net/imagej/ops/geom/GeomNamespace.java @@ -691,6 +691,28 @@ public RandomAccessibleInterval voxelization(final Mesh in, final Inter return result; } + //Two OpMethods below cause a Mismatched inputs error, though they both work fine if built with 'Skip tests' +// @OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class) +// public RandomAccessibleInterval voxelization(final Mesh in, double wallThickness) { +// final RandomAccessibleInterval result = +// (RandomAccessibleInterval) ops().run(net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class, in, null, false, wallThickness); +// return result; +// } +// +// @OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class) +// public RandomAccessibleInterval voxelization(final Mesh in, final Interval dimensions, double wallThickness) { +// final RandomAccessibleInterval result = +// (RandomAccessibleInterval) ops().run(net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class, in, dimensions, false, wallThickness); +// return result; +// } + + @OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class) + public RandomAccessibleInterval voxelization(final Mesh in, final Interval dimensions, boolean scaleMeshToDimesions, double wallThickness) { + final RandomAccessibleInterval result = + (RandomAccessibleInterval) ops().run(net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class, in, dimensions, scaleMeshToDimesions, wallThickness); + return result; + } + @OpMethod(op = net.imagej.ops.geom.geom2d.DefaultVerticesCountPolygon.class) public DoubleType verticesCount(final Polygon2D in) { 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 b45a6894e..e1c1250b2 100644 --- a/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java +++ b/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java @@ -42,13 +42,12 @@ import net.imglib2.util.Intervals; import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; import org.scijava.ItemIO; -import org.scijava.log.LogService; import org.scijava.plugin.Parameter; import org.scijava.plugin.Plugin; /** *

- * This is a voxelizer that produces a binary image with values filled in along + * This is a voxelizer that produces a binary image with values set to true along * the surface of the mesh. *

* @@ -67,7 +66,9 @@ public class DefaultVoxelization3D extends AbstractUnaryFunctionOp calculate(Mesh input) { float[] bounds = Meshes.boundingBox(input); long[] outputInterval = new long[3]; for (int i = 0; i < 3; i++) { - outputInterval[i] = (long)Math.ceil(bounds[i+3]-bounds[i]); + outputInterval[i] = (long)Math.ceil(bounds[i+3]+(2*wallThickness)-bounds[i]); } dimensions = new FinalInterval(outputInterval); - setScale(input, bounds); + setScale(bounds); scaleMeshToDimesions = false; } @@ -99,10 +100,8 @@ public RandomAccessibleInterval calculate(Mesh input) { while (it.hasNext()) { it.fwd(); if(Intervals.contains(dimensions, it.positionAsPoint())) { - if (pointToTriangleDist(new Vector3D(it.getDoublePosition(0), it.getDoublePosition(1), it.getDoublePosition(2)), scaledT) < wallThickness/2) { - synchronized (ra) { + if (pointToTriangleDist(new Vector3D(it.getDoublePosition(0), it.getDoublePosition(1), it.getDoublePosition(2)), scaledT) <= wallThickness/2) { ra.setPositionAndGet(it.positionAsPoint()).set(true); - } } } } @@ -121,14 +120,14 @@ private Vector3D[] scaleTriangleToOutput(Triangle t){ } private void setScale(Mesh input){ - setScale(input, Meshes.boundingBox(input)); + setScale(Meshes.boundingBox(input)); } - private void setScale(Mesh input, float[] bounds) { + private void setScale(float[] bounds) { double[] axisScaling = new double[3]; for (int i = 0; i < 3; i++) { - offset[i] = (long) Math.floor(bounds[i]) - dimensions.min(i); - axisScaling[i] = (dimensions.max(i) - dimensions.min(i)) / ((bounds[i + 3]+2) - bounds[i]); + offset[i] = Math.round(bounds[i] - (dimensions.min(i)+wallThickness)); + axisScaling[i] = ((dimensions.max(i)-2*wallThickness) - (dimensions.min(i))) / (bounds[i + 3] - bounds[i]); } scale = Math.min(axisScaling[0], Math.min(axisScaling[1], axisScaling[2])); } @@ -214,7 +213,6 @@ private Vector3D nearestPointInTriangle3D(Vector3D p, Vector3D[] t) { 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)) { // Nearest orthogonal projection point is in triangle @@ -239,7 +237,6 @@ private Vector3D nearestPointInTriangle3D(Vector3D p, Vector3D[] t) { if (v <=0 && w <= 0){ //this should be redundant, but for some reason isn't return t[0]; } - return new Vector3D(-100,-100,-100); } } From 34f8073cae85f3a19ef06d11f21af7b7fa768f45 Mon Sep 17 00:00:00 2001 From: andmccall Date: Fri, 3 Jul 2026 14:51:29 -0400 Subject: [PATCH 08/13] Moved new algorithm to Voxelize3D Op and Deprecated old Voxelization3D Op. Also created better Feature Test. --- .../net/imagej/ops/geom/GeomNamespace.java | 57 ++- .../geom/geom3d/DefaultVoxelization3D.java | 435 ++++++++++++------ .../ops/geom/geom3d/DefaultVoxelize3D.java | 230 +++++++++ .../net/imagej/ops/geom/MeshFeatureTests.java | 36 +- 4 files changed, 589 insertions(+), 169 deletions(-) create mode 100644 src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelize3D.java diff --git a/src/main/java/net/imagej/ops/geom/GeomNamespace.java b/src/main/java/net/imagej/ops/geom/GeomNamespace.java index ac4f1a6dc..6b0f64a4c 100644 --- a/src/main/java/net/imagej/ops/geom/GeomNamespace.java +++ b/src/main/java/net/imagej/ops/geom/GeomNamespace.java @@ -280,6 +280,14 @@ public List convexHull(final Mesh 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; + } + @OpMethod(op = net.imagej.ops.geom.geom2d.DefaultConvexityPolygon.class) public DoubleType convexity(final Polygon2D in) { final DoubleType result = @@ -669,50 +677,53 @@ public double[] vertexInterpolator(final int[] p1, final int[] p2, p1Value, p2Value); return result; } - + + @Deprecated @OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelization3D.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.geom.geom3d.DefaultVoxelization3D.class, in); return result; } + @Deprecated @OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class) - public RandomAccessibleInterval voxelization(final Mesh in, final Interval dimensions) { + 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, dimensions); + (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 Interval dimensions, boolean scaleMeshToDimesions) { + 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, dimensions, scaleMeshToDimesions); + (RandomAccessibleInterval) ops().run(net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class, in, width, height); return result; } - //Two OpMethods below cause a Mismatched inputs error, though they both work fine if built with 'Skip tests' -// @OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class) -// public RandomAccessibleInterval voxelization(final Mesh in, double wallThickness) { -// final RandomAccessibleInterval result = -// (RandomAccessibleInterval) ops().run(net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class, in, null, false, wallThickness); -// return result; -// } -// -// @OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class) -// public RandomAccessibleInterval voxelization(final Mesh in, final Interval dimensions, double wallThickness) { -// final RandomAccessibleInterval result = -// (RandomAccessibleInterval) ops().run(net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class, in, dimensions, false, wallThickness); -// return result; -// } + @OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelize3D.class) + public RandomAccessibleInterval voxelize(final Mesh in) { + final RandomAccessibleInterval result = + (RandomAccessibleInterval) ops().run(net.imagej.ops.geom.geom3d.DefaultVoxelize3D.class, in); + return result; + } - @OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class) - public RandomAccessibleInterval voxelization(final Mesh in, final Interval dimensions, boolean scaleMeshToDimesions, double wallThickness) { + @OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelize3D.class) + public RandomAccessibleInterval voxelize(final RandomAccessibleInterval out, final Mesh in) { final RandomAccessibleInterval result = - (RandomAccessibleInterval) ops().run(net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class, in, dimensions, scaleMeshToDimesions, wallThickness); + (RandomAccessibleInterval) ops().run(net.imagej.ops.geom.geom3d.DefaultVoxelize3D.class, out, in); return result; } + @OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelize3D.class) + public RandomAccessibleInterval voxelize(final RandomAccessibleInterval out, final Mesh in, final double wallThickness) { + final RandomAccessibleInterval result = + (RandomAccessibleInterval) ops().run(net.imagej.ops.geom.geom3d.DefaultVoxelize3D.class, out, in, wallThickness); + return result; + } @OpMethod(op = net.imagej.ops.geom.geom2d.DefaultVerticesCountPolygon.class) public DoubleType verticesCount(final Polygon2D in) { 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 e1c1250b2..78c360dae 100644 --- a/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java +++ b/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java @@ -6,13 +6,13 @@ * %% * 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 @@ -29,17 +29,19 @@ package net.imagej.ops.geom.geom3d; import net.imagej.mesh.Mesh; -import net.imagej.mesh.Meshes; import net.imagej.mesh.Triangle; +import net.imagej.mesh.Vertices; import net.imagej.ops.OpService; import net.imagej.ops.Ops; import net.imagej.ops.special.function.AbstractUnaryFunctionOp; -import net.imglib2.*; +import net.imglib2.FinalInterval; +import net.imglib2.RandomAccess; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.RealLocalizable; +import net.imglib2.RealPoint; import net.imglib2.img.Img; -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.plugin.Parameter; @@ -47,196 +49,357 @@ /** *

- * This is a voxelizer that produces a binary image with values set to true along + * This is a voxelizer that produces a binary image with values filled in along * the surface of the mesh. *

- * - * @author Andrew McCall (University at Buffalo) + *

+ * Thanks to Tomas Möller for sharing his public domain code: + * http://fileadmin.cs.lth.se/cs/personal/tomas_akenine-moller/code/tribox.txt + *

+ * + * @author Kyle Harrington (University of Idaho) + * @deprecated Use {@link DefaultVoxelize3D} instead. */ @Plugin(type = Ops.Geometric.Voxelization.class) public class DefaultVoxelization3D extends AbstractUnaryFunctionOp> implements Ops.Geometric.Voxelization { @Parameter(type = ItemIO.INPUT, required = false) - private Interval dimensions; + private int width = 10; @Parameter(type = ItemIO.INPUT, required = false) - private boolean scaleMeshToDimesions = false; - - @Parameter - private OpService ops; + private int height = 10; @Parameter(type = ItemIO.INPUT, required = false) - private double wallThickness = 1.0; + private int depth = 10; - private final long[] offset = new long[] {0,0,0}; - private double scale = 1.0; + @Parameter + private OpService ops; @Override public RandomAccessibleInterval calculate(Mesh input) { - if(dimensions == null) { - float[] bounds = Meshes.boundingBox(input); - long[] outputInterval = new long[3]; - for (int i = 0; i < 3; i++) { - outputInterval[i] = (long)Math.ceil(bounds[i+3]+(2*wallThickness)-bounds[i]); - } - dimensions = new FinalInterval(outputInterval); - setScale(bounds); - scaleMeshToDimesions = false; - } + Img outImg = ops.create().img(new FinalInterval(width, height, depth), new BitType()); - if(scaleMeshToDimesions) - setScale(input); + Vertices verts = input.vertices(); - Img outImg = ops.create().img(dimensions, new BitType()); + RealPoint minPoint = new RealPoint(verts.iterator().next()); + RealPoint maxPoint = new RealPoint(verts.iterator().next()); - RandomAccess ra = outImg.randomAccess(); + for (RealLocalizable v : verts) { + if (v.getDoublePosition(0) < minPoint.getDoublePosition(0)) + minPoint.setPosition(v.getDoublePosition(0), 0); + if (v.getDoublePosition(1) < minPoint.getDoublePosition(1)) + minPoint.setPosition(v.getDoublePosition(1), 1); + if (v.getDoublePosition(2) < minPoint.getDoublePosition(2)) + minPoint.setPosition(v.getDoublePosition(2), 2); - input.triangles().forEach((Triangle t) -> { - Vector3D[] scaledT = scaleTriangleToOutput(t); - Interval triangleBox = boundingBox(scaledT); - LocalizingIntervalIterator it = new LocalizingIntervalIterator(triangleBox); - while (it.hasNext()) { - it.fwd(); - if(Intervals.contains(dimensions, it.positionAsPoint())) { - if (pointToTriangleDist(new Vector3D(it.getDoublePosition(0), it.getDoublePosition(1), it.getDoublePosition(2)), scaledT) <= wallThickness/2) { - ra.setPositionAndGet(it.positionAsPoint()).set(true); + if (v.getDoublePosition(0) > maxPoint.getDoublePosition(0)) + maxPoint.setPosition(v.getDoublePosition(0), 0); + if (v.getDoublePosition(1) > maxPoint.getDoublePosition(1)) + maxPoint.setPosition(v.getDoublePosition(1), 1); + if (v.getDoublePosition(2) > maxPoint.getDoublePosition(2)) + maxPoint.setPosition(v.getDoublePosition(2), 2); + } + + RealPoint dimPoint = new RealPoint((maxPoint.getDoublePosition(0) - minPoint.getDoublePosition(0)), + (maxPoint.getDoublePosition(1) - minPoint.getDoublePosition(1)), + (maxPoint.getDoublePosition(2) - minPoint.getDoublePosition(2))); + + double[] stepSizes = new double[3]; + stepSizes[0] = dimPoint.getDoublePosition(0) / width; + stepSizes[1] = dimPoint.getDoublePosition(1) / height; + stepSizes[2] = dimPoint.getDoublePosition(2) / depth; + + double[] voxelHalfsize = new double[3]; + for (int k = 0; k < stepSizes.length; k++) + voxelHalfsize[k] = stepSizes[k] / 2.0; + + for (final Triangle tri : input.triangles()) { + final Vector3D v1 = new Vector3D(tri.v0x(), tri.v0y(), tri.v0z()); + final Vector3D v2 = new Vector3D(tri.v1x(), tri.v1y(), tri.v1z()); + final Vector3D v3 = new Vector3D(tri.v2x(), tri.v2y(), tri.v2z()); + + double[] minSubBoundary = new double[] { + Math.min(Math.min(v1.getX(), v2.getX()), v3.getX()) - minPoint.getDoublePosition(0), + Math.min(Math.min(v1.getY(), v2.getY()), v3.getY()) - minPoint.getDoublePosition(1), + Math.min(Math.min(v1.getZ(), v2.getZ()), v3.getZ()) - minPoint.getDoublePosition(2) }; + double[] maxSubBoundary = new double[] { + Math.max(Math.max(v1.getX(), v2.getX()), v3.getX()) - minPoint.getDoublePosition(0), + Math.max(Math.max(v1.getY(), v2.getY()), v3.getY()) - minPoint.getDoublePosition(1), + Math.max(Math.max(v1.getZ(), v2.getZ()), v3.getZ()) - minPoint.getDoublePosition(2) }; + + RandomAccess ra = outImg.randomAccess();// Should use the + // interval + // implementation + // for speed + + long[] indices = new long[3]; + for (indices[0] = (long) Math.floor(minSubBoundary[0] / stepSizes[0]); indices[0] < Math + .floor(maxSubBoundary[0] / stepSizes[0]); indices[0]++) { + for (indices[1] = (long) Math.floor(minSubBoundary[1] / stepSizes[1]); indices[1] < Math + .floor(maxSubBoundary[1] / stepSizes[1]); indices[1]++) { + for (indices[2] = (long) Math.floor(minSubBoundary[2] / stepSizes[2]); indices[2] < Math + .floor(maxSubBoundary[2] / stepSizes[2]); indices[2]++) { + ra.setPosition(indices); + if (!ra.get().get())// Don't check if voxel is already + // filled + { + double[] voxelCenter = new double[3]; + + for (int k = 0; k < 3; k++) + voxelCenter[k] = indices[k] * stepSizes[k] + voxelHalfsize[k]; + + if (triBoxOverlap(voxelCenter, voxelHalfsize, v1, v2, v3) == 1) { + ra.get().set(true); + } + } } } } - }); + } return outImg; } - private Vector3D[] scaleTriangleToOutput(Triangle t){ - Vector3D[] o = new Vector3D[3]; + private double findMin(double x0, double x1, double x2) { + return Math.min(Math.min(x0, x1), x2); + } - o[0] = new Vector3D((t.v0x()-offset[0])*scale, (t.v0y()-offset[1])*scale, (t.v0z()-offset[2])*scale); - o[1] = new Vector3D((t.v1x()-offset[0])*scale, (t.v1y()-offset[1])*scale, (t.v1z()-offset[2])*scale); - o[2] = new Vector3D((t.v2x()-offset[0])*scale, (t.v2y()-offset[1])*scale, (t.v2z()-offset[2])*scale); - return o; + private double findMax(double x0, double x1, double x2) { + return Math.max(Math.max(x0, x1), x2); } - private void setScale(Mesh input){ - setScale(Meshes.boundingBox(input)); + private double dotArray(double[] v1, double[] v2) { + return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]; } - private void setScale(float[] bounds) { - double[] axisScaling = new double[3]; - for (int i = 0; i < 3; i++) { - offset[i] = Math.round(bounds[i] - (dimensions.min(i)+wallThickness)); - axisScaling[i] = ((dimensions.max(i)-2*wallThickness) - (dimensions.min(i))) / (bounds[i + 3] - bounds[i]); + private int planeBoxOverlap(double[] normalArray, double[] vertArray, double[] maxboxArray) { + double[] vminArray = new double[3]; + double[] vmaxArray = new double[3]; + for (int q = 0; q <= 2; q++) { + double v = vertArray[q]; + if (normalArray[q] > 0.0F) { + vminArray[q] = (-maxboxArray[q] - v); + maxboxArray[q] -= v; + } else { + maxboxArray[q] -= v; + vmaxArray[q] = (-maxboxArray[q] - v); + } + } + if (dotArray(normalArray, vminArray) > 0.0F) { + return 0; } - scale = Math.min(axisScaling[0], Math.min(axisScaling[1], axisScaling[2])); + if (dotArray(normalArray, vmaxArray) >= 0.0F) { + return 1; + } + return 0; } - - 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 int axisTest_x01(double e0, double e02, double fez, double fey, double[] v0, double[] v1, double[] v2, + double[] boxhalfsize) { + double p0 = e0 * v0[1] - e02 * v0[2]; + double p2 = e0 * v2[1] - e02 * v2[2]; + double max; + double min; + + if (p0 < p2) { + min = p0; + max = p2; + } else { + min = p2; + max = p0; + } + double rad = fez * boxhalfsize[1] + fey * boxhalfsize[2]; + if ((min > rad) || (max < -rad)) { + return 0; + } + return 1; } - private double pointToTriangleDist(Vector3D p, Vector3D[] t){ - Vector3D tPoint = nearestPointInTriangle3D(p, t); - return p.distance(tPoint); + private int axisTest_x2(double a, double b, double fa, double fb, double[] v0, double[] v1, double[] v2, + double[] boxhalfsize) { + double p0 = a * v0[1] - b * v0[2]; + double p1 = a * v1[1] - b * v1[2]; + double max; + double min; + + if (p0 < p1) { + min = p0; + max = p1; + } else { + min = p1; + max = p0; + } + double rad = fa * boxhalfsize[1] + fb * boxhalfsize[2]; + if ((min > rad) || (max < -rad)) { + return 0; + } + return 1; } + private int axisTest_y02(double a, double b, double fa, double fb, double[] v0, double[] v1, double[] v2, + double[] boxhalfsize) { + double p0 = -a * v0[0] + b * v0[2]; + double p2 = -a * v2[0] + b * v2[2]; + double max; + double min; + + if (p0 < p2) { + min = p0; + max = p2; + } else { + min = p2; + max = p0; + } + double rad = fa * boxhalfsize[0] + fb * boxhalfsize[2]; + if ((min > rad) || (max < -rad)) { + return 0; + } + return 1; + } + private int axisTest_y1(double a, double b, double fa, double fb, double[] v0, double[] v1, double[] v2, + double[] boxhalfsize) { + double p0 = -a * v0[0] + b * v0[2]; + double p1 = -a * v1[0] + b * v1[2]; + double max; + double min; + + if (p0 < p1) { + min = p0; + max = p1; + } else { + min = p1; + max = p0; + } + double rad = fa * boxhalfsize[0] + fb * boxhalfsize[2]; + if ((min > rad) || (max < -rad)) { + return 0; + } + return 1; + } - private Vector3D nearestPointInTriangle3D(Vector3D p, Vector3D[] t) { - - Vector3D ab = t[1].subtract(t[0]); - Vector3D ac = t[2].subtract(t[0]); - - - //region Obtain projection (p) of origP onto plane of triangle - // Find the normal to the plane: n = (b - a) x (c - a) - Vector3D n = ab.crossProduct(ac); + private int axisTest_z12(double a, double b, double fa, double fb, double[] v0, double[] v1, double[] v2, + double[] boxhalfsize) { + double p1 = a * v1[0] - b * v1[1]; + double p2 = a * v2[0] - b * v2[1]; + double max; + double min; + + if (p2 < p1) { + min = p2; + max = p1; + } else { + min = p1; + max = p2; + } + double rad = fa * boxhalfsize[0] + fb * boxhalfsize[1]; + if ((min > rad) || (max < -rad)) { + return 0; + } + return 1; + } - // Normalize normal vector - try{ - n = n.normalize(); + private int axisTest_z0(double a, double b, double fa, double fb, double[] v0, double[] v1, double[] v2, + double[] boxhalfsize) { + double p0 = a * v0[0] - b * v0[1]; + double p1 = a * v1[0] - b * v1[1]; + double max; + double min; + + if (p0 < p1) { + min = p0; + max = p1; + } else { + min = p1; + max = p0; } - catch(Exception e){ - return new Vector3D(-100,-100,-100); // Triangle is degenerate + double rad = fa * boxhalfsize[0] + fb * boxhalfsize[1]; + if ((min > rad) || (max < -rad)) { + return 0; } + return 1; + } - // Project point origP 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 + private void sub(double[] v0, double[] vert1, double[] boxcenter) { + vert1[0] -= boxcenter[0]; + vert1[1] -= boxcenter[1]; + vert1[2] -= boxcenter[2]; + } - Vector3D ap = projection.subtract(t[0]); + private void cross(double[] dest, double[] v1, double[] v2) { + dest[0] = (v1[1] * v2[2] - v1[2] * v2[1]); + dest[1] = (v1[2] * v2[0] - v1[0] * v2[2]); + dest[2] = (v1[0] * v2[1] - v1[1] * v2[0]); + } - //region nearest point is corners - final double abDOTap = ab.dotProduct(ap); - final double acDOTap = ac.dotProduct(ap); + private int triBoxOverlap(double[] boxcenter, double[] boxhalfsize, Vector3D pf1, Vector3D pf2, Vector3D pf3) { + double[] vert1 = pf1.toArray(); + double[] vert2 = pf2.toArray(); + double[] vert3 = pf3.toArray(); - if (abDOTap <= 0d && acDOTap <= 0d) return t[0]; + double[] v0 = new double[3]; + double[] v1 = new double[3]; + double[] v2 = new double[3]; - final Vector3D bc = t[2].subtract(t[1]); - final Vector3D bp = projection.subtract(t[1]); + double[] normal = new double[3]; + double[] e0 = new double[3]; + double[] e1 = new double[3]; + double[] e2 = new double[3]; - final double baDOTbp = ab.negate().dotProduct(bp); - final double bcDOTbp = bc.dotProduct(bp); - if (baDOTbp <= 0d && bcDOTbp <= 0d) return t[1]; + sub(v0, vert1, boxcenter); + sub(v1, vert2, boxcenter); + sub(v2, vert3, boxcenter); + sub(e0, v1, v0); + sub(e1, v2, v1); + sub(e2, v0, v2); - final Vector3D cp = projection.subtract(t[2]); - final double cbDOTcp = bc.negate().dotProduct(cp); - final double caDOTcp = ac.negate().dotProduct(cp); - if (cbDOTcp <= 0d && caDOTcp <= 0d) return t[2]; - //endregion + double fex = Math.abs(e0[0]); + double fey = Math.abs(e0[1]); + double fez = Math.abs(e0[2]); - double acDOTac = ac.dotProduct(ac); - double abDOTac = ab.dotProduct(ac); - double abDOTab = ab.dotProduct(ab); + axisTest_x01(e0[2], e0[1], fez, fey, v0, v1, v2, boxhalfsize); + axisTest_y02(e0[2], e0[0], fez, fex, v0, v1, v2, boxhalfsize); + axisTest_z12(e0[1], e0[0], fey, fex, v0, v1, v2, boxhalfsize); - // Compute barycentric coordinates (v, w) of projection point - double denom = (acDOTac * abDOTab - abDOTac *abDOTac); - if (Math.abs(denom) < 1.0e-30) { - return new Vector3D(-100,-100,-100); // Triangle is degenerate - } + fex = Math.abs(e1[0]); + fey = Math.abs(e1[1]); + fez = Math.abs(e1[2]); - double w = (acDOTac * abDOTap - abDOTac * acDOTap)/denom; //coordinate towards b from a - double v = (abDOTab * acDOTap - abDOTac * abDOTap)/denom; //coordinate towards c from a + axisTest_x01(e1[2], e1[1], fez, fey, v0, v1, v2, boxhalfsize); + axisTest_y02(e1[2], e1[0], fez, fex, v0, v1, v2, boxhalfsize); + axisTest_z0(e1[1], e1[0], fey, fex, v0, v1, v2, boxhalfsize); - // Check barycentric coordinates - if ((v >= 0) && (w >= 0) && (v + w <= 1)) { - // Nearest orthogonal projection point is in triangle - return projection; - } + fex = Math.abs(e2[0]); + fey = Math.abs(e2[1]); + fez = Math.abs(e2[2]); - //region nearest point is on edge - if(w <= 0 && v > w){ - return t[0].add(ab.scalarMultiply(v)); - } + axisTest_x2(e2[2], e2[1], fez, fey, v0, v1, v2, boxhalfsize); + axisTest_y1(e2[2], e2[0], fez, fex, v0, v1, v2, boxhalfsize); + axisTest_z12(e2[1], e2[0], fey, fex, v0, v1, v2, boxhalfsize); - if(v <= 0 && w > v){ - return t[0].add(ac.scalarMultiply(w)); + double min = findMin(v0[0], v1[0], v2[0]); + double max = findMax(v0[0], v1[0], v2[0]); + if ((min > boxhalfsize[0]) || (max < -boxhalfsize[0])) { + return 0; } - - if(v + w > 1){ - final double scalarValue = bcDOTbp/bc.getNormSq(); - return t[1].add(bc.scalarMultiply(scalarValue)); + min = findMin(v0[1], v1[1], v2[1]); + max = findMax(v0[1], v1[1], v2[1]); + if ((min > boxhalfsize[1]) || (max < -boxhalfsize[1])) { + return 0; } - //endregion - - if (v <=0 && w <= 0){ //this should be redundant, but for some reason isn't - return t[0]; + min = findMin(v0[2], v1[2], v2[2]); + max = findMax(v0[2], v1[2], v2[2]); + if ((min > boxhalfsize[2]) || (max < -boxhalfsize[2])) { + return 0; } - return new Vector3D(-100,-100,-100); + cross(normal, e0, e1); + if (planeBoxOverlap(normal, v0, boxhalfsize) != 1) { + return 0; + } + return 1; } -} + +} \ No newline at end of file diff --git a/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelize3D.java b/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelize3D.java new file mode 100644 index 000000000..aab287595 --- /dev/null +++ b/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelize3D.java @@ -0,0 +1,230 @@ +/* + * #%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.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. + *

+ * + * @author Andrew McCall (University at Buffalo) + */ +@Plugin(type = Ops.Geometric.Voxelize.class) +public class DefaultVoxelize3D > extends AbstractUnaryHybridCF + implements Ops.Geometric.Voxelize { + + @Parameter + private OpService ops; + + @Parameter(type = ItemIO.INPUT, required = false) + 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); + 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) <= wallThickness/2) { + ra.setPositionAndGet(iterator.positionAsPoint()).set(true); + } + } + } + }); + + } + + @Override + public O createOutput(Mesh input){ + + float[] bounds = Meshes.boundingBox(input); + long[] outputInterval = new long[3]; + for (int i = 0; i < 3; i++) { + outputInterval[i] = (long)Math.ceil(bounds[i+3]+(2*wallThickness)-bounds[i]); + } + return (O) ops.create().img(new FinalInterval(outputInterval), 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 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 tPoint = nearestPointInTriangle3D(p, t); + return p.distance(tPoint); + } + + private Vector3D nearestPointInTriangle3D(Vector3D p, Vector3D[] t) { + /* + 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 Obtain projection of point p onto plane of triangle + // 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 + } + + // 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 5d5b98875..febc26bdd 100644 --- a/src/test/java/net/imagej/ops/geom/MeshFeatureTests.java +++ b/src/test/java/net/imagej/ops/geom/MeshFeatureTests.java @@ -31,6 +31,7 @@ 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; @@ -38,11 +39,14 @@ import net.imagej.ops.Ops; import net.imagej.ops.features.AbstractFeatureTest; import net.imagej.ops.geom.geom3d.*; -import net.imglib2.RandomAccessibleInterval; -import net.imglib2.algorithm.neighborhood.DiamondShape; +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; @@ -190,14 +194,26 @@ public void verticesCountMesh() { } @Test - public void voxelization3D() { - /*Value of 184 here corresponds with: - RandomAccessibleInterval result = (RandomAccessibleInterval) ops.run(DefaultVoxelization3D.class, mesh, ROI); - ops.morphology().fillHoles(result, result, new DiamondShape(1)); - assertEquals(Ops.Geometric.Voxelization.NAME, ROI.size(), Regions.countTrue(result)); - */ - assertEquals(Ops.Geometric.Voxelization.NAME, 184, - Regions.countTrue((RandomAccessibleInterval) ops.run(DefaultVoxelization3D.class, mesh, ROI))); + public void voxelize3D(){ + final Img out = new ArrayImgFactory<>(new BitType()).create(getTestImage3D()); + ops.run(DefaultVoxelize3D.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); + } + @Test + public void voxelization3D() { + // https://github.com/imagej/imagej-ops/issues/422 } } From 1bcf0a6af83938bc3621cfb43b6390ef3ffb430a Mon Sep 17 00:00:00 2001 From: andmccall Date: Tue, 7 Jul 2026 14:00:19 -0400 Subject: [PATCH 09/13] Renamed Voxelize3D to EuclideanDistanceVoxelization3D, and implemented it as Ops.Geometric.Voxelization. Also, fixed createOutput to actually generate an appropriately sized output image. Lastly, added comments to test to make it clear why this test works. --- .../net/imagej/ops/geom/GeomNamespace.java | 28 ++++++++----------- .../geom/geom3d/DefaultVoxelization3D.java | 2 +- ...a => EuclideanDistanceVoxelization3D.java} | 16 ++++++----- .../net/imagej/ops/geom/MeshFeatureTests.java | 16 ++++++----- 4 files changed, 30 insertions(+), 32 deletions(-) rename src/main/java/net/imagej/ops/geom/geom3d/{DefaultVoxelize3D.java => EuclideanDistanceVoxelization3D.java} (94%) diff --git a/src/main/java/net/imagej/ops/geom/GeomNamespace.java b/src/main/java/net/imagej/ops/geom/GeomNamespace.java index 6b0f64a4c..029d092e8 100644 --- a/src/main/java/net/imagej/ops/geom/GeomNamespace.java +++ b/src/main/java/net/imagej/ops/geom/GeomNamespace.java @@ -37,7 +37,6 @@ import net.imagej.ops.OpMethod; import net.imagej.ops.Ops.Geometric.Voxelization; import net.imagej.ops.geom.geom3d.mesh.VertexInterpolator; -import net.imglib2.Interval; import net.imglib2.IterableInterval; import net.imglib2.RandomAccessibleInterval; import net.imglib2.RealLocalizable; @@ -678,12 +677,14 @@ public double[] vertexInterpolator(final int[] p1, final int[] p2, return result; } - @Deprecated - @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; } @@ -704,24 +705,17 @@ public RandomAccessibleInterval voxelization(final Mesh in, final int w return result; } - @OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelize3D.class) - public RandomAccessibleInterval voxelize(final Mesh in) { - final RandomAccessibleInterval result = - (RandomAccessibleInterval) ops().run(net.imagej.ops.geom.geom3d.DefaultVoxelize3D.class, in); - return result; - } - - @OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelize3D.class) - public RandomAccessibleInterval voxelize(final RandomAccessibleInterval out, final Mesh in) { + @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.DefaultVoxelize3D.class, out, in); + (RandomAccessibleInterval) ops().run(net.imagej.ops.geom.geom3d.EuclideanDistanceVoxelization3D.class, out, in); return result; } - @OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelize3D.class) - public RandomAccessibleInterval voxelize(final RandomAccessibleInterval out, final Mesh in, final double wallThickness) { + @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.DefaultVoxelize3D.class, out, in, wallThickness); + (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 78c360dae..efea071bb 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,7 @@ *

* * @author Kyle Harrington (University of Idaho) - * @deprecated Use {@link DefaultVoxelize3D} instead. + * @deprecated Use {@link EuclideanDistanceVoxelization3D} instead. */ @Plugin(type = Ops.Geometric.Voxelization.class) public class DefaultVoxelization3D extends AbstractUnaryFunctionOp> diff --git a/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelize3D.java b/src/main/java/net/imagej/ops/geom/geom3d/EuclideanDistanceVoxelization3D.java similarity index 94% rename from src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelize3D.java rename to src/main/java/net/imagej/ops/geom/geom3d/EuclideanDistanceVoxelization3D.java index aab287595..8e3b12a42 100644 --- a/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelize3D.java +++ b/src/main/java/net/imagej/ops/geom/geom3d/EuclideanDistanceVoxelization3D.java @@ -43,6 +43,7 @@ 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; @@ -54,9 +55,9 @@ * * @author Andrew McCall (University at Buffalo) */ -@Plugin(type = Ops.Geometric.Voxelize.class) -public class DefaultVoxelize3D > extends AbstractUnaryHybridCF - implements Ops.Geometric.Voxelize { +@Plugin(type = Ops.Geometric.Voxelization.class, priority = Priority.HIGH) +public class EuclideanDistanceVoxelization3D> extends AbstractUnaryHybridCF + implements Ops.Geometric.Voxelization { @Parameter private OpService ops; @@ -82,18 +83,19 @@ public void compute(final Mesh input, final O output){ } } }); - } @Override public O createOutput(Mesh input){ float[] bounds = Meshes.boundingBox(input); - long[] outputInterval = new long[3]; + long[] min = new long[3]; + long[] max = new long[3]; for (int i = 0; i < 3; i++) { - outputInterval[i] = (long)Math.ceil(bounds[i+3]+(2*wallThickness)-bounds[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(outputInterval), new BitType()); + return (O) ops.create().img(new FinalInterval(min,max), new BitType()); } private Vector3D[] triangleToVector3DTriangle(Triangle t){ diff --git a/src/test/java/net/imagej/ops/geom/MeshFeatureTests.java b/src/test/java/net/imagej/ops/geom/MeshFeatureTests.java index febc26bdd..7f4b3fb34 100644 --- a/src/test/java/net/imagej/ops/geom/MeshFeatureTests.java +++ b/src/test/java/net/imagej/ops/geom/MeshFeatureTests.java @@ -194,9 +194,16 @@ public void verticesCountMesh() { } @Test - public void voxelize3D(){ + 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(DefaultVoxelize3D.class,out, mesh, 1.0); + 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(); @@ -211,9 +218,4 @@ public void voxelize3D(){ matches = false; assertTrue(matches); } - - @Test - public void voxelization3D() { - // https://github.com/imagej/imagej-ops/issues/422 - } } From bb905de2c8323fb926640e423a8cabf3aee7712b Mon Sep 17 00:00:00 2001 From: andmccall Date: Wed, 8 Jul 2026 09:51:53 -0400 Subject: [PATCH 10/13] Moved normal calculation to be once per triangle instead of caluculating it for every pixel analyzed. --- .../EuclideanDistanceVoxelization3D.java | 42 ++++++++++++------- 1 file changed, 26 insertions(+), 16 deletions(-) diff --git a/src/main/java/net/imagej/ops/geom/geom3d/EuclideanDistanceVoxelization3D.java b/src/main/java/net/imagej/ops/geom/geom3d/EuclideanDistanceVoxelization3D.java index 8e3b12a42..485773186 100644 --- a/src/main/java/net/imagej/ops/geom/geom3d/EuclideanDistanceVoxelization3D.java +++ b/src/main/java/net/imagej/ops/geom/geom3d/EuclideanDistanceVoxelization3D.java @@ -71,13 +71,17 @@ public void compute(final Mesh input, final O output){ 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) <= wallThickness/2) { + if (pointToTriangleDist(new Vector3D(iterator.getDoublePosition(0), iterator.getDoublePosition(1), iterator.getDoublePosition(2)), vector3DTriangle, normal) <= wallThickness/2) { ra.setPositionAndGet(iterator.positionAsPoint()).set(true); } } @@ -108,6 +112,23 @@ private Vector3D[] triangleToVector3DTriangle(Triangle t){ 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]; @@ -122,12 +143,12 @@ private Interval boundingBox(Vector3D[] t){ return new FinalInterval(min, max); } - private double pointToTriangleDist(Vector3D p, Vector3D[] t){ - Vector3D tPoint = nearestPointInTriangle3D(p, t); + 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) { + 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.' @@ -135,18 +156,7 @@ private Vector3D nearestPointInTriangle3D(Vector3D p, Vector3D[] t) { Vector3D ab = t[1].subtract(t[0]); Vector3D ac = t[2].subtract(t[0]); - //region Obtain projection of point p onto plane of triangle - // 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 - } - + //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)); From 09bbaaa0e7b3ca00efd9b789279bd32a599794ac Mon Sep 17 00:00:00 2001 From: andmccall Date: Thu, 16 Jul 2026 15:47:18 -0400 Subject: [PATCH 11/13] Added more description in comments. --- .../ops/geom/geom3d/DefaultVoxelization3D.java | 1 + .../geom3d/EuclideanDistanceVoxelization3D.java | 17 +++++++++++++---- 2 files changed, 14 insertions(+), 4 deletions(-) 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 efea071bb..eba0f1db2 100644 --- a/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java +++ b/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java @@ -60,6 +60,7 @@ * @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 index 485773186..c5600f65c 100644 --- a/src/main/java/net/imagej/ops/geom/geom3d/EuclideanDistanceVoxelization3D.java +++ b/src/main/java/net/imagej/ops/geom/geom3d/EuclideanDistanceVoxelization3D.java @@ -49,10 +49,18 @@ /** *

- * This is a voxelizer that produces a binary image with values set to true along - * the surface of the mesh. + * This is a voxelizer that produces a binary image with values set to true along the surface of the mesh. + * Calculates the Euclidean distance in pixel 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, 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 net.imagej.ops.Ops.Geometric.MarchingCubes} 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) @@ -62,7 +70,8 @@ public class EuclideanDistanceVoxelization3D Date: Fri, 17 Jul 2026 08:57:25 -0400 Subject: [PATCH 12/13] Fixed some grammar errors. --- .../geom/geom3d/EuclideanDistanceVoxelization3D.java | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/main/java/net/imagej/ops/geom/geom3d/EuclideanDistanceVoxelization3D.java b/src/main/java/net/imagej/ops/geom/geom3d/EuclideanDistanceVoxelization3D.java index c5600f65c..7a9409dfa 100644 --- a/src/main/java/net/imagej/ops/geom/geom3d/EuclideanDistanceVoxelization3D.java +++ b/src/main/java/net/imagej/ops/geom/geom3d/EuclideanDistanceVoxelization3D.java @@ -50,12 +50,12 @@ /** *

* This is a voxelizer that produces a binary image with values set to true along the surface of the mesh. - * Calculates the Euclidean distance in pixel between each pixel in the output image and the nearest point in 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, and subsequently performing a + * 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 net.imagej.ops.Ops.Geometric.MarchingCubes} for any object in the original data that is successfully converted + * {@link DefaultMarchingCubes} for any object in the original data that is successfully converted * to a mesh via marching cubes. *

* @@ -70,8 +70,8 @@ public class EuclideanDistanceVoxelization3D Date: Fri, 17 Jul 2026 09:06:17 -0400 Subject: [PATCH 13/13] Attempt to fix whitespace changes from IntelliJ idea --- .../geom/geom3d/DefaultVoxelization3D.java | 28 +++++++++---------- 1 file changed, 14 insertions(+), 14 deletions(-) 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 eba0f1db2..73fd13e8b 100644 --- a/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java +++ b/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java @@ -6,13 +6,13 @@ * %% * 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 @@ -56,7 +56,7 @@ * Thanks to Tomas Möller for sharing his public domain code: * http://fileadmin.cs.lth.se/cs/personal/tomas_akenine-moller/code/tribox.txt *

- * + * * @author Kyle Harrington (University of Idaho) * @deprecated Use {@link EuclideanDistanceVoxelization3D} instead. */ @@ -131,9 +131,9 @@ public RandomAccessibleInterval calculate(Mesh input) { Math.max(Math.max(v1.getZ(), v2.getZ()), v3.getZ()) - minPoint.getDoublePosition(2) }; RandomAccess ra = outImg.randomAccess();// Should use the - // interval - // implementation - // for speed + // interval + // implementation + // for speed long[] indices = new long[3]; for (indices[0] = (long) Math.floor(minSubBoundary[0] / stepSizes[0]); indices[0] < Math @@ -144,7 +144,7 @@ public RandomAccessibleInterval calculate(Mesh input) { .floor(maxSubBoundary[2] / stepSizes[2]); indices[2]++) { ra.setPosition(indices); if (!ra.get().get())// Don't check if voxel is already - // filled + // filled { double[] voxelCenter = new double[3]; @@ -198,7 +198,7 @@ private int planeBoxOverlap(double[] normalArray, double[] vertArray, double[] m } private int axisTest_x01(double e0, double e02, double fez, double fey, double[] v0, double[] v1, double[] v2, - double[] boxhalfsize) { + double[] boxhalfsize) { double p0 = e0 * v0[1] - e02 * v0[2]; double p2 = e0 * v2[1] - e02 * v2[2]; double max; @@ -219,7 +219,7 @@ private int axisTest_x01(double e0, double e02, double fez, double fey, double[] } private int axisTest_x2(double a, double b, double fa, double fb, double[] v0, double[] v1, double[] v2, - double[] boxhalfsize) { + double[] boxhalfsize) { double p0 = a * v0[1] - b * v0[2]; double p1 = a * v1[1] - b * v1[2]; double max; @@ -240,7 +240,7 @@ private int axisTest_x2(double a, double b, double fa, double fb, double[] v0, d } private int axisTest_y02(double a, double b, double fa, double fb, double[] v0, double[] v1, double[] v2, - double[] boxhalfsize) { + double[] boxhalfsize) { double p0 = -a * v0[0] + b * v0[2]; double p2 = -a * v2[0] + b * v2[2]; double max; @@ -261,7 +261,7 @@ private int axisTest_y02(double a, double b, double fa, double fb, double[] v0, } private int axisTest_y1(double a, double b, double fa, double fb, double[] v0, double[] v1, double[] v2, - double[] boxhalfsize) { + double[] boxhalfsize) { double p0 = -a * v0[0] + b * v0[2]; double p1 = -a * v1[0] + b * v1[2]; double max; @@ -282,7 +282,7 @@ private int axisTest_y1(double a, double b, double fa, double fb, double[] v0, d } private int axisTest_z12(double a, double b, double fa, double fb, double[] v0, double[] v1, double[] v2, - double[] boxhalfsize) { + double[] boxhalfsize) { double p1 = a * v1[0] - b * v1[1]; double p2 = a * v2[0] - b * v2[1]; double max; @@ -303,7 +303,7 @@ private int axisTest_z12(double a, double b, double fa, double fb, double[] v0, } private int axisTest_z0(double a, double b, double fa, double fb, double[] v0, double[] v1, double[] v2, - double[] boxhalfsize) { + double[] boxhalfsize) { double p0 = a * v0[0] - b * v0[1]; double p1 = a * v1[0] - b * v1[1]; double max; @@ -403,4 +403,4 @@ private int triBoxOverlap(double[] boxcenter, double[] boxhalfsize, Vector3D pf1 return 1; } -} \ No newline at end of file +}