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

}