Введение
Сегодня мне потребовалось нарисовать сферу из треугольников как часть 3D визуализации для научного проекта. До этого я всегда экспортировал данные для визуализации в VTK с помощью замечательной библиотеки pyevtk [3], которая за меня строила довольно сложный XML файл, который я потом рисовал в ParaView.
На удивление, гугление "evtk unstructured grid example" давало только ссылку на открытый баг [2]; Яндекс тоже ничего не нашёл.
Взятый за основу файл с демонстрацией бага послужил достаточно хорошим толчком, спасибо Maxim Grechkin [2]. Ниже будет описан процесс правильного создания корректного VTK-файла с неструктурированными сетками, который годится для ParaView и товарищей.
Воспроизведение бага
Скачиваем файл Максима с исходниками и смотрим код (я привёл код в соответствие PEP8 требованиям и добавил все недостающие импорты):
Видно, что человек собирает файл из кусочков, т.е. сам создаёт и пишет разделы VTK файла, используя низкоуровневые функции.
import numpy as np from evtk.vtk import VtkFile, VtkUnstructuredGrid, VtkTriangle from evtk.hl import pointsToVTK def dump_vtu_evtk(filename, x, y, faces, ux, uy): vec_speed = (ux, uy, np.zeros_like(ux)) vertices = (x, y, np.zeros_like(x)) w = VtkFile(filename, VtkUnstructuredGrid) w.openGrid() w.openPiece(npoints=len(x), ncells=len(faces)) w.openData("Cell", vectors="Velocity") w.addData("Velocity", vec_speed) w.closeData("Cell") w.openElement("Points") w.addData("Points", vertices) w.closeElement("Points") ncells = len(faces) # index of last node in each cell offsets = np.arange(start=3, stop=3*(ncells + 1), step=3, dtype='uint32') print offsets, len(offsets) connectivity = faces.reshape(ncells*3).astype('int32') + 1 print connectivity, len(connectivity) cell_types = np.empty(ncells, dtype='uint8') cell_types[:] = VtkTriangle.tid print cell_types, len(cell_types)
Видно, что человек собирает файл из кусочков, т.е. сам создаёт и пишет разделы VTK файла, используя низкоуровневые функции.
Пробуем запустить файл с простыми данными:
arr = np.array dump_vtu_evtk("demo_maxim", x=arr([-1.0, 0, 1]), y=arr([+0.0, 1, 0]), ux=arr([0, 0, 1.0]), uy=arr([1, 1, 0.0]), faces=arr([[0, 1, 2]]))
Файл создаётся, но в ParaView открываться не хочет:
Вроде да, баг.
Делаем над собой насилие и читаем весь баг-трекер. Другой хороший человек, Panagiotis Mavrogiorgos, тоже боролся с задачей и приложил свой файл points.vtu, сгенерированный переработанной версией функции Максима. Она тоже с багами, правда. Скачиваем, пробуем открыть, читаем сообщение ParaView внимательно.
Видно, что тут уже на синтаксические ошибки не жалуются, но ParaView не хватает информации о connectivity.
Видно, что тут уже на синтаксические ошибки не жалуются, но ParaView не хватает информации о connectivity.
К счастью, evtk поставляется с примерами. Ближайший аналог решаемой задачи (создания неструктурированной сетки) ─ это функция evtk.hl.pointsToVTK(...), которая пишет уже полностью корректный файл:
pointsToVTK("demo_evtk", x=arr([-1.0, 0, 1]), y=arr([0, 1.0, 0]), z=arr([0, 0, 0.0]), data={"rho": arr([1, 1, 1.0])})
Пробуем. Да, генерирует нормальный файл, который спокойно читается и рисуется.
Как настоящие детективы,наконец-то открываем оба vtu-файла в текстовом редакторе и видим, что файл Максима без ног и одной почки содержит незакрытые секции. Это с точки зрения XML (не говоря уже про VTK) весьма печально.
Файл с ошибками:
Как настоящие детективы,
Файл с ошибками:
<?xml version="1.0"?> <VTKFile byte_order="LittleEndian" version="0.1" type="UnstructuredGrid"> <UnstructuredGrid> <Piece NumberOfCells="1" NumberOfPoints="3"> <CellData vectors="Velocity"> <DataArray NumberOfComponents="3" offset="0" type="Float64" Name="Velocity" format="appended"/> </CellData> <Points> <DataArray NumberOfComponents="3" offset="76" type="Float64" Name="Points" format="appended"/> </Points>
Начинаем думать, какие XML разделы нам нужны и как они создаются со стороны evtk. В конечном итоге обращаемся к
Действительно, в секцию Cells надо добавить три DataArray: connectivity, offset и types, которые есть в правильном файле и без которых непонятно, как складывать вершины в грани. Эти разделы показаны на рисунке ниже (... обозначает сокращённые места):
Получается, что нам надо дополнить функцию Максима записью метаинформации о массивах connectivity и offsets в соответствующий XML узел раздела Cells, после чего не забыть записать сами массивы как бинарные данные в монолитном блоке двоичного фарша, живущего пониже XML шапки. Так устроен формат, ничего не поделаешь.
Подглядывая в evtk.hl.pointsToVTK, пробуем:
Работает!#! /usr/bin/env python # -*- coding: UTF-8 -*- from evtk.vtk import VtkFile, VtkUnstructuredGrid, VtkTriangle, VtkVertex import numpy as np def _addDataToFile(vtkFile, cellData, pointData): # Point data if pointData is not None: keys = pointData.keys() vtkFile.openData("Point", scalars=keys[0]) for key in keys: data = pointData[key] vtkFile.addData(key, data) vtkFile.closeData("Point") # Cell data if cellData is not None: keys = cellData.keys() vtkFile.openData("Cell", scalars=keys[0]) for key in keys: data = cellData[key] vtkFile.addData(key, data) vtkFile.closeData("Cell") def _appendDataToFile(vtkFile, cellData, pointData): # Append data to binary section if pointData is not None: keys = pointData.keys() for key in keys: data = pointData[key] vtkFile.appendData(data) if cellData is not None: keys = cellData.keys() for key in keys: data = cellData[key] vtkFile.appendData(data) def dump_vtu_evtk(filename, x, y, faces, vx, vy, point_data, cell_data): w = VtkFile(filename, VtkUnstructuredGrid) w.openGrid() w.openPiece(npoints=len(x), ncells=len(faces)) vertices = (x, y, np.zeros_like(x)) w.openElement("Points") w.addData("Points", vertices) w.closeElement("Points") # Create some temporary arrays to write grid topology ... # ... index of last node in each cell ... ncells = len(faces) # index of last node in each cell offsets = np.arange(start=3, stop=3*(ncells + 1), step=3, dtype='uint32') # ... unwinding our triangles into connectivity list. connectivity = faces.reshape(ncells*3).astype('int32') cell_types = np.ones(ncells, dtype='uint8')*VtkTriangle.tid w.openElement("Cells") w.addData("connectivity", connectivity) w.addData("offsets", offsets) w.addData("types", cell_types) w.closeElement("Cells") _addDataToFile(w, cellData=cell_data, pointData=point_data) w.closePiece() w.closeGrid() w.appendData(vertices) w.appendData(connectivity).appendData(offsets).appendData(cell_types) _appendDataToFile(w, cellData=cell_data, pointData=point_data) w.save() return w.getFileName() a = np.array dump_vtu_evtk("demo", x=a([-1.0, 0, 1]), y=a([0.0, 1, 0]), faces=a([[0.0, 1, 2]]), vx=a([1.0]), vy=a([1.0]), point_data={'T': a([1.0, 2, 3]), 'u': (a([0, 0, 1.0]), a([1, 1, 0.0]), a([0, 1, 1.0]))}, cell_data={'rho': a([1]), 'v': (a([1.0]), a([2.0]), a([3.0]))})
Для демонстрации я написал маленький тесселятор сферы, который берёт базовый многогранник и, последовательно уменьшая треугольные грани в 2 раза, строит всё более пригодное приближение к сфере, как показано ниже на рисунке:
Файл исходника можно скачать [4] или взять из приложения (сделал на случай проблем с хостингом исходника). Нормального человека из всей этой простыни интересует только функция triangle_faces_to_VTK.
После визуализации в ParaView вы получите отличную сферу со скалярными и векторными данными в узлах (point data) и на гранях (cell data).
Замечание: для визуализации векторных cell data помогает фильтр "Cell To Point Data".
Картинки:
Файл исходника можно скачать [4] или взять из приложения (сделал на случай проблем с хостингом исходника). Нормального человека из всей этой простыни интересует только функция triangle_faces_to_VTK.
После визуализации в ParaView вы получите отличную сферу со скалярными и векторными данными в узлах (point data) и на гранях (cell data).
Замечание: для визуализации векторных cell data помогает фильтр "Cell To Point Data".
Картинки:
Заключение
Запись в двух местах с точки зрения программирования "попахивает", но это определяется самим VTK-фундаментом и линейным характером записи в файл ─ сперва пишется вся шапка, потом данные.
Что делать 1: есть соблазн сделать наследника класса с внутренней очередью, в которой собирать массивы на запись. Метод при получении нового массива сразу сможет писать мета-информацию и ставить массив в очередь, а в деструкторе делать саму запись всей очереди в двоичный 'raw' блок в гарантированно правильном порядке.
Что делать 1: есть соблазн сделать наследника класса с внутренней очередью, в которой собирать массивы на запись. Метод при получении нового массива сразу сможет писать мета-информацию и ставить массив в очередь, а в деструкторе делать саму запись всей очереди в двоичный 'raw' блок в гарантированно правильном порядке.
Здесь будет скрытая проблема: массивы при более-менее серьёзном моделировании могут занимать до сотен гигабайт, и рачительный программист может после "отдачи на запись" сгоряча их переиспользовать, что при передаче данных по ссылке испортит ещё не записанные данные внутри класса (так как запись данных фактически отложена до вызова деструктора).
Что делать 2: двоичные данные можно было бы писать в отдельный файл и потом делать некий include, но это не предусмотрено форматом, насколько я его вычитал.
Что делать 3: можно попробовать зарезервировать место на шапку и писать всё сразу по получении с помощью fseek/ftell, проверяя наложения.
Что делать 2: двоичные данные можно было бы писать в отдельный файл и потом делать некий include, но это не предусмотрено форматом, насколько я его вычитал.
Что делать 3: можно попробовать зарезервировать место на шапку и писать всё сразу по получении с помощью fseek/ftell, проверяя наложения.
Как всегда, решение принимать
Удачи!
[2] Описание открытого бага в репозитории pyevtk: https://bitbucket.org/pauloh/pyevtk/issues/1/unstructured-grid-output-doesnt-work
[3] Библиотека evtk для создания VTK файлов из Питона: https://bitbucket.org/pauloh/pyevtk
[4] Демонстрация экспорта в действии на примере тесселированной сферы https://gist.github.com/dromanov/0fb8bacff5342a56a690
Ссылки
[1] Описание формата VTK-файла: http://www.vtk.org/VTK/img/file-formats.pdf[2] Описание открытого бага в репозитории pyevtk: https://bitbucket.org/pauloh/pyevtk/issues/1/unstructured-grid-output-doesnt-work
[3] Библиотека evtk для создания VTK файлов из Питона: https://bitbucket.org/pauloh/pyevtk
[4] Демонстрация экспорта в действии на примере тесселированной сферы https://gist.github.com/dromanov/0fb8bacff5342a56a690
Приложение: полный исходник с демонстрацией экспорта векторных и скалярных point и cell данных.
Единственная ценная функция ─ это triangle_faces_to_VTK. Всё остальное служит задаче тесселяции сферы для генерации неструктурированной сетки для экспорта.
#! /usr/bin/env python # -*- coding: UTF-8 -*- from evtk.vtk import VtkFile, VtkUnstructuredGrid, VtkTriangle, VtkVertex import numpy as np from copy import deepcopy """ UnstructuredGrid ================ Each 'UnstructuredGrid' piece specifies a set of points and cells independently from the other pieces. The points are described explicitly by the 'Points' element. The cells are described explicitly by the 'Cells' element. <VTKFile type=”UnstructuredGrid” ...> <UnstructuredGrid> <Piece NumberOfPoints=”#” NumberOfCells=”#”> <PointData>...</PointData> <CellData>...</CellData> <Points>...</Points> <Cells>...</Cells> </Piece> </UnstructuredGrid> </VTKFile> Every dataset describes the data associated with its points and cells with 'PointData' and 'CellData' XML elements as follows: <PointData Scalars=”Temperature” Vectors=”Velocity”> <DataArray Name=”Velocity” .../> <DataArray Name=”Temperature” .../> <DataArray Name=”Pressure” .../> </PointData> VTK allows an arbitrary number of data arrays to be associated with the points and cells of a dataset. Each data array is described by a 'DataArray' element which, among other things, gives each array a name. The following attributes of 'PointData' and 'CellData' are used to specify the active arrays by name: 'Scalars' — The name of the active scalars array, if any. 'Vectors' — The name of the active vectors array, if any. 'Normals' — The name of the active normals array, if any. 'Tensors' — The name of the active tensors array, if any. 'TCoords' — The name of the active texture coordinates array, if any. Points ====== The 'Points' element explicitly defines coordinates for each point individually. It contains one 'DataArray' element describing an array with three components per value, each specifying the coordinates of one point. <Points> <DataArray NumberOfComponents=”3” .../> </Points> Coordinates =========== The 'Coordinates' element defines point coordinates for an extent by specifying the ordinate along each axis for each integer value in the extent’s range. It contains three 'DataArray' elements describing the ordinates along the x-y-z axes, respectively. <Coordinates> <DataArray .../> <DataArray .../> <DataArray .../> </Coordinates> Verts, Lines, Strips, and Polys =============================== The 'Verts', 'Lines', 'Strips', and 'Polys' elements define cells explicitly by specifying point connectivity. Cell types are implicitly known by the type of element in which they are specified. Each element contains two 'DataArray' elements. The first array specifies the point connectivity. All the cells’ point lists are concatenated together. The second array specifies the offset into the connectivity array for the end of each cell. <Verts> <DataArray type=”Int32” Name=”connectivity” .../> <DataArray type=”Int32” Name=”offsets” .../> </Verts> Cells ===== The 'Cells' element defines cells explicitly by specifying point connectivity and cell types. It contains three 'DataArray' elements. The first array specifies the point connectivity. All the cells’ point lists are concatenated together. The second array specifies the offset into the connectivity array for the end of each cell. The third array specifies the type of each cell. (Note: the cell types are defined in Figure 2 and Figure 3.) <Cells> <DataArray type=”Int32” Name=”connectivity” .../> <DataArray type=”Int32” Name=”offsets” .../> <DataArray type=”UInt8” Name=”types” .../> </Cells> All of the data and geometry specifications use 'DataArray' elements to describe their actual content as follows: DataArray ========= The 'DataArray' element stores a sequence of values of one type. There may be one or more components per value. <DataArray type=”Float32” Name=”vectors” NumberOfComponents=”3” format=”appended” offset=”0”/> <DataArray type=”Float32” Name=”scalars” format=”binary”> bAAAAAAAAAAAAIA/AAAAQAAAQEAAAIBA... </DataArray> <DataArray type=”Int32” Name=”offsets” format=”ascii”> 10 20 30 ... </DataArray> The attributes of the 'DataArray' elements are described as follows: type — The data type of a single component of the array. This is one of Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Float32, Float64. Note: the 64-bit integer types are only supported if VTK_USE_64BIT_IDS is on (a CMake variable—see “CMake” on page 8) or the platform is 64-bit. Name — The name of the array. This is usually a brief description of the data stored in the array. NumberOfComponents — The number of components per value in the array. format — The means by which the data values themselves are stored in the file. This is “ascii”, “binary”, or “appended”. offset — If the format attribute is “appended”, this specifies the offset from the beginning of the appended data section to the beginning of this array’s data. The format attribute chooses among the three ways in which data values can be stored: format=”ascii” — The data are listed in ASCII directly inside the DataArray element. Whitespace is used for separation. format=”binary” — The data are encoded in base64 and listed contiguously inside the DataArray element. Data may also be compressed before encoding in base64. The byte-order of the data matches that specified by the byte_order attribute of the 'VTKFile' element. format=”appended” — The data are stored in the appended data section. Since many 'DataArray' elements may store their data in this section, the 'offset' attribute is used to specify where each DataArray’s data begins. This format ("appended") is the default used by VTK’s writers. The appended data section is stored in an 'AppendedData' element that is nested inside 'VTKFile' after the dataset element: <VTKFile ...> ... <AppendedData encoding=”base64”> _QMwEAAAAAAAAA... </AppendedData> </VTKFile> The appended data section begins with the first character after the underscore inside the AppendedData element. The underscore is not part of the data, but is always present. Data in this section is always in binary form, but can be compressed and/or base64 encoded. The byte-order of the data matches that specified by the byte_order attribute of the 'VTKFile' element. Each DataArray’s data are stored contiguously and appended immediately after the previous DataArray’s data without a separator. The DataArray’s offset attribute indicates the file position offset from the first character after the underscore to the beginning its data. """ # These two functions are taken from original 'evtk.hl' module without changes. def _addDataToFile(vtkFile, cellData, pointData): # Point data if pointData is not None: keys = pointData.keys() vtkFile.openData("Point", scalars=keys[0]) for key in keys: data = pointData[key] vtkFile.addData(key, data) vtkFile.closeData("Point") # Cell data if cellData is not None: keys = cellData.keys() vtkFile.openData("Cell", scalars=keys[0]) for key in keys: data = cellData[key] vtkFile.addData(key, data) vtkFile.closeData("Cell") def _appendDataToFile(vtkFile, cellData, pointData): # Append data to binary section if pointData is not None: keys = pointData.keys() for key in keys: data = pointData[key] vtkFile.appendData(data) if cellData is not None: keys = cellData.keys() for key in keys: data = cellData[key] vtkFile.appendData(data) def triangle_faces_to_VTK(filename, x, y, z, faces, point_data, cell_data): vertices = (x, y, z) w = VtkFile(filename, VtkUnstructuredGrid) w.openGrid() w.openPiece(npoints=len(x), ncells=len(faces)) w.openElement("Points") w.addData("Points", vertices) w.closeElement("Points") # Create some temporary arrays to write grid topology. ncells = len(faces) # Index of last node in each cell. offsets = np.arange(start=3, stop=3*(ncells + 1), step=3, dtype='uint32') # Connectivity as unrolled array. connectivity = faces.reshape(ncells*3).astype('int32') cell_types = np.ones(ncells, dtype='uint8')*VtkTriangle.tid w.openElement("Cells") w.addData("connectivity", connectivity) w.addData("offsets", offsets) w.addData("types", cell_types) w.closeElement("Cells") _addDataToFile(w, cellData=cell_data, pointData=point_data) w.closePiece() w.closeGrid() w.appendData(vertices) w.appendData(connectivity).appendData(offsets).appendData(cell_types) _appendDataToFile(w, cellData=cell_data, pointData=point_data) w.save() return w.getFileName() def pointsToVTK(path, x, y, z, data): """ Export points and associated data as an unstructured grid. PARAMETERS: path: name of the file without extension where data to be saved. x, y, z: 1D arrays with coordinates of the points. data: dictionary {'varname': data_array} of point-data to export. RETURNS: Full path to saved file. """ assert (x.size == y.size == z.size) npoints = x.size # Create some temporary arrays to write grid topology ... # ... index of last node in each cell ... offsets = np.arange(start=1, stop=npoints + 1, dtype='int32') # ... unwinding our triangles into connectivity list. connectivity = np.arange(npoints, dtype='int32') cell_types = np.empty(npoints, dtype='uint8') cell_types[:] = VtkVertex.tid w = VtkFile(path, VtkUnstructuredGrid) w.openGrid() w.openPiece(ncells=npoints, npoints=npoints) w.openElement("Points") w.addData("points", (x, y, z)) w.closeElement("Points") w.openElement("Cells") w.addData("connectivity", connectivity) w.addData("offsets", offsets) w.addData("types", cell_types) w.closeElement("Cells") _addDataToFile(w, cellData=None, pointData=data) w.closePiece() w.closeGrid() w.appendData((x, y, z)) w.appendData(connectivity).appendData(offsets).appendData(cell_types) _appendDataToFile(w, cellData=None, pointData=data) w.save() return w.getFileName() class Level: """Approximation of the spherical surface made from triangles. It support refinement of the surface to increase detalization uniformly by decreasing face size 4 times on each step. Enumeration of vertices within face is shown below together with the numbers of new vertices introduced during refinement: 2 2 * * / \ . / \ / \ |\ / \ / \ +---------+ \ 5 / \ 4 / \ | * *-------* / \ +---------+ / / \ / \ / \ |/ / \ / \ / \ " / \ / \ *---------------* *-------*-------* 0 1 0 3 1 """ def __init__(self, faces, vertices): self.faces = faces self.vertices = vertices def refine(self): """Returns level with better sampling.""" cache = {} vert = deepcopy(self.vertices) def _(i, j): '''Makes cached middle vertex and returns its index.''' ij = tuple(sorted([i, j])) if ij not in cache: pos = len(vert) v = 0.5*(vert[i] + vert[j]) vert.append(v/np.sqrt(np.inner(v, v))) cache[ij] = pos return cache[ij] f_pos = 0 face = np.zeros([4*len(self.faces), 3], int) for f in self.faces: i, j, k = f face[f_pos+0, :] = [ i, _(i, j), _(i, k)] face[f_pos+1, :] = [_(i, j), j, _(j, k)] face[f_pos+2, :] = [_(i, j), _(j, k), _(i, k)] face[f_pos+3, :] = [ k, _(i, k), _(j, k)] f_pos += 4 return Level(face, vert) class Polyhedron: """http://student.ulb.ac.be/~claugero/sphere/index.html""" def __init__(self, level=3, base="tetrahedron"): """Creates set of normals.""" assert level in range(1, 8), "bad level: %s" % level if base == "tetrahedron": s = 1.0/np.sqrt(3.0); verts = [[ s, s, s], [-s, -s, s], [-s, s, -s], [ s, -s, -s]] faces = [[0, 2, 1], [0, 1, 3], [2, 3, 1], [3, 2, 0]] elif base == "octahedron": verts = [[ 0.0, 0.0, -1.0], [ 1.0, 0.0, 0.0], [ 0.0, -1.0, 0.0], [-1.0, 0.0, 0.0], [ 0.0, 1.0, 0.0], [ 0.0, 0.0, 1.0] ] faces = [[0, 1, 2], [0, 2, 3], [0, 3, 4], [0, 4, 1], [5, 2, 1], [5, 3, 2], [5, 4, 3], [5, 1, 4]] elif base == "icosahedron": t = (1 + np.sqrt(5))/2 tau = t/np.sqrt(1 + t*t) one = 1/np.sqrt(1 + t*t) verts = [[ tau, one, 0.0], [-tau, one, 0.0], [-tau, -one, 0.0], [ tau, -one, 0.0], [ one, 0.0, tau], [ one, 0.0, -tau], [-one, 0.0, -tau], [-one, 0.0, tau], [ 0.0, tau, one], [ 0.0, -tau, one], [ 0.0, -tau, -one], [ 0.0, tau, -one]] faces = [[4, 8, 7], [4, 7, 9], [5, 6, 11], [5, 10, 6], [0, 4, 3], [0, 3, 5], [2, 7, 1], [2, 1, 6], [8, 0, 11], [8, 11, 1], [9, 10, 3], [9, 2, 10], [8, 4, 0], [11, 0, 5], [4, 9, 3], [5, 3, 10], [7, 8, 1], [6, 1, 11], [7, 2, 9], [6, 10, 2]] else: raise ValueError("unknown base figure: " + str(base)) self.level = level self.base = base self.tree = [Level(faces=faces, vertices=[np.array(v) for v in verts])] for l in range(1, level + 1): self.tree.append(self.tree[-1].refine()) def get_approximation(self, n): return self.tree[n].faces, self.tree[n].vertices # import ipdb; ipdb.set_trace() arr = np.array sphere = Polyhedron(level=3, base="icosahedron") faces, vertices = sphere.get_approximation(3) x = arr([v[0] for v in vertices]) y = arr([v[1] for v in vertices]) z = arr([v[2] for v in vertices]) cell_scalar = arr([x[i] + x[j] + x[k] for i, j, k in faces]) triangle_faces_to_VTK("demo_sphere", x=x, y=y, z=z, faces=arr(faces), point_data={'s': x*x + y*y, 'n': (x, y, z), 'v': (-y, x, -np.sin(np.arccos(-1.0)*z))}, cell_data={'rho': cell_scalar}) pointsToVTK("demo_stock", x, y, z, {"rho": x*x + y*y})
Комментариев нет:
Отправить комментарий