geosAlgorithm.h

00001 /**********************************************************************
00002  * $Id: geosAlgorithm.h,v 1.9.2.1 2005/05/23 16:39:05 strk Exp $
00003  *
00004  * GEOS - Geometry Engine Open Source
00005  * http://geos.refractions.net
00006  *
00007  * Copyright (C) 2001-2002 Vivid Solutions Inc.
00008  * Copyright (C) 2005 Refractions Research Inc.
00009  *
00010  * This is free software; you can redistribute and/or modify it under
00011  * the terms of the GNU Lesser General Public Licence as published
00012  * by the Free Software Foundation. 
00013  * See the COPYING file for more information.
00014  *
00015  **********************************************************************/
00016 
00017 #ifndef GEOS_ALGORITHM_H
00018 #define GEOS_ALGORITHM_H
00019 
00020 #include <memory>
00021 #include <geos/geom.h>
00022 #include <geos/util.h>
00023 #include <geos/platform.h>
00024 #include <geos/indexBintree.h>
00025 #include <geos/indexStrtree.h>
00026 #include <geos/indexStrtree.h>
00027 #include <geos/indexChain.h>
00028 
00029 namespace geos {
00030 
00031 class Coordinate;
00032 
00033 /*
00034  * \class NotRepresentableException geosAlgorithm.h geos/geosAlgorithm.h
00035  * \brief
00036  * Indicates that a HCoordinate has been computed which is
00037  * not representable on the Cartesian plane.
00038  *
00039  * @version 1.4
00040  * @see HCoordinate
00041  */
00042 class NotRepresentableException: public GEOSException {
00043 public:
00044         NotRepresentableException();
00045         NotRepresentableException(string msg);
00046         ~NotRepresentableException();
00047 };
00048 
00049 class PointInRing{
00050 public:
00051         virtual ~PointInRing(){};
00052         virtual bool isInside(const Coordinate& pt)=0;
00053 };
00054 
00055 class CGAlgorithms {
00056 public:
00057         enum {
00058                 CLOCKWISE=-1,
00059                 COLLINEAR,
00060                 COUNTERCLOCKWISE
00061         };
00062         enum {
00063                 RIGHT=-1,
00064                 LEFT,
00065                 STRAIGHT
00066         };
00067         CGAlgorithms(){};
00068 
00082         static bool isPointInRing(const Coordinate& p, const CoordinateSequence* ring);
00090         static bool isOnLine(const Coordinate& p, const CoordinateSequence* pt);
00091 
00092         /*
00093          * Computes whether a ring defined by an array of Coordinate is
00094          * oriented counter-clockwise.
00095          * 
00096          *  - The list of points is assumed to have the first and last
00097          *    points equal.
00098          *  - This will handle coordinate lists which contain repeated points.
00099          *  - If the ring is invalid, the answer returned may not be correct.
00100          * 
00101          *
00102          * @param ring an array of coordinates forming a ring
00103          * @return <code>true</code> if the ring is oriented counter-clockwise.
00104          */
00105         static bool isCCW(const CoordinateSequence* ring);
00106 
00116         static int computeOrientation(const Coordinate& p1, const Coordinate& p2, const Coordinate& q);
00117         static double distancePointLine(const Coordinate& p,const Coordinate& A,const Coordinate& B);
00127         static double distancePointLinePerpendicular(const Coordinate& p,const Coordinate& A,const Coordinate& B);
00128         static double distanceLineLine(const Coordinate& A, const Coordinate& B, const Coordinate& C, const Coordinate& D);
00129         static double signedArea(const CoordinateSequence* ring);
00136         static double length(const CoordinateSequence* pts);
00150         static int orientationIndex(const Coordinate& p1,const Coordinate& p2,const Coordinate& q);
00151 
00152 };
00153 
00155 class HCoordinate {
00156 public:
00157         static Coordinate* intersection(Coordinate& p1,Coordinate& p2,Coordinate& q1,Coordinate& q2);
00158         double x,y,w;
00159         HCoordinate();
00160         HCoordinate(double _x, double _y, double _w);
00161         HCoordinate(Coordinate& p);
00162         HCoordinate(HCoordinate p1, HCoordinate p2);
00163         double getX();
00164         double getY();
00165         Coordinate* getCoordinate();
00166 };
00167 
00168 class SimplePointInRing: public PointInRing {
00169 public:
00170         SimplePointInRing(LinearRing *ring);
00171         virtual ~SimplePointInRing();
00172         bool isInside(const Coordinate& pt);
00173 private:
00174         const CoordinateSequence* pts;
00175 };
00176 
00177 class LineIntersector{
00178 public: 
00179         // Return a Z value being the interpolation of Z from p0 and p1 at
00180         // the given point p
00181         static double interpolateZ(const Coordinate &p, const Coordinate &p0, const Coordinate &p1);
00182         static double computeEdgeDistance(const Coordinate& p, const Coordinate& p0, const Coordinate& p1);
00183         static double nonRobustComputeEdgeDistance(const Coordinate& p,const Coordinate& p1,const Coordinate& p2);
00184         LineIntersector();
00185         virtual ~LineIntersector();
00186 
00187         /*
00188          * Tests whether either intersection point is an interior point of
00189          * one of the input segments.
00190          *
00191          * @return <code>true</code> if either intersection point is in
00192          * the interior of one of the input segments
00193          */
00194         virtual bool isInteriorIntersection();
00195 
00196         /*
00197          * Tests whether either intersection point is an interior point
00198          * of the specified input segment.
00199          *
00200          * @return <code>true</code> if either intersection point is in
00201          * the interior of the input segment
00202          */
00203         virtual bool isInteriorIntersection(int inputLineIndex);
00204 
00205         virtual void setMakePrecise(const PrecisionModel *newPM);
00206 
00207         virtual void setPrecisionModel(const PrecisionModel *newPM);
00208 
00209         /*
00210          * Compute the intersection of a point p and the line p1-p2.
00211          * This function computes the boolean value of the hasIntersection test.
00212          * The actual value of the intersection (if there is one)
00213          * is equal to the value of <code>p</code>.
00214          */
00215         virtual void computeIntersection(const Coordinate& p,const Coordinate& p1,const Coordinate& p2)=0;
00216 
00217         enum {
00218                 DONT_INTERSECT,
00219                 DO_INTERSECT,
00220                 COLLINEAR
00221         };
00222 
00223         virtual void computeIntersection(const Coordinate& p1,const Coordinate& p2,const Coordinate& p3, const Coordinate& p4);
00224 
00225         virtual string toString() const;
00226 
00227         virtual bool hasIntersection() const;
00228 
00229         virtual int getIntersectionNum() const;
00230 
00231         virtual const Coordinate& getIntersection(int intIndex) const;
00232 
00233         static bool isSameSignAndNonZero(double a,double b);
00234 
00235         virtual bool isIntersection(const Coordinate& pt) const;
00236         virtual bool isProper() const;
00237         virtual const Coordinate& getIntersectionAlongSegment(int segmentIndex,int intIndex);
00238         virtual int getIndexAlongSegment(int segmentIndex,int intIndex);
00239         virtual double getEdgeDistance(int geomIndex,int intIndex) const;
00240 protected:
00245         const PrecisionModel *precisionModel;
00246         int result;
00247         Coordinate inputLines[2][2];
00248         Coordinate intPt[2];
00253         int intLineIndex[2][2];
00254         bool isProperVar;
00255         Coordinate pa;
00256         Coordinate pb;
00257         virtual bool isCollinear() const;
00258         virtual int computeIntersect(const Coordinate& p1,const Coordinate& p2,const Coordinate& q1,const Coordinate& q2)=0;
00259         virtual bool isEndPoint() const;
00260         virtual void computeIntLineIndex();
00261         virtual void computeIntLineIndex(int segmentIndex);
00262 };
00263 
00264 class RobustDeterminant {
00265 public:
00266         static int signOfDet2x2(double x1,double y1,double x2,double y2);
00267 };
00268 
00269 class RobustLineIntersector: public LineIntersector {
00270 public:
00271         RobustLineIntersector();
00272         virtual ~RobustLineIntersector();
00273         void computeIntersection(const Coordinate& p,const Coordinate& p1,const Coordinate& p2);
00274         int computeIntersect(const Coordinate& p1,const Coordinate& p2,const Coordinate& q1,const Coordinate& q2);
00275 private:
00276 //      bool between(Coordinate& p1,Coordinate& p2,Coordinate& q);
00277         int computeCollinearIntersection(const Coordinate& p1,const Coordinate& p2,const Coordinate& q1,const Coordinate& q2);
00278         Coordinate* intersection(const Coordinate& p1,const Coordinate& p2,const Coordinate& q1,const Coordinate& q2) const;
00279         void normalize(Coordinate *n1,Coordinate *n2,Coordinate *n3,Coordinate *n4,Coordinate *normPt) const;
00280         double smallestInAbsValue(double x1,double x2,double x3,double x4) const;
00290         bool isInSegmentEnvelopes(const Coordinate& intPt);
00291 };
00292 
00293 class NonRobustLineIntersector: public LineIntersector {
00294 public:
00295         NonRobustLineIntersector();
00296         void computeIntersection(const Coordinate& p,const Coordinate& p1,const Coordinate& p2);
00297 protected:
00298         int computeIntersect(const Coordinate& p1,const Coordinate& p2,const Coordinate& p3,const Coordinate& p4);
00299 private:
00300         int computeCollinearIntersection(const Coordinate& p1,const Coordinate& p2,const Coordinate& p3,const Coordinate& p4);
00301         double rParameter(const Coordinate& p1,const Coordinate& p2,const Coordinate& p) const;
00302 };
00303 
00304 /*
00305  * Stub version of RobustCGAlgorithms for backwards compatibility.
00306  * Will be deprecated in next release - use CGAlgorithms instead.
00307  */
00308 class RobustCGAlgorithms: public CGAlgorithms {
00309 };
00310 
00311 class NonRobustCGAlgorithms: public CGAlgorithms {
00312 public:
00313         NonRobustCGAlgorithms();
00314         ~NonRobustCGAlgorithms();
00324         static bool isPointInRing(const Coordinate& p, const CoordinateSequence* ring);
00325 //      static bool isOnLine(const Coordinate& p, const CoordinateSequence* pt) const;
00336         static bool isCCW(const CoordinateSequence* ring);
00337         static int computeOrientation(const Coordinate& p1,const Coordinate& p2,const Coordinate& q);
00338 };
00339 
00340 class SimplePointInAreaLocator {
00341 public:
00342         static int locate(const Coordinate& p, const Geometry *geom);
00343         static bool containsPointInPolygon(const Coordinate& p,const Polygon *poly);
00344 private:
00345         static bool containsPoint(const Coordinate& p,const Geometry *geom);
00346 };
00347 
00348 /*
00349  * \class PointLocator geosAlgorithm.h geos/geosAlgorithm.h
00350  *
00351  * \brief
00352  * Computes the topological relationship (Location)
00353  * of a single point to a Geometry.
00354  *
00355  * The algorithm obeys the SFS boundaryDetermination rule to correctly determine
00356  * whether the point lies on the boundary or not.
00357  * Note that instances of this class are not reentrant.
00358  * @version 1.3
00359  */
00360 class PointLocator {
00361 public:
00362         PointLocator();
00363         ~PointLocator();
00364         int locate(const Coordinate& p,const Geometry *geom);
00365         bool intersects(const Coordinate& p,const Geometry *geom);
00366         int locate(const Coordinate& p,const LineString *l);
00367         int locate(const Coordinate& p,const LinearRing *ring);
00368         int locate(const Coordinate& p,const Polygon *poly);
00369 private:
00370         bool isIn;         // true if the point lies in or on any Geometry element
00371         int numBoundaries;    // the number of sub-elements whose boundaries the point lies in
00372         void computeLocation(const Coordinate& p,const Geometry *geom);
00373         void updateLocationInfo(int loc);
00374 };
00375 
00376 
00377 class MCPointInRing: public PointInRing {
00378 public:
00379         MCPointInRing(LinearRing *newRing);
00380         virtual ~MCPointInRing();
00381         bool isInside(const Coordinate& pt);
00382         void testLineSegment(Coordinate& p,LineSegment *seg);
00383         class MCSelecter: public MonotoneChainSelectAction {
00384         private:
00385                 Coordinate p;
00386                 MCPointInRing *parent;
00387         public:
00388                 MCSelecter(const Coordinate& newP,MCPointInRing *prt);
00389             void select(LineSegment *ls);
00390         };
00391 private:
00392         LinearRing *ring;
00393         BinTreeInterval *interval;
00394         CoordinateSequence *pts;
00395         Bintree *tree;
00396         int crossings;  // number of segment/ray crossings
00397         void buildIndex();
00398         void testMonotoneChain(Envelope *rayEnv,MCSelecter *mcSelecter,indexMonotoneChain *mc);
00399 };
00400 
00401 class SIRtreePointInRing: public PointInRing {
00402 private:
00403         LinearRing *ring;
00404         SIRtree *sirTree;
00405         int crossings;  // number of segment/ray crossings
00406         void buildIndex();
00407         void testLineSegment(const Coordinate& p,LineSegment *seg);
00408 public:
00409         SIRtreePointInRing(LinearRing *newRing);
00410         bool isInside(const Coordinate& pt);
00411 };
00412 
00413 class CentroidPoint {
00414 private:
00415         int ptCount;
00416         Coordinate* centSum;
00417 public:
00418         CentroidPoint();
00419         virtual ~CentroidPoint();
00420         void add(const Geometry *geom);
00421         void add(const Coordinate *pt);
00422         Coordinate* getCentroid() const;
00423 };
00424 
00425 class CentroidLine {
00426 private:
00427         Coordinate* centSum;
00428         double totalLength;
00429 public:
00430         CentroidLine();
00431         virtual ~CentroidLine();
00432         void add(const Geometry *geom);
00433         void add(const CoordinateSequence *pts);
00434         Coordinate* getCentroid() const;
00435 };
00436 
00437 /*
00438  * \class CentroidArea geosAlgorithm.h geos/geosAlgorithm.h
00439  *
00440  * \brief Computes the centroid of an area geometry.
00441  *
00442  * Algorithm:
00443  *
00444  * Based on the usual algorithm for calculating
00445  * the centroid as a weighted sum of the centroids
00446  * of a decomposition of the area into (possibly overlapping) triangles.
00447  * The algorithm has been extended to handle holes and multi-polygons.
00448  * See <code>http://www.faqs.org/faqs/graphics/algorithms-faq/</code>
00449  * for further details of the basic approach.
00450  */
00451 class CentroidArea {
00452 public:
00453         CentroidArea();
00454         virtual ~CentroidArea();
00455         void add(const Geometry *geom);
00456         void add(const CoordinateSequence *ring);
00457         Coordinate* getCentroid() const;
00458 private:
00459         CGAlgorithms *cga;
00460         Coordinate* basePt;// the point all triangles are based at
00461         Coordinate* triangleCent3;// temporary variable to hold centroid of triangle
00462         double areasum2;        /* Partial area sum */
00463         Coordinate* cg3; // partial centroid sum
00464         void setBasePoint(const Coordinate *newbasePt);
00465         void add(const Polygon *poly);
00466         void addShell(const CoordinateSequence *pts);
00467         void addHole(const CoordinateSequence *pts);
00468         inline void addTriangle(const Coordinate &p0, const Coordinate &p1, const Coordinate &p2,bool isPositiveArea);
00469         static inline  void centroid3(const Coordinate &p1,const Coordinate &p2,const Coordinate &p3,Coordinate *c);
00470         static inline double area2(const Coordinate &p1,const Coordinate &p2,const Coordinate &p3);
00471 };
00472 
00473 /*
00474  * \class InteriorPointPoint geosAlgorithm.h geos/geosAlgorithm.h
00475  * \brief
00476  * Computes a point in the interior of an point geometry.
00477  *
00478  * Algorithm:
00479  *
00480  * Find a point which is closest to the centroid of the geometry.
00481  */
00482 class InteriorPointPoint {
00483 private:
00484         const Coordinate* centroid;
00485         double minDistance;
00486         Coordinate *interiorPoint;
00487         void add(const Geometry *geom);
00488         void add(const Coordinate *point);
00489 public:
00490         InteriorPointPoint(const Geometry *g);
00491         virtual ~InteriorPointPoint();
00492         Coordinate* getInteriorPoint() const;
00493 };
00494 
00495 /*
00496  * Computes a point in the interior of an linear geometry.
00497  * <h2>Algorithm</h2>
00498  * <ul>
00499  * <li>Find an interior vertex which is closest to
00500  * the centroid of the linestring.
00501  * <li>If there is no interior vertex, find the endpoint which is
00502  * closest to the centroid.
00503  * </ul>
00504  */
00505 class InteriorPointLine {
00506 public:
00507         InteriorPointLine(Geometry *g);
00508         virtual ~InteriorPointLine();
00509         Coordinate* getInteriorPoint() const;
00510 private:
00511         const Coordinate *centroid;
00512         double minDistance;
00513         Coordinate *interiorPoint;
00514         void addInterior(const Geometry *geom);
00515         void addInterior(const CoordinateSequence *pts);
00516         void addEndpoints(const Geometry *geom);
00517         void addEndpoints(const CoordinateSequence *pts);
00518         void add(const Coordinate *point);
00519 };
00520 
00521 /*
00522  * Computes a point in the interior of an area geometry.
00523  *
00524  * <h2>Algorithm</h2>
00525  * <ul>
00526  *   <li>Find the intersections between the geometry
00527  *       and the horizontal bisector of the area's envelope
00528  *   <li>Pick the midpoint of the largest intersection (the intersections
00529  *       will be lines and points)
00530  * </ul>
00531  *
00532  * <b>
00533  * Note: If a fixed precision model is used,
00534  * in some cases this method may return a point
00535  * which does not lie in the interior.
00536  * </b>
00537  */
00538 class InteriorPointArea {
00539 private:
00540         static double avg(double a, double b);
00541         const GeometryFactory *factory;
00542         Coordinate *interiorPoint;
00543         double maxWidth;
00544         void add(const Geometry *geom);
00545 public:
00546         InteriorPointArea(const Geometry *g);
00547         virtual ~InteriorPointArea();
00548         Coordinate* getInteriorPoint() const;
00549         void addPolygon(const Geometry *geometry);
00550         Coordinate* centre(const Envelope *envelope) const;
00551 protected:
00552         const Geometry *widestGeometry(const Geometry *geometry);
00553         const Geometry *widestGeometry(const GeometryCollection *gc);
00554         LineString *horizontalBisector(const Geometry *geometry);
00555 
00556 };
00557 
00558 class BigQuad {
00559 public:
00560         Coordinate northmost;
00561         Coordinate southmost;
00562         Coordinate westmost;
00563         Coordinate eastmost;
00564 };
00565 
00566 class ConvexHull {
00567 private:
00568         PointLocator *pointLocator;
00569         //CGAlgorithms *cgAlgorithms;
00570         const Geometry *geometry;
00571         const GeometryFactory *factory;
00572         CoordinateSequence* reduce(const CoordinateSequence *pts);
00573         CoordinateSequence* preSort(CoordinateSequence *pts);
00574         CoordinateSequence* grahamScan(const CoordinateSequence *c);
00575         void radialSort(CoordinateSequence *p);
00576         int polarCompare(Coordinate o, Coordinate p, Coordinate q);
00577         bool isBetween(Coordinate c1, Coordinate c2, Coordinate c3);
00578     BigQuad* makeBigQuad(const CoordinateSequence *pts);
00579         Geometry* lineOrPolygon(CoordinateSequence *newCoordinates);
00580         CoordinateSequence* cleanRing(CoordinateSequence *original);
00581 public:
00582         ConvexHull(const Geometry *newGeometry);
00583         ~ConvexHull();
00584         Geometry* getConvexHull();
00585 };
00586 
00587 
00588 /*
00589  * Computes the minimum diameter of a {@link Geometry}.
00590  * The minimum diameter is defined to be the
00591  * width of the smallest band that
00592  * contains the geometry,
00593  * where a band is a strip of the plane defined
00594  * by two parallel lines.
00595  * This can be thought of as the smallest hole that the geometry can be
00596  * moved through, with a single rotation.
00597  * <p>
00598  * The first step in the algorithm is computing the convex hull of the Geometry.
00599  * If the input Geometry is known to be convex, a hint can be supplied to
00600  * avoid this computation.
00601  *
00602  * @see ConvexHull
00603  *
00604  */
00605 class MinimumDiameter {
00606 private:
00607         const Geometry* inputGeom;
00608         bool isConvex;
00609         LineSegment* minBaseSeg;
00610         Coordinate* minWidthPt;
00611         int minPtIndex;
00612         double minWidth;
00613         void computeMinimumDiameter();
00614         void computeWidthConvex(const Geometry* geom);
00622         void computeConvexRingMinDiameter(const CoordinateSequence *pts);
00623         int findMaxPerpDistance(const CoordinateSequence* pts, LineSegment* seg, int startIndex);
00624         static int getNextIndex(const CoordinateSequence* pts, int index);
00625 public:
00626         ~MinimumDiameter();
00627 
00633         MinimumDiameter(const Geometry* newInputGeom);
00634 
00645         MinimumDiameter(const Geometry* newInputGeom,const bool newIsConvex);
00646 
00652         double getLength();
00653 
00659         Coordinate* getWidthCoordinate();
00660 
00666         LineString* getSupportingSegment();
00667 
00673         LineString* getDiameter();
00674 };
00675 
00676 } // namespace geos
00677 
00678 #endif
00679 
00680 /**********************************************************************
00681  * $Log: geosAlgorithm.h,v $
00682  * Revision 1.9.2.1  2005/05/23 16:39:05  strk
00683  * log lines moved at bottom file, added Refractions copyright
00684  *
00685  * Revision 1.9  2004/11/23 19:53:07  strk
00686  * Had LineIntersector compute Z by interpolation.
00687  *
00688  * Revision 1.8  2004/11/06 08:16:46  strk
00689  * Fixed CGAlgorithms::isCCW from JTS port.
00690  * Code cleanup in IsValidOp.
00691  *
00692  * Revision 1.7  2004/10/21 22:29:54  strk
00693  * Indentation changes and some more COMPUTE_Z rules
00694  *
00695  * Revision 1.6  2004/09/13 10:12:49  strk
00696  * Added invalid coordinates checks in IsValidOp.
00697  * Cleanups.
00698  *
00699  **********************************************************************/
00700 

Generated on Wed Jan 10 22:13:55 2007 for GEOS by  doxygen 1.5.1