Навеяло вопросом в алгоритмах:
http://www.rsdn.ru/Forum/Message.aspx?mid=1217769&only=1

Автор:
Дата: 10.06.05

Продолжение вопроса:
http://www.rsdn.ru/Forum/Message.aspx?mid=1218153&only=1

Автор:
Дата: 11.06.05

На что я лихо ответил, что типа "все элементарно", но в действительности оказалось что все совсем не так, как на самом деле. Тем не менее, кое-что придумалось.

    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;
            }

            // Calculate the secondary curvature value. To do that
            // we estimate the sum of the angles: a1+a2, where 
            // a1 is the absolute value of the angle between vectors:
            //     (p1,p2)->(p2,p3) 
            // a2 is the absolute value of the angle between vectors:
            //     (p2,p3)->(p3,p4)
            // 
            // Calculation of the angles is expensive, so, the estimation is 
            // calculated as:
            // 2.0 - (d13 / (d12+d23) + d24 / (d23+d34)); 
            //----------------------
            double d12 = calc_distance(x1,y1, x2,y2);
            double d23 = calc_distance(x2,y2, x3,y3);
            double d34 = calc_distance(x3,y3, x4,y4);
            curvature = 2.0 - (calc_distance(x1,y1, x3,y3) / (d12 + d23) +
                               calc_distance(x2,y2, x4,y4) / (d23 + d34));

            if(curvature < 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/250);
        path.line_to(x4, y4);
    }

Пояснения.
Алгоритм аппроксимирует кубическую кривую Безье кусочно-линейным способом. Цель задачи — сохранить "гладкость" кривой при минимальном количестве отрезков. За основу взят классический рекурсивный способ разбиения пополам. Чертеж:


Итак, кривая задана четырьмя точками — 1,2,3,4. 1,4 — это концы, 2,3 — контрольные точки ("усы" в графических редакторах). Деление кривой пополам делается очень просто. Вычисляются средние точки — 12,23,34, после чего — средние точки 123 и 234, и наконец — точка 1234, лежащая на кривой. Таким образом мы получили две более "плоские" кривые — 1,12,123,1234 и 1234,234,34,4. Повторяем операцию рекурсивно до тех пор пока... Вот в этом "пока" вся сложность и заключается. В простейшем случае берем, например, расстояние между точками 1 и 4 и завершаем деление, если оно меньше какого-то значения. Это плохой способ, он выдает много точек на плоских участках, но слишком мало на крутых. К тому же, точки 1 и 4 могут совпадать.

Более хороший способ — это оценка мгновенной "курватуры" (не знаю, как по-русски, но "кризизна кривой" как-то не звучит). В результате деления, мы получаем три треугольника — 12,2,23, 123,23,234 и 23,3,34. Если сумма площадей этих треугольников не превышает некого tolerance, можно считато кривую достаточно плоской. Но площадь — это тоже плохой критерий, поскольку треугольники могут получаться очень длинными и узкими. Гораздо лучше оказалась просто сумма расстояний (высот треугольников) — d123+d1234+d234.
В принципе, на этом можно и закончить, критерий работает вполне приемлемо. Значение "distance_tolerance" приблизительно соответствует абсолютному значению максимального отклонения. Например, 0.25 пиксела на экранном разрешении получается вполне точно.

Дальше — улучшения. На очень крутых поворотах, как здесь, точек получается все-же недостаточно, особенно если нам надо считать эквидистанту к кривой (например, stroke). Бледно-зеленая "пипка" — это тот самый "не-фонтан", который получается при инкрементальном способе с равномерным распределением точек и последующим вычислением stroke. При этом инкрементальный способ выдает почти в 8 раз больше точек, в большинстве своем бесполезных.

Для таких случаев можно посчитать сумму углов между исходными векторами (1,2)->(2,3) и (2,3)->(3,4). Но вычисление углов — операция дорогая, поэтому мы заменяем ее на некую оценочную функцию с расстояниями:
d13/(d12+d23) + d24/(d23+d34).
Здесь уже возможны деления на ноль, но простая проверка чуть выше защищает от всех вырожденных случаев:

            if(curvature < intersection_epsilon)
            {
                path.line_to(x1234, y1234);
                return;
            }

Значение angle_tolerance — нормализовано эмпирически, в функции bezier: angle_tolerance/250. Таким образом, значение 1.0, передаваемое в функцию соответствует добавлению нескольких дополнительных точек на изломах. Чем оно меньше, там больше точек.

Еднственный вырожденный случай, плохой для данного алгоритма, это когда все 4 точки лежат на одной прямой в следующем порядке:
2-----1------4------3
При этом, результат выродится в отрезок 1----4. Не знаю пока как его побороть.
Другая актуальная задача — раскрутить рекурсию в нечто event-driven, чтобы можно было пользоваться так:

bezier c(. . .);
double x, y;
while(c.vertex(&x, &y))
{
   . . .
}

Может кто поможет?
Поиграться можно здесь.