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.framework.util;
038    
039    import static org.deegree.framework.log.LoggerFactory.getLogger;
040    
041    import static org.deegree.framework.util.CollectionUtils.map;
042    import static org.deegree.model.spatialschema.GeometryFactory.createCurve;
043    import static org.deegree.model.spatialschema.GeometryFactory.createCurveSegment;
044    import static org.deegree.model.spatialschema.GeometryFactory.createEnvelope;
045    import static org.deegree.model.spatialschema.GeometryFactory.createPoint;
046    import static org.deegree.model.spatialschema.GeometryFactory.createPosition;
047    import static org.deegree.model.spatialschema.GeometryFactory.createSurface;
048    import static org.deegree.model.spatialschema.GeometryFactory.createSurfacePatch;
049    
050    import java.util.ArrayList;
051    import java.util.List;
052    
053    import org.deegree.framework.log.ILogger;
054    import org.deegree.framework.util.CollectionUtils.Mapper;
055    import org.deegree.model.crs.CoordinateSystem;
056    import org.deegree.model.spatialschema.Aggregate;
057    import org.deegree.model.spatialschema.Curve;
058    import org.deegree.model.spatialschema.CurveSegment;
059    import org.deegree.model.spatialschema.Envelope;
060    import org.deegree.model.spatialschema.Geometry;
061    import org.deegree.model.spatialschema.GeometryException;
062    import org.deegree.model.spatialschema.GeometryFactory;
063    import org.deegree.model.spatialschema.JTSAdapter;
064    import org.deegree.model.spatialschema.MultiSurface;
065    import org.deegree.model.spatialschema.Point;
066    import org.deegree.model.spatialschema.Position;
067    import org.deegree.model.spatialschema.Ring;
068    import org.deegree.model.spatialschema.Surface;
069    import org.deegree.model.spatialschema.SurfaceInterpolationImpl;
070    import org.deegree.model.spatialschema.SurfacePatch;
071    
072    import com.vividsolutions.jts.algorithm.CGAlgorithms;
073    
074    /**
075     * 
076     * 
077     * @author <a href="mailto:poth@lat-lon.de">Andreas Poth</a>
078     * @author last edited by: $Author: poth $
079     * 
080     * @version $Revision: 6251 $, $Date: 2007-03-19 16:59:28 +0100 (Mo, 19 Mrz 2007) $
081     */
082    public class GeometryUtils {
083    
084        private static final ILogger LOG = getLogger( GeometryUtils.class );
085    
086        /**
087         * 
088         * @param p1
089         * @param p2
090         * @return distance between two 2D positions
091         */
092        public static double distance( Position p1, Position p2 ) {
093            return distance( p1.getX(), p1.getY(), p2.getX(), p2.getY() );
094        }
095    
096        /**
097         * 
098         * @param xx
099         * @param yy
100         * @param xx_
101         * @param yy_
102         * @return distance between two points
103         */
104        public static double distance( double xx, double yy, double xx_, double yy_ ) {
105            double dx = xx - xx_;
106            double dy = yy - yy_;
107            return Math.sqrt( dx * dx + dy * dy );
108        }
109    
110        /**
111         * @param <T>
112         * @param geom
113         * @param x
114         * @param y
115         * @return moves 2D geometries only (will discard any z values
116         */
117        @SuppressWarnings("unchecked")
118        public static <T> T move( T geom, double x, double y ) {
119            final Mapper<Position, Position> mover = getMover( x, y );
120            if ( geom instanceof Envelope ) {
121                Envelope e = (Envelope) geom;
122                // dammit, where's the proper typecase?
123                return (T) createEnvelope( move( e.getMin(), x, y ), move( e.getMax(), x, y ), e.getCoordinateSystem() );
124            }
125            if ( geom instanceof Point ) {
126                Point p = (Point) geom;
127                return (T) createPoint( p.getX() + x, p.getY() + y, p.getCoordinateSystem() );
128            }
129            if ( geom instanceof Curve ) {
130                Curve c = (Curve) geom;
131                try {
132                    return (T) createCurve( map( c.getCurveSegments(), GeometryUtils.<CurveSegment> getMover( x, y ) ),
133                                            c.getCoordinateSystem() );
134                } catch ( GeometryException e ) {
135                    LOG.logError( "Unknown error", e );
136                    return null;
137                }
138            }
139            if ( geom instanceof CurveSegment ) {
140                CurveSegment c = (CurveSegment) geom;
141                try {
142                    return (T) createCurveSegment( map( c.getPositions(), mover ), c.getCoordinateSystem() );
143                } catch ( GeometryException e ) {
144                    LOG.logError( "Unknown error", e );
145                    return null;
146                }
147            }
148            if ( geom instanceof Position ) {
149                Position p = (Position) geom;
150                return (T) createPosition( p.getX() + x, p.getY() + y );
151            }
152            if ( geom instanceof Surface ) {
153                Surface s = (Surface) geom;
154                SurfacePatch[] patches = new SurfacePatch[s.getNumberOfSurfacePatches()];
155                for ( int i = 0; i < patches.length; ++i ) {
156                    try {
157                        patches[i] = move( s.getSurfacePatchAt( 0 ), x, y );
158                    } catch ( GeometryException e ) {
159                        LOG.logError( "Unknown error", e );
160                        return null;
161                    }
162                }
163                try {
164                    return (T) createSurface( patches, s.getCoordinateSystem() );
165                } catch ( GeometryException e ) {
166                    LOG.logError( "Unknown error", e );
167                    return null;
168                }
169            }
170            if ( geom instanceof SurfacePatch ) {
171                SurfacePatch p = (SurfacePatch) geom;
172                Position[] ring = p.getExteriorRing();
173                Position[] outer = map( ring, mover ).toArray( new Position[ring.length] );
174                Position[][] inners = p.getInteriorRings();
175                if ( inners != null ) {
176                    for ( int i = 0; i < inners.length; ++i ) {
177                        inners[i] = map( inners[i], mover ).toArray( new Position[inners[i].length] );
178                    }
179                }
180                try {
181                    return (T) createSurfacePatch( outer, inners, p.getInterpolation(), p.getCoordinateSystem() );
182                } catch ( GeometryException e ) {
183                    LOG.logError( "Unknown error", e );
184                    return null;
185                }
186            }
187            if ( geom instanceof Aggregate ) {
188                Aggregate m = (Aggregate) geom;
189                for ( int i = 0; i < m.getSize(); ++i ) {
190                    try {
191                        m.setObjectAt( move( m.getObjectAt( i ), x, y ), i );
192                    } catch ( GeometryException e ) {
193                        LOG.logError( "Unknown error", e );
194                        return null;
195                    }
196                }
197                return (T) m;
198            }
199    
200            throw new UnsupportedOperationException( "Moving geometries of class " + geom.getClass() + " is not supported." );
201        }
202    
203        /**
204         * @param <T>
205         * @param x
206         * @param y
207         * @return a move mapper wrapper
208         */
209        public static <T> Mapper<T, T> getMover( final double x, final double y ) {
210            return new Mapper<T, T>() {
211                public T apply( T u ) {
212                    return move( u, x, y );
213                }
214            };
215        }
216    
217        /**
218         * 
219         * @param surface
220         * @return surface with inverted order of vertices
221         * @throws GeometryException
222         */
223        public static Surface invertOrder( Surface surface )
224                                throws GeometryException {
225            Position[] exring = surface.getSurfaceBoundary().getExteriorRing().getPositions();
226            // invert exterior ring vertices order
227            for ( int i = 0; i < exring.length / 2; i++ ) {
228                Position p = exring[i];
229                exring[i] = exring[exring.length - 1 - i];
230                exring[exring.length - 1 - i] = p;
231            }
232            Ring[] inner = surface.getSurfaceBoundary().getInteriorRings();
233            Position[][] inPos = new Position[inner.length][];
234            // invert interior ring vertices order
235            for ( int i = 0; i < inner.length; i++ ) {
236                Position[] tmp = inner[i].getPositions();
237                for ( int j = 0; j < tmp.length; j++ ) {
238                    Position p = tmp[j];
239                    tmp[j] = tmp[tmp.length - 1 - j];
240                    tmp[tmp.length - 1 - j] = p;
241                }
242                inPos[i] = tmp;
243            }
244            return GeometryFactory.createSurface( exring, inPos, new SurfaceInterpolationImpl(),
245                                                  surface.getCoordinateSystem() );
246        }
247    
248        /**
249         * 
250         * @param curve
251         * @return curve with inverted order of vertices for each segment
252         * @throws GeometryException
253         */
254        public static Curve invertOrder( Curve curve )
255                                throws GeometryException {
256            CurveSegment[] segments = curve.getCurveSegments();
257            for ( int i = 0; i < segments.length; i++ ) {
258                Position[] tmp = segments[i].getPositions();
259                for ( int j = 0; j < tmp.length; j++ ) {
260                    Position p = tmp[j];
261                    tmp[j] = tmp[tmp.length - 1 - j];
262                    tmp[tmp.length - 1 - j] = p;
263                }
264                segments[i] = GeometryFactory.createCurveSegment( tmp, curve.getCoordinateSystem() );
265            }
266            return GeometryFactory.createCurve( segments );
267        }
268    
269        /**
270         * 
271         * @param surface
272         * @return true if an array of passed {@link Position} forms a clockwise orientated ring
273         */
274        public static boolean isClockwise( Surface surface ) {
275            Position[] ring = surface.getSurfaceBoundary().getExteriorRing().getPositions();
276            return !CGAlgorithms.isCCW( JTSAdapter.export( ring ).getCoordinates() );
277        }
278    
279        /**
280         * 
281         * @param geom
282         * @return surface or multi surface with guaranteed clockwise vertices orientation
283         * @throws GeometryException
284         */
285        public static Geometry ensureClockwise( Geometry geom )
286                                throws GeometryException {
287            if ( geom instanceof Surface ) {
288                if ( !isClockwise( (Surface) geom ) ) {
289                    geom = invertOrder( (Surface) geom );
290                }
291            } else if ( geom instanceof MultiSurface ) {
292                Surface[] surfaces = ( (MultiSurface) geom ).getAllSurfaces();
293                for ( int i = 0; i < surfaces.length; i++ ) {
294                    surfaces[i] = invertOrder( surfaces[i] );
295                }
296                geom = GeometryFactory.createMultiSurface( surfaces );
297            }
298            return geom;
299        }
300    
301        /**
302         * 
303         * @param distance
304         *            if distance is < 0 left parallel will be created
305         * @param curve
306         * @return parallel curve with distance.
307         * @throws GeometryException
308         */
309        // TODO
310        // remove arte facts occuring if distance between at least two vertices is less passed distance
311        public static Curve createCurveParallel( double distance, Curve curve )
312                                throws GeometryException {
313            Position[] pos = curve.getAsLineString().getPositions();
314    
315            List<Position[]> posList = new ArrayList<Position[]>( pos.length );
316            // create parallel for each segment
317            for ( int j = 0; j < pos.length - 1; j++ ) {
318                // calculate normal vector
319                // swap x and y and change sign
320                double nx = -( pos[j].getY() - pos[j + 1].getY() );
321                double ny = pos[j].getX() - pos[j + 1].getX();
322                double nl = Math.sqrt( nx * nx + ny * ny );
323                nx = nx / nl * distance;
324                ny = ny / nl * distance;
325                Position[] p = new Position[2];
326                p[0] = GeometryFactory.createPosition( pos[j].getX() + nx, pos[j].getY() + ny );
327                p[1] = GeometryFactory.createPosition( pos[j + 1].getX() + nx, pos[j + 1].getY() + ny );
328                posList.add( p );
329            }
330            List<Position> pList = new ArrayList<Position>( pos.length );
331            // first point
332            pList.add( posList.get( 0 )[0] );
333            for ( int j = 0; j < posList.size() - 1; j++ ) {
334                // find intersection point between j'th and j+1'th segment
335                Position is = intersection( posList.get( j )[0], posList.get( j )[1], posList.get( j + 1 )[0],
336                                            posList.get( j + 1 )[1] );
337                if ( is != null ) {
338                    // is == null if to segments are parallel or part of the same line
339                    pList.add( is );
340                }
341            }
342            // last point
343            pList.add( posList.get( posList.size() - 1 )[1] );
344            curve = GeometryFactory.createCurve( pList.toArray( new Position[pList.size()] ), curve.getCoordinateSystem() );
345            return curve;
346        }
347    
348        /**
349         * 
350         * @param startPoint1
351         * @param endPoint1
352         * @param startPoint2
353         * @param endPoint2
354         * @return intersection coordinates between to lines (not line segments!!!). This means the intersection point may
355         *         not lies between passed start- and end-points
356         */
357        public static Position intersection( Position startPoint1, Position endPoint1, Position startPoint2,
358                                             Position endPoint2 ) {
359            double m1 = ( endPoint1.getY() - startPoint1.getY() ) / ( endPoint1.getX() - startPoint1.getX() );
360            double m2 = ( endPoint2.getY() - startPoint2.getY() ) / ( endPoint2.getX() - startPoint2.getX() );
361    
362            // lines are parallels
363            if ( m1 == m2 ) {
364                return null;
365            }
366    
367            double t1 = -( m1 * startPoint1.getX() - startPoint1.getY() );
368            double t2 = -( m2 * startPoint2.getX() - startPoint2.getY() );
369            double x = ( t2 - t1 ) / ( m1 - m2 );
370            double y = m1 * x + t1;
371            return GeometryFactory.createPosition( x, y );
372        }
373    
374        /**
375         * 
376         * @param a1
377         *            the first point
378         * @param a2
379         *            the second point
380         * @param l
381         *            the length of the segment starting from the second point
382         * @param alpha
383         *            the angle that the a1-a2 segment makes the following segment
384         * @return the point found at length l from a2 and which (connected with a2) forms an angle of alpha to a1a2.
385         */
386        public static Point vectorByAngle( Point a1, Point a2, double l, double alpha, boolean useAbsoluteAngle ) {
387            double beta = Math.atan2( a1.getY() - a2.getY(), a1.getX() - a2.getX() );
388    
389            double absoluteAngle;
390            if ( useAbsoluteAngle ) {
391                if ( alpha >= 0 ) {
392                    absoluteAngle = beta + alpha;
393                } else {
394                    absoluteAngle = beta - alpha;
395                }
396            } else {
397                absoluteAngle = beta - alpha;
398            }
399            return GeometryFactory.createPoint( a2.getX() + l * Math.cos( absoluteAngle ), a2.getY() + l
400                                                                                           * Math.sin( absoluteAngle ),
401                                                a1.getCoordinateSystem() );
402        }
403    
404        /**
405         * 
406         * @param center
407         * @param r
408         * @param noOfPos
409         * @param startPosition
410         * @param endPosition
411         * @return
412         * @throws GeometryException
413         */
414        public static Curve calcCircleCoordinates( Position center, double r, int nSeg, Position startPosition,
415                                                   Position endPosition, CoordinateSystem crs )
416                                throws GeometryException {
417            Curve curve = null;
418            if ( startPosition != null && endPosition != null ) {
419                double startCos = ( startPosition.getX() - center.getX() ) / r;
420                double startSin = ( startPosition.getY() - center.getY() ) / r;
421                double startArc = getArc( startSin, startCos );
422                double endCos = ( endPosition.getX() - center.getX() ) / r;
423                double endSin = ( endPosition.getY() - center.getY() ) / r;
424                double endArc = getArc( endSin, endCos );
425    
426                if ( startArc > endArc ) {
427                    double t = startArc;
428                    startArc = endArc;
429                    endArc = t;
430                }
431                if ( endArc - startArc > 180 ) {
432                    double t = startArc;
433                    startArc = endArc;
434                    endArc = 360 + t;
435                }
436    
437                curve = GeometryFactory.createCurveAsArc( center.getX(), center.getY(), r, r, nSeg, startArc, endArc, crs );
438                List<Position> l = null;
439                // ensure that curve starts/ends at start/end point of lines
440                if ( !curve.getEndPoint().getPosition().equals( startPosition ) ) {
441                    Position[] p = curve.getAsLineString().getPositions();
442                    l = new ArrayList<Position>( p.length + 1 );
443                    for ( Position position : p ) {
444                        l.add( position );
445                    }
446                    if ( GeometryUtils.distance( p[p.length - 1], startPosition ) < GeometryUtils.distance(
447                                                                                                            p[p.length - 1],
448                                                                                                            endPosition ) ) {
449                        l.add( startPosition );
450                    } else {
451                        l.add( endPosition );
452                    }
453                    curve = GeometryFactory.createCurve( l.toArray( new Position[l.size()] ), curve.getCoordinateSystem() );
454                }
455                if ( !curve.getEndPoint().getPosition().equals( endPosition ) ) {
456                    Position[] p = curve.getAsLineString().getPositions();
457                    l = new ArrayList<Position>( p.length + 1 );
458                    for ( Position position : p ) {
459                        l.add( position );
460                    }
461                    if ( GeometryUtils.distance( p[p.length - 1], startPosition ) < GeometryUtils.distance(
462                                                                                                            p[p.length - 1],
463                                                                                                            endPosition ) ) {
464                        l.add( startPosition );
465                    } else {
466                        l.add( endPosition );
467                    }
468                    curve = GeometryFactory.createCurve( l.toArray( new Position[l.size()] ), curve.getCoordinateSystem() );
469                }
470                if ( !curve.getStartPoint().getPosition().equals( startPosition )
471                     && curve.getStartPoint().getPosition().equals( endPosition ) ) {
472                    Position[] p = curve.getAsLineString().getPositions();
473                    l = new ArrayList<Position>( p.length + 1 );
474                    for ( Position position : p ) {
475                        l.add( position );
476                    }
477                    if ( GeometryUtils.distance( p[0], startPosition ) < GeometryUtils.distance( p[0], endPosition ) ) {
478                        l.add( 0, startPosition );
479                    } else {
480                        l.add( 0, endPosition );
481                    }
482                    curve = GeometryFactory.createCurve( l.toArray( new Position[l.size()] ), curve.getCoordinateSystem() );
483                }
484    
485            }
486            return curve;
487        }
488    
489        private static double getArc( double sin, double cos ) {
490            double atan = Math.atan2( sin, cos );
491            return Math.toDegrees( atan ) + 90;
492        }
493    
494        /**
495         * 
496         * @param p0x
497         * @param p0y
498         * @param p1x
499         * @param p1y
500         * @param p2x
501         * @param p2y
502         * @return arc between two line segments where p0x/p0y is common point of both segments
503         */
504        public static double getArc( double p0x, double p0y, double p1x, double p1y, double p2x, double p2y ) {
505            double d1 = GeometryUtils.distance( p0x, p0y, p1x, p1y );
506            double d2 = GeometryUtils.distance( p2x, p2y, p1x, p1y );
507            double d3 = GeometryUtils.distance( p2x, p2y, p0x, p0y );
508            double rad = 180 / Math.PI;
509            double s = ( d1 + d2 + d3 ) / 2d;
510            return rad * 2 * Math.asin( Math.sqrt( ( s - d1 ) * ( s - d3 ) / d1 / d3 ) );
511        }
512    
513        /**
514         * 
515         * @param p0x
516         * @param p0y
517         * @param p1x
518         * @param p1y
519         * @param p2x
520         * @param p2y
521         * @return true if p2x/p2y is left of the line defined by p0x/p0y p1x/p1y
522         */
523        public static boolean isLeft( double p0x, double p0y, double p1x, double p1y, double p2x, double p2y ) {
524            double p = ( p2x - p0x ) * ( p0y - p1y ) + ( p2y - p0y ) * ( p1x - p0x );
525            return ( p > 0 );
526        }
527    
528    }