Phase retrieval example

[ghost]

Hi, I would like to leverage Pylops for the domain of 'phase retrieval', but I don't know how to start on this. A typical phase retrieval problem has the following form:

|| y - | F^-1{ K * F{x}} | ||^2 + lam || |x| ||^1

y is a measured real-valued 2D image x is an unknown complex-valued 2D image F and F^-1 are the forward and backward 2D Fourier transforms K is a complex-valued 2D kernel (predefined dense matrix), defined in Fourier domain

So, intuitively, complex-valued 2D image x is 2D convoluted (by multiplication in Fourier domain) and only the magnitude is fitted to the known real-valued 2D image y. This problem is ill-posed, and needs regularization. I would e.g. like to add a sparsity prior on the magnitude of x and solve this using a Split-bregman solver or ISTA solver.

Would it be possible for you to provide a minimal example script that does just that, to point me in the right direction?

[mrava87]

Hi @tomhaeck0, interesting problem and definitely something pylops can solve :)

Here is a skeleton code to get you going (however I have a couple of questions/comments below):

import numpy as np
import matplotlib.pyplot as plt

from pylops.utils import dottest
from pylops.basicoperators   import *
from pylops.signalprocessing import *
from pylops.optimization.sparsity  import *

# Model
nx, ny = 20, 30
x = np.zeros(nx*ny) + 1j * np.zeros(nx*ny)
ipos = np.random.permutation(np.arange(nx*ny))[:100]
x[ipos] = np.random.normal(0, 1, 100) + 1j * np.random.normal(0, 1, 100)
x = x.reshape(nx, ny)

# Fourier operator
F = FFT2D(dims=(nx,ny), nffts=(nx,ny))

# Filtering operator
h = np.ones((2,2), dtype=np.complex)
K = Convolve2D(nx*ny, h, dims=(nx,ny), dtype=np.complex128)

# or Dense operator
#h = np.random.normal(0., 1., (nx*ny, nx*ny)) + 1j * np.random.normal(0., 1., (nx*ny, nx*ny))
#K = MatrixMult(h, dtype=np.complex128)

# Total operator
Op = F.H * K * F
Op = Op.toreal(forw=True, adj=False)

# Data
y = Op * x.ravel()
y  = y.reshape(nx, ny)

# Inverse
xinv = FISTA(Op, y.ravel(), eps=0.01, niter=1000, show=True)[0]
xinv = xinv.reshape(nx, ny)

plt.figure()
plt.imshow(np.abs(x), vmin=-1, vmax=1)

plt.figure()
plt.imshow(np.abs(y), vmin=-1, vmax=1)

plt.figure()
plt.imshow(np.abs(x)-np.abs(y), vmin=-1, vmax=1)

1. I wasn't sure how to interpret your definition of K, if its a dense matrix representing 2D FFT of a blurring filter in the original domain or if its a blurring filter in that domain (I gave you both options).
2. Regarding the solver I suggest FISTA. Split-Bregman makes only sense for more complex regularizations, say for example you want to minimize the TV(x) but for the L1 norm of x it is much faster to use FISTA

The x and y here are just random numbers, would be good to see how things work with your data, which are likely more meaningful :) Let us know how it goes!