src/terralib/geometry/Utils.cpp

Go to the documentation of this file.

/*  Copyright (C) 2008 National Institute For Space Research (INPE) - Brazil.

    This file is part of the TerraLib - a Framework for building GIS enabled applications.

    TerraLib 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 3 of the License,
    or (at your option) any later version.

    TerraLib 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 TerraLib. See COPYING. If not, write to
    TerraLib Team at <terralib-team@terralib.org>.
 */

/*!
  \file terralib/geometry/Utils.cpp

  \brief Utility functions for the Geometry Module.
*/

// TerraLib
#include "../common/STLUtils.h"
#include "../core/translator/Translator.h"
#include "Curve.h"
#include "Envelope.h"
#include "Exception.h"
#include "Geometry.h"
#include "GeometryCollection.h"
#include "GEOSGeometryFactory.h"
#include "GEOSReader.h"
#include "GEOSWriter.h"
#include "Line.h"
#include "LinearRing.h"
#include "LineString.h"
#include "MultiLineString.h"
#include "MultiPolygon.h"
#include "Point.h"
#include "Polygon.h"
#include "Utils.h"

// STL
#include <cmath>
#include <memory>

#include <geos/algorithm/CentroidArea.h>
#include <geos/algorithm/CGAlgorithms.h>
#include <geos/algorithm/LineIntersector.h>
#include <geos/geom/Coordinate.h>
#include <geos/geom/Geometry.h>
#include <geos/geom/LineSegment.h>
#include <geos/geom/Point.h>
#include <geos/geom/Polygon.h>
#include <geos/operation/polygonize/Polygonizer.h>
#include <geos/operation/union/CascadedPolygonUnion.h>
#include <geos/operation/union/UnaryUnionOp.h>
#include <geos/operation/distance/DistanceOp.h>
#include <geos/operation/overlay/snap/GeometrySnapper.h>
#include <geos/operation/valid/IsValidOp.h>
#include <geos/util/GEOSException.h>

const double TOLERANCE = 0.000000001;

bool te::gm::IsEqual(const te::gm::Point& p1, const te::gm::Point& p2)
{
  if ((p1.getX() == p2.getX()) && (p1.getY() == p2.getY()))
  {
    return true;
  }

  return false;
}

bool te::gm::IsEqual(const te::gm::Point& p1, const te::gm::Point& p2, double tol)
{
  if (std::fabs(p1.getX() - p2.getX()) > tol)
  {
    return false;
  }
  if (std::fabs(p1.getY() - p2.getY()) > tol)
  {
    return false;
  }

  return true;
}

te::gm::Envelope te::gm::CreateEnvelope(double x1, double y1, double x2, double y2)
{
  te::gm::Envelope envelope(std::min(x1, x2)
    , std::min(y1, y2)
    , std::max(x1, x2)
    , std::max(y1, y2));

  return envelope;
}

te::gm::Envelope te::gm::CreateEnvelope(const te::gm::Point& p1, const te::gm::Point& p2)
{
  te::gm::Envelope envelope(std::min(p1.getX(), p2.getX())
    , std::min(p1.getY(), p2.getY())
    , std::max(p1.getX(), p2.getX())
    , std::max(p1.getY(), p2.getY()));

  return envelope;
}

std::unique_ptr<te::gm::Geometry> te::gm::GetGeometryUnion(const std::vector<gm::Geometry*> &geomVec)
{
  std::unique_ptr<te::gm::Geometry> geometry;

  if (geomVec.size() == 1)
  {
    te::dt::AbstractData* abs = geomVec[0]->clone();
    geometry.reset(static_cast<te::gm::Geometry*>(abs));
  }
  else if (geomVec.size() > 1)
  {
    //First clone the vector
    te::sam::rtree::Index<std::size_t, 8> rtree;
    
    std::vector<te::gm::Geometry*> geomVecCopy;
    geomVecCopy.reserve(geomVec.size());

    
    for (std::size_t i = 0; i < geomVec.size(); ++i)
    {
      te::gm::Geometry* geometry = static_cast<te::gm::Geometry*>(geomVec[i]->clone());

      rtree.insert(*geometry->getMBR(), geomVecCopy.size());
      geomVecCopy.push_back(geometry);
    }

    for (std::size_t i = 0; i < geomVecCopy.size(); ++i)
    {
      const te::gm::Envelope* env_A = geomVecCopy[i]->getMBR();

      te::gm::Envelope expandedEnv_A(env_A->m_llx - TOLERANCE, env_A->m_lly - TOLERANCE, env_A->m_urx + TOLERANCE, env_A->m_ury + TOLERANCE);

      std::vector<std::size_t> vecCandidatesIndexes;
      rtree.search(expandedEnv_A, vecCandidatesIndexes);

      if (vecCandidatesIndexes.empty())
      {
        continue;
      }

      //we copy the geom vector so we can remove the current geometry from it
      std::vector<te::gm::Geometry*> geomVecCopyFiltered;
      for (std::size_t j = 0; j < vecCandidatesIndexes.size(); ++j)
      {
        std::size_t index = vecCandidatesIndexes[j];
        if(i == index)
        {
          continue;
        }

        te::gm::Geometry* geometry = geomVecCopy[index];
        geomVecCopyFiltered.push_back(geometry);
      }

      bool wasChanged = false;
      std::unique_ptr<te::gm::Geometry> geometryResult(te::gm::PrepareGeometriesToIntersection(geomVecCopy[i], geomVecCopyFiltered, wasChanged));

     if (wasChanged == true)
      {
        //then we replace the resulted geometry into our base geometry vector
        delete geomVecCopy[i];
        geomVecCopy[i] = geometryResult.release();
      }
    }

    std::unique_ptr<te::gm::GeometryCollection> gc(new te::gm::GeometryCollection(0, te::gm::GeometryCollectionType, geomVec[0]->getSRID()));

    for (std::size_t i = 0; i < geomVecCopy.size(); ++i)
      gc->add(geomVecCopy[i]);

    bool success = false;

    try
    {
      geometry.reset(te::gm::UnaryUnion(gc.get()));
      success = true;
    }
    catch(...)
    {
    }

    if(success == false)
    {
      try
      {
        std::unique_ptr<te::gm::Geometry> bufferedGeometry(gc->buffer(0));
        geometry.reset(te::gm::UnaryUnion(bufferedGeometry.get()));
      }
      catch(geos::util::GEOSException& e)
      {
        throw Exception(e.what());
      }
    }
  }

  return geometry;
}

te::gm::Geometry* te::gm::GetGeomFromEnvelope(const Envelope* const e, int srid)
{
// create an outer ring with the same envelope as our envelope
  LinearRing* r = new LinearRing(5, LineStringType, srid, new Envelope(*e));

  r->setPoint(0, e->m_llx, e->m_lly);
  r->setPoint(1, e->m_urx, e->m_lly);
  r->setPoint(2, e->m_urx, e->m_ury);
  r->setPoint(3, e->m_llx, e->m_ury);
  r->setPoint(4, e->m_llx, e->m_lly);

// create the polygon
  Polygon* p = new Polygon(1, PolygonType, srid, new Envelope(*e));
  p->setRingN(0, r);

  return p;
}

bool te::gm::SatisfySpatialRelation(const Geometry* g1,
                                    const Geometry* g2,
                                    SpatialRelation relation)
{
  switch(relation)
  {
    case CONTAINS:
      return g1->contains(g2);

    case COVEREDBY:
      return g1->coveredBy(g2);

    case COVERS:
      return g1->covers(g2);

    case CROSSES:
      return g1->crosses(g2);

    case DISJOINT:
      return g1->disjoint(g2);

    case EQUALS:
      return g1->equals(g2);

    case INTERSECTS:
      return g1->intersects(g2);

    case OVERLAPS:
      return g1->overlaps(g2);

    case TOUCHES:
      return g1->touches(g2);

    case WITHIN:
      return g1->within(g2);

    default:
      throw Exception(TE_TR("Invalid spatial relation!"));
  }
}

te::gm::Envelope te::gm::AdjustToCut(const Envelope & env, double bWidth, double bHeight)
{
  double auxD;
  int auxI;

  int magicX;
  auxD = env.m_llx/bWidth;
  auxI = static_cast<int>(env.m_llx/bWidth);
  if (env.m_llx < 0 && (auxD - auxI) != 0)
    magicX = static_cast<int>(env.m_llx/bWidth - 1);
  else
    magicX = auxI;

  int magicY;
  auxD = env.m_lly/bHeight;
  auxI = static_cast<int>(env.m_lly/bHeight);
  if (env.m_lly < 0 && (auxD - auxI) != 0)
    magicY  = static_cast<int>(env.m_lly/bHeight - 1);
  else
    magicY  = auxI;

  double xi = magicX*bWidth;
  double yi = magicY*bHeight;

  int magicX2;
  auxD = env.m_urx/bWidth;
  auxI = static_cast<int>(env.m_urx/bWidth);
  if ((env.m_urx < 0) || (auxD - auxI) == 0)
    magicX2 = static_cast<int>(env.m_urx/bWidth);
  else
    magicX2 = static_cast<int>(env.m_urx/bWidth + 1);

  int magicY2;
  auxD = env.m_ury/bHeight;
  auxI = static_cast<int>(env.m_ury/bHeight);
  if ((env.m_ury < 0) || (auxD - auxI) == 0)
    magicY2 = static_cast<int>(env.m_ury/bHeight);
  else
    magicY2 = static_cast<int>(env.m_ury/bHeight + 1);

  double xf = (magicX2)*bWidth;
  double yf = (magicY2)*bHeight;

  return te::gm::Envelope(xi,yi,xf,yf);
}

template<class T1, class T2> bool te::gm::Intersects(const T1& o1, const T2& o2)
{
  return o1->intersects(o2);
}

template<> bool te::gm::Intersects(const te::gm::Point& point, const te::gm::Envelope& e)
{
  // point to the right of envelope
  if(point.getX() > e.m_urx)
    return false;

  // point to the left of envelope
  if(e.m_llx > point.getX())
    return false;

  // point is above envelope
  if(point.getY() > e.m_ury)
    return false;

  // point is below envelope
  if(e.m_lly > point.getY())
    return false;

  return true;
}

te::gm::Coord2D* te::gm::locateAlong(const LineString* line, double initial, double final, double target)
{
  double tTof; // Distance of target to fist point

  std::map<int, double> pointLenghtFromFirst;

  double fullLineLenght = 0;

  pointLenghtFromFirst[0] = 0;
  for(std::size_t i = 1; i < line->getNPoints(); i++)
  {
    std::unique_ptr<Point> pA = line->getPointN(i);
    std::unique_ptr<Point> pB = line->getPointN(i - 1);

    fullLineLenght += pA->distance(pB.get());
    pointLenghtFromFirst[static_cast<int>(i)] = fullLineLenght;
  }

  double diference = final-initial;

  tTof = ((target-initial)*fullLineLenght)/diference;

  int onTheFly = -1;

  int pointBeforeTarget = -1;

  for(std::size_t i = 0; i < pointLenghtFromFirst.size(); i++)
  {
    if(pointLenghtFromFirst[static_cast<int>(i)] == tTof)
    {
      onTheFly = static_cast<int>(i);
      break;
    }
    if(tTof >= pointLenghtFromFirst[static_cast<int>(i)] && tTof < pointLenghtFromFirst[static_cast<int>(i+1)])
    {
      pointBeforeTarget = static_cast<int>(i);
      break;
    }
  }

  te::gm::Coord2D* targetCoord = new te::gm::Coord2D();

  if(onTheFly >= 0)
  {
    targetCoord->x = line->getPointN(static_cast<std::size_t>(onTheFly))->getX();
    targetCoord->y = line->getPointN(static_cast<std::size_t>(onTheFly))->getY();

    return targetCoord;
  }

  std::unique_ptr<Point> p1 = line->getPointN(static_cast<std::size_t>(pointBeforeTarget));
  std::unique_ptr<Point> p2 = line->getPointN(static_cast<std::size_t>(pointBeforeTarget + 1));

  // Partial distance: target to p1
  double pd = tTof - pointLenghtFromFirst[pointBeforeTarget];

  // Total distance: p2 to p1
  double td = pointLenghtFromFirst[pointBeforeTarget+1] - pointLenghtFromFirst[pointBeforeTarget];

  targetCoord->x = ((pd*(p2->getX()-p1->getX()))/td)+p1->getX();

  targetCoord->y = ((pd*(p2->getY()-p1->getY()))/td)+p1->getY();

  return targetCoord;
}

void te::gm::Multi2Single(const te::gm::Geometry* g, std::vector<te::gm::Geometry*>& geoms)
{
  const te::gm::GeometryCollection* gc = dynamic_cast<const te::gm::GeometryCollection*>(g);
  if (gc)
  {
    for (std::size_t i = 0; i < gc->getNumGeometries(); ++i)
    {
      te::gm::Geometry* currentGeom = gc->getGeometryN(i);
      currentGeom->setSRID(gc->getSRID());
      Multi2Single(currentGeom, geoms);
    }
  }
  else
    geoms.push_back(dynamic_cast<te::gm::Geometry*>(g->clone()));
}

te::gm::Geometry* te::gm::CascadedPolygonUnion(const std::vector<Polygon*>& polygonVector)
{
#ifdef TERRALIB_GEOS_ENABLED

  std::vector<geos::geom::Polygon*> geosPolygonVector;

  try
  {
    if(polygonVector.empty())
      return nullptr;

    int srid = polygonVector[0]->getSRID();

    for(auto i : polygonVector)
    {
      if(i->getGeomTypeId() == te::gm::PolygonType)
      {
        geos::geom::Polygon* geosPolygon = GEOSWriter::write(i);
        geosPolygonVector.push_back(geosPolygon);
      }
    }

    std::unique_ptr<geos::geom::Geometry> unionGeom(geos::operation::geounion::CascadedPolygonUnion::Union(&geosPolygonVector));

    te::common::FreeContents(geosPolygonVector);

    if(unionGeom == nullptr)
      return nullptr;

    unionGeom->setSRID(srid);

    return GEOSReader::read(unionGeom.get());
  }
  catch(geos::util::GEOSException& e)
  {
    te::common::FreeContents(geosPolygonVector);
    throw Exception(e.what());
  }

#else
  throw te::common::Exception(TE_TR("CascadedPolygonUnion routine is supported by GEOS! Please, enable the GEOS support."));
#endif
}

te::gm::Geometry* te::gm::UnaryUnion(te::gm::Geometry* geom)
{
#ifdef TERRALIB_GEOS_ENABLED

  try
  {
    std::unique_ptr<geos::geom::Geometry> geomGeos(GEOSWriter::write(geom));

    std::unique_ptr<geos::geom::Geometry> unionGeom = geos::operation::geounion::UnaryUnionOp::Union(*(geomGeos.get()));

    unionGeom->setSRID(geom->getSRID());

    return GEOSReader::read(unionGeom.get());
  }
  catch(geos::util::GEOSException& e)
  {
    throw Exception(e.what());
  }

#else
  throw te::common::Exception(TE_TR("Union routine is supported by GEOS! Please, enable the GEOS support."));
#endif
}

void te::gm::ClosestPoints(te::gm::Geometry* geomA, te::gm::Geometry* geomB,
  te::gm::Coord2D& coordA, te::gm::Coord2D& coordB)
{
#ifdef TERRALIB_GEOS_ENABLED
  std::unique_ptr<geos::geom::Geometry> thisGeomA(GEOSWriter::write(geomA));
  std::unique_ptr<geos::geom::Geometry> thisGeomB(GEOSWriter::write(geomB));

  geos::operation::distance::DistanceOp op(thisGeomA.get(), thisGeomB.get());

  geos::geom::CoordinateSequence* cs = op.closestPoints();

  if (cs->getSize() == 2)
  {
    coordA = te::gm::Coord2D(cs->getAt(0).x, cs->getAt(0).y);
    coordB = te::gm::Coord2D(cs->getAt(1).x, cs->getAt(1).y);
  }

  delete cs;

#else
  throw te::common::Exception(TE_TR("getCentroid routine is supported by GEOS! Please, enable the GEOS support."));
#endif
}

te::gm::Line* te::gm::GetIntersectionLine(te::gm::Geometry* geom, te::gm::Coord2D coord)
{
  std::vector<te::gm::Geometry*> geomVec;
  std::unique_ptr<te::gm::Point> pRef(new te::gm::Point(coord.x, coord.y, geom->getSRID()));

  te::gm::Multi2Single(geom, geomVec);

  for (std::size_t t = 0; t < geomVec.size(); ++t)
  {
    te::gm::Geometry* g = geomVec[t];

    if (g->getGeomTypeId() == te::gm::LineStringType)
    {
      te::gm::LineString* ls = static_cast<te::gm::LineString*>(g);

      if (ls->getNPoints() == 2)
      {
        std::unique_ptr<te::gm::Point> p0 = ls->getPointN(0);
        std::unique_ptr<te::gm::Point> p1 = ls->getPointN(1);

        te::gm::Line* line = new te::gm::Line(te::gm::LineStringType, geom->getSRID());

        line->setCoord(0, p0->getX(), p0->getY());
        line->setCoord(1, p1->getX(), p1->getY());

        return line;
      }
      else
      {
        for (std::size_t p = 0; p < ls->getNPoints() - 1; ++p)
        {
          std::unique_ptr<te::gm::Point> p0 = ls->getPointN(p);
          std::unique_ptr<te::gm::Point> p1 = ls->getPointN(p + 1);

          std::unique_ptr<te::gm::Line> line(new te::gm::Line(te::gm::LineStringType, geom->getSRID()));
          line->setCoord(0, p0->getX(), p0->getY());
          line->setCoord(1, p1->getX(), p1->getY());

          std::unique_ptr<te::gm::Geometry> geomMBR(te::gm::GetGeomFromEnvelope(line->getMBR(), geom->getSRID()));

          if (geomMBR->intersects(pRef.get()))
          {
            return line.release();
          }
        }
      }
    }
  }

  return nullptr;
}

double te::gm::GetAngle(te::gm::Coord2D coordA, te::gm::Coord2D coordB)
{
  if (coordA.x == coordB.x)
  {
    return 1.57079632679489661923;
  }

  if (coordA.y == coordB.y)
  {
    return 0.0;
  }

  double dx = (coordA.x - coordB.x);
  double dy = (coordA.y - coordB.y);

  double rad = std::atan(dy / dx);

  //Convert radians to degrees
  double angle = rad / 0.01745329251994329576;

  return angle;
}

bool te::gm::Rotate(te::gm::Coord2D pr, te::gm::LineString* l, double angle, te::gm::LineString* lOut)
{
  double alfa;
  double dx, dy;
  double x, y;
  double xr, yr;

  try
  {
    if (l->size() < 2)
      return (false);

    //transform Degree to Radius
    alfa = (4.*std::atan(1.)*angle) / 180.;

    dx = pr.getX();
    dy = pr.getY();

    for (std::size_t count = 0; count < l->size(); count++)
    {
      std::unique_ptr<te::gm::Point> curPoint = l->getPointN(count);

      x = curPoint->getX() - dx;
      y = curPoint->getY() - dy;

      xr = x * std::cos(alfa) - y * std::sin(alfa);
      yr = x * std::sin(alfa) + y * std::cos(alfa);

      lOut->setPoint(count, xr + dx, yr + dy);
    }
  }
  catch (...)
  {
    return false;
  }

  return true;
}

bool te::gm::AdjustSegment(te::gm::Point* P0, te::gm::Point* P1, double d0, te::gm::Coord2D& P0out, te::gm::Coord2D& P1out)
{
  double vL_norm1;
  double vL_norm2;

  try
  {
    if (P0->getX() == P1->getX() && P0->getY() == P1->getY())
      return false;

    te::gm::Coord2D vL1((P1->getX() - P0->getX()), (P1->getY() - P0->getY()));
    te::gm::Coord2D vL2(-1 * (P1->getX() - P0->getX()), -1 * (P1->getY() - P0->getY()));
    vL_norm1 = std::sqrt(vL1.getX() * vL1.getX() + vL1.getY() * vL1.getY());
    vL_norm2 = std::sqrt(vL2.getX() * vL2.getX() + vL2.getY() * vL2.getY());

    te::gm::Coord2D uL1((vL1.getX() / vL_norm1), (vL1.getY() / vL_norm1));
    te::gm::Coord2D uL2((vL2.getX() / vL_norm2), (vL2.getY() / vL_norm2));

    te::gm::Point* pFim1 = new te::gm::Point(P0->getX() + uL1.getX() * d0, P0->getY() + uL1.getY() * d0, P0->getSRID());
    te::gm::Point* pFim2 = new te::gm::Point(P0->getX() + uL2.getX() * d0, P0->getY() + uL2.getY() * d0, P1->getSRID());

    if (pFim1->distance(P1) <= pFim2->distance(P1))
    {
      P0out.x = P0->getX();
      P0out.y = P0->getY();
      P1out.x = pFim1->getX();
      P1out.y = pFim1->getY();
    }
    else
    {
      P0out.x = P0->getX();
      P0out.y = P0->getY();
      P1out.x = pFim2->getX();
      P1out.y = pFim2->getY();
    }
  }
  catch (...)
  {
    return false;
  }

  return true;
}

void te::gm::Polygonizer(te::gm::Geometry* g, std::vector<te::gm::Polygon*>& pols)
{
#ifdef TERRALIB_GEOS_ENABLED

  geos::operation::polygonize::Polygonizer polygonizer;

  std::unique_ptr<geos::geom::Geometry> geomGeos(GEOSWriter::write(g));
  polygonizer.add(geomGeos.get());

  std::vector <geos::geom::Polygon*>* vecpolGeos = polygonizer.getPolygons();

  for (std::size_t i = 0; i < vecpolGeos->size(); i++)
  {
    vecpolGeos->at(i)->setSRID(g->getSRID());
    pols.push_back(GEOSReader::read(vecpolGeos->at(i)));
  }

#else
  throw te::common::Exception(TE_TR("Polygonizer routine is supported by GEOS! Please, enable the GEOS support."));
#endif
}

bool te::gm::CheckValidity(const te::gm::Geometry* geom, te::gm::TopologyValidationError& error)
{
#ifdef TERRALIB_GEOS_ENABLED

  try
  {
    std::unique_ptr<geos::geom::Geometry> g(GEOSWriter::write(geom));

    geos::operation::valid::IsValidOp vop(g.get());

    geos::operation::valid::TopologyValidationError* err = vop.getValidationError();

    if (err)
    {
      error.m_coordinate = Coord2D(err->getCoordinate().x, err->getCoordinate().y);
      error.m_message = err->getMessage();

      return false;
    }

    return true;
  }
  catch(const geos::util::GEOSException& e)
  {
    error.m_message = e.what();

    return false;
  }

#else
  throw Exception(TE_TR("isValid routine is supported by GEOS! Please, enable the GEOS support."));
#endif
}

te::gm::Geometry* te::gm::Validate(te::gm::Geometry* geom)
{
#ifdef TERRALIB_GEOS_ENABLED

  std::unique_ptr<geos::geom::Geometry> g(GEOSWriter::write(geom));

  std::vector<te::gm::Geometry*> pAdd;
  geos::operation::polygonize::Polygonizer polygonizer;

  switch (g->getGeometryTypeId())
  {
    case geos::geom::GEOS_POLYGON:
    {
      if (g->isValid())
      {
        g->normalize(); // validate does not pick up rings in the wrong order - this will fix that
        return GEOSReader::read(g.get()); // If the polygon is valid just return it
      }

      AddPolygon(dynamic_cast<te::gm::Polygon*>(GEOSReader::read(g.get())), pAdd);
    }
    break;

    case geos::geom::GEOS_MULTIPOLYGON:
    {
      if (g->isValid()){
        g->normalize(); // validate does not pick up rings in the wrong order - this will fix that
        return GEOSReader::read(g.get()); // If the multipolygon is valid just return it
      }

      for (std::size_t n = g->getNumGeometries(); n-- > 0;)
        AddPolygon(dynamic_cast<te::gm::Polygon*>(GEOSReader::read(g->getGeometryN(n))), pAdd);

    }
    break;

    default:
      return geom;
  }

  for (std::size_t i = 0; i < pAdd.size(); i++)
    polygonizer.add(GEOSWriter::write(pAdd[i]));

  te::common::FreeContents(pAdd);

  std::vector <geos::geom::Polygon*>* vecpolGeos = polygonizer.getPolygons();

  switch (vecpolGeos->size())
  {
    case 1:
      vecpolGeos->at(0)->setSRID(geom->getSRID());
      return GEOSReader::read(vecpolGeos->at(0)); // single polygon - no need to wrap
    default:
    {
      //filter the invalid polygons
      std::size_t i = 0;
      while (i < vecpolGeos->size())
      {
        if (!vecpolGeos->at(i)->isValid())
        {
          vecpolGeos->erase(vecpolGeos->begin() + static_cast<int>(i));
        }
        else
        {
          ++i;
        }
      }

      //polygons may still overlap! Need to sym difference them
      vecpolGeos->at(0)->setSRID(geom->getSRID());
      te::gm::Geometry* ret = GEOSReader::read(vecpolGeos->at(0));

      for (std::size_t p = 1; p < vecpolGeos->size(); p++)
      {
        vecpolGeos->at(p)->setSRID(geom->getSRID());
        ret = ret->symDifference(GEOSReader::read(vecpolGeos->at(p)));
      }

      te::common::FreeContents(*vecpolGeos);

      return ret;
    }
  }

#else
  throw Exception(TE_TR("validate routine is supported by GEOS! Please, enable the GEOS support."));
#endif
}

void te::gm::AddPolygon(te::gm::Polygon* polygon, std::vector<te::gm::Geometry*>& pAdd)
{
#ifdef TERRALIB_GEOS_ENABLED

  AddLineString(dynamic_cast<te::gm::LineString*>(polygon->getExteriorRing()), pAdd);

  for (std::size_t n = polygon->getNumInteriorRings(); n-- > 0;){
    AddLineString(dynamic_cast<te::gm::LineString*>(polygon->getInteriorRingN(n)), pAdd);
  }

#else
  throw Exception(TE_TR("addPolygon routine is supported by GEOS! Please, enable the GEOS support."));
#endif
}

void te::gm::AddLineString(te::gm::LineString* lineString, std::vector<te::gm::Geometry*>& pAdd)
{
#ifdef TERRALIB_GEOS_ENABLED

  std::unique_ptr<geos::geom::LineString> ls(GEOSWriter::write(lineString));

  // unioning the linestring with the point makes any self intersections explicit.
  geos::geom::Point* point = GEOSGeometryFactory::getGeomFactory()->createPoint(ls->getCoordinateN(0));

  geos::geom::Geometry* toAdd = ls->Union(point);

  pAdd.push_back(GEOSReader::read(toAdd));

#else
  throw Exception(TE_TR("addLineString routine is supported by GEOS! Please, enable the GEOS support."));
#endif
}

double te::gm::PerpendicularDistance(const te::gm::Point& first,
                                     const te::gm::Point& last,
                                     te::gm::Point& pin, te::gm::Point& pinter)
{
#ifdef TERRALIB_GEOS_ENABLED
  double d12, xmin, ymin;

  double xi = first.getX();
  double xf = last.getX();
  double yi = first.getY();
  double yf = last.getY();
  double x = pin.getX();
  double y = pin.getY();

  double dx = xf - xi;
  double dy = yf - yi;
  double a2 = (y - yi) * dx - (x - xi)*dy;

  if (dx == 0. && dy == 0.)
  {
    d12 = sqrt(((x - xi) * (x - xi)) + ((y - yi) * (y - yi)));
    d12 *= d12;
  }
  else
    d12 = a2 * a2 / (dx * dx + dy * dy);

  if (dx == 0.)
  {
    xmin = xi;
    ymin = y;
  }
  else if (dy == 0.)
  {
    xmin = x;
    ymin = yi;
  }
  else
  {
    double alfa = dy / dx;
    xmin = (x + alfa * (y - yi) + alfa * alfa * xi) / (1. + alfa * alfa);
    ymin = (x - xmin) / alfa + y;
  }

  pinter = te::gm::Point(xmin, ymin, pin.getSRID());

  return (sqrt(d12));

#else
  throw Exception(TE_TR("PerpendicularDistance routine is supported by GEOS! Please, enable the GEOS support."));
#endif
}

te::gm::Geometry* te::gm::PrepareGeometriesToIntersection(const te::gm::Geometry* geom_A, te::gm::Geometry* geomB, bool& wasChanged)
{
  std::vector<te::gm::Geometry*> vecGeomB;
  vecGeomB.push_back(geomB);

  return te::gm::PrepareGeometriesToIntersection(geom_A, vecGeomB, wasChanged);
}

te::gm::Geometry* te::gm::PrepareGeometriesToIntersection(const te::gm::Geometry* geom_A, const std::vector<te::gm::Geometry*>& vecGeomB, bool& wasChanged)
{
  wasChanged = false;

#ifdef TERRALIB_GEOS_ENABLED

  // Check if MBRs of both LineStrings have intersection.
  const te::gm::Envelope* env_A = geom_A->getMBR();
  te::gm::Envelope expandedEnv_A(env_A->m_llx - TOLERANCE,
    env_A->m_lly - TOLERANCE,
    env_A->m_urx + TOLERANCE,
    env_A->m_ury + TOLERANCE);

  te::gm::Geometry* geomResult = nullptr;

  GeomType geomTypeA = geom_A->getGeomTypeId();

  std::vector<te::gm::MultiLineString*> multiLineStringVectorA = GetAsMultiLineStringVector(*geom_A);


  // Creating the RTree with the ls_B segment boxes intersects ls_A segment box
  te::sam::rtree::Index<std::size_t, 8> rtree;
  std::vector<te::gm::LineString> rtreeSegments;

  //we decompose all geometries in vecGeomB into line segments, and then we index it
  for (std::size_t i = 0; i < vecGeomB.size(); ++i)
  {
    te::gm::Geometry* geometryB = vecGeomB[i];

    const te::gm::Envelope* envGeometryB = geometryB->getMBR();
    te::gm::Envelope expandedGeometryEnv_B(envGeometryB->m_llx - TOLERANCE,
      envGeometryB->m_lly - TOLERANCE,
      envGeometryB->m_urx + TOLERANCE,
      envGeometryB->m_ury + TOLERANCE);

    if (expandedEnv_A.intersects(expandedGeometryEnv_B) == false)
    {
      continue;
    }
    
    std::vector<te::gm::MultiLineString*> multiLineStringVectorB = GetAsMultiLineStringVector(*vecGeomB[i]);
    for (std::size_t j = 0; j < multiLineStringVectorB.size(); ++j)
    {
      te::gm::MultiLineString* multiLineB = multiLineStringVectorB[j];
      const std::vector<te::gm::Geometry*>& vecGeometries = multiLineB->getGeometries();

      for (std::size_t k = 0; k < vecGeometries.size(); ++k)
      {
        te::gm::LineString* lineB = static_cast<te::gm::LineString*>(vecGeometries[k]);
        
        const te::gm::Envelope* envLineB = lineB->getMBR();
        te::gm::Envelope expandedLineEnv_B(envLineB->m_llx - TOLERANCE,
          envLineB->m_lly - TOLERANCE,
          envLineB->m_urx + TOLERANCE,
          envLineB->m_ury + TOLERANCE);

        if (expandedEnv_A.intersects(expandedLineEnv_B) == false)
        {
          continue;
        }
        
        //finally we get all the segments
        std::size_t nSegmentsB = lineB->size() - 1;
        for (std::size_t s = 0; s < nSegmentsB; ++s)
        {

          te::gm::Envelope currentSegmentEnvelope = CreateEnvelope(lineB->getX(s), lineB->getY(s), lineB->getX(s + 1), lineB->getY(s + 1));

          te::gm::Envelope expandedSegmentEnv_B(currentSegmentEnvelope.getLowerLeftX() - TOLERANCE,
            currentSegmentEnvelope.getLowerLeftY() - TOLERANCE,
            currentSegmentEnvelope.getUpperRightX() + TOLERANCE,
            currentSegmentEnvelope.getUpperRightY() + TOLERANCE);

          if (expandedEnv_A.intersects(expandedSegmentEnv_B) == true)
          {
            te::gm::LineString segment_B(2, te::gm::LineStringType, geom_A->getSRID());

            segment_B.setPoint(0, lineB->getX(s), lineB->getY(s));
            segment_B.setPoint(1, lineB->getX(s + 1), lineB->getY(s + 1));

            rtree.insert(currentSegmentEnvelope, rtreeSegments.size());
            rtreeSegments.push_back(segment_B);
          }
        }
      }
    }
    te::common::FreeContents(multiLineStringVectorB);
  }

  for (std::size_t i = 0; i < multiLineStringVectorA.size(); ++i)
  {
    bool currentWasChanged = false;
    std::unique_ptr<te::gm::MultiLineString> currentResult(PrepareGeometriesToIntersection(*multiLineStringVectorA[i], rtree, rtreeSegments, currentWasChanged));

    if (currentWasChanged == true)
    {
      wasChanged = true;
      delete multiLineStringVectorA[i];
      multiLineStringVectorA[i] = currentResult.release();
    }
  }

  if (wasChanged)
  {
    std::unique_ptr<te::gm::Geometry> preparedGeometryResult(MultiLineToDefinedType(multiLineStringVectorA, geomTypeA));

    TopologyValidationError error;
    bool result = CheckValidity(preparedGeometryResult.get(), error);
    if (result == false)
    {
      preparedGeometryResult.reset(preparedGeometryResult->buffer(0));
    }

    geomResult = preparedGeometryResult.release();
  }
  else
  {
    geomResult = dynamic_cast<te::gm::Geometry*>(geom_A->clone());
  }


  te::common::FreeContents(multiLineStringVectorA);

  return geomResult;

#else
  throw Exception(TE_TR("PrepareGeometriesToIntersection routine is supported by GEOS! Please, enable the GEOS support."));
#endif
}

te::gm::MultiLineString* te::gm::PrepareGeometriesToIntersection(const te::gm::MultiLineString& mline_A,
                                                                 const te::sam::rtree::Index<std::size_t, 8>& rtree,
                                                                 std::vector<te::gm::LineString>& vecRtreeSegments,
                                                                 bool& wasChanged)
{
  wasChanged = false;

#ifdef TERRALIB_GEOS_ENABLED

  te::gm::MultiLineString* multiLineStringResult = new te::gm::MultiLineString(0, mline_A.getGeomTypeId(), mline_A.getSRID());

  const std::vector<te::gm::Geometry*>& geometryVectorA = mline_A.getGeometries();

  for (std::size_t i = 0; i < geometryVectorA.size(); ++i)
  {
    bool currentWasChanged = false;

    te::gm::Geometry* singleResult = PrepareGeometriesToIntersection(*dynamic_cast<te::gm::LineString*>(geometryVectorA[i]), rtree, vecRtreeSegments, currentWasChanged);
    multiLineStringResult->add(singleResult);

    if (currentWasChanged == true)
    {
      wasChanged = true;
    }
  }

  return multiLineStringResult;

#else
  throw Exception(TE_TR("PrepareGeometriesToIntersection routine is supported by GEOS! Please, enable the GEOS support."));
#endif
}

te::gm::LineString* te::gm::PrepareGeometriesToIntersection(const te::gm::LineString& ls_A,
                                                            const te::sam::rtree::Index<std::size_t, 8>& rtree,
                                                            std::vector<te::gm::LineString>& vecRtreeSegments, 
                                                            bool& wasChanged)
{
  wasChanged = false;

#ifdef TERRALIB_GEOS_ENABLED

  // First try to snap the points from ls_A in the points of ls_B
  bool wasSnapped = false;

  //if there are none candidate segments, we avoid some processing
  std::vector<std::size_t> report;
  rtree.search(*ls_A.getMBR(), report);

  if (report.empty())
  {
    return static_cast<te::gm::LineString*>(ls_A.clone());
  }

  std::size_t n_pointsA = ls_A.getNPoints();
  std::size_t n_segments_A = n_pointsA - 1;

  std::vector<te::gm::Point> vecOutputPoints;
  vecOutputPoints.reserve(n_pointsA);

  std::set<std::size_t> setIndexesIgnored;
  for (std::size_t i = 0; i < n_segments_A; ++i)
  {
    te::gm::Point segmentP1(ls_A.getX(i), ls_A.getY(i));
    te::gm::Point segmentP2(ls_A.getX(i+1), ls_A.getY(i+1));

    te::gm::Envelope segmentEnvelope = CreateEnvelope(segmentP1.getX(), segmentP1.getY(), segmentP2.getX(), segmentP2.getY());

    //now we create an ignore list based in all candidate segments that are EQUAL to the segments in the prepared line
    std::vector<std::size_t> segmentReport;
    rtree.search(segmentEnvelope, segmentReport);

    for (std::size_t i = 0; i < segmentReport.size(); ++i)
    {
      std::size_t index = segmentReport[i];
      //if it is already ignored, we go to the next one
      if (setIndexesIgnored.find(index) != setIndexesIgnored.end())
      {
        continue;
      }

      const te::gm::LineString& candidate = vecRtreeSegments[index];

      //if the boxes arent igual, the segment must not be ignored
      if(segmentEnvelope.equals(*candidate.getMBR()) == false)
      {
        continue;
      }

      //if the boxes are equal, we must check the coordinates  
      //if one coordinate is equal, necessarely the other one is tool
      //so we can check only  the first coordinate
      te::gm::Point candidateP1(candidate.getX(0), candidate.getY(0));
      te::gm::Point candidateP2(candidate.getX(1), candidate.getY(1));
      
      //if the segment fist point is not equal any point of the candidate line, the segment must not be ignored
      if (IsEqual(segmentP1, candidateP1) == false && IsEqual(segmentP1, candidateP2) == false)
      {
        continue;
      }

      //if we got here, both segments are equal and so they can be added to the ignore list
      setIndexesIgnored.insert(index);
    }
  }

  std::unique_ptr<te::gm::LineString> lr_Result(dynamic_cast<te::gm::LineString*>(ls_A.clone()));
  lr_Result.reset(SnapLineToPoints(rtree, vecRtreeSegments, setIndexesIgnored, *lr_Result.get(), wasSnapped));

  lr_Result.reset(AddIntersectionPoints(rtree, vecRtreeSegments, setIndexesIgnored, *lr_Result.get(), true));
  lr_Result.reset(AddIntersectionPoints(rtree, vecRtreeSegments, setIndexesIgnored, *lr_Result.get()));

  if (wasSnapped == true || lr_Result->size() != ls_A.size())
  {
    wasChanged = true;
  }

  return lr_Result.release();

#else
  throw Exception(TE_TR("PrepareGeometriesToIntersection routine is supported by GEOS! Please, enable the GEOS support."));
#endif
}

te::gm::LineString* te::gm::SnapLineToPoints(const te::sam::rtree::Index<std::size_t, 8>& rtree,
                                             const std::vector<te::gm::LineString>& referenceSegments,
                                             const std::set<std::size_t>& setIndexesIgnored,
                                             const LineString& linearRingToSnap, bool& wasChanged)
{
  wasChanged = false;

#ifdef TERRALIB_GEOS_ENABLED

  te::gm::LineString* lr_Result = dynamic_cast<te::gm::LineString*>(linearRingToSnap.clone());

  std::size_t npoints_A = linearRingToSnap.getNPoints();

  for (std::size_t i = 0; i < npoints_A; ++i)
  {
    std::unique_ptr<Point> currentPoint = lr_Result->getPointN(i);

    te::gm::Envelope expandedEnv;
    expandedEnv.init(
        currentPoint->getX() - TOLERANCE, currentPoint->getY() - TOLERANCE,
        currentPoint->getX() + TOLERANCE, currentPoint->getY() + TOLERANCE);

    std::vector<size_t> report;
    rtree.search(expandedEnv, report);

    if (report.empty())
      continue;

    // Extracts all the points.
    std::vector<te::gm::Point> candidates;
    for (std::size_t ls = 0; ls < report.size(); ++ls)
    {
      std::size_t index = report[ls];

      //we the candidate is in the ignore list, we go to the next one
      if (setIndexesIgnored.find(index) != setIndexesIgnored.end())
      {
        continue;
      }

      std::unique_ptr<Point> pStart = referenceSegments[index].getPointN(0);
      candidates.push_back(*pStart.get());

      std::unique_ptr<Point> pEnd = referenceSegments[index].getPointN(1);
      candidates.push_back(*pEnd.get());
    }

    if (candidates.empty() == true)
    {
      continue;
    }

    //Checks if there is on point that is within the tolerance.
    //If there is more than one, is found the closest.
    DistanceOrderFunctor func(*currentPoint);
    std::sort(candidates.begin(), candidates.end(), func);

    double distance = currentPoint->distance(&candidates[0]);
    if(distance != 0 && distance <= TOLERANCE)
    {
      wasChanged = true;

      lr_Result->setPointN(i, candidates[0]);
      
      if (i == 0)
      {
        lr_Result->setPointN(npoints_A - 1, candidates[0]);
      }
    }
  }

  return lr_Result;

#else
  throw Exception(TE_TR("SnapLineToPoints routine is supported by GEOS! Please, enable the GEOS support."));
#endif
}

te::gm::LineString* te::gm::AddIntersectionPoints(const te::sam::rtree::Index<std::size_t, 8>& rtree,
                                                  const std::vector<te::gm::LineString>& candidateSegments,
                                                  const std::set<std::size_t>& setIndexesIgnored,
                                                  const LineString& lr_Reference,
                                                  const bool& usePerpendicularDistance)
{
#ifdef TERRALIB_GEOS_ENABLED

  bool addedIntersectionPoints = false;

  std::size_t n_pointsA = lr_Reference.getNPoints();
  std::size_t n_segments_A = n_pointsA - 1;

  std::vector<te::gm::Point> vecOutputPoints;
  vecOutputPoints.reserve(n_pointsA);

  for (std::size_t i = 0; i < n_segments_A; ++i)
  {
    // Search for candidates.
    te::gm::Envelope currentSegmentEnvelope = CreateEnvelope(lr_Reference.getX(i), lr_Reference.getY(i), lr_Reference.getX(i + 1), lr_Reference.getY(i + 1));

    std::unique_ptr<Point> pCurrentStartPoint = lr_Reference.getPointN(i);

    vecOutputPoints.push_back(*pCurrentStartPoint.get());

    te::gm::Envelope expandedEnv(currentSegmentEnvelope.getLowerLeftX() - TOLERANCE,
      currentSegmentEnvelope.getLowerLeftY() - TOLERANCE,
      currentSegmentEnvelope.getUpperRightX() + TOLERANCE,
      currentSegmentEnvelope.getUpperRightY() + TOLERANCE);

    std::vector<size_t> report;
    rtree.search(expandedEnv, report);

    if (report.empty())
      continue;

    //creates the curent segment
    te::gm::LineString segment_A(2, te::gm::LineStringType, lr_Reference.getSRID());
    segment_A.setPoint(0, lr_Reference.getX(i), lr_Reference.getY(i));
    segment_A.setPoint(1, lr_Reference.getX(i + 1), lr_Reference.getY(i + 1));

    std::vector<te::gm::LineString> candidates;
    candidates.reserve(report.size());
    for(std::size_t ls = 0; ls < report.size(); ++ls)
    {
      std::size_t index = report[ls];

      //we the candidate is in the ignore list, we go to the next one
      if (setIndexesIgnored.find(index) != setIndexesIgnored.end())
      {
        continue;
      }
      
      const te::gm::LineString& segmentCandidate = candidateSegments[index];
      candidates.push_back(segmentCandidate);
    }

    if (candidates.empty() == true)
    {
      continue;
    }

    // Executes the intersection for each candidate point, and get the coordinates.
    std::vector<te::gm::Point> vecIntersectionPoints;

    if(usePerpendicularDistance)
      vecIntersectionPoints = GetIntersectionPointsByPerpendicularDistance(segment_A, candidates);
    else
      vecIntersectionPoints = GetIntersectionPointsByOperators(segment_A, candidates);

    DistanceOrderFunctor func(*pCurrentStartPoint.get());
    std::sort(vecIntersectionPoints.begin(), vecIntersectionPoints.end(), func);

    for (std::size_t j = 0; j < vecIntersectionPoints.size(); ++j)
    {
      const te::gm::Point& currentIntersectionPoint = vecIntersectionPoints[j];
      const te::gm::Point& lastPoint = vecOutputPoints[vecOutputPoints.size() - 1];

      if (IsEqual(currentIntersectionPoint, lastPoint, TOLERANCE) == true)
      {
        continue;
      }
      if (IsEqual(currentIntersectionPoint, *pCurrentStartPoint.get(), TOLERANCE) == true)
      {
        continue;
      }
      std::unique_ptr<Point> pCurrentEndPoint = lr_Reference.getPointN(i + 1);
      if (IsEqual(currentIntersectionPoint, *pCurrentEndPoint.get(), TOLERANCE) == true)
      {
        continue;
      }

      addedIntersectionPoints = true;
      vecOutputPoints.push_back(vecIntersectionPoints[j]);
    }
  }

  //we onlye create anew line if we detected at least one distinct intersection point
  if (addedIntersectionPoints == false)
  {
    return dynamic_cast<te::gm::LineString*>(lr_Reference.clone());
  }

  // Set the end point in the output vector.
  std::unique_ptr<Point> pStart = lr_Reference.getStartPoint();
  vecOutputPoints.push_back(*pStart.get());

  te::gm::LineString* result = new te::gm::LineString(vecOutputPoints.size(), te::gm::LineStringType, lr_Reference.getSRID());

  for (std::size_t i = 0; i < vecOutputPoints.size(); ++i)
    result->setPointN(i, vecOutputPoints[i]);

  return result;

#else
  throw Exception(TE_TR("GetIntersectionPoints routine is supported by GEOS! Please, enable the GEOS support."));
#endif
}

std::vector<te::gm::Point> te::gm::GetIntersectionPointsByOperators(
    const te::gm::LineString& lsReference,
    const std::vector<te::gm::LineString>& candidates)
{
#ifdef TERRALIB_GEOS_ENABLED

  std::vector<te::gm::Point> intersectionPoints;
  std::unique_ptr<te::gm::Point> startPoint;
  std::unique_ptr<te::gm::Point> endPoint;

  geos::geom::Coordinate referenceP1(lsReference.getX(0), lsReference.getY(0));
  geos::geom::Coordinate referenceP2(lsReference.getX(1), lsReference.getY(1));

  for (std::size_t i = 0; i < candidates.size(); ++i)
  {
    const te::gm::LineString& lineCandidate = candidates[i];

    geos::geom::Coordinate candidateP1(lineCandidate.getX(0), lineCandidate.getY(0));
    geos::geom::Coordinate candidateP2(lineCandidate.getX(1), lineCandidate.getY(1));

    geos::algorithm::LineIntersector intersector;
    intersector.computeIntersection(referenceP1, referenceP2, candidateP1, candidateP2);
    if (intersector.hasIntersection() == false)
    {
      continue;
    }

    std::vector<te::gm::Point> candidatesIntersectionPoints;
    int numberOfIntersections = intersector.getIntersectionNum();
    for (int j = 0; j < numberOfIntersections; ++j)
    {
      const geos::geom::Coordinate& interCoord = intersector.getIntersection(j);

      candidatesIntersectionPoints.push_back(te::gm::Point(interCoord.x, interCoord.y, lsReference.getSRID()));
    }

    for (std::size_t j = 0; j < candidatesIntersectionPoints.size(); ++j)
    {
      if (startPoint.get() == nullptr || endPoint.get() == nullptr)
      {
        startPoint = lsReference.getStartPoint();
        endPoint = lsReference.getEndPoint();
      }      

      const te::gm::Point& candidateIntersectionPoint = candidatesIntersectionPoints[j];
      
      //if the candidate intersection point is "near" to the start or end point of the current segment, we drop it
      if (IsEqual(candidateIntersectionPoint, *startPoint.get(), TOLERANCE) == true || IsEqual(candidateIntersectionPoint, *endPoint.get(), TOLERANCE) == true)
      {
        continue;
      }

      std::unique_ptr<te::gm::Point>candidateStartPoint = candidates[i].getStartPoint();
      std::unique_ptr<te::gm::Point>candidateEndPoint = candidates[i].getEndPoint();

      // If the calculated intersection point is "near" the start or end point of the current segment, just snap it.
      double distanceStartPointB = candidateIntersectionPoint.distance(candidateStartPoint.get());
      double distanceEndPointB = candidateIntersectionPoint.distance(candidateEndPoint.get());

      if (distanceStartPointB <= TOLERANCE)
      {
        intersectionPoints.push_back(*candidateStartPoint.get());
      }
      else if (distanceEndPointB <= TOLERANCE)
      {
        intersectionPoints.push_back(*candidateEndPoint.get());
      }
      else
      {
        intersectionPoints.push_back(candidateIntersectionPoint);
      }
    }
  }

  return intersectionPoints;

#else
  throw Exception(TE_TR("GetIntersectionPointsByOperators routine is supported by GEOS! Please, enable the GEOS support."));
#endif
}

std::vector<te::gm::Point> te::gm::GetIntersectionPointsByPerpendicularDistance(
    const te::gm::LineString& lsReference,
    const std::vector<te::gm::LineString>& candidates)
{
#ifdef TERRALIB_GEOS_ENABLED

  std::vector<te::gm::Point> intersectionPoints;

  std::unique_ptr<te::gm::Point> startPoint = lsReference.getStartPoint();
  std::unique_ptr<te::gm::Point> endPoint = lsReference.getEndPoint();

  for (std::size_t i = 0; i < candidates.size(); ++i)
  {
    std::unique_ptr<te::gm::Point>candidateStartPoint = candidates[i].getStartPoint();
    std::unique_ptr<te::gm::Point>candidateEndPoint = candidates[i].getEndPoint();

    std::unique_ptr<geos::geom::LineSegment> lineSegment(
        new geos::geom::LineSegment(startPoint->getX(), startPoint->getY(),
                                    endPoint->getX(), endPoint->getY()));

    //we try to find the an intersection point touching our current segment
    //we first compare our segment to the candidate start point
    te::gm::Point pInterStart;
    double startDistance = te::gm::PerpendicularDistance(*startPoint.get(), *endPoint.get(), *candidateStartPoint.get(), pInterStart);

    if (IsEqual(*startPoint.get(), pInterStart, TOLERANCE) == false && IsEqual(*endPoint.get(), pInterStart, TOLERANCE) == false)
    {
      std::unique_ptr<geos::geom::Coordinate> coordStart(new geos::geom::Coordinate(pInterStart.getX(), pInterStart.getY()));

      double projFactorStart = lineSegment->projectionFactor(*coordStart.get());
      if (startDistance <= TOLERANCE && (projFactorStart >= 0 && projFactorStart <= 1))
      {
        intersectionPoints.push_back(*candidateStartPoint.get());
      }
    }
     

    //then we compare our segment to the candidate end point
    te::gm::Point pInterEnd;
    double endDistance = te::gm::PerpendicularDistance(*startPoint.get(), *endPoint.get(), *candidateEndPoint.get(), pInterEnd);

    if (IsEqual(*startPoint.get(), pInterEnd, TOLERANCE) == false && IsEqual(*endPoint.get(), pInterEnd, TOLERANCE) == false)
    {
      std::unique_ptr<geos::geom::Coordinate> coordEnd(new geos::geom::Coordinate(pInterEnd.getX(), pInterEnd.getY()));

      double projFactorEnd = lineSegment->projectionFactor(*coordEnd.get());
      if (endDistance <= TOLERANCE && (projFactorEnd >= 0 && projFactorEnd <= 1))
      {
        intersectionPoints.push_back(*candidateEndPoint.get());
      }
    }
  }

  return intersectionPoints;

#else
  throw Exception(TE_TR("GetIntersectionPointsByPerpendicularDistance routine is supported by GEOS! Please, enable the GEOS support."));
#endif
}

std::vector<te::gm::MultiLineString*> te::gm::GetAsMultiLineStringVector(const te::gm::Geometry& geometry)
{
  std::vector<te::gm::MultiLineString*> multiLineStringVector;

  switch(geometry.getGeomTypeId())
  {
    case te::gm::LineStringType:
    case te::gm::LineStringMType:
    case te::gm::LineStringZType:
    case te::gm::LineStringZMType:
    {
      te::gm::MultiLineString* multiLineString = new te::gm::MultiLineString(0, GetMultiLineStringType(geometry), geometry.getSRID());

      te::gm::LineString* lineString = dynamic_cast<te::gm::LineString*>(geometry.clone());

      multiLineString->add(lineString);

      multiLineStringVector.push_back(multiLineString);
    }
    break;

    case te::gm::MultiLineStringType:
    case te::gm::MultiLineStringMType:
    case te::gm::MultiLineStringZType:
    case te::gm::MultiLineStringZMType:
    {
      te::gm::MultiLineString* multiLineString = new te::gm::MultiLineString(0, GetMultiLineStringType(geometry), geometry.getSRID());

      std::vector<te::gm::Geometry*> geometryVector;
      te::gm::Multi2Single(&geometry, geometryVector);

      for(auto g : geometryVector)
      {
        te::gm::LineString* lineString = dynamic_cast<te::gm::LineString*>(g->clone());

        multiLineString->add(lineString);

        multiLineStringVector.push_back(multiLineString);
      }

      te::common::FreeContents(geometryVector);
    }
    break;

    case te::gm::PolygonType:
    case te::gm::PolygonMType:
    case te::gm::PolygonZType:
    case te::gm::PolygonZMType:
    {
      te::gm::MultiLineString* multiLineString = new te::gm::MultiLineString(0, GetMultiLineStringType(geometry), geometry.getSRID());

      const te::gm::Polygon* polygon = dynamic_cast<const te::gm::Polygon*>(&geometry);

      const std::vector<te::gm::Curve*>& curveVector = polygon->getRings();

      for (auto curve : curveVector)
      {
        te::gm::LineString* lineString = dynamic_cast<te::gm::LineString*>(curve->clone());

        multiLineString->add(lineString);
      }

      multiLineStringVector.push_back(multiLineString);
    }
    break;

    case te::gm::MultiPolygonType:
    case te::gm::MultiPolygonMType:
    case te::gm::MultiPolygonZType:
    case te::gm::MultiPolygonZMType:
    {
      std::vector<te::gm::Geometry*> geometryVector;
      Multi2Single(&geometry, geometryVector);

      for(auto g : geometryVector)
      {
        std::vector<te::gm::MultiLineString*> multiLineStringVectorTemp = GetAsMultiLineStringVector(*g);

        multiLineStringVector.insert(multiLineStringVector.end(), multiLineStringVectorTemp.begin(), multiLineStringVectorTemp.end());
      }

      te::common::FreeContents(geometryVector);
    }
    break;

    default:
      throw Exception("This type of geometry is not composed by LineString.");
  }

  return multiLineStringVector;
}

te::gm::Geometry* te::gm::MultiLineToDefinedType(const std::vector<te::gm::MultiLineString*>& multiLineStringVector, te::gm::GeomType geometryType)
{
  te::gm::Geometry* geometry = nullptr;

  if(multiLineStringVector.empty())
    return geometry;

  int srid = multiLineStringVector[0]->getSRID();

  switch (geometryType)
  {
    case te::gm::LineStringType:
    case te::gm::LineStringZType:
    case te::gm::LineStringMType:
    case te::gm::LineStringZMType:
    {
      te::gm::MultiLineString* mls = multiLineStringVector[0];

      if(mls->isEmpty())
        return geometry;

      geometry = dynamic_cast<te::gm::Geometry*>(mls->getGeometryN(0)->clone());
    }
    break;

    case te::gm::PolygonType:
    case te::gm::PolygonZType:
    case te::gm::PolygonMType:
    case te::gm::PolygonZMType:
    {
      te::gm::MultiLineString* mls = multiLineStringVector[0];

      if(mls->isEmpty())
        return geometry;

      te::gm::Polygon* polygon = new te::gm::Polygon(0, geometryType, srid);

      for(auto g : mls->getGeometries())
        polygon->add(dynamic_cast<te::gm::Curve*>(g->clone()));

      geometry = polygon;
    }
    break;

    case te::gm::MultiLineStringType:
    case te::gm::MultiLineStringZType:
    case te::gm::MultiLineStringMType:
    case te::gm::MultiLineStringZMType:
    {
      te::gm::MultiLineString* outMls = new te::gm::MultiLineString(0, geometryType, srid);

      for(auto inMls : multiLineStringVector)
      {
        if(inMls->isEmpty())
          continue;

        for(auto g : inMls->getGeometries())
          outMls->add(dynamic_cast<te::gm::Geometry*>(g->clone()));
      }

      geometry = outMls;
    }
    break;

    case te::gm::MultiPolygonType:
    case te::gm::MultiPolygonZType:
    case te::gm::MultiPolygonMType:
    case te::gm::MultiPolygonZMType:
    {
      te::gm::MultiPolygon* outMPolygon = new te::gm::MultiPolygon(0, geometryType, srid);

      for(auto inMls : multiLineStringVector)
      {
        if(inMls->getGeometries().empty())
          continue;

        std::size_t numRings = inMls->getGeometries().size();

        te::gm::Polygon* polygon = new te::gm::Polygon(numRings, GetSimpleType(geometryType), srid);

        for (std::size_t j = 0; j < numRings; ++j)
        {
          polygon->setRingN(j, dynamic_cast<te::gm::Curve*>(inMls->getGeometryN(j)->clone()));
        }

        outMPolygon->add(polygon);
      }

      geometry = outMPolygon;
    }
    break;

    default:
      throw Exception("This type of geometry is not composed by LineString.");
  }

  return geometry;
}

te::gm::GeomType te::gm::GetMultiLineStringType(const te::gm::Geometry& geometry)
{
  switch (geometry.getGeomTypeId())
  {
    case te::gm::LineStringZType:
    case te::gm::PolygonZType:
    case te::gm::MultiLineStringZType:
    case te::gm::MultiPolygonZType:
      return te::gm::MultiLineStringZType;

    case te::gm::LineStringMType:
    case te::gm::PolygonMType:
    case te::gm::MultiLineStringMType:
    case te::gm::MultiPolygonMType:
      return te::gm::MultiLineStringMType;

    case te::gm::LineStringZMType:
    case te::gm::PolygonZMType:
    case te::gm::MultiLineStringZMType:
    case te::gm::MultiPolygonZMType:
      return te::gm::MultiLineStringZMType;

    default:
      return te::gm::MultiLineStringType;
  }
}

bool te::gm::IsMultiType(te::gm::GeomType geomType)
{
  switch(geomType)
  {
    case te::gm::MultiLineStringType:
    case te::gm::MultiLineStringMType:
    case te::gm::MultiLineStringZType:
    case te::gm::MultiLineStringZMType:
    case te::gm::MultiPointType:
    case te::gm::MultiPointMType:
    case te::gm::MultiPointZType:
    case te::gm::MultiPointZMType:
    case te::gm::MultiPolygonType:
    case te::gm::MultiPolygonMType:
    case te::gm::MultiPolygonZType:
    case te::gm::MultiPolygonZMType:
    case te::gm::GeometryCollectionType:
    case te::gm::GeometryCollectionMType:
    case te::gm::GeometryCollectionZType:
    case te::gm::GeometryCollectionZMType:
      return true;
    default:
      return false;
  }
}

te::gm::GeomType te::gm::GetSimpleType(te::gm::GeomType geomType)
{
  switch(geomType)
  {
    case te::gm::MultiLineStringType:
      return te::gm::LineStringType;
    case te::gm::MultiLineStringMType:
      return te::gm::LineStringMType;
    case te::gm::MultiLineStringZType:
      return te::gm::LineStringZType;
    case te::gm::MultiLineStringZMType:
      return te::gm::LineStringZMType;
    case te::gm::MultiPointType:
      return te::gm::PointType;
    case te::gm::MultiPointMType:
      return te::gm::PointMType;
    case te::gm::MultiPointZType:
      return te::gm::PointZType;
    case te::gm::MultiPointZMType:
      return te::gm::PointZMType;
    case te::gm::MultiPolygonType:
      return te::gm::PolygonType;
    case te::gm::MultiPolygonMType:
      return te::gm::PolygonMType;
    case te::gm::MultiPolygonZType:
      return te::gm::PolygonZType;
    case te::gm::MultiPolygonZMType:
      return te::gm::PolygonZMType;
    case te::gm::GeometryCollectionType:
      return te::gm::GeometryType;
    case te::gm::GeometryCollectionMType:
      return te::gm::GeometryMType;
    case te::gm::GeometryCollectionZType:
      return te::gm::GeometryZType;
    case te::gm::GeometryCollectionZMType:
      return te::gm::GeometryZMType;
    default:
      return te::gm::UnknownGeometryType;
  }
}

te::gm::GeomType te::gm::GetMultiType(te::gm::GeomType geomType)
{
  switch (geomType)
  {
  case te::gm::LineStringType:
    return te::gm::MultiLineStringType;
  case te::gm::LineStringMType:
    return te::gm::MultiLineStringMType;
  case te::gm::LineStringZType:
    return te::gm::MultiLineStringZType;
  case te::gm::LineStringZMType:
    return te::gm::MultiLineStringZMType;
  case te::gm::PointType:
    return te::gm::MultiPointType;
  case te::gm::PointMType:
    return te::gm::MultiPointMType;
  case te::gm::PointZType:
    return te::gm::MultiPointZType;
  case te::gm::PointZMType:
    return te::gm::MultiPointZMType;
  case te::gm::PolygonType:
    return te::gm::MultiPolygonType;
  case te::gm::PolygonMType:
    return te::gm::MultiPolygonMType;
  case te::gm::PolygonZType:
    return te::gm::MultiPolygonZType;
  case te::gm::PolygonZMType:
    return te::gm::MultiPolygonZMType;
  case te::gm::GeometryType:
    return te::gm::GeometryCollectionType;
  case te::gm::GeometryMType:
    return te::gm::GeometryCollectionMType;
  case te::gm::GeometryZType:
    return te::gm::GeometryCollectionZType;
  case te::gm::GeometryZMType:
    return te::gm::GeometryCollectionZMType;
  default:
    return te::gm::UnknownGeometryType;
  }
}

te::gm::LineString te::gm::CreateLine(const te::gm::Coord2D& startCoord, const te::gm::Coord2D& endCoord, const int& srid, const int& numberOfIntermediateCoords)
{
  double gapX = (endCoord.getX() - startCoord.getX()) / static_cast<double>(numberOfIntermediateCoords + 1);
  double gapY = (endCoord.getY() - startCoord.getY()) / static_cast<double>(numberOfIntermediateCoords + 1);

  std::vector<te::gm::Coord2D> coords;

  coords.push_back(startCoord);
  for (std::size_t i = 0; i < static_cast<std::size_t>(numberOfIntermediateCoords); ++i)
  {
    std::size_t coordIndex = i + 1;
    te::gm::Coord2D coord(startCoord.getX() + (coordIndex * gapX), startCoord.getY() + (coordIndex * gapY));
    coords.push_back(coord);
  }
  coords.push_back(endCoord);

  te::gm::LineString resultLine(coords.size(), te::gm::LineStringType, srid);

  for (std::size_t t = 0; t < coords.size(); ++t)
    resultLine.setPoint(t, coords[t].getX(), coords[t].getY());

  return resultLine;
}

te::gm::Polygon te::gm::CreatePolygon(const te::gm::Envelope& box, const int& srid, const int& numberOfIntermediateCoords)
{
  te::gm::Coord2D lowerLeft(box.m_llx, box.m_lly);
  te::gm::Coord2D lowerRight(box.m_urx, box.m_lly);

  te::gm::Coord2D upperLeft(box.m_llx, box.m_ury);
  te::gm::Coord2D upperRight(box.m_urx, box.m_ury);

  te::gm::LineString bottomLine = CreateLine(lowerLeft, lowerRight, srid, numberOfIntermediateCoords);
  te::gm::LineString leftLine = CreateLine(lowerRight, upperRight, srid, numberOfIntermediateCoords);
  te::gm::LineString topLine = CreateLine(upperRight, upperLeft, srid, numberOfIntermediateCoords);
  te::gm::LineString rightLine = CreateLine(upperLeft, lowerLeft, srid, numberOfIntermediateCoords);

  std::vector<te::gm::Coord2D> coords;

  //we do not add the last coord to avoid duplicity
  for (std::size_t i = 0; i < bottomLine.size() - 1; ++i)
    coords.push_back(te::gm::Coord2D(bottomLine.getX(i), bottomLine.getY(i)));

  for (std::size_t i = 0; i < leftLine.size() - 1; ++i)
    coords.push_back(te::gm::Coord2D(leftLine.getX(i), leftLine.getY(i)));

  for (std::size_t i = 0; i < topLine.size() - 1; ++i)
    coords.push_back(te::gm::Coord2D(topLine.getX(i), topLine.getY(i)));

  for (std::size_t i = 0; i < rightLine.size() - 1; ++i)
    coords.push_back(te::gm::Coord2D(rightLine.getX(i), rightLine.getY(i)));

  te::gm::LineString* fullLine = new te::gm::LineString(coords.size(), te::gm::LineStringType);

  for (std::size_t t = 0; t < coords.size(); ++t)
    fullLine->setPoint(t, coords[t].getX(), coords[t].getY());

  te::gm::Polygon polygon(1, te::gm::PolygonType, srid);
  polygon.setRingN(0, fullLine);

  return polygon;
}