нахождение площади замкнутого двумерного однородного кубического 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)

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

  1. Это правильный расчет площади? В идеале аналитическое решение сделало бы мой день лучше!

  2. Как только я найду эту область, как я могу определить, нужно ли мне добавить или вычесть ее из базового многоугольника (красные и зеленые области на первой диаграмме)?

  3. Существуют ли какие-либо модули Python, которые будут выполнять этот расчет за меня? У Numpy есть несколько методов для вычисления кубических B-сплайнов, но ни один из них, похоже, не имеет отношения к площади.

  4. Есть ли более простой способ сделать это? Я думаю о том, чтобы, возможно, оценить P (t) в нескольких точках - что-то вроде t = numpy.arange(0.0, 1.0, 0.05) - и рассматривать все это как многоугольник. Есть идеи, сколько подразделений необходимо, чтобы гарантировать заданный уровень точности (я бы хотел погрешность <1%)?


python numpy math spline
person Hugh Bothwell    schedule 03.06.2012    source источник

Ответы (3)


arrow_upward
3
arrow_downward

  1. Выберите произвольную точку как точку поворота p 0 (например, начало координат (0,0))
  2. Укажите точку на кривой p 1 = (x, y)
  3. Продифференцируйте кривую в этой точке, чтобы получить скорость v = ‹vx, vy>
  4. Сформируйте треугольник из трех точек и вычислить площадь. Легче всего сделать с помощью перекрестного произведения векторов p 0 p 1 и v , и разделив на два.
  5. Интегрируйте эту область по 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 < / i> +2 y i - 31 x i < / i> +2 y i +1) / 720

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

Площадь i = ‹x i -1 x i x i +1 x i +2 ›P‹ y i -1 y i y i +1 y i +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.

person Markus Jarderot    schedule 03.06.2012
comment
Почему P именно эта матрица? - person huon; 03.06.2012
comment
Это коэффициенты из решения. Я не анализировал это дальше. - person Markus Jarderot; 03.06.2012
comment
Есть более глубокое объяснение? Или это просто так? (Извините, пропустил вашу правку, не обращайте на это внимания :)) - person huon; 03.06.2012
comment
@Markus Jarderot: похоже, это легко реализовать, но я все еще не совсем понимаю, как вы получили P. Извините, если я немного медленный - мой последний урок исчисления был пятнадцать лет назад: - / - person Hugh Bothwell; 03.06.2012
comment
-1: Ваша матрица кажется неправильной. У меня есть 3/5 для площади вашего примера, вычислив интеграл вручную. Ваш результат хотя и близок. Кроме того, ваша процедура недостаточно ясна. Шаг 5 мог бы быть намного более подробным, и если бы это было так, читатели могли бы отредактировать ваш ответ с помощью правильной формулы. - person fkorsa; 02.02.2017

arrow_upward
4
arrow_downward

Лично я бы использовал сплайны как можно лучше и переписал бы интеграл по площади как контурный интеграл, используя теорему Грина .

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

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

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

person duffymo    schedule 03.06.2012
comment
Это очень хорошая идея ... Не знаю, как я это пропустил :) Никакой путаницы с многоугольниками, никакого сдвига / поворота .... Просто интегрируйте по кривой. Прохладный! - person BlacKow; 03.06.2012
comment
Обратите внимание, что это в основном то же решение, что и у Маркуса Джардерота, он просто объяснил, как работает теорема Грина, вместо того, чтобы называть ее. - person Vaughn Cato; 03.06.2012
comment
Нет, это не так. Теорема Грина работает совсем не так. Необходим линейный интеграл, а не произведение векторов треугольника. Вам лучше перечитать ту статью, которую я опубликовал. - person duffymo; 03.06.2012
comment
Последнее уравнение (интеграл для площади, заключенной в параметрическую кривую) можно переписать как ∫ (‹x, y› × ‹x ', y'›) / 2 dt. Это то же самое, что и шаги 4 и 5 в моем ответе. - person Markus Jarderot; 03.06.2012
comment
Но реализация, с помощью которой вы этого добьетесь, совершенно другая. Это означало бы квадратуру Гаусса с использованием точек вдоль линии, которые легко вычислить. Попробуйте оба варианта и скажите, что проще. Дополнительные очки за работу в области с дырой. - person duffymo; 03.06.2012
comment
@duffymo: спасибо, что указали мне на теорему Грина - я думаю, она сделает именно то, что мне нужно. Я обновлю свой вопрос, как только у меня заработает код решения. - person Hugh Bothwell; 03.06.2012
comment
Ключевым шагом за пределами формулы Вольфрама является запись контурного интеграла в терминах параметрического расстояния вдоль кривой. Двухмерный интеграл площади по (x, y) становится одномерным интегралом по параметрическому расстоянию t. Вы пишете x (t), x '(t), y (t) и y' (t) в терминах этого, используя уравнения сплайна. Квадратуру Гаусса сделать легко. - person duffymo; 03.06.2012

arrow_upward
1
arrow_downward

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

  1. Да, это правильный способ вычисления площади под сегментом, действительно

    dS=Y*dX=Y*(dX/dt)*dt

  2. Вам не нужно заботиться о сложении или вычитании. Вам нужно добавлять все время, но некоторые интегралы (красные сегменты) будут отрицательными (если вы всегда устанавливаете, ставьте P n в начало координат и P n + 1 < / sub> в положительном направлении оси X). Итак, для каждого сегмента вам нужно рассчитать перемещение и вращение, а затем интеграл (все выполняется аналитически).

  3. Не знаю, как насчет модуля Python, но похоже, что создание такого модуля займет не больше суток.
  4. Думаю, аналитическое решение будет лучше и не так сложно.

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

person BlacKow    schedule 03.06.2012
comment
Хороший ответ, и согласен по поводу графики! (Кроме того, я поставил n и n + 1 в качестве индексов: вы можете использовать буквальный HTML в ответах, например <sub>..</sub>) - person huon; 03.06.2012
comment
Спасибо за подсказку в нижнем индексе - person BlacKow; 03.06.2012