Создание скалярной функции на базе кривой Безье: проблемы
Для стандартной области применения Безье (2D и 3D графика) слабо изученных проблем немного. Когда же мы начинаем пытаться «вылепить» на базе 2D-кривой Безье график скалярной функции, вылезает ворох головоломок, о которых написано весьма ограниченное число научных статей. Плюс к ним привязываются проблемы, относящиеся к удобству интерфейса, и о чем научные статьи, увы, не рассказывают.
Проблема 1: Как гарантировать единственность значения для каждого аргумента?
Требование единственности значения – суть определения функции. 2D кривым Безье на это требование естественно наплевать. Они нередко имеют по 2 и более точек для одного значения x
Проблема 2: Как найти значение функции для заданного аргумента?
Фактически, нам нужно отыскать пересечение кривой Безье и прямой x = x*, где х* - заданный нами аргумент.
Проблема 3: эстетика, или что делать с ограничениями на вершины и направляющие?
В принципе, при решенных проблемах 1 и 2 мы из любого, даже завязанного морским узлом, сплайна сможем выудить график скалярной функции.
Но будет ли понятно пользователю, как он получил результат? Далеко не всегда. Выбрасывание кусков из кривой с целью получить график функции может изменить кривую до неузнаваемости, и не будет даже понятно, за какую ее часть отвечает контрольная точка или направляющая. Для получения понятных результатов необходимо ввести эвристические ограничения на позиции конрольных точек и направляющих.
Решение проблемы единственности значения для функции
Для начала cделаем так, чтобы кривая, состоящая лишь из Corner-точек, не имела несколько Y для 1 значения X.
Для этого просто будем сортировать точки кривой по X.
Кривые могут иметь самопересечения по оси X, т.е. для одному х могут соответствовать более 1 одной точки кривой.
Введем ограничение: пусть имеем график функции, полученный из кривых Безье. Тогда для любых 2 контрольных точек P1P2 значения Y для любого X из интервала [P1.x, P2.x) принадлежат кривой Безье, построенной по точкам P1P2. Это означает, что соседние кривые не могут залезать на полосу [P1.x, P2.x).
На основании этого ограничения самопересечения кривых можно условно разделить на 2 категории:
Если кривая "залезает" на интервал соседних контрольных точек, мы получаем "провисающую" петлю.
У кубической кривой, в отличие от квадратичной, могут быть еще и самопересечения и S-образные петли. Такие петли могут лежать и внутри интервалов конрольных точек.
С петлями желательно бороться. Желательно, но не обязательно потому, что нижеизложенный алгоритм нахождения значения функции справляется и с кривым при наличии петель. При желании проблему с петлями можно решить при помощи ограничений в интерфейсе - этому повящено решение Проблемы 3.
Для удаления провисающих петель достаточно найти параметры кривой, соответствующие точкам пересечения кривой с прямыми x=P1.x и x=P2.x, разрезать кривые в данных точках и отбросить провисающие части кривой. Как найти пересечение кривой и прямой x=x* описано в следующем пункте статьи.
Для удаления петель на кубических кривых я использовал следующую эвристику.
Построим неадаптивную ломаную по кривой - «пройдемся» по ней, изменяя параметр кривой t с некоторым небольшим шагом и получим множество точек
P[i] i=1..N
и соответсвующим им параметрам кривой
T[i] i=1..N
Пойдем по кривой с 2 противоположных концов по индексам i=1..N и j = N..1
На каждой итерации цикла:
Увеличиваем i пока P[i].x >= P[i + 1].x (т.е. пока кривая не начнет двигаться в обратном направлении)
Уменьшаем j пока P[j+1].x <= P[j].x (т.е. пока кривая не начнет двигаться в обратном направлении)
Если P[i].x >= P[j].x, выходим из цикла
Если в итоге цикл завершился не по условию P(i).x >= P(i).x, то кривая «визуально» не имеет самопересечений и S-петель. «Визуально» - потому как алгоритм может не поймать мелкие петли. Но мелкие петли волновать нас не будут – они не вносят проблем в отображение, а алгоритм нахождения значения функции по заданному аргументу, описанный в следующем пункте, может справится даже с кривыми с наличием петель.
Если же такая пара существует, то
находим общий отрезок для проекций отрезков на оси Х. Пусть x* - середина этого отрезка.
Находим si, sj - линейные параметры на отрезках P(i)P(i+1) и P(j+1)P(j), соответствующие x*, т.е
Это – линейные интерполяции параметров кубической кривой, соответствующих точкам по полученным si и sj.
Код вышеописанных пунктов:
// Структура для хранения пар точка кривой-параметр, соответсвующий данной точкеstruct ApproxPoint
{
// Координаты точки
Vec2 pos;
// Параметр кривой, соответствующий данной точкеdouble t;
// Простой конструктор от 2 параметров здесь опущен
};
// Находит линейный параметр для пересечения отрезка x1x2 в точке xdouble get_lerp_param(constdouble x1, constdouble x2, constdouble x)
{
ASSERT(x1 <= x2)
ASSERT(x1 <= x && x <= x2)if(x2 - x1 > std::numeric_limits<double>::epsilon())return(x - x1) / (x2 - x1);
elsereturn0.5;
}
// Находит приближенное значение параметра кривой Безье для заданного x по 2 точкам аппроксимацииdouble get_curve_intersection_t(const ApproxPoint& point1, const ApproxPoint& point2, constdouble x)
{
return lerp(point1.t, point2.t, get_lerp_param(point1.pos.x, point2.pos.x, x)));
}
// Находит ортогональную проекцию вектора vector на ось axisstaticdouble get_orthogonal_projection(const Vec2& vector, const Vec2& axis)
{
return vector.x * axis.y - vector.y * axis.x;
}
// Проверка на отсутствие самопересечений или S петельstaticbool cubic_curve_has_no_inner_loops(const CurbicPoints& points)
{
Vec2 vec1 = points.p1- points.p2
Vec2 middle_vec = points.p3 - points.p2;
Vec2 vec2 = points.p4 - points.p3;
// Если углы p1p2p3 и p2p3p4 - тупые и точки p1 и p4 лежат по одну сторону от p2p3, возращаем truereturn(get_orthogonal_projection(vec1, middle_vec) * get_orthogonal_projection(vec2, middle_vec)>= 0.0 &&
vec1.dot(middle_vec)<= 0.0 && vec2.dot(middle_vec)>= 0.0);
}
// Удаляет петли кубической кривой, если они есть. Кривая задается точками points// В range1 и range2 возращаются параметры конца и начала частей кривой после отбрасывания части, // образующей петлю. // Если петли не были выброшены, range1 и range2 не используется// Функция возращает true, если кривая имела петли и они были удалены, false в противном случаеbool remove_cubic_curve_loops(const CurbicPoints& points, Vec2& range1, Vec2& range2)
{
// Число отрезков, которое содержит аппроксимация кривойstaticconstunsigned SUBDIVISIONS = 100;
staticdouble delta = 1.0 / SUBDIVISIONS;
if(cubic_curve_has_no_inner_loops(point1, point2, point3, point4))returnfalse;
// Создаем набор точек
std::vector<ApproxPoint> points;
double t = 0.0;
while(t <1.0)
{
Vec2 point = get_cubic_bezier_point(points, t);
points.push_back(ApproxPoint(point, t));
t += delta;
}
points.push_back((points.p4, 1.0));
size_t left_i = 0;
size_t right_i = points.size() - 1;
while(left_i < right_i - 2)
{
ApproxPoint min1 = points[left_i];
ApproxPoint max1 = points[left_i + 1];
ApproxPoint min2 = points[right_i];
ApproxPoint max2 = points[right_i - 1];
if(min1.pos.x > max1.pos.x)
std::swap(min1, max1);
if(min2.pos.x > max2.pos.x)
std::swap(min2, max2);
// Если нашли пересечение отрезков min1,max1 и min2,max2 на оси Xif(!(max1.pos.x < min2.pos.x || min1.pos.x > max2.pos.x))
{
constdouble slice_x = (std::max(min1.pos.x, min2.pos.x) + std::min(max1.pos.x, max2.pos.x)) * 0.5;
range1.x = 0.0;
range1.y = get_curve_intersection_t(min1, max1, slice_x);
range2.x = get_curve_intersection_t(min2, max2, slice_x);
range2.y = 1.0;
returntrue;
}
constbool move_left_ind = points[left_i + 2].pos.x >= points[left_i + 1].pos.x;
constbool move_right_ind = points[right_i - 2].pos.x <= points[right_i - 1].pos.x;
// Если проверка в ASSERT не выполнится, мы будем крутится в цикле вечно
ASSERT(move_left_ind || move_right_ind)if(move_left_ind)
++left_i;
if(move_right_ind)
--right_i;
}
returnfalse;
}
Разрезаем кривую в точках ti_cross и tj_cross. Оставляем куски кривой, соответствующие диапазону параметора оригинальной кривой (0, ti_cross) и (tj_cross, 1).
Из-за численных ошибок у нас с большой вероятностью не совпадут значения x для концов частей кривой. Стянем их, просто назначив одинаковые значения x для точек, которые образуют разрыв.
Из-за неточности эвристики возле места стыковки кусков кривой могут появится «провисающие» петли. Избавляемся от них, просто сдвинув контрольные точки P1 и P2 кривой P0P1P2P3 по х, чтобы они поместились в интервале (P0.x P3.x)
if (p1.x < p0.x)
p1.x = p0.x;
if (p2.x > p3.x)
p2.x = p3.x;
Может показаться, что это очень грубое действие. На самом деле, в контексте данного алгоритма если точки P1 и P2 и выходят из интервала (P0.x P3.x), то очень незначительно, и форма кривой от такого сдвига визуально практически не меняется.
После всех описанных выше манипуляций из любой квадратной или кубической кривой получается кривая или набор кривых, у которых визуально нет самопересечений по оси x.
Решение проблемы нахождения значения функции
Для квадратичных кривых найти пересечение прямой x=x* легко - нужно лишь решить квадратное уравнение c ограничениями:
В реализации есть несколько нюансов:
удобно всегда получать лишь 1 допустимое решение для приведенного выше уравнения; чтобы никогда не получать 2 решения в точках P0 и P1, достаточно просто отступить от них внутрь отрезка P0P1
если расстояние между P0 и P1 мало, могут появляться численные ошибки. Я для такого случая можно стягивать вершины в одну или использовать для данной пары точек линейную кривую вместо квадратичной.
Для кубических кривых все сложнее. Для нахождения t нужно решить уже уравнение 3-ей степени. Нахождение пересечения кривой Безье и прямой x = x* можно свести к решению 1-го полинома 3-ей степени:
Посмотрев в wikipedia аналитическое решение, я отказался от него в пользу численного метода. Какие есть варианты? Есть Bezier Clipping – итеративный метод, использующий некоторые особенности кривых Безье для быстрого нахождения интервала, на котором лежит корень.
Однако, метод использует пересечения выпуклого многогранника и прямой. Collision detection для многогранника и прямой не то чтобы сложен, но он имеет неприятные частные случаи плюс демонстрирует нестабильность в случае малых длин сторон многогранника. Кроме того, многогранник в этом методе – это выпуклая оболочка, построенная по 4 точкам. Ее построение требует дополнительных вычислений. После реализации Bezier Clipping в более-менее стабильном виде я сравнил с его обычной дихотомией – выяснилось, что дихотомия работает быстрее.
Как Bezier Clipping, так и простая дихотомия имеют 1 недостаток - они не могут найти множественные корни. Но у нас такой проблемы не возникает, т.к. мы предварительно обрезаем самопересечения.
Казалось бы, мы решили проблему. Но и тут есть один небольшой нюанс. В проекте, для которых разрабатывался редактор, предъявлялись жесткие условия по гладкости кривой, так как предполагалось семплирование точек на кривой с очень малым шагом по x. Для квадратичных кривых все было хорошо – описанные выше квадратное уравнение решалось с достаточной точностью для семплирования. В случае кубической кривой все было плохо – дихотомия была плоха как по производительности (а нам нужно было вычислять по нескольку значений функции на каждый пиксел просчитываемого изображения), так и по гладкости – для нормальной производительности приходилось увеличивать интервал нахождения корня, что давало артефакты – «ступеньки».
Я использовал подход, описанный в данной статье. Суть его – приближение кубических кривых кусками квадратичных Безье. При малой кривизне углов для конрольных точек кубическая прямая почти совпадает к квадратичной кривой. Естественно, мы уже не получаем «честных» кубических кривых, зато с такой аппроксимацией проще работать, и вычисление точки для кубической кривой сводится к линейному поиску по аппроксимирующим отрезкам (а их нужно в большинстве случаев не очень много) + решению одного квадратного уравнения.
Строится приближение следующим образом:
Разбиваем кубическую кривую на куски, руководствуясь критерием останова:
Условие (*) означает, что углы P0P1P2 и P1P2P3 - тупые
Условие (**) заставляет разбивать кривую до тех пор, пока для каждого куска кривой, обозначенного P0-P1-P2-P3, точки P0 и P3 не будут лежать по одну сторону от отрезка P1P2.
Аппроксимируем каждый из кусков кубической P0-P1-P2-P3 кривой квадратичной кривой P0-P*-P3, где P* находится как пересечение прямых P0-P1 и P2-P3.
Существование пересечения обеспечивают условия (*) и (**). В случае, если для каких-либо кусков кривых "схлопываются" контрольные точки, то такие частные случаи решаются так же, как и в случае построения адаптивной ломаной (смотрите выше).
В итоге для нахождения пересечений с кубической кривой я использовал метод аппроксимации квадратичными кривыми. Дихотомию я использовал в процессе разработки отрезания "провисающих" петель у кубических кривых, но введя впоследствии ограничения на положение направляющих, я ликвидировал варианты кубических кривых с "провисающими петлями". Так что дихотомию для кубических кривых привожу лишь для полноты картины - вдруг кому-нибудь понадобится.
Ограничения для точек и направляющих в интерфейсе.
Для точек ограничение заключается в их отсортированности по оси х. Т.е. сдвиге точки так, что она попадает в новый интервал между 2 точками, порядок контрольных точек изменяется.
Для направляющих имеются следующие ограничения:
Для связанных направляющих (тип точек – Smooth, Symmetrical) ограничений нет, однако, при переходе одного из концов направляющей через контрольную точку кривая перестраивается, так для каждой пары контрольных точек кривая строится по точкам «позиция P1»-«правая направляющая P1»-«левая направляющая P2»-«позиция P2»
Для несвязанных направляющих (тип точек – Cusp) одна из направляющих всегда находится слева от контрольной точки, другая – всегда справа.
При наличии данных ограничений на направляющие у нас не будут получатся петли-«провисания» у кубических кривых (у квадратичных они сохранятся).
Можно было бы ввести еще и ограничение:
Т.е. чтобы направляющие не заходили за соседние контрольные точки. В этом случае мы избавимся от всех возможных петель:
"Провисающих петель" - по свойству, что выпуклая оболочка из контрольных точек содержит кривую.
Самопересечений и S-образный петель у кубических кривых - к сожалению, я не нашел и не сделал строгого доказательства данного факта. Визуально это так.
Но данное ограничение достаточно сильно уменьшает возможности пользователя в движении направляющих, что не очень удобно. По этой причине я не применял его.
Заключение
Описанные мной приемы работы с кривыми Безье можно использовать в различных комбинациях как строительные блоки. Возможно, вам не потребуется высокая гладкость кривой или большая свобода положений направляющих - в этом случае можно аппроксимировать кривые ломаными или не волноваться о петлях, введя более жесткие ограничения в интерфейсе. Как всегда, не существует идеального "рецепта приготовления".
Сопроводительного кода, к сожалению, я не могу предоставить - проект коммерческий. Примеры в данной статьи - это упрощенный код. Работоспособность примеров не тестировалась, так что если будете их использовать и найдете ошибки в коде - напишите, пожалуйста, в комментариях.
Сам редактор в действии вы можете посмотреть, зайдя на http://filterforge.com/, и скачав trial свежей версии (на момент написания статьи наиболее свежая версия - 2-я бета версии 3.0)/
