Examples with TriangularFacet net.imagej.ops.geom.geom3d.mesh.TriangularFacet used on opensource projects

Example 6 with TriangularFacet

use of net.imagej.ops.geom.geom3d.mesh.TriangularFacet in project imagej-ops by imagej.

the class DefaultVoxelization3D method calculate.

@Override
public RandomAccessibleInterval<BitType> calculate(Mesh input) {
    Img<BitType> outImg = ops.create().img(new FinalInterval(width, height, depth), new BitType());
    DefaultMesh dMesh = (DefaultMesh) input;
    Set<RealLocalizable> verts = dMesh.getVertices();
    RealPoint minPoint = new RealPoint(verts.iterator().next());
    RealPoint maxPoint = new RealPoint(verts.iterator().next());
    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);
        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 (Facet f : dMesh.getFacets()) {
        TriangularFacet tri = (TriangularFacet) f;
        Vector3D v1 = tri.getP0();
        Vector3D v2 = tri.getP1();
        Vector3D v3 = tri.getP2();
        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) };
        // Should use the
        RandomAccess<BitType> ra = outImg.randomAccess();
        // 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 (// Don't check if voxel is already
                    !ra.get().get()) // 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;
}

Example 7 with TriangularFacet

use of net.imagej.ops.geom.geom3d.mesh.TriangularFacet in project imagej-ops by imagej.

the class DefaultConvexHull3D method assignPointsToFacets.

/**
 * Assigns all points which are not part of the convex hull to a facet. A
 * point is assigned to a facet if the point is in front of this facet. Every
 * point is assigned to only one facet. If a facet has a point in front the
 * facet is added to {@link DefaultConvexHull3D#facetsWithPointInFront}. After
 * this call {@link DefaultConvexHull3D#vertices} is empty. Points which are
 * behind all facets are removed because they are on the inside of the convex
 * hull.
 *
 * @param newFacets which could have a point in front
 * @param facetsWithPointInFront
 */
private void assignPointsToFacets(final double epsilon, final Set<Vertex> vertices, final List<TriangularFacet> newFacets, final List<TriangularFacet> facets, final List<TriangularFacet> facetsWithPointInFront) {
    Iterator<Vertex> vertexIt = vertices.iterator();
    while (vertexIt.hasNext()) {
        Vertex v = vertexIt.next();
        Iterator<TriangularFacet> facetIt = newFacets.iterator();
        TriangularFacet maxFacet = null;
        double maxdis = epsilon;
        while (facetIt.hasNext()) {
            TriangularFacet f = facetIt.next();
            double distanceToPlane = f.distanceToPlane(v);
            // point is assigned to the facet with maximum distance
            if (distanceToPlane > maxdis) {
                maxdis = distanceToPlane;
                maxFacet = f;
            }
        }
        // therefore on the inside of the convex hull.
        if (maxFacet != null) {
            maxFacet.setVertexInFront(v, maxdis);
            if (!facetsWithPointInFront.contains(maxFacet)) {
                facetsWithPointInFront.add(maxFacet);
            }
        }
    }
    facets.addAll(newFacets);
    // All vertices are reassigned or are inside of the convex hull.
    vertices.clear();
}

Example 8 with TriangularFacet

use of net.imagej.ops.geom.geom3d.mesh.TriangularFacet in project imagej-ops by imagej.

the class DefaultConvexHull3D method createSimplex.

/**
 * Computes an initial simplex of four facets. The simplex consists of the
 * four points v0-v3. v0 and v1 have the largest possible distance in one
 * dimension. v2 is the point with the largest distance to v0----v1. v3 is the
 * point with the largest distance to the plane described by v0, v1, v2.
 */
private double createSimplex(final Set<Vertex> vertices, final List<TriangularFacet> facets, final List<TriangularFacet> facetsWithPointInFront) {
    final Pair<Double, Vertex[]> minMax = computeMinMax(vertices);
    final double epsilon = minMax.getA();
    int i = getMaxDistPointIndex(minMax.getB());
    Vertex v0 = minMax.getB()[i];
    Vertex v1 = minMax.getB()[i + 3];
    vertices.remove(v0);
    vertices.remove(v1);
    Vertex v2 = getV2(epsilon, vertices, v0, v1);
    vertices.remove(v2);
    Vertex v3 = getV3(epsilon, vertices, v0, v1, v2);
    vertices.remove(v3);
    TriangularFacet f0 = new TriangularFacet(v0, v1, v2);
    if (f0.distanceToPlane(v3) > epsilon) {
        // change triangle orientation to counter clockwise
        Vertex tmp = v1;
        v1 = v2;
        v2 = tmp;
        f0 = new TriangularFacet(v0, v1, v2);
    }
    // v3 is behind f0
    assert f0.distanceToPlane(v3) < epsilon;
    TriangularFacet f1 = new TriangularFacet(v1, v0, v3);
    TriangularFacet f2 = new TriangularFacet(v2, v1, v3);
    TriangularFacet f3 = new TriangularFacet(v0, v2, v3);
    f0.setNeighbor(0, f3);
    f0.setNeighbor(1, f1);
    f0.setNeighbor(2, f2);
    f1.setNeighbor(0, f2);
    f1.setNeighbor(1, f0);
    f1.setNeighbor(2, f3);
    f2.setNeighbor(0, f3);
    f2.setNeighbor(1, f0);
    f2.setNeighbor(2, f1);
    f3.setNeighbor(0, f1);
    f3.setNeighbor(1, f0);
    f3.setNeighbor(2, f2);
    assert f0.distanceToPlane(v3) < epsilon;
    assert f1.distanceToPlane(v2) < epsilon;
    assert f2.distanceToPlane(v0) < epsilon;
    assert f3.distanceToPlane(v1) < epsilon;
    List<TriangularFacet> newFacets = new ArrayList<>();
    newFacets.add(f0);
    newFacets.add(f1);
    newFacets.add(f2);
    newFacets.add(f3);
    assignPointsToFacets(epsilon, vertices, newFacets, facets, facetsWithPointInFront);
    return epsilon;
}

Example 9 with TriangularFacet

use of net.imagej.ops.geom.geom3d.mesh.TriangularFacet in project imagej-ops by imagej.

the class DefaultMarchingCubes method calculate.

@SuppressWarnings({ "unchecked" })
@Override
public DefaultMesh calculate(final RandomAccessibleInterval<T> input) {
    DefaultMesh output = new DefaultMesh();
    ExtendedRandomAccessibleInterval<T, RandomAccessibleInterval<T>> extended = Views.extendValue(input, (T) new BoolType(false));
    Cursor<T> c = Views.interval(extended, new FinalInterval(new long[] { input.min(0) - 1, input.min(1) - 1, input.min(2) - 1 }, new long[] { input.max(0) + 1, input.max(1) + 1, input.max(2) + 1 })).localizingCursor();
    while (c.hasNext()) {
        c.next();
        int cursorX = c.getIntPosition(0);
        int cursorY = c.getIntPosition(1);
        int cursorZ = c.getIntPosition(2);
        Cursor<T> cu = getCube(extended, cursorX, cursorY, cursorZ);
        int i = 0;
        double[] vertex_values = new double[8];
        while (cu.hasNext()) {
            vertex_values[i++] = (cu.next().get()) ? 1 : 0;
        }
        // 6------7
        // /| /|
        // 2-----3 |
        // | 4---|-5
        // |/ |/
        // 0-----1
        vertex_values = mapFlatIterableToLookUpCube(vertex_values);
        // 4------5
        // /| /|
        // 7-----6 |
        // | 0---|-1
        // |/ |/
        // 3-----2
        int cubeindex = getCubeIndex(vertex_values);
        if (EDGE_TABLE[cubeindex] != 0) {
            int[] p0 = new int[] { 0 + cursorX, 0 + cursorY, 1 + cursorZ };
            int[] p1 = new int[] { 1 + cursorX, 0 + cursorY, 1 + cursorZ };
            int[] p2 = new int[] { 1 + cursorX, 0 + cursorY, 0 + cursorZ };
            int[] p3 = new int[] { 0 + cursorX, 0 + cursorY, 0 + cursorZ };
            int[] p4 = new int[] { 0 + cursorX, 1 + cursorY, 1 + cursorZ };
            int[] p5 = new int[] { 1 + cursorX, 1 + cursorY, 1 + cursorZ };
            int[] p6 = new int[] { 1 + cursorX, 1 + cursorY, 0 + cursorZ };
            int[] p7 = new int[] { 0 + cursorX, 1 + cursorY, 0 + cursorZ };
            double[][] vertlist = new double[12][];
            /* Find the vertices where the surface intersects the cube */
            if (0 != (EDGE_TABLE[cubeindex] & 1)) {
                vertlist[0] = interpolatePoint(p0, p1, vertex_values[0], vertex_values[1]);
            }
            if (0 != (EDGE_TABLE[cubeindex] & 2)) {
                vertlist[1] = interpolatePoint(p1, p2, vertex_values[1], vertex_values[2]);
            }
            if (0 != (EDGE_TABLE[cubeindex] & 4)) {
                vertlist[2] = interpolatePoint(p2, p3, vertex_values[2], vertex_values[3]);
            }
            if (0 != (EDGE_TABLE[cubeindex] & 8)) {
                vertlist[3] = interpolatePoint(p3, p0, vertex_values[3], vertex_values[0]);
            }
            if (0 != (EDGE_TABLE[cubeindex] & 16)) {
                vertlist[4] = interpolatePoint(p4, p5, vertex_values[4], vertex_values[5]);
            }
            if (0 != (EDGE_TABLE[cubeindex] & 32)) {
                vertlist[5] = interpolatePoint(p5, p6, vertex_values[5], vertex_values[6]);
            }
            if (0 != (EDGE_TABLE[cubeindex] & 64)) {
                vertlist[6] = interpolatePoint(p6, p7, vertex_values[6], vertex_values[7]);
            }
            if (0 != (EDGE_TABLE[cubeindex] & 128)) {
                vertlist[7] = interpolatePoint(p7, p4, vertex_values[7], vertex_values[4]);
            }
            if (0 != (EDGE_TABLE[cubeindex] & 256)) {
                vertlist[8] = interpolatePoint(p0, p4, vertex_values[0], vertex_values[4]);
            }
            if (0 != (EDGE_TABLE[cubeindex] & 512)) {
                vertlist[9] = interpolatePoint(p1, p5, vertex_values[1], vertex_values[5]);
            }
            if (0 != (EDGE_TABLE[cubeindex] & 1024)) {
                vertlist[10] = interpolatePoint(p2, p6, vertex_values[2], vertex_values[6]);
            }
            if (0 != (EDGE_TABLE[cubeindex] & 2048)) {
                vertlist[11] = interpolatePoint(p3, p7, vertex_values[3], vertex_values[7]);
            }
            /* Create the triangle */
            for (i = 0; TRIANGLE_TABLE[cubeindex][i] != -1; i += 3) {
                TriangularFacet face = new TriangularFacet(new Vertex(vertlist[TRIANGLE_TABLE[cubeindex][i + 2]][0], vertlist[TRIANGLE_TABLE[cubeindex][i + 2]][1], vertlist[TRIANGLE_TABLE[cubeindex][i + 2]][2]), new Vertex(vertlist[TRIANGLE_TABLE[cubeindex][i + 1]][0], vertlist[TRIANGLE_TABLE[cubeindex][i + 1]][1], vertlist[TRIANGLE_TABLE[cubeindex][i + 1]][2]), new Vertex(vertlist[TRIANGLE_TABLE[cubeindex][i]][0], vertlist[TRIANGLE_TABLE[cubeindex][i]][1], vertlist[TRIANGLE_TABLE[cubeindex][i]][2]));
                face.getArea();
                output.addFace(face);
            }
        }
    }
    return output;
}

Example 10 with TriangularFacet

use of net.imagej.ops.geom.geom3d.mesh.TriangularFacet in project imagej-ops by imagej.

the class DefaultVolumeMesh method calculate.

@Override
public DoubleType calculate(final Mesh input) {
    double volume = 0;
    for (Facet f : input.getFacets()) {
        TriangularFacet tf = (TriangularFacet) f;
        volume += signedVolumeOfTriangle(tf.getP0(), tf.getP1(), tf.getP2());
    }
    return new DoubleType(Math.abs(volume));
}