%%html
<style>
.container-fluid.main-container {
margin-left: 3em;
margin-right: 8em;
max-width: 100vw;
}
</style>
Introduction¶
This notebook showcases different functionalities of the GVS propagator, which is an efficient algorithm for simulating partially coherent light. The model uses the Van Cittert-Zernike theorem to find the complex coherence factor from a light source reaching the object plane, and it then uses the generalized Schell’s theorem to compute the diffraction pattern produced at a given propagation distance, not limited to the far field.
Eq. (1) gives the numerical propagation of GVS, with $\mu_{\mathcal{F}}$ given in Eq. (2). GVS also supports the use of any existing coherent propagator, as shown in Eq. (3), where $I_{\textrm{coherent}}$ is calculated with the object multiplied by a quadratic phase factor, $P(\xi, \eta) e^{i\frac{k}{2 z_{so}}(\xi^2 + \eta ^2)}$.
$$ I \left( \frac{x}{\lambda z_{oc}}, \frac{y}{\lambda z_{oc}} \right) = \frac{I_0}{(\lambda z_{oc})^2} \pmb{F} \left\{ \pmb{F^{-1}} \left\{ \left\| \pmb{F} \left\{ P(\xi, \eta) e^{i\frac{k}{2} \left ( \frac{1}{z_{so}} + \frac{1}{z_{oc}} \right ) (\xi^2 + \eta ^2)} \right\} \right\| ^ 2 \right\} \mu_{\mathcal{F}}(\Delta \xi, \Delta \eta) \right\} \qquad (1) $$$$ \mu_{\mathcal{F}} \left ( \frac{\Delta \xi}{\lambda z_{so}}, \frac{\Delta \eta}{\lambda z_{so}} \right ) = \frac{\pmb{F} \{I(u, v) \}}{\|\pmb{F} \{I(u, v) \}\|_\textrm{MAX}} \qquad (2) $$$$ I (x, y) = \frac{I_0}{(\lambda z_{oc})^2} \pmb{F} \left\{ \pmb{F^{-1}} \left\{ I_\textrm{coherent} \right\} \mu_{\mathcal{F}}(\Delta \xi, \Delta \eta) \right\} \qquad (3) $$import base64
import numpy as np
import pandas as pd
import plotly
import plotly.express as px
import imageio.v3 as imageio
from scipy.ndimage import zoom
from IPython.display import display, HTML, Image, Markdown
from gvs_propagator.constants import nm, um, mm, cm
from gvs_propagator.coherent import compute_angular_spectrum
from gvs_propagator.partially_coherent import (
compute_complex_coherence_factor,
compute_schells_theorem,
compute_gvs,
compute_angular_spectrum_gvs,
compute_fresnel_convolution_gvs,
compute_shifted_fresnel_gvs,
compute_mutual_intensity_integral_cross_section,
compute_intensity_point_sources_dask,
compute_intensity_cpr,
)
from gvs_propagator.utils import (
normalize,
crop_center,
mean_squared_error,
zoom_matrix
)
from gvs_propagator.notebooks.utils import display_image
from gvs_propagator.masks import (
create_circular_mask,
create_double_pinhole_mask,
)
plotly.offline.init_notebook_mode()
def add_margin(fig, margin=20, **kwargs):
fig.update_layout(margin_b=margin)
return Image(fig.to_image(**kwargs))
GVS Diffraction vs Far Field Schell Diffraction¶
def show_diffraction():
dx_object = (wavelength * z_object_camera) / (2 * dx_camera * num_pixels_x)
dx_source = (wavelength * z_source_object) / (2 * dx_object * num_pixels_x)
dy_object = (wavelength * z_object_camera) / (2 * dy_camera * num_pixels_y)
dy_source = (wavelength * z_source_object) / (2 * dy_object * num_pixels_y)
source_intensity = create_circular_mask(num_pixels_x * 2, dx_source, source_radius, num_pixels_y * 2, dy_source)
complex_coherence_factor = compute_complex_coherence_factor(source_intensity)
object = create_double_pinhole_mask(
num_pixels_x,
dx_object,
radius=pinhole_radius,
separation=pinhole_separation,
num_pixels_y=num_pixels_y,
dy=dy_object,
)
intensity_gvs = crop_center(
compute_gvs(
object,
complex_coherence_factor,
z_source_object,
z_object_camera,
wavelength,
dx_object
),
)
intensity_far_field = crop_center(
compute_schells_theorem(object, complex_coherence_factor, z_object_camera, wavelength)
)
x = np.linspace(-num_pixels_x // 2, num_pixels_x // 2, num_pixels_x) * dx_camera
y = np.linspace(num_pixels_y // 2, -num_pixels_y // 2, num_pixels_y) * dy_camera
return add_margin(display_image(
field=[intensity_gvs, intensity_far_field],
x=x,
y=y,
show_ticks=True,
show_as_grid=True,
normalize_all_frames=True,
grid_titles=["GVS", "Far Field Schell"],
title=f"Propagation distance: {int(z_object_camera * 100)} cm",
return_figure=True
), width=512, height=300)
Source Radius: $100 \ \mu m$¶
num_pixels_x = num_pixels_y = 1024
dx_camera = dy_camera = 5 * um
wavelength = 633 * nm
z_object_camera = 10 * cm
z_source_object = 100 * cm
source_radius = 100 * um
pinhole_radius = 100 * um
pinhole_separation = 600 * um
for z_object_camera in np.array([10, 20, 50]) * cm:
display(show_diffraction())
Source Radius: $1 \ mm$¶
source_radius = 1 * mm
for z_object_camera in np.array([10, 20, 50]) * cm:
display(show_diffraction())
Source Radius: $4 \ mm$¶
source_radius = 4 * mm
for z_object_camera in np.array([10, 20, 50]) * cm:
display(show_diffraction())
Model Comparison - Double pinhole diffraction¶
def show_diffraction(integral_1d=False):
dx_object = (wavelength * z_object_camera) / (2 * dx_camera * num_pixels_x)
dx_source = (wavelength * z_source_object) / (2 * dx_object * num_pixels_x)
dy_object = (wavelength * z_object_camera) / (2 * dy_camera * num_pixels_y)
dy_source = (wavelength * z_source_object) / (2 * dy_object * num_pixels_y)
dx_source_alternative = (wavelength * z_source_object) / (2 * dx_camera * num_pixels_x)
dy_source_alternative = (wavelength * z_source_object) / (2 * dy_camera * num_pixels_y)
object = create_double_pinhole_mask(
num_pixels_x,
dx_object,
radius=pinhole_radius,
separation=pinhole_separation,
num_pixels_y=num_pixels_y,
dy=dy_object,
)
object_approximation = create_double_pinhole_mask(
num_pixels_x,
dx_camera,
radius=pinhole_radius,
separation=pinhole_separation,
num_pixels_y=num_pixels_y,
dy=dy_camera,
)
intensities = []
radii = []
names = ["Integral Computation 1", "Integral Computation 2"] * len(source_radii)
if integral_1d:
num_pixels_y_integral = 1
else:
num_pixels_y_integral = num_pixels_y
for source_radius in source_radii:
source_intensity = create_circular_mask(num_pixels_x * 2, dx_source, source_radius, num_pixels_y * 2, dy_source)
complex_coherence_factor = compute_complex_coherence_factor(source_intensity)
intensities.append(normalize(
compute_mutual_intensity_integral_cross_section(
crop_center(object, num_pixels_y_integral, num_pixels_x),
np.fft.ifftshift(crop_center(complex_coherence_factor, num_pixels_y_integral, num_pixels_x)),
z_source_object,
z_object_camera,
wavelength,
dx_camera,
dx_object,
chunk_size=300
)
))
radii.append(f"{round(source_radius * 1e6)} µm")
source_intensity = create_circular_mask(
num_pixels_x * 2, dx_source_alternative, source_radius, num_pixels_y * 2, dy_source_alternative
)
complex_coherence_factor = compute_complex_coherence_factor(source_intensity)
intensities.append(normalize(
compute_mutual_intensity_integral_cross_section(
crop_center(object_approximation, num_pixels_y_integral, num_pixels_x),
np.fft.ifftshift(crop_center(complex_coherence_factor, num_pixels_y_integral, num_pixels_x)),
z_source_object,
z_object_camera,
wavelength,
dx_camera,
chunk_size=300
)
))
radii.append(f"{round(source_radius * 1e6)} µm")
intensities_gvs = []
intensities_far_field_schell = []
for i, source_radius in enumerate(source_radii):
radii += [f"{round(source_radius * 1e6)} µm"] * 7
source_intensity = create_circular_mask(num_pixels_x * 2, dx_source, source_radius, num_pixels_y * 2, dy_source)
cpr_filter = create_circular_mask(num_pixels_x * 2, dx_source, cpr_radii[i], num_pixels_y * 2, dy_source)
complex_coherence_factor = compute_complex_coherence_factor(source_intensity)
intensity = np.roll(normalize(crop_center(
compute_gvs(
object,
complex_coherence_factor,
z_source_object,
z_object_camera,
wavelength,
dx_object
),
num_pixels_y, num_pixels_x
)), -1) # Roll to shift the array to the center
intensities_gvs.append(intensity)
intensities.append(intensity[num_pixels_y // 2, :])
names.append("GVS")
intensities.append(normalize(crop_center(
compute_shifted_fresnel_gvs(
object,
complex_coherence_factor,
z_source_object,
z_object_camera,
wavelength,
dx_camera,
dx_object=dx_object
),
num_pixels_y, num_pixels_x
)[num_pixels_y // 2, :]))
names.append("Shifted Fresnel GVS")
intensities.append(normalize(crop_center(
compute_fresnel_convolution_gvs(
object_approximation,
complex_coherence_factor,
z_source_object,
z_object_camera,
dx_camera,
wavelength,
),
num_pixels_y, num_pixels_x
)[num_pixels_y // 2, :]))
names.append("Fresnel Convolution GVS")
intensities.append(normalize(crop_center(
compute_angular_spectrum_gvs(
object_approximation,
complex_coherence_factor,
z_source_object,
z_object_camera,
dx_camera,
wavelength,
),
num_pixels_y, num_pixels_x
)[num_pixels_y // 2, :]))
names.append("Angular Spectrum GVS")
intensities.append(np.roll(normalize(crop_center(
compute_intensity_point_sources_dask(
object_approximation,
source_intensity,
dx_source,
dx_camera,
z_source_object,
z_object_camera,
wavelength,
sample_size=10000
),
num_pixels_y, num_pixels_x
)[num_pixels_y // 2, :]), 0))
names.append("Finite Point Sources")
intensities.append(np.roll(normalize(crop_center(
compute_intensity_cpr(
cpr_filter,
dx_camera,
object_approximation,
z_object_camera,
wavelength,
sample_size=30
),
num_pixels_y, num_pixels_x
)[num_pixels_y // 2, :]), 0))
names.append("CP-R (M = 30)")
intensities.append(np.roll(normalize(crop_center(
compute_intensity_cpr(
cpr_filter,
dx_camera,
object_approximation,
z_object_camera,
wavelength,
sample_size=300
),
num_pixels_y, num_pixels_x
)[num_pixels_y // 2, :]), 0))
names.append("CP-R (M = 300)")
intensity = np.roll(normalize(crop_center(
compute_schells_theorem(object, complex_coherence_factor, z_object_camera, wavelength),
num_pixels_y, num_pixels_x
)), -1) # Roll to shift the array to the center
intensities_far_field_schell.append(intensity)
x = np.linspace(-num_pixels_x // 2, num_pixels_x // 2, num_pixels_x) * dx_camera
df = pd.DataFrame({
"Normalized Intensity": intensities,
"Source Radius": radii,
"Model": names,
})
results = []
for source_radius in source_radii:
df_mse = df[(df["Source Radius"] == f"{round(source_radius * 1e6)} µm")]
integral_1 = df_mse[df_mse["Model"] == "Integral Computation 1"].iloc[0]["Normalized Intensity"]
integral_2 = df_mse[df_mse["Model"] == "Integral Computation 2"].iloc[0]["Normalized Intensity"]
df_mse = df_mse.loc[~df["Model"].isin(["Integral Computation 1", "Integral Computation 2"])]
df_mse_1 = df_mse.loc[
(df["Model"].isin(["GVS", "Shifted Fresnel GVS"]))
]
df_mse_1 = df_mse_1.groupby("Model")["Normalized Intensity"].apply(
lambda s: mean_squared_error(s.iloc[0], integral_1)
)
df_mse_2 = df_mse.loc[
(df["Model"].isin([
"Fresnel Convolution GVS",
"Angular Spectrum GVS",
"Finite Point Sources",
"CP-R (M = 30)",
"CP-R (M = 300)",
]))
]
df_mse_2 = df_mse_2.groupby("Model")["Normalized Intensity"].apply(
lambda s: mean_squared_error(s.iloc[0], integral_2)
)
df_mse = pd.concat([df_mse_1, df_mse_2])
df_mse.name = f"Radius = {round(source_radius * 1e6)} µm"
results.append(df_mse.reset_index())
df_mse = results[0]
for result in results[1:]:
df_mse = df_mse.merge(result, on="Model", how="outer")
display(Markdown(f"## 2D diffraction"))
display(display_image(
intensities_gvs,
show_as_grid=True,
normalize_all_frames=True,
grid_titles=[f"{round(r * 1e6)} µm" for r in source_radii],
title="GVS",
width=1000,
height=300
))
display(display_image(
intensities_far_field_schell,
show_as_grid=True,
normalize_all_frames=True,
grid_titles=[f"{round(r * 1e6)} µm" for r in source_radii],
title="Schell (far field)",
width=1000,
height=300
))
display(Markdown(f"## Cross sections"))
display(Markdown(f"### Mean Squared Error:"))
display(df_mse)
df = df.explode("Normalized Intensity")
df["x (m)"] = x.tolist() * len(intensities)
fig = px.line(df, x="x (m)", y="Normalized Intensity", color="Model", facet_col="Source Radius")
fig.update_traces(selector=dict(name="Integral Computation 1"), line=dict(width=3))
fig.update_traces(selector=dict(name="Integral Computation 2"), line=dict(width=3))
fig.update_traces(selector=dict(name="GVS"), line=dict(dash="dot"))
fig.update_traces(selector=dict(name="Shifted Fresnel GVS"), line=dict(dash="dot"))
fig.update_traces(selector=dict(name="Fresnel Convolution GVS"), line=dict(dash="dot"))
fig.update_traces(selector=dict(name="Angular Spectrum GVS"), line=dict(dash="dot"))
fig.update_traces(selector=dict(name="Finite Point Sources"), line=dict(dash="longdashdot"))
fig.update_traces(selector=dict(name="CP-R (M = 30)"), line=dict(dash="longdashdot", color="#1f77b4"), opacity=0.5)
fig.update_traces(selector=dict(name="CP-R (M = 300)"), line=dict(dash="longdashdot"), opacity=0.5)
fig.update_layout(margin_t=20, margin_l=0, margin_r=0)
fig.update_xaxes(showticklabels=False, title=None)
fig.update_yaxes(showticklabels=False)
fig.update_legends(title_text=None)
fig.show()
num_pixels_x = num_pixels_y = 256
z_source_object = 10 * cm
source_radii = [100 * um, 250 * um, 500 * um, 1 * mm]
cpr_radii = [50 * um, 200 * um, 300 * um, 500 * um]
pinhole_radius = 25 * um
pinhole_separation = 200 * um
dx_camera = dy_camera = 5 * um
z_object_camera = 10 * mm
show_diffraction()
2D diffraction¶
Cross sections¶
Mean Squared Error:¶
| Model | Radius = 100 µm | Radius = 250 µm | Radius = 500 µm | Radius = 1000 µm | |
|---|---|---|---|---|---|
| 0 | Angular Spectrum GVS | 0.000122 | 0.000053 | 0.000045 | 0.000043 |
| 1 | CP-R (M = 30) | 0.004074 | 0.007874 | 0.027226 | 0.035960 |
| 2 | CP-R (M = 300) | 0.003913 | 0.007687 | 0.004589 | 0.031018 |
| 3 | Finite Point Sources | 0.000122 | 0.000053 | 0.000046 | 0.000051 |
| 4 | Fresnel Convolution GVS | 0.000121 | 0.000053 | 0.000045 | 0.000043 |
| 5 | GVS | 0.000035 | 0.000013 | 0.000013 | 0.000015 |
| 6 | Shifted Fresnel GVS | 0.000008 | 0.000002 | 0.000004 | 0.000013 |
Arbitrary shifts and pixel pitches¶
This section shows the off-axis propagation and arbitrary pixel pitch propagation by using the shifted fresnel GVS diffraction.
Pixel pitch variation¶
def display_three_images():
image_1_ = base64.b64encode(image_1.data).decode("utf-8")
image_2_ = base64.b64encode(image_2.data).decode("utf-8")
image_3_ = base64.b64encode(image_3.data).decode("utf-8")
html = f"""
<div style="display: flex; justify-content: space-around; align-items: center;">
<img src="data:image/png;base64,{image_1_}" style="width: 33%;">
<img src="data:image/png;base64,{image_2_}" style="width: 33%;">
<img src="data:image/png;base64,{image_3_}" style="width: 33%;">
</div>
"""
return HTML(html)
num_pixels_x = num_pixels_y = 1024
dx_camera = 10 * um
dx_object = 10 * um
shift_object_x = 0.0 * um
shift_object_y = 0.0 * um
shift_camera_x = 0.0 * um
shift_camera_y = 0.0 * um
z_object_camera = 40 * cm
z_source_object = 10 * cm
pinhole_radius = 50 * um
pinhole_separation = 300 * um
source_radius = 100 * um
dx_source = z_source_object / z_object_camera * dx_camera
object = create_double_pinhole_mask(
num_pixels_x,
dx_object,
radius=pinhole_radius,
separation=pinhole_separation,
num_pixels_y=num_pixels_y,
)
source_intensity = create_circular_mask(num_pixels_x * 2, dx_source, source_radius)
complex_coherence_factor = compute_complex_coherence_factor(source_intensity)
intensity_1 = crop_center(
compute_shifted_fresnel_gvs(
object,
complex_coherence_factor,
z_source_object,
z_object_camera,
wavelength,
dx_camera,
dx_object=dx_object
),
num_pixels_y, num_pixels_x
)
x = np.linspace(-num_pixels_x // 2, num_pixels_x // 2, num_pixels_x) * dx_camera
image_1 = display_image(intensity_1, x=x, y=-x, title="$dx = 10 \ \mu m$")
dx_camera = 5 * um
dx_source = z_source_object / z_object_camera * dx_camera
source_intensity = create_circular_mask(num_pixels_x * 2, dx_source, source_radius)
complex_coherence_factor = compute_complex_coherence_factor(source_intensity)
intensity_2 = crop_center(
compute_shifted_fresnel_gvs(
object,
complex_coherence_factor,
z_source_object,
z_object_camera,
wavelength,
dx_camera,
dx_object=dx_object
),
num_pixels_y, num_pixels_x
)
x = np.linspace(-num_pixels_x // 2, num_pixels_x // 2, num_pixels_x) * dx_camera
image_2 = display_image(intensity_2, x=x, y=-x, title="$dx = 5 \ \mu m$")
dx_camera = 2.5 * um
dx_source = z_source_object / z_object_camera * dx_camera
source_intensity = create_circular_mask(num_pixels_x * 2, dx_source, source_radius)
complex_coherence_factor = compute_complex_coherence_factor(source_intensity)
intensity_3 = crop_center(
compute_shifted_fresnel_gvs(
object,
complex_coherence_factor,
z_source_object,
z_object_camera,
wavelength,
dx_camera,
dx_object=dx_object
),
num_pixels_y, num_pixels_x
)
x = np.linspace(-num_pixels_x // 2, num_pixels_x // 2, num_pixels_x) * dx_camera
image_3 = display_image(intensity_3, x=x, y=-x, title="$dx = 2.5 \ \mu m$")
display_three_images()
Shifted diffraction¶
Here we compute the diffraction with four different shifts, corresponding to the four quadrants of a cartesian plane.
num_pixels_x = 512
dx_camera = 5 * um
dx_source = z_source_object / z_object_camera * dx_camera
object = create_double_pinhole_mask(
num_pixels_x,
dx_object,
radius=pinhole_radius,
separation=pinhole_separation,
)
source_intensity = create_circular_mask(num_pixels_x * 2, dx_source, source_radius)
complex_coherence_factor = compute_complex_coherence_factor(source_intensity)
shift_object_y = -1.25 * mm
shift_object_x = 1.25 * mm
intensity_1 = crop_center(
compute_shifted_fresnel_gvs(
object,
complex_coherence_factor,
z_source_object,
z_object_camera,
wavelength,
dx_camera,
dx_object=dx_object,
shift_object_x=shift_object_x,
shift_object_y=shift_object_y
),
num_pixels_x, num_pixels_x
)
shift_object_y = 1.25 * mm
intensity_2 = crop_center(
compute_shifted_fresnel_gvs(
object,
complex_coherence_factor,
z_source_object,
z_object_camera,
wavelength,
dx_camera,
dx_object=dx_object,
shift_object_x=shift_object_x,
shift_object_y=shift_object_y
),
num_pixels_x, num_pixels_x
)
shift_object_x = -1.25 * mm
intensity_3 = crop_center(
compute_shifted_fresnel_gvs(
object,
complex_coherence_factor,
z_source_object,
z_object_camera,
wavelength,
dx_camera,
dx_object=dx_object,
shift_object_x=shift_object_x,
shift_object_y=shift_object_y
),
num_pixels_x, num_pixels_x
)
shift_object_y = -1.25 * mm
intensity_4 = crop_center(
compute_shifted_fresnel_gvs(
object,
complex_coherence_factor,
z_source_object,
z_object_camera,
wavelength,
dx_camera,
dx_object=dx_object,
shift_object_x=shift_object_x,
shift_object_y=shift_object_y
),
num_pixels_x, num_pixels_x
)
display_image(
[intensity_2, intensity_3, intensity_1, intensity_4],
show_as_grid=True,
facet_col_wrap=2,
grid_titles=[""] * 4,
width=490,
height=512
)
Random objects¶
Here we compare the fourier-based models with the integral computation using small random objects. We separate the GVS and shifted fresnel models from the angular spectrum and fresnel convolution models for fair comparison.
Random Object¶
num_pixels = 40
dx_camera = 5 * um
z_object_camera = 2 * cm
z_source_object = 10 * cm
source_radius = 50 * um
dx_object = (wavelength * z_object_camera) / (2 * dx_camera * num_pixels)
dx_source = (wavelength * z_source_object) / (2 * dx_object * num_pixels)
random_object = np.pad(np.random.rand(num_pixels, num_pixels), num_pixels // 2)
display_image(random_object)
Models (1)¶
source_intensity = create_circular_mask(num_pixels * 2, dx_source, source_radius)
complex_coherence_factor = compute_complex_coherence_factor(source_intensity)
intensity_integral_1 = normalize(
compute_mutual_intensity_integral_cross_section(
random_object,
np.fft.ifftshift(complex_coherence_factor),
z_source_object,
z_object_camera,
wavelength,
dx_camera,
dx_object,
chunk_size=300
)
)[num_pixels - num_pixels // 2:num_pixels + num_pixels // 2]
intensity_gvs = np.roll(normalize(crop_center(
compute_gvs(
random_object,
complex_coherence_factor,
z_source_object,
z_object_camera,
wavelength,
dx_object
),
num_pixels, num_pixels
)[num_pixels // 2, :]), -1) # Roll to shift the array to the center
intensity_shifted_fresnel = normalize(crop_center(
compute_shifted_fresnel_gvs(
random_object,
complex_coherence_factor,
z_source_object,
z_object_camera,
wavelength,
dx_camera,
dx_object=dx_object
),
num_pixels, num_pixels
)[num_pixels // 2, :])
fig = (
px.line(
y=[
intensity_integral_1,
intensity_gvs,
intensity_shifted_fresnel,
])
.update_yaxes(title="Normalized Intensity")
.update_xaxes(title="Pixel Index")
.update_legends(
title=None,
yanchor="top",
y=0.99,
xanchor="right",
x=0.99
)
)
fig.data[0].update(name="Integral Computation", line=dict(width=3))
fig.data[1].update(name="GVS", line=dict(dash="dot"))
fig.data[2].update(name="Shifted Fresnel GVS", line=dict(dash="dot"))
fig.show()
Models (2)¶
z_object_camera = 0.2 * cm
source_radius = 500 * um
dx_object = (wavelength * z_object_camera) / (2 * dx_camera * num_pixels)
dx_source = (wavelength * z_source_object) / (2 * dx_object * num_pixels)
random_object = np.pad(np.random.rand(num_pixels, num_pixels), num_pixels // 2)
source_intensity = create_circular_mask(num_pixels * 2, dx_source, source_radius)
source_intensity_point_sources = create_circular_mask(num_pixels * 2, dx_camera, source_radius)
complex_coherence_factor = compute_complex_coherence_factor(source_intensity)
complex_coherence_factor_alternative = zoom_matrix(complex_coherence_factor, dx_object / dx_camera)
intensity_integral_2 = normalize(
compute_mutual_intensity_integral_cross_section(
random_object,
np.fft.ifftshift(complex_coherence_factor_alternative),
z_source_object,
z_object_camera,
wavelength,
dx_camera,
chunk_size=300
)
)[num_pixels - num_pixels // 2:num_pixels + num_pixels // 2]
intensity_fresnel_convolution = normalize(crop_center(
compute_fresnel_convolution_gvs(
random_object,
complex_coherence_factor,
z_source_object,
z_object_camera,
dx_camera,
wavelength
),
num_pixels, num_pixels
)[num_pixels // 2, :])
intensity_angular_spectrum_fresnel = normalize(crop_center(
compute_angular_spectrum_gvs(
random_object,
complex_coherence_factor,
z_source_object,
z_object_camera,
dx_camera,
wavelength
),
num_pixels, num_pixels
)[num_pixels // 2, :])
fig = (
px.line(
y=[
intensity_integral_2,
intensity_fresnel_convolution,
intensity_angular_spectrum_fresnel,
]
)
.update_yaxes(title="Normalized Intensity")
.update_xaxes(title="Pixel Index")
.update_legends(title_text="Model")
)
fig.data[0].update(name="Integral Computation", line=dict(width=3))
fig.data[1].update(name="Fresnel Convolution Schell", line=dict(dash="dot"))
fig.data[2].update(name="Angular Spectrum Schell", line=dict(dash="dot"))
fig.show()
USAF animation¶
def zoom_image(image):
return zoom(
image,
(
num_pixels_y / image.shape[0],
num_pixels_x / image.shape[1]
)
)
def show_diffraction():
object = normalize(zoom_image(imageio.imread("img/usaf.png", mode="L")))
object = np.pad(object, ((num_pixels_y // 2, num_pixels_y // 2), (num_pixels_x // 2, num_pixels_x // 2)))
ny_padded, nx_padded = object.shape
all_intensities = []
for z_object_camera in propagation_distances:
dx_object = (wavelength * z_object_camera) / (2 * dx_camera * num_pixels_x)
dx_source = (wavelength * z_source_object) / (2 * dx_object * num_pixels_x)
dy_object = (wavelength * z_object_camera) / (2 * dy_camera * num_pixels_y)
dy_source = (wavelength * z_source_object) / (2 * dy_object * num_pixels_y)
intensities = []
intensities.append(normalize(crop_center(abs(compute_angular_spectrum(
object, z_object_camera, dx_camera, wavelength, spherical=True, z_source_object=z_source_object
) ** 2))))
source_intensity = create_circular_mask(nx_padded, dx_source, source_radii[0], ny_padded, dy_source)
complex_coherence_factor = compute_complex_coherence_factor(source_intensity)
intensities.append(normalize(crop_center(compute_angular_spectrum_gvs(
object, complex_coherence_factor, z_source_object, z_object_camera, dx_camera, wavelength
))))
source_intensity = create_circular_mask(nx_padded, dx_source, source_radii[1], ny_padded, dy_source)
complex_coherence_factor = compute_complex_coherence_factor(source_intensity)
intensities.append(normalize(crop_center(compute_angular_spectrum_gvs(
object, complex_coherence_factor, z_source_object, z_object_camera, dx_camera, wavelength
))))
source_intensity = create_circular_mask(nx_padded, dx_source, source_radii[2], ny_padded, dy_source)
complex_coherence_factor = compute_complex_coherence_factor(source_intensity)
intensities.append(normalize(crop_center(compute_angular_spectrum_gvs(
object, complex_coherence_factor, z_source_object, z_object_camera, dx_camera, wavelength
))))
all_intensities.append(intensities)
fig = display_image(
all_intensities,
convert_to=None,
show_as_grid=True,
show_animation=True,
propagation_distances=propagation_distances.tolist(),
grid_titles=[
"Angular Spectrum (Coherent)",
"Angular Spectrum GVS (High Coherence)",
"Angular Spectrum GVS (Medium Coherence)",
"Angular Spectrum GVS (Low Coherence)",
],
return_figure=True,
width=None,
height=430
)
fig.update_layout(
sliders=[
{
"currentvalue": {
"font": {"size": 13},
"xanchor": "center",
"offset": 20,
},
"pad": {"t": 0}
}
],
)
for annotation in fig.layout.annotations:
annotation.font.size = 9
fig.frames = fig.frames[::2]
fig.show()
num_pixels_x = 1000
num_pixels_y = 750
z_source_object = 10 * cm
source_radii = [80 * um, 250 * um, 1000 * um]
dx_camera = dy_camera = 2 * um
propagation_distances = np.linspace(2 * mm, 20 * mm, 20).tolist()
propagation_distances = np.concat([propagation_distances, propagation_distances[::-1]])
show_diffraction()