separated removal of linear redundancy from other refinement =>

satisfactory result for all the wallpaper.

pub/lib/plot.cpp

@@ -288,25 +288,37 @@ int plot_curve_refine(CPlot* plot, const COlc* fc, int k, int j, int cstyle)
         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;
-        vector<CPt> ret(x.size());
-        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}};
-            else if (y[i]>plot->Y.sup)
-                ret[i] = {x[i], 0., 1, {xp, NAN}};
-            else if (!inRange[i])
-                ret[i] = {x[i], 0., 1, {xp, NAN}};
-            else
-                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; };
+    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;
+            vector<CPt> ret(x.size());
+            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}};
+                else if (y[i]>plot->Y.sup)
+                    ret[i] = {x[i], 0., 1, {xp, NAN}};
+                else if (!inRange[i])
+                    ret[i] = {x[i], 0., 1, {xp, NAN}};
+                else
+                    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; };
+
+    auto nonlinearity = [](const vector<CPt>::iterator b) -> double
+        {
+            auto a = b-1;
+            auto c = b+1;
+            CP2 BA = a->P-b->P;
+            CP2 AC = c->P-a->P;
+            // squared  deviation from linearity:
+            return BA*BA - pow(BA*AC,2)/(AC*AC);
+        };
 
     // start with equidistant grid:
     vector<double> xn;
@@ -314,45 +326,47 @@ int plot_curve_refine(CPlot* plot, const COlc* fc, int k, int j, int cstyle)
 
     // refinement loop:
     vector<CPt> pp;
-    for (int iref=0; iref<47 && pp.size() < plot->maxpoints; ++iref) {
+    for (int iref=0; iref<49 && pp.size() < plot->maxpoints; ++iref) {
         vector<CPt> pn = x2p(xn);
         pp = triv::merge_sorted(pp, pn, pt_order);
         xn.clear();
 
         auto a = pp.begin();
         auto b = pp.begin()+1;
-        vector<CPt> po = { *a };
         for (; b<pp.end(); ++a, ++b) {
             if (a->code!=0 || b->code!=0) {
-                if (a->code!=b->code)
+                if (a->code!=b->code) {
                     xn.push_back( (b->xd + a->xd)/2 );
+                }
             } else { // valid data
                 auto c = b+1;
                 if (c<pp.end() && c->code==0) {
-                    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
+                    double delta2 = nonlinearity(b);
                     if ( delta2 > 1e-8 ) {
-                        xn.push_back( (b->xd + a->xd)/2 );
-                        xn.push_back( (c->xd + b->xd)/2 );
+                        double dx = std::min(b->xd - a->xd, c->xd - b->xd)/2;
+                        xn.push_back( b->xd - dx );
+                        xn.push_back( b->xd + dx );
                     }
                 }
             }
-            po.push_back( *b );
         }
-        pp = po;
         if (!xn.size())
             break;
         xn = std::vector<double>( xn.begin(), std::unique(xn.begin(), xn.end()) );
     }
 
+    // remove points in straight lines
+    vector<CPt> po { *pp.begin() };
+    for (auto b = pp.begin()+1; b<pp.end()-1; ++b)
+        if (b->code!=0 || (b-1)->code!=0 | (b+1)->code!=0 || nonlinearity(b)>1e-9)
+            po.push_back(*b);
+    po.push_back( *(pp.end()-1) );
+
     // divide into segments
     vector<vector<CPt>> segments;
     vector<CPt> seg;
-    for (CPt p: pp) {
+    for (CPt p: po) {
         if (p.code==0)
             seg.push_back(p);
         else if (seg.size()) {