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 }