PTI example jupyter notebook to script and .html (#160)

* first-pass script and .html

* html -> pdf

* correct a docstring

* remove obsolete waveOrderWriter calls

* stub for iohub writer + correcting names of quantities

* Save reconstruction to OME zarr for viewing in napari

* fix mean/differential permittivity mistake

* fix comments

* don't track jupyter notebook

* clean path

* update pdf and readme

---------

Co-authored-by: Shalin Mehta <[email protected]>

examples/documentation/PTI_experiment/README.md

@@ -1,6 +1,11 @@
-The notebooks in this folder require data from this source:
+# PTI reconstructions
+
+The demos in this folder illustrate reconstruction algorithm for permittivity tensor imaging (PTI). The main demo script is [PTI_Experiment_Recon3D_anisotropic_target_small.py](PTI_Experiment_Recon3D_anisotropic_target_small.py), which illustrates reconstruction of the permittivity tensor of an anisotropic glass target. The reconstructed results are rendered as PDF file [PTI_Experiment_Recon3D_anisotropic_target_small.pdf](PTI_Experiment_Recon3D_anisotropic_target_small.pdf).
+
+The remaining notebooks illustrate the reconstruction of permittivity tensor of cells and tissues. However they are not maintained and require an old version of this repository from https://github.com/mehta-lab/waveorder/tree/1.0.0b0.
+
+The scripts and notebooks in this folder require data from the image data repository:
 https://www.ebi.ac.uk/biostudies/bioimages/studies/S-BIAD1063
 
-One notebook `PTI_Experiment_Recon3D_anisotropic_target_small.ipynb` will run as is with the latest version of the code and a ~500 MB download from https://www.ebi.ac.uk/biostudies/files/S-BIAD1063/PTI-BIA/Anisotropic_target_small.zip.
+The demo script [PTI_Experiment_Recon3D_anisotropic_target_small.py](PTI_Experiment_Recon3D_anisotropic_target_small.py) will run as is with the latest version of the code and a ~500 MB download from https://www.ebi.ac.uk/biostudies/files/S-BIAD1063/PTI-BIA/Anisotropic_target_small.zip.
 
-The remaining notebooks are not maintained, so they will require a download and an old version of this repository from https://github.com/mehta-lab/waveorder/tree/1.0.0b0.

examples/models/isotropic_thin_3d.py

@@ -16,7 +16,7 @@
     "wavelength_illumination": 0.532,
     "index_of_refraction_media": 1.3,
 }
-phantom_arguments = {"index_of_refraction_sample": 1.50, "sphere_radius": 5}
+phantom_arguments = {"index_of_refraction_sample": 1.33, "sphere_radius": 5}
 z_shape = 100
 z_pixel_size = 0.25
 transfer_function_arguments = {

examples/models/isotropic_thin_3d_resolution.py

@@ -0,0 +1,88 @@
+# 3D partially coherent optical diffraction tomography (ODT) simulation
+# J. M. Soto, J. A. Rodrigo, and T. Alieva, "Label-free quantitative
+# 3D tomographic imaging for partially coherent light microscopy," Opt. Express
+# 25, 15699-15712 (2017)
+
+import napari
+import numpy as np
+from waveorder import util
+from waveorder.models import isotropic_thin_3d
+
+# Parameters
+# all lengths must use consistent units e.g. um
+simulation_arguments = {
+    "yx_shape": (256, 256),
+    "yx_pixel_size": 6.5 / 63,
+    "wavelength_illumination": 0.532,
+    "index_of_refraction_media": 1.3,
+}
+phantom_arguments = {"index_of_refraction_sample": 1.4, "sphere_radius": 0.05}
+z_shape = 100
+z_pixel_size = 0.25
+transfer_function_arguments = {
+    "z_position_list": (np.arange(z_shape) - z_shape // 2) * z_pixel_size,
+    "numerical_aperture_illumination": 0.9,
+    "numerical_aperture_detection": 1.2,
+}
+
+# Create a phantom
+yx_absorption, yx_phase = isotropic_thin_3d.generate_test_phantom(
+    **simulation_arguments, **phantom_arguments
+)
+
+# Calculate transfer function
+(
+    absorption_2d_to_3d_transfer_function,
+    phase_2d_to_3d_transfer_function,
+) = isotropic_thin_3d.calculate_transfer_function(
+    **simulation_arguments, **transfer_function_arguments
+)
+
+# Display transfer function
+viewer = napari.Viewer()
+zyx_scale = np.array(
+    [
+        z_pixel_size,
+        simulation_arguments["yx_pixel_size"],
+        simulation_arguments["yx_pixel_size"],
+    ]
+)
+isotropic_thin_3d.visualize_transfer_function(
+    viewer,
+    absorption_2d_to_3d_transfer_function,
+    phase_2d_to_3d_transfer_function,
+)
+input("Showing OTFs. Press <enter> to continue...")
+viewer.layers.select_all()
+viewer.layers.remove_selected()
+
+# Simulate
+zyx_data = isotropic_thin_3d.apply_transfer_function(
+    yx_absorption,
+    yx_phase,
+    absorption_2d_to_3d_transfer_function,
+    phase_2d_to_3d_transfer_function,
+)
+
+# Reconstruct
+(
+    yx_absorption_recon,
+    yx_phase_recon,
+) = isotropic_thin_3d.apply_inverse_transfer_function(
+    zyx_data,
+    absorption_2d_to_3d_transfer_function,
+    phase_2d_to_3d_transfer_function,
+)
+
+# Display
+arrays = [
+    (yx_absorption, "Phantom - absorption"),
+    (yx_phase, "Phantom - phase"),
+    (zyx_data, "Data"),
+    (yx_absorption_recon, "Reconstruction - absorption"),
+    (yx_phase_recon, "Reconstruction - phase"),
+]
+
+for array in arrays:
+    viewer.add_image(array[0].cpu().numpy(), name=array[1])
+input("Showing object, data, and recon. Press <enter> to quit...")

waveorder/models/isotropic_thin_3d.py

@@ -28,7 +28,7 @@ def generate_test_phantom(
         / wavelength_illumination
     )  # phase in radians
 
-    yx_absorption = 0.99 * sphere[1]
+    yx_absorption = 0.02 * sphere[1]
 
     return yx_absorption, yx_phase
 
@@ -113,6 +113,28 @@ def visualize_transfer_function(
     viewer.dims.order = (0, 1, 2)
 
 
+def visualize_point_spread_function(
+    viewer,
+    absorption_2d_to_3d_transfer_function,
+    phase_2d_to_3d_transfer_function,
+):
+    arrays = [
+        (torch.fft.ifftn(absorption_2d_to_3d_transfer_function), "absorb PSF"),
+        (torch.fft.ifftn(phase_2d_to_3d_transfer_function), "phase PSF"),
+    ]
+
+    for array in arrays:
+        lim = 0.5 * torch.max(torch.abs(array[0]))
+        viewer.add_image(
+            torch.fft.ifftshift(array[0], dim=(1, 2)).cpu().numpy(),
+            name=array[1],
+            colormap="bwr",
+            contrast_limits=(-lim, lim),
+            scale=(1, 1, 1),
+        )
+    viewer.dims.order = (0, 1, 2)
+
+
 def apply_transfer_function(
     yx_absorption,
     yx_phase,

waveorder/optics.py

@@ -1147,7 +1147,7 @@ def unit_conversion_from_scattering_potential_to_permittivity(
 ):
     """
 
-    compute the optic sign probability from the reconstructed material tendancy
+    Convert the scattering potential of the specimen to the relative permittivity.
 
     Parameters
     ----------