finish getting rid of boost

* Assume dash patterns are all on-off(-on-off...) even though this is
not necessarily true

 * now unroll the uvdraw loop 1 time so that there's no test of whether
   to draw or not, and half as many 'maybe reset i' checks

Makefile

@@ -1,13 +1,21 @@
 .PHONY: default
-default: dashing
+default: dashing dashing.omp
 
 .PHONY: bench
 bench: dashing
 	perf stat --log-fd=2 ./dashing -b -s .1 data/HWOOD6E1.pat data/sf.seg
 
 dashing: main.cc dashing.cc dashing.hh parse_numbers.hh contours_and_segments.hh
-	g++ -W -Wall -O2 -g -std=c++11 $(filter %.cc, $^) -o $@ -lboost_random
+	+g++ -W -Wall -O2 -g -std=c++20 $(filter %.cc, $^) -o $@ -flto -DNDEBUG
 dashing-noopt: main.cc dashing.cc dashing.hh parse_numbers.hh contours_and_segments.hh
-	g++ -W -Wall -O0 -g -std=c++11 $(filter %.cc, $^) -o $@ -lboost_random
+	+g++ -W -Wall -O0 -g -std=c++20 $(filter %.cc, $^) -o $@
 dashing-clang: main.cc dashing.cc dashing.hh parse_numbers.hh
-	clang++ -W -Wall -O2 -g -std=c++11 $(filter %.cc, $^) -o $@ -lboost_random
+	+clang++ -W -Wall -O2 -g -std=c++20 $(filter %.cc, $^) -o $@ -flto -DNDEBUG
+dashing.pgo.0: main.cc dashing.cc dashing.hh parse_numbers.hh contours_and_segments.hh
+	+g++ -W -Wall -O2 -g -std=c++20 $(filter %.cc, $^) -o $@ -flto -DNDEBUG -fprofile-generate
+dashing.pgo.1: dashing.pgo.0 main.cc dashing.cc dashing.hh parse_numbers.hh contours_and_segments.hh
+	./dashing.pgo.0 -b -s .002 data/HWOOD6E1.pat  data/sf.seg
+	+g++ -W -Wall -O2 -g -std=c++20 $(filter %.cc, $^) -o $@ -flto -DNDEBUG -fprofile-use
+
+dashing.omp: main.cc dashing.cc dashing.hh parse_numbers.hh contours_and_segments.hh
+	+g++ -W -Wall -O2 -g -std=c++20 $(filter %.cc, $^) -o $@ -flto -DNDEBUG -fopenmp -DDASHING_OMP

README.md

@@ -2,11 +2,11 @@
 
 License: permissive (zlib); see source files for additional details.
 
-On a i5-3320M in benchmark mode, it runs at over 175 million dashes per second:
+On a i5-1235U in multithreaded benchmark mode, it runs at over 3.3 billion dashes per second:
 ```
-$ time ./dashing -b -s .01 data/HWOOD6E1.pat  data/sf.seg
-111430525
-user    0m0.616s
+$ time ./dashing.omp -b -s .002 data/HWOOD6E1.pat  data/sf.seg
+2822685873
+real	0m0.851s
 ```
 
 ![Example image](data/sf.png)
@@ -31,7 +31,7 @@ I would be interested in passing this project to a new maintainer.
     The winding rule `wr` defines which regions are in the interior of the contours.
     For each dash or dot in the resulting hatch, `cb` is called with the output segment.
 
-`xyhatch(const HatchPattern&, const C &segments, Cb cb, const char \*wr)`:
+`xyhatch(const HatchPattern&, const C &segments, Cb cb, const char *wr)`:
     The container C holds segments which must define a set of closed
     contours.
     The winding rule `wr` defines which regions are in the interior of the contours.
@@ -39,7 +39,7 @@ I would be interested in passing this project to a new maintainer.
 
 `HatchPattern::FromFile`: Read a hatch pattern from a file.
 
-`parse\_numbers(std::string line)`: Read a comma and/or space-separated
+`parse_numbers(std::string line)`: Read a comma and/or space-separated
     sequence of numbers into a vector
 
 `SegmentsFromFile`, `ContoursFromFile`, `ContourToSegments`, `ContoursToSegments`: Read and convert segments and contours
@@ -51,6 +51,22 @@ Useful winding rules include:
  
 but any predicate of a single integer may be used.
 
+### Parallel API
+These APIs are available if built with `-fopenmp -DDASHING_OMP`
+
+`xyhatch_omp(const HatchPattern&, It start, It end, Cb cb, Wr wr)`:
+    Iterators `start`..`end` define a range of segments, which must define a set of closed contours.
+    The winding rule `wr` defines which regions are in the interior of the contours.
+    For each dash or dot in the resulting hatch, `cb` is called with the output segment and the thread ID.
+
+`xyhatch_omp(const HatchPattern&, const C &segments, Cb cb, Wr wr)`:
+    The container C holds segments which must define a set of closed
+    contours.
+    The winding rule `wr` defines which regions are in the interior of the contours.
+    For each dash or dot in the resulting hatch, `cb` is called with the output segment and the thread ID.
+
+### Other APIs
+
 Other items in the header files are implementation details.
 
 ## Demo program

contours_and_segments.hh

@@ -12,7 +12,7 @@ typedef std::vector<Contour> Contours;
 
 void ContourToSegments(Segments &dest,
         /* EXPLICIT COPY */ Contour src,
-        double jitter = 0)
+        F jitter = 0)
 {
     if(jitter)
     {
@@ -32,7 +32,7 @@ void ContourToSegments(Segments &dest,
     dest.emplace_back(Segment{src[i-1], src[0], false});
 }
 
-void ContoursToSegments(Segments &dest, Contours const &src, double jitter=0)
+void ContoursToSegments(Segments &dest, Contours const &src, F jitter=0)
 {
     dest.clear();
 
@@ -62,14 +62,14 @@ Contours ContoursFromFile(std::istream &fi)
     return result;
 }
 
-Segments SegmentsFromFile(std::istream &fi, double jitter) {
+Segments SegmentsFromFile(std::istream &fi, F jitter) {
     auto && contours = ContoursFromFile(fi);
     auto result = Segments{};
     ContoursToSegments(result, contours, jitter);
     return result;
 }
 
-std::vector<Segment> SegmentsFromFile(const char *filename, double jitter) {
+std::vector<Segment> SegmentsFromFile(const char *filename, F jitter) {
     std::fstream fi(filename);
     return SegmentsFromFile(fi, jitter);
 }

dashing.cc

@@ -2,37 +2,41 @@
 #include "dashing.hh"
 namespace dashing
 {
-std::vector<double> parse_numbers(std::string line) {
-    boost::algorithm::replace_all(line, ",", " ");
+std::vector<F> parse_numbers(std::string line) {
     std::istringstream fi(line);
-    std::vector<double> result;
-    double d;
-    while((fi >> d)) result.push_back(d);
+    std::vector<F> result;
+    F d;
+    while((fi >> d)) {
+        result.push_back(d);
+        for (auto p = fi.peek(); p == ',' || isspace(p); p = fi.peek()) {
+            fi.get();
+        }
+    }
     return result;
 }
 
 PSMatrix PSMatrix::inverse() const {
-    auto i = 1. / determinant();
+    auto i = F(1.) / determinant();
     return PSMatrix{d*i, -b*i, -c*i, a*i, i*(c*f-e*d), i*(b*e-a*f)};
 }
 
-PSMatrix Translation(double x, double y) {
+PSMatrix Translation(F x, F y) {
     PSMatrix r{1.,0.,0.,1.,x,y};
     return r;
 }
 
-PSMatrix Rotation(double theta) {
-    double c = cos(theta), s = sin(theta);
+PSMatrix Rotation(F theta) {
+    F c = cos(theta), s = sin(theta);
     PSMatrix r{c,s,-s,c,0.,0.};
     return r;
 }
 
-PSMatrix XSkew(double xk) {
+PSMatrix XSkew(F xk) {
     PSMatrix r{1.,0.,xk,1.,0.,0.};
     return r;
 }
 
-PSMatrix YScale(double ys) {
+PSMatrix YScale(F ys) {
     PSMatrix r{1.,0.,0.,ys,0.,0.};
     return r;
 }
@@ -48,21 +52,31 @@ PSMatrix operator*(const PSMatrix &m1, const PSMatrix m2) {
     return r;
 }
 
-static double radians(double degrees) { return degrees * acos(0) / 90.; }
+static F radians(F degrees) { return degrees * acos(0) / 90.; }
 
-Dash::Dash(double th, double x0, double y0, double dx, double dy,
-        const std::vector<double>::const_iterator dbegin,
-        const std::vector<double>::const_iterator dend) : dash(dbegin, dend) {
+Dash::Dash(F th, F x0, F y0, F dx, F dy,
+        const std::vector<F>::const_iterator dbegin,
+        const std::vector<F>::const_iterator dend)
+    : dash{dbegin, dend} {
     auto s = 0.;
-    for(auto d : dash) { sum.push_back(s); s += fabs(d); }
+    for(size_t i = 0; i < dash.size(); i++) {
+        bool is_negative = std::signbit(dash[i]);
+        bool index_is_odd = (i % 2) != 0;
+        if(is_negative != index_is_odd) {
+            throw std::invalid_argument("not a supported dash specification (but probably valid)");
+        }
+        dash[i] = std::abs(dash[i]);
+    }
+    if (dash.size() % 2) dash.push_back(0);
+    for(auto d : dash) { sum.push_back(s); s += d; }
     sum.push_back(s);
 
     tr = Translation(x0, y0) * Rotation(th) * XSkew(dx / dy) * YScale(dy);
     tf = tr.inverse();
 }
 
-Dash Dash::FromString(const std::string &line, double scale) {
-    std::vector<double> words = parse_numbers(line);
+Dash Dash::FromString(const std::string &line, F scale) {
+    std::vector<F> words = parse_numbers(line);
     if(words.size() < 5)
         throw std::invalid_argument("not a valid dash specification");
     for(auto i = words.begin() + 1; i != words.end(); i++) *i *= scale;

dashing.hh

@@ -18,43 +18,47 @@ freely, subject to the following restrictions:
 3. This notice may not be removed or altered from any source distribution.
 */
 
-#ifndef DASHING_H
-#define DASHING_H
+#pragma once
 
+#include <algorithm>
 #include <cmath>
 #include <cassert>
 #include <vector>
 #include <string>
 #include <fstream>
 #include <stdexcept>
-#include <boost/algorithm/string/replace.hpp>
-#include <boost/algorithm/string/trim.hpp>
+#include <span>
+#if defined(DASHING_OMP)
+#include <omp.h>
+#endif
+
+#include "dashing_F.hh"
 
 namespace dashing
 {
 
 struct PSMatrix {
-    double a, b, c, d, e, f;
+    F a, b, c, d, e, f;
 
     PSMatrix inverse() const;
-    double determinant() const { return a * d - b * c; }
+    F determinant() const { return a * d - b * c; }
 };
 
-PSMatrix Translation(double x, double y);
-PSMatrix Rotation(double theta);
-PSMatrix XSkew(double xk);
-PSMatrix YScale(double ys);
+PSMatrix Translation(F x, F y);
+PSMatrix Rotation(F theta);
+PSMatrix XSkew(F xk);
+PSMatrix YScale(F ys);
 
-struct Point { double x, y; };
+struct Point { F x, y; };
 inline Point operator*(const Point &p, const PSMatrix &m) {
     return Point{ p.x*m.a + p.y*m.c + m.e,
                   p.x*m.b + p.y*m.d + m.f };
 }
-inline Point operator*(const Point &p, double d)
+inline Point operator*(const Point &p, F d)
 {
     return Point{ p.x * d, p.y * d };
 }
-inline Point operator*(double d, const Point &p)
+inline Point operator*(F d, const Point &p)
 {
     return Point{ p.x * d, p.y * d };
 }
@@ -67,17 +71,17 @@ PSMatrix operator*(const PSMatrix &m1, const PSMatrix m2);
 
 struct Dash {
     PSMatrix tr, tf;
-    std::vector<double> dash, sum;
+    std::vector<F> dash, sum;
 
-    Dash(double th, double x0, double y0, double dx, double dy,
-            const std::vector<double>::const_iterator dbegin,
-            const std::vector<double>::const_iterator dend);
+    Dash(F th, F x0, F y0, F dx, F dy,
+            const std::vector<F>::const_iterator dbegin,
+            const std::vector<F>::const_iterator dend);
 
-    static Dash FromString(const std::string &line, double scale);
+    static Dash FromString(const std::string &line, F scale);
 };
 
 struct Segment { Point p, q; bool swapped; };
-struct Intersection { double u; bool positive; };
+struct Intersection { F u; bool positive; };
 inline bool operator<(const Intersection &a, const Intersection &b)
 {
     return a.u < b.u;
@@ -90,16 +94,16 @@ inline void ysort(Segment &s) {
     std::swap(s.p, s.q);
 }
 
-inline double intceil(double x) { return int(ceil(x)); }
-inline double intfloor(double x) { return int(floor(x)); }
+inline int intceil(F x) { return int(ceil(x)); }
+inline int intfloor(F x) { return int(floor(x)); }
 
-inline double pythonmod(double a, double b) {
+inline F pythonmod(F a, F b) {
     auto r = a - floor(a / b) * b;
     if(r == b) return 0;
     return r;
 }
 
-inline size_t utoidx(const Dash &d, double u, double &o) {
+inline size_t utoidx(const Dash &d, F u, F &o) {
     u = pythonmod(u, d.sum.back());
     for(size_t i = 1; i != d.sum.size(); i++) {
         if(u < d.sum[i]) { o = u - d.sum[i-1]; return i-1; }
@@ -108,26 +112,30 @@ inline size_t utoidx(const Dash &d, double u, double &o) {
 }
 
 template<class Cb>
-void uvdraw(const Dash &pattern, double v, double u1, double u2, Cb cb) {
+void uvdraw(const Dash &pattern, F v, F u1, F u2, Cb cb) {
     if(pattern.dash.empty()) { cb(v, u1, u2); return;  }
-    double o;
+    F o;
     auto i = utoidx(pattern, u1, o);
-    const auto &pi = pattern.dash[i];
-    if(pi >= 0) { cb(v, u1, std::min(u2, u1+pi-o)); u1 += pi-o; }
+    const auto pi = pattern.dash[i];
+    if(i % 2 == 0) { cb(v, u1, std::min(u2, u1+pi-o)); u1 += pi-o; }
     else { u1 -= pi+o; }
-    i = i + 1;
-    if(i == pattern.dash.size()) i = 0;
+    i++;
+    if(i % 2) {
+        u1 += pattern.dash[i];
+        i++;
+    }
     for(auto u = u1; u < u2;) {
-        const auto &pi = pattern.dash[i];
-        if(pi >= 0) { cb(v, u, std::min(u2, u+pi)); u += pi; }
-        else { u -= pi; }
-        i = i + 1;
-        if(i == pattern.dash.size()) i = 0;
+        if(i >= pattern.dash.size()) i = 0;
+        const auto pi = pattern.dash[i];
+        cb(v, u, std::min(u2, u+pi));
+        u += pi;
+        u += pattern.dash[i+1];
+        i+=2;
     }
 }
 
 template<class Cb, class Wr>
-void uvspans(const Dash &pattern, std::vector<Segment> && segments, Cb cb, std::vector<Intersection> &uu, Wr wr) {
+void uvspans(const Dash &pattern, std::vector<Segment> & segments, Cb cb, std::vector<Intersection> &uu, Wr wr) {
     if(segments.empty()) return; // no segments
 
     for(auto &s : segments) ysort(s);
@@ -173,16 +181,16 @@ void uvspans(const Dash &pattern, std::vector<Segment> && segments, Cb cb, std::
             segments_begin ++;
         }
 
-        for(const auto &s : boost::make_iterator_range(heap_begin, heap_end)) {
+        for(const auto &s : std::span(heap_begin, heap_end)) {
             auto du = s.q.x - s.p.x;
             auto dv = s.q.y - s.p.y;
             assert(dv);
-            if(dv) uu.push_back(
+            uu.push_back(
                     Intersection{s.p.x + du * (v - s.p.y) / dv,s.swapped});
         }
         std::sort(uu.begin(), uu.end());
         int winding = 0;
-        double old_u = -std::numeric_limits<double>::infinity();
+        F old_u = -std::numeric_limits<F>::infinity();
         for(const auto &isect : uu) {
             if(wr(winding)) uvdraw(pattern, v, old_u, isect.u, cb);
             winding += 2*isect.positive - 1;
@@ -193,14 +201,16 @@ void uvspans(const Dash &pattern, std::vector<Segment> && segments, Cb cb, std::
 
 struct HatchPattern {
     std::vector<Dash> d;
-    static HatchPattern FromFile(std::istream &fi, double scale) {
+    static HatchPattern FromFile(std::istream &fi, F scale) {
         HatchPattern result;
 
         std::string line;
         while(getline(fi, line)) {
             auto i = line.find(";");
             if(i != line.npos) line.erase(i, line.npos);
-            boost::algorithm::trim(line);
+            while(!line.empty() && isspace(line.back())) {
+                line.pop_back();
+            }
             if(line.empty()) continue;
             if(line[0] == '*') continue;
             result.d.push_back(Dash::FromString(line, scale));
@@ -208,34 +218,34 @@ struct HatchPattern {
         return result;
     }
 
-    static HatchPattern FromFile(const char *filename, double scale) {
+    static HatchPattern FromFile(const char *filename, F scale) {
         std::ifstream fi(filename);
         return FromFile(fi, scale);
     }
 };
 
 template<class It, class Cb, class Wr>
-void xyhatch(const Dash &pattern, It start, It end, Cb cb, std::vector<Segment> &uvsegments, std::vector<Intersection> &uu, Wr wr) {
-    uvsegments.clear();
+void xyhatch(const Dash &pattern, It start, It end, Cb cb, Wr wr) {
+    std::vector<Segment> uvsegments;
+    uvsegments.reserve(end-start);
+
+    std::vector<Intersection> uu;
+    uu.reserve(8);
+
     bool swapped = pattern.tf.determinant() < 0;
     std::transform(start, end, std::back_inserter(uvsegments),
         [&](const Segment &s)
         { return Segment{s.p * pattern.tf, s.q * pattern.tf, swapped != s.swapped };
     });
-    uvspans(pattern, std::move(uvsegments), [&](double v, double u1, double u2) {
+    uvspans(pattern, uvsegments, [&](F v, F u1, F u2) {
         Point p{u1, v}, q{u2, v};
-        Segment xy{ p * pattern.tr, q * pattern.tr, false };
-        cb(xy);
+        cb(Segment{ p * pattern.tr, q * pattern.tr, false });
     }, uu, wr);
 }
 
 template<class It, class Cb, class Wr>
 void xyhatch(const HatchPattern &pattern, It start, It end, Cb cb, Wr wr) {
-    std::vector<Segment> uvsegments;
-    uvsegments.reserve(end-start);
-    std::vector<Intersection> uu;
-    uu.reserve(8);
-    for(const auto &i : pattern.d) xyhatch(i, start, end, cb, uvsegments, uu, wr);
+    for(const auto &i : pattern.d) xyhatch(i, start, end, cb, wr);
 }
 
 template<class C, class Cb, class Wr>
@@ -243,5 +253,20 @@ void xyhatch(const HatchPattern &pattern, const C &c, Cb cb, Wr wr) {
     xyhatch(pattern, c.begin(), c.end(), cb, wr);
 }
 
+#if defined(DASHING_OMP)
+template<class It, class Cb, class Wr>
+void xyhatch_omp(const HatchPattern &pattern, It start, It end, Cb cb, Wr wr) {
+    #pragma omp parallel for
+    for(const auto &i : pattern.d) {
+        int iam = omp_get_thread_num();
+        xyhatch(i, start, end, [iam, &cb](Segment s) { cb(s, iam); }, wr);
+    }
 }
+template<class C, class Cb, class Wr>
+void xyhatch_omp(const HatchPattern &pattern, const C &c, Cb cb, Wr wr) {
+    xyhatch_omp(pattern, c.begin(), c.end(), cb, wr);
+}
+
 #endif
+
+}

dashing_F.hh

@@ -0,0 +1,26 @@
+/*
+Copyright (c) 2015 Jeff Epler
+
+This software is provided 'as-is', without any express or implied
+warranty. In no event will the authors be held liable for any damages
+arising from the use of this software.
+
+Permission is granted to anyone to use this software for any purpose,
+including commercial applications, and to alter it and redistribute it
+freely, subject to the following restrictions:
+
+1. The origin of this software must not be misrepresented; you must not
+   claim that you wrote the original software. If you use this software
+   in a product, an acknowledgement in the product documentation would be
+   appreciated but is not required.
+2. Altered source versions must be plainly marked as such, and must not be
+   misrepresented as being the original software.
+3. This notice may not be removed or altered from any source distribution.
+*/
+
+#pragma once
+namespace dashing
+{
+
+typedef float F;
+}

main.cc

@@ -20,6 +20,12 @@ freely, subject to the following restrictions:
 
 // demo and benchmark program for dashing
 
+#include <new>
+#include <atomic>
+#include <cstring>
+#if defined(DASHING_OMP)
+#include <omp.h>
+#endif
 #include "dashing.hh"
 #include "contours_and_segments.hh"
 
@@ -85,9 +91,20 @@ int main(int argc, char **argv) {
     if(xit) return 0;
 
     if(bench) {
-        int nseg = 0;
-        auto print_seg = [&nseg](const Segment &s) { (void)s; nseg ++; };
-        xyhatch(h, s, print_seg, [](int i) { return i != 0; } );
+#if defined(DASHING_OMP)
+        struct per_thread_count_type {
+            alignas(std::hardware_destructive_interference_size) size_t n;
+        };
+        auto max_threads = omp_get_max_threads();
+        std::vector<per_thread_count_type> thread_nseg(max_threads, per_thread_count_type{0});
+        auto count_seg = [&thread_nseg](const Segment &s, int iam) { (void)s; thread_nseg[iam].n++; };
+        xyhatch_omp(h, s, count_seg, [](int i) { return i != 0; } );
+        size_t nseg = std::accumulate(thread_nseg.begin(), thread_nseg.end(), size_t{}, [](size_t a, const per_thread_count_type &b) { return a + b.n; });
+#else
+        size_t nseg{};
+        auto count_seg = [&nseg](const Segment &s) { (void)s; nseg ++; };
+        xyhatch(h, s, count_seg, [](int i) { return i != 0; } );
+#endif
         std::cout << nseg << "\n";
         return 0;
     }

parse_numbers.hh

@@ -3,9 +3,9 @@
 #include <vector>
 #include <string>
 #include <sstream>
-#include <boost/algorithm/string/replace.hpp>
+#include "dashing_F.hh"
 namespace dashing
 {
-std::vector<double> parse_numbers(std::string line);
+std::vector<F> parse_numbers(std::string line);
 }
 #endif