Все-таки есть один крайне неприятный вырожденный случай (точнее сказать, тип вырожденных случаев), на котором алгоритм ведет себя грубо. Например, такой:

p1=(100,100), p2=(100,200), p3=(100,100), p4=(110,100)

То есть, такой "сапог".
Получается не фонтан, в критическом месте деление прекращается:


Должно быть вот это:


Виновато упрощение для оценки угла. Поэтому, я не мудрствуя лукаво стал вычислять сумму углов между первым-вторым и вторым-третьим отрезками через atan2. По эффективности получается примерно то же самое, как я выяснил. Вычисление миллиона корней занимает 34 миллисекунды, вычисление миллиона арктангенсов — 90 (на Centrino-1.7). Но в исходном варианте вычисляется 5 корней, а в модифицированном — 3 арктангенса. К тому же, это все теряется среди прочих операций.

Теперь angle_tolerance задается непосредственно в радианах. Значение 0.3 дает вполне хороший результат. Чем меньше значение, тем больше появляется точек в местах изломов.

Скачать здесь: http://antigrain.com/stuff/bezier_div.zip

    const double pi = 3.14159265358979323846;
    const double intersection_epsilon = 1.0e-30;

    inline double calc_distance(double x1, double y1, double x2, double y2)
    {
        double dx = x2-x1;
        double dy = y2-y1;
        return sqrt(dx * dx + dy * dy);
    }

    inline double calc_point_line_distance(double x1, double y1, 
                                           double x2, double y2, 
                                           double x,  double y)
    {
        double dx = x2-x1;
        double dy = y2-y1;
        double d = sqrt(dx * dx + dy * dy);
        if(d < intersection_epsilon)
        {
            return calc_distance(x1, y1, x, y);
        }
        return ((x - x2) * dy - (y - y2) * dx) / d;
    }


    template<class PathStorage>
    void recursive_bezier(PathStorage& path,
                          double x1, double y1, 
                          double x2, double y2, 
                          double x3, double y3, 
                          double x4, double y4,
                          double distance_tolerance,
                          double angle_tolerance)
    {
        // Calculate all the mid-points of the line segments
        //----------------------
        double x12   = (x1 + x2) / 2;                
        double y12   = (y1 + y2) / 2;
        double x23   = (x2 + x3) / 2;
        double y23   = (y2 + y3) / 2;
        double x34   = (x3 + x4) / 2;
        double y34   = (y3 + y4) / 2;
        double x123  = (x12 + x23) / 2;
        double y123  = (y12 + y23) / 2;
        double x234  = (x23 + x34) / 2;
        double y234  = (y23 + y34) / 2;
        double x1234 = (x123 + x234) / 2;
        double y1234 = (y123 + y234) / 2;

        // Calculate the primary curvature as the sum of distances: 
        // d123 + d234 + d1234
        //----------------------
        double curvature = 
            fabs(calc_point_line_distance(x12,  y12,  x23,  y23,  x2,  y2)) + 
            fabs(calc_point_line_distance(x23,  y23,  x34,  y34,  x3,  y3)) +
            fabs(calc_point_line_distance(x123, y123, x234, y234, x23, y23));

        if(curvature < distance_tolerance)
        {
            // If the curvature doesn't exceed the distance_tolerance value
            // we tend to finish subdivisions.
            //----------------------
            if(curvature < intersection_epsilon)
            {
                // This condition protects us from degenerate cases, such as
                // (p1 == p4 && p2 == p3) and others.
                //----------------------
                path.line_to(x1234, y1234);
                return;
            }

            double a12 = atan2(y2 - y1, x2 - x1);
            double a23 = atan2(y3 - y2, x3 - x2);
            double a34 = atan2(y4 - y3, x4 - x3);

            double da1 = fabs(a23 - a12);
            double da2 = fabs(a34 - a23);
            if(da1 >= pi) da1 = 2*pi - da1;
            if(da2 >= pi) da2 = 2*pi - da2;

            if(da1 + da2 < angle_tolerance)
            {
                // Finally we can stop the recursion
                //----------------------
                path.line_to(x1234, y1234);
                return;                 
            }
        }

        // Continue subdivision
        //----------------------
        recursive_bezier(path, x1, y1, x12, y12, x123, y123, x1234, y1234, 
                         distance_tolerance, angle_tolerance);

        recursive_bezier(path, x1234, y1234, x234, y234, x34, y34, x4, y4, 
                         distance_tolerance, angle_tolerance);
    }

    template<class PathStorage>
    void bezier(PathStorage& path,
                double x1, double y1, 
                double x2, double y2, 
                double x3, double y3, 
                double x4, double y4,
                double distance_tolerance,
                double angle_tolerance)
    {
        path.move_to(x1, y1);
        recursive_bezier(path, x1, y1, x2, y2, x3, y3, x4, y4, 
                         distance_tolerance, angle_tolerance);
        path.line_to(x4, y4);
    }