Reflected light

[golovart]

How does one simulate the reflected light (light scattered backwards)? I guess, for the cases when scatterer is out of focal plane, if it's on one side of focal plane - the output is forward scattered light, opposite side would give backward scattered, right? But if you want to use Lens theory to simulate a scatterer in the focal plane, what modifications could be done to obtain the reflected light?

[briandleahy]

@golovart I'm not quite sure I understand your question.

The Lens and MieLens theories work for scatterers above, in, and below the focal plane of a microscope objective lens. So if you want to simulate a hologram of a particle in the focus of an objective lens, then you need no modifications -- just make the object's z-position 0. That hologram is still calculated using the forward scattered light, because that is what physically enters the objective lens.

Or are you asking how to calculate backscattered light, when the imaging setup is not an in-line holography configuration?

It might help if you uploaded or linked a schematic of the imaging set up you're interested in.

[golovart]

@briandleahy in the microscope setup I am interested in, we are looking at reflected light (not for hologram, just optical image)

The schematic microscope configuration is the following:

Therefore, my goal is to obtain backscattered light and propagate it through the lens as usual.

I was thinking about simply changing theta to pi - theta only in the part that is passed to the scattering matrices calculation, like only in lens.py:

    def _calc_scattering_matrix(self, scatterer, medium_wavevec, medium_index):
        if not self.refl_light:
            theta, phi = np.meshgrid(self._theta_pts, self._phi_pts)
        else:
            theta, phi = np.meshgrid(np.pi - self._theta_pts, self._phi_pts)
        pts = detector_points(theta=theta.ravel(), phi=phi.ravel())
        illum_wavelen = 2 * np.pi * medium_index / medium_wavevec

        pts = update_metadata(pts, medium_index=medium_index,
                              illum_wavelen=illum_wavelen)
        S = self.theory.calculate_scattering_matrix(scatterer, pts)
        S = np.conj(S.values.reshape(self.quad_npts_theta,
                                     self.quad_npts_phi, 2, 2))
        S = np.swapaxes(S, 0, 1)
        S1 = S[:, :, 1, 1].reshape(self.quad_npts_theta, self.quad_npts_phi, 1)
        S2 = S[:, :, 0, 0].reshape(self.quad_npts_theta, self.quad_npts_phi, 1)
        S3 = S[:, :, 0, 1].reshape(self.quad_npts_theta, self.quad_npts_phi, 1)
        S4 = S[:, :, 1, 0].reshape(self.quad_npts_theta, self.quad_npts_phi, 1)
        return S1, S2, S3, S4

But it leads to weird and unintuitive results for the output image.

[briandleahy]

Hi @golovart ,

As you guessed, holopy does not currently support the calculation of holograms with backscattered light. However, it should be possible to implement this with a little bit of effort.

I don't see any problems with what you wrote for calculating the scattered beam. (I would do it through composition rather than changing the Lens class itself, but the physics looks OK. If you want, I can append a code snippet.)

There are also some subtleties in how the Lens theory treats the scattered beam. First, the Lens theory treats the scattered beam as undergoing 1 Gouy phase shift as it propagates through the optics. I think this should be the same for your setup, but you should think about it. Second, depending on where your scatterer is relative to the reflective sample plane, the phase difference between the incident and scattered field may change with z, as compared to the case for an inline holographic microscope. For instance, if your particle is a distance z above the reflecting surface, then the unscattered beam will travel an additional distance of 2z and accumulate an additional phase of exp(i*k*2z) as compared to the in-line case. Third, for your scattering geometry, I would think about whether the phase of the incident beam changes after it reflects off the sample.

What are the unintuitive results for the output image?

[golovart]

Here is a code-example for what I mean (with that redefinition in Lens class):

s1 = Sphere(r = 0.05, center = (0, -0.3, 0))
s2 = Sphere(r = 0.05, center = (0, 0.3, 0))
detector = hp.detector_grid(100, .01)
bisph = Scatterer(lambda point: np.logical_or(s1.contains(point), s2.contains(point)), n = 0.05+1j*2.1, center = (0.5, 0.5, 0))
len_dda = Lens(1.38, DDA(use_gpu=True))

# changing bounds by hand, since there is nothing in the center
# and default algorithm will fail
bond = np.array(bisph.bounds)
c = bisph.center
for i,b in enumerate(bond):
    bond[i,0] = min(s1.bounds[i][0], s2.bounds[i][0])
    bond[i,1] = max(s1.bounds[i][1], s2.bounds[i][1])
bisph.indicators.bound = bond

fig, axes = plt.subplots(1, 2, figsize=(15, 6))
axes = np.ravel(axes)
for refl in [False, True]:
    len_dda = Lens(1.38, DDA(use_gpu=True), refl_light=refl)
    holo = calc_field(detector, bisph, medium_index=1., illum_wavelen=0.4, illum_polarization=(1,0), theory=len_dda)
    holo = (np.abs(holo)**2).sum(dim='vector')
    holo.plot(ax=axes[int(refl)])
    axes[int(refl)].set_title('refl' if refl else 'scat')
plt.show()

And the output is:

We are experimentally looking at intensities on the camera pixels, not holograms. The problem with reflected light is that we expect two resolved objects (small spheres and distance is bigger than wavelength), as on the forward scattered image, but we get intensity minima at the spheres' locations and peak in the center...

[briandleahy]

Hi @golovart ,

Thanks for the image. That does look weird.....

One subtlety with this is that it's not obvious to me if the problem is with the code or with my expectations for the physics. I don't really have a good intuition for the image of two strong scatterers separated about 1 wavelength apart, as imaged in reflection mode through a high-aperture lens, so it's not obvious to me that this is a problem with the implementation. (This would not be the first time that the effect of a lens on the imaging has surprised me....)

Usually when I see something like this, I try to find an edge case where I am certain of the physics. For example, I might check how the hologram changes as I move the spheres farther and farther apart. Or I might check if this persists for a single particle rather than two.

How sure are you that the image on the right is wrong here? Could it be possible that there is a scattering resonance that gets propagated back through the lens differently than intuition might suggest?

[golovart]

@briandleahy I checked for the far apart (3 wavelength) spheres, and we get two clear, almost independent images:

And this pattern remains for a single sphere as well:

I can surely say that experimentally in the microscope (of a scheme above) small objects are seen as "blobs", not as circles. We see something similar to the "scattered" image from this simulation.

Could it be due to a phase change of pi from the reflection? And where should that be fixed? At the moment when the fields are calculated from scattering matrices (and before passing to lens)?

[briandleahy]

Hi @golovart,

I think I've figured this out. It is an issue with correctly modeling the optics. It took me a while, so humor me as I walk through my reasoning.

First, I have intuition for Rayleigh scatterers, so I took a look at a Rayleigh scatterer in forward and backscattered light. Rather than changing the lens code, I did this via composition, by defining a new Backscattered theory:

# in a separate file, backscatter.py
import numpy as np

from holopy.scattering.theory.scatteringtheory import ScatteringTheory

class Backscattered(ScatteringTheory):
    def __init__(self, theory):
        self.theory = theory

    def _can_handle(self, scatterer):
        return self.theory._can_handle(scatterer)

    def _raw_scat_matrs(self, scatterer, pos, medium_wavevec, medium_index):
        r, theta, phi = pos
        theta_backscattered = np.pi - theta
        pos_backscattered = np.array([r, theta_backscattered, phi])

        args = (scatterer, pos_backscattered, medium_wavevec, medium_index)
        return self.theory._raw_scat_matrs(*args)

Then, I generate a hologram using the backscattered theory and the Mie calculation instead of DDA (to omit one more possible source of problems):

import holopy as hp
from holopy.scattering.interface import calc_intensity
from holopy.scattering import Sphere
from holopy.scattering.theory import Lens, Mie

import matplotlib.pyplot as plt; plt.close('all')

from backscatter import Backscattered

THEORIES = {
    'normal': Lens(1.38,  Mie()),
    'backscattered': Lens(1.38, Backscattered(Mie()))}

def get_rayleigh_intensity(theory):
    scatterer = Sphere(r=0.05, center=(0.5, 0.5, 0), n=1.01)
    detector = hp.detector_grid(100, .01)
    kwargs = {
        'medium_index': 1.0, 'illum_wavelen': 0.4, 'illum_polarization': (1, 0)}
    return calc_intensity(detector, scatterer, theory=theory, **kwargs)

def make_unscaled_comparison():
    f = plt.figure()
    for i, name in enumerate(THEORIES):
        intensity = get_rayleigh_intensity(THEORIES[name])
        ax = f.add_subplot(1, 2, i + 1)
        ax.imshow(
            intensity.values.squeeze(),
            interpolation='bicubic',
            vmin=0,
            cmap='inferno')
        ax.set_title("Rayleigh scatterer,\nin focus, {}".format(name))

This gives the following image, which is pretty similar to the one you show above:

However the scale in the two images is very different. Plotting on the same scale gives something like this:

def make_scaled_comparison():
    f = plt.figure()
    for i, name in enumerate(THEORIES):
        intensity = get_rayleigh_intensity(THEORIES[name])
        ax = f.add_subplot(1, 2, i + 1)
        ax.imshow(
            intensity.values.squeeze(),
            interpolation='bicubic',
            vmin=0,
            vmax=2.44e-6,
            cmap='inferno')
        ax.set_title("Rayleigh scatterer,\nin focus, {}".format(name))

The backscattered image is barely visible. Here is why:

For a spherically symmetric scatterer imaged through a lens (regardless of backscattered light), the scattered fields can be written in terms of two integrals, I_0 and I_2 (see this paper: https://aip.scitation.org/doi/full/10.1063/5.0015976). The I_0 term corresponds to the rotationally symmetric portion of the field in the detector plane; the integrand of I_0 depends on S_perp(theta) + S_parallel(theta). The I_2 term corresponds to the portion of the field with a cos(2phi) variation in the detector plane; that integrand depends on S_perp(theta) - S_parallel(theta). Other terms in the integrand tend to make the I_2 term's magnitude much smaller than the I_0 term's.

For a Rayleigh scatterer, S_parallel ~ cos(theta), and S_perp ~ 1. Thus, for forward scattered light, cos(theta) > 0, and S_parallel and S_perp add in the I_0 term, giving the typical bright, mostly symmetric spot that you see on the left. But for backscattered light, cos(theta) < 0, and S_parallell + S_perp is basically zero. This gives the drastically reduced intensity that you see in my second figure, on the left. Finally, with the I_0 term drastically reduced for backscattered light, the I_2 term is the dominant one. Since the I_2 term has a cos(2phi) structure, that gives the funny wings that are visible in the first figure on the right, and in your figures. In other words, the problem is that the sign of the scattered light changed when it really shouldn't have.

The reason is that the transformation theta -> pi - theta needs to happen in two places: the calculation of the scattering matrices and the decomposition of the incoming light into theta and phi components. This gives an additional factor of -1 to two of the scattering matrices. If I add that in to the Backscattered theory I defined above, like so:

class Backscattered(ScatteringTheory):
    def __init__(self, theory):
        self.theory = theory

    def _can_handle(self, scatterer):
        return self.theory._can_handle(scatterer)

    def _raw_scat_matrs(self, scatterer, pos, medium_wavevec, medium_index):
        r, theta, phi = pos
        theta_backscattered = np.pi - theta
        pos_backscattered = np.array([r, theta_backscattered, phi])

        args = (scatterer, pos_backscattered, medium_wavevec, medium_index)
        scat_matrices = np.array(self.theory._raw_scat_matrs(*args))

        # account for coordinate change in decomposition of incoming
        # light into theta, phi components:
        scat_matrices[:, 0, 0] *= -1
        scat_matrices[:, 1, 0] *= -1  # FIXME is it this or the transpose?

        return scat_matrices

and I re-run the make_scaled_comparison() above, I get this, which looks great: Note that the backscattered image is ever so slightly fainter than the forward scattered, presumably due to the scatterer not being a perfect Rayleigh scatterer.

I think that solves the problem you were having. However, you should think about this more, since it is possible I made a physics mistake here.

Let me know if that helps, but if it doesn't, I think you're on your own :smile:

[golovart]

@briandleahy thank you for your help and explanation!

One last thing about scattering matrix vs transposed. In Bohren&Huffman (and ADDA manual) the scattering matrix is defined as:

So I would expect elements S[0,0] and S[0,1] to correspond to S-parallel.

While on line 233 in dda.py

scat_matr = np.array([[s[:,1], s[:,2]], [s[:,3], s[:,0]]]).transpose()

I looks like the definition from the formula, but then transposed. Why? Would that mean now S[0,0] and S[1,0] describe S-parallel?

[briandleahy]

Hi @golovart ,

Glad to hear that my explanation helped.

Regarding the scattering matrices, I don't know recall which definition holopy uses for the off-diagonal scattering matrices. Are you asking as part of this reflected light question though, or as a separate question? if it's a separate question, can I ask you to close this issue and open a new one?

[golovart]

I got a bit confused in your code with the line

scat_matrices[:, 1, 0] *= -1 # FIXME is it this or the transpose?

So I tried to understand from dda.py code is this off-diagonal the one we need and got the question above

[briandleahy]

@golovart I think that you're right; it should be the other way.

Closing this issue then.