Polygon3D.java

/*
 * Copyright (C) 2012 Alberto Irurueta Carro (alberto@irurueta.com)
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *         http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package com.irurueta.geometry;

import java.io.Serializable;
import java.util.ArrayList;
import java.util.List;

/**
 * This class defines a polygon in 3D space.
 */
@SuppressWarnings("DuplicatedCode")
public class Polygon3D implements Serializable {

    /**
     * Default threshold value. Thresholds are used to determine whether a point
     * lies inside the polygon or not, or if it's locus or not, etc.
     */
    public static final double DEFAULT_THRESHOLD = 1e-9;

    /**
     * Minimum allowed threshold value.
     */
    public static final double MIN_THRESHOLD = 0.0;

    /**
     * Minimum number of vertices that a polygon is allowed to have.
     */
    public static final int MIN_VERTICES = 3;

    /**
     * Constant defining inhomogeneous coordinates.
     */
    public static final int INHOM_COORDS = 3;

    /**
     * Default method for triangulation.
     */
    public static final TriangulatorMethod DEFAULT_TRIANGULATOR_METHOD = TriangulatorMethod.VAN_GOGH_TRIANGULATOR;

    /**
     * List containing vertices of this polygon. Each vertex is a 3D point.
     */
    private List<Point3D> vertices;

    /**
     * Boolean indicating whether polygon has already been triangulated.
     */
    private boolean triangulated;

    /**
     * List containing triangles found after triangulating this polygon.
     * Initially this list will be null until triangulation is done.
     */
    private List<Triangle3D> triangles;

    /**
     * Method to do triangulation.
     */
    private TriangulatorMethod triangulatorMethod;

    /**
     * Constructor.
     *
     * @param vertices List of vertices forming this polygon.
     * @throws NotEnoughVerticesException Raised if list does not contain enough
     *                                    vertices.
     * @see #MIN_VERTICES
     */
    public Polygon3D(final List<Point3D> vertices) throws NotEnoughVerticesException {
        setVertices(vertices);
        triangulatorMethod = DEFAULT_TRIANGULATOR_METHOD;
    }

    /**
     * Returns triangulator method. Triangulator method determines the way a
     * polygon is divided into triangles.
     * If none has been provided DEFAULT_TRIANGULATOR_METHOD will be returned.
     *
     * @return Triangulator method.
     */
    public TriangulatorMethod getTriangulatorMethod() {
        return triangulatorMethod;
    }

    /**
     * Sets triangulator method. A triangulator method determines the way a
     * polygon is divided into triangles.
     *
     * @param triangulatorMethod A triangulator method.
     */
    public void setTriangulatorMethod(final TriangulatorMethod triangulatorMethod) {
        this.triangulatorMethod = triangulatorMethod;
    }

    /**
     * Returns the list of vertices forming this polygon.
     *
     * @return List of vertices.
     */
    public List<Point3D> getVertices() {
        return vertices;
    }

    /**
     * Sets list of vertices forming this polygon.
     *
     * @param vertices List of vertices.
     * @throws NotEnoughVerticesException Raised if provided list does not have
     *                                    enough vertices.
     * @see #MIN_VERTICES
     */
    public final void setVertices(final List<Point3D> vertices) throws NotEnoughVerticesException {
        if (vertices.size() < MIN_VERTICES) {
            throw new NotEnoughVerticesException();
        }

        if (vertices instanceof Serializable) {
            this.vertices = vertices;
        } else {
            this.vertices = new ArrayList<>(vertices);
        }
        triangulated = false;
        triangles = null;
    }

    /**
     * Determines whether this polygon has already been triangulated.
     * A polygon will only need to be triangulated once, unless the list of
     * vertices is reset.
     *
     * @return True if polygon has already been triangulated, false otherwise.
     */
    public boolean isTriangulated() {
        return triangulated;
    }

    /**
     * Returns a list of triangles forming this polygon.
     * This method checks whether this polygon has already been triangulated,
     * if not, it performs triangulation first.
     *
     * @return A list of triangles forming this polygon.
     * @throws TriangulatorException Raised if triangulation was needed and
     *                               failed.
     */
    public List<Triangle3D> getTriangles() throws TriangulatorException {
        if (!isTriangulated()) {
            triangulate();
        }
        return triangles;
    }

    /**
     * Returns signed area of this polygon.
     * The sign of the area determines whether vertices of the polygon are
     * provided in clockwise (negative sign) or clockwise (positive sign) order.
     *
     * @return Signed area of this polygon.
     */
    public double getArea() {
        final var origin = vertices.get(0);
        final var inhomX0 = origin.getInhomX();
        final var inhomY0 = origin.getInhomY();
        final var inhomZ0 = origin.getInhomZ();

        final var iterator = vertices.iterator();

        // because there are at least 3
        // vertices
        var prevPoint = iterator.next();
        Point3D curPoint;
        var avgX = 0.0;
        var avgY = 0.0;
        var avgZ = 0.0;
        while (iterator.hasNext()) {
            curPoint = iterator.next();

            final var inhomX1 = prevPoint.getInhomX();
            final var inhomY1 = prevPoint.getInhomY();
            final var inhomZ1 = prevPoint.getInhomZ();

            final var inhomX2 = curPoint.getInhomX();
            final var inhomY2 = curPoint.getInhomY();
            final var inhomZ2 = curPoint.getInhomZ();

            // compute cross product of ab = (prevPoint - origin) and
            // ac = (curPoint - origin)
            final var abX = inhomX1 - inhomX0;
            final var abY = inhomY1 - inhomY0;
            final var abZ = inhomZ1 - inhomZ0;

            final var acX = inhomX2 - inhomX0;
            final var acY = inhomY2 - inhomY0;
            final var acZ = inhomZ2 - inhomZ0;

            final var crossX = abY * acZ - abZ * acY;
            final var crossY = abZ * acX - abX * acZ;
            final var crossZ = abX * acY - abY * acX;

            avgX += crossX;
            avgY += crossY;
            avgZ += crossZ;

            prevPoint = curPoint;
        }

        return 0.5 * Math.sqrt(avgX * avgX + avgY * avgY + avgZ * avgZ);
    }

    /**
     * Returns perimeter of this polygon.
     * The perimeter is computed as the sum of the distances between consecutive
     * pairs of vertices.
     *
     * @return Perimeter of this polygon.
     */
    public double getPerimeter() {
        // iterate over all vertices and compute their distance
        final var iterator = vertices.iterator();
        var prevPoint = iterator.next();
        Point3D point;
        var perimeter = 0.0;
        while (iterator.hasNext()) {
            point = iterator.next();
            perimeter += prevPoint.distanceTo(point);
            prevPoint = point;
        }
        // get distance from last point with first one
        perimeter += prevPoint.distanceTo(vertices.get(0));
        return perimeter;
    }

    /**
     * Determines if provided point lies within the region defined by this
     * polygon.
     * Notice that this method is only ensured to work for polygons having no
     * holes or crossing borders. It will safely work on any other polygon,
     * no matter if it is regular, non-regular, convex or concave.
     *
     * @param point Point to be checked.
     * @return True if point lies within the area defined by this polygon, false
     * otherwise.
     * @throws TriangulatorException Raised if triangulation was required but
     *                               failed.
     */
    public boolean isInside(final Point3D point) throws TriangulatorException {
        return isInside(point, DEFAULT_THRESHOLD);
    }

    /**
     * Determines if provided point lies within the region defined by this
     * polygon.
     * Notice that this method is only ensured to work for polygons having no
     * holes or crossing borders. It will safely work on any other polygon,
     * no matter if it is regular, non-regular, convex or concave.
     *
     * @param point     Point to be checked.
     * @param threshold Threshold to determine whether point lies inside this
     *                  polygon. Usually this value should be small.
     * @return True if point lies within the area defined by this polygon, false
     * otherwise.
     * @throws IllegalArgumentException Raised if provided threshold is negative.
     * @throws TriangulatorException    Raised if triangulation was required but
     *                                  failed.
     */
    public boolean isInside(final Point3D point, final double threshold)
            throws TriangulatorException {
        if (threshold < MIN_THRESHOLD) {
            throw new IllegalArgumentException();
        }

        for (final var triangle : getTriangles()) {
            if (triangle.isInside(point, threshold)) {
                return true;
            }
        }
        return false;
    }

    /**
     * Returns the center of this polygon.
     * The center is the average point among all the vertices of this polygon.
     * The center is not ensure to lie within the area formed by this polygon.
     *
     * @return Center of this polygon.
     */
    public Point3D getCenter() {
        final var result = Point3D.create();
        center(result);
        return result;
    }

    /**
     * Computes the center of this polygon.
     * The center is the average point among all the vertices of this polygon.
     * The center is not ensured to lie within the area formed by this polygon.
     *
     * @param result Instance where the computed center will be stored.
     */
    public void center(final Point3D result) {
        // compute average location of all vertices
        var inhomX = 0.0;
        var inhomY = 0.0;
        var inhomZ = 0.0;
        final var total = vertices.size();

        for (final var point : vertices) {
            inhomX += point.getInhomX() / total;
            inhomY += point.getInhomY() / total;
            inhomZ += point.getInhomZ() / total;
        }
        result.setInhomogeneousCoordinates(inhomX, inhomY, inhomZ);
    }

    /**
     * Determines whether provided point is locus of the borders defined by
     * the vertices of this polygon. A point will be locus if it lies in the
     * line defined by two consecutive vertices up to a certain threshold of
     * error.
     *
     * @param point     Point to be checked.
     * @param threshold Threshold of allowed error. This should usually be a
     *                  small value.
     * @return True if provided point lies in a border of this polygon, false
     * otherwise.
     * @throws IllegalArgumentException Raised if provided threshold is negative.
     */
    public boolean isLocus(final Point3D point, final double threshold) {
        if (threshold < MIN_THRESHOLD) {
            throw new IllegalArgumentException();
        }

        // normalize point to increase accuracy
        point.normalize();

        final var iterator = vertices.iterator();
        // it's ok because there are at
        // least 3 vertices
        var prevPoint = iterator.next();
        // normalize to increase accuracy
        prevPoint.normalize();

        Point3D curPoint;
        while (iterator.hasNext()) {
            curPoint = iterator.next();
            // normalize to increase accuracy
            curPoint.normalize();
            if (point.isBetween(prevPoint, curPoint, threshold)) {
                return true;
            }
            prevPoint = curPoint;
        }

        // check last point with first
        return point.isBetween(prevPoint, vertices.get(0), threshold);
    }

    /**
     * Determines whether provided point is locus of the borders defined by the
     * vertices of this polygon. A point will be locus if it lies in the line
     * defined by two consecutive vertices.
     *
     * @param point Point to be checked.
     * @return True if provided point lies in a border of this polygon, false
     * otherwise.
     */
    public boolean isLocus(final Point3D point) {
        return isLocus(point, DEFAULT_THRESHOLD);
    }

    /**
     * Returns the shortest distance from provided point to a border of this
     * polygon. Note that borders are segments defined by consecutive vertices
     *
     * @param point Point to be checked.
     * @return Shortest distance from provided point to this polygon.
     * @throws CoincidentPointsException Raised if points in a polygon are too
     *                                   close. This usually indicates numerical instability or polygon degeneracy.
     */
    public double getShortestDistance(final Point3D point) throws CoincidentPointsException {
        // iterate over all vertices and compute their distance
        var iterator = vertices.iterator();
        var prevPoint = iterator.next();
        // to increase accuracy
        prevPoint.normalize();
        Point3D curPoint;
        var bestDist = Double.MAX_VALUE;
        double dist;
        var found = false;
        Line3D line = null;
        final var pointInLine = Point3D.create();

        while (iterator.hasNext()) {
            curPoint = iterator.next();
            // to increase accuracy
            curPoint.normalize();

            // check if point lies in the segment of the boundary of this polygon
            if (point.isBetween(curPoint, prevPoint)) {
                return 0.0;
            }

            if (line == null) {
                line = new Line3D(curPoint, prevPoint);
            } else {
                line.setPlanesFromPoints(curPoint, prevPoint);
            }
            // to increase accuracy
            line.normalize();

            // find the closest point to line
            line.closestPoint(point, pointInLine);
            // to increase accuracy
            pointInLine.normalize();

            if (pointInLine.isBetween(curPoint, prevPoint)) {
                // closest point lies within segment of polygon boundary, so we
                // keep distance
                dist = point.distanceTo(pointInLine);
                if (dist < bestDist) {
                    // a better point has been found
                    bestDist = dist;
                    found = true;
                }
            }

            prevPoint = curPoint;
        }

        // try last vertex with first
        // check if point lies in the segment of the boundary of this polygon
        final var first = vertices.get(0);
        if (point.isBetween(prevPoint, first)) {
            return 0.0;
        }

        if (line == null) {
            line = new Line3D(prevPoint, first);
        } else {
            line.setPlanesFromPoints(prevPoint, first);
        }
        // to increase accuracy
        line.normalize();

        // find the closest point to line
        line.closestPoint(point, pointInLine);
        // to increase accuracy
        pointInLine.normalize();

        if (pointInLine.isBetween(prevPoint, first)) {
            // closest point lies within segment of polygon boundary, so we
            // keep distance
            dist = point.distanceTo(pointInLine);
            if (dist < bestDist) {
                // a better point has been found
                bestDist = dist;
                found = true;
            }
        }

        if (!found) {
            // no closest point was found on a segment belonging to polygon
            // boundary, so we search for the closest vertex
            iterator = vertices.iterator();
            while (iterator.hasNext()) {
                // a better vertex has been found
                curPoint = iterator.next();
                dist = point.distanceTo(curPoint);
                if (dist < bestDist) {
                    bestDist = dist;
                }
            }
        }

        return bestDist;
    }

    /**
     * Returns the closest point to provided point that is locus of this
     * polygon (i.e. lies on a border of this polygon).
     *
     * @param point Point to be checked.
     * @return Closest point being locus of this polygon.
     * @throws CoincidentPointsException Raised if points in a polygon are too
     *                                   close. This usually indicates numerical instability or polygon degeneracy.
     */
    public Point3D getClosestPoint(final Point3D point) throws CoincidentPointsException {
        final var result = Point3D.create();
        closestPoint(point, result);
        return result;
    }

    /**
     * Computes the closes point to provided point that is locus of this
     * polygon (i.e. lies on a border of this polygon).
     *
     * @param point  Point to be checked.
     * @param result Instance where the closest point will be stored.
     * @throws CoincidentPointsException Raised if points in a polygon are too
     *                                   close. This usually indicates numerical instability or polygon degeneracy.
     */
    public void closestPoint(final Point3D point, final Point3D result) throws CoincidentPointsException {
        // iterate over all vertices and compute their distance
        var iterator = vertices.iterator();
        var prevPoint = iterator.next();
        // to increase accuracy
        prevPoint.normalize();

        Point3D curPoint;
        var bestDist = Double.MAX_VALUE;
        double dist;
        var found = false;
        Line3D line = null;
        final var pointInLine = Point3D.create();

        while (iterator.hasNext()) {
            curPoint = iterator.next();
            // to increase accuracy
            curPoint.normalize();

            // check if point lies in the segment of the boundary of this polygon
            if (point.isBetween(curPoint, prevPoint)) {
                result.setCoordinates(point);
                return;
            }

            if (line == null) {
                line = new Line3D(curPoint, prevPoint);
            } else {
                line.setPlanesFromPoints(curPoint, prevPoint);
            }
            // to increase accuracy
            line.normalize();

            // find the closest point to line
            line.closestPoint(point, pointInLine);
            // to increase accuracy
            pointInLine.normalize();

            if (pointInLine.isBetween(curPoint, prevPoint)) {
                // closest point lies within segment of polygon boundary, so we
                // keep distance and point
                dist = point.distanceTo(pointInLine);
                if (dist < bestDist) {
                    // a better point has been found
                    bestDist = dist;
                    result.setCoordinates(pointInLine);
                    found = true;
                }
            }

            prevPoint = curPoint;
        }

        // try last vertex with first
        // check if point lies in the segment of the boundary of this polygon
        final var first = vertices.get(0);
        if (point.isBetween(prevPoint, first)) {
            result.setCoordinates(point);
            return;
        }

        if (line == null) {
            line = new Line3D(prevPoint, first);
        } else {
            line.setPlanesFromPoints(prevPoint, first);
        }
        // to increase accuracy
        line.normalize();

        // find the closest point to line
        line.closestPoint(point, pointInLine);
        // to increase accuracy
        pointInLine.normalize();

        if (pointInLine.isBetween(prevPoint, first)) {
            // closest point lies within segment of polygon boundary, so we
            // keep distance
            dist = point.distanceTo(pointInLine);
            if (dist < bestDist) {
                // a better point has been found
                bestDist = dist;
                result.setCoordinates(pointInLine);
                found = true;
            }
        }


        if (!found) {
            // no closest point was found on a segment belonging to polygon
            // boundary, so we search for the closest vertex
            iterator = vertices.iterator();
            while (iterator.hasNext()) {
                curPoint = iterator.next();
                dist = point.distanceTo(curPoint);
                if (dist < bestDist) {
                    // a better vertex has been found
                    bestDist = dist;
                    result.setCoordinates(curPoint);
                }
            }
        }
    }

    /**
     * Triangulates this polygon using this polygon's triangulator method.
     * A polygon only will be triangulated once when required or this method is
     * called.
     * This method will make no action if a polygon is already triangulated
     * unless it's vertices are reset.
     *
     * @throws TriangulatorException Raised if triangulation failed.
     * @see #getTriangulatorMethod
     * @see #setTriangulatorMethod(TriangulatorMethod)
     */
    public void triangulate() throws TriangulatorException {
        if (!triangulated) {
            final var triangulator = Triangulator3D.create(triangulatorMethod);
            triangles = triangulator.triangulate(vertices);
            triangulated = true;
        }
    }

    /**
     * Computes the average orientation of a 3D polygon as an array
     * containing the vector coordinates of the orientation.
     * Notice that if all the vertices of the polygon lie in the same plane,
     * then the estimated orientation will be exact, otherwise an average
     * orientation will be estimated.
     * Estimated orientation will be normalized (the norm of the estimated
     * vector will be 1).
     *
     * @param vertices  List of vertices forming a polygon in 3D.
     * @param result    Array where the estimated orientation will be stored.
     * @param threshold Threshold to determine whether the orientation can be
     *                  estimated or not. Because the estimated orientation needs to be
     *                  normalized, it will only be possible to do so when norm has a reasonable
     *                  value. If estimated norm happens to be smaller than provided value, then
     *                  it will be assumed that polygon contains a degeneracy or that a
     *                  consecutive pair of vertices are coincident.
     * @throws IllegalArgumentException  Raised if provided result array does not
     *                                   have length 3, or if the list of provided vertices does not contain at
     *                                   least three vertices, or if provided threshold is negative.
     * @throws CoincidentPointsException Raised usually when it is determined
     *                                   that consecutive vertices of the polygon are too close (i.e. coincident)
     *                                   or when there are polygon degeneracies or numerical instabilities.
     */
    public static void orientation(
            final List<Point3D> vertices, final double[] result, final double threshold)
            throws CoincidentPointsException {
        final var numVertices = vertices.size();
        if (numVertices < MIN_VERTICES) {
            throw new IllegalArgumentException();
        }
        if (result.length != INHOM_COORDS) {
            throw new IllegalArgumentException();
        }
        if (threshold < MIN_THRESHOLD) {
            throw new IllegalArgumentException();
        }

        var avgX = 0.0;
        var avgY = 0.0;
        var avgZ = 0.0;

        final var iterator = vertices.iterator();
        var prevPoint = iterator.next();
        var origin = prevPoint;
        Point3D curPoint;
        final var inhomX0 = origin.getInhomX();
        final var inhomY0 = origin.getInhomY();
        final var inhomZ0 = origin.getInhomZ();

        while (iterator.hasNext()) {
            curPoint = iterator.next();
            final var inhomX1 = prevPoint.getInhomX();
            final var inhomY1 = prevPoint.getInhomY();
            final var inhomZ1 = prevPoint.getInhomZ();

            final var inhomX2 = curPoint.getInhomX();
            final var inhomY2 = curPoint.getInhomY();
            final var inhomZ2 = curPoint.getInhomZ();

            // compute cross product of ab = (vertex1 - origin) and ac = (vertex2
            // - origin)
            final var abX = inhomX1 - inhomX0;
            final var abY = inhomY1 - inhomY0;
            final var abZ = inhomZ1 - inhomZ0;

            final var acX = inhomX2 - inhomX0;
            final var acY = inhomY2 - inhomY0;
            final var acZ = inhomZ2 - inhomZ0;

            final var crossX = abY * acZ - abZ * acY;
            final var crossY = abZ * acX - abX * acZ;
            final var crossZ = abX * acY - abY * acX;

            avgX += crossX;
            avgY += crossY;
            avgZ += crossZ;

            prevPoint = curPoint;
        }

        final var norm = Math.sqrt(avgX * avgX + avgY * avgY + avgZ * avgZ);

        if (norm < threshold) {
            throw new CoincidentPointsException();
        }

        avgX /= norm;
        avgY /= norm;
        avgZ /= norm;

        result[0] = avgX;
        result[1] = avgY;
        result[2] = avgZ;
    }

    /**
     * Computes the average orientation of a 3D polygon as an array
     * containing the vector coordinates of the orientation.
     * Notice that if all the vertices of the polygon lie in the same plane,
     * then the estimated orientation will be exact, otherwise an average
     * orientation will be estimated.
     * Estimated orientation will be normalized (the norm of the estimated
     * vector will be 1).
     *
     * @param vertices List of vertices forming a polygon in 3D.
     * @param result   Array where the estimated orientation will be stored.
     * @throws IllegalArgumentException  Raised if provided result array does not
     *                                   have length 3, or if the list of provided vertices does not contain at
     *                                   least three vertices.
     * @throws CoincidentPointsException Raised usually when it is determined
     *                                   that consecutive vertices of the polygon are too close (i.e. coincident)
     *                                   or when there are polygon degeneracies or numerical instabilities.
     */
    public static void orientation(
            final List<Point3D> vertices, final double[] result) throws CoincidentPointsException {
        orientation(vertices, result, DEFAULT_THRESHOLD);
    }

    /**
     * Returns the average orientation of a 3D polygon as an array containing
     * the vector coordinates of the orientation.
     * Notice that if all the vertices of the polygon lie in the same plane,
     * then the estimated orientation will be exact, otherwise an average
     * orientation will be estimated.
     * Estimated orientation will be normalized (the norm of the estimated
     * vector will be 1).
     *
     * @param vertices  List of vertices forming a polygon in 3D.
     * @param threshold Threshold to determine whether the orientation can be
     *                  estimated or not. Because the estimated orientation needs to be
     *                  normalized, it will only be possible to do so when norm has a reasonable
     *                  value. If estimated norm happens to be smaller than provided value, then
     *                  it will be assumed that polygon contains a degeneracy or that a
     *                  consecutive pair of vertices are coincident.
     * @return Array containing estimated orientation.
     * @throws IllegalArgumentException  Raised if the list of provided vertices
     *                                   does not contain at least three vertices, or if provided threshold is
     *                                   negative.
     * @throws CoincidentPointsException Raised usually when it is determined
     *                                   that consecutive vertices of the polygon are too close (i.e. coincident)
     *                                   or when there are polygon degeneracies or numerical instabilities.
     */
    public static double[] orientation(
            final List<Point3D> vertices, final double threshold) throws CoincidentPointsException {
        final var out = new double[INHOM_COORDS];
        orientation(vertices, out, threshold);
        return out;
    }

    /**
     * Returns the average orientation of a 3D polygon as an array containing
     * the vector coordinates of the orientation.
     * Notice that if all the vertices of the polygon lie in the same plane,
     * then the estimated orientation will be exact, otherwise an average
     * orientation will be estimated.
     * Estimated orientation will be normalized (the norm of the estimated
     * vector will be 1).
     *
     * @param vertices List of vertices forming a polygon in 3D
     * @return Array containing estimated orientation.
     * @throws IllegalArgumentException  Raised if the list of provided vertices
     *                                   does not contain at least three vertices.
     * @throws CoincidentPointsException Raised usually when it is determined
     *                                   that consecutive vertices of the polygon are too close (i.e. coincident)
     *                                   or when there are polygon degeneracies or numerical instabilities.
     */
    public static double[] orientation(final List<Point3D> vertices) throws CoincidentPointsException {
        return orientation(vertices, DEFAULT_THRESHOLD);
    }

    /**
     * Computes the average orientation of provided 3D polygon as an array
     * containing the vector coordinates of the orientation.
     * Notice that if all the vertices of the polygon lie in the same plane,
     * then the estimated orientation will be exact, otherwise an average
     * orientation will be estimated.
     * Estimated orientation will be normalized (the norm of the estimated
     * vector will be 1).
     *
     * @param polygon   A polygon in 3D.
     * @param result    Array where the estimated orientation will be stored.
     * @param threshold Threshold to determine whether the orientation can be
     *                  estimated or not. Because the estimated orientation needs to be
     *                  normalized, it will only be possible to do so when norm has a reasonable
     *                  value. If estimated norm happens to be smaller than provided value, then
     *                  it will be assumed that polygon contains a degeneracy or that a
     *                  consecutive pair of vertices are coincident.
     * @throws IllegalArgumentException  Raised if provided result array does not
     *                                   have length 3 or if provided threshold is negative.
     * @throws CoincidentPointsException Raised usually when it is determined
     *                                   that consecutive vertices of the polygon are too close (i.e. coincident)
     *                                   or when there are polygon degeneracies or numerical instabilities.
     */
    public static void orientation(
            final Polygon3D polygon, final double[] result, final double threshold) throws CoincidentPointsException {
        orientation(polygon.getVertices(), result, threshold);
    }

    /**
     * Computes the average orientation of provided 3D polygon as an array
     * containing the vector coordinates of the orientation.
     * Notice that if all the vertices of the polygon lie in the same plane,
     * then the estimated orientation will be exact, otherwise an average
     * orientation will be estimated.
     * Estimated orientation will be normalized (the norm of the estimated
     * vector will be 1).
     *
     * @param polygon A polygon in 3D.
     * @param result  Array where the estimated orientation will be stored.
     * @throws IllegalArgumentException  Raised if provided result array does not
     *                                   have length 3.
     * @throws CoincidentPointsException Raised usually when it is determined
     *                                   that consecutive vertices of the polygon are too close (i.e. coincident)
     *                                   or when there are polygon degeneracies or numerical instabilities.
     */
    public static void orientation(
            final Polygon3D polygon, final double[] result) throws CoincidentPointsException {
        orientation(polygon, result, DEFAULT_THRESHOLD);
    }

    /**
     * Returns the average orientation of provided 3D polygon as an array
     * containing the vector coordinates of the orientation.
     * Notice that if all the vertices of the polygon lie in the same plane,
     * then the estimated orientation will be exact, otherwise an average
     * orientation will be estimated.
     * Estimated orientation will be normalized (the norm of the estimated
     * vector will be 1).
     *
     * @param polygon   A polygon in 3D.
     * @param threshold Threshold to determine whether the orientation can be
     *                  estimated or not. Because the estimated orientation needs to be
     *                  normalized, it will only be possible to do so when norm has a reasonable
     *                  value. If estimated norm happens to be smaller than provided value, then
     *                  it will be assumed that polygon contains a degeneracy or that a
     *                  consecutive pair of vertices are coincident.
     * @return Array containing estimated orientation.
     * @throws IllegalArgumentException  Raised if provided threshold is
     *                                   negative.
     * @throws CoincidentPointsException Raised usually when it is determined
     *                                   that consecutive vertices of the polygon are too close (i.e. coincident)
     *                                   or when there are polygon degeneracies or numerical instabilities.
     */
    public static double[] orientation(
            final Polygon3D polygon, final double threshold) throws CoincidentPointsException {
        return orientation(polygon.getVertices(), threshold);
    }

    /**
     * Returns the average orientation of provided 3D polygon as an array
     * containing the vector coordinates of the orientation.
     * Notice that if all the vertices of the polygon lie in the same plane,
     * then the estimated orientation will be exact, otherwise an average
     * orientation will be estimated.
     * Estimated orientation will be normalized (the norm of the estimated
     * vector will be 1).
     *
     * @param polygon A polygon in 3D.
     * @return Array containing estimated orientation.
     * @throws CoincidentPointsException Raised usually when it is determined
     *                                   that consecutive vertices of the polygon are too close (i.e. coincident)
     *                                   or when there are polygon degeneracies or numerical instabilities.
     */
    public static double[] orientation(final Polygon3D polygon) throws CoincidentPointsException {
        return orientation(polygon, DEFAULT_THRESHOLD);
    }

    /**
     * Computes the average orientation of this polygon as an array containing
     * the vector coordinates of the orientation.
     * Notice that if all the vertices of the polygon lie in the same plane,
     * then the estimated orientation will be exact, otherwise an average
     * orientation will be estimated.
     * Estimated orientation will be normalized (the norm of the estimated
     * vector will be 1).
     *
     * @param result    Array where the estimated orientation will be stored.
     * @param threshold Threshold to determine whether the orientation can be
     *                  estimated or not. Because the estimated orientation needs to be
     *                  normalized, it will only be possible to do so when norm has a reasonable
     *                  value. If estimated norm happens to be smaller than provided value, then
     *                  it will be assumed that polygon contains a degeneracy or that a
     *                  consecutive pair of vertices are coincident.
     * @throws IllegalArgumentException  Raised if provided result array does not
     *                                   have length 3 or if provided threshold is negative.
     * @throws CoincidentPointsException Raised usually when it is determined
     *                                   that consecutive vertices of the polygon are too close (i.e. coincident)
     *                                   or when there are polygon degeneracies or numerical instabilities.
     */
    public void orientation(final double[] result, final double threshold) throws CoincidentPointsException {
        orientation(vertices, result, threshold);
    }

    /**
     * Computes the average orientation of provided 3D polygon as an array
     * containing the vector coordinates of the orientation.
     * Notice that if all the vertices of the polygon lie in the same plane,
     * then the estimated orientation will be exact, otherwise an average
     * orientation will be estimated.
     * Estimated orientation will be normalized (the norm of the estimated
     * vector will be 1).
     *
     * @param result Array where the estimated orientation will be stored.
     * @throws IllegalArgumentException  Raised if provided result array does not
     *                                   have length 3.
     * @throws CoincidentPointsException Raised usually when it is determined
     *                                   that consecutive vertices of the polygon are too close (i.e. coincident)
     *                                   or when there are polygon degeneracies or numerical instabilities.
     */
    public void orientation(final double[] result) throws CoincidentPointsException {
        orientation(result, DEFAULT_THRESHOLD);
    }

    /**
     * Returns the average orientation of provided 3D polygon as an array
     * containing the vector coordinates of the orientation.
     * Notice that if all the vertices of the polygon lie in the same plane,
     * then the estimated orientation will be exact, otherwise an average
     * orientation will be estimated.
     * Estimated orientation will be normalized (the norm of the estimated
     * vector will be 1).
     *
     * @param threshold Threshold to determine whether the orientation can be
     *                  estimated or not. Because the estimated orientation needs to be
     *                  normalized, it will only be possible to do so when norm has a reasonable
     *                  value. If estimated norm happens to be smaller than provided value, then
     *                  it will be assumed that polygon contains a degeneracy or that a
     *                  consecutive pair of vertices are coincident.
     * @return Array containing estimated orientation.
     * @throws IllegalArgumentException  Raised if provided threshold is
     *                                   negative.
     * @throws CoincidentPointsException Raised usually when it is determined
     *                                   that consecutive vertices of the polygon are too close (i.e. coincident)
     *                                   or when there are polygon degeneracies or numerical instabilities.
     */
    public double[] getOrientation(final double threshold) throws CoincidentPointsException {
        return orientation(vertices, threshold);
    }

    /**
     * Returns the average orientation of provided 3D polygon as an array
     * containing the vector coordinates of the orientation.
     * Notice that if all the vertices of the polygon lie in the same plane,
     * then the estimated orientation will be exact, otherwise an average
     * orientation will be estimated.
     * Estimated orientation will be normalized (the norm of the estimated
     * vector will be 1).
     *
     * @return Array containing estimated orientation.
     * @throws CoincidentPointsException Raised usually when it is determined
     *                                   that consecutive vertices of the polygon are too close (i.e. coincident)
     *                                   or when there are polygon degeneracies or numerical instabilities.
     */
    public double[] getOrientation() throws CoincidentPointsException {
        return getOrientation(DEFAULT_THRESHOLD);
    }

    /**
     * Returns the angle between two polygons, assuming that all vertices of
     * each polygon lie on a given plane. Hence, this is equivalent to
     * estimating the angle between the planes formed by two polygons.
     * The angle between two polygons is estimated by first estimating their
     * orientation.
     *
     * @param polygon1  1st polygon.
     * @param polygon2  2nd polygon.
     * @param threshold Threshold to determine when polygon orientation can be
     *                  estimated.
     * @return Angle between two polygons in radians.
     * @throws IllegalArgumentException  Raised if provided threshold is negative.
     * @throws CoincidentPointsException Raised usually when it is determined
     *                                   that consecutive vertices of the polygon are too close (i.e. coincident)
     *                                   or when there are polygon degeneracies or numerical instabilities.
     * @see #orientation(Polygon3D, double[])
     */
    public static double getAngleBetweenPolygons(
            final Polygon3D polygon1, final Polygon3D polygon2, final double threshold)
            throws CoincidentPointsException {
        return getAngleBetweenPolygons(polygon1.getVertices(), polygon2.getVertices(), threshold);
    }

    /**
     * Returns the angle between two polygons, assuming that all vertices of
     * each polygon lie on a given plane. Hence, this is equivalent to
     * estimating the angle between the planes formed by two polygons.
     * The angle between two polygons is estimated by first estimating their
     * orientation.
     *
     * @param polygon1 1st polygon.
     * @param polygon2 2nd polygon.
     * @return Angle between two polygons in radians.
     * @throws CoincidentPointsException Raised usually when it is determined
     *                                   that consecutive vertices of the polygon are too close (i.e. coincident)
     *                                   or when there are polygon degeneracies or numerical instabilities.
     * @see #orientation(Polygon3D, double[])
     */
    public static double getAngleBetweenPolygons(final Polygon3D polygon1, final Polygon3D polygon2)
            throws CoincidentPointsException {
        return getAngleBetweenPolygons(polygon1, polygon2, DEFAULT_THRESHOLD);
    }

    /**
     * Returns the angle between two polygons formed each of them by the
     * corresponding list of provided vertices and assuming that all vertices of
     * each polygon lie on a given plane. Hence, this is equivalent to
     * estimating the angle between the planes formed by two polygons.
     * The angle between two polygons is estimated by first estimating their
     * orientation.
     *
     * @param vertices1 Vertices of 1st polygon.
     * @param vertices2 2nd polygon.
     * @param threshold Threshold to determine when polygon orientation can be
     *                  estimated.
     * @return Angle between two polygons in radians.
     * @throws IllegalArgumentException  Raised if provided threshold is negative
     *                                   or if list of vertices do not contain at least 3 vertices for each
     *                                   polygon.
     * @throws CoincidentPointsException Raised usually when it is determined
     *                                   that consecutive vertices of the polygon are too close (i.e. coincident)
     *                                   or when there are polygon degeneracies or numerical instabilities.
     * @see #orientation(Polygon3D, double[])
     */
    public static double getAngleBetweenPolygons(
            final List<Point3D> vertices1, final List<Point3D> vertices2, final double threshold)
            throws CoincidentPointsException {
        return getAngleBetweenOrientations(Polygon3D.orientation(vertices1, threshold),
                Polygon3D.orientation(vertices2, threshold));
    }

    /**
     * Returns the angle between two polygons formed each of them by the
     * corresponding list of provided vertices and assuming that all vertices of
     * each polygon lie on a given plane. Hence, this is equivalent to
     * estimating the angle between the planes formed by two polygons.
     * The angle between two polygons is estimated by first estimating their
     * orientation.
     *
     * @param vertices1 Vertices of 1st polygon.
     * @param vertices2 2nd polygon.
     * @return Angle between two polygons in radians.
     * @throws IllegalArgumentException  Raised if list of vertices do not
     *                                   contain at least 3 vertices for each polygon.
     * @throws CoincidentPointsException Raised usually when it is determined
     *                                   that consecutive vertices of the polygon are too close (i.e. coincident)
     *                                   or when there are polygon degeneracies or numerical instabilities.
     * @see #orientation(Polygon3D, double[])
     */
    public static double getAngleBetweenPolygons(
            final List<Point3D> vertices1, final List<Point3D> vertices2) throws CoincidentPointsException {
        return getAngleBetweenPolygons(vertices1, vertices2, DEFAULT_THRESHOLD);
    }

    /**
     * Internal method to compute polygon orientation from their respective
     * orientation vectors.
     * Orientation vectors are provided as arrays and can be obtained by calling
     * orientation(Polygon3D, double[]) among other methods.
     *
     * @param orientation1 Orientation of 1st polygon.
     * @param orientation2 Orientation of 2nd polygon.
     * @return Angle between two polygons in radians.
     * @throws IllegalArgumentException Raised if any of the orientation arrays
     *                                  do not have length 3.
     */
    private static double getAngleBetweenOrientations(final double[] orientation1, final double[] orientation2) {
        if (orientation1.length != INHOM_COORDS || orientation2.length != INHOM_COORDS) {
            throw new IllegalArgumentException();
        }

        final var x1 = orientation1[0];
        final var y1 = orientation1[1];
        final var z1 = orientation1[2];

        final var x2 = orientation2[0];
        final var y2 = orientation2[1];
        final var z2 = orientation2[2];

        final var norm1 = Math.sqrt(x1 * x1 + y1 * y1 + z1 * z1);
        final var norm2 = Math.sqrt(x2 * x2 + y2 * y2 + z2 * z2);

        final var dotProduct = (x1 * x2 + y1 * y2 + z1 * z2) / (norm1 * norm2);

        return Math.acos(dotProduct);
    }
}