Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "PointReduction.h"
      11             : #include "MooseUtils.h"
      12             : #include "MooseError.h"
      13             : 
      14             : #include <algorithm>
      15             : #include <cmath>
      16             : 
      17             : namespace PointReduction
      18             : {
      19             : 
      20             : Real
      21           0 : sqr(Real a)
      22             : {
      23           0 :   mooseDeprecated("use PointReduction::square() instead");
      24           0 :   return a * a;
      25             : }
      26             : 
      27             : Real
      28     1676220 : square(Real a)
      29             : {
      30     1676220 :   return a * a;
      31             : }
      32             : 
      33             : Real
      34      838110 : perpendicularDistance(const FunctionNode & node,
      35             :                       const FunctionNode & begin,
      36             :                       const FunctionNode & end)
      37             : {
      38      838110 :   const Real x0 = node.first;
      39      838110 :   const Real y0 = node.second;
      40      838110 :   const Real x1 = begin.first;
      41      838110 :   const Real y1 = begin.second;
      42      838110 :   const Real x2 = end.first;
      43      838110 :   const Real y2 = end.second;
      44             : 
      45      838110 :   const Real denom = std::sqrt(square(y2 - y1) + square(x2 - x1));
      46             :   mooseAssert(MooseUtils::absoluteFuzzyGreaterThan(denom, 0.0),
      47             :               "Line begin and end points must not be the same");
      48             : 
      49      838110 :   return std::abs((y2 - y1) * x0 - (x2 - x1) * y0 + x2 * y1 - y2 * x1) / denom;
      50             : }
      51             : 
      52             : void
      53         247 : douglasPeuckerRecurse(const FunctionNodeList & list,
      54             :                       const Real epsilon,
      55             :                       std::vector<bool> & keep,
      56             :                       std::size_t begin,
      57             :                       std::size_t end)
      58             : {
      59             :   // Find the point with the maximum distance from the line defined by begin and end
      60         247 :   Real dmax = 0.0;
      61         247 :   std::size_t index = 0;
      62             : 
      63      838357 :   for (std::size_t i = begin; i <= end; ++i)
      64      838110 :     if (keep[i])
      65             :     {
      66      838110 :       const Real d = perpendicularDistance(list[i], list[begin], list[end]);
      67      838110 :       if (d > dmax)
      68             :       {
      69      141531 :         index = i;
      70      141531 :         dmax = d;
      71             :       }
      72             :     }
      73             : 
      74             :   // If max distance is greater than epsilon, recursively simplify
      75         247 :   if (dmax > epsilon)
      76             :   {
      77             :     // Recursive call
      78         117 :     douglasPeuckerRecurse(list, epsilon, keep, begin, index);
      79         117 :     douglasPeuckerRecurse(list, epsilon, keep, index, end);
      80             :   }
      81             :   else
      82             :   {
      83             :     // remove all points between begin and end
      84      129987 :     for (std::size_t i = begin + 1; i < end; ++i)
      85      129857 :       keep[i] = false;
      86             :   }
      87         247 : }
      88             : 
      89             : FunctionNodeList
      90          13 : douglasPeucker(const FunctionNodeList & list, const Real epsilon)
      91             : {
      92             :   // set up keep list for function nodes
      93          13 :   std::vector<bool> keep(list.size(), true);
      94          13 :   douglasPeuckerRecurse(list, epsilon, keep, 0, list.size() - 1);
      95             : 
      96          13 :   FunctionNodeList result;
      97      130013 :   result.reserve(std::count_if(keep.begin(), keep.end(), [](bool k) { return k; }));
      98             : 
      99             :   /// filter result
     100      130013 :   for (std::size_t i = 0; i < list.size(); ++i)
     101      130000 :     if (keep[i])
     102         143 :       result.push_back(list[i]);
     103             : 
     104          26 :   return result;
     105          13 : }
     106             : 
     107             : } // namespace PointReduction