![](https://habrastorage.org/webt/7o/la/oc/7olaocvqltjf9qgmxe-rpbhftva.png)
Допустим, вам потребовалось на языке программирования 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: Обзор библиотеки
Строение полигональной модели:
![Изображение из Википедии Изображение из Википедии](https://habrastorage.org/webt/vw/kl/rp/vwklrpj8lp8bomlojmf5qcyq6qs.png)
Vertices (вершины) — список точек. Каждая точка описана тремя числами — координатами в 3-х мерном пространстве.
Если вы не знакомы с Jupyter notebook
Пример: numpy_stl_example_01.ipynb
import numpy as np from myplot import plot_verticles vertices = np.array( [ [-3, -3, 0], [+3, -3, 0], [+3, +3, 0], [-3, +3, 0], [+0, +0, +3] ] ) plot_verticles(vertices = vertices, isosurf = False)
![Вершины пирамиды Вершины пирамиды](https://habrastorage.org/webt/k3/dz/uv/k3dzuvmjw1tj_ukifsgez_fr-jc.png)
Несмотря на то, что пока описаны только вершины, уже можно взглянуть как будет выглядеть модель, если соединить их треугольниками:
plot_verticles(vertices = vertices, isosurf = True)
![Изоповерхность пирамиды Изоповерхность пирамиды](https://habrastorage.org/webt/bn/tf/ys/bntfysalvm3z0s9jakzpygv-rjw.png)
Выглядит так, будто грани уже существуют. Но пока у нас есть только вершины. Чтобы сформировать 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:
![Принцип описания граней через позиции вершин Принцип описания граней через позиции вершин](https://habrastorage.org/webt/zx/1i/xs/zx1ixspulzbw1ug-adfcl1hcgno.png)
Mesh (сетка) — Набор вершин и граней, которые определяют форму многогранного объекта.
myramid_mesh = mesh.Mesh(np.zeros(faces.shape[0], dtype=mesh.Mesh.dtype)) for i, f in enumerate(faces): for j in range(3): myramid_mesh.vectors[i][j] = vertices[f[j],:] plot_mesh(myramid_mesh)
![Грани пирамиды Грани пирамиды](https://habrastorage.org/webt/ia/aj/ix/iaajixa6kvglnaa3ncdtgnvgxhg.png)
Как видно из рисунка, одна грань пирамиды оказалась перевернута. В последующих примерах, при построении фракталов, метод ConvexHull применяться не будет, так как он располагает точки грани в произвольном порядке, что приводит к переворачиванию некоторых граней.
myramid_mesh.save('numpy_stl_example_02.stl')
Для просмотра stl файлов разработано довольно много программ. Одна из них называется Blender, доступна для скачивания и не требует оплаты за использование
![Снимок экрана. Blender: numpy_stl_example_02.stl Снимок экрана. Blender: numpy_stl_example_02.stl](https://habrastorage.org/webt/gb/dp/hp/gbdphpl7lsqk88chsnjgay0ygtq.png)
Метод 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-')
![Центральные точки не связаны с гранями Центральные точки не связаны с гранями](https://habrastorage.org/webt/ys/mf/yr/ysmfyr94j59btacq1wuhd0gsxqa.png)
Для таких случаев придется найти альтернативу 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-')
![Добавлено две грани: одна для нижней точки, одна для верхней Добавлено две грани: одна для нижней точки, одна для верхней](https://habrastorage.org/webt/iz/ye/ky/izyekyl0v3c0jwthdh_xscuokhc.png)
Numpy-stl: Построение фрактала
Пришло время построить фрактал. В numpy-stl нет функции булева вычитания. Для построения фрактала Menger Sponge пойдем от обратного. В нашем распоряжении два метода:
-
Построение элементарного mesh куба. Назовем его voxel.
-
Объединение нескольких voxel в единый mesh.
Построим фрактал из кубиков, как из конструктора.
Описание логики построения фрактала
Договоимся что длина стороны грани фрактала — единица. Глубина фрактала это, грубо говоря, количество уникальных размеров отверстий. Длина вокселя зависит от глубины фрактала, она делится на 3 с каждым новым уровнем глубины.
Найдем сторону вокселя на глубинах 1 и 2. Упростим задачу, свернув фрактал с 3-х до 1-мерного случая:
![](https://habrastorage.org/webt/xa/lb/rr/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[i][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](https://habrastorage.org/webt/wd/cx/-l/wdcx-ljjdtbvjsaln2u9974uonq.png)
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[i]) imageio_images.append(imageio.imread(filenames[i])) i-=1 print('reverse rendering..') i = gif_half_size-1 while i>0: print(datetime.datetime.now(), filenames[i]) imageio_images.append(imageio.imread(filenames[i])) i-=1 # save gif imageio.mimsave('solidpython_example_03.gif', imageio_images) # remove images for filename in filenames: os.remove(filename) print('job complete!')
![Один из кадров анимации Один из кадров анимации](https://habrastorage.org/webt/uy/9y/xg/uy9yxg3j1eqbxgwhiohzuklvliu.png)
Показать анимацию
![Анимация: solidpython_example_03.gif Анимация: solidpython_example_03.gif](https://habrastorage.org/webt/kz/fh/q6/kzfhq6y_wnkd4lxmuhmxsfw-30c.gif)
PyMesh: Обзор библиотеки
Установка
К сожалению, библиотека PyMesh не установилась у меня ни через менеджер пакетов PIP (несмотря на упоминание этого способа в документации), ни через Anaconda. Есть два способа установки.
-
Следуя инструкции, собрать из исходников.
-
Используя 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](https://habrastorage.org/webt/9y/tv/wb/9ytvwbrdazkgpfkpfvtixdrvd7q.png)
Квадратное отверстие сделаем при помощи булева вычитания. Строим параллелепипед и вычитаем его из основного куба.
Пример: 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](https://habrastorage.org/webt/jv/l5/by/jvl5byacrheodcwux0dcg5g_2oy.png)
Функция boolean проста. Первым параметром передаётся объект из которого вычитаем. Вторым — который вычитаем. Следом передаются операция и движок.
Здесь boolean применяется не только для вычитания. Всего доступно 4 операции:
-
Intersection: A∩B
-
Union: A∪B
-
Difference: A∖B (два последних примера)
-
Symmetric difference: A XOR B (изображение не представлено)
![Иллюстрация применения операции boolean на сфере и кубе Иллюстрация применения операции boolean на сфере и кубе](https://habrastorage.org/webt/qf/x_/wj/qfx_wj6kuzjvk_y1gtdao7rmft0.png)
Чтобы лучше понять как перемещать и вращать объект, бывает удобно временно заменить операцию 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](https://habrastorage.org/webt/iz/yv/bp/izyvbpevgh5gi9awayb6krlwhto.png)
Перечисленных инструментов достаточно для построения фрактала.
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](https://habrastorage.org/webt/g_/ae/l3/g_ael35hczajf025c9hno8cntm0.png)
Если запросить фрактал 4-го уровня
Построение займет около 4-х часов, а размер файла составит 73 мб:
![Снимок экрана. Blender: pymesh_example_04_4.stl Снимок экрана. Blender: pymesh_example_04_4.stl](https://habrastorage.org/webt/4u/if/dd/4uifddiphckpmtrgmd4a1e2qk6i.png)
![Menger sponge 4-го уровня, напечатанный на 3D принтере Menger sponge 4-го уровня, напечатанный на 3D принтере](https://habrastorage.org/webt/pi/ro/rd/pirordukfelh---4mfkmqtnecx8.jpeg)
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[i]) imageio_images.append(imageio.imread(filenames[i])) i-=1 print('reverse rendering..') i = gif_half_size-1 while i>0: print(datetime.datetime.now(), filenames[i]) imageio_images.append(imageio.imread(filenames[i])) 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
![Один из кадров анимации Один из кадров анимации](https://habrastorage.org/webt/7i/l2/fo/7il2fow-jhalj2fmhqjv3aqbz1c.png)
Показать анимацию
![Анимация: pymesh_example_05_3.gif Анимация: pymesh_example_05_3.gif](https://habrastorage.org/webt/-c/hr/zq/-chrzqyavbe-xepwt6vtwdth0u8.gif)
PyTorch3d: Обзор библиотеки
Построение пирамиды. Пример: 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](https://habrastorage.org/webt/fr/c1/p_/frc1p_hvt-fft5df-h-vthny36a.png)
Обратите внимание на команду 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](https://habrastorage.org/webt/y1/of/ec/y1ofec-d-h34f9uvm6h73haypw0.png)
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](https://habrastorage.org/webt/6u/62/yb/6u62yb4tpzst1wd1ofvttvhzjxa.png)
Cкорость расчета увеличилась на порядок, что позволило примерно за 5 часов построить фрактал 5 уровня:
![Снимок экрана. Blender Menger Sponge 5 lvl Снимок экрана. Blender Menger Sponge 5 lvl](https://habrastorage.org/webt/yw/oy/41/ywoy41nayiuh2pd1p5ujnhd_guq.png)
Размер 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[i]) imageio_images.append(imageio.imread(filenames[i])) i-=1 # save gif imageio.mimsave('pytorch3d_example_03.gif', imageio_images) # remove images for filename in filenames: os.remove(filename)
![Один из кадров анимации Один из кадров анимации](https://habrastorage.org/webt/sp/uy/n_/spuyn_vtca3jk0jorz29hzr-dpq.png)
Показать анимацию
![](https://habrastorage.org/webt/lu/hb/rn/luhbrnspu0ep1saufkpl7servzw.gif)
SolidPython: Обзор библиотеки
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](https://habrastorage.org/webt/ln/il/hf/lnilhfl--6f5ib_df8mm0sxw0ie.png)
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](https://habrastorage.org/webt/nb/26/zu/nb26zulkx-swb3rvsk-tv6qkx7g.png)
Общий принцип такой: Любая функция возвращает объект. Если над объектом необходимо произвести некоторое действие, объект (или список объектов) передается в круглых скобках после вызова соответствующей функции.
Пример: 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](https://habrastorage.org/webt/tj/bj/ak/tjbjakjrhqsg8f7jahx2-i7xnw8.png)
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](https://habrastorage.org/webt/p3/dg/l6/p3dgl62d4xqpm82akshpusnbymo.png)
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[i]) imageio_images.append(imageio.imread(filenames[i])) i-=1 # save gif imageio.mimsave('solidpython_example_03.gif', imageio_images) # remove images for filename in filenames: os.remove(filename) print('job complete!')
![Один из кадров анимации Один из кадров анимации](https://habrastorage.org/webt/gi/c-/4v/gic-4voe2nhe4zudn2kzqbvnayk.png)
Показать анимацию
![](https://habrastorage.org/webt/m9/f3/oj/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
![Результат Результат](https://habrastorage.org/webt/6g/pe/an/6gpeans9u0ifehnafwu6kbdjtdi.png)
Сравнение библиотек
Сравнение производительности не совсем объективно, так как имеются значительные различия в алгоритмах. В 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) |
Большой выбор методов |
Характеристики моей машины
-
Ubuntu 20.04.2 LTS
-
Python 3.7.3
-
Intel(R) Core(TM) i3-7350K CPU @ 4.20GHz
-
RAM 35.2Gb
-
GPU GeForce GTX 1080 Ti 11175 MB
ссылка на оригинал статьи https://habr.com/ru/articles/572760/
Добавить комментарий