Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  Marching square algorithm
       4             :  * Purpose:  Core algorithm implementation for contour line generation.
       5             :  * Author:   Oslandia <infos at oslandia dot com>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2018, Oslandia <infos at oslandia dot com>
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : #ifndef MARCHING_SQUARE_UTILITY_H
      13             : #define MARCHING_SQUARE_UTILITY_H
      14             : 
      15             : #include <limits>
      16             : #include <cmath>
      17             : 
      18             : namespace marching_squares
      19             : {
      20             : 
      21             : // This is used to determine the maximum level value for polygons,
      22             : // the one that spans all the remaining plane
      23             : constexpr double Inf = std::numeric_limits<double>::infinity();
      24             : 
      25             : constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
      26             : 
      27             : #define debug(format, ...) CPLDebug("MarchingSquare", format, ##__VA_ARGS__)
      28             : 
      29             : // Perturb a value if it is too close to a level value
      30     2806867 : inline double fudge(double value, double minLevel, double level)
      31             : {
      32             :     // FIXME
      33             :     // This is too "hard coded". The perturbation to apply really depend on
      34             :     // values between which we have to interpolate, so that the result of
      35             :     // interpolation should give coordinates that are "numerically" stable for
      36             :     // classical algorithms to work (on polygons for instance).
      37             :     //
      38             :     // Ideally we should probably use snap rounding to ensure no contour lines
      39             :     // are within a user-provided minimum distance.
      40             : 
      41     2806867 :     const double absTol = 1e-6;
      42             :     // Do not fudge the level that would correspond to the absolute minimum
      43             :     // level of the raster, so it gets included.
      44             :     // Cf scenario of https://github.com/OSGeo/gdal/issues/10167
      45     2806867 :     if (level == minLevel)
      46       92202 :         return value;
      47     2714667 :     return std::abs(level - value) < absTol ? value + absTol : value;
      48             : }
      49             : 
      50             : }  // namespace marching_squares
      51             : #endif