Допустим, вам потребовалось на языке программирования python, построить трёхмерную модель некоторого объекта, затем визуализировать его, или подготовить файл для печати на 3D принтере. Существует несколько библиотек, помогающих в решении этих задач. Поговорим о том, как строить трёхмерные модели из точек, граней и примитивов в python. Как выполнять элементарные приемы 3D моделирования: перемещение, поворот, объединение, вычитание и другие.
Рассмотрим библиотеки:
Примеры кода можно найти в репозитории на GitHub.
Подготовка среды
sudo sh get-docker.sh
sudo apt install g++
cd python3d
pip install -r requirements.txt
pip install "git+https://github.com/facebookresearch/pytorch3d.git"
Документация
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
Метод 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-')
Добавлено две грани: одна для нижней точки, одна для верхней
Описание логики построения фрактала
Договоимся что длина стороны грани фрактала - единица. Глубина фрактала это, грубо говоря, количество уникальных размеров отверстий. Длина вокселя зависит от глубины фрактала, она делится на 3 с каждым новым уровнем глубины.
Найдем сторону вокселя на глубинах 1 и 2. Упростим задачу, свернув фрактал с 3-х до 1-мерного случая:
Если фрактал второго уровня, длина стороны куба составит 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
Пример: 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
Github
Установка
К сожалению, библиотека PyMesh не установилась у меня ни через менеджер пакетов PIP (несмотря на упоминание этого способа в документации), ни через Anaconda. Есть два способа установки.
Начнем с простого куба. В папке 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
Что тут происходит?
Снимок экрана. 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
Функция boolean проста. Первым параметром передаётся объект из которого вычитаем. Вторым - который вычитаем. Следом передаются операция и движок.
Здесь boolean применяется не только для вычитания. Всего доступно 4 операции:
Иллюстрация применения операции 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
Перечисленных инструментов достаточно для построения фрактала.
Показать код
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Если запросить фрактал 4-го уровня
Построение займет около 4-х часов, а размер файла составит 73 мб:
Снимок экрана. Blender: pymesh_example_04_4.stl
Menger sponge 4-го уровня, напечатанный на 3D принтере
Пример: 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
Документация
Установка
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
Обратите внимание на команду 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
Пример: 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
Cкорость расчета увеличилась на порядок, что позволило примерно за 5 часов построить фрактал 5 уровня:
Снимок экрана. Blender Menger Sponge 5 lvl
Размер stl файла составил 1.9 ГБ. При построении фрактала 5-го уровня, программа останавливалась из-за переполнения памяти видеокарты. Пришлось сборку объекта выполнять пакетами. Создавалось по 10 слоев “двумерных” фракталов, затем они присоединялись к основному объекту, до тех пор, пока не построился полный фрактал.
Пример: 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)
Один из кадров анимацииПоказать анимацию
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
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
Общий принцип такой: Любая функция возвращает объект. Если над объектом необходимо произвести некоторое действие, объект (или список объектов) передается в круглых скобках после вызова соответствующей функции.
Пример: 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
Показать код
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
Пример: 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!')
Один из кадров анимацииПоказать анимацию
Кроме того, 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
Результат
Характеристики моей машиныПолезные ссылки
habr.com
Рассмотрим библиотеки:
- numpy-stl
- pymesh
- pytorch3d
- SolidPython
Примеры кода можно найти в репозитории на GitHub.
Подготовка среды
- Docker для исполнения примеров Pymesh:
sudo sh get-docker.sh
- Компилятор g++ для PyTorch3d:
sudo apt install g++
- Клонирование репозитория и установка библиотек:
cd python3d
pip install -r requirements.txt
pip install "git+https://github.com/facebookresearch/pytorch3d.git"
- Установка 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
Метод 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 пойдем от обратного. В нашем распоряжении два метода:- Построение элементарного mesh куба. Назовем его voxel.
- Объединение нескольких voxel в единый mesh.
Описание логики построения фрактала
Договоимся что длина стороны грани фрактала - единица. Глубина фрактала это, грубо говоря, количество уникальных размеров отверстий. Длина вокселя зависит от глубины фрактала, она делится на 3 с каждым новым уровнем глубины.
Найдем сторону вокселя на глубинах 1 и 2. Упростим задачу, свернув фрактал с 3-х до 1-мерного случая:
Если фрактал второго уровня, длина стороны куба составит 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
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
PyMesh: Обзор библиотеки
ДокументацияGithub
Установка
К сожалению, библиотека PyMesh не установилась у меня ни через менеджер пакетов PIP (несмотря на упоминание этого способа в документации), ни через Anaconda. Есть два способа установки.
- Следуя инструкции, собрать из исходников.
- Используя Docker контейнер. Выберем этот вариант, как более интересный. Запуск контейнера будет инициирован с параметрами. При помощи параметров запуска смонтируем в контейнер папку скриптов. Передадим необходимые параметры в скрипт. После завершения работы скрипта, контейнер будет удален. Если Docker вы уже установили следуя инструкции в начале статьи, более ничего устанавливать не нужно.
Начнем с простого куба. В папке 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
Квадратное отверстие сделаем при помощи булева вычитания. Строим параллелепипед и вычитаем его из основного куба.
Пример: 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
Функция boolean проста. Первым параметром передаётся объект из которого вычитаем. Вторым - который вычитаем. Следом передаются операция и движок.
Здесь boolean применяется не только для вычитания. Всего доступно 4 операции:
- Intersection: A∩B
- Union: A∪B
- Difference: A∖B (два последних примера)
- Symmetric difference: A XOR B (изображение не представлено)
Чтобы лучше понять как перемещать и вращать объект, бывает удобно временно заменить операцию 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
Перечисленных инструментов достаточно для построения фрактала.
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Если запросить фрактал 4-го уровня
Построение займет около 4-х часов, а размер файла составит 73 мб:
Снимок экрана. Blender: pymesh_example_04_4.stl
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
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
Обратите внимание на команду 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
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
Cкорость расчета увеличилась на порядок, что позволило примерно за 5 часов построить фрактал 5 уровня:
Снимок экрана. 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: Обзор библиотеки
Документация SolidPythonGithub
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
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
Общий принцип такой: Любая функция возвращает объект. Если над объектом необходимо произвести некоторое действие, объект (или список объектов) передается в круглых скобках после вызова соответствующей функции.
Пример: 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
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
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!')
Один из кадров анимацииПоказать анимацию
Кроме того, 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) | Большой выбор методов |
- 3D ML. Часть 1: формы представления 3D-данных
- 3D ML. Часть 2: функции потерь в задачах 3D ML
- 3D ML. Часть 3: датасеты и фреймворки в 3D ML
- 3D ML. Часть 4: дифференциальный рендеринг
- 3D ML. Часть 5: Свертки на графах
- 3D ML. Часть 6: Обзор алгоритмов семантической сегментации облака точек
- PyTorch3D : RoboLove
- trimesh - Страница проекта
- Работа с 3D моделями в Python с использованием библиотеки OpenMesh
- openscadpy - Страница проекта
- Система скриптового 3д моделирования ZenCad
- zenscad - Документация
- Hands-on Guide to PyTorch 3D – A Library for Deep Learning with 3D Data
- Scripting the Unreal Editor using Python
- Nvidia kaolin
- Guide To NVIDIA’s Kaolin: A 3D Deep Learning Library
- 3D Рендеринг в matplotlib
3D моделирование в Python
Допустим, вам потребовалось на языке программирования python, построить трёхмерную модель некоторого объекта, затем визуализировать его, или подготовить файл для печати на 3D принтере. Существует...
habr.com