Пример использования
//для хранения управляющих точек кривой
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.
Здравствуйте, 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 — скалярное произведение векторов. Еще я не делал здесь проверку на деление на ноль. Углы вычисленные по этому методу не отличались от углов по твоему варианту.
