5

i have the following problem: I want to integrate a 2D array, so basically reversing a gradient operator.

Assuming i have a very simple array as follows:

shape = (60, 60)
sampling = 1
k_mesh = np.meshgrid(np.fft.fftfreq(shape[0], sampling), np.fft.fftfreq(shape[1], sampling))

Then i construct my vectorfield as a complex-valued arreay (x-vector = real part, y-vector = imaginary part):

k = k_mesh[0] + 1j * k_mesh[1]

So the real part for example looks like this enter image description here

Now i take the gradient:

k_grad = np.gradient(k, sampling)

I then use Fourier transforms to reverse it, using the following function:

def freq_array(shape, sampling):

    f_freq_1d_y = np.fft.fftfreq(shape[0], sampling[0])
    f_freq_1d_x = np.fft.fftfreq(shape[1], sampling[1])
    f_freq_mesh = np.meshgrid(f_freq_1d_x, f_freq_1d_y)
    f_freq = np.hypot(f_freq_mesh[0], f_freq_mesh[1])

    return f_freq


def int_2d_fourier(arr, sampling):
    freqs = freq_array(arr.shape, sampling)

    k_sq = np.where(freqs != 0, freqs**2, 0.0001)
    k = np.meshgrid(np.fft.fftfreq(arr.shape[0], sampling), np.fft.fftfreq(arr.shape[1], sampling))

    v_int_x = np.real(np.fft.ifft2((np.fft.fft2(arr[1]) * k[0]) / (2*np.pi * 1j * k_sq)))
    v_int_y = np.real(np.fft.ifft2((np.fft.fft2(arr[0]) * k[0]) / (2*np.pi * 1j * k_sq)))

    v_int_fs = v_int_x + v_int_y
    return v_int_fs


k_int = int_2d_fourier(k, sampling)

Unfortunately, the result is not very accurate at the position where k has an abrupt change, as can be seen in the plot below, which displayes a horizontal line profile of k and k_int.

enter image description here

Any ideas how to improve the accuracy? Is there a way to make it exactly the same?

F. Win
  • 390
  • 5
  • 17
  • Well it crosses 0 at the same place right? which is the place of the transition from white to black –  Nov 27 '18 at 11:41
  • yes, which is good news, but at the maximum and minimum next to it, the deviation is significant. – F. Win Nov 27 '18 at 11:43
  • If your image is artificial make it bigger so you will be sampling more pixels, so you will have more frequencies. The problem comes from abrupt changes, which cannot be modelled unless you have infinite number of harmonics. –  Nov 27 '18 at 11:50
  • Of course i could resample this one, but for other problems my sampling is fixed. – F. Win Nov 27 '18 at 11:50
  • So your answer is that this is an artifact that cannot be resolved, as it stems from the discrete sampling of the data? – F. Win Nov 27 '18 at 11:51
  • Yup, with a finite number of harmonics you will never have a perfect transition. If you know how the transition should be then you can try to "learn" the correction, but in general doing a fft and then reversing it will always loose information. –  Nov 27 '18 at 11:54
  • Do you think one could do better with a real space integration? If yes, do you have an example code? – F. Win Nov 27 '18 at 11:56

1 Answers1

1

I actually found a solution. The integration itself yields very accurate results. However, the gradient function from numpy calculates second order accurate central differences, which means that the gradient itself already is an approximation.

When you replace the problem above with an analytical formula such as a 2D Gaussian, one can calculate the derivative analytically. When integrating this analytically derived function, the error is on the order of 10^-10 (depending on the width of the Gaussian, which can lead to aliasing effects).

So long story short: The integration function proposed above works as intended!

F. Win
  • 390
  • 5
  • 17