3d-интерполяция в Scipy - сетка плотности - PullRequest
Новая
       33

3d-интерполяция в Scipy - сетка плотности

4 голосов
/ 15 июня 2011

У меня есть набор координат, представляющих 3d-позиции ряда объектов. Эти точки взяты из симуляции трехмерного куба. Мне нужно сделать трехмерную сетку на кубе, а затем назначить эти точки в их правильных местах на сетке, чтобы я мог затем найти плотность объектов в каждом разделе сетки. Я искал документацию по интерполяции и сетке (например, http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html#scipy.interpolate.griddata),), но не уверен, что с этим делать, потому что у меня нет функции, связанной с этими точками данных. Я думал о том, чтобы использовать оператор if : дан массив точек = [x1, y1, z1], [x2, y2, z2] и т. д., если points [i] [0] -gridpoint1 [0] <1: if points [i] [1] - gridpoint1 [1] <1: если points [i] [2] -gridpoint1 [2] <1, points [i] = bin1 [i], где bin1 - это готовый массив нулей. Однако я думаю, что Я должен выполнить это для каждой точки сетки на сетке (точка сетки будет в центре каждой ячейки), а затем выяснить, сколько ненулевых элементов было в каждой ячейке, что я также не уверен, как это сделать. У меня есть чувство, что я могу использовать какую-то функцию в scipy, чтобы выполнить эту задачу легче, но я все еще не уверен, как туда добраться. Большое спасибо заранее за вашу помощь! </p>

1 Ответ

4 голосов
/ 16 июня 2011

Если я правильно понимаю, у вас есть координаты (x [i], y [i], z [i]), i = 0, ..., N-1, и вы хотите посчитать, сколько из них заканчивается в заданной ячейке сетки в трехмерном кубе?

Это можно сделать с помощью numpy.histogramdd: 

import numpy as np

# some random points in a cube
x = np.random.rand(100)
y = np.random.rand(100)
z = np.random.rand(100);

# specify grid in a cube
xgrid = np.linspace(0.0, 1.0, 5)
ygrid = np.linspace(0.0, 1.0, 6)
zgrid = np.linspace(0.0, 1.0, 7)

# compute counts
counts, edges = np.histogramdd(np.c_[x, y, z], bins=(xgrid, ygrid, zgrid))

print counts.shape # -> (4, 5, 6)

# number of points in the box 
# xgrid[1] <= x < xgrid[2] && ygrid[2] <= y < ygrid[3] && zgrid[0] <= z < zgrid[1]
print counts[1, 2, 0]

Если вы хотите выяснить, в какой ячейке сетки находится каждая из точек, это можно сделать с помощью searchsorted: 

ix = np.searchsorted(xgrid, x) - 1
iy = np.searchsorted(ygrid, y) - 1
iz = np.searchsorted(zgrid, z) - 1

# point (x[3], y[3], z[3]) is in the following grid cell
print ix[3], iy[3], iz[3]
...