Невероятно быстрое вычисление объёма с помощью теоремы о дивергенции

Представлен быстрый алгоритм вычисления объёма простой, замкнутой, триангулированной 3D-сетки. Это предположение — прямое следствие теоремы о дивергенции. Дальнейшие обобщения на другие типы сеток возможны, но пока выходят за рамки рассмотрения.

Отправная точка — определение объёма как тройного интеграла по области от константы, равной единице:

$V={\iiint }_{R}1dV$

Пусть $𝐅$ — функция в ${ℝ}^3$ такая, что её дивергенция равна единице. Для целей данного вывода выбирается:

Легко проверить, что

$div𝐅=\frac{\partial F}{\partial x}+\frac{\partial F}{\partial y}+\frac{\partial F}{\partial z}=1+0+0=1$

Следовательно,

$V={\iiint }_{R}1dV={\iiint }_{R}div𝐅(x,y,z)dV$

По теореме о дивергенции это равно поверхностному интегралу:

$V={\iint }_S𝐅(x,y,z)d𝐒$

Этот поверхностный интеграл, определённый по поверхности S 3D-сетки, равен сумме своих кусочных треугольных частей. Пусть $T_i$ обозначает поверхность $i$-го треугольника сетки. Тогда,

$V=\sum _{i=0}{\iint }_{T_i}𝐅(x,y,z)d𝐒$

Пусть $T_{in}$ представляет $n$-ю вершину $i$-го треугольника. Пусть ${\Delta }_1$ равно векторной разности между $T_{i1}$ и $T_{i0}$, а ${\Delta }_2$ аналогично равно $T_{i2}-Ti0$. Каждый отдельный треугольник $T_i$ таким образом может быть параметризован как:

$𝐫(u,v)=T_{i0}+u{\Delta }_1+v{\Delta }_2$

Тогда простое дифференцирование даёт:

${𝐫}_u={\Delta }_1$ ${𝐫}_v={\Delta }_2$

Следовательно,

${𝐫}_u\times {𝐫}_v={\Delta }_1\times {\Delta }_2$

Таким образом, поверхностный интеграл можно переписать в терминах этой параметризации, подставив определение $𝐅$ где необходимо:

$V=\sum _{i=0}{\iint }_{T_i}𝐅(x,y,z)({𝐫}_u\times {𝐫}_v)dA$ $=\sum _{i=0}{\iint }_{T_i}𝐅(x,y,z)\dot{(}{\Delta }_{i1}\times {\Delta }_{i2})dA$

Это векторное произведение постоянно на всём треугольнике и легко вычисляется по данным вершин. Достаточно вычислить только X-компоненту векторного произведения; остальные равны нулю из-за скалярного произведения с нулевыми компонентами $𝐅$. $V$ таким образом можно переписать как:

$V=\sum _{i=0}{({\Delta }_{i1}\times {\Delta }_{i2})}_x{\iint }_{T_i}xdA$

Далее рассматривается поверхностный интеграл ${\iint }_{T_i}xdA$. Раскрытие с помощью параметризации даёт:

${\iint }_{T_i}xdA={\int }_0^1{\int }_0^{u}xdvdu={\int }_0^1{\int }_0^u(T_{i0x}+u{\Delta }_{i1x}+v{\Delta }_{i2x})dvdu$

Этот интеграл можно вычислить напрямую, рассматривая данные вершин как константы:

${\int }_0^1{\int }_0^{1-u}(T_{i0x}+u{\Delta }_{i1x}+v{\Delta }_{i2x})dvdu$ $=T_{i0x}{\int }_0^1{\int }_0^{1-u}dvdu+{\Delta }_{i1x}{\int }_0^1{\int }_0^{1-u}udvdu+{\Delta }_{i2x}){\int }_0^1{\int }_0^{1-u}vdvdu$ $=T_{i0x}\left(\frac{1}{2}\right)+{\Delta }_{i1x}\left(\frac{1}{6}\right)+{\Delta }_{i2x}\left(\frac{1}{6}\right)$ $=T_{i0x}\left(\frac{1}{2}\right)+(T_{i1x}-T_{i0x})\left(\frac{1}{6}\right)+(T_{i2x}-T_{i0x})\left(\frac{1}{6}\right)$ $=T_{i0x}\left(\frac{1}{6}\right)+\left(T_{i1x}\right)\left(\frac{1}{6}\right)+\left(T_{i2x}\right)\left(\frac{1}{6}\right)$ $=\frac{1}{6}(T_{i0x}+T_{i1x}+T_{i2x})$

Подставляя это в исходную сумму и вынося постоянный множитель $\frac{1}{6}$ за пределы внутреннего цикла (чтобы избежать деления в цикле), получается следующая компактная формула для объёма:

$V=\frac{1}{6}\sum _{i=0}{({\Delta }_{i1}\times {\Delta }_{i2})}_x(T_{i0x}+T_{i1x}+T_{i2x})$

Анализ производительности

Итоговый алгоритм не содержит ни численного интегрирования, ни дифференцирования. В отличие от распространённых наивных алгоритмов вычисления объёма, которые по сути эквивалентны рендерингу сетки с последующей выборкой из результата рендера — дорогостоящей операции, — здесь присутствует только один цикл, проходящий по треугольникам. Таким образом, алгоритм вычисления объёма имеет сложность O(n) относительно числа треугольников. Более того, вычисление для каждого отдельного треугольника столь же эффективно: при естественном раскрытии векторного произведения внутренняя часть содержит семь сложений и три умножения. Вне цикла требуется всего одно умножение. Таким образом, для сетки из $n$ треугольников алгоритму требуется $8n-1$ сложений и $3n+1$ умножений, то есть $11n$ операций с плавающей точкой. Это очень быстро.

Для примерной оценки: если объём нужно вычислять каждый кадр в высокопроизводительном приложении на 60 кадров в секунду, без помощи GPU, используя только возможности CPU Raspberry Pi за $35, за один кадр можно обработать порядка 30 миллионов треугольников.

Мотивация

Скоро экзамен по векторному исчислению, и требуется к нему подготовиться. К тому же, кто не любит 3D-графику?!

Было бы (приятной) неожиданностью, если бы этот алгоритм оказался новым. Дополнительное исследование, проведённое после публикации, показало, что статья Efficient Feature Extraction for 2D/3D Objects in Mesh Representation авторства Ча Чжэн и Цухана Чена, судя по всему, описывает тот же алгоритм, хотя и с другим выводом. Было весело, пока это длилось!