src/algorithm/Centroid.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  *
   8  * This is free software; you can redistribute and/or modify it under
   9  * the terms of the GNU Lesser General Public Licence as published
  10  * by the Free Software Foundation.
  11  * See the COPYING file for more information.
  12  *
  13  **********************************************************************
  14  *
  15  * Last port: algorithm/Centroid.java r728 (JTS-0.13+)
  16  *
  17  **********************************************************************/
  18
  19 #include <geos/algorithm/Centroid.h>
  20 #include <geos/algorithm/CGAlgorithms.h>
  21 #include <geos/geom/CoordinateSequence.h>
  22 #include <geos/geom/Geometry.h>
  23 #include <geos/geom/Point.h>
  24 #include <geos/geom/Polygon.h>
  25 #include <geos/geom/GeometryCollection.h>
  26 #include <geos/geom/LineString.h>
  27
  28
  29 #include <cmath> // for std::abs
  30
  31 using namespace geos::geom;
  32
  33 namespace geos {
  34 namespace algorithm { // geos.algorithm
  35
  36 /* static public */
  37 bool
  38 Centroid::getCentroid(const Geometry& geom, Coordinate& pt)
  39 {
  40   Centroid cent(geom);
  41   return cent.getCentroid(pt);
  42 }
  43
  44 /* public */
  45 bool
  46 Centroid::getCentroid(Coordinate& cent) const
  47 {
  48   if (std::abs(areasum2) > 0.0) {
  49     cent.x = cg3.x / 3 / areasum2;
  50     cent.y = cg3.y / 3 / areasum2;
  51   }
  52   else if (totalLength > 0.0) {
  53     // if polygon was degenerate, compute linear centroid instead
  54     cent.x = lineCentSum.x / totalLength;
  55     cent.y = lineCentSum.y / totalLength;
  56   }
  57   else if (ptCount > 0){
  58     cent.x = ptCentSum.x / ptCount;
  59     cent.y = ptCentSum.y / ptCount;
  60   }
  61   else {
  62     return false;
  63   }
  64   return true;
  65 }
  66
  67 /* private */
  68 void
  69 Centroid::add(const Geometry& geom)
  70 {
  71   if (geom.isEmpty()) return;
  72
  73   if (const Point* g = dynamic_cast<const Point*>(&geom)) {
  74     addPoint(*g->getCoordinate());
  75   }
  76   else if (const LineString* g = dynamic_cast<const LineString*>(&geom)) {
  77     addLineSegments(*g->getCoordinatesRO());
  78   }
  79   else if (const Polygon* g = dynamic_cast<const Polygon*>(&geom)) {
  80     add(*g);
  81   }
  82   else if (const GeometryCollection* g = dynamic_cast<const GeometryCollection*>(&geom)) {
  83     for (size_t i = 0; i < g->getNumGeometries(); i++) {
  84       add(*g->getGeometryN(i));
  85     }
  86   }
  87 }
  88
  89 /* private */
  90 void
  91 Centroid::setBasePoint(const Coordinate& basePt)
  92 {
  93   if ( ! areaBasePt.get() )
  94       areaBasePt.reset( new Coordinate(basePt) );
  95 }
  96
  97 /* private */
  98 void
  99 Centroid::add(const Polygon& poly)
 100 {
 101   addShell(*poly.getExteriorRing()->getCoordinatesRO());
 102   for (size_t i = 0; i < poly.getNumInteriorRing(); i++) {
 103     addHole(*poly.getInteriorRingN(i)->getCoordinatesRO());
 104   }
 105 }
 106
 107 /* private */
 108 void
 109 Centroid::addShell(const CoordinateSequence& pts)
 110 {
 111   size_t len = pts.size();
 112   if (len > 0)
 113     setBasePoint(pts[0]);
 114   bool isPositiveArea = ! CGAlgorithms::isCCW(&pts);
 115   for (size_t i = 0; i < len - 1; ++i) {
 116     addTriangle(*areaBasePt, pts[i], pts[i+1], isPositiveArea);
 117   }
 118   addLineSegments(pts);
 119 }
 120
 121 /* private */
 122 void
 123 Centroid::addHole(const CoordinateSequence& pts)
 124 {
 125   bool isPositiveArea = CGAlgorithms::isCCW(&pts);
 126   for (size_t i = 0, e = pts.size() - 1; i < e; ++i) {
 127     addTriangle(*areaBasePt, pts[i], pts[i+1], isPositiveArea);
 128   }
 129   addLineSegments(pts);
 130 }
 131
 132 /* private */
 133 void
 134 Centroid::addTriangle(const Coordinate& p0, const Coordinate& p1, const Coordinate& p2, bool isPositiveArea)
 135 {
 136   double sign = (isPositiveArea) ? 1.0 : -1.0;
 137   centroid3( p0, p1, p2, triangleCent3 );
 138   double a2 = area2( p0, p1, p2 );
 139   cg3.x += sign * a2 * triangleCent3.x;
 140   cg3.y += sign * a2 * triangleCent3.y;
 141   areasum2 += sign * a2;
 142 }
 143
 144 /* static private */
 145 void
 146 Centroid::centroid3(const Coordinate& p1, const Coordinate& p2, const Coordinate& p3, Coordinate& c )
 147 {
 148   c.x = p1.x + p2.x + p3.x;
 149   c.y = p1.y + p2.y + p3.y;
 150   return;
 151 }
 152
 153 /* static private */
 154 double
 155 Centroid::area2(const Coordinate& p1, const Coordinate& p2, const Coordinate& p3 )
 156 {
 157   return
 158     (p2.x - p1.x) * (p3.y - p1.y) -
 159     (p3.x - p1.x) * (p2.y - p1.y);
 160 }
 161
 162 /* private */
 163 void
 164 Centroid::addLineSegments(const CoordinateSequence& pts)
 165 {
 166   size_t npts = pts.size();
 167   double lineLen = 0.0;
 168   for (size_t i = 0; i < npts - 1; i++) {
 169     double segmentLen = pts[i].distance(pts[i + 1]);
 170     if (segmentLen == 0.0)
 171       continue;
 172
 173     lineLen += segmentLen;
 174
 175     double midx = (pts[i].x + pts[i + 1].x) / 2;
 176     lineCentSum.x += segmentLen * midx;
 177     double midy = (pts[i].y + pts[i + 1].y) / 2;
 178     lineCentSum.y += segmentLen * midy;
 179   }
 180   totalLength += lineLen;
 181   if (lineLen == 0.0 && npts > 0)
 182     addPoint(pts[0]);
 183 }
 184
 185 /* private */
 186 void
 187 Centroid::addPoint(const Coordinate& pt)
 188 {
 189   ptCount += 1;
 190   ptCentSum.x += pt.x;
 191   ptCentSum.y += pt.y;
 192 }
 193
 194 } // namespace geos.algorithm
 195 } // namespace geos