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

Kate

Administrator
Команда форума
Допустим, вам потребовалось на языке программирования python, построить трёхмерную модель некоторого объекта, затем визуализировать его, или подготовить файл для печати на 3D принтере. Существует несколько библиотек, помогающих в решении этих задач. Поговорим о том, как строить трёхмерные модели из точек, граней и примитивов в python. Как выполнять элементарные приемы 3D моделирования: перемещение, поворот, объединение, вычитание и другие.
Рассмотрим библиотеки:
  • numpy-stl
  • pymesh
  • pytorch3d
  • SolidPython
Применяя каждую библиотеку, построим фрактал - Губку Менгера (Menger sponge), сохраним модель в stl файл и выполним рендеринг изображения. Попутно, кратко ознакомимся с используемыми структурами данных и терминами.
Примеры кода можно найти в репозитории на GitHub.
Подготовка среды
  • Docker для исполнения примеров Pymesh:
curl -fsSL https://get.docker.com -o get-docker.sh
sudo sh get-docker.sh
  • Компилятор g++ для PyTorch3d:
sudo apt update
sudo apt install g++
  • Клонирование репозитория и установка библиотек:
git clone https://github.com/format37/python3d.git
cd python3d
pip install -r requirements.txt
pip install "git+https://github.com/facebookresearch/pytorch3d.git"
  • Установка OpenScad:
sudo apt-get install openscad

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

Страница проекта
Документация
Github
Строение полигональной модели:
Изображение из Википедии
Изображение из Википедии
Vertices (вершины) - список точек. Каждая точка описана тремя числами - координатами в 3-х мерном пространстве.
Если вы не знакомы с Jupyter notebook
Что такое jupyter-ноутбук и зачем он нужен
Пример: 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[j] = vertices[f[j],:]
plot_mesh(myramid_mesh)
Грани пирамиды


Грани пирамиды
Как видно из рисунка, одна грань пирамиды оказалась перевернута. В последующих примерах, при построении фракталов, метод ConvexHull применяться не будет, так как он располагает точки грани в произвольном порядке, что приводит к переворачиванию некоторых граней.
myramid_mesh.save('numpy_stl_example_02.stl')
Для просмотра stl файлов разработано довольно много программ. Одна из них называется Blender, доступна для скачивания и не требует оплаты за использование
Снимок экрана. Blender: numpy_stl_example_02.stl


Снимок экрана. Blender: numpy_stl_example_02.stl
Метод 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.
Построим фрактал из кубиков, как из конструктора.
Описание логики построения фрактала
Договоимся что длина стороны грани фрактала - единица. Глубина фрактала это, грубо говоря, количество уникальных размеров отверстий. Длина вокселя зависит от глубины фрактала, она делится на 3 с каждым новым уровнем глубины.
Найдем сторону вокселя на глубинах 1 и 2. Упростим задачу, свернув фрактал с 3-х до 1-мерного случая:
xalbrr-2y0ukyaaaqpic-lchdsq.png

Если фрактал второго уровня, длина стороны куба составит 1/(3**2) или 1/9. Составим набор кубов так, чтобы своим расположением они заполнили исходный куб. Получится воксельный куб. Вычислим области отверстий. Исключим воксели, которые входят в области отверстий. В завершение, объединим оставшиеся воксели в один объект и сохраним.
Пример: numpy_stl_example_04.ipynb
Показать код
from stl import mesh
import numpy as np
import datetime

time_start = datetime.datetime.now()


def create_voxel(side, tx, ty, tz):
vertices = np.array([
[tx + 0, ty + 0, tz + 0],
[tx + side, ty + 0, tz + 0],
[tx + side, ty + side, tz + 0],
[tx + 0, ty + side, tz + 0],
[tx + 0, ty + 0, tz + side],
[tx + side, ty + 0, tz + side],
[tx + side, ty + side, tz + side],
[tx + 0, ty + side, tz + side]
])

faces = np.array([
[3, 2, 1],
[3, 1, 0],
[0, 1, 5],#
[5, 4, 0],
[7, 3, 0],
[0, 4, 7],#
[1, 2, 6],#
[6, 5, 1],
[2, 3, 6],#
[3, 7, 6],#
[4, 5, 6],#
[6, 7, 4]
])
voxel = mesh.Mesh(np.zeros(faces.shape[0], dtype=mesh.Mesh.dtype))
for i, f in enumerate(faces):
for j in range(3):
voxel.vectors[j] = vertices[f[j],:]
return voxel


def voxel_in_hole(holes, x, y):
comp_error = 0
for hole in holes:
if x + comp_error >= hole['x'][0] and x - comp_error < hole['x'][1] and \
y + comp_error >= hole['y'][0] and y - comp_error < hole['y'][1]:
return True
return False


depth = 3
side = 1/(3**depth)

# holes
holes = []
side = 1
round_factor = 14
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):
holes.append({
'x':[round(side*x,round_factor), round(side*x+side,round_factor)],
'y':[round(side*y,round_factor), round(side*y+side,round_factor)]
})

sponge = None
for x in range(3**depth):
for y in range(3**depth):
for z in range(3**depth):
rx = round(x/(3**depth),round_factor)
ry = round(y/(3**depth),round_factor)
rz = round(z/(3**depth),round_factor)
if voxel_in_hole(holes, ry, rz) or \
voxel_in_hole(holes, rx, ry) or \
voxel_in_hole(holes, rx, rz):
continue
if sponge is None:
sponge = create_voxel(side, x/(3**depth), y/(3**depth), z/(3**depth))
else:
new_voxel = create_voxel(side, x/(3**depth), y/(3**depth), z/(3**depth))
sponge = mesh.Mesh(np.concatenate([sponge.data, new_voxel.data]))

sponge.save('numpy_stl_example_04.stl')


time_end = datetime.datetime.now()
print('spent', (time_end - time_start).seconds, 'seconds')

Снимок экрана. Blender: numpy_stl_example_04.stl


Снимок экрана. Blender: numpy_stl_example_04.stl

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

Для рендеринга, в функцию вывода plot_mesh, будем передавать mesh, загруженный из stl файла.
Пример: numpy_stl_example_05.ipynb
Показать код
from myplot import plot_mesh
from stl import mesh
import datetime
import imageio
import math
import os

sponge = mesh.Mesh.from_file('numpy_stl_example_04.stl')

gif_half_size = 14
filenames = []
imageio_images = []
for i in range(gif_half_size):
sponge.rotate([1, 0.0, 1], math.radians(1))
filename = 'solidpython_example_03_'+str(i)+'.png'
filenames.append(filename)
plot_mesh(sponge, filename = filename)
print(datetime.datetime.now(), filename)

print('forward rendering..')
for i in range(gif_half_size):
print(datetime.datetime.now(), filenames)
imageio_images.append(imageio.imread(filenames))
i-=1

print('reverse rendering..')
i = gif_half_size-1
while i>0:
print(datetime.datetime.now(), filenames)
imageio_images.append(imageio.imread(filenames))
i-=1

# save gif
imageio.mimsave('solidpython_example_03.gif', imageio_images)

# remove images
for filename in filenames:
os.remove(filename)

print('job complete!')
Один из кадров анимации


Один из кадров анимации. Показать анимацию
Анимация: solidpython_example_03.gif


Анимация: solidpython_example_03.gif

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

Документация
Github
Установка
К сожалению, библиотека PyMesh не установилась у меня ни через менеджер пакетов PIP (несмотря на упоминание этого способа в документации), ни через Anaconda. Есть два способа установки.
  1. Следуя инструкции, собрать из исходников.
  2. Используя Docker контейнер. Выберем этот вариант, как более интересный. Запуск контейнера будет инициирован с параметрами. При помощи параметров запуска смонтируем в контейнер папку скриптов. Передадим необходимые параметры в скрипт. После завершения работы скрипта, контейнер будет удален. Если Docker вы уже установили следуя инструкции в начале статьи, более ничего устанавливать не нужно.
Если Docker вам не подходит, скомпилируйте PyMesh, следуя инструкции в документации. Такой вариант тоже будет работать.
Начнем с простого куба. В папке 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. При первом запуске будет загружен образ, это займёт некоторое время.
  • Папка pymesh_examples монтируется в контейнер
  • В контейнере запускается python скрипт /pymesh_examples/pymesh_example_01.py
  • Импортируется библиотека pymesh
  • функция generate_box_mesh генерирует куб на основании двух противоположных вершин в точках [0,0,0] и [1,1,1]
  • функция save_mesh сохраняет объект в stl файл.
  • После исполнения, в папке pymesh_examples появляется файл pymesh_example_01.stl
Снимок экрана. Blender: pymesh_example_01.stl


Снимок экрана. Blender: pymesh_example_01.stl
Квадратное отверстие сделаем при помощи булева вычитания. Строим параллелепипед и вычитаем его из основного куба.
Пример: 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
Снимок экрана. Blender: pymesh_example_02.stl


Снимок экрана. Blender: pymesh_example_02.stl
Функция boolean проста. Первым параметром передаётся объект из которого вычитаем. Вторым - который вычитаем. Следом передаются операция и движок.
Здесь boolean применяется не только для вычитания. Всего доступно 4 операции:
  • Intersection: AB
  • Union: AB
  • Difference: AB (два последних примера)
  • Symmetric difference: A XOR B (изображение не представлено)
Иллюстрация применения операции boolean на сфере и кубе
Иллюстрация применения операции boolean на сфере и кубе
Чтобы лучше понять как перемещать и вращать объект, бывает удобно временно заменить операцию Difference на Union.
Сделаем второе отверстие, переместим и повернем его.
Пример: pymesh_example_03.py
Показать код
import pymesh
import numpy as np
import math

def mesh_trans(mesh, x, y, z):
return pymesh.form_mesh(mesh.vertices + [[x, y, z]], mesh.faces)

def mesh_rotate(mesh, x, y, z, angle):
rot = pymesh.Quaternion.fromAxisAngle(np.array([x, y, z]), math.pi*2*angle/360)
return pymesh.form_mesh(np.dot(rot.to_matrix(), mesh.vertices.T).T, mesh.faces)

sponge = pymesh.generate_box_mesh([0,0,0], [1,1,1])
bool_box = pymesh.generate_box_mesh([0.4,0.4,-0.1], [0.6,0.6,1.1])
sponge = pymesh.boolean(sponge, bool_box, operation='difference', engine="igl")

bool_box = mesh_rotate(mesh = bool_box, x = 1, y = 0, z = 0, angle = 90)
bool_box = mesh_trans(mesh = bool_box, x = 0, y = 1, z = 0)
sponge = pymesh.boolean(sponge, bool_box, operation='difference', engine="igl")

filename = "/pymesh_examples/pymesh_example_03.stl"
pymesh.save_mesh(filename, sponge, ascii=False)
Запуск:
sh pymesh_example_03.sh
В этом скрипте добавлены функции перемещения и поворота. При перемещении создается новый mesh объект на основании измененных вершин и граней исходного объекта. Для поворота, сначала при помощи класса Quaternion, описывается поворот, а затем, аналогично случаю с перемещением, используется создание нового, повернутого объекта, на основании вершин и граней исходного объекта, а также описания поворота.
Кватернионы - система гиперкомплексных чисел, образующая векторное пространство размерностью четыре над полем вещественных чисел.
О кватернионах есть довольно подробная статья:
Доступно о кватернионах и их преимуществах
В результате исполнения скрипта, получается куб с двумя пересекающимися отверстиями:
Снимок экрана. Blender: pymesh_example_03.stl


Снимок экрана. Blender: pymesh_example_03.stl
Перечисленных инструментов достаточно для построения фрактала.

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

Пример: pymesh_example_04.py
Показать код
import numpy as np
import datetime
import pymesh
import math
import sys

time_start = datetime.datetime.now()

depth = int(sys.argv[1])

def mesh_trans(mesh, x, y, z):
return pymesh.form_mesh(mesh.vertices + [[x, y, z]], mesh.faces)

def mesh_rotate(mesh, x, y, z, angle):
rot = pymesh.Quaternion.fromAxisAngle(np.array([x, y, z]), math.pi*2*angle/360)
return pymesh.form_mesh(np.dot(rot.to_matrix(), mesh.vertices.T).T, mesh.faces)

def menger_sponge(depth):
mesh_fractal = pymesh.generate_box_mesh([0,0,0], [1,1,1])
z = [-0.1,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):
print(datetime.datetime.now(), d, x, y, '/', depth, ':', 3**d)
box_a = pymesh.generate_box_mesh([side*x,side*y,z[0]], [side*x+side,side*y+side,z[1]])

box_b = mesh_rotate(mesh = box_a, x = 1, y = 0, z = 0, angle = 90)
box_b = mesh_trans(mesh = box_b, x = 0, y = 1, z = 0)

box_c = mesh_rotate(mesh = box_a, x = 0, y = 1, z = 0, angle = 90)
box_c = mesh_trans(mesh = box_c, x = 0, y = 0, z = 1)

mesh_fractal = pymesh.boolean(mesh_fractal, box_a, operation='difference', engine="igl")
mesh_fractal = pymesh.boolean(mesh_fractal, box_b, operation='difference', engine="igl")
mesh_fractal = pymesh.boolean(mesh_fractal, box_c, operation='difference', engine="igl")

return mesh_fractal

mesh_fractal = menger_sponge(depth)

pymesh.save_mesh("/pymesh_examples/pymesh_example_04_"+str(depth)+".stl", mesh_fractal, ascii=False)

time_end = datetime.datetime.now()
print('spent', (time_end - time_start).seconds, 'seconds')
В этом скрипте добавлен входящий параметр для передачи глубины вычисления фрактала. Для каждой глубины создается параллелепипед, который затем дважды копируется, поворачивается и смещается. Получается всего 3 параллелепипеда, которые вычитаются из основного куба. По одному на каждую грань. Операция повторяется x и y раз, чтобы заполнить все строки и колонки грани. Проверка на вычитание из пустого пространства не выполняется.
На этот раз при запуске необходимо явно указать глубину фрактала:
sh pymesh_example_04.sh 3
Исполняться он будет 5-15 минут. После исполнения, в папке pymesh_examples появляется stl файл:
Снимок экрана. Blender: pymesh_example_04_3.stl


Снимок экрана. Blender: pymesh_example_04_3.stlЕсли запросить фрактал 4-го уровня
Построение займет около 4-х часов, а размер файла составит 73 мб:
Снимок экрана. Blender: pymesh_example_04_4.stl


Снимок экрана. Blender: pymesh_example_04_4.stl
Menger sponge 4-го уровня, напечатанный на 3D принтере
Menger sponge 4-го уровня, напечатанный на 3D принтере

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

Mesh мы уже поворачивали, на этот раз повернём камеру.
Пример: pymesh_example_05.py
Показать код
import pymesh
from myplot_vf import plot_vf
import imageio
import os
import datetime

gif_half_size = 16

sponge = pymesh.meshio.load_mesh('/pymesh_examples/pymesh_example_04_3.stl')

filenames = []
imageio_images = []

# rendering
for i in range(gif_half_size):
vertices = sponge.vertices
faces = sponge.faces
filename = '/pymesh_examples/pymesh_example_05_3_'+str(i)+'.png'
print(filename)
filenames.append(filename)
plot_vf(vertices, faces, cam_x=i, cam_y=i, filename = filename)

print('forward rendering..')
for i in range(gif_half_size):
print(datetime.datetime.now(), filenames)
imageio_images.append(imageio.imread(filenames))
i-=1

print('reverse rendering..')
i = gif_half_size-1
while i>0:
print(datetime.datetime.now(), filenames)
imageio_images.append(imageio.imread(filenames))
i-=1

# save gif
imageio.mimsave('/pymesh_examples/pymesh_example_05_3.gif', imageio_images)

# remove images
for filename in filenames:
os.remove(filename)
Запуск:
sh pymesh_example_05.sh
Один из кадров анимации


Один из кадров анимацииПоказать анимацию
Анимация: pymesh_example_05_3.gif


Анимация: pymesh_example_05_3.gif

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

Страница проекта
Документация
Установка
Github
Построение пирамиды. Пример: pytorch3d_example_01.py
Показать код
import torch
from pytorch3d.io import save_obj
from scipy import spatial
import numpy as np

if torch.cuda.is_available():
device = torch.device("cuda:0")
torch.cuda.set_device(device)
else:
device = torch.device("cpu")

with_grad = False

verts = torch.tensor(
[
[-3, -3, 0],
[+3, -3, 0],
[+3, +3, 0],
[-3, +3, 0],
[+0, +0, +3]
],
device=device,
dtype=torch.float32,
requires_grad=with_grad,
)

hull = spatial.ConvexHull(verts.cpu().numpy())
faces_data = hull.simplices

faces = torch.tensor(
faces_data,
device=device,
dtype=torch.int64,
)

save_obj('pytorch3d_example_01.obj', verts, faces)
Подход очень похож на аналогичный в numpy-stl. Но так, как предполагается работа на GPU, будем разделять понятия хоста и устройства. Хост - наш компьютер. Устройство - GPU карта. Если карты нет, пользоваться библиотекой тоже можно, тогда в роли GPU будет выступать CPU. У хоста и устройства своя память. Чтобы передать объект с хоста на устройство и обратно, необходимо произвести небольшой ритуал.
В примере ниже, сразу на устройстве описываем вершины, копируем их с устройства на хост. На основании вершин вычисляем грани. Сохраняем объект. Файл в формате obj можно импортировать в blender:
Снимок экрана. Blender: pytorch3d_example_01.obj


Снимок экрана. Blender: pytorch3d_example_01.obj
Обратите внимание на команду verts.cpu().numpy()
Вершины копируются с устройства на хост. Если вы работаете с GPU, каждое копирование будет замедлять работу алгоритма. В планировании архитектуры программы, количество операций копирования между хостом и устройством, по возможности, лучше свести к минимуму. Например, изначально имея на хосте список вершин, можно вычислить грани, не прибегая к копированию вершин с устройства на хост, как это будет сделано в следующем примере.
Пример: pytorch3d_example_02.py
Показать код
import torch
from pytorch3d.io import IO
from pytorch3d.structures.meshes import Meshes
from pytorch3d.structures import join_meshes_as_scene
from pytorch3d.transforms import RotateAxisAngle
from pytorch3d.transforms import Translate
from scipy import spatial
import numpy as np

if torch.cuda.is_available():
device = torch.device("cuda:0")
torch.cuda.set_device(device)
else:
device = torch.device("cpu")

with_grad = False

host_verts = [
[-3, -3, 0],
[+3, -3, 0],
[+3, +3, 0],
[-3, +3, 0],
[+0, +0, +3]
]

verts = torch.tensor(
host_verts,
device=device,
dtype=torch.float32,
requires_grad=with_grad,
)

hull = spatial.ConvexHull(host_verts)
faces = torch.tensor(
hull.simplices,
device=device,
dtype=torch.int64,
)

mesh_a = Meshes(verts=[verts], faces=[faces])

rot = RotateAxisAngle(180,'X', device=device)
verts_b = rot.transform_points(mesh_a.verts_list()[0])

trans = Translate(0,0,5, device=device)
verts_b = trans.transform_points(verts_b)

mesh_b = Meshes(verts=[verts_b], faces=[faces])

mesh = join_meshes_as_scene([mesh_a, mesh_b])
IO().save_mesh(mesh, 'pytorch3d_example_02.obj')
Снимок экрана. Blender: pytorch3d_example_02.obj


Снимок экрана. Blender: pytorch3d_example_02.obj

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

Использование GPU даёт некоторый прирост производительности.
Пример: pytorch3d_example_03.py
Показать код
import torch
from pytorch3d.io import IO
from pytorch3d.structures.meshes import Meshes
from pytorch3d.structures import join_meshes_as_scene
from pytorch3d.transforms import Translate
from scipy import spatial
import numpy as np
import sys
import datetime

batch_size = 10

time_start = datetime.datetime.now()

def log(message):
print(datetime.datetime.now().strftime('%Y.%m.%d %H:%M:%S'), message)

if torch.cuda.is_available():
device = torch.device("cuda:0")
torch.cuda.set_device(device)
else:
device = torch.device("cpu")

depth = 3
side = 1/(3**depth)
with_grad = False

# vertices
vertices_base = [
[0, 0, 0],
[side, 0, 0],
[side, side, 0],
[0, side, 0],
[0, 0, side],
[side, 0, side],
[side, side, side],
[0, side, side]
]
vertices = torch.tensor(
vertices_base,
device=device,
dtype=torch.float32,
requires_grad=with_grad,
)

# faces
#hull = spatial.ConvexHull(vertices_base)
faces_data = [
[3, 2, 1],
[3, 1, 0],
[0, 1, 5],#
[5, 4, 0],
[7, 3, 0],
[0, 4, 7],#
[1, 2, 6],#
[6, 5, 1],
[2, 3, 6],#
[3, 7, 6],#
[4, 5, 6],#
[6, 7, 4]
]
faces = torch.tensor(
faces_data,
device=device,
dtype=torch.int64,
)

# holes
holes = []
side = 1
round_factor = 14
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):
holes.append({
'x':[round(side*x,round_factor), round(side*x+side,round_factor)],
'y':[round(side*y,round_factor), round(side*y+side,round_factor)]
})

def voxel_in_hole(holes, x, y):
comp_error = 0
for hole in holes:
if x + comp_error >= hole['x'][0] and x - comp_error < hole['x'][1] and \
y + comp_error >= hole['y'][0] and y - comp_error < hole['y'][1]:
return True
return False

# voxels
voxels = []
batch = 0
mesh = None
for x in range(3**depth):
log('x '+ str(x)+' in '+str(3**depth))
for y in range(3**depth):
for z in range(3**depth):
if voxel_in_hole(holes, round(y/(3**depth),round_factor), round(z/(3**depth),round_factor)) or \
voxel_in_hole(holes, round(x/(3**depth),round_factor), round(y/(3**depth),round_factor)) or \
voxel_in_hole(holes, round(x/(3**depth),round_factor), round(z/(3**depth),round_factor)):
continue
voxel_translate = Translate(x/(3**depth), y/(3**depth), z/(3**depth), device=device)
voxel_vertices_translated = voxel_translate.transform_points(vertices)
voxels.append(Meshes(verts=[voxel_vertices_translated], faces=[faces]))
batch += 1
if batch > batch_size:
if not mesh is None:
log('append mesh')
voxels.append(mesh)
log('join')
mesh = join_meshes_as_scene(voxels)
voxels = []
batch = 0


log('join')
if len(voxels):
if not mesh is None:
voxels.append(mesh)
mesh = join_meshes_as_scene(voxels)
log('save')
IO().save_mesh(mesh, 'pytorch3d_example_03.obj')
log('end')

time_end = datetime.datetime.now()
print('spent', (time_end - time_start).seconds, 'seconds')
В этом скрипте объявляем вершины минимального для указанной глубины вокселя. По знакомому из прошлого примера алгоритму вычисляем координаты отверстий в двух измерениях. Заполняем первичный куб вокселями, которые не попадают в пределы отверстий.
Снимок экрана. Blender: pytorch3d_example_03.obj


Снимок экрана. Blender: pytorch3d_example_03.obj
Cкорость расчета увеличилась на порядок, что позволило примерно за 5 часов построить фрактал 5 уровня:
Снимок экрана. Blender Menger Sponge 5 lvl


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

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

Помимо plotly визуализаций, pytorch3d отдельно выделяет рендеринг и подход нему тут довольно основательный, с текстурами и шейдерами.
Пример: pytorch3d_example_04.py
Показать код
import torch
from pytorch3d.structures.meshes import Meshes
from pytorch3d.io import IO, load_obj
from pytorch3d.renderer import (
FoVPerspectiveCameras, look_at_view_transform,
RasterizationSettings, BlendParams,
MeshRenderer, MeshRasterizer, HardPhongShader, TexturesVertex
)
import matplotlib.pyplot as plt
import imageio
import os

gif_half_size = 30

if torch.cuda.is_available():
device = torch.device("cuda:0")
torch.cuda.set_device(device)
else:
device = torch.device("cpu")

# Load object
verts, faces, aux = load_obj('pytorch3d_example_03.obj')

# init texture
verts_list = []
faces_list = []
texts_list = []

verts = torch.as_tensor(verts, dtype=torch.float32, device=device)
faces = torch.as_tensor(faces.verts_idx, dtype=torch.float32, device=device)
texts = torch.rand((len(verts), 3), dtype=torch.float32, device=device)
verts_list.append(verts)
faces_list.append(faces)
texts_list.append(texts)

# create textures
tex = TexturesVertex(texts_list)

# Initialise the mesh with textures
mesh = Meshes(verts=[verts], faces=[faces], textures=tex)[0]

filenames = []
imageio_images = []

# rendering
for i in range(gif_half_size):

filename = 'pytorch3d_example_'+str(i)+'.png'
print(filename)
filenames.append(filename)

# Select the viewpoint using spherical angles
distance = 3.5 # distance from camera to the object
elevation = 0.0 + i # angle of elevation in degrees
azimuth = 0.0 + i # No rotation so the camera is positioned on the +Z axis.

# Get the position of the camera based on the spherical angles
R, T = look_at_view_transform(distance, elevation, azimuth, device=device)

# Initialize an OpenGL perspective camera.
cameras = FoVPerspectiveCameras(device=device, R=R, T=T)

# Define the settings for rasterization and shading. Here we set the output image to be of size
# 512x512. As we are rendering images for visualization purposes only we will set faces_per_pixel=1
# and blur_radius=0.0. Refer to rasterize_meshes.py for explanations of these parameters.
raster_settings = RasterizationSettings(
image_size=512,
blur_radius=0.0,
faces_per_pixel=1,
)

# Create a phong renderer by composing a rasterizer and a shader. Here we can use a predefined
# PhongShader, passing in the device on which to initialize the default parameters
renderer = MeshRenderer(
rasterizer=MeshRasterizer(cameras=cameras, raster_settings=raster_settings),
shader=HardPhongShader(device=device, cameras=cameras)
)

images = renderer(mesh)
plt.figure(figsize=(10, 10))
plt.imshow(images[0, ..., :3].cpu().numpy())
plt.savefig(filename)
plt.axis("off")
imageio_images.append(imageio.imread(filename))


# reverse rendering
i = gif_half_size-1
while i>0:
print(filenames)
imageio_images.append(imageio.imread(filenames))
i-=1

# save gif
imageio.mimsave('pytorch3d_example_03.gif', imageio_images)

# remove images
for filename in filenames:
os.remove(filename)
Один из кадров анимации


Один из кадров анимацииПоказать анимацию

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.
Снимок экрана. Openscad: solidpython_example_01.scad


Снимок экрана. Openscad: solidpython_example_01.scad
solidpython удобно отлаживать. С одной стороны экрана открыт scad файл, с другой jupyter notebook. При исполнении scad_render_to_file картинка в openscad автоматически обновляется.
Если нужен stl, openscad умеет рендерить файлы этого формата через команды консоли. Пример вызова из jupyter notebook:
!openscad solidpython_example_01.scad -o solidpython_example_01.stl
Снимок экрана. Blender: solidpython_example_01.stl


Снимок экрана. Blender: 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.
Снимок экрана. Openscad: solidpython_example_02.scad


Снимок экрана. Openscad: solidpython_example_02.scad

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

Пример: solidpython_example_03.ipynb
Показать код
from solid import *
from IPython.display import Image
import datetime

depth = 3
sponge = 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 = cube([side,side,z], center = True)
box_a = translate([side*x-0.5+side/2, side*y-0.5+side/2, 0]) (box_a)
box_b = rotate([90,0,0]) (box_a)
box_c = rotate([0,90,0]) (box_a)
sponge -= box_a
sponge -= box_b
sponge -= box_c
path = scad_render_to_file(sponge, 'solidpython_example_03.scad')

# export to stl
!openscad solidpython_example_03.scad -o solidpython_example_03.stl

time_end = datetime.datetime.now()
print('spent', (time_end - time_start).seconds, 'seconds')
Снимок экрана. Openscad: solidpython_example_03.scad


Снимок экрана. Openscad: solidpython_example_03.scad

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

Сформируем серию изображений последней сцены, поворачивая камеру в каждом изображении.
Пример: solidpython_example_04.ipynb
Показать код
from solid import *
import subprocess
import datetime
import imageio
import os

gif_half_size = 30
filenames = []
imageio_images = []
for i in range(gif_half_size):
filename = 'solidpython_example_03_'+str(i)+'.png'
filenames.append(filename)
# export to png with camera params
# camera parameters when exporting png:
# =translate_x,y,z,rot_x,y,z,dist or
# =eye_x,y,z,center_x,y,z
rot_x = str(i-20)
rot_y = str(i-20)
rot_z = str(i-20)
dist = str(5+i/10)
cmd = "openscad solidpython_example_03.scad --camera '0,0,0,"+\
rot_x+","+rot_y+","+rot_z+","+dist+"' -o "+filename
MyOut = subprocess.Popen(
cmd,
shell=True,
stdout=subprocess.PIPE,
stderr=subprocess.PIPE,
stdin=subprocess.PIPE
)
stdout,stderr = MyOut.communicate()
imageio_images.append(imageio.imread(filename))
print(datetime.datetime.now(), filename)
print('stdout', stdout)
print('stderr', stderr)

print('reverse rendering..')
i = gif_half_size-1
while i>0:
print(datetime.datetime.now(), filenames)
imageio_images.append(imageio.imread(filenames))
i-=1

# save gif
imageio.mimsave('solidpython_example_03.gif', imageio_images)

# remove images
for filename in filenames:
os.remove(filename)

print('job complete!')
Один из кадров анимации


Один из кадров анимацииПоказать анимацию
m9f3ojaysevcjpqoohud7drifzs.gif

Кроме того, 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.
Библиотека​
Производительность
(Время вычисления фрактала 3-го уровня, в секунду)
Форматы импортируемых файлов​
Форматы экспортируемых файлов​
Отличия​
Numpy-stl​
24​
stl​
stl​
Pymesh​
823​
obj, off, ply, stl, mesh, msh, node, face, ele​
obj, off, ply, stl, mesh, msh, node, face, ele​
Pytorch3d​
12​
obj, ply, npz​
obj, ply​
GPU ускорение, Загрузка облака точек
SolidPython​
347​
scad (через консоль: stl, 3mf, off, amf, dxf, svg, csg)​
scad (через консоль: stl, 3mf, off, amf, dxf, svg, csg, png, stl, off, echo, ast, term, nef3, nefdbg)​
Большой выбор методов​
Характеристики моей машиныПолезные ссылки

 
Сверху