diff --git a/src/main/java/net/imagej/ops/geom/GeomNamespace.java b/src/main/java/net/imagej/ops/geom/GeomNamespace.java
index 6e0d449f0..ac4f1a6dc 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<BitType> voxelization(final Mesh in, final int width, final int height, final int depth ) {
-		final RandomAccessibleInterval<BitType> result = (RandomAccessibleInterval<BitType>) 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,47 @@ public double[] vertexInterpolator(final int[] p1, final int[] p2,
 	
 	@OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class)
 	public RandomAccessibleInterval<BitType> voxelization(final Mesh in) {
-		@SuppressWarnings("unchecked")
 		final RandomAccessibleInterval<BitType> result =
 			(RandomAccessibleInterval<BitType>) ops().run(net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class, in);
 		return result;
 	}
 
 	@OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class)
-	public RandomAccessibleInterval<BitType> voxelization(final Mesh in, final int width) {
-		@SuppressWarnings("unchecked")
+	public RandomAccessibleInterval<BitType> voxelization(final Mesh in, final Interval dimensions) {
+		final RandomAccessibleInterval<BitType> result =
+			(RandomAccessibleInterval<BitType>) ops().run(net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class, in, dimensions);
+		return result;
+	}
+
+	@OpMethod(op = net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class)
+	public RandomAccessibleInterval<BitType> voxelization(final Mesh in, final Interval dimensions, boolean scaleMeshToDimesions) {
 		final RandomAccessibleInterval<BitType> result =
-			(RandomAccessibleInterval<BitType>) ops().run(net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class, in, width);
+				(RandomAccessibleInterval<BitType>) ops().run(net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class, in, dimensions, scaleMeshToDimesions);
 		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<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)
-	public RandomAccessibleInterval<BitType> voxelization(final Mesh in, final int width, final int height) {
+	public RandomAccessibleInterval<BitType> voxelization(final Mesh in, final Interval dimensions, boolean scaleMeshToDimesions, double wallThickness) {
 		final RandomAccessibleInterval<BitType> result =
-			(RandomAccessibleInterval<BitType>) ops().run(net.imagej.ops.geom.geom3d.DefaultVoxelization3D.class, in, width, height);
+				(RandomAccessibleInterval<BitType>) 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 8837875db..e1c1250b2 100644
--- a/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java
+++ b/src/main/java/net/imagej/ops/geom/geom3d/DefaultVoxelization3D.java
@@ -29,19 +29,17 @@
 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.plugin.Parameter;
@@ -49,356 +47,196 @@
 
 /**
  * <p>
- * 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.
  * </p>
- * <p>
- * Thanks to Tomas Möller for sharing his public domain code:
- * http://fileadmin.cs.lth.se/cs/personal/tomas_akenine-moller/code/tribox.txt
- * </p>
  * 
- * @author Kyle Harrington (University of Idaho)
+ * @author Andrew McCall (University at Buffalo)
  */
 @Plugin(type = Ops.Geometric.Voxelization.class)
 public class DefaultVoxelization3D extends AbstractUnaryFunctionOp<Mesh, RandomAccessibleInterval<BitType>>
 		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(type = ItemIO.INPUT, required = false)
+	private double wallThickness = 1.0;
+
+	private final long[] offset = new long[] {0,0,0};
+	private double scale = 1.0;
+
 	@Override
 	public RandomAccessibleInterval<BitType> calculate(Mesh input) {
 
-		Img<BitType> 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]+(2*wallThickness)-bounds[i]);
+			}
+			dimensions = new FinalInterval(outputInterval);
+			setScale(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<BitType> 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<BitType> 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<BitType> 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) {
+							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(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(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]);
 		}
-		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];
-	}
 
-	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]);
-	}
+		//    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 int triBoxOverlap(double[] boxcenter, double[] boxhalfsize, Vector3D pf1, Vector3D pf2, Vector3D pf3) {
-		double[] vert1 = pf1.toArray();
-		double[] vert2 = pf2.toArray();
-		double[] vert3 = pf3.toArray();
+		Vector3D ap = projection.subtract(t[0]);
 
-		double[] v0 = new double[3];
-		double[] v1 = new double[3];
-		double[] v2 = new double[3];
+		//region nearest point is corners
+		final double abDOTap = ab.dotProduct(ap);
+		final double acDOTap = ac.dotProduct(ap);
 
-		double[] normal = new double[3];
-		double[] e0 = new double[3];
-		double[] e1 = new double[3];
-		double[] e2 = new double[3];
+		if (abDOTap <= 0d && acDOTap <= 0d) return t[0];
 
-		sub(v0, vert1, boxcenter);
-		sub(v1, vert2, boxcenter);
-		sub(v2, vert3, boxcenter);
+		final Vector3D bc = t[2].subtract(t[1]);
+		final Vector3D bp = projection.subtract(t[1]);
 
-		sub(e0, v1, v0);
-		sub(e1, v2, v1);
-		sub(e2, v0, v2);
+		final double baDOTbp = ab.negate().dotProduct(bp);
+		final double bcDOTbp = bc.dotProduct(bp);
+		if (baDOTbp <= 0d && bcDOTbp <= 0d) return t[1];
 
-		double fex = Math.abs(e0[0]);
-		double fey = Math.abs(e0[1]);
-		double fez = Math.abs(e0[2]);
 
-		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);
+		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
 
-		fex = Math.abs(e1[0]);
-		fey = Math.abs(e1[1]);
-		fez = Math.abs(e1[2]);
+		double acDOTac = ac.dotProduct(ac);
+		double abDOTac =  ab.dotProduct(ac);
+		double abDOTab = ab.dotProduct(ab);
 
-		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);
+		// 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(e2[0]);
-		fey = Math.abs(e2[1]);
-		fez = Math.abs(e2[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_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<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)));
+
 	}
 }