El filtro en el espacio de Fourier no se comporta como se supone que debe

Este es un seguimiento de una pregunta respondida que hice y que se puede encontrar aquí .

Tengo varios puntos (coordenadas x, y, z) en un cuadro 3D con masas asociadas. Quiero dibujar un histogtwig de la densidad de masa que se encuentra en las esferas de un radio R dado. La idea es calcular un histogtwig 3D de mi caja (con un intervalo mucho más pequeño que el radio), tomar su FFT, multiplicar por el filtro (una bola en el espacio real) e invertir FFT el resultado. A partir de ahí, simplemente calculo el histogtwig 1D de los valores obtenidos en cada bandeja 3D.

Siguiendo el problema que tuve al usar una expresión analítica del filtro en el espacio de Fourier, ahora estoy generando la bola en el espacio real y tomo su FFT para obtener mi filtro. Sin embargo, el histogtwig que obtengo de este método es realmente extraño, donde esperaría que un gaussiano obtuviera esto: Histograma normalizado de la densidad de masa en mi caja

Mi código es el siguiente:

 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= 0.1 and a = 0.2 and a = 0.3 and a = 0.4 and a = 0.5 and a = 0.6 and a = 0.7 and a = 0.8 and a = 0.9 and 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() 

¿Entiendes por qué sucede esto? Muchas gracias por su ayuda.

PD: la larga lista de sentencias if cuando defino el filtro en el espacio real proviene de cómo estoy dibujando la esfera en la cuadrícula. Asigno el valor 1 a todos los contenedores que son 100% dentro de la esfera, y luego el valor disminuye a medida que el volumen ocupado por la esfera en el contenedor disminuye. Comprobé que me da una esfera del radio deseado. Los detalles sobre el tema se pueden encontrar aquí (parte 2.5 y figura 8 para mayor precisión).

–EDITAR–

El código solo parece comportarse así cuando todas las masas de partículas son idénticas

Mi problema viene de cómo estoy generando mi filtro. En mi código, la forma en que asocio el peso a los voxeles que no están completamente en la esfera es discontinua: por ejemplo, le doy el peso 0.1 a un voxel cuya relación de volumen está entre 0.1 y 0.2.

Por lo tanto, lo que sucede cuando todos los puntos tienen la misma masa es: tengo múltiplos de 1 en mi cuadrícula que multiplico con un número finito de coeficientes, por lo que hay un número finito de valores posibles que mi cuadrícula puede tomar, y por lo tanto, algunos contenedores son Vacío o al menos ‘menos lleno’. Esto es menos probable que ocurra cuando las masas de mi partícula se distribuyen de manera más continua.

Por lo tanto, una solución es asignar el peso correcto a los 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