    001    //$HeadURL$
    002    /*----------------------------------------------------------------------------
    003     This file is part of deegree, http://deegree.org/
    004     Copyright (C) 2001-2009 by:
    005     Department of Geography, University of Bonn
    006     and
    007     lat/lon GmbH
    008
    009     This library is free software; you can redistribute it and/or modify it under
    010     the terms of the GNU Lesser General Public License as published by the Free
    011     Software Foundation; either version 2.1 of the License, or (at your option)
    012     any later version.
    013     This library is distributed in the hope that it will be useful, but WITHOUT
    014     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
    015     FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
    016     details.
    017     You should have received a copy of the GNU Lesser General Public License
    018     along with this library; if not, write to the Free Software Foundation, Inc.,
    019     59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
    020
    021     Contact information:
    022
    023     lat/lon GmbH
    024     Aennchenstr. 19, 53177 Bonn
    025     Germany
    026     http://lat-lon.de/
    027
    028     Department of Geography, University of Bonn
    029     Prof. Dr. Klaus Greve
    030     Postfach 1147, 53001 Bonn
    031     Germany
    032     http://www.geographie.uni-bonn.de/deegree/
    033
    034     e-mail: info@deegree.org
    035     ----------------------------------------------------------------------------*/
    036
    037    package org.deegree.model.spatialschema;
    038
    039    import javax.vecmath.Point3d;
    040    import javax.vecmath.Vector2d;
    041    import javax.vecmath.Vector3d;
    042
    043    import org.deegree.framework.log.ILogger;
    044    import org.deegree.framework.log.LoggerFactory;
    045
    046    /**
    047     * Utility class for the linearization of arcs and circles.
    048     *
    049     * @author <a href="mailto:elmasry@lat-lon.de">Moataz Elmasry</a>
    050     * @author last edited by: $Author: elmasri$
    051     *
    052     * @version $Revision: $, $Date: 9 May 2008 13:09:29$
    053     */
    054    public class LinearizationUtil {
    055
    056        private static ILogger LOG = LoggerFactory.getLogger( LinearizationUtil.class );
    057
    058        private static double PI2 = Math.PI * 2;
    059
    060        private static double EPSILON = 1E-11;
    061
    062        /**
    063         * Takes in a three Position of an arc and computes a linearized Arc (linestring) using 150 positions/vertices
    064         *
    065         * @param p0
    066         * @param p1
    067         * @param p2
    068         * @return An array of linearized Positions
    069         */
    070        public static Position[] linearizeArc( Position p0, Position p1, Position p2 ) {
    071            return linearizeArc( p0, p1, p2, 150 );
    072        }
    073
    074        /**
    075         * Takes in a three Positions of an Arc and computes a linearized Arc (linestring).
    076         *
    077         * @param p0
    078         * @param p1
    079         * @param p2
    080         *
    081         * @param numberOfOpPositions
    082         *            the number of linearized Positions to produce
    083         * @return An array of linearized Positions
    084         */
    085        public static Position[] linearizeArc( Position p0, Position p1, Position p2, int numberOfOpPositions ) {
    086
    087            // shift the points down (to reduce the occurrence of floating point errors), independently on the x and y axes
    088            double shiftx = findShiftX( p0, p1, p2 );
    089            double shifty = findShiftY( p0, p1, p2 );
    090            // if the points are already shifted, this does no harm!
    091            Position p0Shifted = new PositionImpl( p0.getX() - shiftx, p0.getY() - shifty );
    092            Position p1Shifted = new PositionImpl( p1.getX() - shiftx, p1.getY() - shifty );
    093            Position p2Shifted = new PositionImpl( p2.getX() - shiftx, p2.getY() - shifty );
    094
    095            if ( areCollinear( p0Shifted, p1Shifted, p2Shifted ) ) {
    096                Position[] Positions = new Position[3];
    097                Positions[0] = p0;
    098                Positions[1] = p1;
    099                Positions[2] = p2;
    100                return Positions;
    101            }
    102
    103            int j = numberOfOpPositions - 1;
    104            if ( j < 1 ) {
    105                return null;
    106            }
    107
    108            double[] positionsDouble = computeArc( p0Shifted.getX(), p0Shifted.getY(), p1Shifted.getX(), p1Shifted.getY(),
    109                                                   p2Shifted.getX(), p2Shifted.getY() );
    110            double d11 = computeArcOrientation( p0Shifted.getX(), p0Shifted.getY(), p1Shifted.getX(), p1Shifted.getY(),
    111                                                p2Shifted.getX(), p2Shifted.getY() );
    112            double d12 = positionsDouble[3];
    113            double d13 = positionsDouble[5];
    114            if ( d11 > 0 && d13 < d12 ) {
    115                d13 += PI2;
    116            } else if ( d11 < 0 && d13 > d12 ) {
    117                d13 -= PI2;
    118            }
    119            Position outPositions[] = new Position[numberOfOpPositions];
    120            for ( int k = 0; k <= j; k++ ) {
    121                double d15 = d12 + ( ( d13 - d12 ) * k ) / j;
    122                double xCoord = positionsDouble[0] + positionsDouble[2] * Math.cos( d15 );
    123                double yCoord = positionsDouble[1] + positionsDouble[2] * Math.sin( d15 );
    124                outPositions[k] = GeometryFactory.createPosition( xCoord, yCoord );
    125            }
    126
    127            // ensure numerical stability for start and end points
    128            outPositions[0] = GeometryFactory.createPosition( p0Shifted.getX(), p0Shifted.getY() );
    129            outPositions[outPositions.length - 1] = GeometryFactory.createPosition( p2Shifted.getX(), p2Shifted.getY() );
    130
    131            // shift the points back up
    132            for ( int i = 0; i < outPositions.length; i++ ) {
    133                outPositions[i] = new PositionImpl( outPositions[i].getX() + shiftx, outPositions[i].getY() + shifty );
    134            }
    135
    136            return outPositions;
    137        }
    138
    139        private static double findShiftY( Position p0, Position p1, Position p2 ) {
    140            double miny = p0.getY();
    141            if ( p1.getY() < miny ) {
    142                miny = p1.getY();
    143            }
    144            if ( p2.getY() < miny ) {
    145                miny = p2.getY();
    146            }
    147
    148            double maxy = p0.getY();
    149            if ( p1.getX() > maxy ) {
    150                maxy = p1.getY();
    151            }
    152            if ( p2.getY() > maxy ) {
    153                maxy = p2.getY();
    154            }
    155            return ( miny + maxy ) / 2;
    156        }
    157
    158        private static double findShiftX( Position p0, Position p1, Position p2 ) {
    159            double minx = p0.getX();
    160            if ( p1.getX() < minx ) {
    161                minx = p1.getX();
    162            }
    163            if ( p2.getX() < minx ) {
    164                minx = p2.getX();
    165            }
    166
    167            double maxx = p0.getX();
    168            if ( p1.getX() > maxx ) {
    169                maxx = p1.getX();
    170            }
    171            if ( p2.getX() > maxx ) {
    172                maxx = p2.getX();
    173            }
    174            return ( minx + maxx ) / 2;
    175        }
    176
    177        /**
    178         * Calculates a sequence of {@link Position}s on the arc of a circle.
    179         * <p>
    180         * The circle is constructed from the input points: All three points belong to the arc of the circle. They must be
    181         * distinct and non-colinear. To form a complete circle, the arc is extended past the third control point until the
    182         * first control point is encountered.
    183         *
    184         * @param p0
    185         *            start point on the arc of the circle
    186         * @param p1
    187         *            second point on the arc of the circle
    188         * @param p2
    189         *            third point on the arc of the circle
    190         * @param numberOfPositions
    191         *            number of interpolation points, the returned array contains numberOfPositions+1 entries, to ensure
    192         *            that the first point is identical to the last
    193         * @return Position array with numberOfPositions+1 entries, first entry and last entry are identical to p0
    194         * @throws IllegalArgumentException
    195         *             if no order can be determined, because the points are identical or co-linear
    196         */
    197        public static Position[] linearizeCircle( Position p0, Position p1, Position p2, int numberOfPositions ) {
    198
    199            // shift the points down (to reduce the occurrence of floating point errors), independently on the x and y axes
    200            double shiftx = findShiftX( p0, p1, p2 );
    201            double shifty = findShiftY( p0, p1, p2 );
    202            // if the points are already shifted, this does no harm!
    203            Position p0Shifted = new PositionImpl( p0.getX() - shiftx, p0.getY() - shifty );
    204            Position p1Shifted = new PositionImpl( p1.getX() - shiftx, p1.getY() - shifty );
    205            Position p2Shifted = new PositionImpl( p2.getX() - shiftx, p2.getY() - shifty );
    206
    207            Position[] positions = new Position[numberOfPositions + 1];
    208
    209            Position center = findCircleCenter( p0Shifted, p1Shifted, p2Shifted );
    210
    211            double centerX = center.getX();
    212            double centerY = center.getY();
    213            double dx = p0Shifted.getX() - centerX;
    214            double dy = p0Shifted.getY() - centerY;
    215            double radius = Math.sqrt( dx * dx + dy * dy );
    216
    217            double angleStep = 0.0;
    218            if ( isClockwise( p0Shifted, p1Shifted, p2Shifted ) ) {
    219                angleStep = -Math.PI * 2 / numberOfPositions;
    220            } else {
    221                angleStep = Math.PI * 2 / numberOfPositions;
    222            }
    223            double angle = Math.atan2( dy, dx );
    224
    225            for ( int i = 0; i < numberOfPositions; i++ ) {
    226                double x = centerX + Math.cos( angle ) * radius;
    227                double y = centerY + Math.sin( angle ) * radius;
    228                positions[i] = GeometryFactory.createPosition( x, y );
    229                angle += angleStep;
    230            }
    231            positions[numberOfPositions] = positions[0];
    232
    233            // shift the points back up
    234            for ( int i = 0; i <= numberOfPositions; i++ ) {
    235                positions[i] = new PositionImpl( positions[i].getX() + shiftx, positions[i].getY() + shifty );
    236            }
    237
    238            return positions;
    239        }
    240
    241        /**
    242         * Returns whether the order of the given three points is clockwise or counterclockwise. Uses the (signed) area of a
    243         * planar triangle to get to know about the order of the points.
    244         *
    245         * @param p0
    246         *            first point
    247         * @param p1
    248         *            second point
    249         * @param p2
    250         *            third point
    251         * @return true, if order is clockwise, otherwise false
    252         * @throws IllegalArgumentException
    253         *             if no order can be determined, because the points are identical or co-linear
    254         */
    255        static final boolean isClockwise( Position p0, Position p1, Position p2 )
    256                                throws IllegalArgumentException {
    257
    258            if ( areCollinear( p0, p1, p2 ) ) {
    259                throw new IllegalArgumentException( "Cannot evaluate isClockwise(). The three points are colinear." );
    260            }
    261            double res = ( p2.getX() - p0.getX() ) * ( ( p2.getY() + p0.getY() ) / 2 ) + ( p1.getX() - p2.getX() )
    262                         * ( ( p1.getY() + p2.getY() ) / 2 ) + ( p0.getX() - p1.getX() ) * ( ( p0.getY() + p1.getY() ) / 2 );
    263            return res < 0.0 ? true : false;
    264        }
    265
    266        /**
    267         * Tests if the given three points are collinear.
    268         * <p>
    269         * NOTE: Only this method should be used throughout the whole linearization process for testing collinearity to
    270         * avoid inconsistent results (the necessary EPSILON would differ).
    271         *
    272         * @param p0
    273         * @param p1
    274         * @param p2
    275         * @return true if the points are collinear, false otherwise
    276         */
    277        public static boolean areCollinear( Position p0, Position p1, Position p2 ) {
    278            double res = ( p2.getX() - p0.getX() ) * ( ( p2.getY() + p0.getY() ) / 2 ) + ( p1.getX() - p2.getX() )
    279                         * ( ( p1.getY() + p2.getY() ) / 2 ) + ( p0.getX() - p1.getX() ) * ( ( p0.getY() + p1.getY() ) / 2 );
    280            return Math.abs( res ) < EPSILON;
    281        }
    282
    283        /**
    284         * Computes an Arc from of three points (wich must not be collinear).
    285         *
    286         * @return array of points representing the computed arc
    287         */
    288        private static final double[] computeArc( double d1, double d2, double d3, double d4, double d5, double d6 ) {
    289
    290            double d13 = d1 - d3;
    291            double d14 = d3 - d5;
    292            double d15 = d4 - d6;
    293            double d16 = d2 - d4;
    294            double d17 = d1 + d3;
    295            double d18 = d3 + d5;
    296            double d19 = d2 + d4;
    297            double d20 = d4 + d6;
    298            double d21 = d13 * d15 - d14 * d16;
    299
    300            double d7 = ( ( d15 * d17 * d13 - d16 * d18 * d14 ) + d16 * d15 * ( d16 + d15 ) ) / d21;
    301            double d8;
    302            if ( Math.abs( d4 - d6 ) < 5E-008D ) {
    303                d8 = ( d19 * d16 + ( d17 - d7 ) * d13 ) / d16;
    304            } else {
    305                d8 = ( d20 * d15 + ( d18 - d7 ) * d14 ) / d15;
    306            }
    307            d7 *= 0.5;
    308            d8 *= 0.5;
    309            double d9 = Math.sqrt( ( d7 - d1 ) * ( d7 - d1 ) + ( d8 - d2 ) * ( d8 - d2 ) );
    310            double d10 = Math.atan2( d2 - d8, d1 - d7 );
    311            if ( d10 < 0 ) {
    312                d10 += PI2;
    313            }
    314            double d11 = Math.atan2( d4 - d8, d3 - d7 );
    315            if ( d11 < 0 ) {
    316                d11 += PI2;
    317            }
    318            double d12 = Math.atan2( d6 - d8, d5 - d7 );
    319            if ( d12 < 0 ) {
    320                d12 += PI2;
    321            }
    322
    323            double ad[] = new double[6];
    324            ad[0] = d7;
    325            ad[1] = d8;
    326            ad[2] = d9;
    327            ad[3] = d10;
    328            ad[4] = d11;
    329            ad[5] = d12;
    330            return ad;
    331
    332        }
    333
    334        private static final double computeArcOrientation( double d, double d1, double d2, double d3, double d4, double d5 ) {
    335            return ( d * d3 + d2 * d5 + d4 * d1 ) - ( d4 * d3 + d2 * d1 + d * d5 );
    336        }
    337
    338        /**
    339         * Transforms an array of Position to an array of double
    340         *
    341         * @param positions
    342         * @return array of double
    343         */
    344        protected static double[] positionToDouble( Position[] positions ) {
    345            double[] ad = new double[positions.length * 2];
    346            for ( int i = 0; i < positions.length; i++ ) {
    347                ad[i] = positions[i].getX();
    348                ad[i + 1] = positions[i].getY();
    349            }
    350            return ad;
    351        }
    352
    353        /**
    354         * Transforms an array of double to an array of Position
    355         *
    356         * @param doubles
    357         * @return generated array of Position or null if doubles%2 != 0
    358         */
    359        protected static Position[] doubleToPosition( double[] doubles ) {
    360            // The number of doubles has to be even
    361            if ( doubles.length % 2 != 0 ) {
    362                return null;
    363            }
    364            Position[] positions = new Position[doubles.length / 2];
    365            for ( int i = 0; i < positions.length; i++ ) {
    366                positions[i] = GeometryFactory.createPosition( doubles[i], doubles[i + 1] );
    367            }
    368            return positions;
    369        }
    370
    371        /**
    372         * Finds the center of a cirlce/arc using three points that lie on the circle. Thanks to wikipedia:
    373         * http://en.wikipedia.org/wiki/Circumradius#Coordinates_of_circumcenter.
    374         *
    375         * @param p0
    376         * @param p1
    377         * @param p2
    378         * @return center of the circle
    379         * @throws IllegalArgumentException
    380         *             if the points are co linear, e.g. on a line.
    381         */
    382        static Position findCircleCenter( Position p0, Position p1, Position p2 )
    383                                throws IllegalArgumentException {
    384
    385            // shift the points down (to reduce the occurrence of floating point errors), independently on the x and y axes
    386            double shiftx = findShiftX( p0, p1, p2 );
    387            double shifty = findShiftY( p0, p1, p2 );
    388            // if the points are already shifted, this does no harm!
    389            Position p0Shifted = new PositionImpl( p0.getX() - shiftx, p0.getY() - shifty );
    390            Position p1Shifted = new PositionImpl( p1.getX() - shiftx, p1.getY() - shifty );
    391            Position p2Shifted = new PositionImpl( p2.getX() - shiftx, p2.getY() - shifty );
    392
    393            if ( areCollinear( p0Shifted, p1Shifted, p2Shifted ) ) {
    394                throw new IllegalArgumentException( "The given points are co linear, no circum center can be calculated." );
    395            }
    396
    397            Vector3d a = new Vector3d( p0Shifted.getAsPoint3d() );
    398            Vector3d b = new Vector3d( p1Shifted.getAsPoint3d() );
    399            Vector3d c = new Vector3d( p2Shifted.getAsPoint3d() );
    400
    401            if ( Double.isNaN( a.z ) ) {
    402                a.z = 0;
    403            }
    404            if ( Double.isNaN( b.z ) ) {
    405                b.z = 0;
    406            }
    407            if ( Double.isNaN( c.z ) ) {
    408                c.z = 0;
    409            }
    410
    411            Vector3d ab = new Vector3d( a );
    412            Vector3d ac = new Vector3d( a );
    413            Vector3d bc = new Vector3d( b );
    414            Vector3d ba = new Vector3d( b );
    415            Vector3d ca = new Vector3d( c );
    416            Vector3d cb = new Vector3d( c );
    417
    418            ab.sub( b );
    419            ac.sub( c );
    420            bc.sub( c );
    421            ba.sub( a );
    422            ca.sub( a );
    423            cb.sub( b );
    424
    425            Vector3d cros = new Vector3d();
    426            cros.cross( ab, bc );
    427            double crosSquare = 2 * cros.length() * cros.length();
    428
    429            if ( LOG.isDebug() ) {
    430                double radius = ( ab.length() * bc.length() * ca.length() ) / ( 2 * cros.length() );
    431                LOG.logDebug( "Radius: " + radius );
    432            }
    433
    434            a.scale( ( ( bc.length() * bc.length() ) * ab.dot( ac ) ) / crosSquare );
    435            b.scale( ( ( ac.length() * ac.length() ) * ba.dot( bc ) ) / crosSquare );
    436            c.scale( ( ( ab.length() * ab.length() ) * ca.dot( cb ) ) / crosSquare );
    437
    438            Point3d circle = new Point3d( a );
    439            circle.add( b );
    440            circle.add( c );
    441
    442            // shift the center circle back up
    443            circle.x += shiftx;
    444            circle.y += shifty;
    445
    446            return GeometryFactory.createPosition( circle );
    447        }
    448
    449        /**
    450         * Finds the angle between three points, p0,p1 and a center, where the angle to be measured. In other words, its the
    451         * angle between the segments p0-center and p1-center.
    452         *
    453         * @param p0
    454         * @param p1
    455         * @param centre
    456         * @return calculated angle in degree
    457         */
    458        private static double findAngle( Position p0, Position p1, Position centre ) {
    459            Vector2d vec1 = calcVector( centre, p0 );
    460            Vector2d vec2 = calcVector( centre, p1 );
    461            double cosTheta = ( vec1.dot( vec2 ) ) / ( vec1.length() * vec2.length() );
    462            int round = 1000;
    463            // Using the 1000th decimal to round the cosTheta.
    464            cosTheta = Math.ceil( cosTheta * round ) / round;
    465            // convert from radian to degree
    466            return Math.acos( cosTheta ) / Math.PI * 180;
    467
    468        }
    469
    470        /**
    471         * Finds the vector between two points. which is the difference between p1 and p0 and the direction is from p0 to
    472         * p1.
    473         *
    474         * @param p0
    475         * @param p1
    476         * @return calculated vector
    477         */
    478        protected static Vector2d calcVector( Position p0, Position p1 ) {
    479            return new Vector2d( p1.getX() - p0.getX(), p1.getY() - p0.getY() );
    480        }
    481    }