Phase Reconstruction via Gerchberg–Saxton Algorithm

Light propagates as a wave and waves have not only amplitude, but also phase.  We typically aren't aware of this as the majority of the light we see is highly incoherent and spectrally broad.  However, if we have a scene illuminated with a laser, the speckle we observe will be a direct result of the coherent nature of the light and the phase information it contains.  I wanted to know what that phase was.

Suppose that we have light propagating through a series of planes and there are no scatterers in between them.  If we know the complex field (ie amplitude and phase) on one plane, we can calculate the complex field on another plane through a variety of mathematical techniques.  My favorite uses FFTs to perform a plane wave expansion.  However, in practice, we don't know the complex field on a plane.  A standard camera records only the power of the wave which is related to abs(E).  Put more simply, when you snap an image of a plane, you get one number per pixel, but the complex field requires two.

However, suppose we measure the power at two planes.  This can be achieved simply by defocusing the camera and snapping a second "blurry" image of the scene.  That "blurry" image contains the exact same information as the focused image because the one is a simple mathematical transformation of the other.  In some sense, you have just measured the same thing twice and therefore have two measurements per pixel.

The challenge is to then combine the images take at two different planes into a single complex image on a single plane.  In theory, you have all the information you need, however, the quantities you know are abs(E(x,y,z)) on z = z1 and z = z2.  There is no direct algebraic way of solving this.  One way to get there numerically, however, is to employ what is known as the Gershberg-Saxton Method.  

The algorithm is pretty simple.  You begin on "plane-A" and make the gross assumption that the phase is uniform.  You then mathematically propagate that field to "plane-B".  You have measured data for the magnitude on "plane-B" and so it will be obvious that the magnitude of the field you just calculated doesn't match well with the one you measured.  No matter.  You keep the phase from the field you just calculated and take the magnitude of the field you know to be there.  You mathematically propagate this field back to "plane-A".  Again, you keep the phase, but substitue in the experimental magnitude.  You continue to bounce the wave back and forth between the two planes and with each pass the difference between the calculated magnitude and the known magnitude will become less and less.

There are variations one could do on this.  We could use images from many planes.  

In order to verify this actually works, I tested it out on some manufactured data from COMSOL.  You can see the starting data in Fig 1 and the results of the code in Fig 2 and Fig 3.