lemon/bits/bezier.h
author Peter Kovacs <kpeter@inf.elte.hu>
Fri, 03 Apr 2009 18:59:15 +0200
changeset 600 6ac5d9ae1d3d
parent 314 2cc60866a0c9
child 792 761fe0846f49
permissions -rw-r--r--
Support real types + numerical stability fix in NS (#254)

- Real types are supported by appropriate inicialization.
- A feature of the XTI spanning tree structure is removed to ensure
numerical stability (could cause problems using integer types).
The node potentials are updated always on the lower subtree,
in order to prevent overflow problems.
The former method isn't notably faster during to our tests.
     1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library.
     4  *
     5  * Copyright (C) 2003-2009
     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7  * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8  *
     9  * Permission to use, modify and distribute this software is granted
    10  * provided that this copyright notice appears in all copies. For
    11  * precise terms see the accompanying LICENSE file.
    12  *
    13  * This software is provided "AS IS" with no warranty of any kind,
    14  * express or implied, and with no claim as to its suitability for any
    15  * purpose.
    16  *
    17  */
    18 
    19 #ifndef LEMON_BEZIER_H
    20 #define LEMON_BEZIER_H
    21 
    22 //\ingroup misc
    23 //\file
    24 //\brief Classes to compute with Bezier curves.
    25 //
    26 //Up to now this file is used internally by \ref graph_to_eps.h
    27 
    28 #include<lemon/dim2.h>
    29 
    30 namespace lemon {
    31   namespace dim2 {
    32 
    33 class BezierBase {
    34 public:
    35   typedef lemon::dim2::Point<double> Point;
    36 protected:
    37   static Point conv(Point x,Point y,double t) {return (1-t)*x+t*y;}
    38 };
    39 
    40 class Bezier1 : public BezierBase
    41 {
    42 public:
    43   Point p1,p2;
    44 
    45   Bezier1() {}
    46   Bezier1(Point _p1, Point _p2) :p1(_p1), p2(_p2) {}
    47 
    48   Point operator()(double t) const
    49   {
    50     //    return conv(conv(p1,p2,t),conv(p2,p3,t),t);
    51     return conv(p1,p2,t);
    52   }
    53   Bezier1 before(double t) const
    54   {
    55     return Bezier1(p1,conv(p1,p2,t));
    56   }
    57 
    58   Bezier1 after(double t) const
    59   {
    60     return Bezier1(conv(p1,p2,t),p2);
    61   }
    62 
    63   Bezier1 revert() const { return Bezier1(p2,p1);}
    64   Bezier1 operator()(double a,double b) const { return before(b).after(a/b); }
    65   Point grad() const { return p2-p1; }
    66   Point norm() const { return rot90(p2-p1); }
    67   Point grad(double) const { return grad(); }
    68   Point norm(double t) const { return rot90(grad(t)); }
    69 };
    70 
    71 class Bezier2 : public BezierBase
    72 {
    73 public:
    74   Point p1,p2,p3;
    75 
    76   Bezier2() {}
    77   Bezier2(Point _p1, Point _p2, Point _p3) :p1(_p1), p2(_p2), p3(_p3) {}
    78   Bezier2(const Bezier1 &b) : p1(b.p1), p2(conv(b.p1,b.p2,.5)), p3(b.p2) {}
    79   Point operator()(double t) const
    80   {
    81     //    return conv(conv(p1,p2,t),conv(p2,p3,t),t);
    82     return ((1-t)*(1-t))*p1+(2*(1-t)*t)*p2+(t*t)*p3;
    83   }
    84   Bezier2 before(double t) const
    85   {
    86     Point q(conv(p1,p2,t));
    87     Point r(conv(p2,p3,t));
    88     return Bezier2(p1,q,conv(q,r,t));
    89   }
    90 
    91   Bezier2 after(double t) const
    92   {
    93     Point q(conv(p1,p2,t));
    94     Point r(conv(p2,p3,t));
    95     return Bezier2(conv(q,r,t),r,p3);
    96   }
    97   Bezier2 revert() const { return Bezier2(p3,p2,p1);}
    98   Bezier2 operator()(double a,double b) const { return before(b).after(a/b); }
    99   Bezier1 grad() const { return Bezier1(2.0*(p2-p1),2.0*(p3-p2)); }
   100   Bezier1 norm() const { return Bezier1(2.0*rot90(p2-p1),2.0*rot90(p3-p2)); }
   101   Point grad(double t) const { return grad()(t); }
   102   Point norm(double t) const { return rot90(grad(t)); }
   103 };
   104 
   105 class Bezier3 : public BezierBase
   106 {
   107 public:
   108   Point p1,p2,p3,p4;
   109 
   110   Bezier3() {}
   111   Bezier3(Point _p1, Point _p2, Point _p3, Point _p4)
   112     : p1(_p1), p2(_p2), p3(_p3), p4(_p4) {}
   113   Bezier3(const Bezier1 &b) : p1(b.p1), p2(conv(b.p1,b.p2,1.0/3.0)),
   114                               p3(conv(b.p1,b.p2,2.0/3.0)), p4(b.p2) {}
   115   Bezier3(const Bezier2 &b) : p1(b.p1), p2(conv(b.p1,b.p2,2.0/3.0)),
   116                               p3(conv(b.p2,b.p3,1.0/3.0)), p4(b.p3) {}
   117 
   118   Point operator()(double t) const
   119     {
   120       //    return Bezier2(conv(p1,p2,t),conv(p2,p3,t),conv(p3,p4,t))(t);
   121       return ((1-t)*(1-t)*(1-t))*p1+(3*t*(1-t)*(1-t))*p2+
   122         (3*t*t*(1-t))*p3+(t*t*t)*p4;
   123     }
   124   Bezier3 before(double t) const
   125     {
   126       Point p(conv(p1,p2,t));
   127       Point q(conv(p2,p3,t));
   128       Point r(conv(p3,p4,t));
   129       Point a(conv(p,q,t));
   130       Point b(conv(q,r,t));
   131       Point c(conv(a,b,t));
   132       return Bezier3(p1,p,a,c);
   133     }
   134 
   135   Bezier3 after(double t) const
   136     {
   137       Point p(conv(p1,p2,t));
   138       Point q(conv(p2,p3,t));
   139       Point r(conv(p3,p4,t));
   140       Point a(conv(p,q,t));
   141       Point b(conv(q,r,t));
   142       Point c(conv(a,b,t));
   143       return Bezier3(c,b,r,p4);
   144     }
   145   Bezier3 revert() const { return Bezier3(p4,p3,p2,p1);}
   146   Bezier3 operator()(double a,double b) const { return before(b).after(a/b); }
   147   Bezier2 grad() const { return Bezier2(3.0*(p2-p1),3.0*(p3-p2),3.0*(p4-p3)); }
   148   Bezier2 norm() const { return Bezier2(3.0*rot90(p2-p1),
   149                                   3.0*rot90(p3-p2),
   150                                   3.0*rot90(p4-p3)); }
   151   Point grad(double t) const { return grad()(t); }
   152   Point norm(double t) const { return rot90(grad(t)); }
   153 
   154   template<class R,class F,class S,class D>
   155   R recSplit(F &_f,const S &_s,D _d) const
   156   {
   157     const Point a=(p1+p2)/2;
   158     const Point b=(p2+p3)/2;
   159     const Point c=(p3+p4)/2;
   160     const Point d=(a+b)/2;
   161     const Point e=(b+c)/2;
   162     const Point f=(d+e)/2;
   163     R f1=_f(Bezier3(p1,a,d,e),_d);
   164     R f2=_f(Bezier3(e,d,c,p4),_d);
   165     return _s(f1,f2);
   166   }
   167 
   168 };
   169 
   170 
   171 } //END OF NAMESPACE dim2
   172 } //END OF NAMESPACE lemon
   173 
   174 #endif // LEMON_BEZIER_H