Нахождение области замкнутого 2d-однородного кубического B-сплайна

У меня есть список 2d точек, которые являются управляющими вершинами (Dx) для замкнутого однородного кубического B-сплайна. Я предполагаю простую кривую (несамопересекающиеся, все контрольные точки различны).

Я пытаюсь найти область, заключенную кривой:

Если я вычислил узловые точки (Px), я могу рассматривать кривую как многоугольник; то мне "просто" нужно найти оставшиеся области дельты для каждого сегмента, между фактической кривой и прямой, соединяющей узловые точки.

Я понимаю, что форма (и, следовательно, площадь) Bspline инвариантна относительно вращения и трансляции, поэтому для каждого сегмента я могу найти перевод, чтобы положить узел t = 0 в начало координат и поворот, чтобы положить t = 1 узел на оси + x:

Я могу найти уравнение для кривой, подключив точки и перегруппировав:

P(t) = (
    (t**3)*(-Dm1 + 3*D0 - 3*D1 + D2)
    + (t**2)*(3*Dm1 - 6*D0 + 3*D1)
    + t*(-3*Dm1 + 3*D1)
    + (Dm1 + 4*D0 + D1)
) / 6.

но я отрываю свои волосы, пытаясь их интегрировать - я могу сделать

 1
/
| Py(t) dt
/
t=0

но это не дает мне области. Я думаю, что мне нужно

 Px(t=1)
/
| Py(t) (dPx(t) / dt) dt
/
x = Px(t=0)

но прежде, чем я пойду дальше, мне очень хотелось бы знать:

Является ли это правильным вычислением для области? В идеале, аналитическое решение сделает мой день!
Как только я нахожу эту область, как я могу определить, нужно ли добавлять или вычитать ее из базового поля (красные или зеленые области на первой диаграмме)?
Существуют ли какие-либо модули Python, которые будут делать это для меня? У Numpy есть несколько методов для оценки кубических Bsplines, но нет, которые, похоже, имеют дело с областью.
Есть ли более простой способ сделать это? Я думаю о том, чтобы оценить P (t) в кучу точек - что-то вроде t = numpy.arange(0.0, 1.0, 0.05) - и рассматривать все это как многоугольник. Любая идея, сколько подразделений необходимо для обеспечения заданного уровня точности (я бы хотел, чтобы ошибка была меньше 1%)?

Ответ 1

Выберите произвольную точку как pivot p ₀ (например, начало координат (0,0))
Выберите некоторую точку вдоль кривой p ₁= (x, y)
Дифференцируем кривую в этой точке, чтобы получить скорость v = < vx, vy >
Создайте треугольник из трех точек и вычислить площадь. Проще всего сделать с помощью перекрестного произведения между векторами p ₀ p ₁ и v и делить на два.
Интегрируйте эту область по t, от 0 до 1.

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

Результат:

Область _i= (31 x _{i -1} y _i + 28 x _{i -1} y _{i +1} + x _{i - 1} y _{i +2} - 31 x _i y _{i -1} + 183 x _i y _{i +1} + 28 x _i y _{i +2} - 28 x _{i +1} y _{i -1} - 183 x _{i +1} y _i + 31 x _{i +1} y _{i +2} - x _{i +2} y _{i -1} - 28 x _{i +2} y _i - 31 x _{i +2} y _{i +1})/720

Если вы преобразуете его в матричную форму, вы получите:

Область _i= < x _{i -1} x _i x _{i +1} <Я > х <суб > + 2суб → </Мидот; P & middot; < y _{i -1} y _i y _{i +1} <Я > у <суб > + 2суб → ^T

где P -

[    0   31   28    1]
[  -31    0  183   28] / 720
[  -28 -183    0   31]
[   -1  -28  -31    0]

Если контрольные точки [(0,0) (1,0) (1,1) (0,1)], результирующие области становятся:

[(0,0), (1,0), (1,1), (0,1)] -> 242/720
[(1,0), (1,1), (0,1), (0,0)] -> 242/720
[(1,1), (0,1), (0,0), (1,0)] ->   2/720
[(0,1), (0,0), (1,0), (1,1)] ->   2/720

Сумма составляет 488/720 = 61/90.

Ответ 2

Лично я использовал сплайны в своих лучших интересах и переписал интеграл области как контурный интеграл, используя "Зеленая теорема" .

Поскольку вы уже знаете кривую, это будет непросто сделать интеграцию с использованием гауссовой квадратуры достаточного порядка. Никаких махинаций для оценки дополнительных областей. Я буду держать пари, что это будет также эффективным с точки зрения вычислений, потому что гауссова квадратура над многочленами так хорошо себя ведет. Кубические B-сплайны будут хорошо интегрироваться.

Я бы написал ваш код таким образом, чтобы первая и последняя точки должны совпадать. Это инвариант проблемы.

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

Ответ 3

Ну... похоже, вы уже знаете, что делать.

Да, это правильный способ вычисления области под сегментом, действительно
```
dS=Y*dX=Y*(dX/dt)*dt
```
Вам не нужно заботиться о добавлении или вычитании. Вы должны добавить все время, но некоторые из интегралов (красные сегменты) будут отрицательными (если вы всегда устанавливаете свой P _n в начале координат и P _{n + 1} вдоль положительного направления оси X). Поэтому для каждого сегмента вам нужно вычислить перевод и вращение, а затем интеграл (все сделано аналитически).
Я не знаю о модуле Python, но кажется, что для создания такого модуля потребуется максимум дня.
Я думаю, что аналитическое решение будет лучше и не так сложно.

В любом случае графика абсолютно потрясающая