Уже год как я сменил работу на новую. В этой статье я хочу поделиться опытом, накопленным на прошлом месте. Здесь рассмотрены методы аппроксимации кривых Безье, а также обработка исключительных случаев, при которых простые алгоритмы показывают себя не очень хорошо. Все, кому близка тема векторной графики прошу под кат.

Введение

На предыдущей работе я был оператором ЧПУ-станка, делал мебельные фасады из МДФ. Работа разнообразная: то сидишь за компом и чертишь в CAD/CAM-программе, то стоишь контролируешь резку, то приносишь заготовки и уносишь готовые детали.

До того, как я пришёл в фирму, там был довольно скудный каталог рисунков, которого явно не хватало клиентам, поэтому они всё время просили сделать как на фото из интернета. Все специфичные рисунки я рисовал вручную, стандартные же были довольно просты, а так как я увлекался программированием, то решил генерировать их G-код для станка напрямую без участия CAD/CAM-программы.

В интернете существует пара вариантов в виде макросов для CorelDRAW и что-то ещё в этом духе. Я же решил придумать более лёгкое, быстрое и кроссплатформенное решение, которое решило бы мои текущие задачи максимально эффективно. Выбор был сделан в пользу консольной программы с выводом результата в stdout.

В данной статье я не буду затрагивать саму программу, поскольку это тема отдельной статьи. Я расскажу о том, с чем мне пришлось столкнуться в самом начале, а именно как представить кривую Безье в виде прямых отрезков с точностью, при которой будет создаваться минимальное количество промежуточных точек при сохранении визуальной плавности на изделии. Станок очень хорошо передаёт издаваемыми звуками то, насколько плавно он движется по траектории. Я ориентировался на звук при резке чертежей, созданных профессиональной CAD/CAM-программой.

На тот момент самой дельной информацией по этому вопросу была вот эта статья на RSDN. Статья заслуживает внимания. В ней много скриншотов демонстрационной программы и её исходный код. Т.е. можно было не только почитать статью, но и запустить всё это у себя на компьютере. Перед написанием статьи я зашёл на страницу и обнаружил, что ссылки на программу и её исходный код мёртвые, поэтому теперь можно только почитать.

С чем мы собираемся бороться?

Итак, начнём. На более прямых участках нужно минимальное количество точек, но в определённых местах, таких как крутые повороты, перегибы и узкие петли, кривая очень сильно изгибается, и чтобы сохранить плавность, нужно генертровать больше точек. Потому разбиение и называется адаптивным. Ниже примеры кривых, которые обычно сложно аппроксимировать.

Рис. 1 Крутой поворот

Рис. 2 Перегиб

Рис. 3 Узкая петля

В общем-то внешне это примерно одно и то же. И тут мой подход отличается от подхода автора той статьи. Он решил обрабатывать все исключительные ситуации в отдельности. Как итог нагромождение кода, который тем не менее не особо-то эффективно решает все проблемы. Я начал искать решение, которое было бы простым и универсальным. Лень двигатель прогресса. В итоге все проверки сводятся к одной, а степень детализации задаётся одним параметром.

Кривые 2-го порядка

Первый пример будет с кривой Безье 2-го порядка. Начальная и конечная точки эталонной кривой совпадают с таковыми у исходной. Опорная точка находится посередине, т.к. кривая 2-го порядка, а значит мы делим её на две части.

Рис. 4 Кривая Безье 2-го порядка

Расстояние d1 на рисунке показывает, на сколько исходная опорная точка удалена от эталонной. Если d1 меньше заданной нами величины, чертим прямую линию из точки 0 в точку 2 и заканчиваем аппроксимацию, если нет делим кривую пополам и рекурсивно для каждой половины повторяем проверку. Извините frontend-разработчики, джаваскрипта не будет. Ниже пример кода на питоне.

from math import sqrtdef b2(x0, y0, x1, y1, x2, y2, d):    mx1 = x1 - (x0 + x2) / 2    my1 = y1 - (y0 + y2) / 2    d1 = sqrt(mx1 ** 2 + my1 ** 2)    if d1 < d:        print(x2, y2)    else:        x01 = (x0 + x1) / 2        y01 = (y0 + y1) / 2        x12 = (x1 + x2) / 2        y12 = (y1 + y2) / 2        x012 = (x01 + x12) / 2        y012 = (y01 + y12) / 2        b2(x0, y0, x01, y01, x012, y012, d)        b2(x012, y012, x12, y12, x2, y2, d)

Кривые 3-го порядка

Второй пример будет с кривой Безье 3-го порядка. Начальная и конечная точки эталонной кривой совпадают с таковыми у исходной. Опорные точки находятся на 1/3 и 2/3 её длины, т.к. кривая 3-го порядка, а значит мы делим её на три части.

Рис. 5 Кривая Безье 3-го порядка.

Аналогично предыдущему примеру, расстояния d1 и d2 на рисунке показывают отклонение опорных точек от их эталонных значений. Если оба расстояния d1 и d2 меньше заданной нами величины, чертим прямую линию из точки 0 в точку 3 и заканчиваем аппроксимацию, если нет делим кривую пополам и рекурсивно для каждой половины повторяем проверку. Пример кода прилагается.

from math import sqrtdef b3(x0, y0, x1, y1, x2, y2, x3, y3, d):    px = (x3 - x0) / 3    py = (y3 - y0) / 3    mx1 = x1 - x0 - px    my1 = y1 - y0 - py    mx2 = x2 - x3 + px    my2 = y2 - y3 + py    d1 = sqrt(mx1 ** 2 + my1 ** 2)    d2 = sqrt(mx2 ** 2 + my2 ** 2)    if d1 < d and d2 < d:        print(x3, y3)    else:        x01 = (x0 + x1) / 2        y01 = (y0 + y1) / 2        x12 = (x1 + x2) / 2        y12 = (y1 + y2) / 2        x23 = (x2 + x3) / 2        y23 = (y2 + y3) / 2        x012 = (x01 + x12) / 2        y012 = (y01 + y12) / 2        x123 = (x12 + x23) / 2        y123 = (y12 + y23) / 2        x0123 = (x012 + x123) / 2        y0123 = (y012 + y123) / 2        b3(x0, y0, x01, y01, x012, y012, x0123, y0123, d)        b3(x0123, y0123, x123, y123, x23, y23, x3, y3, d)

Выведение

Поскольку массово используются только кривые этих порядков, копать дальше считаю лишним. Этот код несколько лет помогал выполнять мне мою работу. Требования к плавности на станке высокие, вибрация его убивает. Станок жив до сих пор, поэтому код прошёл боевое крещение. При работе на станке все мои размеры были в миллиметрах, переменная d в моём случае была 0.025 мм.

Также есть материал по рациональным кривым тех же порядков (это те, что с весами). Если будет интерес, про это тоже напишу. Благодарю за прочтение.

источник изображения

Возможно, вы видели предыдущий пост, где были предоставлены визуализации первых 1000 цифр и . Он возник в результате небольшого спора о том, лучше ли , чем . По этому поводу идут бесконечные дебаты, и я подумал, что могу пошутить по этому поводу. В этом посте я хочу показать, как создать визуализации, и надеюсь, что вы захотите попробовать удивительный пакет Luxor.jl после прочтения. Вчера я начал читать туториал, и это потрясающе! В прошлый раз визуализация делалась на Javascript, и я подумал, что этот аккуратный маленький проект сойдет, чтобы начать изучать Луксор. Как уже упоминалось в let me be your mentor: я думаю, что очень важно иметь такие маленькие проекты, чтобы освоить новый инструмент.

Основная идея

Я хотел воссоздать визуализацию, которую видел в Numberphile от Мартина Крживинского.

Там был круг (который, вполне ассоциируется и с и с ) разделенный на 10 сегментов, по одному для каждой цифры. Цифры нашего иррационального числа представляются кривыми внутри этого круга, так что 3.1415 (я начинаю с 14) это кривая от сегмента 1 до сегмента 4, а затем обратно к 1, потом до 5 и так далее. Каждый раз мы перемещаемся немного по часовой стрелке в сегменте так, что 14 создает различные кривые (в зависимости от текущего положения, в котором мы находимся).

Потом надобавляем всякие фичи. Мы должны начать чувствовать себя комфортно с Луксором. Важно: не надо искать математическую интерпретацию это просто небольшой проект визуализации ;)

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

Начинаем

using Luxorfunction vis()    @png begin        sethue("black")        circle(O, 50, :stroke)        setdash("dot")        circle(O, 70, :stroke)        sethue("darkblue")        circle(O, 10, :fill)    end 500 200 "./start.png"end

вызываем vis() и создаем файл start.png который будет выглядеть как-то так:

Давайте быстренько пройдемся по командам:

@png beginend width height "filename.png"

просто хороший макрос. :)

sethue задает цвет и принимает либо строку, как показано выше или цвет пакета из Colors. Он устанавливает цвет для следующих команд рисования до тех пор, пока вы не выберете другой. То же самое верно и при установке ширины линии с помощью setline, или при установке размера шрифта, или при других общих настройках.

Команды рисования, такие как circle, обычно принимают некоторые параметры и заканчиваются параметром действия, таким как :stroke или :fill.

О это буква "О", а не число "0". :) Она представляет собой начало координат и является краткой формой для Point(0, 0). В Луксоре начало находится в центре полотна. В качестве второго параметра должен быть задан радиус.

Давайте сначала нарисуем внешний круг и добавим цифры:

radius = 100@png begin    background("black")    sethue("white")    circle(O, radius, :stroke)    for i in 0:9         = 2*0.1*i+0.1*        mid = Point(            radius*sin(),            -radius*cos(),        )        label(string(i), :N, mid)    endend 700 300 "./first_step.png"

Первая часть должна быть достаточно простой.

 = 2*0.1*i+0.1*

возможно, это не идеально написано (кроме того, я мог бы использовать :D). 2*0.1*i начинает с северного положения, а затем для следующего i происходит перемещение на . Я добавляю "0.1 ", потому что хочу переходить к середине каждого сегмента. Может быть, следует написать 0.5/10*2. Затем мы просто поворачиваем наш холст и двигаясь чуть выше радиуса, рисуем метки. На самом деле такое можно проделать в Luxor, используя rotate и translate. Но я решил сделать вручную, так как мне все равно это пригодится позже. В общем формула такова:

Такое преобразование поворачивает плоскость на и производит трансляцию на x,y. Поскольку я перевожу только на y, мне не нужно первое тождество. Помните, что y увеличивается, когда идет вниз.

В настоящее время есть две проблемы:

• на самом деле нам не нужен круг, нам нужны дуги (сегменты) для каждой цифры
• подписи не читаются

Команда label принимает три значения: текст, вращение и положение, где вращение может быть записано как :N,: E,: S,: W для севера, востока, юга, запада или как угол (в радианах). :N есть . Поэтому мы хотим начать с , а потом добавлять текущий угол поворота. Кроме того, смещение было бы здорово, если бы оно не доставало непосредственно до окружности или не подходило слишком близко к ней. Здесь мы могли бы увеличить радиус или использовать ;offset в команде label.

Для первой задачи нам нужна функция arc2r, которая принимает три аргумента
c1, p1, p2 + действие: c1 это центр окружности, а p1 и p2 точки на окружности, между которыми должен быть показан сегмент. По умолчанию выбрано направление по часовой стрелке.

Мы определяем следующую функцию, чтобы получить и соответствующую точку более простым способом:

function get_coord(val, radius)     = 2*0.1*val    return Point(        radius*sin(),        -radius*cos(),    )end

а потом:

background("black")for i in 0:9    from = get_coord(i, radius)    to = get_coord(i+1, radius)    randomhue()     = 2*0.1*i+0.1*    mid = Point(        radius*sin(),        -radius*cos(),    )    label(string(i), -/2+, mid; offset=15)    move(from)    arc2r(O, from, to, :stroke)end

Я использовал randomhue, чтобы получить случайный цвет. Мы исправим это в следующий раз :)
Также я переставлял порядок Label и arc2r и поставил move, так как в противном случае линии рисуются от метки дуги. Это происходит потому, что arc продолжает текущий путь.

Выглядит намного лучше! Давайте возьмем несколько хороших цветов из Colorschemes.jl.

Я использовал схему rainbow, начиная с 7-го цвета :D. Вы, возможно, захотите испытать другие цветовые схемы, так как здесь цвета не так легко различить, но мне все равно почему-то нравится именно она.

using ColorSchemescolors = ColorSchemes.rainbow[7:end]

и затем

sethue(colors[i+1])

помните, что индексация массивов в Julia начинается с единицы.

Каковы следующие шаги?

• Добавление строк
• Рефакторинг кода
• Оживление процесса
• Добавление точек
• Добавление гистограммы сверху

Я думаю, что визуально привлекательно иметь круг посередине, где мы можем добавить символ (или ) позже.
Поэтому мы не можем провести прямые линии от одного сегмента к другому. Для этого я использую квадратичные кривые Безье.

Давайте сначала получим цифры числа Пи:

max_digits = 10digits = setprecision(BigFloat, Int(ceil(log2(10) * max_digits+10))) do    return parse.(Int, collect(string(BigFloat(pi))[3:max_digits+2]))end

это дает нам первые 10 цифр после десятичной точки числа Пи. Для этого мне нужно установить точность BigFloat. Довольно интересно, что пи не является жестко закодированной константой в Джулии. Оно вычислено таким образом, что я в принципе могу получить любую точность, какую захочу. Точность должна быть задана в количестве битов, так что необходимо выполнить небольшое вычисление. Я добавил +10 в конце, чтобы быть уверенным :D

Чтобы нарисовать квадратичную кривую Безье, нам нужны три точки. Начало, конец и контрольная точка. В качестве контрольной точки я выбираю точку на внутреннем круге, который просто также разделен на десять сегментов, и выбираю сегмент, который находится посередине между текущей цифрой from_val и следующей цифрой to_val.

Я должен уточнить, что я имею в виду под серединой: средняя точка между 0 и 4 должна быть 2, но между 8 и 0 она должна быть 9. Она определяется кратчайшим путем от одного сегмента к другому, а потом берется середина.

Кроме того, у меня на самом деле нет 10 дискретных сегментов, это просто для понимания. Я могу использовать среднюю точку 1,23 или что-то в этом роде. Это используется, потому что мы меняем нашу начальную и конечную позиции на основе текущей позиции, которую мы находимся в нашем массиве цифр.

Я надеюсь, что все станет яснее, ели взглянуть на код:

small_radius = 70for i in 1:max_digits-1    from_val = digits[i]    to_val = digits[i+1]    sethue(colors[from_val+1])    f = from_val+(i-1)/max_digits    t = to_val+i/max_digits    from = get_coord(f, radius)    to = get_coord(t, radius)    # get the correct mid point for example for 0-9 it should be 9.5 and not 4.5    mid_val = (f+t)/2    mid_control = get_coord(mid_val, small_radius)    if abs(f-t) >= 5        mid_control = get_coord(mid_val+5, small_radius)    end    pts = Point[from, mid_control, mid_control, to]    bezpath = BezierPathSegment(pts...)    drawbezierpath(bezpath, :stroke, close=false)end

Думаю, уже выглядит достаточно хорошо. Цвета линий подгоняются под цвета из под цифр. Итак, в какой-то момент мы переходим от 9 к 2. Вместо этого я хотел бы посмотреть, куда мы идем и откуда идем. Это можно сделать с помощью blend и setblend. Это линейная смена цвета "от" и "до", так что на самом деле не по кривой, но я думаю, что она достаточно хороша.

setblend(blend(from, to, colors[to_val+1], colors[from_val+1]))

Это похоже на sethue поэтому нам нужно задать его в какой-то момент, прежде чем мы вызовем drawbezierpath.

Давайте добавим еще несколько цифр и немного уменьшим ширину линии: setline(0.1)

Ладно я думаю что внутренний радиус немного велик:

small_radius = 40

Затем мы можем добавить в середине, прежде чем немного очистить код, чтобы создать нашу первую анимацию.

Luxor.jl не поддерживает латексные стринги LaTeXStrings.jl это облом, но мы можем использовать UnicodeFun.jl.

using UnicodeFuncenter_text = to_latex("\\pi")

и промеж циклов ставим:

sethue("white")fontsize(60)text(center_text, Point(-2, 0), valign=:middle, halign=:center)

Мне кажется Point(-2, 0) более центральная, чем Point(0, 0) или O.

Анимация

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

В Луксоре это можно сделать с помощью функции animate, которая берет несколько сцен и их номера кадров. Это также обеспечит немного большую структуру кода.

У нас может быть сцена для устойчивого фона и одна для линий.

Прежде чем мы напишем функцию, давайте определим очень короткую анимацию, чтобы увидеть, как это делается.

function draw_background(scene, framenumber)    background("black")endfunction circ(scene, framenumber)    setdash("dot")    sethue("white")    translate(-200, 0)    @layer begin         translate(framenumber*2, 0)        circle(O, 50, :fill)    endendfunction anim()    anim = Movie(600, 200, "test")    animate(anim, [        Scene(anim, draw_background, 0:200),        Scene(anim, circ, 0:200),    ],    creategif = true,    pathname = "./test.gif"    )end

Сначала мы создаем Movie с width, height и name.
Затем мы вызываем animate с помощью созданного Movie и списка scenes, а затем функции и диапазон кадров, начинающихся с 0.

Происходит вызов draw_background(сцена, 0) и circ(scene, 0) для первого кадра. Сцена может содержать некоторые аргументы, которые мы будем использовать для нашей анимации. Остальное в основном так же, как и раньше, просто мы можем, конечно, использовать переменную framenumber.

Теперь я разделю все это дело на функции и определю переменные, такие как цифры, которые мы хотим визуализировать, чтобы нам было легче визуализировать или другие вещи.

using Luxor, ColorSchemesusing UnicodeFunfunction get_coord(val, radius)     = 2*0.1*val    return Point(        radius*sin(),        -radius*cos(),    )endfunction draw_background(scene, framenumber)    background("black")    radius = scene.opts[:radius]    colors = scene.opts[:colors]    center_text = scene.opts[:center_text]    for i in 0:9        from = get_coord(i, radius)        to = get_coord(i+1, radius)        sethue(colors[i+1])         = 2*0.1*i+0.1*        mid = Point(            radius*sin(),            -radius*cos(),        )        label(string(i), -/2+, mid; offset=15)        move(from)        arc2r(O, from, to, :stroke)    end    sethue("white")    fontsize(60)    text(center_text, Point(-2, 0), valign=:middle, halign=:center)endfunction dig_line(scene, framenumber)    radius = scene.opts[:radius]    colors = scene.opts[:colors]    center_text = scene.opts[:center_text]    bezier_radius = scene.opts[:bezier_radius]    max_digits = scene.opts[:max_digits]    digits = scene.opts[:digits]    setline(0.1)    for i in 1:min(framenumber, max_digits-1)        from_val = digits[i]        to_val = digits[i+1]        f = from_val+(i-1)/max_digits        t = to_val+i/max_digits        from = get_coord(f, radius)        to = get_coord(t, radius)        # get the correct mid point for example for 0-9 it should be 9.5 and not 4.5        mid_val = (f+t)/2        mid_control = get_coord(mid_val, bezier_radius)        if abs(f-t) >= 5            mid_control = get_coord(mid_val+5, bezier_radius)        end        pts = Point[from, mid_control, mid_control, to]        bezpath = BezierPathSegment(pts...)        # reverse the color to see where it is going        setblend(blend(from, to, colors[to_val+1], colors[from_val+1]))        drawbezierpath(bezpath, :stroke, close=false)    endendfunction anim()    anim = Movie(700, 300, "test")    radius = 100    bezier_radius = 40    colors = ColorSchemes.rainbow[7:end]    max_digits = 1000    center_text = to_latex("\\pi")    digits_arr = setprecision(BigFloat, Int(ceil(log2(10) * max_digits+10))) do        return parse.(Int, collect(string(BigFloat(pi))[3:max_digits+2]))    end    args = Dict(:radius => radius,        :bezier_radius => bezier_radius,        :colors => colors, :max_digits => max_digits,        :digits => digits_arr, :center_text => center_text    )    animate(anim, [        Scene(anim, draw_background, 0:max_digits+50, optarg=args),        Scene(anim, dig_line, 0:max_digits+50, optarg=args),    ],    creategif = true,    pathname = "./pi_first.gif"    )end

Единственное, что я еще не объяснил, это optarg в функции Scene и получение его с помощью radius = scene.opts[:radius].

Мы как бы потеряли возможность создавать простые образы. Поэтому я создал структуру

struct PNGScene    opts::Dict{Symbol, Any}end

и использую некоторые аргументы в функции anim, которую я переименую в viz :D

Тогда я могу использовать что-то вроде:

scene = PNGScene(args)@png begin    draw_background(scene, max_digits)    dig_line(scene, max_digits)end 700 300 "./$fname.png"

Не волнуйтесь, в конце есть репка, где вы можете увидеть весь код целиком. Просто немного сложно описать здесь изменения.

Может, мне стоило снять видео? :D

Добавление точки Фейнмана

Мы визуализировали соединение цифр с цифрами с помощью кривых, но если бы у нас встретилось что-то вроде 555 в цифрах, мы бы видели только линию, идущую в направлении центра и обратно (или, может быть, мы видим две в зависимости от наших максимальных цифр, разрешения и т. д.)

Вместо этого мы можем показать дополнительную точку всякий раз, когда это происходит. Это можно получить благодаря аргументу функции show_dots, что вы можете найти в моем коде ;)

Я только что проверил длину последовательности, и когда она больше 1, я рисую круг, где это происходит, и цвет это цифра после этой последовательности. Большой круг в сегменте 9 это так называемая точка Фейнмана, где цифра 9 появляется 6 раз в позиции 762.

Добавление гистограмм

Последняя вещь в моем списке получить гистограмму на каждом сегменте, чтобы показать, случаются ли некоторые комбинации пар чаще, чем другие.

Для этого я использую функцию poly с четырьмя точками. В идеале, она должна быть ограничена двумя дугами, а не двумя линиями, но я оставляю это читателю :)

Тау

Да, можно было бы в принципе сгенерировать случайное число с 1000 цифрами и получить аналогичный результат...

Простое число

В двух словах: использование нашей функции для визуализации большилства элементов не так разумно, но так или иначе может получится что-то интересное.

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

Тем не менее, гистограммы становятся все интереснее, я думаю:

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

Давайте сосредоточимся на гистограммах и визуализируем простые числа под 10 000 000:

Узор сохраняется.

Код

Окай, тут у нас репка

Я хотел бы создать что-то вроде штучек, из 3b1b.

По крайней мере, небольшие простые версии с некоторыми удобными функциями визуализации :)

Спасибо за чтение и особая благодарность моим 10 покровителям!

Я буду держать вас в курсе событий на Twitter OpenSourcES и на более личном:
Twitter Wikunia_de

Article

Вы сталкивались когда-нибудь с построением (непрерывного) пути обхода кривой на плоскости, заданной отрезками и кривыми Безье?

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

Всё было хорошо, пока из-под рук дизайнеров не стали вылезать монструозные пути, где отдельные кривые могли пересекаться или не точно состыковываться. Объяснение было предельно простым визуально они все лежат как надо, а для станка, который этот путь будет обходить, такие отклонения незаметны.

Вооружившись знанием о величине максимально допустимого отклонения, я приступил к исследованию, результатами которого хочу поделиться.

Первое что я сделал это разобрался как на сегодняшний день (октябрь 2020) обстоят дела с поиском точек пересечения кривых. То ли я не там искал, то ли не то спрашивал, но найти простого решения не получилось. Хотя, идея с результантом пары полиномов довольно занимательна. Много разных алгоритмов связанных с кривыми Безье собрано здесь.

Что не понравилось в известных способах и что точно не хочется делать, так это численно искать корни полиномов, или даже решать квадратичные уравнения. Очень не хочется исследовать кривые на экстремумы. Да и вообще, хотелось бы избежать деления, возведения в степень и всего того, что может привести к неопределённому поведению.

Итак, с чем придётся работать.

• Точки задаются типом Point, например так:

  
  

  using Point = std::array<double, 2>;
  

  
  

  Для Point определены операторы сложения, вычитания, умножения на скаляр, скалярного умножения.

  
  

• Задана величина R допустимого отклонения точек.

  
  

• Кривые заданы массивами опорных (контрольных) точек std::vector<Point>.

  
  

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

  
  

Первое, что точно понадобиться, это простой способ вычисления значения параметрической кривой. (Простой в смысле реализации и читабельности):

Point point(const std::vector<Point> &curve, double t, int n, int i){    return n == 0 ? curve[i] : (point(curve, t, n - 1, i - 1) * (1 - t) + point(curve, t, n - 1, i) * t);}

Оставлять функцию в таком виде для постоянного использования не стоит лучше спрятать её подальше, а пользоваться такой:

Point point(const std::vector<Point> &curve, double t){    return point(curve, t, curve.size() - 1, curve.size() - 1);}

Здесь, curve контейнер для опорных точек: для отрезка их две, для кривой Безье три или четыре или более.

Второе точки надо как-то сравнивать, с учётом R:

template <>struct std::less<Point>{    bool operator()(const Point &a, const Point &b, const double edge = R) const    {        for (auto i = a.size(); i-- > 0;) {            if (a[i] + edge < b[i])                return true;            if (a[i] > b[i] + edge)                return true;        }        return false;    }};

Для поиска точек пересечения пары кривых я воспользовался идеей деления каждой кривой на две до тех пор пока есть пересечение области значений этих кривых. А чтобы не заморачиваться с экстремумами и интервалами монотонности, использовал тот факт, что кривая Безье ограничена выпуклым многоугольником, вершинами которого являются опорные точки кривой. Или даже ещё проще достаточно ограничивающего эти точки прямоугольника:

struct Rect{    Point topLeft, bottomRight;    Rect(const Point &point);    Rect(const std::vector<Point> &curve);    bool isCross(const Rect &rect, const double edge) const    {        for (auto i = topLeft.size(); i-- > 0;) {            if (topLeft[i] > rect.bottomRight[i] + edge                || bottomRight[i] + edge < rect.topLeft[i])                return false;        }        return true;    }};

Алгоритм поиска рекурсивный и достаточно простой

void find(const std::vector<Point> &curveA, const std::vector<Point> &curveB,          double tA, double tB, double dt){

    if (m_isSimilarCurve)        return;

    Rect aRect(curveA);    Rect bRect(curveB);    if (!aRect.isCross(bRect, R))        return;

    if (isNear(aRect.tl, aRect.br, R / 2) && isNear(bRect.tl, bRect.br, R / 2)) {        // 3.1 Для найденного пересечения сохранить наиболее близкие концы кривых        addBest(curveA.front(), curveA.back(), curveB.front(), curveB.back(), tA, tB, dt);        m_isSimilarCurve = (m_result.size() > curveA.size() * curveB.size());        return;    }

    const auto curveALeft = subCurve(curveA, 0, 0.5);    const auto curveARight = subCurve(curveA, 0.5, 1.0);    const auto curveBLeft = subCurve(curveB, 0, 0.5);    const auto curveBRight = subCurve(curveB, 0.5, 1.0);

    const auto dtHalf = dt / 2;    find(curveALeft, curveBLeft, tA, tB, dtHalf);    find(curveALeft, curveBRight, tA, tB + dtHalf, dtHalf);    find(curveARight, curveBLeft, tA + dtHalf, tB, dtHalf);    find(curveARight, curveBRight, tA + dtHalf, tB + dtHalf, dtHalf);

}

Вот тут-то и выполз самый главный вопрос: как найти опорные точки кривой, которая является частью исходной кривой в интервале t от t1 до t2?

После исследования к чему приводит подстановка t = (t2 - t1) t' + t1 я обнаружил простую закономерность. Первая опорная точка вычисляется по исходной кривой при t = t1, последняя при t = t2. Это логично, так как по свойствам кривых Безье (полиномов Бернштейна) крайние точки лежат на кривой. Для остальных точек нашлось простое правило: если в процессе вычисления на шаге k вместо t2 подставить t1, то в результате получим опорную точку под номером k:

Point point(const std::vector<Point> &curve, double t1, int n, int i, double t2, int k){    return n > k ? (point(curve, t1, n - 1, i - 1, t2, k) * (1 - t2) +                    point(curve, t1, n - 1, i, t2, k) * t2)                 : point(curve, t1, n, i);}

Эту функцию тоже лучше спрятать куда подальше, а для получения всех опорных точек воспользоваться этой:

std::vector<Point> subCurve(const std::vector<Point> &curve, double t1, double t2){    std::vector<Point> result(curve.size());    for (auto k = result.size(); k-- > 0;) {        result[k] = point(curve, t1, curve.size() - 1, curve.size() - 1, t2, curve.size() - 1 - k);    }    return result;}

Вот, собственно, и всё.

Примечания.

1. t1 и t2 могут быть любыми:
   
   • subCurve(curve, 1, 0) даст кривую, которая "движется" от конечной точки curve к начальной,
   • subCurve(curve, 1, 2) экстраполирует curve за пределами последней опорной точки.
2. Реализации некоторых методов опущены намеренно, так как не содержат ничего особенно интересного.
3. Функция point(curve, t) не подходит для вычисления множества точек на кривой (например для растрезации), для этого лучше подойдёт вычисление с помощью треугольной матрицы.