The problem of retrieving a signal given the magnitude of its Furier transform is ubiquitous across all of applied sciences, due to the interpretation of such transform as the near field (Fresnel) diffraction. In essence, phase retrieval corresponds to finding out what object in space casts a given shadow.

In the 1950s, the shape of DNA was established using what is now called X-ray crystallography and variations of this technique are still widely employed in biophysical chemistry to study the shape and behaviour of different molecules. This approach to understanding the word at the nano scale culminated with coherent diffraction imaging, and transmission electron microscopy, both ideas based on phase retrieval. Not just our understanding of the very small has benefited from the efficient solving of this problem, but also the understanding of the vast space. phase retrieval being the key ingredient in correcting the defects of the Hubbel Space Telescope.

The problem of phase retrieval consists of finding an $x\in {ℂ}^n$ such that $$|{(ℱx)}_i|=b_i\forall i\in 1..n,$$ where $b\in {ℝ}^n$ is the given data.

Clearly, in this general form the problem is ill posed because there are infinitely many solutions (in fact the solution correspond to ${𝕋}^n$ ), so further constraints have to be imposed. Standard across the literature are either real ( $x\in {ℝ}^n$ ), positivity ( $x_i\ge 0$ , $\forall i\in 1..n$ ), bounded support ( $x_i=0$ , $\forall u\notin m_1\mathrm{\dots }{m}_2$ ), and sparsity ( ${\Vert x\Vert }_0\le s$ ) or a combination of the above.

In this work, we will concentrate on the real and positivity constraints, because they are illustrative of what happens with the broader class. As a remark, the real constraint is a consequence of a certain symmetry of the given data, so imposing this can yield an inconsistent problem. In this setting, the problem can be relaxed to $$\underset{x\in {ℝ}_{+}^n}{\min }\sum _{i=1}^n{\left|{|ℱx_i|}^2-b_i^2\right|}^2,$$ or $$\underset{x\in {ℝ}_{+}^n}{\min }\sum _{i=1}^n{\left||ℱx_i|-b_i\right|}^2.$$

The advantage of the first formulation is that it avoids the non smoothness of the absolute value, while the second formulation does not suffer from the slow decrease of a quartic.

Now that we are more familiar with the problem we can ask how well posed it is. For this, a simple numerical experiment may provide a crucial insight.

First, let's define a utility function that provides the phase and amplitude of the Furier transform of a given image.

using FFTW
function decomposed_furier(x)
    ℱx = fft(x)
    ϕ_ℱx = ℱx ./ abs.(ℱx)
    abs_ℱx = ℱx ./ ϕ_ℱx
    ϕ_ℱx, abs_ℱx
end

The inverse function will also provide a certain utility.

function recomposed_furier(ϕ_ℱx, abx_ℱx)
    ifft(ϕ_ℱx .* abx_ℱx)
end

We can now read two test images and decompose them into the phase and amplitude of their respective Furier transforms.

using Images
using TestImages

image₁ = Float32.(Gray.(testimage("barbara_gray_512.bmp")))
image₂ = Float32.(Gray.(testimage("cameraman.tif")))

Gray scale images are simpler to handle, but this experiment maps nicely to full color, by doing it once per RGB channel.

Now, we simply decompose both images, exchange the phase information between the two, and construct two new images.

ϕ_image₁, abs_image₁ = decomposed_furier(image₁)
ϕ_image₂, abs_image₂ = decomposed_furier(image₂)

new_image₁ = recomposed_furier(ϕ_image₁, abs_image₂)
new_image₂ = recomposed_furier(ϕ_image₂, abs_image₁)

using Plots
p₁ = plot(Gray.(image₁))
p₂ = plot(Gray.(image₂))
new_p₁ = plot(Gray.(new_image₁))
new_p₂ = plot(Gray.(new_image₂))

p = plot(p₁, p₂, new_p₁, new_p₂, layout=(2, 2))
savefig("results.png")

This makes it painfully clear that most of the information is stored in the phase of an object, since even random amplitudes can generate a visually recognizable image.

Soon to come