Это продолжение вопроса, на который я ответил, и который можно найти здесь .
У меня есть несколько точек (координаты x, y, z) в 3D-окне со связанными массами. Я хочу нарисовать гистограмму плотности массы, которая находится в сферах данного радиуса R. Идея состоит в том, чтобы вычислить трехмерную гистограмму моего блока (с биннингом, намного меньшим, чем радиус), взять его БПФ, умножить на фильтр (шар в реальном пространстве) и инвертировать БПФ результат. Оттуда я просто вычисляю 1D гистограмму значений, полученных в каждой 3D-корзине.
Следуя проблеме, возникшей у меня с использованием аналитического выражения фильтра в пространстве Фурье, я теперь генерирую шар в реальном пространстве и беру его БПФ, чтобы получить мой фильтр. Тем не менее, гистограмма, которую я получаю из этого метода, действительно странная, где я ожидал бы получить гауссовскую форму:
Мой код следующий:
import numpy as np
import matplotlib.pyplot as plt
import random
from numba import njit
# 1. Generate a bunch of points with masses from 1 to 3 separated by a radius of 1 cm
size = 100
radius = 1
rangeX = (0, size)
rangeY = (0, size)
rangeZ = (0, size)
rangem = (1,3)
qty = 300000 # or however many points you want
deltas = set()
for x in range(-radius, radius+1):
for y in range(-radius, radius+1):
for z in range(-radius, radius+1):
if x*x + y*y + z*z<= radius*radius:
deltas.add((x,y,z))
X = []
Y = []
Z = []
M = []
excluded = set()
for i in range(qty):
x = random.randrange(*rangeX)
y = random.randrange(*rangeY)
z = random.randrange(*rangeZ)
m = random.uniform(*rangem)
if (x,y,z) in excluded: continue
X.append(x)
Y.append(y)
Z.append(z)
M.append(1)
excluded.update((x+dx, y+dy, z+dz) for (dx,dy,dz) in deltas)
#print("There is ",len(X)," points in the box")
# Compute the 3D histogram
a = np.vstack((X, Y, Z)).T
b = 200
R = 10
H, edges = np.histogramdd(a, weights=M, bins = b)
Fh = np.fft.fftn(H, axes=(-3,-2, -1))
# Generate the filter in real space
Kreal = np.zeros((b,b,b))
X = edges[0]
Y = edges[1]
Z = edges[2]
mid = int(b/2)
s = (X.max()-X.min()+Y.max()-Y.min()+Z.max()-Z.min())/(3*b)
cst = 1/2 + (1/12 - (R/s)**2)*np.arctan((0.5*np.sqrt((R/s)**2-0.5))/(0.5-(R/s)**2)) + 1/3*np.sqrt((R/s)**2-0.5) + ((R/s)**2 - 1/12)*np.arctan(0.5/(np.sqrt((R/s)**2-0.5))) - 4/3*(R/s)**3*np.arctan(0.25/((R/s)*np.sqrt((R/s)**2-0.5)))
@njit(parallel=True)
def remp(Kreal):
for i in range(b):
for j in range(b):
for k in range(b):
a = cst - np.sqrt((X[i]-X[mid])**2 + (Y[j]-Y[mid])**2 + (Z[k]-Z[mid])**2)/s
if a >= 0.1 and a < 0.2:
Kreal[i][j][k] = 0.1
elif a >= 0.2 and a < 0.3:
Kreal[i][j][k] = 0.2
elif a >= 0.3 and a < 0.4:
Kreal[i][j][k] = 0.3
elif a >= 0.4 and a < 0.5:
Kreal[i][j][k] = 0.4
elif a >= 0.5 and a < 0.6:
Kreal[i][j][k] = 0.5
elif a >= 0.6 and a < 0.7:
Kreal[i][j][k] = 0.6
elif a >= 0.7 and a < 0.8:
Kreal[i][j][k] = 0.7
elif a >= 0.8 and a < 0.9:
Kreal[i][j][k] = 0.8
elif a >= 0.9 and a < 0.99:
Kreal[i][j][k] = 0.9
elif a >= 0.99:
Kreal[i][j][k] = 1
return Kreal
Kreal = remp(Kreal)
Kreal = np.fft.ifftshift(Kreal)
Kh = np.fft.fftn(Kreal, axes=(-3,-2, -1))
Gh = np.multiply(Fh, Kh)
Density = np.real(np.fft.ifftn(Gh,axes=(-3,-2, -1)))
# Generate the filter in fourier space using its analytic expression
kx = 2*np.pi*np.fft.fftfreq(len(edges[0][:-1]))*len(edges[0][:-1])/(np.amax(X)-np.amin(X))
ky = 2*np.pi*np.fft.fftfreq(len(edges[1][:-1]))*len(edges[1][:-1])/(np.amax(Y)-np.amin(Y))
kz = 2*np.pi*np.fft.fftfreq(len(edges[2][:-1]))*len(edges[2][:-1])/(np.amax(Z)-np.amin(Z))
kr = np.sqrt(kx[:,None,None]**2 + ky[None,:,None]**2 + kz[None,None,:]**2)
kr *= R
Kh = (np.sin(kr)-kr*np.cos(kr))*3/(kr)**3
Kh[0,0,0] = 1
Gh = np.multiply(Fh, Kh)
Density2 = np.real(np.fft.ifftn(Gh,axes=(-3,-2, -1)))
D = Density.flatten()
N = np.mean(D)
D2 = Density2.flatten()
N2 = np.mean(D2)
# I then compute the histogram I want
hist, bins = np.histogram(D/N, bins='auto', density=True)
bin_centers = (bins[1:]+bins[:-1])*0.5
plt.plot(bin_centers, hist,'.',label = "Defining the Filter in real space")
hist, bins = np.histogram(D2/N2, bins='auto', density=True)
bin_centers = (bins[1:]+bins[:-1])*0.5
plt.plot(bin_centers, hist,'.',label = "Using analytic expression")
plt.xlabel('Normalised Density')
plt.ylabel('Probability density')
plt.legend()
plt.show()
Вы понимаете, почему это происходит? Большое спасибо за вашу помощь.
PS: длинный список if операторов, когда я определяю Фильтр в реальном пространстве, основан на том, как я рисую сферу на сетке. Я присваиваю значение 1 всем ячейкам, которые на 100% находятся внутри сферы, а затем значение уменьшается по мере уменьшения объема, занимаемого сферой в ячейке. Я проверил, что это дает мне сферу нужного радиуса. Подробности по этому вопросу можно найти здесь (часть 2.5 и рисунок 8 для точности).
- EDIT -
Кажется, что код ведет себя так, только когда все массы частиц идентичны