more stable computation of deviation from linearity.

pub/lib/plot.cpp

@@ -7,8 +7,9 @@
 //! \file  plot.cpp
 //! \brief NPlot: plot data and curves.
 
-#include "defs.hpp"
+#include <algorithm>
 
+#include "defs.hpp"
 
 #include "../plot/dualplot.hpp"
 #include "../trivia/vector_ops.hpp"
@@ -264,20 +265,29 @@ int plot_curve_equidist(CPlot* plot, const COlc* fc, int k, int j, int cstyle)
     return xp.size();
 }
 
-class CPt {
-public:
-    double xd; // data point
-    double yd;
-    int code; // 0=valid, 1=below, 2=above, 3=off_range
-    double xp; // plot point
-    double yp;
-};
-
 //! Plots slice j of non-convolved curve file fc, refining the grid where necessary;
 //! return number of points plotted.
 
 int plot_curve_refine(CPlot* plot, const COlc* fc, int k, int j, int cstyle)
 {
+    class CP2 {
+    public:
+        double x;
+        double y;
+        CP2 operator-(const CP2 other) {
+            return { x-other.x, y-other.y }; };
+        double operator*(const CP2 other) {
+            return x*other.x + y*other.y; };
+    };
+
+    class CPt {
+    public:
+        double xd; // data point
+        double yd;
+        int code; // 0=valid, 1=below, 2=above, 3=off_range
+        CP2 P; // plot point
+    };
+
     auto x2p = [&](const vector<double>& x) -> vector<CPt> {
         vector<double> y = PCAST<const CObjVecDbl>(fc->eval_curve(x, k, j, false))->v;
         vector<int> inRange = fc->eval_prange(x, k, j)->v;
@@ -285,16 +295,17 @@ int plot_curve_refine(CPlot* plot, const COlc* fc, int k, int j, int cstyle)
         for (size_t i=0; i<x.size(); ++i) {
             double xp = plot->X.pc(x[i]);
             if (y[i]<plot->Y.inf)
-                ret[i] = {x[i], 0., 1, xp, NAN};
+                ret[i] = {x[i], 0., 1, {xp, NAN}};
             else if (y[i]>plot->Y.sup)
-                ret[i] = {x[i], 0., 1, xp, NAN};
+                ret[i] = {x[i], 0., 1, {xp, NAN}};
             else if (!inRange[i])
-                ret[i] = {x[i], 0., 1, xp, NAN};
+                ret[i] = {x[i], 0., 1, {xp, NAN}};
             else
-                ret[i] = {x[i], y[i], 0, xp, plot->Y.pc(y[i])};
+                ret[i] = {x[i], y[i], 0, {xp, plot->Y.pc(y[i])}};
         }
         return ret;
     };
+
     auto pt_order = [](const CPt a, const CPt b) -> bool { return a.xd<=b.xd; };
 
     // start with equidistant grid:
@@ -318,19 +329,24 @@ int plot_curve_refine(CPlot* plot, const COlc* fc, int k, int j, int cstyle)
             } else { // valid data
                 auto c = b+1;
                 if (c<pp.end() && c->code==0) {
-                    if ( std::abs( ((b->xp-a->xp)*c->yp
-                                    +(c->xp-b->xp)*a->yp)
-                                   /(c->xp-a->xp) - b->yp ) < 1e-6 ) // it's linear
+                    CP2 BA = a->P-b->P;
+                    CP2 AC = c->P-a->P;
+                    // squared  deviation from linearity:
+                    double delta2 = BA*BA - pow(BA*AC,2)/(AC*AC);
+                    if ( delta2 < 1e-12 ) // it's linear
                         continue; // don't copy *b to po
+                    if ( delta2 > 1e-8 ) {
+                        xn.push_back( (b->xd + a->xd)/2 );
+                        xn.push_back( (c->xd + b->xd)/2 );
+                    }
                 }
-                if (hypot(b->xp - a->xp, b->yp - a->yp)>1e-1)
-                    xn.push_back( (b->xd + a->xd)/2 );
             }
             po.push_back( *b );
         }
         pp = po;
         if (!xn.size())
             break;
+        xn = std::vector<double>( xn.begin(), std::unique(xn.begin(), xn.end()) );
     }
 
     // divide into segments