Delaunay2D.cxx

Go to the documentation of this file.

// @(#)root/mathcore:$Id: Delaunay2D.h,v 1.00
// Author: Daniel Funke, Lorenzo Moneta
 
/*************************************************************************
 * Copyright (C) 2015 ROOT Math Team                                     *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/
 
// Implementation file for class Delaunay2D
 
#include "Math/Delaunay2D.h"
#include "Rtypes.h"
 
//#include <thread>
 
// use the triangle library if we do not use CGAL
#ifndef HAS_CGAL
#include "triangle.h"
#endif
 
#include <algorithm>
#include <stdlib.h>
 
#include <iostream>
#include <limits>
 
namespace ROOT {
 
   namespace Math {
 
 
/// class constructor from array of data points
Delaunay2D::Delaunay2D(int n, const double * x, const double * y, const double * z,
                       double xmin, double xmax, double ymin, double ymax)
{
   // Delaunay2D normal constructor
 
   fX            = x;
   fY            = y;
   fZ            = z;
   fZout         = 0.;
   fNpoints      = n;
   fOffsetX      = 0;
   fOffsetY      = 0;
   fScaleFactorX = 0;
   fScaleFactorY = 0;
   fNdt          = 0;
   fXNmax        = 0;
   fXNmin        = 0;
   fYNmin        = 0;
   fYNmax        = 0;
 
   SetInputPoints(n,x,y,z,xmin,xmax,ymin,ymax);
 
}
 
/// set the input points
void Delaunay2D::SetInputPoints(int n, const double * x, const double * y, const double * z,
                           double xmin, double xmax, double ymin, double ymax) {
 
 
   fInit         = kFALSE;
 
   if (n == 0 || !x || !y || !z ) return;
 
   if (xmin >= xmax) {
      xmin = *std::min_element(x, x+n);
      xmax = *std::max_element(x, x+n);
 
      ymin = *std::min_element(y, y+n);
      ymax = *std::max_element(y, y+n);
   }
 
   fOffsetX      = -(xmax+xmin)/2.;
   fOffsetY      = -(ymax+ymin)/2.;
 
   if ( (xmax-xmin) != 0 ) {
      fScaleFactorX = 1./(xmax-xmin);
      fXNmax        = Linear_transform(xmax, fOffsetX, fScaleFactorX); //xTransformer(xmax);
      fXNmin        = Linear_transform(xmin, fOffsetX, fScaleFactorX); //xTransformer(xmin);
   } else {
      // we can't find triangle in this case
      fInit = true;
   }
 
   if (ymax-ymin != 0) {
      fScaleFactorY = 1./(ymax-ymin);
      fYNmax        = Linear_transform(ymax, fOffsetY, fScaleFactorY); //yTransformer(ymax);
      fYNmin        = Linear_transform(ymin, fOffsetY, fScaleFactorY); //yTransformer(ymin);
   } else {
      fInit = true;
   }
   //printf("Normalized space extends from (%f,%f) to (%f,%f)\n", fXNmin, fYNmin, fXNmax, fYNmax);
 
 
#ifndef HAS_CGAL
   fXCellStep    = 0.;
   fYCellStep    = 0.;
#endif
}
 
//______________________________________________________________________________
double Delaunay2D::Interpolate(double x, double y)
{
   // Return the interpolated z value corresponding to the given (x,y) point
 
   // Initialise the Delaunay algorithm if needed.
   // CreateTrianglesDataStructure computes fXoffset, fYoffset,
   // fXScaleFactor and fYScaleFactor;
   // needed in this function.
   FindAllTriangles();
 
   // handle case no triangles are found  return default value (i.e. 0)
   // to do: if points are aligned in one direction we could return in
   // some case the 1d interpolated value
   if (fNdt == 0) {
      return fZout;
   }
 
   // Find the z value corresponding to the point (x,y).
   double xx, yy;
   xx = Linear_transform(x, fOffsetX, fScaleFactorX); //xx = xTransformer(x);
   yy = Linear_transform(y, fOffsetY, fScaleFactorY); //yy = yTransformer(y);
   double zz = DoInterpolateNormalized(xx, yy);
 
   // the case of points on a regular grid (i.e. points on triangle edges) it is now handles in
   // DoInterpolateNormalized
 
   return zz;
}
 
//______________________________________________________________________________
void Delaunay2D::FindAllTriangles()
{
 
   if (fInit)
      return;
   else
      fInit = kTRUE;
 
   // Function used internally only. It creates the data structures needed to
   // compute the Delaunay triangles.
 
   // Offset fX and fY so they average zero, and scale so the average
   // of the X and Y ranges is one. The normalized version of fX and fY used
   // in Interpolate.
 
   DoNormalizePoints(); // call backend specific point normalization
 
   DoFindTriangles(); // call backend specific triangle finding
 
   fNdt = fTriangles.size();
}
 
// backend specific implementations
 
#ifndef HAS_CGAL
 
// Triangle implementation (default case)
 
/// Triangle implementation for points normalization
void Delaunay2D::DoNormalizePoints() {
   for (Int_t n = 0; n < fNpoints; n++) {
      fXN.push_back(Linear_transform(fX[n], fOffsetX, fScaleFactorX));
      fYN.push_back(Linear_transform(fY[n], fOffsetY, fScaleFactorY));
   }
 
   //also initialize fXCellStep and FYCellStep
   fXCellStep = fNCells / (fXNmax - fXNmin);
   fYCellStep = fNCells / (fYNmax - fYNmin);
}
 
/// Triangle implementation for finding all the triangles
void Delaunay2D::DoFindTriangles() {
 
   auto initStruct = [] (triangulateio & s) {
      s.pointlist = nullptr;              /* In / out */
      s.pointattributelist = nullptr;     /* In / out */
      s.pointmarkerlist = nullptr;        /* In / out */
      s.numberofpoints = 0;               /* In / out */
      s.numberofpointattributes = 0;      /* In / out */
 
      s.trianglelist = nullptr;           /* In / out */
      s.triangleattributelist = nullptr;  /* In / out */
      s.trianglearealist = nullptr;       /* In only */
      s.neighborlist = nullptr;           /* Out only */
      s.numberoftriangles = 0;            /* In / out */
      s.numberofcorners = 0;              /* In / out */
      s.numberoftriangleattributes = 0;   /* In / out */
 
      s.segmentlist = nullptr;            /* In / out */
      s.segmentmarkerlist = nullptr;      /* In / out */
      s.numberofsegments = 0;             /* In / out */
 
      s.holelist = nullptr;               /* In / pointer to array copied out */
      s.numberofholes = 0;                /* In / copied out */
 
      s.regionlist = nullptr;             /* In / pointer to array copied out */
      s.numberofregions = 0;              /* In / copied out */
 
      s.edgelist = nullptr;               /* Out only */
      s.edgemarkerlist = nullptr;         /* Not used with Voronoi diagram; out only */
      s.normlist = nullptr;               /* Used only with Voronoi diagram; out only */
      s.numberofedges = 0;                /* Out only */
   };
 
   auto freeStruct = [] (triangulateio & s) {
      if(s.pointlist != nullptr) free(s.pointlist);                         /* In / out */
      if(s.pointattributelist != nullptr) free(s.pointattributelist);       /* In / out */
      if(s.pointmarkerlist != nullptr) free(s.pointmarkerlist);             /* In / out */
 
      if(s.trianglelist != nullptr) free(s.trianglelist);                   /* In / out */
      if(s.triangleattributelist != nullptr) free(s.triangleattributelist); /* In / out */
      if(s.trianglearealist != nullptr) free(s.trianglearealist);           /* In only */
      if(s.neighborlist != nullptr) free(s.neighborlist);                   /* Out only */
 
      if(s.segmentlist != nullptr) free(s.segmentlist);                     /* In / out */
      if(s.segmentmarkerlist != nullptr) free(s.segmentmarkerlist);         /* In / out */
 
      if(s.holelist != nullptr) free(s.holelist);             /* In / pointer to array copied out */
 
      if(s.regionlist != nullptr) free(s.regionlist);         /* In / pointer to array copied out */
 
      if(s.edgelist != nullptr) free(s.edgelist);             /* Out only */
      if(s.edgemarkerlist != nullptr) free(s.edgemarkerlist); /* Not used with Voronoi diagram; out only */
      if(s.normlist != nullptr) free(s.normlist);             /* Used only with Voronoi diagram; out only */
   };
 
   struct triangulateio in, out;
   initStruct(in); initStruct(out);
 
   /* Define input points. */
 
   in.numberofpoints = fNpoints;
   in.pointlist = (REAL *) malloc(in.numberofpoints * 2 * sizeof(REAL));
 
   for (Int_t i = 0; i < fNpoints; ++i) {
      in.pointlist[2 * i] = fXN[i];
      in.pointlist[2 * i + 1] = fYN[i];
   }
 
   triangulate((char *) "zQN", &in, &out, nullptr);
 
   fTriangles.resize(out.numberoftriangles);
   for(int t = 0; t < out.numberoftriangles; ++t){
      Triangle tri;
 
      auto transform = [&] (const unsigned int v) {
         //each triangle as numberofcorners vertices ( = 3)
         tri.idx[v] = out.trianglelist[t*out.numberofcorners + v];
 
         //printf("triangle %u vertex %u: point %u/%i\n", t, v, tri.idx[v], out.numberofpoints);
 
         //pointlist is [x0 y0 x1 y1 ...]
         tri.x[v] = in.pointlist[tri.idx[v] * 2 + 0];
 
         //printf("\t x: %f\n", tri.x[v]);
 
         tri.y[v] = in.pointlist[tri.idx[v] * 2 + 1];
 
         //printf("\t y: %f\n", tri.y[v]);
      };
 
      transform(0);
      transform(1);
      transform(2);
 
      //see comment in header for CGAL fallback section
      tri.invDenom = 1 / ( (tri.y[1] - tri.y[2])*(tri.x[0] - tri.x[2]) + (tri.x[2] - tri.x[1])*(tri.y[0] - tri.y[2]) );
 
      fTriangles[t] = tri;
 
      auto bx = std::minmax({tri.x[0], tri.x[1], tri.x[2]});
      auto by = std::minmax({tri.y[0], tri.y[1], tri.y[2]});
 
      unsigned int cellXmin = CellX(bx.first);
      unsigned int cellXmax = CellX(bx.second);
 
      unsigned int cellYmin = CellY(by.first);
      unsigned int cellYmax = CellY(by.second);
 
      for(unsigned int i = cellXmin; i <= cellXmax; ++i) {
         for(unsigned int j = cellYmin; j <= cellYmax; ++j) {
            //printf("(%u,%u) = %u\n", i, j, Cell(i,j));
            fCells[Cell(i,j)].insert(t);
         }
      }
   }
 
   freeStruct(in); freeStruct(out);
}
 
/// Triangle implementation for interpolation
/// Finds the Delaunay triangle that the point (xi,yi) sits in (if any) and
/// calculate a z-value for it by linearly interpolating the z-values that
/// make up that triangle.
/// Relay that all the triangles have been found before
/// see comment in class description (in Delaunay2D.h) for implementation details:
/// finding barycentric coordinates and computing the interpolation
double Delaunay2D::DoInterpolateNormalized(double xx, double yy)
{
 
   // compute barycentric coordinates of a point P(xx,yy,zz)
   auto bayCoords = [&](const unsigned int t) -> std::tuple<double, double, double> {
      double la = ((fTriangles[t].y[1] - fTriangles[t].y[2]) * (xx - fTriangles[t].x[2]) +
                   (fTriangles[t].x[2] - fTriangles[t].x[1]) * (yy - fTriangles[t].y[2])) *
                  fTriangles[t].invDenom;
      double lb = ((fTriangles[t].y[2] - fTriangles[t].y[0]) * (xx - fTriangles[t].x[2]) +
                   (fTriangles[t].x[0] - fTriangles[t].x[2]) * (yy - fTriangles[t].y[2])) *
                  fTriangles[t].invDenom;
 
      return std::make_tuple(la, lb, (1 - la - lb));
   };
 
   // function to test if a point with barycentric coordinates (a,b,c) is inside the triangle
   // If the point is outside one or more of the coordinate are negative.
   // If the point is on a triangle edge, one of the coordinate (the one not part of the edge) is zero.
   // Due to numerical error, it can happen that if the point is at the edge the result is a small negative value.
   // Use then a tolerance (of - eps) to still consider the point within the triangle
   auto inTriangle = [](const std::tuple<double, double, double> &coords) -> bool {
      constexpr double eps = -4 * std::numeric_limits<double>::epsilon();
      return std::get<0>(coords) > eps && std::get<1>(coords) > eps && std::get<2>(coords) > eps;
   };
 
   int cX = CellX(xx);
   int cY = CellY(yy);
 
   if (cX < 0 || cX > fNCells || cY < 0 || cY > fNCells)
      return fZout; // TODO some more fancy interpolation here
 
   for (unsigned int t : fCells[Cell(cX, cY)]) {
 
      auto coords = bayCoords(t);
 
      // std::cout << "result of bayCoords " << std::get<0>(coords) <<
      //     "  " << std::get<1>(coords)  << "   " << std::get<2>(coords) << std::endl;
 
      if (inTriangle(coords)) {
 
         // we found the triangle -> interpolate using the barycentric interpolation
 
         return std::get<0>(coords) * fZ[fTriangles[t].idx[0]] + std::get<1>(coords) * fZ[fTriangles[t].idx[1]] +
                std::get<2>(coords) * fZ[fTriangles[t].idx[2]];
      }
   }
 
   // printf("Could not find a triangle for point (%f,%f)\n", xx, yy);
 
   // no triangle found return standard value
   return fZout;
}
 
#else //HAS_CGAL: case of using GCAL
 
/// CGAL implementation of normalize points
void Delaunay2D::DonormalizePoints() {
   for (Int_t n = 0; n < fNpoints; n++) {
      //Point p(xTransformer(fX[n]), yTransformer(fY[n]));
      Point p(linear_transform(fX[n], fOffsetX, fScaleFactorX),
              linear_transform(fY[n], fOffsetY, fScaleFactorY));
 
      fNormalizedPoints.insert(std::make_pair(p, n));
   }
}
 
/// CGAL implementation for finding triangles
void Delaunay2D::DoFindTriangles() {
   fCGALdelaunay.insert(fNormalizedPoints.begin(), fNormalizedPoints.end());
 
   std::transform(fCGALdelaunay.finite_faces_begin(),
                  fCGALdelaunay.finite_faces_end(), std::back_inserter(fTriangles),
                  [] (const Delaunay::Face face) -> Triangle {
 
                     Triangle tri;
 
                     auto transform = [&] (const unsigned int i) {
                        tri.x[i] = face.vertex(i)->point().x();
                        tri.y[i] = face.vertex(i)->point().y();
                        tri.idx[i] = face.vertex(i)->info();
                     };
 
                     transform(0);
                     transform(1);
                     transform(2);
 
                     return tri;
 
                  });
}
 
/// CGAL implementation for interpolation
double Delaunay2D::DoInterpolateNormalized(double xx, double yy)
{
   // Finds the Delaunay triangle that the point (xi,yi) sits in (if any) and
   // calculate a z-value for it by linearly interpolating the z-values that
   // make up that triangle.
 
   // initialise the Delaunay algorithm if needed
    FindAllTriangles();
 
   //coordinate computation
   Point p(xx, yy);
 
   std::vector<std::pair<Point, Coord_type> > coords;
   auto nn = CGAL::natural_neighbor_coordinates_2(fCGALdelaunay, p,
                                                  std::back_inserter(coords));
 
   //std::cout << std::this_thread::get_id() << ": Found " << coords.size() << " points" << std::endl;
 
   if(!nn.third) //neighbor finding was NOT successful, return standard value
      return fZout;
 
   //printf("found neighbors %u\n", coords.size());
 
   Coord_type res = CGAL::linear_interpolation(coords.begin(), coords.end(),
                                               nn.second, Value_access(fNormalizedPoints, fZ));
 
   //std::cout << std::this_thread::get_id() << ": Result " << res << std::endl;
 
   return res;
}
#endif // HAS_GCAL
 
} // namespace Math
} // namespace ROOT