Add support for expressing Mesh as RealMask

[ctrueden]

The imglib2-roi library has a RealMask API for working with ROIs.

If we have a Mesh that is closed and convex, we can express this as a 3D RealMask by triangulating the convex hull and then testing points for inclusion by walking the triangles.

[ctrueden]

Here is a proof-of-concept implementation that uses VisAD. You can add the org.scijava:visad:2.0.0-scijava-1 artifact to your Maven dependencies.

import java.util.function.BiConsumer;

import net.imagej.Dataset;
import net.imagej.ImageJ;
import net.imagej.mesh.Mesh;
import net.imagej.mesh.Vertex;
import net.imglib2.Localizable;
import net.imglib2.RandomAccessible;
import net.imglib2.RandomAccessibleInterval;
import net.imglib2.RealLocalizable;
import net.imglib2.RealPoint;
import net.imglib2.RealRandomAccessible;
import net.imglib2.position.FunctionRandomAccessible;
import net.imglib2.roi.Masks;
import net.imglib2.roi.RealMask;
import net.imglib2.type.logic.BitType;
import net.imglib2.type.logic.BoolType;
import net.imglib2.view.Views;

import visad.Irregular3DSet;
import visad.RealTupleType;
import visad.RealType;
import visad.VisADException;

public class MeshToRealMask {

	public static void main(String[] args) throws Throwable {
		ImageJ ij = new ImageJ();

		// Generate a silly binary image.
		BiConsumer<Localizable, BitType> function = (l, t) -> {
			long x = l.getLongPosition(0);
			long y = l.getLongPosition(1);
			long z = l.getLongPosition(2);
			t.set(
				x + y > 50 && x + y < 80 ||
				x + z > 70 && x + z < 110
			);
		};
		FunctionRandomAccessible<BitType> fra = //
			new FunctionRandomAccessible<>(3, function, () -> new BitType());
		long[] min = {0, 0, 0};
		long[] max = {256, 256, 64};
		RandomAccessibleInterval<BitType> image = Views.interval(fra, min, max);
		Dataset dataset = ij.dataset().create(image);
		System.out.println("Dataset created: " + dataset);

		// Compute mesh.
		Mesh mesh = ij.op().geom().marchingCubes(image);
		System.out.println("Mesh computed");

		// Compute convex hull.
		Mesh hull = (Mesh) ij.op().geom().convexHull(mesh).get(0);
		System.out.println("Hull computed");

		// Compute Delaunay triangulation.
		long vertexCount = hull.vertices().size();
		if (vertexCount > Integer.MAX_VALUE) throw new RuntimeException("Hull too large");
		float[][] samples = new float[3][(int) vertexCount];
		int i = 0;
		for (Vertex vertex : hull.vertices()) {
			samples[0][i] = vertex.xf();
			samples[1][i] = vertex.yf();
			samples[2][i] = vertex.zf();
			i++;
		}
		RealType xType = RealType.getRealType("X");
		RealType yType = RealType.getRealType("Y");
		RealType zType = RealType.getRealType("Z");
		RealType[] xyz = {xType, yType, zType};
		RealTupleType xyzType = new RealTupleType(xyz);
		Irregular3DSet set = new Irregular3DSet(xyzType, samples);
		System.out.println("Triangulation computed");

		// Wrap as a RealMask.
		RealMask roi = new RealMask() {
			private ThreadLocal<float[][]> value = new ThreadLocal<float[][]>() {
				@Override
				public float[][] initialValue() {
					return new float[3][1];
				}
			};

			@Override
			public boolean test(RealLocalizable t) {
				float[][] v = value.get();
				try {
					v[0][0] = t.getFloatPosition(0);
					v[1][0] = t.getFloatPosition(1);
					v[2][0] = t.getFloatPosition(2);
					int[] tri = set.valueToTri(v);
					return tri[0] >= 0;
				}
				catch (VisADException exc) {
					throw new RuntimeException(exc);
				}
			}

			@Override
			public int numDimensions() { return 3; }
		};
		System.out.println("ROI created: " + roi);

		class Tester {
			private RealPoint point = new RealPoint(3);
			public void inside(double... p) {
				point.setPosition(p[0], 0);
				point.setPosition(p[1], 1);
				point.setPosition(p[2], 2);
				System.out.println(point + " -> " + roi.test(point));
			}
		}
		Tester t = new Tester();
		System.out.println("Testing some points:");
		t.inside(1, 2, 3);
		t.inside(50, 45, 40);
		t.inside(100, 55, 10);

		// Display the results.
		ij.ui().showUI();
//		ij.get(ROIService.class).add(roi, dataset);
		ij.ui().show("Dataset", dataset);

		// Use ImageJ1 to display the ROI. :-(
//		ImagePlus imp = IJ.getImage();
//		Roi ij1Roi = ij.convert().convert(roi, Roi.class);
//		System.out.println("IJ1 roi = " + ij1Roi);
//		imp.setRoi(ij1Roi);

		// Use BDV to display the ROI. :-/
//		RealRandomAccessible<BoolType> maskRA = Masks.toRealRandomAccessible(roi);
//		BdvFunctions.show(maskRA, dataset, "ROI");

		// Wrap the ROI in a RandomAccessibleInterval and display that.
		RealRandomAccessible<BoolType> roiRA = Masks.toRealRandomAccessible(roi);
		RandomAccessible<BoolType> roiRaster = Views.raster(roiRA);
		RandomAccessibleInterval<BoolType> roiRAI = Views.interval(roiRaster, dataset);
		ij.ui().show("ROI", roiRAI);
	}
}