BornAgain ≻ git-develop ≻ Documentation ≻ Reference ≻ Data types ≻ Class Datafield ≻ Fourier transform of 2D data

Fourier transform of 2D scattering data

The Fourier transform of an intensity map can expose real-space length scales that underlie a scattering pattern. An FFT requires values on an equidistant grid in the coordinate conjugate to the requested output. A spherical detector is equidistant in exit angles, not in scattering-vector coordinates. First create the equidistant $q_y,q_z$ approximation described in Axis coordinates and transformations:

transformation = ba.FrameTrafo.ScatteringToQ(wavelength, alpha_i)
q_result = transformation.transformedDatafield(angular_result)

Applying an FFT directly to the angular result is mathematically possible, but its conjugate axes are not real-space lengths, and the nonlinear angle-to-q relation distorts inferred periods. For real-space interpretation, use the reciprocal-space result. NumPy then provides the transform:

intensity = q_result.intensities()
fourier_values = np.fft.fft2(intensity)
fourier_magnitude = np.abs(np.fft.fftshift(fourier_values))

The conjugate coordinate values follow from the reciprocal-space bin widths. NumPy uses cycles in its frequency coordinates, whereas the scattering vector enters the Fourier phase as $q\cdot r$. Therefore the spatial coordinates contain a factor of $2\pi$:

q_y_step = np.diff(q_result.xCenters()).mean()
q_z_step = np.diff(q_result.yCenters()).mean()
n_z, n_y = intensity.shape
y = 2*np.pi*np.fft.fftshift(np.fft.fftfreq(n_y, d=q_y_step))
z = 2*np.pi*np.fft.fftshift(np.fft.fftfreq(n_z, d=q_z_step))

Attach these coordinates to the Fourier magnitude to retain Datafield plotting and cropping operations:

real_space_frame = ba.Frame(
    ba.ListScan("y (nm)", y.tolist()),
    ba.ListScan("z (nm)", z.tolist()))
fourier_result = ba.Datafield(
    real_space_frame, fourier_magnitude.ravel().tolist())

The complete example below simulates scattering from a square lattice, converts the angular detector axes to $q_y,q_z$, and plots the Fourier magnitude on real-space axes. The maxima repeat with the lattice period. The result is an autocorrelation-like map, not a direct reconstruction of the sample: scattering intensities contain no phase information, and the finite detector window can introduce artifacts.

Both heatmaps in the example use the default logarithmic color normalization. Changing to a linear normalization would only change the display; the plotted quantity would still be Fourier magnitude, not scattering-length density.

In the example, the complete fourier_result is kept; only a separate fourier_crop is plotted. Its horizontal range spans six lattice periods on each side of zero, and its more compact vertical range spans three. These limits show several repeats clearly in the stacked layout. They can be changed or the crop can be removed without recomputing the transform.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123

#!/usr/bin/env python3
# /// script
# requires-python = ">=3.10"
# dependencies = ["bornagain>=25,<26"]
# ///
"""
Fourier transform of a simulated 2D scattering pattern.
"""
import bornagain as ba
ba.require_versions("bornagain>=25,<26")
from bornagain import deg, nm
import numpy as np


wavelength = 0.04*nm
alpha = 0.2*deg
lattice_period = 10*nm


def get_sample():
    """
    Creates spheres on a square lattice.
    """
    # Materials
    red = (0.86, 0.24, 0.18)
    particle_material = ba.RefractiveMaterial("Particle", red, 6e-5, 2e-8)
    blue = (0.28, 0.57, 0.82)
    substrate = ba.RefractiveMaterial("Substrate", blue, 6e-6, 2e-8)

    # Particle arrangement
    particle = ba.Particle(particle_material, ba.Sphere(2.5*nm))
    lattice = ba.SquareLattice2D(lattice_period, 2*deg)
    layout = ba.Crystal2D(particle, lattice)
    layout.setDecayFunction(ba.Profile2DCauchy(50*nm, 50*nm, 0))

    # Layers
    vacuum_layer = ba.Layer(ba.Vacuum())
    vacuum_layer.deposit2D(layout)
    substrate_layer = ba.Layer(substrate)

    # Sample
    sample = ba.Sample()
    sample.addLayer(vacuum_layer)
    sample.addLayer(substrate_layer)
    return sample


def get_simulation(sample):
    """
    Creates a GISAS simulation with a two-dimensional detector.
    """
    beam = ba.Beam(1e9, wavelength, alpha)
    n = 200
    detector = ba.SphericalDetector(n, -1*deg, 1*deg, n, 0, 1*deg)
    return ba.ScatteringSimulation(beam, sample, detector)


def fourier_transform(q_result):
    """
    Transforms an equidistant q_y/q_z Datafield to real-space axes.
    """
    # Intensity and reciprocal-space bin widths
    intensity_values = q_result.intensities()
    q_y_step = np.diff(q_result.xCenters()).mean()
    q_z_step = np.diff(q_result.yCenters()).mean()

    # Fourier magnitude centered at zero frequency
    fourier_values = np.fft.fft2(intensity_values)
    centered_fourier_values = np.fft.fftshift(fourier_values)
    fourier_magnitude = np.abs(centered_fourier_values)

    # Conjugate real-space coordinates
    y_frequencies = np.fft.fftfreq(intensity_values.shape[1], d=q_y_step)
    z_frequencies = np.fft.fftfreq(intensity_values.shape[0], d=q_z_step)
    y_coordinates = 2*np.pi*np.fft.fftshift(y_frequencies)
    z_coordinates = 2*np.pi*np.fft.fftshift(z_frequencies)

    # Real-space Datafield
    real_space_frame = ba.Frame(
        ba.ListScan("y (nm)", y_coordinates.tolist()),
        ba.ListScan("z (nm)", z_coordinates.tolist()))
    flat_magnitude = fourier_magnitude.ravel().tolist()
    return ba.Datafield(real_space_frame, flat_magnitude)


if __name__ == '__main__':
    # Simulate intensity on angular detector axes
    sample = get_sample()
    angular_result = get_simulation(sample).simulate()

    # Convert the detector axes to reciprocal-space coordinates
    transformation = ba.FrameTrafo.ScatteringToQ(wavelength, alpha)
    q_result = transformation.transformedDatafield(angular_result)

    # Transform the intensity map to real-space coordinates
    fourier_result = fourier_transform(q_result)
    ba.showSample3D(sample, sample_size=100*nm, seed=0)

    # Crop only the displayed result
    y_min = -6*lattice_period
    y_max = 6*lattice_period
    z_min = -3*lattice_period
    z_max = 3*lattice_period
    fourier_crop = fourier_result.crop(y_min, z_min, y_max, z_max)

    # Plot reciprocal-space intensity and its Fourier magnitude
    figure = ba.plt.figure(figsize=(7, 7.5), layout="constrained")
    plot_axes = figure.subplots(2, 1)
    ba.plot_heatmap(
        q_result,
        ax=plot_axes[0],
        with_cb=True,
        unit_aspect=1,
        title="Scattering intensity",
        zlabel="Intensity")
    ba.plot_heatmap(
        fourier_crop,
        ax=plot_axes[1],
        with_cb=True,
        unit_aspect=1,
        title="Magnitude of Fourier transform",
        zlabel="Fourier magnitude (a.u.)")
    ba.plt.show()
auto/Examples/gisas/methods/FourierTransform.py