Навеяло вопросом в алгоритмах:
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))
{
   . . .
}

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

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

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

Хорошее исследование. Меня оно многому научило. Я попробовал написать свой класс для нерекурсивного преобразования кривой Безье в ломаную. Вот что у меня вышло. Кратко опишу используемые классы

TVector2D — класс вектора, в нем реализованы только те методы, которые необходимы для вычислений

TBezier2D — класс хранит управляющие точки кривой Безье

TBezierCurve2D — класс позволяет последовательно вычислять точки кривой
основной метод

bool TBezierCurve2D::vertex(TVector2D* point)

возвращает точку кривой.

getDistance — функция вычисления расстояния от точки до прямой

Я разделил код на несколько функций. Скорее всего это замедляет работу, но позволяет анализировать алгоритм. В основном код позаимствован у McSeem2. Надеюсь он на меня за это не обидится.

#include "math.h"
#include <stack>

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

//класс вектора
class TVector2D
{
public:
double x,y;
//операция вычитания векторов
friend TVector2D inline operator - (const TVector2D& Vector1,const TVector2D& Vector2);
//операция сложения векторов
friend TVector2D inline operator + (const TVector2D& Vector1,const TVector2D& Vector2);
//операция умножения на скаляр
friend TVector2D inline operator * (const TVector2D& Vector,double Scalar);
//векторное произведение векторов возвращается модуль со знаком
double VectorCrossProduct(const TVector2D& Vector) const
{return x*Vector.y-y*Vector.x;}
TVector2D(){}
TVector2D(double x,double y):x(x),y(y){}
};

//операция вычитания векторов
TVector2D inline operator - (const TVector2D& Vector1,const TVector2D& Vector2)
{
return TVector2D(Vector1.x-Vector2.x,Vector1.y-Vector2.y);
}

//операция сложения векторов
TVector2D inline operator + (const TVector2D& Vector1,const TVector2D& Vector2)
{
return TVector2D(Vector1.x+Vector2.x,Vector1.y+Vector2.y);
}

//операция умножения на скаляр
TVector2D inline operator * (const TVector2D& Vector,double Scalar)
{
return TVector2D(Scalar*Vector.x,Scalar*Vector.y);
}

//вспомогательная функция для определения расстояния
//от точки P до отрезка P0,P1
double inline getDistance(const TVector2D& P,
                          const TVector2D& P0,
                          const TVector2D& P1)
{
//расстояние между прямой и точкой определяется по формуле
//|(P-P0)x(P1-P0)|/|P1-P0|
TVector2D P10=P1-P0;

//знак нам не важен поэтому вычисляем векторное произведение наоборот
double cross=P10.VectorCrossProduct(P-P0);
return fabs(cross)/sqrt(P10.x*P10.x+P10.y*P10.y);
}

//класс предназначен для хранения 4 точек кривой Безье
//кривая может быть представлена в виде
//C(u)=B0(u)*P0+B1(u)*P1+B2(u)*P2+B3(u)*P3
//где B0(u),B1(u),B2(u),B3(u) - полиномы Берштейна
//так же кривая может быть представлена в виде
//C(u)=a+b*u+c*u^2+d*u^3 хотя данную форму мы не используем
class TBezier2D
{
public:
TBezier2D(){}
TVector2D P0,P1,P2,P3;
};

//класс является чем-то вроде итератора по точкам кривой Безье
class TBezierCurve2D
{
public:
TBezierCurve2D(const TBezier2D& bezier,
               double distance_tolerance=0.25,
               double angle_tolerance=0.3):
               distance_tolerance_(distance_tolerance),
               angle_tolerance_(angle_tolerance)
{
stack_.push(bezier);
firstPoint_=bezier.P0;lastPoint_=bezier.P3;
state_=START;
}

//заносит точку на кривой в point возвращает false если точек нет
bool vertex(TVector2D* point);

private:
//разбиение кривой на левую и правую подкривые в точке 0.5
void subdivide(const TBezier2D& bezier,TBezier2D* left,TBezier2D* right) const;

//проверка насколько точно ломаная интерполирует кривую bezier
//bezier - исходная кривая left,right - ее подкривые
//возвращает true если ломаная приближает кривую в пределах заданной точности
bool checkTolerance(const TBezier2D& bezier,
                    const TBezier2D& left,
                    const TBezier2D& right) const;

//стек для проведения рекурсивного разбиения
std::stack<TBezier2D> stack_;
//точность для расстояния и угла
double distance_tolerance_,angle_tolerance_;
//первая и последняя точка кривой
TVector2D firstPoint_,lastPoint_;
enum TState {START,MOVE,STOP};
//для определения начальной точки и завершения работы
TState state_;
};

//разбиение кривой на левую и правую подкривые в точке 0.5
void TBezierCurve2D::subdivide(const TBezier2D& bezier,
                                     TBezier2D* left,
                                     TBezier2D* right) const
{
//исходная кривая задается точками P0,P1,P2,P3
//кривая left точками P0,P10,P20,P30
//кривая right точками P30,P21,P12,P3
//сначала делим каждый отрезок control polyline исходной кривой пополам
TVector2D P10,P11,P12;
P10=(bezier.P0+bezier.P1)*0.5;
P11=(bezier.P1+bezier.P2)*0.5;
P12=(bezier.P2+bezier.P3)*0.5;

//отрезки соединяющие точки вычисленные на предыдущем шаге опять делим пополам
TVector2D P20,P21;
P20=(P10+P11)*0.5;
P21=(P11+P12)*0.5;

//наконец считаем точку P30 в этой точке будут соединяться левая и правая
//части исходной кривой, так же эта точка соответствует точке C(0.5) исходной кривой
TVector2D P30;
P30=(P20+P21)*0.5;

//разнесем полученные точки по левой и правой части кривой
//левая часть
left->P0=bezier.P0;
left->P1=P10;
left->P2=P20;
left->P3=P30;
//правая часть
right->P0=P30;
right->P1=P21;
right->P2=P12;
right->P3=bezier.P3;
}

//проверка насколько точно ломаная интерполирует кривую bezier
//bezier - исходная кривая left,right - ее подкривые
//возвращает true если ломаная приближает кривую в пределах заданной точности
bool TBezierCurve2D::checkTolerance(const TBezier2D& bezier,
                                    const TBezier2D& left,
                                    const TBezier2D& right) const
{
        TVector2D P11=(bezier.P1+bezier.P2)*0.5;

        //содрано у McSeem2
        // Calculate the primary curvature as the sum of distances:
        // d123 + d234 + d1234
       double curvature =
            getDistance(bezier.P1,left.P1,P11)+getDistance(P11,left.P2,right.P1)+
            getDistance(bezier.P2,P11,left.P2);

        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.
                //----------------------
                return true;
            }
            TVector2D P01=bezier.P1-bezier.P0;
            TVector2D P12=bezier.P2-bezier.P1;
            TVector2D P23=bezier.P3-bezier.P2;

            double a12 = atan2(P01.y, P01.x);
            double a23 = atan2(P12.y, P12.x);
            double a34 = atan2(P23.y, P23.x);

            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
                //----------------------
                return true;
            }
        }
        return false;
}

//заносит точку на кривой в point возвращает false если точек нет
bool TBezierCurve2D::vertex(TVector2D* point)
{
if(state_==STOP) return false;
//возвращаем первую точку вначале
if(state_==START)
{
*point=firstPoint_;
state_=MOVE;
return true;
}
//последняя точка
if(stack_.empty())
{
*point=lastPoint_;
state_=STOP;
return true;
}

TBezier2D left,right;
while(1)
{
subdivide(stack_.top(),&left,&right);
if(checkTolerance(stack_.top(),left,right))
{
*point=left.P3;
stack_.pop();
return true;
}
stack_.pop();
stack_.push(right);
stack_.push(left);
}
}

Пример использования

//управляющие точки кривой
TBezier2D bezier;

//заносим точки
bezier.P0 = TVector2D(10, 10);
bezier.P1 = TVector2D(70, 70);
bezier.P2 = TVector2D(107, 108);
bezier.P3 = TVector2D(207, 211);

//вычисление точек
TBezierCurve2D iter(bezier);
TVector2D point;
while( iter.vertex(&point) )
{
    //какие-то действия над точкой
    foo(point.x, point.y);
}

Здравствуйте, runtime2, Вы писали:

R>Я разделил код на несколько функций. Скорее всего это замедляет работу, но позволяет анализировать алгоритм. В основном код позаимствован у McSeem2. Надеюсь он на меня за это не обидится.

Замечательно! Нисколько не обижусь.
Но зря ты убрал проверку из getDistance. Две точки отрезка могут совпадать и тогда — деление на ноль. В данном случае это важно, поскольку ситуация может запросто возникнуть в вырожденных случаях. Вообще, числовая стабильность в подобных задачах — это основная проблема.

Здравствуйте, McSeem2, Вы писали:

MS>Замечательно! Нисколько не обижусь.
MS>Но зря ты убрал проверку из getDistance. Две точки отрезка могут совпадать и тогда — деление на ноль. В данном случае это важно, поскольку ситуация может запросто возникнуть в вырожденных случаях. Вообще, числовая стабильность в подобных задачах — это основная проблема.

Огромное спасибо за отзыв! Спасибо за уточнение функции getDistance. Я действительно не учел деление на 0 (у меня нет такого опыта в написание подобного кода, как у тебя). Я ее писал по аналогии с функцией agg::calc_point_line_distance, а надо было брать функцию с исходника на форуме.

Я компилировал пример bezier_div.cpp на C++ Builder 6.0 SP4 с полной оптимизацией. Программа работает примерно в 2 раза медленее, чем откомпилированная в архиве. Я несколько разочарован в своем компиляторе.

Посмотрел твой пример. Зашел на www.antigrain.com. Скачал дистрибутив библиотеки. Откомпилировал примеры.
Нет слов — одни эмоции. Круто!

KK>Посмотрел твой пример. Зашел на www.antigrain.com. Скачал дистрибутив библиотеки. Откомпилировал примеры.
KK>Нет слов — одни эмоции. Круто!

Спасибо на добром слове! Но как кто-то заметил, нам еще расти и расти.

Здравствуйте, runtime2, Вы писали:

R>Огромное спасибо за отзыв! Спасибо за уточнение функции getDistance. Я действительно не учел деление на 0 (у меня нет такого опыта в написание подобного кода, как у тебя). Я ее писал по аналогии с функцией agg::calc_point_line_distance, а надо было брать функцию с исходника на форуме.

R>Я компилировал пример bezier_div.cpp на C++ Builder 6.0 SP4 с полной оптимизацией. Программа работает примерно в 2 раза медленее, чем откомпилированная в архиве. Я несколько разочарован в своем компиляторе.

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

r=1000, num_points=517
r=10K,  num_points=2053
r=100K, num_points=8197
r=1M,   num_points=16389
r=10M,  num_points=65541

При этом точность аппроксимации остается той же!

Я там еще добавил кое-что, в частности, cusp_limit. И глубину рекурсии надо все-таки ограничивать. Полагаться на одинаковость поведения чисел с плавающей точкой на разных платформах нельзя. Скоро сделаю более подробный отчет.

Сколько человек пишут библиотеку ?

Здравствуйте, Ka3a4oK, Вы писали:

KK>Сколько человек пишут библиотеку ?

Я один и пишу. Люди иногда помогают идеями, алгоритмами, математикой.

Здравствуйте, McSeem2, Вы писали:

MS>Здравствуйте, Ka3a4oK, Вы писали:

KK>>Сколько человек пишут библиотеку ?

MS>Я один и пишу. Люди иногда помогают идеями, алгоритмами, математикой.

Почему один ? Проект вроде опенсоурс ? Неужели желающих нет ?

Здравствуйте, Ka3a4oK, Вы писали:

KK>Почему один ? Проект вроде опенсоурс ? Неужели желающих нет ?

Как ни странно, но желающих что-то добавить в основную часть библиотеки не много. Например, Tony Juricic сделал NURBS и ими вполне можно пользоваться. Но я пока что не включил их в библиотеку, поскольку не знаю как интегрировать в общий конвейер — там кроме координат нужны еще и веса.

Другая часть, непосредственно не относящаяся к библиотеке — это agg_platform_support.cpp. Это обертка над неким API для компиляции демо-примеров. Я написал варианты для WinAPI и X11. Потом Hansruedi Baer написал вариант для CarbonAPI, Stephan Assmus — для BeOS, Steven Solie — для AmigaOS. Кто бы вот еще помог с QNX Proton API?

Nickolas Zimmermann (AKA wildfox) добавил файлы для automake, таким образом AGG попал в дистрибутив SuSE, чем я очень горжусь.

А вообще, архитектура позволяет заменять и/или переписывать части библиотеки, не трогая исходников из основного пакета. Например, народ из Liberty Technology Systems дописал свои специализированные трансформеры изображений. Они так же спонсировали доработки — 48- и 64-битовые цвета и экранное пространство коодинат в 32 бита (раньше было 16, я даже не предполагал, что кому-то когда-то может понадобиться буфер шириной больше 32K пикселов).
Carl Sassenrath использует AGG в своем REBOL, но исторически там alpha интерпретируется наоборот (255=прозрачно, 0=opaque). Поэтому они используют свои версии некоторых файлов.
Это я к тому, что включать подобные частные решения в библтотеку может и не иметь смысла, но завсегда можно использовать их как некие собственные расширения.

Ну а вообще, есть некоторые действительно нетривиальные вещи, требующие решения, вот, например:
http://www.rsdn.ru/Forum/Message.aspx?mid=1093589&only=1

Автор: McSeem2
Дата: 27.03.05

Но люди (и я в том числе) как правило заняты более насущными проблемами. Увы, но такова жизнь.

У меня рекурсивный алгоритм обгоняет нерекурсивный тоже где-то на 200%. Я провел следующие улучшения алгоритма.

1. Часть замедления связана с использованием функции TBezierCurve2D::vertex. Для того чтобы уравнять рекурсивный и нерекурсивный алгоритм я заменил эту функцию на TBezierCurve2D::getPolyline, которая сразу записывает все точки в массив. Так же я немного оптимизировал функцию TBezierCurve2D::subdivide. В результате отставание снизилось процентов на 20%.

2. Так же, я переписал алгоритм работы со стеком. Стек представляет собой глобальный массив на 1000 элементов. Ты уже говорил, в примере с окружностями, что количество точек растет медленно по сравнению с размером кривой. Вообще, если я правильно понял твою идею, наибольшее количество точек получается если кривая близка к дуге окружности. Глубина стека тоже растет очень медленно. Для 1000 точек глубина стека равна 10. Вообще глубина стека скорее всего зависит от расстояния между начальной и конечной точкой кривой и будет максимальной для управляющей ломаной с каждым отрезком равным этому расстоянию. Для маленьких кривых можно пользоваться статическим стеком. Для больших выделять динамически стек в зависимости от расстояния. Хотя я думаю, что 1000 для простых задач достаточно.

Переписанный алгоритм работы со стеком уменьшает разницу на 30%. Это не обязательно из-за медленности std:stack, в данном алгоритме просто используется знание о структуре данных в стеке.

Остальное скорее всего вызвано использованием TVector2D, если переписать нерекурсивный алгоритм так же как в твоем рекурсивном, то я думаю что скрость не будет сильно отличаться. Я не проверял просто мысль.

3. Я выделил функцию checkTolerance в отдельный шаблонный параметр. Так же туда я поместил distance_tolerance_ и angle_tolerance_. Я думаю, что так будет проще исследовать класс алгоритмов, которые я опишу ниже.

4. Что-бы избавиться от операций вычисления квадратных корней и операций деления я изменил алгоритм вычисления предварительной кривизны.

        // 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));

Я решил для оценки кривизны использовать расстояние точек P1 и P2 до прямой P0P3. Можно брать максимальное из этих расстояний. Но вроде бы кривая при этом выглядит не очень хорошо. Поэтому я, так же как и ты, взял сумму этих расстояний. Это позволяет избавиться от квадратных корней и операций деления. Так же на надо вычислять лишний раз left и right части кривой для оценки кривизны.

double cross1, cross2;
        //расстояние между прямой P30 и точкой P определяется по формуле
        //|(P-P0)x(P30)|/|P30|
        //вктор соединяющий точки P3 и P0 кривой
        TVector2D P30 = bezier.P3 - bezier.P0;
        //квадрат длинны вектора P30
        double quadricLength = P30.x * P30.x + P30.y * P30.y;
        //вектор соединяющий точки P1 и P0
        TVector2D P01 = bezier.P1 - bezier.P0;
        //вектор соединяющий точки P3 и P2
        TVector2D P23 = bezier.P3 - bezier.P2;

        //модуль векторного произведения вкторов P30 и P01
        //знак нам не важен поэтому вычисляем векторное произведение наоборот
        cross1=fabs(P30.VectorCrossProduct(P01));
        //модуль вкторного произведения векторов P30 и P23
        cross2=fabs(P30.VectorCrossProduct(P23));

        double curvature=cross1+cross2;

        if(curvature*curvature < distance_tolerance_*quadricLength)
        {
        //далее определение углов
        }

Очень интересное наблюдение! Твой алгритм и моя модификация дают одинаковое количество точек при distance_length=0.5 для твоего и distance_length=1 для моей модификации. Скорее всего это вызывано тем, что твои расстояния можно пересчитать через мои расстояиния и отрезок соединяющий начальную и конечную точку кривой. Причем эта зависимость получается приблизительно линейная с постоянным коэффициентом (или мало меняющимся).

В исходнике для вычислений я использую distance_length=1, что дает неплохие результаты. Как я уже говорил, у меня мало опыта в серьезном исследовании таких вещей. Поэтому моя оценка неплохие результаты может быть сильно тобой скорректирована.

У меня модифицированный нерекурсивный алгоритм рабтает примерно в 2 раза быстрее, чем исходный нерекурсивный.

На данный момент вычисления связанные с арктангенсами составляют 25% всего времени выполнения функции checkTolerance. Не знаю можно ли что-нибудь улучшить с арктангенсами. Для двух углов вычисляется три арктангенса.

Если точку P0 двигать возле точки P1, то возникает ошибка при вычислении арктангенсов. Поэтому мне пришлось ввести проверку.

if(P01.y < intersection_epsilon ||
   P12.y < intersection_epsilon ||
   P23.y < intersection_epsilon)
  {
       return;
  }

У тебя в исходнике bezier_div.cpp так же возникает такая ошибка. Для исправления я добавлял

if(fabs(y2 - y1) < intersection_epsilon ||
   fabs(y3 - y2) < intersection_epsilon ||
   fabs(y4 - y3) < intersection_epsilon)
  {
       path->line_to(x1234, y1234);
       return;
  }

Я использую класс TPolyline2D в качестве массива векторов.

//стек для разбиения кривой
static TBezier2D stack[1000];

//класс предназначен для вычисления точек на кривой Безье
//TAlgorithm - алгоритм проверки окончания разбиения должен содержать
//функцию checkTolerance
template<class TAlgorithm>
class TBezierCurve2D:public TAlgorithm
{
public:
//заносит точки кривой bezier в массив polyline
void getSplinePolyline(const TBezier2D& bezier, TPolyline2D* polyline);
private:
//разбиение кривой на левую и правую подкривые в точке 0.5
void subdivide(const TBezier2D& bezier,TBezier2D* left,TBezier2D* right) const;
};

//разбиение кривой на левую и правую подкривые в точке 0.5
template<class TAlgorithm>
void TBezierCurve2D<TAlgorithm>::subdivide(const TBezier2D& bezier,
                                           TBezier2D* left,
                                           TBezier2D* right) const
{
//исходная кривая задается точками P0,P1,P2,P3
//средняя точка отрезка P1P2
TVector2D P11;
P11=(bezier.P1+bezier.P2)*0.5;
//левая часть
left->P0=bezier.P0;
left->P1=(left->P0+bezier.P1)*0.5;;
left->P2=(left->P1+P11)*0.5;
//правая часть
right->P3=bezier.P3;
right->P2=(right->P3+bezier.P2)*0.5;
right->P1=(P11+right->P2)*0.5;
right->P0=(left->P2+right->P1)*0.5;
left->P3=right->P0;
}

//заносит точки кривой bezier в массив polyline
template<class TAlgorithm>
void TBezierCurve2D<TAlgorithm>::getSplinePolyline(const TBezier2D& bezier, TPolyline2D* polyline)
{
//указатель на вершину стека это первый свободный элемент
TBezier2D* top=stack;
//начальный адрес стека
TBezier2D* base=stack;
//первая точка
polyline->Add(bezier.P0);
*top=bezier;
while(1)
{
//последняя точка
if(top<base)
{
polyline->Add(bezier.P3);
return;
}
//вызваем функцию унаследованую из класса TAlgorithm
if(checkTolerance(*top))
{
polyline->Add(top->P3);
top--;
continue;
}
//top bezier
//top+1 right
//top+2 left
subdivide(*top,top+2,top+1);
memcpy(top,top+1,2*sizeof(TBezier2D));
top++;
}
}

//алгоритм для проверки точности
class TQuadricDistanceAngle
{
public:
TQuadricDistanceAngle()
{
//начальные значения точности
init(1,0.3);
}
void init(double distance_tolerance, double angle_tolerance)
{
//используем distance_tolerance в квадрате
distance_tolerance_=distance_tolerance*distance_tolerance;
angle_tolerance_=angle_tolerance;
}
bool checkTolerance(const TBezier2D& bezier);
private:
//точность для расстояния и угла
double distance_tolerance_,angle_tolerance_;
};


bool TQuadricDistanceAngle::checkTolerance(const TBezier2D& bezier)
{

        double cross1, cross2;
        //расстояние между прямой P30 и точкой P определяется по формуле
        //|(P-P0)x(P30)|/|P30|
        //вктор соединяющий точки P3 и P0 кривой
        TVector2D P30 = bezier.P3 - bezier.P0;
        //квадрат длинны вектора P30
        double quadricLength = P30.x * P30.x + P30.y * P30.y;
        //вектор соединяющий точки P1 и P0
        TVector2D P01 = bezier.P1 - bezier.P0;
        //вектор соединяющий точки P3 и P2
        TVector2D P23 = bezier.P3 - bezier.P2;

        //модуль векторного произведения вкторов P30 и P01
        //знак нам не важен поэтому вычисляем векторное произведение наоборот
        cross1=fabs(P30.VectorCrossProduct(P01));
        //модуль вкторного произведения векторов P30 и P23
        cross2=fabs(P30.VectorCrossProduct(P23));

        double curvature=cross1+cross2;

        if(curvature*curvature < distance_tolerance_*quadricLength)
        {
            // If the curvature doesn't exceed the distance_tolerance value
            // we tend to finish subdivisions.
            //----------------------
            TVector2D P12=bezier.P2-bezier.P1;

            if(P01.y < intersection_epsilon ||
               P12.y < intersection_epsilon ||
               P23.y < intersection_epsilon)
            {
                // This condition protects us from degenerate cases, such as
                // (p1 == p4 && p2 == p3) and others.
                //----------------------
                return true;
            }


            double a12 = atan2(P01.y, P01.x);
            double a23 = atan2(P12.y, P12.x);
            double a34 = atan2(P23.y, P23.x);

            double da = fabs(a23 - a12) + fabs(a34 - a23);
            if(da >= 2*pi) da = 4*pi - da;

            if(da < angle_tolerance_)
            {
                // Finally we can stop the recursion
                //----------------------
                return true;
            }



        }
        return false;

}

Пример использования

//для хранения управляющих точек кривой
TBezier2D bezier;
//для хранения точек кривой
TPolyline2D p;
//управляющие точки кривой
bezier.P0 = TVector2D(7, 17);
bezier.P1 = TVector2D(7, 109);
bezier.P2 = TVector2D(107, 109);
bezier.P3 = TVector2D(107, 9);
//вычисление точек кривой
TBezierCurve2D< TQuadricDistanceAngle > curve;
curve.getSplinePolyline(bezier, p);

Здравствуйте, runtime2, Вы писали:

R>Я решил для оценки кривизны использовать расстояние точек P1 и P2 до прямой P0P3. Можно брать максимальное из этих расстояний. Но вроде бы кривая при этом выглядит не очень хорошо. Поэтому я, так же как и ты, взял сумму этих расстояний. Это позволяет избавиться от квадратных корней и операций деления. Так же на надо вычислять лишний раз left и right части кривой для оценки кривизны.

Вот за эту идею — огромное спасибо! Похоже, что эта оценка работает лучше всего, я ее и оставлю. Она работает очень похоже на мой метод с тремя расстояниями, но оценивает ошибку более точно. Тем не менее, твой код не справится со случаем, когда начальнная и конечная точки (P0,P3 в твоих терминах) совпадают (петля). Казалось бы, это решается элементарно — надо принудительно заставить кривую делиться первый раз. Даже есть доказательство, что кривая Безье третьего порядка после одного деления пополам не может иметь петли. Но не все так просто. Ни моя, ни твоя оценки не спасают от случая P2----P0----P3----p2. Кроме того, есть вырожденные случаи с изломами, например, такой: (475, 157, 200, 100, 453, 100, 222, 157). Тем действительно есть бесконечная кривизна в некой точке, то есть, производная имеет разрыв.

Довольно интересный и простой метод предложил Timothee Groleau (http://www.timotheegroleau.com). Оценка
кривизны считается как расстояние между точкой на кривой (в моих терминах — 1234) и серединой отрезка 1-4.

            double x14 = (x1 + x4) / 2;
            double y14 = (y1 + y4) / 2;
            double dx  = x1234 - x14;
            double dy  = y1234 - y14;
            double distance = dx*dx + dy*dy;

При этом кривую надо тоже принудительно поделить один раз, иначе не сработает случай формы "Z" — точка на кривой попадает в точности на середину отрезка 1-4. Но этот способ оценивает расстояние между точками, в то время как реальная ошибка аппроксимации оценивается расстоянием между точкой и прямой.

Таким образом, я пришел к следующему коду:

    //------------------------------------------------------------------------
    const double curve_distance_epsilon         = 1e-20;
    const double curve_collinearity_epsilon     = 1e-30;
    const double curve_angle_tolerance_epsilon  = 0.01;
    enum       { curve_recursion_limit          = 32 };

    //------------------------------------------------------------------------
    bool curve4_div::angle_condition(double x1,    double y1, 
                                     double x2,    double y2, 
                                     double x3,    double y3, 
                                     double x4,    double y4,
                                     double x1234, double y1234)
    {
        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 < m_angle_tolerance)
        {
            // Finally we can stop the recursion
            //----------------------
            m_points.add(point_type(x1234, y1234));
            return true;
        }

        if(m_cusp_limit != 0.0)
        {
            if(da1 > m_cusp_limit)
            {
                m_points.add(point_type(x2, y2));
                return true;
            }

            if(da2 > m_cusp_limit)
            {
                m_points.add(point_type(x3, y3));
                return true;
            }
        }
        return false;
    }

    void curve4_div::recursive_bezier_line_point(double x1, double y1, 
                                                 double x2, double y2, 
                                                 double x3, double y3, 
                                                 double x4, double y4,
                                                 unsigned level)
    {
        if(level > curve_recursion_limit) return;

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

        if(level > 0) // Enforce subdivision first time
        {
            // Try to approximate the full cubic curve by a single straight line
            //------------------
            double dx = x4-x1;
            double dy = y4-y1;
            double d14 = sqrt(dx*dx + dy*dy);
            double distance = 0;

            if(d14 > curve_distance_epsilon)
            {
                double e = fabs(((x2 - x4) * dy - (y2 - y4) * dx)) + 
                           fabs(((x3 - x4) * dy - (y3 - y4) * dx));
                if(e > curve_collinearity_epsilon)
                {
                    distance = e / d14;
                }
                else
                {
                    double x14 = (x1 + x4) / 2;
                    double y14 = (y1 + y4) / 2;
                    dx = x1234 - x14;
                    dy = y1234 - y14;
                    distance = dx*dx + dy*dy;
                }
            }
            else
            {
                distance = calc_distance(x1, y1, x2, y2) + calc_distance(x4, y4, x3, y3);
            }

            // Another method. Works slightly worse and doesn't handle 
            // the case of 2-----1-----4-----3
            //double distance = 
            //    fabs(calc_line_point_distance(x12,  y12,  x23,  y23,  x2,  y2)) + 
            //    fabs(calc_line_point_distance(x23,  y23,  x34,  y34,  x3,  y3)) +
            //    fabs(calc_line_point_distance(x123, y123, x234, y234, x23, y23));

            if(distance <= m_distance_tolerance)
            {
                // If the curvature doesn't exceed the distance_tolerance value
                // we tend to finish subdivisions.
                //----------------------
                if(m_angle_tolerance < curve_angle_tolerance_epsilon || 
                   distance <= curve_distance_epsilon)
                {
                    m_points.add(point_type(x1234, y1234));
                    return;
                }
                if(angle_condition(x1, y1, x2, y2, x3, y3, x4, y4, x1234, y1234)) return;
            }
        }

        // Continue subdivision
        //----------------------
        recursive_bezier_line_point(x1, y1, x12, y12, x123, y123, x1234, y1234, level + 1); 
        recursive_bezier_line_point(x1234, y1234, x234, y234, x34, y34, x4, y4, level + 1); 
    }

    void curve4_div::bezier_line_point(double x1, double y1, 
                                       double x2, double y2, 
                                       double x3, double y3, 
                                       double x4, double y4)
    {
        m_points.add(point_type(x1, y1));
        recursive_bezier_line_point(x1, y1, x2, y2, x3, y3, x4, y4, 0);
        m_points.add(point_type(x4, y4));
    }

Здравствуйте, runtime2, Вы писали:

R>Если точку P0 двигать возле точки P1, то возникает ошибка при вычислении арктангенсов. Поэтому мне пришлось ввести проверку.
R>

R>if(P01.y < intersection_epsilon ||
R>   P12.y < intersection_epsilon ||
R>   P23.y < intersection_epsilon)
R>  {
R>       return;
R>  }
R>

Быть такого не может. Какая именно ошибка возникает?
Функция atan2 устойчива и dy=0 для нее — нормальное явление. Очень малая величина dy — тоже нормально. Даже atan2(0,0) выдает строгий 0.

MS>Здравствуйте, runtime2, Вы писали:

R>>Если точку P0 двигать возле точки P1, то возникает ошибка при вычислении арктангенсов. Поэтому мне пришлось ввести проверку.

MS>Быть такого не может. Какая именно ошибка возникает?
MS>Функция atan2 устойчива и dy=0 для нее — нормальное явление. Очень малая величина dy — тоже нормально. Даже atan2(0,0) выдает строгий 0.

На самом деле, ты прав. Если точки совпадают, то atan2 выдает 0. А угол между векторами вычисляется как разность углов наклона каждого из векторов. Таким образом, в случае одного вырожденного вектора получается полная фигня. Но твоя проверка все равно не спасает, это надо обрабатывать по-другому. Как конкретно — скоро напишу подробный отчет и приведу окончательный вариант. Непростое это дело оказалось...

Здравствуйте, McSeem2, Вы писали:

MS>Вот за эту идею — огромное спасибо! Похоже, что эта оценка работает лучше всего, я ее и оставлю. Она работает очень похоже на мой метод с тремя расстояниями, но оценивает ошибку более точно. Тем не менее, твой код не справится со случаем, когда начальнная и конечная точки (P0,P3 в твоих терминах) совпадают (петля).

Мой код справится с случаем, когда начальная и конечная точки совпадают (P0 и P3).

    
    //вектор соединяющий точки P3 и P0 кривой
    TVector2D P30 = bezier.P3 - bezier.P0;
    //квадрат длинны вектора P30
    double quadricLength = P30.x * P30.x + P30.y * P30.y;

    ...
    
    if(curvature*curvature < distance_tolerance_*quadricLength)
    {
        //разбиение подходит с точки зрения расстояний
    }

Для случая P0=P3 quadricLength=0 и соответственно условие if не выполняется.

MS>Ни моя, ни твоя оценки не спасают от случая P2----P0----P3----p2.

Ты имеешь ввиду случай P1-P0-P3-P2? Да, в этом случае у меня в результирующей ломаной оказывается только первая и последняя точка.

MS>Кроме того, есть вырожденные случаи с изломами, например, такой: (475, 157, 200, 100, 453, 100, 222, 157). Тем действительно есть бесконечная кривизна в некой точке, то есть, производная имеет разрыв.

Да, твой последний алгоритм выдает в центре закругление. Мой алгоритм выдает острый угол в этом месте. За счет каких строк в твоем алгоритме это происходит?

Здравствуйте, McSeem2, Вы писали:

MS>>Здравствуйте, runtime2, Вы писали:

R>>>Если точку P0 двигать возле точки P1, то возникает ошибка при вычислении арктангенсов. Поэтому мне пришлось ввести проверку.

MS>>Быть такого не может. Какая именно ошибка возникает?
MS>>Функция atan2 устойчива и dy=0 для нее — нормальное явление. Очень малая величина dy — тоже нормально. Даже atan2(0,0) выдает строгий 0.

MS>На самом деле, ты прав. Если точки совпадают, то atan2 выдает 0. А угол между векторами вычисляется как разность углов наклона каждого из векторов. Таким образом, в случае одного вырожденного вектора получается полная фигня. Но твоя проверка все равно не спасает, это надо обрабатывать по-другому. Как конкретно — скоро напишу подробный отчет и приведу окончательный вариант. Непростое это дело оказалось...

Моя проверка не подходит, так как считает точным любое разбиение с управляющей ломаной, содержащей отрезок параллельный оси Y.
Я попробовал использовать код

            double da;
            //если одновременно x и y окажется равным 0
            try
            {
            double a12 = atan2(P01.y, P01.x);
            double a23 = atan2(P12.y, P12.x);
            double a34 = atan2(P23.y, P23.x);

            da = fabs(a23 - a12) + fabs(a34 - a23);
            if(da >= 2*pi) da = 4*pi - da;
            }
            catch(...)
            {
            return true;
            }

У меня он работает на 30% медленне чем исходный с if. Как я понимаю исключения сильно замедляют работу (хотя ни одного исключения при измерении времени не возникало). Или я что-то напутал. Да и вообще судя по твоему сообщению нельзя из за ошибки в одном atan2 считать разбиение точным.

У меня есть еще несколько вопросов, если можно.

if(m_angle_tolerance < curve_angle_tolerance_epsilon || 
   distance <= curve_distance_epsilon)
{
    m_points.add(point_type(x1234, y1234));
    return;
}

Зачем в этом коде проверка m_angle_tolerance < curve_angle_tolerance_epsilon?

Что такое m_cusp_limit и чему он должен быть равен?

Еще раньше я пробовал ускорить алгоритм вычисляя сумму углов через углы между P0P1 и P0P3, P2P3 и P0P3. Но потом отказался от этого вот кусок кода

//вычисление угла в нерекурсивной реализации
double testAngle2(const TBezier2D& bezier)
{

            TVector2D P01=bezier.P1-bezier.P0;
            TVector2D P12=bezier.P2-bezier.P1;
            TVector2D P23=bezier.P3-bezier.P2;
            //отрезок соединяющий начальную и конечную точку кривой
            TVector2D a=bezier.P3-bezier.P0;

            double da1;
            double da2;

            da1=acos(a.VectorDotProduct(P01)/(sqrt(P01.x*P01.x+P01.y*P01.y)*
                                                   sqrt(a.x*a.x+a.y*a.y)));

            da2=acos(a.VectorDotProduct(P23)/(sqrt(P23.x*P23.x+P23.y*P23.y)*
                                                   sqrt(a.x*a.x+a.y*a.y)));

            da1 = fabs(da1);
            da2 = fabs(da2);
            double da=da1+da2;
            if(da>2*pi) da=4*pi-da;
            return da;
}

VectorDotProduct — скалярное произведение векторов. Еще я не делал здесь проверку на деление на ноль. Углы вычисленные по этому методу не отличались от углов по твоему варианту.