3D-контур на основе данных с использованием Mayavi/Python

Я хотел бы сделать трехмерный контурный график, используя Mayavi, точно так же, как третий рисунок на этой странице (модель электронного облака водорода):

http://www.sethanil.com/python-for-reseach/5

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

XYZV = [[1, 2, 3, 4],
        [6, 7, 8, 9],
        ...
        [4, 5, 6, 7]]

Точки данных не распределены равномерно в пространстве XYZ и не хранятся в каком-либо определенном порядке. Я думаю, что в примере используется сетка для создания точек данных - я просмотрел это, но совершенно не понимаю. Любая помощь приветствуется?

H
(источник: sethanil.com)


person joshlk    schedule 23.02.2012    source источник
comment
Покажите нам, что вы пробовали до сих пор. Мы будем помогать с этого момента.   -  person Don Question    schedule 23.02.2012
comment
Просто для справки в будущем, такие вопросы отлично подходят для вычислительных наук.   -  person David Z    schedule 24.02.2012


Ответы (2)


Хитрость заключается в интерполяции по сетке перед построением графика — я бы использовал для этого scipy. Ниже R находится массив (500,3) значений XYZ, а V — «величина» в каждой точке XYZ.

from scipy.interpolate import griddata
import numpy as np

# Create some test data, 3D gaussian, 200 points
dx, pts = 2, 100j

N = 500
R = np.random.random((N,3))*2*dx - dx
V = np.exp(-( (R**2).sum(axis=1)) )

# Create the grid to interpolate on
X,Y,Z = np.mgrid[-dx:dx:pts, -dx:dx:pts, -dx:dx:pts]

# Interpolate the data
F = griddata(R, V, (X,Y,Z))

Отсюда совсем несложно отобразить наши данные:

from mayavi.mlab import *
contour3d(F,contours=8,opacity=.2 )

Это дает красивую (неровную) гауссиану.

введите здесь описание изображения

Взгляните на документы для griddata, обратите внимание, что вы можете изменить метод интерполяции. Если у вас больше точек (как на интерполированной сетке, так и в наборе данных), интерполяция становится лучше и лучше представляет базовую функцию, которую вы пытаетесь проиллюстрировать. Вот приведенный выше пример с 10 000 точек и более мелкой сеткой:

введите здесь описание изображения

person Hooked    schedule 24.02.2012
comment
Большое спасибо. Работает как часы! Только один вопрос: если бы я хотел удвоить количество точек на аппроксимирующей сетке, что бы я изменил в строке «dx, pts = 2, 100j»? - person joshlk; 24.02.2012
comment
@ Джош, вы бы изменили его на dx, pts = 2, 200j, однако это удвоило бы количество точек в каждом измерении, поэтому у вас было бы в 2 ^ 3 = 8 раз больше точек для интерполяции. dx управляет размером сетки графика. Для более точного контроля достаточно mgrid для каждого линейного размера. - person Hooked; 24.02.2012
comment
Я считаю, что вычисление функции griddata может занять очень много времени. Знаете ли вы какие-нибудь советы по ускорению процесса? - person joshlk; 28.02.2012
comment
@ Джош, ты мог бы попробовать изменить интерполяцию на ближайшую, а не на линейную, я не пробовал тест времени, поэтому я не уверен, насколько это изменит. Кроме того, вы можете разделить график на секции: области с высокой плотностью могут иметь больше точек/вокселей, а области с низкой плотностью могут быть более разреженными. Однако для этого вам понадобятся отдельные команды griddata. - person Hooked; 28.02.2012

Вы можете использовать фильтр delaunay3d для создания ячеек из точек. Затем вы можете создать iso_surface() для вывода UnstructuredGrid delaunay3d. Если вам нужны ImageData, вы можете использовать фильтр image_data_probe.

import numpy as np
from tvtk.api import tvtk
from mayavi import mlab

points = np.random.normal(0, 1, (1000, 3))
ug = tvtk.UnstructuredGrid(points=points)
ug.point_data.scalars = np.sqrt(np.sum(points**2, axis=1))
ug.point_data.scalars.name = "value"
ds = mlab.pipeline.add_dataset(ug)
delaunay = mlab.pipeline.delaunay3d(ds)
iso = mlab.pipeline.iso_surface(delaunay)
iso.actor.property.opacity = 0.1
iso.contour.number_of_contours = 10
mlab.show()

введите здесь описание изображения

person HYRY    schedule 24.02.2012
comment
Привет, когда я использую этот метод со своими собственными данными, я продолжаю получать множество ошибок, появляющихся в новых окнах с последним высказыванием: Предупреждение: в C:\pisi\tmp\VTK-5.6.0-2\work\VTK\Graphics \vtkDelaunay3D.cxx, строка 488 vtkDelaunay3D (03AD2250): 27 вырожденных треугольников, подозрительное качество сетки. Вы знаете, что происходит? - person joshlk; 25.02.2012