傅立叶空间中的过滤器的行为不符合预期 [英] Filter in Fourier space does not behave like it's supposed to
问题描述
This is a follow-up to an answered question that I asked and that can be found here.
我在3D盒子中有几个点(x，y，z坐标)，并带有相关的质量.我想绘制在给定半径R的球体中发现的质量密度的直方图.想法是计算我的盒子的3D直方图(装仓比半径小得多)，进行FFT，乘以滤波器(真实空间中的球)，然后对FFT进行逆变换.从那里，我只计算每个3D仓中获得的值的1D直方图.
I have several points (x,y,z coordinates) in a 3D box with associated masses. I want to draw an histogram of the mass-density that is found in spheres of a given radius R. The idea is to compute a 3D histogram of my box (with binning much smaller than the radius), take its FFT, multiply by the filter (a ball in real space) and inverse FFT the result. From there, I just compute the 1D histogram of the values obtained in each 3D-bin.
在傅里叶空间中使用滤波器的解析表达式遇到了我遇到的问题之后，我现在在真实空间中生成球，并对其进行FFT以获得我的滤波器.但是，我从这种方法中得到的直方图确实很奇怪，在这里我期望得到一个高斯:
Following the issue I had by using an analytic expression of the filter in Fourier space, I am now generating the ball in real space and take its FFT to obtain my filter. However, the histogram I get out of this method is really strange, where I would expect a Gaussian I am getting this:
我的代码如下:
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()
您知道为什么会这样吗?非常感谢您的帮助.
Do you understand why this happens ? Thank you very much for your help.
PS:当我在实际空间中定义过滤器"时，if语句的长长列表来自于我在网格上绘制球体的方式.我将值1分配给球体中所有100％的分箱，然后随着分箱中球体所占体积的减小，该值减小.我检查它是否给了我想要的半径范围.可以在此处找到(第2.5部分和图8)以确保准确性).
PS: the long list of if statements when I define the Filter in real space comes from how i'm drawing the sphere on the grid. I assign the value 1 to all the bins that are 100% within the sphere, and then the value decreases as the volume occupied by the sphere in the bin decreases. I checked that it gives me a sphere of the radius wanted. Details on the subject can be found here(part 2.5 and figure 8 for accuracy).
-编辑-
仅当所有粒子质量相同时，代码才看起来像这样
The code only seems to behave like this when all the particle masses are identical
推荐答案
我的问题来自我如何生成过滤器.在我的代码中，将权重不完全与球体中的体素相关联的方式是不连续的:例如，将权重0.1赋予体积比在0.1到0.2之间的体素.
My problem comes from how I am generating my filter. In my code, the way I associate weight to voxels not entirely in the sphere is discontinuous: For example I give the weight 0.1 to a voxel whose volume ratio is between 0.1 et 0.2.
因此，当所有点都具有相同质量时会发生以下情况:我在网格中乘以1的倍数，然后乘以有限数量的系数，因此网格可以取有限数量的可能值，因此一些垃圾箱是空的，或者至少是装满"的.当我的粒子质量更连续地分布时，这种情况不太可能发生.
Thus what happens when all points have the same mass is: I have multiples of 1 in my grid that I multiply with a finite number of coefficients, thus there is a finite nummber of possible values that my grid can take, and thus some bins are empty or at least 'less full'. This is less likely to happen when the masses of my particle are more continuously distributed.
因此，一种解决方法是为体素指定合适的权重.
A fix is thus to appoint the right weight to the voxels.
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.99:
Kreal[i][j][k] = a
elif a >= 0.99:
Kreal[i][j][k] = 1
这篇关于傅立叶空间中的过滤器的行为不符合预期的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
