Geometry.h

Go to the documentation of this file.

#ifndef OSMSCOUT_UTIL_GEOMETRY_H
#define OSMSCOUT_UTIL_GEOMETRY_H
 
/*
  This source is part of the libosmscout library
  Copyright (C) 2009  Tim Teulings
 
  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
 
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
 
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA
*/
 
#include <algorithm>
#include <array>
#include <functional>
#include <list>
#include <unordered_map>
#include <utility>
#include <vector>
#include <tuple>
 
#include <osmscout/lib/CoreImportExport.h>
 
#include <osmscout/system/Assert.h>
#include <osmscout/system/Math.h>
#include <osmscout/system/SystemTypes.h>
 
#include <osmscout/GeoCoord.h>
#include <osmscout/Point.h>
#include <osmscout/OSMScoutTypes.h>
 
#include <osmscout/util/GeoBox.h>
#include <osmscout/util/Distance.h>
#include <osmscout/util/Bearing.h>
 
namespace osmscout {

 

  inline double DegToRad(double deg)
  {
    return deg*M_PI/180.0;
  }

  inline double RadToDeg(double rad)
  {
    return rad*180.0/M_PI;
  }

  inline double AngleDiff(double a, double b)
  {
    double diff = std::abs(a - b);
 
    if (diff > M_PI) {
      return 2 * M_PI - diff;
    }
 
    return diff;
  }

  template< class InputIt >
  void GetBoundingBox(const InputIt first,
                      const InputIt last,
                      GeoBox& boundingBox)
  {
    assert(first!=last);
 
    double minLon=first->GetLon();
    double maxLon=minLon;
    double minLat=first->GetLat();
    double maxLat=minLat;
 
    for (InputIt i=first; i!=last; i++) {
      minLon=std::min(minLon,i->GetLon());
      maxLon=std::max(maxLon,i->GetLon());
      minLat=std::min(minLat,i->GetLat());
      maxLat=std::max(maxLat,i->GetLat());
    }
 
    boundingBox.Set(GeoCoord(minLat,minLon),
                    GeoCoord(maxLat,maxLon));
  }

  template<typename N>
  GeoBox GetBoundingBox(const std::vector<N>& nodes)
  {
    assert(!nodes.empty());
 
    GeoBox boundingBox;
 
    for (const auto& node : nodes) {
      boundingBox.Include(node);
    }
 
    return boundingBox;
  }

  inline GeoBox GetBoundingBox(const std::vector<Point>& nodes)
  {
    assert(!nodes.empty());
 
    GeoBox boundingBox;
 
    for (const auto& node : nodes) {
      boundingBox.Include(node.GetCoord());
    }
 
    return boundingBox;
  }

  template<typename N>
  void GetBoundingBox(const std::vector<N>& nodes,
                      GeoBox& boundingBox)
  {
    assert(!nodes.empty());
 
    boundingBox.Invalidate();
    for (const auto& node : nodes) {
      boundingBox.Include(node);
    }
  }

  inline void GetBoundingBox(const std::vector<Point>& nodes,
                             GeoBox& boundingBox)
  {
    assert(!nodes.empty());
 
    boundingBox.Invalidate();
    for (const auto& node : nodes) {
      boundingBox.Include(node.GetCoord());
    }
  }

  template<typename N>
  bool LinesIntersect(const N& a1,
                      const N& a2,
                      const N& b1,
                      const N& b2)
  {
    if (a1.IsEqual(b1) ||
        a1.IsEqual(b2) ||
        a2.IsEqual(b1) ||
        a2.IsEqual(b2)) {
      return true;
    }
    if (a1.IsEqual(a2) &&
        b1.IsEqual(b2)){
      // two different zero size vectors can't intersects
      return false;
    }
 
    // denominator, see GetLineIntersection for more details
    double denr=(b2.GetLat()-b1.GetLat())*(a2.GetLon()-a1.GetLon())-
                (b2.GetLon()-b1.GetLon())*(a2.GetLat()-a1.GetLat());
 
    double ua_numr=(b2.GetLon()-b1.GetLon())*(a1.GetLat()-b1.GetLat())-
                   (b2.GetLat()-b1.GetLat())*(a1.GetLon()-b1.GetLon());
    double ub_numr=(a2.GetLon()-a1.GetLon())*(a1.GetLat()-b1.GetLat())-
                   (a2.GetLat()-a1.GetLat())*(a1.GetLon()-b1.GetLon());
 
    if (denr==0.0) { // parallels
      if (ua_numr==0.0 && ub_numr==0.0) {
        // two lines are on one straight line, check the bounds
        GeoBox aBox(GeoCoord(a1.GetLat(),a1.GetLon()),
                    GeoCoord(a2.GetLat(),a2.GetLon()));
        GeoBox bBox(GeoCoord(b1.GetLat(),b1.GetLon()),
                    GeoCoord(b2.GetLat(),b2.GetLon()));
 
        return bBox.Includes(a1,false) ||
               bBox.Includes(a2,false) ||
               aBox.Includes(b1,false) ||
               aBox.Includes(b2,false);
      }
 
      return false;
    }
 
    double ua=ua_numr/denr;
    double ub=ub_numr/denr;
 
    return ua>=0.0 &&
           ua<=1.0 &&
           ub>=0.0 &&
           ub<=1.0;
  }

  template<typename N,typename I>
  bool GetLineIntersection(const N& a1,
                           const N& a2,
                           const N& b1,
                           const N& b2,
                           I& intersection)
  {
    if (a1.IsEqual(b1) ||
        a1.IsEqual(b2)){
      intersection.Set(a1.GetLat(),a1.GetLon());
      return true;
    }
    if (a2.IsEqual(b1) ||
        a2.IsEqual(b2)) {
      intersection.Set(a2.GetLat(),a2.GetLon());
      return true;
    }
    if (a1.IsEqual(a2) &&
        b1.IsEqual(b2)){
      // two different zero size vectors can't intersects
      return false;
    }
 
    // for geographic, expect point.x=lon and point.y=lat
 
    // see https://en.wikipedia.org/wiki/Line%E2%80%93line_intersection
    // and http://www.geeksforgeeks.org/orientation-3-ordered-points/
 
    // denominator < 0  : left angle
    // denominator == 0 : parallels
    // denominator > 0  : right angle
    double denr=(b2.GetLat()-b1.GetLat())*(a2.GetLon()-a1.GetLon())-
                (b2.GetLon()-b1.GetLon())*(a2.GetLat()-a1.GetLat());
 
    double ua_numr=(b2.GetLon()-b1.GetLon())*(a1.GetLat()-b1.GetLat())-
                   (b2.GetLat()-b1.GetLat())*(a1.GetLon()-b1.GetLon());
    double ub_numr=(a2.GetLon()-a1.GetLon())*(a1.GetLat()-b1.GetLat())-
                   (a2.GetLat()-a1.GetLat())*(a1.GetLon()-b1.GetLon());
 
    if (denr==0.0) {
      if (ua_numr==0.0 && ub_numr==0.0) {
        // two lines are on one straight line, check the bounds
        GeoBox aBox(GeoCoord(a1.GetLat(),a1.GetLon()),
                    GeoCoord(a2.GetLat(),a2.GetLon()));
        GeoBox bBox(GeoCoord(b1.GetLat(),b1.GetLon()),
                    GeoCoord(b2.GetLat(),b2.GetLon()));
        if (bBox.Includes(a1,false)){
          intersection.Set(a1.GetLat(),a1.GetLon());
          return true;
        }
        if (bBox.Includes(a2,false)){
          intersection.Set(a2.GetLat(),a2.GetLon());
          return true;
        }
        if (aBox.Includes(b1,false)){
          intersection.Set(b1.GetLat(),b1.GetLon());
          return true;
        }
        if (aBox.Includes(b2,false)){
          intersection.Set(b2.GetLat(),b2.GetLon());
          return true;
        }
 
        return false;
      }
 
      return false;
    }
 
    double ua=ua_numr/denr;
    double ub=ub_numr/denr;
 
    if (ua>=0.0 &&
        ua<=1.0 &&
        ub>=0.0 &&
        ub<=1.0) {
      intersection.Set(a1.GetLat()+ua*(a2.GetLat()-a1.GetLat()),
                       a1.GetLon()+ua*(a2.GetLon()-a1.GetLon()));
      return true;
    }
 
    return false;
  }

  template<typename N>
  bool GetLineIntersectionPixel(const N& a1,
                                const N& a2,
                                const N& b1,
                                const N& b2,
                                N& intersection)
  {
    if (a1.IsEqual(b1) ||
        a1.IsEqual(b2)){
      intersection.Set(a1.GetX(),a1.GetY());
      return true;
    }
    if (a2.IsEqual(b1) ||
        a2.IsEqual(b2)) {
      intersection.Set(a2.GetX(),a2.GetY());
      return true;
    }
    if (a1.IsEqual(a2) &&
        b1.IsEqual(b2)){
      // two different zero size vectors can't intersects
      return false;
    }
 
    double denr=(b2.GetY()-b1.GetY())*(a2.GetX()-a1.GetX())-
                (b2.GetX()-b1.GetX())*(a2.GetY()-a1.GetY());
 
    double ua_numr=(b2.GetX()-b1.GetX())*(a1.GetY()-b1.GetY())-
                   (b2.GetY()-b1.GetY())*(a1.GetX()-b1.GetX());
    double ub_numr=(a2.GetX()-a1.GetX())*(a1.GetY()-b1.GetY())-
                   (a2.GetY()-a1.GetY())*(a1.GetX()-b1.GetX());
 
    if (denr==0.0) {
      // This gives currently false hits because of number resolution problems, if two lines are very
      // close together and for example are part of a very details node curve intersections are detected.
 
      // FIXME: setup intersection
      return ua_numr==0.0 && ub_numr==0.0;
    }
 
    double ua=ua_numr/denr;
    double ub=ub_numr/denr;
 
    if (ua>=0.0 &&
        ua<=1.0 &&
        ub>=0.0 &&
        ub<=1.0) {
      intersection.Set(a1.GetX()+ua*(a2.GetX()-a1.GetX()),
                       a1.GetY()+ua*(a2.GetY()-a1.GetY()));
      return true;
    }
 
    return false;
  }
 
  template<typename N>
  double DistanceSquare(const N& a,
                        const N& b)
  {
    double lonDelta=a.GetLon()-b.GetLon();
    double latDelta=a.GetLat()-b.GetLat();
 
    return lonDelta*lonDelta+latDelta*latDelta;
  }

/*
  template<typename N, typename M>
  inline double IsLeft(const N& p0,
                       const N& p1,
                       const M& p2)
  {
    if ((p2.GetLat()==p0.GetLat() &&
         p2.GetLon()==p0.GetLon()) ||
        (p2.GetLat()==p1.GetLat() &&
         p2.GetLon()==p1.GetLon())) {
      return 0;
    }
 
    return (p1.GetLon()-p0.GetLon())*(p2.GetLat()-p0.GetLat())-
           (p2.GetLon()-p0.GetLon())*(p1.GetLat()-p0.GetLat());
  }
 
  template<typename N, typename M>
  inline bool IsCoordInArea(const N& point,
                            const std::vector<M>& nodes)
  {
    for (int i=0; i<(int)nodes.size()-1; i++) {
      if (point.GetLat()==nodes[i].GetLat() &&
          point.GetLon()==nodes[i].GetLon()) {
        return true;
      }
    }
 
    int wn=0;    // the winding number counter
 
    // loop through all edges of the polygon
    for (int i=0; i<(int)nodes.size()-1; i++) {   // edge from V[i] to V[i+1]
      if (nodes[i].GetLat()<=point.GetLat()) {         // start y <= P.y
        if (nodes[i+1].GetLat()>point.GetLat()) {     // an upward crossing
          if (IsLeft(nodes[i],nodes[i+1],point) > 0) { // P left of edge
            ++wn;            // have a valid up intersect
          }
        }
      }
      else {                       // start y > P.y (no test needed)
        if (nodes[i+1].GetLat()<=point.GetLat()) {    // a downward crossing
          if (IsLeft(nodes[i],nodes[i+1],point) < 0) { // P right of edge
            --wn;            // have a valid down intersect
          }
        }
      }
    }
 
    return wn!=0;
  }*/

  template<typename N, typename M>
  inline bool IsCoordInArea(const N& point,
                            const std::vector<M>& nodes)
  {
    size_t i;
    size_t j;
    bool   c=false;
 
    for (i=0, j=nodes.size()-1; i<nodes.size(); j=i++) {
      if (point.GetLat()==nodes[i].GetLat() &&
          point.GetLon()==nodes[i].GetLon()) {
        return true;
      }
 
      if ((((nodes[i].GetLat()<=point.GetLat()) && (point.GetLat()<nodes[j].GetLat())) ||
           ((nodes[j].GetLat()<=point.GetLat()) && (point.GetLat()<nodes[i].GetLat()))) &&
          (point.GetLon()<(nodes[j].GetLon()-nodes[i].GetLon())*(point.GetLat()-nodes[i].GetLat())/(nodes[j].GetLat()-nodes[i].GetLat())+
           nodes[i].GetLon())) {
        c=!c;
      }
    }
 
    return c;
  }

  template<typename N,typename M>
  inline int GetRelationOfPointToArea(const N& point,
                                      const std::vector<M>& nodes)
  {
    size_t i;
    size_t j;
    bool   c=false;
 
    for (i=0, j=nodes.size()-1; i<nodes.size(); j=i++) {
      if (point.GetLat()==nodes[i].GetLat() &&
          point.GetLon()==nodes[i].GetLon()) {
        return 0;
      }
 
      if ((((nodes[i].GetLat()<=point.GetLat()) && (point.GetLat()<nodes[j].GetLat())) ||
           ((nodes[j].GetLat()<=point.GetLat()) && (point.GetLat()<nodes[i].GetLat()))) &&
          (point.GetLon()<(nodes[j].GetLon()-nodes[i].GetLon())*(point.GetLat()-nodes[i].GetLat())/(nodes[j].GetLat()-nodes[i].GetLat())+
           nodes[i].GetLon())) {
        c=!c;
      }
    }
 
    return c ? 1 : -1;
  }

  template<typename N,typename M>
  inline bool IsAreaCompletelyInArea(const std::vector<N>& a,
                                     const std::vector<M>& b)
  {
    for (const auto& node : a) {
      if (GetRelationOfPointToArea(node,b)<0) {
        return false;
      }
    }
 
    return true;
  }

  template<typename N,typename M>
  inline bool IsAreaAtLeastPartlyInArea(const std::vector<N>& a,
                                        const std::vector<M>& b,
                                        const GeoBox& aBox,
                                        const GeoBox& bBox)
   {
    if (!aBox.Intersects(bBox)){
      return false;
    }
 
    for (const auto& node : a) {
      if (bBox.Includes(node, /*openInterval*/ false) && GetRelationOfPointToArea(node,b)>=0) {
        return true;
      }
    }
 
    return false;
  }

  template<typename N,typename M>
  inline bool IsAreaAtLeastPartlyInArea(const std::vector<N>& a,
                                        const std::vector<M>& b)
  {
    GeoBox aBox;
    GeoBox bBox;
    GetBoundingBox(a, aBox);
    GetBoundingBox(b, bBox);
 
    return IsAreaAtLeastPartlyInArea(a,b,aBox,bBox);
  }

  template<typename N,typename M>
  inline bool IsAreaSubOfArea(const std::vector<N>& a,
                              const std::vector<M>& b)
  {
    for (const auto& node : a) {
      int relPos=GetRelationOfPointToArea(node,b);
 
      if (relPos>0) {
        return true;
      }
 
      if (relPos<0) {
        return false;
      }
    }
 
    return false;
  }

  template<typename N,typename M>
  inline bool IsAreaSubOfAreaQuorum(const std::vector<N>& a,
                                    const std::vector<M>& b)
  {
    size_t pro=0;
    size_t contra=0;
    size_t count=0;
 
    for (const auto& node : a) {
      int relPos=GetRelationOfPointToArea(node,b);
 
      ++count;
 
      if (relPos>0) {
        ++pro;
      }
      else if (relPos<0) {
        ++contra;
      }
 
      if (count>=100 && double(pro)/20.0>double(contra)) {
        return true;
      }
    }
 
    return double(pro)/20.0>double(contra);
  }

  template<typename N,typename M>
  inline bool IsAreaSubOfAreaOrSame(const std::vector<N>& a,
                                    const std::vector<M>& b)
  {
    size_t pro=0;
    size_t contra=0;
    size_t count=0;
 
    pro=0;
    contra=0;
    count=0;
 
    for (const auto& node : a) {
      int relPos=GetRelationOfPointToArea(node,b);
 
      ++count;
 
      if (relPos>0) {
        ++pro;
      }
      else if (relPos<0) {
        ++contra;
      }
 
      if (count>=100 && pro/20>contra) {
        return true;
      }
    }
 
    if (pro == 0 && contra == 0 && count > 0) {
      return true;
    }
 
    return pro/20.0>contra;
  }

  template<typename N>
  bool AreaIsSimple(std::vector<N> points)
  {
    if (points.size()<3) {
      return false;
    }
 
    points.push_back(points[0]);
 
    size_t edgesIntersect=0;
 
    for (size_t i=0; i<points.size()-1; i++) {
      edgesIntersect=0;
 
      for (size_t j=i+1; j<points.size()-1; j++) {
        if (LinesIntersect(points[i],
                           points[i+1],
                           points[j],
                           points[j+1])) {
          edgesIntersect++;
 
          if (i==0) {
            if (edgesIntersect>2) {
              return false;
            }
          }
          else {
            if (edgesIntersect>1) {
              return false;
            }
          }
        }
      }
    }
 
    return true;
  }

  template<typename N>
  bool AreaIsSimple(const std::vector<std::pair<N,N> >& edges,
                    const std::vector<bool>& edgeStartsNewPoly)
  {
    size_t edgesIntersect=0;
 
    for (size_t i=0; i<edges.size(); i++) {
      edgesIntersect=0;
 
      for (size_t j=i+1; j<edges.size(); j++) {
        if (LinesIntersect(edges[i].first,
                           edges[i].second,
                           edges[j].first,
                           edges[j].second)) {
          edgesIntersect++;
 
          // we check the first edge of a sole polygon against every
          // other edge and expect to see 2 intersections for
          // adjacent edges; polygon is complex if there are more
          // intersections
          if (edgeStartsNewPoly[i]) {
            if (edgesIntersect>2) {
              return false;
            }
          }
 
          // otherwise we check an edge that isn't the first
          // edge against every other edge excluding those that
          // have already been tested (this means one adjacent
          // edge); polygon is complex if there is more than one
          // intersection
          else {
            if (edgesIntersect>1) {
              return false;
            }
          }
        }
      }
    }
 
    return true;
  }

  template<typename N>
  bool AreaIsCCW(const std::vector<N> &edges)
  {
    // based on http://en.wikipedia.org/wiki/Curve_orientation
    // and http://local.wasp.uwa.edu.au/~pbourke/geometry/clockwise/
 
    // note: polygon must be simple
    // note: this assumes 2d cartesian coordinate space is used;
    // for geographic, expect Vec2.x=lon and Vec2.y=lat!
 
    size_t ptIdx=0;
 
    for (size_t i=1; i<edges.size(); i++) {
      // find the point with the smallest y value,
      if (edges[i].y<edges[ptIdx].y) {
        ptIdx=i;
      }
      // if y values are equal save the point with greatest x
      else if (edges[i].y==edges[ptIdx].y) {
        if (edges[i].x<edges[ptIdx].x) {
          ptIdx=i;
        }
      }
    }
 
    size_t prevIdx=(ptIdx==0) ? edges.size()-1 : ptIdx-1;
    size_t nextIdx=(ptIdx==edges.size()-1) ? 0 : ptIdx+1;
 
    double signedArea=(edges[ptIdx].x-edges[prevIdx].x)*
                      (edges[nextIdx].y-edges[ptIdx].y)-
                      (edges[ptIdx].y-edges[prevIdx].y)*
                      (edges[nextIdx].x-edges[ptIdx].x);
 
    return signedArea>0.0;
  }

  template<typename N>
  bool AreaIsValid(std::vector<N>& outerPoints,
                   std::vector<std::vector<N> >& innerPoints)
  {
    if (outerPoints.size()<3) {
      return false;
    }
 
    size_t numEdges=outerPoints.size();
 
    for (size_t i=0; i<innerPoints.size(); i++) {
      numEdges+=innerPoints[i].size();
    }
 
    std::vector<bool>            edgeStartsNewPoly(numEdges,false);
    std::vector<std::pair<N,N> > listEdges(numEdges);
    size_t                       cEdge=0;
 
    // temporarily wrap around vertices
    // (first == last) to generate edge lists
    outerPoints.push_back(outerPoints[0]);
    for (size_t i=0; i<innerPoints.size(); i++) {
      innerPoints[i].push_back(innerPoints[i][0]);
    }
 
    // outer poly
    edgeStartsNewPoly[0]=true;
    for (size_t i=1; i<outerPoints.size(); i++) {
      std::pair<N, N> outerEdge;
 
      outerEdge.first=outerPoints[i-1];
      outerEdge.second=outerPoints[i];
      listEdges[cEdge]=outerEdge;
 
      cEdge++;
    }
 
    // inner polys
    for (size_t i=0; i<innerPoints.size(); i++) {
      edgeStartsNewPoly[cEdge]=true;
 
      for (size_t j=1; j<innerPoints[i].size(); j++) {
        std::pair<N, N> innerEdge;
 
        innerEdge.first=innerPoints[i][j-1];
        innerEdge.second=innerPoints[i][j];
        listEdges[cEdge]=innerEdge;
 
        cEdge++;
      }
    }
 
    // revert vertex list modifications (not
    // really the 'nicest' way of doing this)
    outerPoints.pop_back();
    for (size_t i=0; i<innerPoints.size(); i++) {
      innerPoints[i].pop_back();
    }
 
    if (AreaIsSimple(listEdges,edgeStartsNewPoly))
        {
      // expect listOuterPts to be CCW and innerPts
      // to be CW, if not then reverse point order
 
      if (!AreaIsCCW(outerPoints)) {
        std::reverse(outerPoints.begin(),
                     outerPoints.end());
      }
 
      for (int i=0; i<innerPoints.size(); i++) {
        if (AreaIsCCW(innerPoints[i])) {
          std::reverse(innerPoints[i].begin(),
                       innerPoints[i].end());
        }
      }
    }
    else {
      return false;
    }
 
    return true;
  }

  inline void Normalize(double x,
                        double y,
                        double& nx,
                        double& ny)
  {
    double length=sqrt(x*x+y*y);
 
    if (length!=0) {
      nx=x/length;
      ny=y/length;
    }
  }

  inline double Det(double x1,
                    double y1,
                    double x2,
                    double y2)
  {
    return x1*y2-y1*x2;
  }
 

  extern OSMSCOUT_API size_t Pow(size_t a, size_t b);

  template<typename N>
  bool AreaIsClockwise(const std::vector<N>& edges)
  {
    assert(edges.size()>=3);
    // based on http://en.wikipedia.org/wiki/Curve_orientation
    // and http://local.wasp.uwa.edu.au/~pbourke/geometry/clockwise/
 
    // note: polygon must be simple
 
    size_t ptIdx=0;
 
    for (size_t i=1; i<edges.size(); i++) {
      // find the point with the smallest y value,
      if (edges[i].GetLat()<edges[ptIdx].GetLat()) {
        ptIdx=i;
      }
        // if y values are equal save the point with greatest x
      else if (edges[i].GetLat()==edges[ptIdx].GetLat()) {
        if (edges[i].GetLon()<edges[ptIdx].GetLon()) {
          ptIdx=i;
        }
      }
    }
 
    size_t prevIdx=(ptIdx==0) ? edges.size()-1 : ptIdx-1;
    size_t nextIdx=(ptIdx==edges.size()-1) ? 0 : ptIdx+1;
 
    double signedArea=(edges[ptIdx].GetLon()-edges[prevIdx].GetLon())*
                      (edges[nextIdx].GetLat()-edges[ptIdx].GetLat())-
                      (edges[ptIdx].GetLat()-edges[prevIdx].GetLat())*
                      (edges[nextIdx].GetLon()-edges[ptIdx].GetLon());
 
    return signedArea<0.0;
  }
 

  template<typename N>
  double CalculateDistancePointToLineSegment(const N& p,
                                             const N& a,
                                             const N& b)
  {
    double xdelta=b.GetLon()-a.GetLon();
    double ydelta=b.GetLat()-a.GetLat();
 
    if (xdelta==0 && ydelta==0) {
      return std::numeric_limits<double>::infinity();
    }
 
    double u=((p.GetLon()-a.GetLon())*xdelta+(p.GetLat()-a.GetLat())*ydelta)/(xdelta*xdelta+ydelta*ydelta);
 
    double cx;
    double cy;
 
    if (u<0) {
      cx=a.GetLon();
      cy=a.GetLat();
    }
    else if (u>1) {
      cx=b.GetLon();
      cy=b.GetLat();
    }
    else {
      cx=a.GetLon()+u*xdelta;
      cy=a.GetLat()+u*ydelta;
    }
 
    double dx=cx-p.GetLon();
    double dy=cy-p.GetLat();
 
    return sqrt(dx*dx+dy*dy);
  }
 
  extern OSMSCOUT_API double CalculateDistancePointToLineSegment(const GeoCoord& p,
                                                                 const GeoCoord& a,
                                                                 const GeoCoord& b,
                                                                 GeoCoord& intersection);

  extern OSMSCOUT_API Distance GetSphericalDistance(const GeoCoord& a,
                                                    const GeoCoord& b);

  extern OSMSCOUT_API Distance GetEllipsoidalDistance(double aLon, double aLat,
                                                      double bLon, double bLat);

  extern OSMSCOUT_API Distance GetEllipsoidalDistance(const GeoCoord& a,
                                                      const GeoCoord& b);

  extern OSMSCOUT_API void GetEllipsoidalDistance(double lat1, double lon1,
                                                  const Bearing &bearing,
                                                  const Distance &distance,
                                                  double& lat2, double& lon2);

  extern OSMSCOUT_API GeoCoord GetEllipsoidalDistance(const GeoCoord& position,
                                                      const Bearing &bearing,
                                                      const Distance &distance);

  extern OSMSCOUT_API Bearing GetSphericalBearingInitial(const GeoCoord& a,
                                                         const GeoCoord& b);

  extern OSMSCOUT_API Bearing GetSphericalBearingFinal(const GeoCoord& a,
                                                       const GeoCoord& b);

  extern OSMSCOUT_API double GetDistanceInLonDegrees(const Distance &d,
                                                     double latitude=0);

  extern OSMSCOUT_API double NormalizeRelativeAngle(double angle);
 
  struct OSMSCOUT_API ScanCell
  {
    int x;
    int y;
 
    ScanCell(int x, int y);
 
    void Set(int x, int y)
    {
      this->x=x;
      this->y=y;
    }
 
    bool operator==(const ScanCell& other) const
    {
      return x==other.x && y==other.y;
    }
 
    bool operator<(const ScanCell& other) const
    {
      return std::tie(y, x) < std::tie(other.y, other.x);
    }
 
    bool IsEqual(const ScanCell& other) const
    {
      return x==other.x && y==other.y;
    }
 
    int GetX() const
    {
      return x;
    }
 
    int GetY() const
    {
      return y;
    }
  };
 
  struct OSMSCOUT_API SegmentGeoBox
  {
    size_t from;
    size_t to; 
    GeoBox bbox;
  };

  void OSMSCOUT_API ScanConvertLine(int x1, int y1,
                                    int x2, int y2,
                                    std::vector<ScanCell>& cells);

  extern OSMSCOUT_API double DistanceToSegment(double px, double py,
                                               double p1x, double p1y,
                                               double p2x, double p2y,
                                               double& r,
                                               double& qx, double& qy);
 
  extern OSMSCOUT_API double DistanceToSegment(const GeoCoord& point,
                                               const GeoCoord& segmentStart,
                                               const GeoCoord& segmentEnd,
                                               double& r,
                                               GeoCoord& intersection);
 
  extern OSMSCOUT_API double DistanceToSegment(const std::vector<Point>& points,
                                               const GeoCoord& segmentStart,
                                               const GeoCoord& segmentEnd,
                                               GeoCoord& location,
                                               GeoCoord& intersection);
 
  extern OSMSCOUT_API double DistanceToSegment(const GeoBox& boundingBox,
                                               const GeoCoord& segmentStart,
                                               const GeoCoord& segmentEnd,
                                               GeoCoord& location,
                                               GeoCoord& intersection);

  struct OSMSCOUT_API PathIntersection
  {
    GeoCoord point;           
    size_t   aIndex;          
    size_t   bIndex;          
    double   orientation;     
    double   aDistanceSquare; 
    double   bDistanceSquare; 
  };

  template<typename N>
  void GetSegmentBoundingBox(const std::vector<N>& path,
                             size_t from, size_t to,
                             GeoBox& boundingBox)
  {
    if (path.empty() || from>=to){
      boundingBox.Invalidate();
      return;
    }
 
    double minLon=path[from%path.size()].GetLon();
    double maxLon=minLon;
    double minLat=path[from%path.size()].GetLat();
    double maxLat=minLat;
 
    for (size_t i=from; i<to; i++) {
      minLon=std::min(minLon,path[i%path.size()].GetLon());
      maxLon=std::max(maxLon,path[i%path.size()].GetLon());
      minLat=std::min(minLat,path[i%path.size()].GetLat());
      maxLat=std::max(maxLat,path[i%path.size()].GetLat());
    }
 
    boundingBox.Set(GeoCoord(minLat,minLon),
                    GeoCoord(maxLat,maxLon));
  }

  template<typename N>
  void ComputeSegmentBoxes(const std::vector<N>& path,
                           std::vector<SegmentGeoBox> &segmentBoxes,
                           size_t bound,
                           size_t segmentSize = 1000)
  {
    assert(segmentSize>0);
    for (size_t i=0;i<bound;i+=segmentSize){
      SegmentGeoBox box;
      box.from = i;
      box.to = std::min(i+segmentSize,bound);
      GetSegmentBoundingBox(path, box.from, box.to, box.bbox);
      segmentBoxes.push_back(box);
    }
  }

  template<typename N>
  void FindPathIntersections(const std::vector<N> &aPath,
                             const std::vector<N> &bPath,
                             bool aClosed,
                             bool bClosed,
                             std::vector<PathIntersection> &intersections,
                             size_t aStartIndex=0,
                             size_t bStartIndex=0)
  {
    size_t aIndex=aStartIndex;
    size_t bIndex=bStartIndex;
    // bounds are inclusive
    size_t aBound=aClosed?aIndex+aPath.size():aPath.size()-1;
    size_t bBound=bClosed?bIndex+bPath.size():bPath.size()-1;
    size_t bStart=bIndex;
 
    if (bStart>=bBound || aIndex>=aBound){
      return;
    }
 
    GeoBox bBox;
    GetSegmentBoundingBox(bPath,bIndex,bBound+1,bBox);
    GeoBox aLineBox;
    GeoBox bLineBox;
 
    // compute b-boxes for B path, each 1000 point-long segment
    std::vector<SegmentGeoBox> bSegmentBoxes;
    ComputeSegmentBoxes(bPath,bSegmentBoxes,bBound+1, 1000);
 
    for (;aIndex<aBound;aIndex++){
      N a1=aPath[aIndex%aPath.size()];
      N a2=aPath[(aIndex+1)%aPath.size()];
      aLineBox.Set(GeoCoord(a1.GetLat(),a1.GetLon()),
                   GeoCoord(a2.GetLat(),a2.GetLon()));
      if (!bBox.Intersects(aLineBox,/*openInterval*/false)){
        continue;
      }
      for (bIndex=bStart;bIndex<bBound;bIndex++){
        N b1=bPath[bIndex%bPath.size()];
        N b2=bPath[(bIndex+1)%bPath.size()];
        bLineBox.Set(GeoCoord(b1.GetLat(),b1.GetLon()),
                     GeoCoord(b2.GetLat(),b2.GetLon()));
 
        if (!bSegmentBoxes[bIndex/1000].bbox.Intersects(aLineBox,/*openInterval*/false) &&
            !aLineBox.Intersects(bLineBox,/*openInterval*/false)){
          // round up
          bIndex+=std::max(0, 998-(int)(bIndex%1000));
          continue;
        }
        if (!aLineBox.Intersects(bLineBox,/*openInterval*/false)){
          continue;
        }
        PathIntersection pi;
        if (GetLineIntersection(a1,a2,
                                b1,b2,
                                pi.point)){
 
          // if b2==pi or a1==pi, orientation is zero then
          // for that reason, we prolong both lines for computation
          // of intersection orientation
          GeoCoord pointBefore(a1.GetLat()-(a2.GetLat()-a1.GetLat()),
                               a1.GetLon()-(a2.GetLon()-a1.GetLon()));
          GeoCoord pointAfter(b2.GetLat()+(b2.GetLat()-b1.GetLat()),
                              b2.GetLon()+(b2.GetLon()-b1.GetLon()));
          pi.orientation=(pi.point.GetLon()-pointBefore.GetLon())*
                         (pointAfter.GetLat()-pi.point.GetLat())-
                         (pi.point.GetLat()-pointBefore.GetLat())*
                         (pointAfter.GetLon()-pi.point.GetLon());
 
          pi.aIndex=aIndex;
          pi.bIndex=bIndex;
          pi.aDistanceSquare=DistanceSquare(GeoCoord(a1.GetLat(),a1.GetLon()),pi.point);
          pi.bDistanceSquare=DistanceSquare(GeoCoord(b1.GetLat(),b1.GetLon()),pi.point);
 
          intersections.push_back(pi);
        }
      }
    }
  }

  template<typename N>
  bool FindIntersection(const std::vector<N> &way,
                        size_t &i,
                        size_t &j)
  {
    GeoBox lineBox;
 
    // compute b-boxes for path, each 1000 point-long segment
    std::vector<SegmentGeoBox> segmentBoxes;
    ComputeSegmentBoxes(way,segmentBoxes,way.size(),1000);
 
    for (; i<way.size()-1; i++) {
      N i1=way[i];
      N i2=way[i+1];
      lineBox.Set(GeoCoord(i1.GetLat(),i1.GetLon()),
                  GeoCoord(i2.GetLat(),i2.GetLon()));
 
      for (j=i+1; j<way.size()-1; j++) {
        if (!segmentBoxes[j/1000].bbox.Intersects(lineBox,/*openInterval*/false) &&
            !segmentBoxes[(j+1)/1000].bbox.Intersects(lineBox,/*openInterval*/false)){
          // round up
          j+=std::max(0, 998-(int)(j%1000));
          continue;
        }
        if (LinesIntersect(way[i],
                           way[i+1],
                           way[j],
                           way[j+1])) {
 
          if (way[i].IsEqual(way[j]) ||
              way[i].IsEqual(way[j+1]) ||
              way[i+1].IsEqual(way[j]) ||
              way[i+1].IsEqual(way[j+1])) {
            continue; // ignore sibling edges
          }
 
          return true;
        }
      }
    }
    return false;
  }
 
  class OSMSCOUT_API PolygonMerger
  {
  private:
    struct Edge
    {
      size_t fromIndex;
      size_t toIndex;
    };
 
  public:
    struct Polygon
    {
      std::vector<Point> coords;
    };
 
  private:
    std::unordered_map<Id,size_t>                               nodeIdIndexMap;
    std::vector<Point>                                          nodes;
    std::list<Edge>                                             edges;
 
  private:
    void RemoveEliminatingEdges();
 
  public:
    void AddPolygon(const std::vector<Point>& polygonsCoords);
 
    bool Merge(std::list<Polygon>& result);
  };
 
  struct OSMSCOUT_API CellDimension
  {
    double width;
    double height;
  };
 
  const size_t CELL_DIMENSION_MAX   = 25;
  const size_t CELL_DIMENSION_COUNT = CELL_DIMENSION_MAX+1;
 
  extern OSMSCOUT_API const std::array<CellDimension,CELL_DIMENSION_COUNT> cellDimension;

  class OSMSCOUT_API GeoBoxPartitioner
  {
  public:
    enum class Direction
    {
      HORIZONTAL,
      VERTICAL
    };
 
  private:
    GeoBox    box;
    Direction direction;
    double    parts;
    size_t    currentIndex=0;
    GeoBox    currentBox;
 
  private:
    void CalculateBox();
 
  public:
    GeoBoxPartitioner(const GeoBox& box,
                     Direction direction,
                     size_t parts)
     : box(box),
       direction(direction),
       parts((double)parts)
    {
      assert(currentIndex<parts);
      CalculateBox();
    }
 
    void Advance()
    {
      currentIndex++;
 
      if (currentIndex<parts) {
        CalculateBox();
      }
    }
 
    GeoBox GetCurrentGeoBox() const
    {
      assert(currentIndex<parts);
      return currentBox;
    }
 
    size_t GetCurrentIndex() const
    {
      assert(currentIndex<parts);
      return currentIndex;
    }
  };
}
 
namespace std {
  template <>
  struct hash<osmscout::ScanCell>
  {
    size_t operator()(const osmscout::ScanCell& cell) const
    {
      return hash<int>{}(cell.GetX()) ^ (hash<int>{}(cell.GetY()) << 1);
    }
  };
}
 
#endif