src/terralib/mnt/core/Utils.cpp

Go to the documentation of this file.

/*!
\file terralib/mnt/core/Utils.cpp

\brief Utility functions for MNT support.
*/

#include "Utils.h"
#include "CalculateGrid.h"
#include "SplineGrass.h"

#include "../../common/progress/TaskProgress.h"
#include "../../dataaccess/datasource/DataSourceTransactor.h"
#include "../../dataaccess/utils/Utils.h"
#include "../../datatype/SimpleProperty.h"
#include "../../geometry/GeometryProperty.h"
#include "../../geometry/Point.h"
#include "../../memory/DataSet.h"
#include "../../memory/DataSetItem.h"
#include "../../raster.h"
#include "../../vp/Utils.h"

//#include <geos.h>
#include <geos/geom/Coordinate.h>
#include <geos/simplify/DouglasPeuckerLineSimplifier.h>

#include <cmath>
#include <iostream>
#include <limits>
#include <stdint.h>


size_t te::mnt::ReadPoints(std::string &inDsetName, te::da::DataSourcePtr &inDsrc, std::string &atrZ, double tol, 
  te::gm::MultiPoint &mpt, std::string &geostype, te::gm::Envelope &env)
{
  if (inDsetName.empty())
    return 0;

  size_t nsamples = 0;

   std::unique_ptr<te::da::DataSet> inDset(inDsrc->getDataSet(inDsetName));

  std::size_t geo_pos = te::da::GetFirstPropertyPos(inDset.get(), te::dt::GEOMETRY_TYPE);

  te::common::TaskProgress task("Reading Samples...", te::common::TaskProgress::UNDEFINED, (int)inDset->size());

  inDset->moveBeforeFirst();
  double value=std::numeric_limits<double>::max();
  while (inDset->moveNext())
  {
    if (!task.isActive())
      return false;
    task.pulse();

    std::unique_ptr<te::gm::Geometry> gin(inDset->getGeometry(geo_pos));
    geostype = gin.get()->getGeometryType();

    if (geostype == "MultiPoint")
    {
      te::gm::MultiPoint *g = dynamic_cast<te::gm::MultiPoint*>(gin.get());
      std::size_t np = g->getNumGeometries();
      if (!atrZ.empty())
        value = inDset->getDouble(atrZ);
      for (std::size_t i = 0; i < np; ++i)
      {
        te::gm::Point *p = dynamic_cast<te::gm::Point*>(g->getGeometryN(i));
        if (atrZ.empty())
          value = p->getZ();

        te::gm::Point pz(0, te::gm::PointZType);
        pz.setX(p->getX());
        pz.setY(p->getY());
        pz.setZ(value);

        mpt.add(dynamic_cast<te::gm::Geometry*>(pz.clone()));
        nsamples++;
      }
    }
    else if (geostype == "Point")
    {
      te::gm::Point *p = dynamic_cast<te::gm::Point*>(gin.get());
      if (atrZ.empty())
        value = p->getZ();
      else
        value = inDset->getDouble(atrZ);

      te::gm::Point pz(0, te::gm::PointZType);
      pz.setX(p->getX());
      pz.setY(p->getY());
      pz.setZ(value);

      mpt.add(dynamic_cast<te::gm::Geometry*>(pz.clone()));
      nsamples++;
    }
    else
      return 0;
  }

  std::unique_ptr<te::gm::Envelope> e(inDset->getExtent(geo_pos));
  env.init((e->getLowerLeftX() - tol), (e->getLowerLeftY() - tol), (e->getUpperRightX() + tol), (e->getUpperRightY() + tol));

  return nsamples;
}

size_t te::mnt::ReadSamples(std::string &inDsetName, te::da::DataSourcePtr &inDsrc, std::string &atrZ, double tol, double max, Simplify alg,
  te::gm::MultiPoint &mpt, te::gm::MultiLineString &isolines, std::string &geostype, te::gm::Envelope &env, int srid)
{

  if (inDsetName.empty())
    return 0;

  std::unique_ptr<te::da::DataSet> inDset;
  size_t nsamples = mpt.getNumGeometries();

  inDset = inDsrc->getDataSet(inDsetName);

  const std::size_t npr = inDset->getNumProperties();

  //Read attributes
  std::vector<std::string>pnames;
  std::vector<int> ptypes;
  for (std::size_t i = 0; i != npr; ++i)
  {
    pnames.push_back(inDset->getPropertyName(i));
    ptypes.push_back(inDset->getPropertyDataType(i));
  }

  te::common::TaskProgress task("Reading Isolines...", te::common::TaskProgress::UNDEFINED, (int)inDset->size());

  std::size_t geo_pos = te::da::GetFirstPropertyPos(inDset.get(), te::dt::GEOMETRY_TYPE);

  inDset->moveBeforeFirst();
  double value = std::numeric_limits<double>::max();

  while (inDset->moveNext())
  {
    if (!task.isActive())
      return false;
    task.pulse();

    std::unique_ptr<te::gm::Geometry> gin(inDset->getGeometry(geo_pos));
    gin->setSRID(srid);
    geostype = gin.get()->getGeometryType();

    if (geostype == "LineString")
    {
      te::gm::LineString *l = dynamic_cast<te::gm::LineString*>(gin.get());
      if (atrZ.empty())
      {
        if (l->is3D())
        value = *l->getZ();
      }
      else
        value = inDset->getDouble(atrZ);

      te::gm::LineString *ls;
      if (alg == Spline)
        ls = te::mnt::SplineInterpolationGrass::pointListSimplify(l, tol, value);
      else if (alg == DouglasPeucker)
        ls = GEOS_DouglasPeucker(l, tol, value);
      else if (alg == None)
      {
        ls = new te::gm::LineString(l->size(), te::gm::LineStringZType, isolines.getSRID());
        for (std::size_t il = 0; il < l->size(); il++)
          ls->setPointZ(il, l->getX(il), l->getY(il), value);
      }
      else
        ls = pointListSimplify(l, tol, max, value);

      isolines.add(dynamic_cast<te::gm::Geometry*>(ls));
      for (std::size_t p = 0; p < ls->size(); p++)
      {
        mpt.add(ls->getPointN(p).release());
      }
      nsamples += ls->size();
    }
    if (geostype == "MultiLineString")
    {
      te::gm::MultiLineString *g = dynamic_cast<te::gm::MultiLineString*>(gin.get());
      std::size_t np = g->getNumGeometries();
      if (!atrZ.empty())
        value = inDset->getDouble(atrZ);
      for (std::size_t i = 0; i < np; ++i)
      {
        te::gm::LineString *l = dynamic_cast<te::gm::LineString*>(g->getGeometryN(i));
        std::unique_ptr<te::gm::LineString> lz(new te::gm::LineString(l->size(), te::gm::LineStringZType, isolines.getSRID()));
        if (atrZ.empty() && l->is3D())
          value = *l->getZ();

        for (std::size_t il = 0; il < l->size(); il++)
          lz->setPointZ(il, l->getX(il), l->getY(il), value);
        l->setSRID(isolines.getSRID());
        std::unique_ptr<te::gm::LineString> ls;
        if (alg == Spline)
          ls.reset(te::mnt::SplineInterpolationGrass::pointListSimplify(l, tol, value));
        else if (alg == DouglasPeucker)
          ls.reset(GEOS_DouglasPeucker(l, tol, value));
        else if (alg == None)
          ls.reset(lz.release());
        else
          ls.reset(pointListSimplify(l, tol, max, value));
        if (ls->size())
        {
          isolines.add(dynamic_cast<te::gm::Geometry*>(ls.get()));
          nsamples += ls->size();

          for (std::size_t p = 0; p < ls->size(); p++)
          {
            mpt.add(ls->getPointN(p).release());
          }
          ls.release();
        }
      }
    }
  }

  std::unique_ptr<te::gm::Envelope> e = inDset->getExtent(geo_pos);
  env.init((e->getLowerLeftX() - tol), (e->getLowerLeftY() - tol), (e->getUpperRightX() + tol), (e->getUpperRightY() + tol));

  return nsamples;
}

double te::mnt::Distance(const te::gm::Coord2D &pt1, const te::gm::Coord2D &pt2)
{
  return(Distance(pt1.getX(), pt1.getY(), pt2.getX(), pt2.getY()));
}

double te::mnt::Distance(const double pt1x, const double pt1y, const double pt2x, const double pt2y)
{
  double  dx,
          dy,
          dist; 

  dx = pt1x - pt2x;
  dx = dx * dx;
  dy = pt1y - pt2y;
  dy = dy * dy;
  dist = sqrt(dy + dx);
  return(dist);
}


te::gm::Point te::mnt::Max(te::gm::Point& p1, te::gm::Point& p2) 
{
  return te::gm::Point(MAX(p1.getX(), p2.getX()), MAX(p1.getY(), p2.getY()));
}

te::gm::Point te::mnt::Min(te::gm::Point& p1, te::gm::Point& p2)
{
  return te::gm::Point(MIN(p1.getX(), p2.getX()), MIN(p1.getY(), p2.getY()));
}


bool te::mnt::Equal(te::gm::Point &p1, te::gm::Point &p2, double &tol)
{
  return (std::fabs(p1.getX() - p2.getX()) < tol && std::fabs(p1.getY() - p2.getY()) < tol && std::fabs(p1.getZ() - p2.getZ()) < tol);
}

bool te::mnt::equalFptSpt(te::gm::Point & fpt, te::gm::Point &spt, double scale)
{
  double  delta = (double)(scale/* * 0.04 / 1000.*/); // point one percent of map tolerance

  if ((fpt.getX() == spt.getX()) && (fpt.getY() == spt.getY())) return true;
  if (fabs(fpt.getX() - spt.getX()) > delta) return false;
  if (fabs(fpt.getY() - spt.getY()) < delta) return true;
  return false;
}

double TePerpendicularDistance(const te::gm::Coord2D& first, const te::gm::Coord2D& last, const te::gm::Coord2D& pin, te::gm::Coord2D &pinter)
{
  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.x = xmin;
  pinter.y = ymin;

  return (sqrt(d12));
}

te::gm::LineString* te::mnt::DouglasPeuckerTA(te::gm::LineString *lineIn, double simpFactor, double Zvalue)
{
  te::gm::LineString *lineOut;

  if (lineIn->size() <= 3)
  {
    lineOut = new te::gm::LineString(*lineIn);
    return lineOut;
  }

  te::gm::LineString lringIn(*lineIn);

  std::vector<te::gm::Point*> lringOut;
  lringOut.push_back(lringIn.getPointN(0).release());

  unsigned int initial = 0;
  unsigned int final = 2;

  while (initial < (lringIn.size() - 2))
  {
    bool distIsGreater = false;

    while (!distIsGreater)
    {
      for (unsigned i = initial + 1; i < final; i++)
      {
        te::gm::Coord2D pInter;
        if (TePerpendicularDistance(te::gm::Coord2D(lringIn.getX(initial), lringIn.getY(initial)), te::gm::Coord2D(lringIn.getX(final), lringIn.getY(final)), 
          te::gm::Coord2D(lringIn.getX(i), lringIn.getY(i)), pInter) > simpFactor)
        {//if distance is greater than maximum distance
          lringOut.push_back(lringIn.getPointN(final - 1).release());

          initial = final - 1;
          final = initial + 2;

          distIsGreater = true;

          break;
        }
      }

      if (!distIsGreater)
      {
        final++;

        if ((final) >= lringIn.size())
        {
          lringOut.push_back(lringIn.getPointN(final - 2).release());

          initial = final;
          distIsGreater = true;

          break;
        }
      }
      else
        break;

      if (final >= lringIn.size())
      {
        initial = final;
        distIsGreater = true;

        break;
      }
    }
  }

  if (lringOut.size() <= 3)
  {
    lineOut = new te::gm::LineString(*lineIn);
  }
  else
  {
    lineOut = new te::gm::LineString(lringOut.size(), te::gm::LineStringZType);
    for (std::size_t i = 0; i < lringOut.size(); i++)
    {
      lineOut->setPointN(i, *lringOut[i]);
      lineOut->setZ(i, Zvalue);
    }
  }

  te::common::FreeContents(lringOut);

  return lineOut;
}

te::gm::LineString* te::mnt::GEOS_DouglasPeucker(te::gm::LineString *ls, double snap, double Zvalue)
{
  std::vector<geos::geom::Coordinate> coordpts;
  for (std::size_t j = 0; j < ls->size(); ++j)
  {
    std::unique_ptr<te::gm::Point> lpt = ls->getPointN(j);
    geos::geom::Coordinate coo(lpt->getX(), lpt->getY());
    coordpts.push_back(coo);
  }
  geos::simplify::DouglasPeuckerLineSimplifier douglas(coordpts);
  douglas.setDistanceTolerance(snap); 
  geos::simplify::DouglasPeuckerLineSimplifier::CoordsVectAutoPtr simplified(douglas.simplify());

  te::gm::LineString* lsout = new te::gm::LineString(simplified->size(), te::gm::LineStringZType);

  for (std::size_t j = 0; j < simplified->size(); ++j)
  {
    geos::geom::Coordinate c1 = simplified.get()->at(j);
    lsout->setPointZ(j, c1.x, c1.y, Zvalue);
  }

  return lsout;
}

te::gm::LineString* te::mnt::pointListSimplify(te::gm::LineString *ls, double snap, double maxdist, double Zvalue)
{
  size_t npts = ls->getNPoints();
  size_t npte = npts;

  te::gm::Coord2D* vxy = ls->getCoordinates();

  //  If line is too short do nothing
  if (npts <= 3 || snap == 0.)
  {
    std::unique_ptr<te::gm::LineString > ptlist(new te::gm::LineString(npts, te::gm::LineStringZType));
    for (size_t i = 0; i < npts; i++)
    {
      ptlist->setPointZ(i, vxy[i].getX(), vxy[i].getY(), Zvalue);
    }
    return ptlist.release();
  }

  //  snap2 = 400 * snap*snap;
  double snap2 = maxdist*maxdist;

  size_t npt;
  bool island = false;
  //  Check for islands before defining number of points to be used
  if (vxy[0] == vxy[npte - 1])
  {
    npt = npte - 1;
    island = true;
  }
  else
    npt = npte;

  //  initialize variables
  size_t i = 0;
  size_t numa = 0;
  size_t numpf = npt - 1;

  //  define anchor
  te::gm::Coord2D axy = vxy[numa];

  //  define floating point
  te::gm::Coord2D pfxy = vxy[numpf];

  bool retv;
  double dmax, d,
    a = 0, b = 0,
    aa1 = 0;
  size_t numdmax;

  while (numa < (npt - 1))
  {
    //    Compute coeficients of straight line y=ax+b
    if (pfxy.getX() == axy.getX())
      retv = true;
    else
    {
      retv = false;
      a = (pfxy.getY() - axy.getY()) / (pfxy.getX() - axy.getX());
      b = pfxy.getY() - a * pfxy.getX();
      aa1 = sqrt(a * a + 1.);
    }

    dmax = 0;
    numdmax = numpf;
    size_t k;

    for (k = numa + 1; k < numpf; k++)
    {
      //      Distance between point and line
      if (retv)
        d = fabs(axy.getX() - vxy[k].getX());
      else
        d = fabs(vxy[k].getY() - a*vxy[k].getX() - b) / aa1;

      if (d > dmax)
      {
        dmax = d;
        numdmax = k;
      }
    }

    if (dmax <= snap)
    {
      //      Store selected point
      if (i > (npt - 1))
      {
        throw;
      }
      vxy[i++] = axy;

      double axbx = pfxy.getX() - axy.getX();
      double ayby = pfxy.getY() - axy.getY();
      double dist2 = axbx*axbx + ayby*ayby;
      if (dist2 > snap2)
      {
        for (k = numpf; k > numa + 1; k--)
        {
          axbx = vxy[k].getX() - axy.getX();
          ayby = vxy[k].getY() - axy.getY();
          dist2 = axbx*axbx + ayby*ayby;
          if (dist2 < snap2)
            break;
        }
        //        Shift anchor
        numa = k;
        axy = vxy[k];
      }
      else
      {
        //        Shift anchor
        numa = numpf;
        axy = vxy[numpf];
      }
      numpf = npt - 1;
    }
    else
      //      Shift floating point
      numpf = numdmax;

    pfxy = vxy[numpf];
  }

  //  Store results
  vxy[i] = vxy[numa];
  npts = i + 1;

  if ((Distance(vxy[i], vxy[i - 1]) < snap) || (Distance(vxy[i], vxy[0]) < snap))
    npts = i;

  size_t newpts;
  if (island)
    newpts = npts + 1;
  else
    newpts = npts;

  std::unique_ptr< te::gm::LineString> ptlist(new te::gm::LineString(newpts, te::gm::LineStringZType, ls->getSRID()));

  for (i = 0; i < npts; i++)
  {
    ptlist->setPointZ(i, vxy[i].getX(), vxy[i].getY(), Zvalue);
  }

  if (island)
    ptlist->setPointZ(npts, vxy[0].getX(), vxy[0].getY(), Zvalue);

  return ptlist.release();
}

/*!
\brief Function that evaluates the distance between a point and a segment
\param fseg is a pointer to the segment first point
\param lseg is a pointer to the segment last point
\param pt is a pointer to the considered point
\param pti is a pointer to the intersection point
\return the distance between the point and the segment
*/
double te::mnt::SegmentDistance(double fx, double fy, double lx, double ly, double ptx, double pty, double *pix, double *piy)
{
  double ux = lx - fx;
  double uy = ly - fy;
  double vx = ptx - fx;
  double vy = pty - fy;
  double uv = (ux*vx + uy*vy);
  if (uv < 0.)
  {
    if (pix && piy)
    {
      *pix = lx;
      *piy = ly;
    }
    return (sqrt(vx*vx + vy*vy));
  }
  else
  {
    ux = fx - lx;
    uy = fy - ly;
    vx = ptx - lx;
    vy = pty - ly;
    uv = (ux*vx + uy*vy);
    if (uv < 0.)
    {
      if (pix && piy)
      {
        *pix = lx;
        *piy = ly;
      }
      return (sqrt(vx*vx + vy*vy));
    }
  }

  double  a, b, c, k, dist, aabb;

  a = ly - fy;
  b = fx - lx;
  c = lx*fy - fx*ly;
  aabb = a*a + b*b;

  dist = fabs((a*ptx + b*pty + c) / (sqrt(aabb)));
  if (pix && piy)
  {
    k = b*ptx - a*pty;
    *pix = (b*k - a*c) / aabb;
    *piy = (-a*k - b*c) / aabb;
  }

  return dist;

}

double te::mnt::coordToSegmentDistance(te::gm::Coord2D &fseg, te::gm::Coord2D &lseg, te::gm::Coord2D &pt, te::gm::Coord2D *pti)
{
  double pix, piy;
  double dist = SegmentDistance(fseg.getX(), fseg.getY(), lseg.getX(), lseg.getY(), pt.getX(), pt.getY(), &pix, &piy);
  if (pti)
  {
    pti->x = pix;
    pti->y = piy;
  }
  return dist;
}

double te::mnt::pointToSegmentDistance(te::gm::Point &fseg, te::gm::Point &lseg, te::gm::Point &pt, te::gm::Point *pti)
{
  double pix, piy;
  double dist = SegmentDistance(fseg.getX(), fseg.getY(), lseg.getX(), lseg.getY(), pt.getX(), pt.getY(), &pix, &piy);
  if (pti)
  {
    pti->setX(pix);
    pti->setY(piy);
  }
  return dist;
}

/*!
\brief Function that finds the center of the vertices of a triangle
\param vert is pointer to the triangle vertices
\param pcx and pcy are pointers to coordinates x and y of the vertex center points
\return TRUE if the centers are evaluated with no errors or FALSE otherwise
*/

short te::mnt::findCenter(te::gm::Point* vert, double* pcx, double* pcy)
{
  double  m1 = 10., m2 = 10.;  // normais aos segmentos 1 e 2 
  double  xm1, ym1,     // coords do ponto medio do segm. 1 
    xm2, ym2,    // coords do ponto medio do segm. 2 
    xaux, yaux,
    px[3], py[3];
  bool  perpvert1 = false,  // perpendicular vertical ao segmento 1 
    perpvert2 = false;  // perpendicular vertical ao segmento 2 
  short  i;

  for (i = 0; i < 3; i++)
  {
    px[i] = vert[i].getX();
    py[i] = vert[i].getY();
  }

  //  calcula o coeficiente angular perpendicular a reta 1
  if ((float)py[0] == (float)py[1])
  {
    if ((float)px[0] == (float)px[1])
    {
      //      pontos sao iguais-> retorna
      return false;
    }
    perpvert1 = true;
  }
  else
    m1 = -(px[1] - px[0]) / (py[1] - py[0]);

  //  calcula o coeficiente angular da perpendicular a reta 2
  if ((float)py[1] == (float)py[2])
  {
    if ((float)px[1] == (float)px[2])
    {
      //      pontos sao iguais-> retorna
      return false;
    }
    perpvert2 = false;
  }
  else
    m2 = -(px[2] - px[1]) / (py[2] - py[1]);

  //  Caso as retas sejam coincidentes, uma circunferencia 
  //   nao eh definida
  if ((float)m1 == (float)m2)
  {
    return false;
  }

  //  passou pelos testes: os pontos definem uma circunferencia
  //  calculo do ponto medio do segmento ponto0-ponto1 (segmento 1)
  xm1 = (px[0] + px[1]) / 2.;
  ym1 = (py[0] + py[1]) / 2.;

  //  calculo do ponto medio do segmento ponto1-ponto2 (segmento 2)
  xm2 = (px[1] + px[2]) / 2.;
  ym2 = (py[1] + py[2]) / 2.;

  //  calculo das coordenadas do centro: ponto de interseccao das mediatrizes
  //   dos segmentos ponto0-ponto1 e ponto1-ponto2
  if (perpvert1 == true)
  {
    xaux = xm1;
    yaux = (ym2 + m2 * (xaux - xm2));
  }
  else if (perpvert2 == true)
  {
    xaux = xm2;
    yaux = (ym1 + m1 * (xaux - xm1));
  }
  else
  {
    xaux = (m1*xm1 - m2*xm2 - ym1 + ym2) / (m1 - m2);
    yaux = (ym1 + m1 * (xaux - xm1));
  }
  *pcx = xaux;
  *pcy = yaux;

  return true;
}
/*!
\brief Function for check segment intersection
\param pfr, pto, lfr, lto are pointers to Point objects
\return point of intersection
*/
short te::mnt::segIntersect(te::gm::Point &pfr, te::gm::Point &pto, te::gm::Point &lfr, te::gm::Point &lto)
{
  te::gm::Point  pt(0, te::gm::PointZType);

  return segInterPoint(pfr, pto, lfr, lto, &pt);
}

/*!
\brief Function for determines the point of two segment intersections
\param pfr, pto, lfr, lto and pt are pointers to Point objects
\return TRUE if the intersection was calculated with no errors or FALSE otherwise
*/

bool te::mnt::segInterPoint(te::gm::Point &pfr, te::gm::Point &pto, te::gm::Point &lfr, te::gm::Point &lto, te::gm::Point *pt)
{
  //  Adapted from TWO-DIMENSIONAL CLIPPING: A VECTOR-BASED APPROACH
  //  Hans J.W. Spoelder, Fons H. Ullings in:
  ////  Graphics Gems I, pp.701, 

  double  a, b, c,
    px, py, lx, ly, lpx, lpy, s;

  px = pto.getX() - pfr.getX();
  py = pto.getY() - pfr.getY();
  lx = lto.getX() - lfr.getX();
  ly = lto.getY() - lfr.getY();
  lpx = pfr.getX() - lfr.getX();
  lpy = pfr.getY() - lfr.getY();

  a = py * lx - px * ly;
  b = lpx * ly - lpy * lx;
  c = lpx * py - lpy * px;

  if (a == 0) // Linhas paralelas
    return false;
  else
  {
    if (a > 0)
    {
      if ((b < 0) || (b > a) || (c < 0) || (c > a))
        return false;
    }
    else
    {
      if ((b > 0) || (b < a) || (c > 0) || (c < a))
        return false;
    }
    if (pt)
    {
      s = b / a;
      pt->setX(pfr.getX() + (px*s));
      pt->setY(pfr.getY() + (py*s));
    }
  }
  return true;
}


/*!
\brief Function that test the vertex values
\param isolin is the isolin value of a given isolin
\param p3da is a pointer to a vector with Point3d objects
\return TRUE if the vertex values are in the same isolin or FALSE otherwise
*/
bool te::mnt::testVertexValues(double isolin, te::gm::Point *p3da)
{
  double  delta = .9999;
  size_t  i;

  for (i = 0; i < 3; i++)
    if (p3da[i].getZ() == isolin)
      p3da[i].setZ(isolin*delta);


  if ((p3da[0].getZ() == std::numeric_limits<double>::max()) ||
    (p3da[1].getZ() == std::numeric_limits<double>::max()) ||
    (p3da[2].getZ() == std::numeric_limits<double>::max()))
    return false;

  if ((p3da[0].getZ() < isolin) &&
    (p3da[1].getZ() < isolin) &&
    (p3da[2].getZ() < isolin))
    return false;

  if ((p3da[0].getZ() > isolin) &&
    (p3da[1].getZ() > isolin) &&
    (p3da[2].getZ() > isolin))
    return false;

  for (i = 0; i < 3; i++)
  {
    if ((p3da[i].getZ() == isolin) &&
      (p3da[(i + 1) % 3].getZ() > isolin) &&
      (p3da[(i + 2) % 3].getZ() > isolin))
      return false;
    if ((p3da[i].getZ() == isolin) &&
      (p3da[(i + 1) % 3].getZ() < isolin) &&
      (p3da[(i + 2) % 3].getZ() < isolin))
      return false;
  }

  return true;
}


/*!
\brief Function that defines the intersections between an isoline and a triangle
\param ptline is a pointer to a Line object
\param p3da is a pointer to a Point3d object
\return TRUE if the intersections are defined or FALSE otherwise
*/

bool te::mnt::defineInterTriangle(std::vector<te::gm::Point> &ptline, te::gm::Point *p3da, double isolin)
{
  size_t contrl = 0;

  // Verify if contour intercepts triangle
  if ((isolin < p3da[0].getZ()) && (isolin < p3da[1].getZ()) &&
    (isolin < p3da[2].getZ()))
    return false;

  if ((isolin > p3da[0].getZ()) && (isolin > p3da[1].getZ()) &&
    (isolin > p3da[2].getZ()))
    return false;

  if (((isolin <= p3da[0].getZ()) && (isolin > p3da[1].getZ())) ||
    ((isolin > p3da[0].getZ()) && (isolin <= p3da[1].getZ())))
  {// Calculate intersection with edge 0
    if (!defineInterEdge(ptline, p3da[0], p3da[1], isolin))
      return false;
    contrl++;
  }

  if (((isolin <= p3da[1].getZ()) && (isolin > p3da[2].getZ())) ||
    ((isolin > p3da[1].getZ()) && (isolin <= p3da[2].getZ())))
  {// Calculate intersection with edge 0
    if (!defineInterEdge(ptline, p3da[1], p3da[2], isolin))
      return false;
    contrl++;
  }

  if (((isolin <= p3da[2].getZ()) && (isolin > p3da[0].getZ())) ||
    ((isolin > p3da[2].getZ()) && (isolin <= p3da[0].getZ())))
  {// Calculate intersection with edge 0
    if (!defineInterEdge(ptline, p3da[2], p3da[0], isolin))
      return false;
    contrl++;
  }
  if (contrl != 2)
    return false;

  return true;
}


/*!
\brief Function that defines the intersection with an edge of a triangle
\param ptline is a pointer to a Line object
\param fpt and spt are the first and second point of the edge
\return TRUE if the intersection is defined or FALSE otherwise
*/

bool te::mnt::defineInterEdge(std::vector<te::gm::Point> &ptline, te::gm::Point &fpt, te::gm::Point &spt, double isolin)
{
  te::gm::Point dvect(0, te::gm::PointZType); // Vector from fpt to spt
  te::gm::Point pt(0, te::gm::PointZType);
  double  paux;// Parametric variable

  if (spt.getZ() == fpt.getZ())
    return false;

  dvect.setX(spt.getX() - fpt.getX());
  dvect.setY(spt.getY() - fpt.getY());
  dvect.setZ(spt.getZ() - fpt.getZ());

  paux = (isolin - fpt.getZ()) / dvect.getZ(); // t = (Z - Z0) / C

  pt.setX((fpt.getX() + dvect.getX() * paux)); // X = X0 + A * t
  pt.setY((fpt.getY() + dvect.getY() * paux)); // Y = Y0 + B * t
  pt.setZ(isolin);

  ptline.push_back(pt);
  return true;
}

bool te::mnt::extractLines(std::vector<te::gm::Point> &pline, std::vector< std::unique_ptr< te::gm::LineString > > &clinlist, double scale)
{
  std::vector<te::gm::Point> velin;
  te::gm::Point pt(0, te::gm::PointZType);

  std::vector<te::gm::Point>::iterator it, itf, itn;
  short borda; // borda = 1 no boundary reached, 2 one boundary reached.

 
  initLineVector(pline, velin);
  borda = 1;

  //index pindex to the current point of old line
  it = itf = pline.begin();
  itn = ++it;
  te::gm::Point ptf(*itf);
  te::gm::Point ptn(*itn);
  while (it != pline.end())
  {
    ptf = *itf;
    ptn = *itn;
    pt = velin[velin.size()-1];
    if (equalFptSpt(ptf, pt, scale))
    {// conection on ptf
      velin.push_back(ptn);
      --it;
      it = itf = pline.erase(it);
      itn = pline.erase(it);
    }
    else if (equalFptSpt(ptn, pt, scale))
    {// conection on ptn
      velin.push_back(ptf);
      --it;
      it = itf = pline.erase(it);
      itn = pline.erase(it);
    }
    else
    {
      itf = ++it;;
      if (it != pline.end())
        itn = ++it;
      if (itf == pline.end() || itn == pline.end())
      {// There are no points to match
        if (borda == 2)
        {// If first point is on boundary -> open contour
          assembLine(clinlist, velin);
          initLineVector(pline, velin); //SSL0296
          borda = 1;
        }
        else
        {// If not on boundary, invert line points
          std::reverse(velin.begin(), velin.end());
          borda = 2;
        }
        it = itf = pline.begin();
        if (it == pline.end())
          break;
        itn = ++it;
      }
      continue;
    }

    ptf = velin[0];
    ptn = velin[velin.size()-1];
    if (Equal(ptf, ptn, scale))
    {// If closed contour
      assembLine(clinlist, velin);
      initLineVector(pline, velin);
    }
    it= itf = pline.begin();
    if (it == pline.end())
      break;
    itn = ++it;
  }

  assembLine(clinlist, velin);

  return true;
}

/*************************************************************
NAME    : initLineVector
AUTHOR    : Laercio Namikawa      aug-95
RESUME    : Initialize Line Vector with points of a line.
INPUT    : pline    line
OUTPUT    : vect    line with first and second points
RETURN    : TRUE    successful
FALSE    failure
*************************************************************/
bool te::mnt::initLineVector(std::vector<te::gm::Point> &pline, std::vector<te::gm::Point> &vect)
{
  short  c;
  std::vector<te::gm::Point>::iterator it;

  if (vect.size() != 0)
    vect.clear();

  if (pline.size() < 2)
    return true;

  for (c = 0; c < 2; c++)
  {
    vect.push_back(*pline.begin());
    pline.erase(pline.begin());
  }

  if (vect.size() != 2)
    return false;
  return true;
}


/*************************************************************
NAME    : assembLine
AUTHOR    : Laercio Namikawa      aug-95
RESUME    : Assemble a line with Line points and
save in LineList.
INPUT    : vect    line with points
OUTPUT    : linlout  output linelist
RETURN    : TRUE    successful
FALSE    failure
*************************************************************/
bool te::mnt::assembLine(std::vector< std::unique_ptr< te::gm::LineString > > &linlout, std::vector<te::gm::Point> &vect)
{
  if (vect.size() <= 1)
    return true;

  std::unique_ptr< te::gm::LineString > linout(new te::gm::LineString(vect.size(), te::gm::LineStringZType));

  for (size_t i = 0; i < vect.size(); i++)
    linout->setPointN(i, vect[i]);

  linlout.push_back(std::move(linout));

  return true;
}

/*!
\brief Function that defines the triangle normal vector
\param p3da are pointers to objets Point3d representing the triangle vertices
\param nvector are the coordinates of the normal vector
\return TRUE always
*/

bool te::mnt::triangleNormalVector(te::gm::Point *p3da, double *nvector)
{
  double  ux, vx, uy, vy, uz, vz;

  if ((p3da == NULL) || (nvector == NULL)){
    return false;
  }

  nvector[0] = 0.;
  nvector[1] = 0.;

  // Define normal vector (uvx,uvy,uvz)
  ux = p3da[1].getX() - p3da[0].getX();
  vx = p3da[2].getX() - p3da[0].getX();
  uy = p3da[1].getY() - p3da[0].getY();
  vy = p3da[2].getY() - p3da[0].getY();
  uz = p3da[1].getZ() - p3da[0].getZ();
  vz = p3da[2].getZ() - p3da[0].getZ();

  if ((ux == 0) && (vx == 0))
  {
    nvector[0] = 1.;
    nvector[1] = 0.;
    nvector[2] = 0.;
    return true;
  }
  if ((uy == 0) && (vy == 0))
  {
    nvector[0] = 0.;
    nvector[1] = 1.;
    nvector[2] = 0.;
    return true;
  }
  if ((uz == 0) && (vz == 0))
  {
    nvector[0] = 0.;
    nvector[1] = 0.;
    nvector[2] = 1.;
    return true;
  }

  nvector[2] = ux * vy - vx * uy;
  if (nvector[2] < 0.)
  {
    // Make sure that normal vector is always positive
    nvector[2] = -nvector[2];
    nvector[0] = vy * uz - uy * vz;
    nvector[1] = ux * vz - vx * uz;
  }
  else
  {
    nvector[0] = uy * vz - vy * uz;
    nvector[1] = vx * uz - ux * vz;
  }
  return true;
}

/*!
\brief Function that normalize a vector by its size
\param nvector is a pointer to the vector coordinates x, y and z
\return TRUE always
*/

bool te::mnt::normalizeVector(double *nvector)
{
  double  vectorSize;

  if (nvector == NULL){
    return false;
  }

  vectorSize = sqrt(nvector[0] * nvector[0] + nvector[1] * nvector[1] + nvector[2] * nvector[2]);

  if (vectorSize != 0.)
  {
    nvector[0] = nvector[0] / vectorSize;
    nvector[1] = nvector[1] / vectorSize;
    nvector[2] = nvector[2] / vectorSize;
  }
  return true;
}

/*!
\brief Function that checks if a point pt is on a segment
\param pt, fseg, lseg are pointers to Point objects
\param tol is the tolerance distance between the point and the segment
\return true if the point is on the segment or false otherwise
*/
bool te::mnt::onSegment(te::gm::Point& pt, te::gm::Point& fseg, te::gm::Point& lseg, double tol)
{
  double  pax, pay, bax, bay;
  double  area2, minxy, daux;

  if ((pt.getX() < fseg.getX()) && (pt.getX() < lseg.getX()))
    return false;
  if ((pt.getY() < fseg.getY()) && (pt.getY() < lseg.getY()))
    return false;
  if ((pt.getX() > fseg.getX()) && (pt.getX() > lseg.getX()))
    return false;
  if ((pt.getY() > fseg.getY()) && (pt.getY() > lseg.getY()))
    return false;
  pax = pt.getX() - fseg.getX();
  bax = lseg.getX() - fseg.getX();
  if ((pax == 0.) && (bax == 0.))
    // On segment
    return true;

  pay = pt.getY() - fseg.getY();
  bay = lseg.getY() - fseg.getY();
  if ((pay == 0.) && (bay == 0.))
    // On segment
    return true;

  area2 = fabs((pay * bax) - (pax * bay));
  bax = fabs(bax);
  bay = fabs(bay);
  minxy = bay;
  if (bax < bay)
    minxy = bax;
  daux = area2 / (bax + bay - (minxy / 2));

  if (daux < tol)
    // On segment
    return true;

  return false;
}


int te::mnt::onSameSide(te::gm::Coord2D pt1, te::gm::Coord2D pt2, te::gm::Coord2D fseg, te::gm::Coord2D lseg)
{
  double  a, b, c, ip, ipt;

  a = lseg.getY() - fseg.getY();
  b = fseg.getX() - lseg.getX();
  c = lseg.getX()*fseg.getY() - fseg.getX()*lseg.getY();
  ip = a*pt1.getX() + b*pt1.getY() + c;
  if (ip == 0.) // On segment
    return -1;
  ipt = a*pt2.getX() + b*pt2.getY() + c;
  if ((ip > 0.) && (ipt < 0.))
    return 0;
  if ((ip < 0.) && (ipt > 0.))
    return 0;

  // On same side
  return 1;
}

/*!
\brief Function that filters the points that are closer considering a given tolerance
\param p3dl is a pointer to a list of Point3dList objects
\param tol is a tolerance value
\return TRUE always
*/

bool te::mnt::point3dListFilter(std::vector<te::gm::Point> &p3dl, std::vector<bool> &fixed, double tol)
{
  te::gm::Point p3d(0, te::gm::PointZType);
  int32_t npts, i, j, maxdiffindex = 0;
  double x0 = 0, y0 = 0, coef[5];
  double maxdiff;
  double maxfixptdiff;
  short degree;

  npts = 0;
 
  try
  {
    double *vectd = new double[200];
    double *vectz = new double[200];
    double *vectx = new double[200];
    double *vecty = new double[200];
    double *fvectz = new double[200];
    short *ptype = new short[200];

    te::gm::Point *p3dlaux[200];
    for(i = 0; i < 200; ++i)
      p3dlaux[i] = new te::gm::Point(0, te::gm::PointZType);

    npts = 0;
    for (size_t pp = 0; pp < p3dl.size(); pp++)
    {
      p3d = p3dl[pp];
      if ((p3d.getX() >= std::numeric_limits< float >::max()) || (npts > 199))
      {
        // If last point of a line
        maxdiff = std::numeric_limits< float >::max();
        maxfixptdiff = std::numeric_limits< float >::max();
        j = 4;
        while ((maxdiff > tol) && (maxfixptdiff > tol) &&
          (j > 3) && (npts > 3))
        {
          j = 0;
          for (i = 0; i < npts; i++)
          {
            if (ptype[i] == 0)
              continue;
            vectx[j] = vectd[i];
            vecty[j++] = vectz[i];
          }
          if (j < 3)
            break;
          if (j < 5)
            degree = 1;
          else if (j < 8)
            degree = 2; //6 is minimum npts for 3rd degree
          else
            degree = 3; //8 is minimum npts for 3rd degree
          Least_square_fitting(vectx, vecty, (short)j, degree, coef);
          for (i = 0; i < npts; i++)
          {
            fvectz[i] = coef[0];
            for (j = 1; j <= degree; j++)
              fvectz[i] += coef[j] * pow(vectd[i], (double)j);
          }
          maxfixptdiff = 0.;
          for (i = 0; i < npts; i++)
          {
            if ((ptype[i] == 0) || (ptype[i] == 1))
              continue;
            if (fabs(fvectz[i] - vectz[i]) > maxfixptdiff)
              maxfixptdiff = (float)fabs(fvectz[i] - vectz[i]);
          }
          maxdiff = 0.;
          for (i = 0; i < npts; i++)
          {
            if ((ptype[i] == 0) || (ptype[i] == 2))
              continue;
            if (fabs(fvectz[i] - vectz[i]) > maxdiff)
            {
              maxdiff = (float)fabs(fvectz[i] - vectz[i]);
              maxdiffindex = i;
            }
          }
          if (i == 0)
            break;
          ptype[maxdiffindex] = 0;
        }

        for (i = 0; i < npts; i++)
        {
          p3d = *p3dlaux[i];
          p3d.setZ(fvectz[i]);
        }
        npts = 0;
        continue;
      }
      if (npts == 0)
      {
        x0 = p3d.getX();
        y0 = p3d.getY();
      }

      p3dlaux[npts]->setX(p3d.getX());
      p3dlaux[npts]->setY(p3d.getY());

      vectd[npts] = sqrt((p3d.getX() - x0)*(p3d.getX() - x0) +
        (p3d.getY() - y0)*(p3d.getY() - y0));
      if (npts > 0)
        vectd[npts] += vectd[npts - 1];
      if (p3d.getZ() < std::numeric_limits< float >::max())
        ptype[npts] = 1;
      else
        ptype[npts] = 0;
      if (!fixed[pp])
        ptype[npts] = 2;
      vectz[npts] = p3d.getZ();
      fvectz[npts++] = p3d.getZ();
      x0 = p3d.getX();
      y0 = p3d.getY();
    }

    delete[] vectd;
    delete[] vectz;
    delete[] vectx;
    delete[] vecty;
    delete[] fvectz;
    delete[] ptype;

    for(size_t ii = 0; ii < 200; ++ii)
      delete p3dlaux[ii];

  }
  catch (std::bad_alloc& ba)
  {
    std::cerr << "bad_alloc caught: " << ba.what() << '\n';
    return false;
  }

  return true;
}

/*!
\brief Function that performs a leat square fitting in a set of points
\param vectx and vecty are the point coordinates
\param np is the number of points
\param deg is degree of the polynomium used in the fitting
\param coef is a pointer to the the polynomiun coefficients
\return TRUE if the fitting was performed with no errors or FALSE otherwise
*/

bool te::mnt::Least_square_fitting(double *vectx, double *vecty, short np, short deg, double *coef)
{
  short  i, j, m, n;
  double matx[6][6];
  double *powx,
    *sumpow,
    *fx;

  if (np > 200)
    return false;

  try
  {
    powx = new double[(unsigned)np];
    sumpow = new double[200];
    fx = new double[(unsigned)np];
  }
  catch (std::bad_alloc& ba)
  {
    std::cerr << "bad_alloc caught: " << ba.what() << '\n';
    return false;
  }

  // initialization 
  for (i = 0; i<np; i++)
  {
    powx[i] = 1.;
    fx[i] = vecty[i];
  }

  // compute the rigth hand side of the first normal equation 
  matx[0][deg + 1] = 0;
  for (i = 0; i<np; i++)
    matx[0][deg + 1] = matx[0][deg + 1] + vecty[i];

  // compute the rigth hand side of the other normal equation 
  // and also the sums  
  sumpow[0] = np;
  for (i = 1; i <= deg; i++)
  {
    sumpow[i] = 0;
    matx[i][deg + 1] = 0;

    for (j = 0; j< np; j++)
    {
      powx[j] = powx[j] * vectx[j];
      fx[j] = fx[j] * vectx[j];
      sumpow[i] = sumpow[i] + powx[j];
      matx[i][deg + 1] = matx[i][deg + 1] + fx[j];
    }
  }

  for (i = deg + 1; i< 2 * (deg + 1); i++)
  {
    sumpow[i] = 0;
    for (j = 0; j<np; j++)
    {
      powx[j] = powx[j] * vectx[j];
      sumpow[i] = sumpow[i] + powx[j];
    }
  }

  for (i = 0; i <= deg; i++)
    for (j = 0; j <= deg; j++)
      matx[i][j] = sumpow[i + j];

  // perform Gaussian elimination
  m = deg + 1;
  n = 1;

  delete[] powx;
  delete[] sumpow;
  delete[] fx;

  if (!Gauss_elimination(m, n, matx))
  {
    return false;
  }

  for (i = 0; i<deg + 1; i++)
    coef[i] = matx[i][m];

  return false;
}


/*********************************************************************
NAME    : Gauss_elimination
AUTHOR    : Jose Claudio Mura     oct - 95
RESUME    : Solve a system of linear equation throught the Gauss
elimination method.
INPUT         : Number of the variables (m), number of equations (n) and
array with the variables (mat).
OUTPUT        : Solution in the last column of the (a).
RETURN    : TRUE    successful
FALSE   failure
*********************************************************************/

bool te::mnt::Gauss_elimination(short m, short n, double mat[6][6])
{
  short i, j, k, l, ii, jj, kk;
  double aux, piv, eps;

  /*-- estimated error bound (machine epsilon) --*/
  eps = .000000001;

  if (m > 1)
  {
    for (i = 0; i<m - 1; i++)
    {
      piv = fabs(mat[i][i]);
      ii = i + 1;
      jj = i;

      /*-- search for the index in of maximum pivot value -*/
      for (j = ii; j< m; j++)
      {
        if (fabs(mat[j][i]) > piv)
        {
          piv = fabs(mat[j][i]);
          jj = j;
        }
      }

      if (i != jj)
      {
        /*-- interchange rows i and index jj --*/
        for (j = i; j< m + n; j++)
        {
          aux = mat[i][j];
          mat[i][j] = mat[jj][j];
          mat[jj][j] = aux;
        }
      }

      /*-- check if pivot is too smal --*/
      if (piv < eps)
      {
        return false;
      }

      /*-- forward elimination step --*/
      for (j = ii; j< m; j++)
      {
        for (k = ii; k< m + n; k++)
          mat[j][k] = mat[j][k] - mat[j][i] * mat[i][k] / mat[i][i];
      }
    }

    if (fabs(mat[m - 1][m - 1]) < eps)
    {
      return false;
    }

    /*-- back substitution --*/
    for (i = 0; i<n; i++)
    {
      mat[m - 1][i + m] = mat[m - 1][i + m] / mat[m - 1][m - 1];

      for (j = 0; j<m - 1; j++)
      {
        k = m - j - 2;
        kk = k + 1;

        for (l = kk; l<m; l++)
          mat[k][i + m] = mat[k][i + m] - mat[l][i + m] * mat[k][l];

        mat[k][i + m] = mat[k][i + m] / mat[k][k];
      }
    }
  }
  else
  {
    if (fabs(mat[0][0]) < eps)
    {
      return false;
    }

    for (i = 0; i<n; i++)
      mat[0][m + i] = mat[0][m + i] / mat[0][0];
  }
  
  return true;
}


bool te::mnt::SaveIso(std::string& outDsetName, te::da::DataSourcePtr &outDsrc, std::vector< std::unique_ptr< te::gm::LineString > > &isolist, std::vector<double> &guidevalues, int srid)
{
  std::unique_ptr<te::da::DataSetType> dt(new te::da::DataSetType(outDsetName));

  //Primary key
  std::unique_ptr<te::dt::SimpleProperty> prop0(new te::dt::SimpleProperty("FID", te::dt::INT32_TYPE));
  prop0->setAutoNumber(true);
  std::unique_ptr<te::dt::SimpleProperty> prop1 (new te::dt::SimpleProperty("Z", te::dt::DOUBLE_TYPE));
  std::unique_ptr<te::dt::SimpleProperty> prop11(new te::dt::SimpleProperty("type", te::dt::STRING_TYPE));
  std::unique_ptr<te::gm::GeometryProperty> prop2(new te::gm::GeometryProperty("iso", 0, te::gm::LineStringZType, true));
  prop2->setSRID(srid);
  dt->add(prop0.release());
  dt->add(prop1.release());
  dt->add(prop11.release());
  dt->add(prop2.release());

  std::unique_ptr<te::mem::DataSet> ds(new te::mem::DataSet(dt.get()));

  int id = 0;

  te::common::TaskProgress task("Saving Isolines...", te::common::TaskProgress::UNDEFINED, (int)isolist.size());

  for (unsigned int Idx = 0; Idx < isolist.size(); ++Idx)
  {
    if (!task.isActive())
      return false;
    task.pulse();

    te::mem::DataSetItem* dataSetItem =new te::mem::DataSetItem(ds.get());
    isolist[Idx]->setSRID(srid);
    double *zvalue = isolist[Idx]->getZ();
    dataSetItem->setInt32(0, id++);
    if (zvalue){
      dataSetItem->setDouble(1, zvalue[0]);
      if (std::find(guidevalues.begin(), guidevalues.end(), zvalue[0]) != guidevalues.end())
        dataSetItem->setString(2, "GUIDELINE");
      else
        dataSetItem->setString(2, "NORMAL");
    }
    dataSetItem->setGeometry(3, isolist[Idx].release());

    ds->add(dataSetItem);
  }
  te::vp::Save(outDsrc.get(), ds.get(), dt.get());

  return true;

}

bool te::mnt::convertAngleToPlanar(double& val, te::common::UnitOfMeasurePtr planar, int type)
{
  if (type == 0) //lenght
  {
    switch (planar->getId())
    {
    case te::common::UOM_Metre:
      val *= 111000.;            // 1 degree = 111.000 meters
      break;
    case te::common::UOM_Kilometre:
      val *= 111.;                  // 1 degree = 111 kilometers
      break;
    case te::common::UOM_Foot:
      val *= 364173.24;        //  1 feet  = 3.28084 meters
      break;
    default:
      return false;
    }
  }

  if (type == 1) //area
  {
    switch (planar->getId())
    {
    case te::common::UOM_Metre:
      val *= 111000. * 111000.;            // 1 degree = 111.000 meters
      break;
    case te::common::UOM_Kilometre:
      val *= 111. * 111.;                  // 1 degree = 111 kilometers
      break;
    case te::common::UOM_Foot:
      val *= 364173.24 * 364173.24;        //  1 feet  = 3.28084 meters
      break;
    default:
      return false;
    }
  }
  return true;
}


bool te::mnt::convertPlanarToAngle(double& val, te::common::UnitOfMeasurePtr unit)
{
  switch (unit->getId())
  {
  case te::common::UOM_Metre:
    val /= 111000;            // 1 degree = 111.000 meters
    break;
  case te::common::UOM_Kilometre:
    val /= 111;               // 1 degree = 111 kilometers
    break;
  case te::common::UOM_Foot:
    val /= 364173.24;        //  1 feet  = 3.28084 meters
    break;
  default:
    return false;
  }
  return true;
}

void te::mnt::getMinMax(te::rst::Raster *inputRst, double &vmin, double &vmax)
{
  double min = std::numeric_limits<double>::max();
  double max = std::numeric_limits<double>::min();

  std::complex<double> pixel;
  unsigned int rf = inputRst->getNumberOfRows() - 1;
  unsigned int cf = inputRst->getNumberOfColumns() - 1;
  double no_data = inputRst->getBand(0)->getProperty()->m_noDataValue;

  for (unsigned r = 0; r <= rf; r++)
    for (unsigned c = 0; c <= cf; c++)
    {
      inputRst->getValue(c, r, pixel);

      if (pixel.real() == no_data)
        continue;

      if (pixel.real() < min)
        min = pixel.real();

      if (pixel.real() > max)
        max = pixel.real();
    }

  vmin = min;
  vmax = max;
}

double te::mnt::PerpendicularDistance(te::gm::Coord2D& first, te::gm::Coord2D& last, te::gm::Coord2D& pin, te::gm::Coord2D& pinter)
{
  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.x = xmin;
  pinter.y = ymin;
  return (sqrt(d12));
}


te::mnt::mntType te::mnt::getMNTType(const te::da::DataSetType* dt)
{
  if (dt->hasGeom())
  {
    te::gm::GeomType gmType = te::da::GetFirstGeomProperty(dt)->getGeometryType();
    if (gmType == te::gm::PointType || gmType == te::gm::MultiPointType ||
      gmType == te::gm::PointZType || gmType == te::gm::MultiPointZType ||
      gmType == te::gm::PointMType || gmType == te::gm::MultiPointMType ||
      gmType == te::gm::PointZMType || gmType == te::gm::MultiPointZMType) //samples
      return SAMPLE;

    if (gmType == te::gm::LineStringType || gmType == te::gm::MultiLineStringType ||
      gmType == te::gm::LineStringZType || gmType == te::gm::MultiLineStringZType ||
      gmType == te::gm::LineStringMType || gmType == te::gm::MultiLineStringMType ||
      gmType == te::gm::LineStringZMType || gmType == te::gm::MultiLineStringZMType) //isolines
      return ISOLINE;

    if (gmType == te::gm::TINZType || gmType == te::gm::MultiPolygonZType || gmType == te::gm::PolyhedralSurfaceZType || gmType == te::gm::PolygonZType ||
      gmType == te::gm::GeometryType)//TIN
      return TIN;

  }
  if (dt->hasRaster()) //GRID
  {
    return GRID;
  }

  return OTHER;

}