diff --git a/src/main/java/net/imagej/mesh/alg/TaubinSmoothing.java b/src/main/java/net/imagej/mesh/alg/TaubinSmoothing.java new file mode 100644 index 0000000..8b51095 --- /dev/null +++ b/src/main/java/net/imagej/mesh/alg/TaubinSmoothing.java @@ -0,0 +1,271 @@ +package net.imagej.mesh.alg; + +import java.util.Arrays; + +import net.imagej.mesh.Mesh; +import net.imagej.mesh.Meshes; +import net.imagej.mesh.Triangles; +import net.imagej.mesh.Vertices; +import net.imagej.mesh.nio.BufferMesh; + +/** + * Taubin's mesh smoothing algorithm. + *

+ * Adapted by the Javascript code of mykolalysenko, MIT license. + * https://github.com/mikolalysenko/taubin-smooth/blob/master/smooth.js + * + * @author Jean-Yves Tinevez + * @see Taubin, G. “Curve and + * Surface Smoothing without Shrinkage.” In Proceedings of IEEE + * International Conference on Computer Vision, 852–57, 1995. + * + */ +public class TaubinSmoothing { + + /** + * Smooth the specified mesh using the Taubin smoothing algorithm, with default + * parameters. + * + * @param mesh the mesh to smooth. + * @return a new smoothed mesh. + */ + public static final BufferMesh smooth(final Mesh mesh) { + return smooth(mesh, 10, 0.1); + } + + /** + * Smooth the specified mesh using the Taubin smoothing algorithm, using + * cotangent weights. + * + * @param mesh the mesh to smooth. + * @param iters the number of iterations to apply. Typical values: 10. + * @param passBand the spatial frequency to use for smoothing. For instance 0.1. + * @return a new smoothed mesh. + */ + public static final BufferMesh smooth(final Mesh mesh, final int iters, final double passBand) { + final double A = -1.; + final double B = passBand; + final double C = 2.; + + final double discr = Math.sqrt(B * B - 4. * A * C); + final double r0 = (-B + discr) / (2. * A * C); + final double r1 = (-B - discr) / (2. * A * C); + + final double lambda = Math.max(r0, r1); + final double mu = Math.min(r0, r1); + + return smooth(mesh, iters, lambda, mu, TaubinWeightType.COTANGENT); + } + + /** + * Smooth the specified mesh using the Taubin smoothing algorithm. + * + * @param mesh the mesh to smooth. + * @param iters the number of iterations to apply. Typical values: 10. + * @param lambda the value of the lambda step in Taubin smoothing. Positive, + * typical values are around 0.5; + * @param mu the value of the mu step in Taubin smoothing. Negative, + * typical values are around -0.5. + * @param weightType the type of weights to use when summing contribution from + * neighbor edges. + * @return a new smoothed mesh. + */ + public static final BufferMesh smooth(final Mesh mesh, final int iters, final double lambda, final double mu, + final TaubinWeightType weightType) { + final int nvs = (int) mesh.vertices().size(); + final int nts = (int) mesh.triangles().size(); + final BufferMesh meshA = new BufferMesh(nvs, nts); + Meshes.copy(mesh, meshA); + final BufferMesh meshB = new BufferMesh(nvs, nts); + Meshes.copy(mesh, meshB); + + final double[] trace = new double[nvs]; + + switch (weightType) { + case COTANGENT: + for (int i = 0; i < iters; ++i) { + smoothStepCotangent(mesh.triangles(), meshA, meshB, trace, lambda); + smoothStepCotangent(mesh.triangles(), meshB, meshA, trace, mu); + } + break; + case NAIVE: + for (int i = 0; i < iters; ++i) { + smoothStepNaive(mesh.triangles(), meshA, meshB, trace, lambda); + smoothStepNaive(mesh.triangles(), meshB, meshA, trace, mu); + } + break; + default: + throw new IllegalArgumentException("Unhandled weight type: " + weightType); + } + + return meshA; + } + + private static void smoothStepCotangent(final Triangles triangles, final BufferMesh source, final BufferMesh target, + final double[] trace, final double weigth) { + + final int nvs = (int) source.vertices().size(); + final int nts = (int) source.triangles().size(); + + // Zero target. + for (int i = 0; i < nvs; i++) + target.vertices().set(i, 0., 0., 0.); + + // Zero trace. + Arrays.fill(trace, 0.); + + for (int i = 0; i < nts; ++i) { + + final int ia = (int) source.triangles().vertex0(i); + final double ax = source.vertices().x(ia); + final double ay = source.vertices().y(ia); + final double az = source.vertices().z(ia); + + final int ib = (int) source.triangles().vertex1(i); + final double bx = source.vertices().x(ib); + final double by = source.vertices().y(ib); + final double bz = source.vertices().z(ib); + + final int ic = (int) source.triangles().vertex2(i); + final double cx = source.vertices().x(ic); + final double cy = source.vertices().y(ic); + final double cz = source.vertices().z(ic); + + final double abx = ax - bx; + final double aby = ay - by; + final double abz = az - bz; + + final double bcx = bx - cx; + final double bcy = by - cy; + final double bcz = bz - cz; + + final double cax = cx - ax; + final double cay = cy - ay; + final double caz = cz - az; + + final double area = 0.5 * hypot(aby * caz - abz * cay, abz * cax - abx * caz, abx * cay - aby * cax); + if (area < 1e-8) + continue; + + final double w = -0.5 / area; + final double wa = w * (abx * cax + aby * cay + abz * caz); + final double wb = w * (bcx * abx + bcy * aby + bcz * abz); + final double wc = w * (cax * bcx + cay * bcy + caz * bcz); + + trace[ia] += wb + wc; + trace[ib] += wc + wa; + trace[ic] += wa + wb; + + accumulate(target.vertices(), ib, source.vertices(), ic, wa); + accumulate(target.vertices(), ic, source.vertices(), ib, wa); + accumulate(target.vertices(), ic, source.vertices(), ia, wb); + accumulate(target.vertices(), ia, source.vertices(), ic, wb); + accumulate(target.vertices(), ia, source.vertices(), ib, wc); + accumulate(target.vertices(), ib, source.vertices(), ia, wc); + } + + for (int i = 0; i < nvs; i++) { + final double ox = target.vertices().x(i); + final double oy = target.vertices().y(i); + final double oz = target.vertices().z(i); + final double ix = source.vertices().x(i); + final double iy = source.vertices().y(i); + final double iz = source.vertices().z(i); + final double tr = trace[i]; + + target.vertices().set(i, + ix + weigth * (ox / tr - ix), + iy + weigth * (oy / tr - iy), + iz + weigth * (oz / tr - iz)); + } + } + + private static void smoothStepNaive(final Triangles triangles, final BufferMesh source, final BufferMesh target, + final double[] trace, final double weigth) { + + final int nvs = (int) source.vertices().size(); + final int nts = (int) source.triangles().size(); + + // Zero target. + for (int i = 0; i < nvs; i++) + target.vertices().set(i, 0., 0., 0.); + + // Zero trace. + Arrays.fill(trace, 0.); + + for (int i = 0; i < nts; ++i) { + + final int ia = (int) source.triangles().vertex0(i); + final int ib = (int) source.triangles().vertex1(i); + final int ic = (int) source.triangles().vertex2(i); + + final double w = 0.5; + final double wa = w; + final double wb = w; + final double wc = w; + + trace[ia] += wb + wc; + trace[ib] += wc + wa; + trace[ic] += wa + wb; + + accumulate(target.vertices(), ib, source.vertices(), ic, wa); + accumulate(target.vertices(), ic, source.vertices(), ib, wa); + accumulate(target.vertices(), ic, source.vertices(), ia, wb); + accumulate(target.vertices(), ia, source.vertices(), ic, wb); + accumulate(target.vertices(), ia, source.vertices(), ib, wc); + accumulate(target.vertices(), ib, source.vertices(), ia, wc); + } + + for (int i = 0; i < nvs; i++) { + final double ox = target.vertices().x(i); + final double oy = target.vertices().y(i); + final double oz = target.vertices().z(i); + final double ix = source.vertices().x(i); + final double iy = source.vertices().y(i); + final double iz = source.vertices().z(i); + final double tr = trace[i]; + + target.vertices().set(i, + ix + weigth * (ox / tr - ix), + iy + weigth * (oy / tr - iy), + iz + weigth * (oz / tr - iz)); + } + } + + private final static void accumulate(final Vertices out, final int ok, final Vertices in, final int ik, + final double w) { + final double xo = out.x(ok); + final double yo = out.y(ok); + final double zo = out.z(ok); + final double xi = in.x(ik); + final double yi = in.y(ik); + final double zi = in.z(ik); + out.set(ok, + xo + xi * w, + yo + yi * w, + zo + zi * w); + } + + private static final double hypot(final double x, final double y, final double z) { + return Math.sqrt(x * x + y * y + z * z); + } + + /** + * The weight to use when accumlating contribution for neighborhood edges in + * Taubin smoothing. + */ + public enum TaubinWeightType { + + /** + * Use cotangent of the 2 triangles around an edge when adding the contribution + * of one edge. + */ + COTANGENT, + /** + * Use identical weight for all edges. + */ + NAIVE; + + } + +} diff --git a/src/test/java/net/imagej/mesh/alg/TaubinSmoothingDemo.java b/src/test/java/net/imagej/mesh/alg/TaubinSmoothingDemo.java new file mode 100644 index 0000000..faa0746 --- /dev/null +++ b/src/test/java/net/imagej/mesh/alg/TaubinSmoothingDemo.java @@ -0,0 +1,48 @@ +package net.imagej.mesh.alg; + +import java.io.IOException; +import java.util.Random; + +import net.imagej.mesh.Mesh; +import net.imagej.mesh.Meshes; +import net.imagej.mesh.alg.TaubinSmoothing.TaubinWeightType; +import net.imagej.mesh.io.ply.PLYMeshIO; +import net.imglib2.FinalInterval; +import net.imglib2.img.Img; +import net.imglib2.img.array.ArrayImgs; +import net.imglib2.type.numeric.integer.UnsignedByteType; +import net.imglib2.util.Intervals; +import net.imglib2.view.Views; + +public class TaubinSmoothingDemo { + + public static void main(final String[] args) throws IOException { + + // Make a cube + final Img img = ArrayImgs.unsignedBytes(100, 100, 100); + final FinalInterval cube = Intervals.createMinMax(25, 25, 25, 75, 75, 75); + Views.interval(img, cube).forEach(p -> p.setOne()); + + final Mesh m1 = Meshes.marchingCubes(img, 0.5); + final Mesh mesh = Meshes.removeDuplicateVertices(m1, 2); + + // Add noise. + final Random ran = new Random(45l); + for (int i = 0; i < mesh.vertices().size(); i++) { + final double x = mesh.vertices().x(i); + final double y = mesh.vertices().y(i); + final double z = mesh.vertices().z(i); + + mesh.vertices().set(i, + x + ran.nextDouble() - 0.5, + y + ran.nextDouble() - 0.5, + z + ran.nextDouble() - 0.5 ); + } + + new PLYMeshIO().save(mesh, "samples/BeforeSmooth.ply"); + new PLYMeshIO().save(TaubinSmoothing.smooth(mesh), "samples/SmoothedBit.ply"); + new PLYMeshIO().save(TaubinSmoothing.smooth(mesh, 10, 0.5, -0.53, TaubinWeightType.NAIVE), + "samples/SmoothedMore.ply"); + } + +}