| gvSIG desktop 1

package org.gvsig.jts.voronoi.chew;

/*

 * Copyright (c) 2005 by L. Paul Chew.

 * 

 * Permission is hereby granted, without written agreement and without

 * license or royalty fees, to use, copy, modify, and distribute this

 * software and its documentation for any purpose, subject to the following 

 * conditions:

 *

 * The above copyright notice and this permission notice shall be included 

 * in all copies or substantial portions of the Software.

 * 

 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 

 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 

 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 

 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 

 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING 

 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER 

 * DEALINGS IN THE SOFTWARE.

 */

/**

 * Points in Euclidean space, implemented as double[].

 * 

 * Includes simple geometric operations.

 * Uses matrices; a matrix is represented as an array of Pnts.

 * Uses simplices; a simplex is represented as an array of Pnts.

 * 

 * @author Paul Chew

 * 

 * Created July 2005.  Derived from an earlier, messier version.

 */

public class Pnt {

    private double[] coordinates;          // The point's coordinates

    /**

     * Constructor.

     * @param coords the coordinates

     */

    public Pnt (double[] coords) {

        // Copying is done here to ensure that Pnt's coords cannot be altered.

        coordinates = new double[coords.length];

        System.arraycopy(coords, 0, coordinates, 0, coords.length);

    }

    /**

     * Constructor.

     * @param coordA

     * @param coordB

     */

    public Pnt (double coordA, double coordB) {

        this(new double[] {coordA, coordB});

    }

    /**

     * Constructor.

     * @param coordA

     * @param coordB

     * @param coordC

     */

    public Pnt (double coordA, double coordB, double coordC) {

        this(new double[] {coordA, coordB, coordC});

    }

    /**

     * Create a String for this Pnt.

     * @return a String representation of this Pnt.

     */

    public String toString () {

        if (coordinates.length == 0) return "()";

        String result = "Pnt(" + coordinates[0];

        for (int i = 1; i < coordinates.length; i++)

            result = result + "," + coordinates[i];

        result = result + ")";

        return result;

    }

    /**

     * Equality.

     * @param other the other Object to compare to

     * @return true iff the Pnts have the same coordinates

     */

    public boolean equals (Object other) {

        if (!(other instanceof Pnt)) return false;

        Pnt p = (Pnt) other;

        if (this.coordinates.length != p.coordinates.length) return false;

        for (int i = 0; i < this.coordinates.length; i++)

            if (this.coordinates[i] != p.coordinates[i]) return false;

        return true;

    }

    /**

     * HashCode.

     * @return the hashCode for this Pnt

     */

    public int hashCode () {

        int hash = 0;

        for (int i = 0; i < this.coordinates.length; i++) {

            long bits = Double.doubleToLongBits(this.coordinates[i]);

            hash = (31*hash) ^ (int)(bits ^ (bits >> 32));

        }

        return hash;

    }

    /* Pnts as vectors */

    /**

     * @return the specified coordinate of this Pnt

     * @throws ArrayIndexOutOfBoundsException for bad coordinate

     */

    public double coord (int i) {

        return this.coordinates[i];

    }

    /**

     * @return this Pnt's dimension.

     */

    public int dimension () {

        return coordinates.length;

    }

    /**

     * Check that dimensions match.

     * @param p the Pnt to check (against this Pnt)

     * @return the dimension of the Pnts

     * @throws IllegalArgumentException if dimension fail to match

     */

    public int dimCheck (Pnt p) {

        int len = this.coordinates.length;

        if (len != p.coordinates.length)

            throw new IllegalArgumentException("Dimension mismatch");

        return len;

    }

    /**

     * Create a new Pnt by adding additional coordinates to this Pnt.

     * @param coords the new coordinates (added on the right end)

     * @return a new Pnt with the additional coordinates

     */

    public Pnt extend (double[] coords) {

        double[] result = new double[coordinates.length + coords.length];

        System.arraycopy(coordinates, 0, result, 0, coordinates.length);

        System.arraycopy(coords, 0, result, coordinates.length, coords.length);

        return new Pnt(result);

    }

    /**

     * Dot product.

     * @param p the other Pnt

     * @return dot product of this Pnt and p

     */

    public double dot (Pnt p) {

        int len = dimCheck(p);

        double sum = 0;

        for (int i = 0; i < len; i++)

            sum += this.coordinates[i] * p.coordinates[i];

        return sum;

    }

    /**

     * Magnitude (as a vector).

     * @return the Euclidean length of this vector

     */

    public double magnitude () {

        return Math.sqrt(this.dot(this));

    }

    /**

     * Subtract.

     * @param p the other Pnt

     * @return a new Pnt = this - p

     */

    public Pnt subtract (Pnt p) {

        int len = dimCheck(p);

        double[] coords = new double[len];

        for (int i = 0; i < len; i++)

            coords[i] = this.coordinates[i] - p.coordinates[i];

        return new Pnt(coords);

    }

    /**

     * Add.

     * @param p the other Pnt

     * @return a new Pnt = this + p

     */

    public Pnt add (Pnt p) {

        int len = dimCheck(p);

        double[] coords = new double[len];

        for (int i = 0; i < len; i++)

            coords[i] = this.coordinates[i] + p.coordinates[i];

        return new Pnt(coords);

    }

    /**

     * Angle (in radians) between two Pnts (treated as vectors).

     * @param p the other Pnt

     * @return the angle (in radians) between the two Pnts

     */

    public double angle (Pnt p) {

        return Math.acos(this.dot(p) / (this.magnitude() * p.magnitude()));

    }

    /**

     * Perpendicular bisector of two Pnts.

     * Works in any dimension.  The coefficients are returned as a Pnt of one

     * higher dimension (e.g., (A,B,C,D) for an equation of the form

     * Ax + By + Cz + D = 0).

     * @param point the other point

     * @return the coefficients of the perpendicular bisector

     */

    public Pnt bisector (Pnt point) {

        int dim = dimCheck(point);

        Pnt diff = this.subtract(point);

        Pnt sum = this.add(point);

        double dot = diff.dot(sum);

        return diff.extend(new double[] {-dot / 2});

    }

    /* Pnts as matrices */

    /**

     * Create a String for a matrix.

     * @param matrix the matrix (an array of Pnts)

     * @return a String represenation of the matrix

     */

    public static String toString (Pnt[] matrix) {

        StringBuffer buf = new StringBuffer("{");

        for (int i = 0; i < matrix.length; i++) buf.append(" " + matrix[i]);

        buf.append(" }");

        return buf.toString();

    }

    /**

     * Compute the determinant of a matrix (array of Pnts).

     * This is not an efficient implementation, but should be adequate 

     * for low dimension.

     * @param matrix the matrix as an array of Pnts

     * @return the determinnant of the input matrix

     * @throws IllegalArgumentException if dimensions are wrong

     */

    public static double determinant (Pnt[] matrix) {

        if (matrix.length != matrix[0].dimension())

            throw new IllegalArgumentException("Matrix is not square");

        boolean[] columns = new boolean[matrix.length];

        for (int i = 0; i < matrix.length; i++) columns[i] = true;

        try {return determinant(matrix, 0, columns);}

        catch (ArrayIndexOutOfBoundsException e) {

            throw new IllegalArgumentException("Matrix is wrong shape");

        }

    }

    /**

     * Compute the determinant of a submatrix specified by starting row

     * and by "active" columns.

     * @param matrix the matrix as an array of Pnts

     * @param row the starting row

     * @param columns a boolean array indicating the "active" columns

     * @return the determinant of the specified submatrix

     * @throws ArrayIndexOutOfBoundsException if dimensions are wrong

     */

    private static double determinant(Pnt[] matrix, int row, boolean[] columns) {

        if (row == matrix.length) return 1;

        double sum = 0;

        int sign = 1;

        for (int col = 0; col < columns.length; col++) {

            if (!columns[col]) continue;

            columns[col] = false;

            sum += sign * matrix[row].coordinates[col] *

                   determinant(matrix, row+1, columns);

            columns[col] = true;

            sign = -sign;

        }

        return sum;

    }

    /**

     * Compute generalized cross-product of the rows of a matrix.

     * The result is a Pnt perpendicular (as a vector) to each row of

     * the matrix.  This is not an efficient implementation, but should 

     * be adequate for low dimension.

     * @param matrix the matrix of Pnts (one less row than the Pnt dimension)

     * @return a Pnt perpendicular to each row Pnt

     * @throws IllegalArgumentException if matrix is wrong shape

     */

    public static Pnt cross (Pnt[] matrix) {

        int len = matrix.length + 1;

        if (len != matrix[0].dimension())

            throw new IllegalArgumentException("Dimension mismatch");

        boolean[] columns = new boolean[len];

        for (int i = 0; i < len; i++) columns[i] = true;

        double[] result = new double[len];

        int sign = 1;

        try {

            for (int i = 0; i < len; i++) {

                columns[i] = false;

                result[i] = sign * determinant(matrix, 0, columns);

                columns[i] = true;

                sign = -sign;

            }

        } catch (ArrayIndexOutOfBoundsException e) {

            throw new IllegalArgumentException("Matrix is wrong shape");

        }

        return new Pnt(result);

    }

    /* Pnts as simplices */

    /**

     * Determine the signed content (i.e., area or volume, etc.) of a simplex.

     * @param simplex the simplex (as an array of Pnts)

     * @return the signed content of the simplex

     */

    public static double content (Pnt[] simplex) {

        Pnt[] matrix = new Pnt[simplex.length];

        for (int i = 0; i < matrix.length; i++)

            matrix[i] = simplex[i].extend(new double[] {1});

        int fact = 1;

        for (int i = 1; i < matrix.length; i++) fact = fact*i;

        return determinant(matrix) / fact;

    }

    /**

     * Relation between this Pnt and a simplex (represented as an array of Pnts).

     * Result is an array of signs, one for each vertex of the simplex, indicating

     * the relation between the vertex, the vertex's opposite facet, and this

     * Pnt. <pre>

     *   -1 means Pnt is on same side of facet

     *    0 means Pnt is on the facet

     *   +1 means Pnt is on opposite side of facet</pre>

     * @param simplex an array of Pnts representing a simplex

     * @return an array of signs showing relation between this Pnt and the simplex

     * @throws IllegalArgumentExcpetion if the simplex is degenerate

     */

    public int[] relation (Pnt[] simplex) {

        /* In 2D, we compute the cross of this matrix:

         *    1   1   1   1

         *    p0  a0  b0  c0

         *    p1  a1  b1  c1

         * where (a, b, c) is the simplex and p is this Pnt.  The result

         * is a vector in which the first coordinate is the signed area

         * (all signed areas are off by the same constant factor) of

         * the simplex and the remaining coordinates are the *negated*

         * signed areas for the simplices in which p is substituted for

         * each of the vertices. Analogous results occur in higher dimensions.

         */

        int dim = simplex.length - 1;

        if (this.dimension() != dim)

            throw new IllegalArgumentException("Dimension mismatch");

        /* Create and load the matrix */

        Pnt[] matrix = new Pnt[dim+1];

        /* First row */

        double[] coords = new double[dim+2];

        for (int j = 0; j < coords.length; j++) coords[j] = 1;

        matrix[0] = new Pnt(coords);

        /* Other rows */

        for (int i = 0; i < dim; i++) {

            coords[0] = this.coordinates[i];

            for (int j = 0; j < simplex.length; j++) 

                coords[j+1] = simplex[j].coordinates[i];

            matrix[i+1] = new Pnt(coords);

        }

        /* Compute and analyze the vector of areas/volumes/contents */

        Pnt vector = cross(matrix);

        double content = vector.coordinates[0];

        int[] result = new int[dim+1];

        for (int i = 0; i < result.length; i++) {

            double value = vector.coordinates[i+1];

            if (Math.abs(value) <= 1.0e-6 * Math.abs(content)) result[i] = 0;

            else if (value < 0) result[i] = -1;

            else result[i] = 1;

        }

        if (content < 0) {

            for (int i = 0; i < result.length; i++) result[i] = -result[i];

        }

        if (content == 0) {

            for (int i = 0; i < result.length; i++) result[i] = Math.abs(result[i]);

        }

        return result;

    }

    /**

     * Test if this Pnt is outside of simplex.

     * @param simplex the simplex (an array of Pnts)

     * @return the simplex Pnt that "witnesses" outsideness (or null if not outside)

     */

    public Pnt isOutside (Pnt[] simplex) {

        int[] result = this.relation(simplex);

        for (int i = 0; i < result.length; i++) {

            if (result[i] > 0) return simplex[i];

        }

        return null;

    }

    /**

     * Test if this Pnt is on a simplex.

     * @param simplex the simplex (an array of Pnts)

     * @return the simplex Pnt that "witnesses" on-ness (or null if not on)

     */

    public Pnt isOn (Pnt[] simplex) {

        int[] result = this.relation(simplex);

        Pnt witness = null;

        for (int i = 0; i < result.length; i++) {

            if (result[i] == 0) witness = simplex[i];

            else if (result[i] > 0) return null;

        }

        return witness;

    }

    /**

     * Test if this Pnt is inside a simplex.

     * @param simplex the simplex (an arary of Pnts)

     * @return true iff this Pnt is inside simplex.

     */

    public boolean isInside (Pnt[] simplex) {

        int[] result = this.relation(simplex);

        for (int i = 0; i < result.length; i++) if (result[i] >= 0) return false;

        return true;

    }

    /**

     * Test relation between this Pnt and circumcircle of a simplex.

     * @param simplex the simplex (as an array of Pnts)

     * @return -1, 0, or +1 for inside, on, or outside of circumcircle

     */

    public int vsCircumcircle (Pnt[] simplex) {

        Pnt[] matrix = new Pnt[simplex.length + 1];

        for (int i = 0; i < simplex.length; i++)

            matrix[i] = simplex[i].extend(new double[] {1, simplex[i].dot(simplex[i])});

        matrix[simplex.length] = this.extend(new double[] {1, this.dot(this)});

        double d = determinant(matrix);

        int result = (d < 0)? -1 : ((d > 0)? +1 : 0);

        if (content(simplex) < 0) result = - result;

        return result;

    }

    /**

     * Circumcenter of a simplex.

     * @param simplex the simplex (as an array of Pnts)

     * @return the circumcenter (a Pnt) of simplex

     */

    public static Pnt circumcenter (Pnt[] simplex) {

        int dim = simplex[0].dimension();

        if (simplex.length - 1 != dim)

            throw new IllegalArgumentException("Dimension mismatch");

        Pnt[] matrix = new Pnt[dim];

        for (int i = 0; i < dim; i++) 

            matrix[i] = simplex[i].bisector(simplex[i+1]);

        Pnt hCenter = cross(matrix);      // Center in homogeneous coordinates

        double last = hCenter.coordinates[dim];

        double[] result = new double[dim];

        for (int i = 0; i < dim; i++) result[i] = hCenter.coordinates[i] / last;

        return new Pnt(result);

    }

    /**

     * Main program (used for testing).

     */

    public static void main (String[] args) {

        Pnt p = new Pnt(1, 2, 3);

        System.out.println("Pnt created: " + p);

        Pnt[] matrix1 = {new Pnt(1,2), new Pnt(3,4)};

        Pnt[] matrix2 = {new Pnt(7,0,5), new Pnt(2,4,6), new Pnt(3,8,1)};

        System.out.print("Results should be -2 and -288: ");

        System.out.println(determinant(matrix1) + " " + determinant(matrix2));

        Pnt p1 = new Pnt(1,1); Pnt p2 = new Pnt(-1,1);

        System.out.println("Angle between " + p1 + " and " + p2 + ": " + p1.angle(p2));

        System.out.println(p1 + " subtract " + p2 + ": " + p1.subtract(p2));

        Pnt v0 = new Pnt(0,0), v1 = new Pnt(1,1), v2 = new Pnt(2,2);

        Pnt[] vs = {v0, new Pnt(0,1), new Pnt(1,0)};

        Pnt vp = new Pnt(.1, .1);

        System.out.println(vp + " isInside " + toString(vs) + ": " + vp.isInside(vs));

        System.out.println(v1 + " isInside " + toString(vs) + ": " + v1.isInside(vs));

        System.out.println(vp + " vsCircumcircle " + toString(vs) + ": " +

                           vp.vsCircumcircle(vs));

        System.out.println(v1 + " vsCircumcircle " + toString(vs) + ": " +

                           v1.vsCircumcircle(vs));

        System.out.println(v2 + " vsCircumcircle " + toString(vs) + ": " +

                           v2.vsCircumcircle(vs));

        System.out.println("Circumcenter of " + toString(vs) + " is " + circumcenter(vs));

    }

}