Osmium  0.1
include/osmium/osm/area.hpp
Go to the documentation of this file.
00001 #ifndef OSMIUM_OSM_AREA_HPP
00002 #define OSMIUM_OSM_AREA_HPP
00003 
00004 /*
00005 
00006 Copyright 2011 Jochen Topf <jochen@topf.org> and others (see README).
00007 
00008 This file is part of Osmium (https://github.com/joto/osmium).
00009 
00010 Osmium is free software: you can redistribute it and/or modify it under the
00011 terms of the GNU Lesser General Public License or (at your option) the GNU
00012 General Public License as published by the Free Software Foundation, either
00013 version 3 of the Licenses, or (at your option) any later version.
00014 
00015 Osmium is distributed in the hope that it will be useful, but WITHOUT ANY
00016 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
00017 PARTICULAR PURPOSE. See the GNU Lesser General Public License and the GNU
00018 General Public License for more details.
00019 
00020 You should have received a copy of the Licenses along with Osmium. If not, see
00021 <http://www.gnu.org/licenses/>.
00022 
00023 */
00024 
00025 #include <sys/types.h>
00026 #include <string.h>
00027 #include <stdio.h>
00028 #include <stdlib.h>
00029 #include <stdarg.h>
00030 #include <vector>
00031 #include <sstream>
00032 #include <iomanip>
00033 #include <map>
00034 
00035 #ifdef OSMIUM_WITH_MULTIPOLYGON_PROFILING
00036 # include <osmium/utils/timer.h>
00037 # define START_TIMER(x) x##_timer.start();
00038 # define STOP_TIMER(x) x##_timer.stop();
00039 #else
00040 # define START_TIMER(x)
00041 # define STOP_TIMER(x)
00042 #endif // OSMIUM_WITH_MULTIPOLYGON_PROFILING
00043 
00044 #ifdef OSMIUM_WITH_GEOS
00045 # include <geos/version.h>
00046 # include <geos/geom/Geometry.h>
00047 # include <geos/geom/Point.h>
00048 # include <geos/geom/LineString.h>
00049 # include <geos/geom/Polygon.h>
00050 # include <geos/geom/PrecisionModel.h>
00051 # include <geos/geom/CoordinateSequence.h>
00052 # include <geos/geom/CoordinateArraySequenceFactory.h>
00053 # include <geos/geom/LinearRing.h>
00054 # include <geos/geom/Polygon.h>
00055 # include <geos/geom/MultiPolygon.h>
00056 # include <geos/geom/MultiLineString.h>
00057 # include <geos/geom/prep/PreparedPolygon.h>
00058 # include <geos/io/WKTWriter.h>
00059 # include <geos/util/GEOSException.h>
00060 # include <geos/opLinemerge.h>
00061 # include <geos/operation/polygonize/Polygonizer.h>
00062 # include <geos/operation/distance/DistanceOp.h>
00063 # include <geos/opPolygonize.h>
00064 # include <geos/algorithm/LineIntersector.h>
00065 # include <geos/geomgraph/GeometryGraph.h>
00066 # include <geos/geomgraph/index/SegmentIntersector.h>
00067 
00068 // this should come from /usr/include/geos/algorithm, but its missing there in some Ubuntu versions
00069 # include "../CGAlgorithms.h"
00070 #endif // OSMIUM_WITH_GEOS
00071 
00072 #ifdef OSMIUM_WITH_SHPLIB
00073 # include <shapefil.h>
00074 #endif // OSMIUM_WITH_SHPLIB
00075 
00076 #include <osmium/osm/way.hpp>
00077 #include <osmium/osm/relation.hpp>
00078 #include <osmium/geometry.hpp>
00079 
00080 namespace Osmium {
00081 
00082     namespace OSM {
00083 
00084         enum innerouter_t { UNSET, INNER, OUTER };
00085         enum direction_t { NO_DIRECTION, CLOCKWISE, COUNTERCLOCKWISE };
00086 
00087 #ifdef OSMIUM_WITH_GEOS
00088         class WayInfo {
00089 
00090             friend class AreaFromRelation;
00091 
00092             Osmium::OSM::Way *way;
00093             int used;
00094             int sequence;
00095             bool invert;
00096             bool duplicate;
00097             std::string errorhint;
00098             innerouter_t innerouter;
00099             innerouter_t orig_innerouter;
00100             geos::geom::Geometry *way_geom;
00101             int firstnode;
00102             int lastnode;
00103             bool tried;
00104 
00105             WayInfo() {
00106                 way = NULL;
00107                 used = -1;
00108                 innerouter = UNSET;
00109                 orig_innerouter = UNSET;
00110                 sequence = 0;
00111                 invert = false;
00112                 duplicate = false;
00113                 way_geom = NULL;
00114                 firstnode = -1;
00115                 lastnode = -1;
00116                 tried = false;
00117             }
00118 
00119             WayInfo(Osmium::OSM::Way *w, innerouter_t io) {
00120                 way = w;
00121                 way_geom = w->create_geos_geometry();
00122                 orig_innerouter = io;
00123                 used = -1;
00124                 innerouter = UNSET;
00125                 sequence = 0;
00126                 invert = false;
00127                 duplicate = false;
00128                 tried = false;
00129                 firstnode = w->get_first_node_id();
00130                 lastnode = w->get_last_node_id();
00131             }
00132 
00134             WayInfo(geos::geom::Geometry *g, int first, int last, innerouter_t io) {
00135                 way = NULL;
00136                 way_geom = g;
00137                 orig_innerouter = io;
00138                 used = -1;
00139                 innerouter = UNSET;
00140                 sequence = 0;
00141                 invert = false;
00142                 duplicate = false;
00143                 tried = false;
00144                 firstnode = first;
00145                 lastnode = last;
00146             }
00147 
00148             ~WayInfo() {
00149                 delete way_geom;
00150             }
00151 
00152             geos::geom::Point *get_firstnode_geom() {
00153                 return (way ? way->get_first_node_geometry() : NULL);
00154             }
00155 
00156             geos::geom::Point *get_lastnode_geom() {
00157                 return (way ? way->get_last_node_geometry() : NULL);
00158             }
00159 
00160         };
00161 
00162         class RingInfo {
00163 
00164             friend class AreaFromRelation;
00165 
00166             geos::geom::Polygon *polygon;
00167             direction_t direction;
00168             std::vector<WayInfo *> ways;
00169             std::vector<RingInfo *> inner_rings;
00170             bool nested;
00171             RingInfo *contained_by;
00172             int ring_id;
00173 
00174             RingInfo() {
00175                 polygon = NULL;
00176                 direction = NO_DIRECTION;
00177                 nested = false;
00178                 contained_by = NULL;
00179                 ring_id = -1;
00180             }
00181         };
00182 #endif // OSMIUM_WITH_GEOS
00183 
00185         class Area : public Object {
00186 
00187         protected:
00188 
00189 #ifdef OSMIUM_WITH_GEOS
00190             geos::geom::Geometry *geometry;
00191 
00192             Area(geos::geom::Geometry *geom = NULL) : geometry(geom) {
00193 #else
00194             Area() {
00195 #endif // OSMIUM_WITH_GEOS
00196 # ifdef OSMIUM_WITH_JAVASCRIPT
00197                 v8::HandleScope handle_scope;
00198                 js_object_instance = v8::Persistent<v8::Object>::New(JavascriptTemplate::get<JavascriptTemplate>().create_instance(this));
00199 # endif // OSMIUM_WITH_JAVASCRIPT
00200             }
00201 
00202             ~Area() {
00203 #ifdef OSMIUM_WITH_GEOS
00204                 delete(geometry);
00205 #endif // OSMIUM_WITH_GEOS
00206             }
00207 
00208         public:
00209 
00210 #ifdef OSMIUM_WITH_GEOS
00211             geos::geom::Geometry* get_geometry() const {
00212                 return geometry;
00213             }
00214 #endif // OSMIUM_WITH_GEOS
00215 
00216 #ifdef OSMIUM_WITH_JAVASCRIPT
00217             v8::Handle<v8::Value> js_from() const {
00218                 const char *value = (get_type() == AREA_FROM_WAY) ? "way" : "relation";
00219                 return v8::String::NewSymbol(value);
00220             }
00221 
00222             v8::Handle<v8::Value> js_geom() const;
00223 
00224             struct JavascriptTemplate : public Osmium::OSM::Object::JavascriptTemplate {
00225 
00226                 JavascriptTemplate() : Osmium::OSM::Object::JavascriptTemplate() {
00227                     js_template->SetAccessor(v8::String::NewSymbol("from"), accessor_getter<Area, &Area::js_from>);
00228                     js_template->SetAccessor(v8::String::NewSymbol("geom"), accessor_getter<Area, &Area::js_geom>);
00229                 }
00230 
00231             };
00232 
00233 #endif // OSMIUM_WITH_JAVASCRIPT
00234 
00235         }; // class Area
00236 
00237         /***
00238         * Area created from a way.
00239         * The way pointer given to the constructor will not be stored, all
00240         * needed attributes are copied.
00241         */
00242         class AreaFromWay : public Area {
00243 
00244         public:
00245 
00246             WayNodeList m_node_list;
00247 
00248 #ifdef OSMIUM_WITH_GEOS
00249             AreaFromWay(Way* way, geos::geom::Geometry* geom) : Area(geom) {
00250 #else
00251             AreaFromWay(Way* way) : Area() {
00252 #endif // OSMIUM_WITH_GEOS
00253                 id(way->id());
00254                 version(way->version());
00255                 changeset(way->changeset());
00256                 timestamp(way->timestamp());
00257                 uid(way->uid());
00258                 user(way->user());
00259 
00260                 tags(way->tags());
00261                 m_node_list = way->nodes();
00262             }
00263 
00264             ~AreaFromWay() {
00265             }
00266 
00267             const WayNodeList& nodes() const {
00268                 return m_node_list;
00269             }
00270 
00271             WayNodeList& nodes() {
00272                 return m_node_list;
00273             }
00274 
00275             osm_object_type_t get_type() const {
00276                 return AREA_FROM_WAY;
00277             }
00278 
00279         }; // class AreaFromWay
00280 
00281         /***
00282         * Area created from a relation with tag type=multipolygon or type=boundary
00283         */
00284         class AreaFromRelation : public Area {
00285 
00286             bool boundary; 
00287 
00289             Relation *relation;
00290 
00292             std::vector<Osmium::OSM::Way> member_ways;
00293 
00295             int num_ways;
00296 
00298             int missing_ways;
00299 
00300             std::string geometry_error_message;
00301 
00303             void (*callback)(Osmium::OSM::Area*);
00304 
00306             bool attempt_repair;
00307 
00308             bool ignore_tag(const std::string &s) {
00309                 if (s=="type") return true;
00310                 if (s=="created_by") return true;
00311                 if (s=="source") return true;
00312                 if (s=="note") return true;
00313                 return false;
00314             }
00315 
00316             bool same_tags(const Object* a, const Object* b) {
00317                 if ((a == NULL) || (b == NULL)) return false;
00318                 const TagList& at = a->tags();
00319                 const TagList& bt = b->tags();
00320                 std::map<std::string, std::string> aTags;
00321 
00322                 TagList::const_iterator end = at.end();
00323                 for (TagList::const_iterator it = at.begin(); it != end; ++it) {
00324                     if (ignore_tag(it->key())) continue;
00325                     aTags[it->key()] = it->value();
00326                 }
00327                 end = bt.end();
00328                 for (TagList::const_iterator it = bt.begin(); it != end; ++it) {
00329                     if (ignore_tag(it->key())) continue;
00330                     if (aTags[it->key()] != it->value()) return false;
00331                     aTags.erase(it->key());
00332                 }
00333                 if (!aTags.empty()) return false;
00334                 return true;
00335             }
00336 
00338             bool merge_tags(Object* a, const Object* b) {
00339                 bool rv = true;
00340                 TagList& at = a->tags();
00341                 const TagList& bt = b->tags();
00342                 std::map<std::string, std::string> aTags;
00343                 TagList::const_iterator end = at.end();
00344                 for (TagList::const_iterator it = at.begin(); it != end; ++it) {
00345                     if (ignore_tag(it->key())) continue;
00346                     aTags[it->key()] = it->value();
00347                 }
00348                 end = bt.end();
00349                 for (TagList::const_iterator it = bt.begin(); it != end; ++it) {
00350                     if (ignore_tag(it->key())) continue;
00351                     if (aTags.find(it->key()) != aTags.end()) {
00352                         if (aTags[it->key()] != it->value()) rv = false;
00353                     } else {
00354                         at.add(it->key(), it->value());
00355                         aTags[it->key()] = it->value();
00356                     }
00357                 }
00358                 return rv;
00359             }
00360 
00361             bool untagged(const Object* r) {
00362                 if (r == NULL) return true;
00363                 const TagList& tags = r->tags();
00364                 if (tags.empty()) return true;
00365                 TagList::const_iterator end = tags.end();
00366                 for (TagList::const_iterator it = tags.begin(); it != end; ++it) {
00367                     if (! ignore_tag(it->key()) ) {
00368                         return false;
00369                     }
00370                 }
00371                 return true;
00372             }
00373 
00374         public:
00375 
00376             AreaFromRelation(Relation* r, bool b, int n, void (*callback)(Osmium::OSM::Area*), bool repair) : Area(), boundary(b), relation(r), callback(callback) {
00377                 num_ways = n;
00378                 missing_ways = n;
00379 #ifdef OSMIUM_WITH_GEOS
00380                 geometry = NULL;
00381 #endif // OSMIUM_WITH_GEOS
00382                 id(r->id());
00383                 attempt_repair = repair;
00384             }
00385 
00386 #ifdef OSMIUM_WITH_MULTIPOLYGON_PROFILING
00387             static std::vector<std::pair<std::string, timer *> > timers;
00388 
00389             static timer write_complex_poly_timer;
00390             static timer assemble_ways_timer;
00391             static timer assemble_nodes_timer;
00392             static timer make_one_ring_timer;
00393             static timer mor_polygonizer_timer;
00394             static timer mor_union_timer;
00395             static timer contains_timer;
00396             static timer extra_polygons_timer;
00397             static timer polygon_build_timer;
00398             static timer inner_ring_touch_timer;
00399             static timer polygon_intersection_timer;
00400             static timer multipolygon_build_timer;
00401             static timer multipolygon_write_timer;
00402             static timer error_write_timer;
00403 
00404             static void init_timings() {
00405                 timers.push_back(std::pair<std::string, timer *> ("   thereof assemble_ways", &assemble_ways_timer));
00406                 timers.push_back(std::pair<std::string, timer *> ("   thereof make_one_ring", &make_one_ring_timer));
00407                 timers.push_back(std::pair<std::string, timer *> ("      thereof union", &mor_union_timer));
00408                 timers.push_back(std::pair<std::string, timer *> ("      thereof polygonizer", &mor_polygonizer_timer));
00409                 timers.push_back(std::pair<std::string, timer *> ("   thereof contains", &contains_timer));
00410                 timers.push_back(std::pair<std::string, timer *> ("   thereof extra_polygons", &extra_polygons_timer));
00411                 timers.push_back(std::pair<std::string, timer *> ("   thereof polygon_build", &polygon_build_timer));
00412                 timers.push_back(std::pair<std::string, timer *> ("      thereof inner_ring_touch", &inner_ring_touch_timer));
00413                 timers.push_back(std::pair<std::string, timer *> ("      thereof intersections", &polygon_intersection_timer));
00414                 timers.push_back(std::pair<std::string, timer *> ("   thereof multipolygon_build", &multipolygon_build_timer));
00415                 timers.push_back(std::pair<std::string, timer *> ("   thereof multipolygon_write", &multipolygon_write_timer));
00416                 timers.push_back(std::pair<std::string, timer *> ("   thereof error_write", &error_write_timer));
00417             }
00418 
00419             static void print_timings() {
00420                 for (unsigned int i=0; i<timers.size(); i++) {
00421                     std::cerr << timers[i].first << ": " << *(timers[i].second) << std::endl;
00422                 }
00423             }
00424 #endif // OSMIUM_WITH_MULTIPOLYGON_PROFILING
00425 
00426             ~AreaFromRelation() {
00427                 delete relation;
00428                 member_ways.erase(member_ways.begin(), member_ways.end());
00429             }
00430 
00431             osm_object_type_t get_type() const {
00432                 return AREA_FROM_RELATION;
00433             }
00434 
00436             void add_member_way(Osmium::OSM::Way *way) {
00437                 member_ways.push_back(*way);
00438                 missing_ways--;
00439             }
00440 
00442             bool is_complete() {
00443                 return missing_ways == 0;
00444             }
00445 
00446             void handle_complete_multipolygon() {
00447 #ifdef OSMIUM_WITH_GEOS
00448                 if (! build_geometry()) {
00449                     std::cerr << "  geom build error: " << geometry_error_message << "\n";
00450                 }
00451 #endif // OSMIUM_WITH_GEOS
00452                 callback(this);
00453             }
00454 
00455         private:
00456 
00468 #ifdef OSMIUM_WITH_GEOS
00469             geos::geom::LinearRing *create_non_intersecting_linear_ring(geos::geom::CoordinateSequence *orig_cs) {
00470                 const std::vector<geos::geom::Coordinate>* coords = orig_cs->toVector();
00471                 int inv = coords->size();
00472                 int val = 0;
00473                 int current = (inv + val) / 2;
00474                 bool simple;
00475 
00476                 // find the longest non-intersecting stretch from the beginning
00477                 // of the way.
00478                 while (1) {
00479                     std::vector<geos::geom::Coordinate> *vv = new std::vector<geos::geom::Coordinate>(coords->begin(), coords->begin() + current);
00480                     geos::geom::CoordinateSequence *cs = geos::geom::CoordinateArraySequenceFactory::instance()->create(vv);
00481                     geos::geom::LineString *a = Osmium::Geometry::geos_geometry_factory()->createLineString(cs);
00482                     if (!(simple = a->isSimple())) {
00483                         inv = current;
00484                     } else {
00485                         val = current;
00486                     }
00487                     delete a;
00488                     if (current == (inv+val)/2) break;
00489                     current = (inv + val) / 2;
00490                 }
00491 
00492                 if (!simple) current--;
00493 
00494                 unsigned int cutoutstart = current;
00495 
00496                 inv = 0;
00497                 val = coords->size();
00498                 current = (inv + val) / 2;
00499 
00500                 // find the longest non-intersecting stretch from the end
00501                 // of the way. Note that this is likely to overlap with the
00502                 // stretch found above - assume a 10-node way where nodes 3
00503                 // and 7 are identical, then we will find the sequence 0..6
00504                 // above, and 4..9 here!
00505 
00506                 while (1) {
00507                     std::vector<geos::geom::Coordinate> *vv = new std::vector<geos::geom::Coordinate>(coords->begin() + current, coords->end());
00508                     geos::geom::CoordinateSequence *cs = geos::geom::CoordinateArraySequenceFactory::instance()->create(vv);
00509                     geos::geom::LineString *a = Osmium::Geometry::geos_geometry_factory()->createLineString(cs);
00510                     if (!(simple = a->isSimple())) {
00511                         inv = current;
00512                     } else {
00513                         val = current;
00514                     }
00515                     delete a;
00516                     if (current == (inv+val)/2) break;
00517                     current = (inv + val) / 2;
00518                 }
00519                 if (!simple) current++;
00520                 unsigned int cutoutend = current;
00521 
00522                 // assemble a new linear ring by cutting out the problematic bit.
00523                 // if the "problematic bit" however is longer than half the way,
00524                 // then try using the "problematic bit" by itself.
00525 
00526                 std::vector<geos::geom::Coordinate> *vv = new std::vector<geos::geom::Coordinate>();
00527                 if (cutoutstart < cutoutend) {
00528                     unsigned int t = cutoutstart;
00529                     cutoutstart = cutoutend;
00530                     cutoutend = t;
00531                 }
00532                 if (cutoutstart-cutoutend > coords->size() / 2) {
00533                     vv->insert(vv->end(), coords->begin() + cutoutend, coords->begin() + cutoutstart);
00534                     vv->insert(vv->end(), vv->at(0));
00535                 } else {
00536                     vv->insert(vv->end(), coords->begin(), coords->begin() + cutoutend);
00537                     vv->insert(vv->end(), coords->begin() + cutoutstart, coords->end());
00538                 }
00539                 geos::geom::CoordinateSequence *cs = geos::geom::CoordinateArraySequenceFactory::instance()->create(vv);
00540                 geos::geom::LinearRing *a = Osmium::Geometry::geos_geometry_factory()->createLinearRing(cs);
00541 
00542                 // if this results in a valid ring, return it; else return NULL.
00543 
00544                 if (!a->isValid()) return NULL;
00545                 geos::geom::LinearRing *b = dynamic_cast<geos::geom::LinearRing *>(a->clone());
00546                 //delete a;
00547                 return b;
00548             }
00549 
00557             RingInfo *make_one_ring(std::vector<WayInfo *> &ways, osm_object_id_t first, osm_object_id_t last, int ringcount, int sequence) {
00558 
00559                 // have we found a loop already?
00560                 if (first && first == last) {
00561                     geos::geom::CoordinateSequence *cs = geos::geom::CoordinateArraySequenceFactory::instance()->create((size_t)0, (size_t)0);
00562                     geos::geom::LinearRing *lr = NULL;
00563                     try {
00564                         START_TIMER(mor_polygonizer);
00565                         WayInfo **sorted_ways = new WayInfo*[sequence];
00566                         for (unsigned int i=0; i<ways.size(); i++) {
00567                             if (ways[i]->used == ringcount) {
00568                                 sorted_ways[ways[i]->sequence] = ways[i];
00569                             }
00570                         }
00571                         for (int i=0; i<sequence; i++) {
00572                             cs->add(dynamic_cast<geos::geom::LineString *>(sorted_ways[i]->way_geom)->getCoordinatesRO(), false, !sorted_ways[i]->invert);
00573                         }
00574                         delete[] sorted_ways;
00575                         lr = Osmium::Geometry::geos_geometry_factory()->createLinearRing(cs);
00576                         STOP_TIMER(mor_polygonizer);
00577                         if (!lr->isSimple() || !lr->isValid()) {
00578                             //delete lr;
00579                             lr = NULL;
00580                             if (attempt_repair) {
00581                                 lr = create_non_intersecting_linear_ring(cs);
00582                                 if (lr) {
00583                                     if (Osmium::debug())
00584                                         std::cerr << "successfully repaired an invalid ring" << std::endl;
00585                                 }
00586                             }
00587                             if (!lr) return NULL;
00588                         }
00589                         bool ccw = geos::algorithm::CGAlgorithms::isCCW(lr->getCoordinatesRO());
00590                         RingInfo *rl = new RingInfo();
00591                         rl->direction = ccw ? COUNTERCLOCKWISE : CLOCKWISE;
00592                         rl->polygon = Osmium::Geometry::geos_geometry_factory()->createPolygon(lr, NULL);
00593                         return rl;
00594                     } catch (const geos::util::GEOSException& exc) {
00595                         if (Osmium::debug())
00596                             std::cerr << "Exception: " << exc.what() << std::endl;
00597                         return NULL;
00598                     }
00599                 }
00600 
00601                 // have we not allocated anything yet, then simply start with first available way,
00602                 // or return NULL if all are taken.
00603                 if (!first) {
00604                     for (unsigned int i=0; i<ways.size(); i++) {
00605                         if (ways[i]->used != -1) continue;
00606                         ways[i]->used = ringcount;
00607                         ways[i]->sequence = 0;
00608                         ways[i]->invert = false;
00609                         RingInfo *rl = make_one_ring(ways, ways[i]->firstnode, ways[i]->lastnode, ringcount, 1);
00610                         if (rl) {
00611                             rl->ways.push_back(ways[i]);
00612                             return rl;
00613                         }
00614                         ways[i]->used = -2;
00615                         break;
00616                     }
00617                     return NULL;
00618                 }
00619 
00620                 // try extending our current line at the rear end
00621                 // since we are looking for a LOOP, no sense to try extending it at both ends
00622                 // as we'll eventually get there anyway!
00623 
00624                 for (unsigned int i=0; i<ways.size(); i++) {
00625                     if (ways[i]->used < 0) ways[i]->tried = false;
00626                 }
00627 
00628                 for (unsigned int i=0; i<ways.size(); i++) {
00629                     // ignore used ways
00630                     if (ways[i]->used >= 0) continue;
00631                     if (ways[i]->tried) continue;
00632                     ways[i]->tried = true;
00633 
00634                     int old_used = ways[i]->used;
00635                     if (ways[i]->firstnode == last) {
00636                         // add way to end
00637                         ways[i]->used = ringcount;
00638                         ways[i]->sequence = sequence;
00639                         ways[i]->invert = false;
00640                         RingInfo *result = make_one_ring(ways, first, ways[i]->lastnode, ringcount, sequence+1);
00641                         if (result) {
00642                             result->ways.push_back(ways[i]);
00643                             return result;
00644                         }
00645                         ways[i]->used = old_used;
00646                     } else if (ways[i]->lastnode == last) {
00647                         // add way to end, but turn it around
00648                         ways[i]->used = ringcount;
00649                         ways[i]->sequence = sequence;
00650                         ways[i]->invert = true;
00651                         RingInfo *result = make_one_ring(ways, first, ways[i]->firstnode, ringcount, sequence+1);
00652                         if (result) {
00653                             result->ways.push_back(ways[i]);
00654                             return result;
00655                         }
00656                         ways[i]->used = old_used;
00657                     }
00658                 }
00659                 // we have exhausted all combinations.
00660                 return NULL;
00661             }
00662 
00674             bool find_and_repair_holes_in_rings(std::vector<WayInfo *> *ways) {
00675                 // collect the remaining debris (=unused ways) and find dangling nodes.
00676 
00677                 std::map<int,geos::geom::Point *> dangling_node_map;
00678                 for (std::vector<WayInfo *>::iterator i(ways->begin()); i != ways->end(); i++) {
00679                     if ((*i)->used < 0) {
00680                         (*i)->innerouter = UNSET;
00681                         (*i)->used = -1;
00682                         for (int j=0; j<2; j++) {
00683                             int nid = j ? (*i)->firstnode : (*i)->lastnode;
00684                             if (dangling_node_map[nid]) {
00685                                 delete dangling_node_map[nid];
00686                                 dangling_node_map[nid] = NULL;
00687                             } else {
00688                                 dangling_node_map[nid] = j ? (*i)->get_firstnode_geom() : (*i)->get_lastnode_geom();
00689                             }
00690                         }
00691                     }
00692                 }
00693 
00694                 do {
00695                     int mindist_id = 0;
00696                     double mindist = -1;
00697                     int node1_id = 0;
00698                     geos::geom::Point *node1 = NULL;
00699                     geos::geom::Point *node2 = NULL;
00700 
00701                     // find one pair consisting of a random node from the list (node1)
00702                     // plus the node that lies closest to it.
00703                     for (std::map<int,geos::geom::Point *>::iterator i(dangling_node_map.begin()); i!= dangling_node_map.end(); i++) {
00704                         if (!i->second) continue;
00705                         if (node1 == NULL) {
00706                             node1 = i->second;
00707                             node1_id = i->first;
00708                             i->second = NULL;
00709                             mindist = -1;
00710                         } else {
00711 # if GEOS_VERSION_MAJOR < 3 || (GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR <= 2)
00712                             double dist = geos::operation::distance::DistanceOp::distance(node1, (i->second)); // deprecated in newer version of GEOS
00713 # else
00714                             double dist = geos::operation::distance::DistanceOp::distance(*node1, *(i->second));
00715 # endif
00716                             if ((dist < mindist) || (mindist < 0)) {
00717                                 mindist = dist;
00718                                 mindist_id = i->first;
00719                             }
00720                         }
00721                     }
00722 
00723                     // if such a pair has been found, synthesize a connecting way.
00724                     if (node1 && mindist > -1) {
00725                         // if we find that there are dangling nodes but aren't
00726                         // repairing - break out.
00727                         if (!attempt_repair) return false;
00728 
00729                         // drop node2 from dangling map
00730                         node2 = dangling_node_map[mindist_id];
00731                         dangling_node_map[mindist_id] = NULL;
00732 
00733                         std::vector<geos::geom::Coordinate> *c = new std::vector<geos::geom::Coordinate>;
00734                         c->push_back(*(node1->getCoordinate()));
00735                         c->push_back(*(node2->getCoordinate()));
00736                         geos::geom::CoordinateSequence *cs = Osmium::Geometry::geos_geometry_factory()->getCoordinateSequenceFactory()->create(c);
00737                         geos::geom::Geometry *geometry = (geos::geom::Geometry *) Osmium::Geometry::geos_geometry_factory()->createLineString(cs);
00738                         ways->push_back(new WayInfo(geometry, node1_id, mindist_id, UNSET));
00739                         if (Osmium::debug())
00740                             std::cerr << "fill gap between nodes " << node1_id << " and " << mindist_id << std::endl;
00741                     } else {
00742                         break;
00743                     }
00744                 } while (1);
00745 
00746                 return true;
00747             }
00748 
00749 
00754             bool build_geometry() {
00755                 std::vector<WayInfo *> ways;
00756 
00757                 // the timestamp of the multipolygon will be the maximum of the timestamp from the relation and from all member ways
00758                 timestamp(relation->timestamp());
00759 
00760                 // assemble all ways which are members of this relation into a
00761                 // vector of WayInfo elements. this holds room for the way pointer
00762                 // and some extra flags.
00763 
00764                 START_TIMER(assemble_ways);
00765                 for (std::vector<Way>::iterator i(member_ways.begin()); i != member_ways.end(); i++) {
00766                     if (i->timestamp() > timestamp()) {
00767                         timestamp(i->timestamp());
00768                     }
00769                     WayInfo *wi = new WayInfo(&(*i), UNSET);
00770                     if (wi->way_geom) {
00771                         geos::io::WKTWriter wkt;
00772                     } else {
00773                         delete wi;
00774                         return geometry_error("invalid way geometry in multipolygon relation member");
00775                     }
00776                     ways.push_back(wi);
00777                     // TODO drop duplicate ways automatically in repair mode?
00778                     // TODO maybe add INNER/OUTER instead of UNSET to enable later warnings on role mismatch
00779                 }
00780                 STOP_TIMER(assemble_ways);
00781 
00782                 std::vector<RingInfo *> ringlist;
00783 
00784                 // convenience defines to aid in clearing up on error return.
00785 #define clear_ringlist() \
00786                 for (std::vector<RingInfo *>::const_iterator rli(ringlist.begin()); rli != ringlist.end(); rli++) delete *rli;
00787 #define clear_wayinfo() \
00788                 for (std::vector<WayInfo *>::const_iterator win(ways.begin()); win != ways.end(); win++) delete *win;
00789 
00790                 // try and create as many closed rings as possible from the assortment
00791                 // of ways. make_one_ring will automatically flag those that have been
00792                 // used so they are not used again.
00793 
00794                 do {
00795                     START_TIMER(make_one_ring);
00796                     RingInfo *r = make_one_ring(ways, 0, 0, ringlist.size(), 0);
00797                     STOP_TIMER(make_one_ring);
00798                     if (r == NULL) break;
00799                     r->ring_id = ringlist.size();
00800                     ringlist.push_back(r);
00801                 } while (1);
00802 
00803                 if (ringlist.empty()) {
00804                     // FIXME return geometry_error("no rings");
00805                 }
00806 
00807                 if (!find_and_repair_holes_in_rings(&ways)) {
00808                     clear_ringlist();
00809                     clear_wayinfo();
00810                     return geometry_error("un-connectable dangling ends");
00811                 }
00812 
00813                 // re-run ring building, taking into account the newly created "repair" bits.
00814                 // (in case there were no dangling bits, make_one_ring terminates quickly.)
00815                 do {
00816                     START_TIMER(make_one_ring);
00817                     RingInfo *r = make_one_ring(ways, 0, 0, ringlist.size(), 0);
00818                     STOP_TIMER(make_one_ring);
00819                     if (r == NULL) break;
00820                     r->ring_id = ringlist.size();
00821                     ringlist.push_back(r);
00822                 } while (1);
00823 
00824                 if (ringlist.empty()) {
00825                     clear_ringlist();
00826                     clear_wayinfo();
00827                     return geometry_error("no rings");
00828                 }
00829 
00830                 std::vector<geos::geom::Geometry *> *polygons = new std::vector<geos::geom::Geometry *>();
00831 
00832                 geos::geom::MultiPolygon *mp = NULL;
00833 
00834                 // find out which ring contains which other ring, so we know
00835                 // which are inner rings and which outer. don't trust the "role"
00836                 // specifications.
00837 
00838                 START_TIMER(contains);
00839 
00840                 bool **contains = new bool*[ringlist.size()];
00841                 bool *contained_by_even_number = new bool[ringlist.size()];
00842 
00843                 // reset array
00844                 for (unsigned int i=0; i<ringlist.size(); i++) {
00845                     contains[i] = new bool[ringlist.size()];
00846                     contained_by_even_number[i] = true;
00847                     for (unsigned int j=0; j<ringlist.size(); j++) {
00848                         contains[i][j] = false;
00849                     }
00850                 }
00851 
00852                 // build contains relationships.
00853                 // we use contained_by_even_number as a helper for us to determine
00854                 // whether something is an inner (false) or outer (true) ring.
00855 
00856                 for (unsigned int i=0; i<ringlist.size(); i++) {
00857                     const geos::geom::prep::PreparedPolygon *pp = new geos::geom::prep::PreparedPolygon(ringlist[i]->polygon);
00858                     for (unsigned int j=0; j<ringlist.size(); j++) {
00859                         if (i==j) continue;
00860                         if (contains[j][i]) continue;
00861                         contains[i][j] = pp->contains(ringlist[j]->polygon);
00862                         contained_by_even_number[j] ^= contains[i][j];
00863                     }
00864                     delete pp;
00865                 }
00866 
00867                 // we now have an array that has a true value whenever something is
00868                 // contained by something else; if a contains b and b contains c, then
00869                 // our array says that a contains b, b contains c, and a contains c.
00870                 // thin out the array so that only direct relationships remain (and
00871                 // the "a contains c" is dropped).
00872 
00873                 for (unsigned int i=0; i<ringlist.size(); i++) {
00874                     for (unsigned j=0; j<ringlist.size(); j++) {
00875                         if (contains[i][j]) {
00876                             // see if there is an intermediary relationship
00877                             for (unsigned int k=0; k<ringlist.size(); k++) {
00878                                 if (k==i) continue;
00879                                 if (k==j) continue;
00880                                 if (contains[i][k] && contains[k][j]) {
00881                                     // intermediary relationship exists; break this
00882                                     // one up.
00883                                     contains[i][j] = false;
00884                                     ringlist[j]->nested = true;
00885                                     break;
00886                                 }
00887                             }
00888                         }
00889                     }
00890                 }
00891 
00892                 // populate the "inner_rings" list and the "contained_by" pointer
00893                 // in the ring list based on the data collected. the "contains"
00894                 // array can be thrown away afterwards.
00895 
00896                 for (unsigned int i=0; i<ringlist.size(); i++) {
00897                     for (unsigned int j=0; j<ringlist.size(); j++) {
00898                         if (contains[i][j] && !contained_by_even_number[j]) {
00899                             ringlist[j]->contained_by = ringlist[i];
00900                             ringlist[i]->inner_rings.push_back(ringlist[j]);
00901                         }
00902                     }
00903                     delete[] contains[i];
00904                 }
00905 
00906                 delete[] contains;
00907                 delete[] contained_by_even_number;
00908                 STOP_TIMER(contains);
00909 
00910                 // now look at all enclosed (inner) rings that consist of only one way.
00911                 // if such an inner ring has way tags, do the following:
00912                 // * emit an extra polygon for the inner ring if the tags are different
00913                 //   from the relation's
00914                 // * emit a warning, and ignore the inner ring, if the tags are the same
00915                 //   as for the relation
00916 
00917                 START_TIMER(extra_polygons);
00918                 for (unsigned int i=0; i<ringlist.size(); i++) {
00919                     if (ringlist[i]->contained_by) {
00920                         if (ringlist[i]->ways.size() == 1 && !untagged(ringlist[i]->ways[0]->way)) {
00921                             std::vector<geos::geom::Geometry *> *g = new std::vector<geos::geom::Geometry *>;
00922                             if (ringlist[i]->direction == CLOCKWISE) {
00923                                 g->push_back(ringlist[i]->polygon->clone());
00924                             } else {
00925                                 geos::geom::LineString *tmp = dynamic_cast<geos::geom::LineString *>(ringlist[i]->polygon->getExteriorRing()->reverse());
00926                                 geos::geom::LinearRing *reversed_ring =
00927                                     Osmium::Geometry::geos_geometry_factory()->createLinearRing(tmp->getCoordinates());
00928                                 delete tmp;
00929                                 g->push_back(Osmium::Geometry::geos_geometry_factory()->createPolygon(reversed_ring, NULL));
00930                             }
00931 
00932                             geos::geom::MultiPolygon *special_mp = Osmium::Geometry::geos_geometry_factory()->createMultiPolygon(g);
00933 
00934                             if (same_tags(ringlist[i]->ways[0]->way, relation)) {
00935                                 // warning
00936                                 // warnings.insert("duplicate_tags_on_inner");
00937                             } else if (ringlist[i]->contained_by->ways.size() == 1 && same_tags(ringlist[i]->ways[0]->way, ringlist[i]->contained_by->ways[0]->way)) {
00938                                 // warning
00939                                 // warnings.insert("duplicate_tags_on_inner");
00940                             } else {
00941                                 Osmium::OSM::AreaFromWay *internal_mp =
00942                                     new Osmium::OSM::AreaFromWay(ringlist[i]->ways[0]->way, special_mp);
00943                                 callback(internal_mp);
00944                                 delete internal_mp;
00945                                 // AreaFromWay destructor deletes the
00946                                 // geometry, so avoid to delete it again.
00947                                 special_mp = NULL;
00948                             }
00949                             delete special_mp;
00950                         }
00951                     }
00952                 }
00953                 STOP_TIMER(extra_polygons);
00954 
00955                 // for all non-enclosed rings, assemble holes and build polygon.
00956 
00957                 START_TIMER(polygon_build)
00958                 for (unsigned int i=0; i<ringlist.size(); i++) {
00959                     // look only at outer, i.e. non-contained rings. each ends up as one polygon.
00960                     if (ringlist[i] == NULL) continue; // can happen if ring has been deleted
00961                     if (ringlist[i]->contained_by) continue;
00962 
00963                     std::vector<geos::geom::Geometry *> *holes = new std::vector<geos::geom::Geometry *>(); // ownership is later transferred to polygon
00964 
00965                     START_TIMER(inner_ring_touch)
00966                     for (int j=0; j<((int)ringlist[i]->inner_rings.size()-1); j++) {
00967                         if (!ringlist[i]->inner_rings[j]->polygon) continue;
00968                         geos::geom::LinearRing *ring = (geos::geom::LinearRing *) ringlist[i]->inner_rings[j]->polygon->getExteriorRing();
00969 
00970                         // check if some of the rings touch another ring.
00971 
00972                         for (unsigned int k=j + 1; k<ringlist[i]->inner_rings.size(); k++) {
00973                             if (!ringlist[i]->inner_rings[k]->polygon) continue;
00974                             const geos::geom::Geometry *compare = ringlist[i]->inner_rings[k]->polygon->getExteriorRing();
00975                             geos::geom::Geometry *inter = NULL;
00976                             try {
00977                                 if (!ring->intersects(compare)) continue;
00978                                 inter = ring->intersection(compare);
00979                             } catch (const geos::util::GEOSException& exc) {
00980                                 // nop;
00981                             }
00982                             if (inter && (inter->getGeometryTypeId() == geos::geom::GEOS_LINESTRING || inter->getGeometryTypeId() == geos::geom::GEOS_MULTILINESTRING)) {
00983                                 // touching inner rings
00984                                 // this is allowed, but we must fix them up into a valid
00985                                 // geometry
00986                                 geos::geom::Geometry *diff = ring->symDifference(compare);
00987                                 geos::operation::polygonize::Polygonizer *p = new geos::operation::polygonize::Polygonizer();
00988                                 p->add(diff);
00989                                 std::vector<geos::geom::Polygon*>* polys = p->getPolygons();
00990                                 if (polys && polys->size() == 1) {
00991                                     ringlist[i]->inner_rings[j]->polygon = polys->at(0);
00992                                     bool ccw = geos::algorithm::CGAlgorithms::isCCW(polys->at(0)->getExteriorRing()->getCoordinatesRO());
00993                                     ringlist[i]->inner_rings[j]->direction = ccw ? COUNTERCLOCKWISE : CLOCKWISE;
00994                                     ringlist[i]->inner_rings[k]->polygon = NULL;
00995                                     j=-1;
00996                                     break;
00997                                 }
00998                             } else {
00999                                 // other kind of intersect between inner rings; this is
01000                                 // not allwoed and will lead to an exception later when
01001                                 // building the MP
01002                             }
01003                         }
01004                     }
01005                     STOP_TIMER(inner_ring_touch)
01006 
01007                     for (unsigned int j=0; j<ringlist[i]->inner_rings.size(); j++) {
01008                         if (!ringlist[i]->inner_rings[j]->polygon) continue;
01009                         geos::geom::LinearRing *ring = (geos::geom::LinearRing *) ringlist[i]->inner_rings[j]->polygon->getExteriorRing();
01010 
01011                         if (ringlist[i]->inner_rings[j]->direction == CLOCKWISE) {
01012                             // reverse ring
01013                             geos::geom::LineString *tmp = dynamic_cast<geos::geom::LineString *>(ring->reverse());
01014                             geos::geom::LinearRing *reversed_ring =
01015                                 Osmium::Geometry::geos_geometry_factory()->createLinearRing(tmp->getCoordinates());
01016                             delete tmp;
01017                             holes->push_back(reversed_ring);
01018                         } else {
01019                             holes->push_back(ring);
01020                         }
01021                     }
01022 
01023                     geos::geom::LinearRing *ring = (geos::geom::LinearRing *) ringlist[i]->polygon->getExteriorRing();
01024                     if (ringlist[i]->direction == COUNTERCLOCKWISE) {
01025                         geos::geom::LineString *tmp = dynamic_cast<geos::geom::LineString *>(ring->reverse());
01026                         geos::geom::LinearRing *reversed_ring = Osmium::Geometry::geos_geometry_factory()->createLinearRing(tmp->getCoordinates());
01027                         ring = reversed_ring;
01028                         delete tmp;
01029                     } else {
01030                         ring = dynamic_cast<geos::geom::LinearRing *>(ring->clone());
01031                     }
01032                     delete ringlist[i]->polygon;
01033                     ringlist[i]->polygon = NULL;
01034                     geos::geom::Polygon *p = NULL;
01035                     bool valid = false;
01036 
01037                     try {
01038                         p = Osmium::Geometry::geos_geometry_factory()->createPolygon(ring, holes);
01039                         if (p) valid = p->isValid();
01040                     } catch (const geos::util::GEOSException& exc) {
01041                         // nop
01042                         if (Osmium::debug())
01043                             std::cerr << "Exception during creation of polygon for relation #" << relation->id() << ": " << exc.what() << " (treating as invalid polygon)" << std::endl;
01044                     }
01045                     if (!valid) {
01046                         // polygon is invalid.
01047                         clear_ringlist();
01048                         clear_wayinfo();
01049                         if (p) delete p;
01050                         else delete ring;
01051                         return geometry_error("invalid ring");
01052                     } else {
01053                         polygons->push_back(p);
01054                         for (unsigned int k=0; k<ringlist[i]->ways.size(); k++) {
01055                             WayInfo *wi = ringlist[i]->ways[k];
01056                             // may have "hole filler" ways in there, not backed by
01057                             // proper way and thus no tags:
01058                             if (wi->way == NULL) continue;
01059                             if (untagged(wi->way)) {
01060                                 // way not tagged - ok
01061                             } else if (same_tags(relation, wi->way)) {
01062                                 // way tagged the same as relation/previous ways, ok
01063                             } else if (untagged(relation)) {
01064                                 // relation untagged; use tags from way; ok
01065                                 merge_tags(relation, wi->way);
01066                             }
01067 
01068                             wi->innerouter = OUTER;
01069                             if (wi->orig_innerouter == INNER && wi->errorhint.empty()) {
01070                                 // warning: inner/outer mismatch
01071                             }
01072                         }
01073                         // copy tags from relation into area
01074                         tags(relation->tags());
01075                     }
01076                     // later delete ringlist[i];
01077                     // ringlist[i] = NULL;
01078                 }
01079                 STOP_TIMER(polygon_build);
01080 
01081                 clear_ringlist();
01082                 clear_wayinfo();
01083                 if (polygons->empty()) {
01084                     return geometry_error("no rings");
01085                 }
01086 
01087                 START_TIMER(multipolygon_build);
01088                 bool valid = false;
01089                 try {
01090                     mp = Osmium::Geometry::geos_geometry_factory()->createMultiPolygon(polygons);
01091                     valid = mp->isValid();
01092                 } catch (const geos::util::GEOSException& exc) {
01093                     // nop
01094                 };
01095                 STOP_TIMER(multipolygon_build);
01096                 if (valid) {
01097                     geometry = mp;
01098                     return true;
01099                 }
01100                 return geometry_error("multipolygon invalid");
01101             }
01102 
01103             bool geometry_error(const char *message) {
01104                 geometry_error_message = message;
01105                 if (Osmium::debug()) {
01106                     std::cerr << "building mp failed: " << geometry_error_message << std::endl;
01107                 }
01108                 geometry = NULL;
01109                 return false;
01110             }
01111 
01112 #endif // OSMIUM_WITH_GEOS
01113 
01114         }; // class AreaFromRelation
01115 
01116     } // namespace OSM
01117 
01118 } // namespace Osmium
01119 
01120 #endif // OSMIUM_OSM_AREA_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines