Overhaul of Voxelization ops with new algorithm. Also included test

[andmccall]

Wanted to address a major issue with the voxelization Op not working appropriately. The old algorithm was returning every pixel of every triangle's entire bounding box as true. This new method only sets the pixels that the surface goes through to true.

In general, this method should now be more symmetric with marching cubes, so that after:

mesh = ops.geom().marchingCubes(image);
result = ops.geom().voxelization(mesh, image);
ops.morphology().fillHoles(result, result, new DiamondShape(1));

Image and result should match in the majority of cases. The notable exception would be images with holes themselves, generating a surface encompassed entirely by a surface in marching cubes. Voxelization would generate both surfaces, but I believe all encompassed holes would be filled after the fillHoles call, though I haven't tested this yet.

I did add some comments in the code for some optional changes that could be made. Originally I had an additional parameter for the wall thickness, but wasn't sure if this was really beneficial. I also wasn't sure if I should be using another Op (fill holes) during the voxelization Op testing, so I left it out and just compared the result to a static number.

Edit: was just looking through and noticed I forgot to delete the LogService parameter, which I was using for troubleshooting an issue.

[andmccall]

Just discovered a bug, working on fixing it and should have a new commit soon.
Edit: done. Sorry for the changes after the pull request, didn't notice these issues in initial testing.

[andmccall]

@gselzer and @ctrueden,

Think I have all the kinks worked out and it's ready for review. I did find that there were a couple @OpMethods I put in the GeomNamespace that I had to comment out for the JUnit tests to pass (Input Mismatch error), even though the function calls worked fine when I compiled them without tests. Since it sounds like everything is gonna move to SciJava Ops, I figured it wasn't worthwhile for me to figure out the issue as I'm not very familiar with JUnit testing.

Anyway, I'll for SciJava-ops and work on converting this to the SciJava-Ops format.

[gselzer]

@andmccall thanks so much for doing this work! It's great to see contributions!

I really love your code - I think it's a massive improvement over the existing functionality. Just some things I was thinking about as I was looking through.

(Also, for those who aren't aware, there was a bit of discussion here

[gselzer]

This comment has absolutely nothing to do with your work, since it's clear you're just copying the precedent, but man, I'm remembering how much I disliked these tests 😆 ; if the goal is to assert proper voxelization, it is of course not sufficient to count the number of trues - you'll also need to know where they are.

I would love a test like this for the eventual SciJava Ops version, but up to you whether you think we should do something like this now (probably shouldn't go here though):

Suggested change

	public void voxelization3D() {
	// https://github.com/imagej/imagej-ops/issues/422
	/*Value of 184 here corresponds with:
	RandomAccessibleInterval<BitType> result = (RandomAccessibleInterval<BitType>) 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<BitType>) ops.run(DefaultVoxelization3D.class, mesh, ROI)));
	}
	public void voxelization3D() {
	// Create a "square" mesh from [1, 1, 1] to [2, 2, 1]
	Mesh m = new NaiveDoubleMesh();
	m.vertices().add(1, 1, 1);
	m.vertices().add(1, 2, 1);
	m.vertices().add(2, 1, 1);
	m.vertices().add(2, 2, 1);
	m.triangles().add(0, 1, 2);
	m.triangles().add(2, 3, 1);
	// Voxelize it
	Interval interval = new FinalInterval(new long[] {0, 0, 0}, new long[] {3, 3, 3});
	RandomAccessibleInterval<BitType> voxelized = (RandomAccessibleInterval<BitType>) ops.run(DefaultVoxelization3D.class, m, interval);
	// Check the results
	Cursor<BitType> cursor = voxelized.cursor();
	while (cursor.hasNext()) {
	boolean value = cursor.next().get();
	long[] pos = cursor.positionAsLongArray();
	if (pos[0] < 1 || pos[0] > 2) {
	// x dimension outside the square
	assertFalse(value);
	}
	else if (pos[1] < 1 || pos[1] > 2) {
	// y dimension outside the square
	assertFalse(value);
	}
	else if (pos[2] != 1) {
	// z dimension outside the square
	assertFalse(value);
	}
	else {
	// "In" the square
	assertTrue(value);
	}
	}
	}

[andmccall]

Hi @gselzer,

Another option that I just discovered yesterday is that we can compare the result of voxelization(mesh, originalBinary) to the result of net.imglib2.roi.boundary.Boundary(originalBinary, new DiamondShape(1)). It turns out that this produces identical results for nearly all meshes. This would allow the continued use of the slightly more complex ROI image with the pre-calculated corresponding mesh that was present in ImageJ ops (if this carried over to SciJava Ops).

[andmccall]

Correction: not DiamondShape(1), but Boundary.StructuringElement.FOUR_CONNECTED.

[gselzer]

These do not work because optional Op signatures are only "valid" if parameters are left off right to left. In this case, valid signatures would be:

• voxelization(final Mesh in, final Interval dimensions, final boolean scaleMeshToDimensions, final double wallThickness)
• voxelization(final Mesh in, final Interval dimensions, final boolean scaleMeshToDimensions)
• voxelization(final Mesh in, final Interval dimensions)
• voxelization(final Mesh in)

Of course these two signatures are not included. The context for this right-to-left dropping is to remove ambiguity when you have two optional parameters that are of the same type e.g. for some

voxelization(final Mesh in, final double a, final Double b)

you could not have both

voxelization(final Mesh in, final double a)
voxelization(final Mesh in, final double b)

However, your use case does not have this problem. You could consider skipping the type check (unclear if this is what you meant by skip tests) - do you think this is an acceptable use case of that parameter @ctrueden?

Suggested change

	//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<BitType> voxelization(final Mesh in, double wallThickness) {
	// final RandomAccessibleInterval<BitType> result =
	// (RandomAccessibleInterval<BitType>) 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<BitType> voxelization(final Mesh in, final Interval dimensions, double wallThickness) {
	// final RandomAccessibleInterval<BitType> result =
	// (RandomAccessibleInterval<BitType>) ops().run(net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class, in, dimensions, false, wallThickness);
	// return result;
	// }
	@OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class, skipTypeCheck = true)
	public RandomAccessibleInterval<BitType> voxelization(final Mesh in, double wallThickness) {
	final RandomAccessibleInterval<BitType> result =
	(RandomAccessibleInterval<BitType>) ops().run(net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class, in, null, false, wallThickness);
	return result;
	}
	@OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class, skipTypeCheck = true)
	public RandomAccessibleInterval<BitType> voxelization(final Mesh in, final Interval dimensions, double wallThickness) {
	final RandomAccessibleInterval<BitType> result =
	(RandomAccessibleInterval<BitType>) ops().run(net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class, in, dimensions, false, wallThickness);
	return result;
	}

[gselzer]

Would be great to either get these working or delete them from the code

[andmccall]

I'm fine with it either way. Wasn't strongly invested in having these signatures, I just wasn't really sure what the issue was, or how to do the SkipTypeCheck.

[gselzer]

I also do not feel strongly - I just want to make sure that you get the code you want in, and it's easier to add these things in later than to remove them now. So if you don't have a use let's just remove them?

[andmccall]

Sounds good to me.

[gselzer]

This is a breaking change 🛠️

Since your Op has no overlap in parameter types with the existing Op, we could consider offering the two Ops side-by-side? Could also consider some sort of deprecation for the old Op...

[andmccall]

This would change the algorithm used, so I'm not sure you would want to do this, but you could do:

final RandomAccessibleInterval<BitType> result = (RandomAccessibleInterval<BitType>) ops().run( Voxelization.class, in, new FinalInterval(width, height, depth));

This means that the output from any old code would change, but the function call wouldn't.

I imagine that is not the preferred option, but it doesn't really matter to me, can definitely keep the old algorithm in there in whatever manner you deem best. I simply forgot about the whole "can't make breaking changes" thing, never done any library coding before.

[gselzer]

This would change the algorithm used, so I'm not sure you would want to do this, but you could do:

final RandomAccessibleInterval<BitType> result = (RandomAccessibleInterval<BitType>) ops().run( Voxelization.class, in, new FinalInterval(width, height, depth));

This means that the output from any old code would change, but the function call wouldn't.

Yes, this option is certainly available, however I think I'd prefer avoiding breakages if possible. I wouldn't say that the previous implementation was a bug, it's just a different algorithm (although certainly a strange one)

I simply forgot about the whole "can't make breaking changes" thing, never done any library coding before.

It happens to the best of us, don't worry 🙂

[gselzer]

Let's think about the tasks a user might need to accomplish using this Op. I can see the user wanting to:

• Obtain a voxelized image across the entire mesh bounding box
• Obtain a voxelized image across some custom interval (e.g. if the mesh is massive maybe they don't want to voxelize the whole thing).

Unless you can think of additional tasks, I think that this collection of goals would map nicely to a UnaryHybridCF Op, instead of a UnaryFunctionOp. This would allow users to call the Op with any preallocated output image defined on whatever interval they want, which would make the second goal easy. Additionally, if they called the Op as a Function, then you could use the logic you have for the bounding box to create a new image over the entire mesh interval.

What do you think about restructuring the Op to be a UnaryHybridCF? Do you think you could get rid of the Interval parameter that way?

(As an aside, in SciJava Ops hybrid Ops aren't a thing, we'd just write two separate Ops for that. But we'll cross that bridge later 😄)

[gselzer]

I would also wonder whether you could get rid of the scale and offset state within the Op if you did this. One key aspect of SciJava Ops is that Ops are stateless, so it'd be great to remove any state if possible if we want to eventually port this work over there!

[andmccall]

I'm not sure as to the difference between these 😅, or how it would work without the interval. I'll have to take a closer look, but will be quite busy until next week at the earliest. Just wanted to reply cause my work schedule tends to slingshot between totally free and a completely packed schedule, and I think it's swinging back to packed for a little while.

[gselzer]

Just wanted to reply cause my work schedule tends to slingshot between totally free and a completely packed schedule, and I think it's swinging back to packed for a little while.

I definitely understand, there is no rush!

[gselzer]

You said here:

Additionally, if the mesh wasn't previously an image the vertex coordinates in the mesh could be negative, which doesn't work very well with converting to images. Also, it just wasn't that hard to program.

If the goal is to describe the voxels a mesh passes through, setting this parameter to true would almost certainly obfuscate that information, since the scale and the offset are not reflected in the result.

What do you think about changing the type of the Op as I mentioned above? If that were the case, we would no longer have to worry about the scale and offset being necessary to transfer meaning between the mesh and the voxelized image...

[gselzer]

I really appreciate the comments here! This kind of context will ease future maintenance.

Would you be willing to add similar comments to the other functions as well?

[andmccall]

Absolutely. I love adding comments like this to math, as it takes me forever to wrap my head around this stuff and I get lost without the comments.

I did have a related question about SciJava, as I just want to clarify that I'm understanding everything correctly. It seems like a lot of the SciJava Ops are intentionally kept to one function, and these functions are chained together into a computer. So would the ideal be to break apart this Voxelization class into functions: ProjectionOntoPlane?, FindNearestPointInTriangle, GetDistanceToTriangle, and maybe others; then compile all these individual functions into a Voxelization computer? Or would it be better to do this as an Inplace? I'm not really sure if the default is to create more functions, or to create functions only if there's reason to believe they'd be useful to others (and I have no idea if they would be).

[gselzer]

I love adding comments like this to math, as it takes me forever to wrap my head around this stuff and I get lost without the comments.

Me too 😁

It seems like a lot of the SciJava Ops are intentionally kept to one function, and these functions are chained together into a computer.

There's definitely some terminology here. Functions, Computers, and Inplaces are all first-class algorithm types in SciJava Ops (and here in ImageJ Ops), but they all handle their outputs differently:

• Functions create a new object every time.
• Computers fill up an output buffer given to them as a parameter.
• Inplaces overwrite one of their inputs with the output data (this is uncommon in our libraries but exists out in the wild).

So would the ideal be to break apart this Voxelization class into functions: ProjectionOntoPlane?, FindNearestPointInTriangle, GetDistanceToTriangle, and maybe others; then compile all these individual functions into a Voxelization computer?

I don't think it's worth making these functions into Ops right now - we could always do that later if there's a use for others!