3D моделирование в Python

Допустим, вам потребовалось на языке программирования python, построить трёхмерную модель некоторого объекта, затем визуализировать его, или подготовить файл для печати на 3D принтере. Существует несколько библиотек, помогающих в решении этих задач. Поговорим о том, как строить трёхмерные модели из точек, граней и примитивов в python. Как выполнять элементарные приемы 3D моделирования: перемещение, поворот, объединение, вычитание и другие.

Рассмотрим библиотеки:

• numpy-stl
• pymesh
• pytorch3d
• SolidPython

Применяя каждую библиотеку, построим фрактал - Губку Менгера (Menger sponge), сохраним модель в stl файл и выполним рендеринг изображения. Попутно, кратко ознакомимся с используемыми структурами данных и терминами.

Примеры кода можно найти в репозитории на GitHub.

Подготовка среды

Numpy-stl: Обзор библиотеки

Страница проекта

Документация

Github

Строение полигональной модели:

Vertices (вершины) - список точек. Каждая точка описана тремя числами - координатами в 3-х мерном пространстве.

Если вы не знакомы с Jupyter notebook

Пример: numpy_stl_example_01.ipynb

import numpy as np
from myplot import plot_verticles

vertices = np.array(
    [
        [-3, -3, 0],
        [+3, -3, 0],
        [+3, +3, 0],
        [-3, +3, 0],
        [+0, +0, +3]
    ]
)

plot_verticles(vertices = vertices, isosurf = False)

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

plot_verticles(vertices = vertices, isosurf = True)

Выглядит так, будто грани уже существуют. Но пока у нас есть только вершины. Чтобы сформировать stl файл, опишем грани, что можно сделать вручную, или предоставить эту работу функции spatial.ConvexHull из библиотеки scipy

Пример: numpy_stl_example_02.ipynb

import numpy as np
from scipy import spatial
from stl import mesh
from myplot import plot_mesh

vertices = np.array(
    [
        [-3, -3, 0],
        [+3, -3, 0],
        [+3, +3, 0],
        [-3, +3, 0],
        [+0, +0, +3]
    ]
)

hull = spatial.ConvexHull(vertices)
faces = hull.simplices

В результате, массив faces содержит описание граней:

array([[4, 1, 0],
       [4, 2, 1],
       [3, 4, 0],
       [3, 4, 2],
       [3, 2, 1],
       [3, 1, 0]], dtype=int32)

faces (грани) - список граней. Каждая треугольная грань описана тремя вершинами, или точками. А точнее, позицией точек в массиве vertices.

Например, последняя грань содержит числа 3, 1, 0. Значит грань собрана из точек 0-го, 1-го и 3-го элементов массива vertices:

Mesh (сетка) - Набор вершин и граней, которые определяют форму многогранного объекта.

myramid_mesh = mesh.Mesh(np.zeros(faces.shape[0], dtype=mesh.Mesh.dtype))
for i, f in enumerate(faces):
    for j in range(3):
        myramid_mesh.vectors[i][j] = vertices[f[j],:]
plot_mesh(myramid_mesh)

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

myramid_mesh.save('numpy_stl_example_02.stl')

Для просмотра stl файлов разработано довольно много программ. Одна из них называется Blender, доступна для скачивания и не требует оплаты за использование

Метод spatial.ConvexHull предназначен для вычисления выпуклой оболочки и хорошо справляется с пирамидой и кубом. Но в объектах, имеющих впадины, часть точек будет потеряна и при сборке stl произойдет ошибка из-за несоответствия количества граней количеству точек.

Это хорошо видно на примере в двух измерениях: numpy_stl_example_03.ipynb

import matplotlib.pyplot as plt
from scipy import spatial
import numpy as np


points = np.array([
    [0,0],
    [-2,0],
    [-2,2],
    [0,1.5],
    [2,2],
    [2,0]
])

hull = spatial.ConvexHull(points)

В hull.simplices содержится описание граней:

array([[2, 1],
       [2, 4],
       [5, 1],
       [5, 4]], dtype=int32)

Отобразим вершины и грани на графике:

plt.plot(points[:,0], points[:,1], 'o')
for simplex in hull.simplices:
    plt.plot(points[simplex, 0], points[simplex, 1], 'k-')

Для таких случаев придется найти альтернативу ConvexHull, или описать грани вручную:

faces = np.array([
    [0, 1],
    [1, 2],
    [2, 3],
    [3, 4],
    [4, 5],
    [5, 0]
])

plt.plot(points[:,0], points[:,1], 'o')
for simplex in faces:
    plt.plot(points[simplex, 0], points[simplex, 1], 'k-')

Numpy-stl: Построение фрактала

Пришло время построить фрактал. В numpy-stl нет функции булева вычитания. Для построения фрактала Menger Sponge пойдем от обратного. В нашем распоряжении два метода:

1. Построение элементарного mesh куба. Назовем его voxel.
2. Объединение нескольких voxel в единый mesh.

Построим фрактал из кубиков, как из конструктора.

Описание логики построения фрактала

Пример: numpy_stl_example_04.ipynb

Показать код

Numpy-stl: Рендеринг изображения

Для рендеринга, в функцию вывода plot_mesh, будем передавать mesh, загруженный из stl файла.

Пример: numpy_stl_example_05.ipynb

Показать код

Показать анимацию

PyMesh: Обзор библиотеки

Документация

Github

Установка

Начнем с простого куба. В папке pymesh_examples находится скрипт pymesh_example_01.py. В дальнейшем, контейнер будет брать файлы именно из этой папки.

Пример: pymesh_example_01.py

import pymesh

box_a = pymesh.generate_box_mesh([0,0,0], [1,1,1])
filename = "/pymesh_examples/pymesh_example_01.stl"
pymesh.save_mesh(filename, box_a, ascii=False)

Из корня проекта запустим контейнер:

sh pymesh_example_01.sh

Что тут происходит?

Квадратное отверстие сделаем при помощи булева вычитания. Строим параллелепипед и вычитаем его из основного куба.

Пример: pymesh_example_02.py

import pymesh

box_a = pymesh.generate_box_mesh([0,0,0], [1,1,1])
box_b = pymesh.generate_box_mesh([0.4,0.4,0], [0.6,0.6,1])
box_c = pymesh.boolean(box_a, box_b, operation='difference', engine="igl")

filename = "/pymesh_examples/pymesh_example_02.stl"
pymesh.save_mesh(filename, box_c, ascii=False)

Запуск:

sh pymesh_example_02.sh

Функция boolean проста. Первым параметром передаётся объект из которого вычитаем. Вторым - который вычитаем. Следом передаются операция и движок. 

Здесь boolean применяется не только для вычитания. Всего доступно 4 операции:

• Intersection: A∩B
• Union: A∪B
• Difference: A∖B (два последних примера)
• Symmetric difference: A XOR B (изображение не представлено)

Чтобы лучше понять как перемещать и вращать объект, бывает удобно временно заменить операцию Difference на Union.

Сделаем второе отверстие, переместим и повернем его.

Пример: pymesh_example_03.py

Показать код

Запуск:

sh pymesh_example_03.sh

В этом скрипте добавлены функции перемещения и поворота. При перемещении создается новый mesh объект на основании измененных вершин и граней исходного объекта. Для поворота, сначала при помощи класса Quaternion, описывается поворот, а затем, аналогично случаю с перемещением, используется создание нового, повернутого объекта, на основании вершин и граней исходного объекта, а также описания поворота.

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

О кватернионах есть довольно подробная статья:

Доступно о кватернионах и их преимуществах

В результате исполнения скрипта, получается куб с двумя пересекающимися отверстиями:

Перечисленных инструментов достаточно для построения фрактала.

PyMesh: Построение фрактала

Пример: pymesh_example_04.py

Показать код

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

На этот раз при запуске необходимо явно указать глубину фрактала:

sh pymesh_example_04.sh 3

Исполняться он будет 5-15 минут. После исполнения, в папке pymesh_examples появляется stl файл:

Если запросить фрактал 4-го уровня

PyMesh: Рендеринг изображения

Mesh мы уже поворачивали, на этот раз повернём камеру.

Пример: pymesh_example_05.py

Показать код

Запуск:

sh pymesh_example_05.sh

Показать анимацию

PyTorch3d: Обзор библиотеки

Страница проекта

Документация

Установка

Github

Построение пирамиды. Пример: pytorch3d_example_01.py

Показать код

Подход очень похож на аналогичный в numpy-stl. Но так, как предполагается работа на GPU, будем разделять понятия хоста и устройства. Хост - наш компьютер. Устройство - GPU карта. Если карты нет, пользоваться библиотекой тоже можно, тогда в роли GPU будет выступать CPU. У хоста и устройства своя память. Чтобы передать объект с хоста на устройство и обратно, необходимо произвести небольшой ритуал. 

В примере ниже, сразу на устройстве описываем вершины, копируем их с устройства на хост. На основании вершин вычисляем грани. Сохраняем объект. Файл в формате obj можно импортировать в blender:

Обратите внимание на команду verts.cpu().numpy()

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

Пример: pytorch3d_example_02.py

Показать код

PyTorch3d: Построение фрактала

Использование GPU даёт некоторый прирост производительности.

Пример: pytorch3d_example_03.py

Показать код

Cкорость расчета увеличилась на порядок, что позволило примерно за 5 часов построить фрактал 5 уровня:

Размер stl файла составил 1.9 ГБ. При построении фрактала 5-го уровня, программа останавливалась из-за переполнения памяти видеокарты. Пришлось сборку объекта выполнять пакетами. Создавалось по 10 слоев “двумерных” фракталов, затем они присоединялись к основному объекту, до тех пор, пока не построился полный фрактал.

PyTorch3d: Рендеринг изображения

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

Пример: pytorch3d_example_04.py

Показать код

Показать анимацию

SolidPython: Обзор библиотеки

Документация SolidPython

Github

Openscad cheatsheet

Wiki

SolidPython Самая богатая методами моделирования библиотека, среди перечисленных. 3D сцена описывается на python, в формате, очень похожем на openscad, генерируется openscad код, который пишется в scad файл и далее его можно редактировать в openscad или сразу сохранить в stl.

Пример: solidpython_example_01.ipynb

from solid import *

d = difference()(
    cube(size = 10, center = True),
    sphere(r = 6.5, segments=300)
)
path = scad_render_to_file(d, 'solidpython_example_01.scad')

Для указания разрешения сферы, вместо привычного для openscad параметра $fn, используем слово segments.

solidpython удобно отлаживать. С одной стороны экрана открыт scad файл, с другой jupyter notebook. При исполнении scad_render_to_file картинка в openscad автоматически обновляется.

Если нужен stl, openscad умеет рендерить файлы этого формата через команды консоли. Пример вызова из jupyter notebook:

!openscad solidpython_example_01.scad -o solidpython_example_01.stl

Общий принцип такой: Любая функция возвращает объект. Если над объектом необходимо произвести некоторое действие, объект (или список объектов) передается в круглых скобках после вызова соответствующей функции.

Пример: solidpython_example_02.ipynb

from solid import *

c = circle(r = 1)
t = translate([2, 0, 0]) (c)
e = linear_extrude(height = 10, center = True, convexity = 10, twist = -500, slices = 500) (t)
col = color('lightgreen') (e)
path = scad_render_to_file(col, 'solidpython_example_02.scad')

Тут разрешение регулируется параметром slices.

SolidPython: Построение фрактала

Пример: solidpython_example_03.ipynb

Показать код

SolidPython: Рендеринг изображения

Сформируем серию изображений последней сцены, поворачивая камеру в каждом изображении.

Пример: solidpython_example_04.ipynb

Показать код

Кроме того, solidpython предлагает формирование анимации средствами openscad. В документации об этом есть небольшой раздел с примером.

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

Пример: solidpython_example_05.ipynb

КОД

import solid as scad
import subprocess
import os

# sponge
depth = 3
sponge = scad.cube([1,1,1], center = True)
z = 1.1
side = 1
for d in range(1, depth+1):
    side /= 3
    for x in range(1,3**d,3):
        for y in range(1,3**d,3):
            box_a = scad.cube([side,side,z], center = True)
            box_a = scad.translate([side*x-0.5+side/2, side*y-0.5+side/2, 0]) (box_a)            
            box_b = scad.rotate([90,0,0]) (box_a)
            box_c = scad.rotate([0,90,0]) (box_a)
            sponge -= box_a
            sponge -= box_b
            sponge -= box_c

sponge = scad.translate([3.3, 0.5, 0.5]) (sponge)

# text
text = scad.text('Pyth   n 3d', size = 1)
text_gray = scad.linear_extrude(height = 0.9) (text)
text_gray = scad.color('gray', 0.5) (text_gray)

text_white = scad.linear_extrude(height = 0.1) (text)
text_white = scad.color('white', 0.8) (text_white)
text_white = scad.translate([0, 0, 0.9]) (text_white)

text_blue = scad.linear_extrude(height = 0.1) (text)
text_blue = scad.color('blue', 0.2) (text_blue)
text_blue = scad.translate([0, 0, 0.9]) (text_blue)

# merge
objects = text_gray + text_white + text_blue
for i in range(0,10):
    sponge_entity = scad.color('magenta', 0.5/(i**2+1)) (sponge)
    sponge_entity = sponge_entity + scad.color('blue', 0.5/(i**2+1)) (sponge)
    objects += scad.translate([0, 0, i*2]) (sponge_entity)
    if i>0:
        objects += scad.translate([0, 0, -i*2]) (sponge_entity)
        objects += scad.translate([0, i*2, 0]) (sponge_entity)
        objects += scad.translate([0, -i*2, 0]) (sponge_entity)
        
# render
filename = 'solidpython_example_05'
render = scad.scad_render_to_file(objects, filename+'.scad')
# camera parameters when exporting png: 
#                                 =translate_x,y,z,rot_x,y,z,dist or 
#                                 =eye_x,y,z,center_x,y,z

#  =colorscheme: *Cornfield | Metallic | Sunset | 
#                                 Starnight | BeforeDawn | Nature | DeepOcean | 
#                                 Solarized | Tomorrow | Tomorrow 2 | Tomorrow 
#                                 Night | Monotone
cmd = "openscad "+filename+".scad --camera '3,2,2,20,5,-70,30' --imgsize=10000,3000 --colorscheme='Tomorrow' -o "+filename+".png"
MyOut = subprocess.Popen(
    cmd, 
    shell=True, 
    stdout=subprocess.PIPE, 
    stderr=subprocess.PIPE, 
    stdin=subprocess.PIPE
)
MyOut.stdout, MyOut.stderr

Сравнение библиотек

Сравнение производительности не совсем объективно, так как имеются значительные различия в алгоритмах. В Pymesh и SolidPython применялось вычитание, тогда как в Numpy-stl и Pytorch3d объединение mesh.