Ticket #310: 1121710e18c9.patch

lemon/bits/bezier.h

@@ -19 +19 @@
 #ifndef LEMON_BEZIER_H
 #define LEMON_BEZIER_H
 
+#include <cmath>
+
 //\ingroup misc
 //\file
 //\brief Classes to compute with Bezier curves.
@@ -66 +68 @@
   Point norm() const { return rot90(p2-p1); }
   Point grad(double) const { return grad(); }
   Point norm(double t) const { return rot90(grad(t)); }
+  Box<double> boundingBox() const {
+    double minx, miny;
+    double maxx, maxy;
+    if (p1.x < p2.x) {
+      minx = p1.x; maxx = p2.x;
+    } else {
+      minx = p2.x; maxx = p1.x;
+    }
+    if (p1.y < p2.y) {
+      miny = p1.y; maxy = p2.y;
+    } else {
+      miny = p2.y; maxy = p1.y;
+    }
+    return Box<double>(minx, miny, maxx, maxy);
+  }
 };
 
 class Bezier2 : public BezierBase
@@ -100 +117 @@
   Bezier1 norm() const { return Bezier1(2.0*rot90(p2-p1),2.0*rot90(p3-p2)); }
   Point grad(double t) const { return grad()(t); }
   Point norm(double t) const { return rot90(grad(t)); }
+  Box<double> boundingBox() const {
+    double minx, miny;
+    double maxx, maxy;
+    if (p1.x < p3.x) {
+      minx = p1.x; maxx = p3.x;
+    } else {
+      minx = p3.x; maxx = p1.x;
+    }
+    if (p1.y < p3.y) {
+      miny = p1.y; maxy = p3.y;
+    } else {
+      miny = p3.y; maxy = p1.y;
+    }
+    if (p1.x + p3.x - 2 * p2.x != 0.0) {
+      double t = (p1.x - p2.x) / (p1.x + p3.x - 2 * p2.x);
+      if (t >= 0.0 && t <= 1.0) {
+        double x = ((1-t)*(1-t))*p1.x+(2*(1-t)*t)*p2.x+(t*t)*p3.x;
+        if (x < minx) minx = x;
+        if (x > maxx) maxx = x;
+      }
+    }
+    if (p1.y + p3.y - 2 * p2.y != 0.0) {
+      double t = (p1.y - p2.y) / (p1.y + p3.y - 2 * p2.y);
+      if (t >= 0.0 && t <= 1.0) {
+        double y = ((1-t)*(1-t))*p1.y+(2*(1-t)*t)*p2.y+(t*t)*p3.y;
+        if (y < miny) miny = y;
+        if (y > maxy) maxy = y;
+      }
+    }
+    return Box<double>(minx, miny, maxx, maxy);
+  }
 };
 
 class Bezier3 : public BezierBase
@@ -150 +198 @@
                                   3.0*rot90(p4-p3)); }
   Point grad(double t) const { return grad()(t); }
   Point norm(double t) const { return rot90(grad(t)); }
+  Box<double> boundingBox() const {
+    double minx, miny;
+    double maxx, maxy;
+    if (p1.x < p4.x) {
+      minx = p1.x; maxx = p4.x;
+    } else {
+      minx = p4.x; maxx = p1.x;
+    }
+    if (p1.y < p4.y) {
+      miny = p1.y; maxy = p4.y;
+    } else {
+      miny = p4.y; maxy = p1.y;
+    }
+    if (p1.x - 3 * p2.x + 3 * p3.x - p4.x != 0.0) {
+      double disc = p2.x * p2.x + p3.x * p3.x - p1.x * p3.x 
+        + p1.x * p4.x - p2.x * p3.x - p2.x * p4.x;
+      if (disc >= 0.0) {
+        double t = (p1.x + p3.x - 2 * p2.x + std::sqrt(disc)) / 
+          (p1.x - 3 * p2.x + 3 * p3.x - p4.x);
+        if (t >= 0.0 && t <= 1.0) {
+          double x = ((1-t)*(1-t)*(1-t))*p1.x+(3*t*(1-t)*(1-t))*p2.x+
+            (3*t*t*(1-t))*p3.x+(t*t*t)*p4.x;
+          if (x < minx) minx = x;
+          if (x > maxx) maxx = x;
+        }
+        t = (p1.x + p3.x - 2 * p2.x - std::sqrt(disc)) / 
+          (p1.x - 3 * p2.x + 3 * p3.x - p4.x);
+        if (t >= 0.0 && t <= 1.0) {
+          double x = ((1-t)*(1-t)*(1-t))*p1.x+(3*t*(1-t)*(1-t))*p2.x+
+            (3*t*t*(1-t))*p3.x+(t*t*t)*p4.x;
+          if (x < minx) minx = x;
+          if (x > maxx) maxx = x;
+        }
+      }
+    } else if (p1.x + p3.x - 2 * p2.x != 0.0) {
+      double t = (p1.x - p2.x) / (2 * p1.x + 2 * p3.x - 4 * p2.x);
+      if (t >= 0.0 && t <= 1.0) {
+        double x = ((1-t)*(1-t)*(1-t))*p1.x+(3*t*(1-t)*(1-t))*p2.x+
+          (3*t*t*(1-t))*p3.x+(t*t*t)*p4.x;
+        if (x < minx) minx = x;
+        if (x > maxx) maxx = x;
+      }      
+    }
+    if (p1.y - 3 * p2.y + 3 * p3.y - p4.y != 0.0) {
+      double disc = p2.y * p2.y + p3.y * p3.y - p1.y * p3.y 
+        + p1.y * p4.y - p2.y * p3.y - p2.y * p4.y;
+      if (disc >= 0.0) {
+        double t = (p1.y + p3.y - 2 * p2.y + std::sqrt(disc)) / 
+          (p1.y - 3 * p2.y + 3 * p3.y - p4.y);
+        if (t >= 0.0 && t <= 1.0) {
+          double y = ((1-t)*(1-t)*(1-t))*p1.y+(3*t*(1-t)*(1-t))*p2.y+
+            (3*t*t*(1-t))*p3.y+(t*t*t)*p4.y;
+          if (y < miny) miny = y;
+          if (y > maxy) maxy = y;
+        }
+        t = (p1.y + p3.y - 2 * p2.y - std::sqrt(disc)) / 
+          (p1.y - 3 * p2.y + 3 * p3.y - p4.y);
+        if (t >= 0.0 && t <= 1.0) {
+          double y = ((1-t)*(1-t)*(1-t))*p1.y+(3*t*(1-t)*(1-t))*p2.y+
+            (3*t*t*(1-t))*p3.y+(t*t*t)*p4.y;
+          if (y < miny) miny = y;
+          if (y > maxy) maxy = y;
+        }
+      }
+    } else if (p1.y + p3.y - 2 * p2.y != 0.0) {
+      double t = (p1.y - p2.y) / (2 * p1.y + 2 * p3.y - 4 * p2.y);
+      if (t >= 0.0 && t <= 1.0) {
+        double y = ((1-t)*(1-t)*(1-t))*p1.y+(3*t*(1-t)*(1-t))*p2.y+
+          (3*t*t*(1-t))*p3.y+(t*t*t)*p4.y;
+        if (y < miny) miny = y;
+        if (y > maxy) maxy = y;
+      }      
+    }
+    return Box<double>(minx, miny, maxx, maxy);
+  }
 
   template<class R,class F,class S,class D>
   R recSplit(F &_f,const S &_s,D _d) const