К вычислению базисных сплайнов
"The mathematics is basically right or wrong; software is a large gray area with no boundaries." (L. Piegl and W. Tiller. The NURBS Book. p. 594)
Базовые алгоритмы библиотеки OpenCascade разрабатывались довольно давно. С одной стороны, классическая математика не устаревает, и надежные вычислительные схемы ядра работают по-прежнему хорошо. С другой стороны, иногда очень непросто разобраться в том, какова идея конкретно взятого алгоритма. Мы уже как-то писали, что хороший алгоритм должен комплектоваться грамотной документацией, чтобы однажды не стать реликтом — эдаким мертвым символом прежних эпох, — который легче переписать, чем исправить. Задача переписывания ядра OpenCascade на сегодняшний день не стоит, да и ни к чему здесь радикальные методы, коль скоро есть разумная стратегия замещения.
Сделаем важное отступление, прежде чем двинуться дальше. Для разработчиков инженерного ПО ценность представляет не геометрическое ядро как таковое, а наличие достаточного набора операторов моделирования. Условно, мы разделяем все операторы на ГЕОМЕТРИЧЕСКИЕ (CAGD) и ТОПОЛОГИЧЕСКИЕ (CAD).
CAGD = Computer-Aided Geometric Design
CAD = Computer-Aided Design
Эти аббревиатуры, кстати, соответствуют двум «главным» научным журналам, посвященным инженерному геометрическому моделированию. «The word geometric distinguishes this field from computer aided logical design» [R.E. Barnhill, R.F. Riesenfeld, Computer Aided Geometric Design, 1974].
Какие-то операторы доступны в библиотеке OpenCascade, какие-то — нет. Любопытно в этой связи отметить пару фактов. Во-первых, не секрет, что компания OPEN CASCADE занимается разработкой коммерческих алгоритмов, которые обогащают бесплатный набор операторов. То есть расширение алгоритмической палитры — часть коммерческой стратегии компании. Во-вторых, в свое время компания Matra Datavision поглотила компанию CISIGRAPH как раз для того, чтобы расширить доступный набор операторов CAGD. Например, именно так в ядре появился алгоритм построения поверхностей класса А, остающийся востребованным по сей день (пакет GeomPlate).
Основой геометрических операторов ядра OpenCascade является пакет вычисления базисных сплайнов и смежные фундаментальные алгоритмы. Эта функциональность реализована в пакетах BSplCLib и BSplSLib, датированных 1991 годом. С тех пор исследования в области CAGD не топтались на месте. В частности, была издана Главная Книга CAGD — Книга NURBS [L. Piegl, W. Tiller, The NURBS Book, Springer, New York, 1997].
Сегодня, если вы желаете реализовать собственные операторы моделирования, у вас есть выбор: базироваться на функциональности OpenCascade, либо обойтись без нее. В частности, возникает вопрос, следует ли использовать пакеты BSplCLib/BSplSLib? Не устарели ли они? Давайте попробуем в этом разобраться.
Вычисление кривой
Поставим вычислительный эксперимент. Возьмем B-кривую ядра OpenCascade и B-кривую, созданную независимо, то есть без использования сторонних библиотек. Для независимой кривой условимся, что вычисление базисных сплайнов осуществляет алгоритм A2.2 Книги NURBS (он идет в связке с алгоритмом A2.1). Этот алгоритм расчитывает значения по рекуррентной формуле Кокса-де Бура и считается в некотором смысле «стандартным». Интересно посмотреть, насколько быстро работает вычисление точки кривой в обоих случаях.
Реализация алгоритмов A2.1 и A2.2 из Книги NURBS не представляет трудности, поскольку исходный код можно, грубо говоря, скопировать из текста монографии, и он будет работать. Сделаем лишь пару оговорок относительно тех граблей, на которые легко наступить при недостаточно аккуратном программировании. Во-первых, избегайте динамического выделения памяти. Во-вторых, следите за тем, чтобы в вашем коде не происходило излишнее копирование структур данных, таких, например, как узловой вектор. В противном случае вся заявленная теоретическая эффективность алгоритмов A2.1 и A2.2 останется на бумаге.
Таблица ниже дает замеры времени для алгоритма вычисления точки на кривой. Тестировалось на железе Intel Core i5-3470 CPU @ 3.20GHz. Каждый замер отвечает одному миллиону повторов.
| Степень сплайна | Количество узлов | OpenCascade [сек] | Алгоритм A3.1 из Книги NURBS [сек] |
| 3 | 73 | 0.12 | 0.11 |
| 3 | 24 | 0.11 | 0.10 |
| 8 | 1492 | 0.29 | 0.28 |
| 14 | 93 | 0.6 | 0.62 |
Таким образом, при аккуратной реализации вы получите производительность сопоставимую с OpenCascade. Отсюда вывод: вычисление B-сплайнов в ядре OpenCascade сделано как минимум неплохо, если сравнивать его эффективность со «стандартным» методом из Книги NURBS. О том же, кстати, свидетельствует упомянутая выше статья С. Гатилова.
Разработчики геометрического ядра NLib, являющегося воплощением теории из Книги NURBS, решили проблему с динамическим выделением памяти на архитектурном уровне. Алгоритмы оперируют т.н. «стэками» указателей на области памяти, которые могут использоваться локально, либо передаваться между несколькими процедурами.
Экстраполяция
Алгоритмы BSplCLib позволяют вычислять значения сплайнов ВНЕ области их определения. Достигается это тривиальным способом: в случае выпадения параметра из области определения, алгоритм BSplCLib::LocateParameter() — аналог A2.1 из Книги NURBS — возвращает индекс ближайшего интервала узлового вектора. Поскольку все базисные сплайны внутри интервалов суть полиномы, их можно вычислять для любых значений параметров.
Конечно, такой подход к экстраполяции едва ли хорош. Тем не менее, он часто используется на практике. Этот подход называют «натуральной» экстраполяцией, поскольку он не требует привлечения дополнительных условий для продления кривых и поверхностей. Подробнее о методах экстраполяции можно почитать в работе [Wolters H.J. Extensions: extrapolation methods for CAD. 1999], а также, например, в [Hu S.-M., Tai C.-L., Zhang S.-H. An extension algorithm for B-splines by curve unclamping // Computer-Aided Design. 2002. N 5 (34). C. 415–419].
Вычисление поверхности
Таблица ниже дает замеры времени для алгоритма вычисления точки на поверхности (среднее за миллион повторов).
| Степень сплайнов (u, v) | Количество контрольных точек | OpenCascade [сек] | Алгоритм A3.5 из Книги NURBS [сек] |
| 1, 1 | 9 | 0.13 | 0.13 |
| 3, 3 | 1752 | 0.28 | 0.32 |
| 4, 4 | 3973 | 0.38 | 0.4 |
| 8, 3 | 108 | 0.51 | 0.63 |
Таким образом, мы снова наблюдаем хорошую производительность ядра OpenCascade.
Вычисление производных
Если вычислять не только значения базисных сплайнов, но и их производные, то скорость работы «стандартного» алгоритма (A2.3 Книги NURBS) окажется хуже OpenCascade, поскольку он, вообще говоря, требует динамического выделения памяти. Арифметические вычисления теряются на фоне чисто утилитарных операций, таких как зануление элементов массива, аллокация и деаллокация.
Для повышения производительности здесь не обойтись без использования специальных техник управления памятью.