src/algorithm/InteriorPointArea.cpp

   1 /**********************************************************************
   2  *
   3  * GEOS - Geometry Engine Open Source
   4  * http://geos.osgeo.org
   5  *
   6  * Copyright (C) 2013 Sandro Santilli <strk@keybit.net>
   7  * Copyright (C) 2005-2006 Refractions Research Inc.
   8  * Copyright (C) 2001-2002 Vivid Solutions Inc.
   9  *
  10  * This is free software; you can redistribute and/or modify it under
  11  * the terms of the GNU Lesser General Public Licence as published
  12  * by the Free Software Foundation.
  13  * See the COPYING file for more information.
  14  *
  15  **********************************************************************
  16  *
  17  * Last port: algorithm/InteriorPointArea.java r728 (JTS-1.13+)
  18  *
  19  **********************************************************************/
  20
  21 #include <geos/algorithm/InteriorPointArea.h>
  22 #include <geos/geom/Coordinate.h>
  23 #include <geos/geom/Geometry.h>
  24 #include <geos/geom/GeometryCollection.h>
  25 #include <geos/geom/Polygon.h>
  26 #include <geos/geom/LineString.h>
  27 #include <geos/geom/Envelope.h>
  28 #include <geos/geom/GeometryFactory.h>
  29 #include <geos/geom/CoordinateSequenceFactory.h>
  30
  31 #include <vector>
  32 #include <typeinfo>
  33 #include <memory> // for auto_ptr
  34
  35 using namespace std;
  36 using namespace geos::geom;
  37
  38 namespace geos {
  39 namespace algorithm { // geos.algorithm
  40
  41 // file-statics
  42 namespace {
  43
  44   double avg(double a, double b){return (a+b)/2.0;}
  45
  46   /**
  47    * Finds a safe bisector Y ordinate
  48    * by projecting to the Y axis
  49    * and finding the Y-ordinate interval
  50    * which contains the centre of the Y extent.
  51    * The centre of this interval is returned as the bisector Y-ordinate.
  52    *
  53    * @author mdavis
  54    *
  55    */
  56   class SafeBisectorFinder
  57   {
  58   public:
  59           static double getBisectorY(const Polygon& poly)
  60           {
  61                   SafeBisectorFinder finder(poly);
  62                   return finder.getBisectorY();
  63           }
  64           SafeBisectorFinder(const Polygon& nPoly)
  65       : poly(nPoly)
  66     {
  67                   // initialize using extremal values
  68                   hiY = poly.getEnvelopeInternal()->getMaxY();
  69                   loY = poly.getEnvelopeInternal()->getMinY();
  70                   centreY = avg(loY, hiY);
  71           }
  72
  73           double getBisectorY()
  74           {
  75                   process(*poly.getExteriorRing());
  76                   for (size_t i = 0; i < poly.getNumInteriorRing(); i++) {
  77                           process(*poly.getInteriorRingN(i));
  78                   }
  79                   double bisectY = avg(hiY, loY);
  80                   return bisectY;
  81           }
  82
  83
  84         private:
  85           const Polygon& poly;
  86
  87           double centreY;
  88           double hiY;
  89           double loY;
  90
  91           void process(const LineString& line) {
  92       const CoordinateSequence* seq = line.getCoordinatesRO();
  93       for (std::size_t i = 0, s = seq->size(); i < s; i++) {
  94         double y = seq->getY(i);
  95         updateInterval(y);
  96       }
  97     }
  98
  99     void updateInterval(double y) {
 100       if (y <= centreY) {
 101         if (y > loY)
 102           loY = y;
 103       }
 104       else if (y > centreY) {
 105         if (y < hiY) {
 106           hiY = y;
 107         }
 108       }
 109     }
 110
 111     SafeBisectorFinder(SafeBisectorFinder const&); /*= delete*/
 112     SafeBisectorFinder& operator=(SafeBisectorFinder const&); /*= delete*/
 113   };
 114
 115 } // anonymous namespace
 116
 117
 118 /*public*/
 119 InteriorPointArea::InteriorPointArea(const Geometry *g)
 120 {
 121         foundInterior=false;
 122         maxWidth=0.0;
 123         factory=g->getFactory();
 124         add(g);
 125 }
 126
 127 /*public*/
 128 InteriorPointArea::~InteriorPointArea()
 129 {
 130 }
 131
 132 /*public*/
 133 bool
 134 InteriorPointArea::getInteriorPoint(Coordinate& ret) const
 135 {
 136         if ( ! foundInterior ) return false;
 137
 138         ret=interiorPoint;
 139         return true;
 140 }
 141
 142 /*public*/
 143 void
 144 InteriorPointArea::add(const Geometry *geom)
 145 {
 146         const Polygon *poly = dynamic_cast<const Polygon*>(geom);
 147         if ( poly ) {
 148                 addPolygon(geom);
 149                 return;
 150         }
 151
 152         const GeometryCollection *gc = dynamic_cast<const GeometryCollection*>(geom);
 153         if ( gc )
 154         {
 155         for(std::size_t i=0, n=gc->getNumGeometries(); i<n; i++) {
 156                         add(gc->getGeometryN(i));
 157                 }
 158         }
 159 }
 160
 161 /*private*/
 162 void
 163 InteriorPointArea::addPolygon(const Geometry *geometry)
 164 {
 165   if (geometry->isEmpty()) return;
 166
 167   Coordinate intPt;
 168   double width;
 169
 170   auto_ptr<LineString> bisector ( horizontalBisector(geometry) );
 171   if ( bisector->getLength() == 0.0 ) {
 172     width = 0;
 173     intPt = bisector->getCoordinateN(0);
 174   }
 175   else {
 176     auto_ptr<Geometry> intersections ( bisector->intersection(geometry) );
 177     const Geometry *widestIntersection = widestGeometry(intersections.get());
 178     const Envelope *env = widestIntersection->getEnvelopeInternal();
 179     width=env->getWidth();
 180     env->centre(intPt);
 181   }
 182   if (!foundInterior || width>maxWidth) {
 183     interiorPoint = intPt;
 184     maxWidth = width;
 185     foundInterior=true;
 186   }
 187 }
 188
 189 //@return if geometry is a collection, the widest sub-geometry; otherwise,
 190 //the geometry itself
 191 const Geometry*
 192 InteriorPointArea::widestGeometry(const Geometry *geometry)
 193 {
 194         const GeometryCollection *gc = dynamic_cast<const GeometryCollection*>(geometry);
 195         if ( gc ) {
 196                 return widestGeometry(gc);
 197         } else {
 198                 return geometry;
 199         }
 200 }
 201
 202 const Geometry*
 203 InteriorPointArea::widestGeometry(const GeometryCollection* gc) {
 204         if (gc->isEmpty()) {
 205                 return gc;
 206         }
 207         const Geometry* widestGeometry=gc->getGeometryN(0);
 208
 209         // scan remaining geom components to see if any are wider
 210         for(std::size_t i=1, n=gc->getNumGeometries(); i<n; i++) // start at 1
 211         {
 212                 const Envelope *env1(gc->getGeometryN(i)->getEnvelopeInternal());
 213                 const Envelope *env2(widestGeometry->getEnvelopeInternal());
 214                 if (env1->getWidth()>env2->getWidth()) {
 215                                 widestGeometry=gc->getGeometryN(i);
 216                 }
 217         }
 218         return widestGeometry;
 219 }
 220
 221 /* private */
 222 LineString*
 223 InteriorPointArea::horizontalBisector(const Geometry *geometry)
 224 {
 225         const Envelope *envelope=geometry->getEnvelopeInternal();
 226
 227         /**
 228          * Original algorithm.  Fails when geometry contains a horizontal
 229          * segment at the Y midpoint.
 230          */
 231         // Assert: for areas, minx <> maxx
 232         //double avgY=avg(envelope->getMinY(),envelope->getMaxY());
 233
 234         double bisectY = SafeBisectorFinder::getBisectorY(*dynamic_cast<const Polygon *>(geometry));
 235         vector<Coordinate>*cv=new vector<Coordinate>(2);
 236         (*cv)[0].x = envelope->getMinX();
 237         (*cv)[0].y = bisectY;
 238         (*cv)[1].x = envelope->getMaxX();
 239         (*cv)[1].y = bisectY;
 240
 241         CoordinateSequence *cl = factory->getCoordinateSequenceFactory()->create(cv);
 242
 243         LineString *ret = factory->createLineString(cl);
 244         return ret;
 245 }
 246
 247 } // namespace geos.algorithm
 248 } // namespace geos