From be20c8fab041ac24c8e70b85a361ca33405dde15 Mon Sep 17 00:00:00 2001 From: Talon Chandler Date: Mon, 27 May 2024 19:37:06 -0400 Subject: [PATCH 1/4] Authorship and release updates (#165) * author and release updates * nits --------- Co-authored-by: Shalin Mehta <2934183+mattersoflight@users.noreply.github.com> --- CITATION.cff | 29 ++++++++++++++++------------- 1 file changed, 16 insertions(+), 13 deletions(-) diff --git a/CITATION.cff b/CITATION.cff index bbd684fe..dc624506 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -4,30 +4,32 @@ cff-version: 1.2.0 title: waveorder message: >- - A generalist library of wave optical models and inverse - algorithms for label-free imaging of density & - orientation. + Wave optical models and inverse algorithms for label-agnostic imaging of density & orientation. type: software authors: - - given-names: Li-Hao - family-names: Yeh - affiliation: CZ Biohub - orcid: 'https://orcid.org/0000-0003-2803-5996' - given-names: Talon family-names: Chandler - affiliation: CZ Biohub + affiliation: Chan Zuckerberg Biohub San Francisco orcid: 'https://orcid.org/0000-0002-3033-674X' + - given-names: Li-Hao + family-names: Yeh + affiliation: Chan Zuckerberg Biohub San Francisco + orcid: 'https://orcid.org/0000-0003-2803-5996' - given-names: Ivan family-names: Ivanov - affiliation: CZ Biohub + affiliation: Chan Zuckerberg Biohub San Francisco orcid: 'https://orcid.org/0000-0002-4675-5287' - given-names: Cameron family-names: Foltz - affiliation: CZ Biohub + affiliation: Chan Zuckerberg Biohub San Francisco orcid: 'https://orcid.org/0000-0001-8933-2172' + - given-names: Ziwen + family-names: Liu + affiliation: Chan Zuckerberg Biohub San Francisco + orcid: 'https://orcid.org/0000-0001-7482-1299' - given-names: Shalin family-names: Mehta - affiliation: CZ Biohub + affiliation: Chan Zuckerberg Biohub San Francisco orcid: 'https://orcid.org/0000-0002-2542-3582' identifiers: - type: doi @@ -67,6 +69,7 @@ keywords: - qlipp - mipolscope - simulation + - reconstruction license: BSD-3-Clause -version: 1.0.0rc1 -date-released: '2023-02-22' +version: 2.1.0 +date-released: '2024-03-12' From c691f2f0411ead4c5e6f1ed3fa44831d8763c3a4 Mon Sep 17 00:00:00 2001 From: Ziwen Liu <67518483+ziw-liu@users.noreply.github.com> Date: Thu, 13 Jun 2024 13:51:28 -0700 Subject: [PATCH 2/4] Bump Python (#168) * migrate to pyproject.toml * bump python * remove inaccurate topics from pyproject.toml * Revert "remove inaccurate topics from pyproject.toml" This reverts commit 46a5e72500501b9bff3b8d5fb1fe467d44b3bd44. * remove inaccurate topic metadata --- .github/workflows/pytests.yml | 2 +- README.md | 2 +- pyproject.toml | 59 ++++++++++++++++++++++++++++++++++- setup.cfg | 52 ------------------------------ setup.py | 3 -- 5 files changed, 60 insertions(+), 58 deletions(-) delete mode 100644 setup.cfg delete mode 100644 setup.py diff --git a/.github/workflows/pytests.yml b/.github/workflows/pytests.yml index deb6faa6..8028c79d 100644 --- a/.github/workflows/pytests.yml +++ b/.github/workflows/pytests.yml @@ -11,7 +11,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: ['3.8', '3.9', '3.10'] + python-version: ['3.10', '3.11', '3.12'] steps: - name: Checkout repo diff --git a/README.md b/README.md index 5bbb9530..07674395 100644 --- a/README.md +++ b/README.md @@ -51,7 +51,7 @@ Please cite this repository, along with the relevant preprint or paper, if you u (Optional but recommended) install [anaconda](https://www.anaconda.com/products/distribution) and create a virtual environment: ```sh -conda create -y -n waveorder python=3.9 +conda create -y -n waveorder python=3.11 conda activate waveorder ``` diff --git a/pyproject.toml b/pyproject.toml index 7caead92..1a57bb79 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,9 +5,66 @@ build-backend = "setuptools.build_meta" [tool.setuptools_scm] write_to = "waveorder/_version.py" +[project] +name = "waveorder" +description = "Wave-optical simulations and deconvolution of optical properties" +readme = "README.md" +requires-python = ">=3.10" +license = { file = "LICENSE" } +authors = [{ name = "CZ Biohub SF", email = "compmicro@czbiohub.org" }] +maintainers = [ + { name = "Talon Chandler", email = "talon.chandler@czbiohub.org" }, + { name = "Shalin Mehta", email = "shalin.mehta@czbiohub.org" }, +] +keywords = [ + "simulation", + 'optics', + 'phase', + 'scattering', + 'polarization', + 'label-free', + 'permittivity', + "reconstruction-algorithm", + 'qlipp', + 'mipolscope', + 'permittivity-tensor-imaging', +] +classifiers = [ + "Development Status :: 4 - Beta", + "Intended Audience :: Science/Research", + "License :: OSI Approved :: BSD License", + "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3.10", + "Programming Language :: Python :: 3.11", + "Programming Language :: Python :: 3.12", + "Topic :: Scientific/Engineering", + "Topic :: Scientific/Engineering :: Image Processing", + "Topic :: Scientific/Engineering :: Bio-Informatics", + "Operating System :: POSIX :: Linux", + "Operating System :: Microsoft :: Windows", + "Operating System :: MacOS", +] +dependencies = [ + "numpy>=1.21", + "matplotlib>=3.1.1", + "scipy>=1.3.0", + "pywavelets>=1.1.1", + "ipywidgets>=7.5.1", + "torch>=2.2.1", +] +dynamic = ["version"] + +[project.optional-dependencies] +dev = ["pytest", "pytest-cov"] + +[project.urls] +Homepage = "https://github.com/mehta-lab/waveorder" +Repository = "https://github.com/mehta-lab/waveorder" +Issues = "https://github.com/mehta-lab/waveorder/issues" + [tool.black] line-length = 79 [tool.isort] profile = "black" -line_length = 79 \ No newline at end of file +line_length = 79 diff --git a/setup.cfg b/setup.cfg deleted file mode 100644 index 8c6a13ef..00000000 --- a/setup.cfg +++ /dev/null @@ -1,52 +0,0 @@ -[metadata] -name = waveorder -author = Computational Microscopy Platform, CZ Biohub -author_email = shalin.mehta@czbiohub.org -url = https://github.com/mehta-lab/waveorder -license = BSD 3-Clause License -description = Wave optical simulations and deconvolution of optical properties -long_description = file: README.md -long_description_content_type = text/markdown -classifiers = - Development Status :: 4 - Beta - Intended Audience :: Science/Research - License :: OSI Approved :: BSD License - Programming Language :: Python :: 3 :: Only - Programming Language :: Python :: 3.8 - Programming Language :: Python :: 3.9 - Programming Language :: Python :: 3.10 - Topic :: Scientific/Engineering - Topic :: Scientific/Engineering :: Image Processing - Topic :: Scientific/Engineering :: Visualization - Topic :: Scientific/Engineering :: Information Analysis - Topic :: Scientific/Engineering :: Bio-Informatics - Operating System :: Microsoft :: Windows - Operating System :: POSIX - Operating System :: Unix - Operating System :: MacOS -project_urls = - Bug Tracker = https://github.com/mehta-lab/waveorder/issues - Documentation = https://github.com/mehta-lab/waveorder - Source Code = https://github.com/mehta-lab/waveorder - User Support = https://github.com/mehta-lab/waveorder/issues - -[options] -package = find: -include_package_data = True -python_requires = >=3.8 -setup_requires = setuptools_scm -# add your package requirements here -install_requires = - numpy>=1.21 - matplotlib>=3.1.1 - scipy>=1.3.0 - pywavelets>=1.1.1 - ipywidgets>=7.5.1 - torch>=2.0.0 - -[options.extras_require] -dev = - flake8 - pytest>=5.0.0 - pytest-cov - \ No newline at end of file diff --git a/setup.py b/setup.py deleted file mode 100644 index 60684932..00000000 --- a/setup.py +++ /dev/null @@ -1,3 +0,0 @@ -from setuptools import setup - -setup() From 45514b8963f1bb446220eb8701ca0680e0c8ff9d Mon Sep 17 00:00:00 2001 From: Talon Chandler Date: Fri, 5 Jul 2024 13:52:51 -0700 Subject: [PATCH 3/4] Merge `alg-dev` into main ahead of release (#171) * Phase reconstruction is invariant to voxel-size (#164) * fix bug finding focus in stack with only one slice * refactor for clarify * formatting * print -> warnings.warn * test single-slice case * fix test bugs * z-scale-invariant test object * no rescaling on output * forward simulation takes a "brightness" - simulating real microscope * fix example script * add background parameter for fluorescence forward model * test voxel-size invariance * rename I_norm -> direct_intensity * refactor to clarify discretization factor * remove comment * fix fluorescence example bug * improved docsring --------- Co-authored-by: Ivan Ivanov * np.product -> np.prod --------- Co-authored-by: Ivan Ivanov --- .../models/isotropic_fluorescent_thick_3d.py | 4 +- examples/models/phase_thick_3d.py | 5 +- tests/models/test_phase_thick_3d.py | 81 ++++++++++++++++++- .../models/isotropic_fluorescent_thick_3d.py | 21 ++++- waveorder/models/phase_thick_3d.py | 35 ++++---- waveorder/optics.py | 16 ++-- 6 files changed, 127 insertions(+), 35 deletions(-) diff --git a/examples/models/isotropic_fluorescent_thick_3d.py b/examples/models/isotropic_fluorescent_thick_3d.py index 21a68169..01840798 100644 --- a/examples/models/isotropic_fluorescent_thick_3d.py +++ b/examples/models/isotropic_fluorescent_thick_3d.py @@ -6,13 +6,13 @@ # Parameters # all lengths must use consistent units e.g. um simulation_arguments = { - "zyx_shape": (100, 256, 256), + "zyx_shape": (200, 256, 256), "yx_pixel_size": 6.5 / 63, "z_pixel_size": 0.25, } phantom_arguments = {"sphere_radius": 5} transfer_function_arguments = { - "wavelength_illumination": 0.532, + "wavelength_emission": 0.532, "z_padding": 0, "index_of_refraction_media": 1.3, "numerical_aperture_detection": 1.2, diff --git a/examples/models/phase_thick_3d.py b/examples/models/phase_thick_3d.py index aa9bd4aa..54d7e43d 100644 --- a/examples/models/phase_thick_3d.py +++ b/examples/models/phase_thick_3d.py @@ -14,12 +14,12 @@ "zyx_shape": (100, 256, 256), "yx_pixel_size": 6.5 / 63, "z_pixel_size": 0.25, - "wavelength_illumination": 0.532, "index_of_refraction_media": 1.3, } phantom_arguments = {"index_of_refraction_sample": 1.50, "sphere_radius": 5} transfer_function_arguments = { "z_padding": 0, + "wavelength_illumination": 0.532, "numerical_aperture_illumination": 0.9, "numerical_aperture_detection": 1.2, } @@ -61,6 +61,7 @@ zyx_phase, real_potential_transfer_function, transfer_function_arguments["z_padding"], + brightness=1e3, ) # Reconstruct @@ -69,8 +70,6 @@ real_potential_transfer_function, imag_potential_transfer_function, transfer_function_arguments["z_padding"], - simulation_arguments["z_pixel_size"], - simulation_arguments["wavelength_illumination"], ) # Display diff --git a/tests/models/test_phase_thick_3d.py b/tests/models/test_phase_thick_3d.py index 224c96c6..d60c7fa6 100644 --- a/tests/models/test_phase_thick_3d.py +++ b/tests/models/test_phase_thick_3d.py @@ -1,5 +1,5 @@ import pytest - +import numpy as np from waveorder.models import phase_thick_3d @@ -20,3 +20,82 @@ def test_calculate_transfer_function(invert_phase_contrast): assert H_re.shape == (20 + 2 * z_padding, 100, 101) assert H_im.shape == (20 + 2 * z_padding, 100, 101) + + +# Helper function for testing reconstruction invariances +def simulate_phase_recon( + z_pixel_size_um=0.1, + yx_pixel_size_um=6.5 / 63, +): + + z_fov_um = 50 + yx_fov_um = 50 + + n_z = np.int32(z_fov_um / z_pixel_size_um) + n_yx = np.int32(yx_fov_um / yx_pixel_size_um) + + # Parameters + # all lengths must use consistent units e.g. um + simulation_arguments = { + "zyx_shape": (n_z, n_yx, n_yx), + "yx_pixel_size": yx_pixel_size_um, + "z_pixel_size": z_pixel_size_um, + "index_of_refraction_media": 1.3, + } + phantom_arguments = { + "index_of_refraction_sample": 1.40, + "sphere_radius": 5, + } + transfer_function_arguments = { + "z_padding": 0, + "wavelength_illumination": 0.532, + "numerical_aperture_illumination": 0.9, + "numerical_aperture_detection": 1.3, + } + + # Create a phantom + zyx_phase = phase_thick_3d.generate_test_phantom( + **simulation_arguments, **phantom_arguments + ) + + # Calculate transfer function + ( + real_potential_transfer_function, + imag_potential_transfer_function, + ) = phase_thick_3d.calculate_transfer_function( + **simulation_arguments, **transfer_function_arguments + ) + + # Simulate + zyx_data = phase_thick_3d.apply_transfer_function( + zyx_phase, + real_potential_transfer_function, + transfer_function_arguments["z_padding"], + brightness=1000, + ) + + # Reconstruct + zyx_recon = phase_thick_3d.apply_inverse_transfer_function( + zyx_data, + real_potential_transfer_function, + imag_potential_transfer_function, + transfer_function_arguments["z_padding"], + regularization_strength=1e-3, + ) + + Z, Y, X = zyx_phase.shape + recon_center = zyx_recon[Z // 2, Y // 2, X // 2].numpy() + + return recon_center + + +def test_phase_invariance(): + recon = simulate_phase_recon() + + # test z pixel size invariance + recon1 = simulate_phase_recon(z_pixel_size_um=0.3) + assert np.abs((recon1 - recon) / recon) < 0.02 + + # test yx pixel size invariance + recon2 = simulate_phase_recon(yx_pixel_size_um=0.7 * 6.5 / 63) + assert np.abs((recon2 - recon) / recon) < 0.02 \ No newline at end of file diff --git a/waveorder/models/isotropic_fluorescent_thick_3d.py b/waveorder/models/isotropic_fluorescent_thick_3d.py index 58234852..b52e71b5 100644 --- a/waveorder/models/isotropic_fluorescent_thick_3d.py +++ b/waveorder/models/isotropic_fluorescent_thick_3d.py @@ -81,7 +81,24 @@ def visualize_transfer_function(viewer, optical_transfer_function, zyx_scale): viewer.dims.order = (0, 1, 2) -def apply_transfer_function(zyx_object, optical_transfer_function, z_padding): +def apply_transfer_function( + zyx_object, optical_transfer_function, z_padding, background=10 +): + """Simulate imaging by applying a transfer function + + Parameters + ---------- + zyx_object : torch.Tensor + optical_transfer_function : torch.Tensor + z_padding : int + background : int, optional + constant background counts added to each voxel, by default 10 + + Returns + ------- + Simulated data : torch.Tensor + + """ if ( zyx_object.shape[0] + 2 * z_padding != optical_transfer_function.shape[0] @@ -99,7 +116,7 @@ def apply_transfer_function(zyx_object, optical_transfer_function, z_padding): zyx_data = zyx_obj_hat * optical_transfer_function data = torch.real(torch.fft.ifftn(zyx_data)) - data += 10 # Add a direct background + data += background # Add a direct background return data diff --git a/waveorder/models/phase_thick_3d.py b/waveorder/models/phase_thick_3d.py index 39555b50..ac29d356 100644 --- a/waveorder/models/phase_thick_3d.py +++ b/waveorder/models/phase_thick_3d.py @@ -12,7 +12,6 @@ def generate_test_phantom( zyx_shape, yx_pixel_size, z_pixel_size, - wavelength_illumination, index_of_refraction_media, index_of_refraction_sample, sphere_radius, @@ -24,12 +23,9 @@ def generate_test_phantom( radius=sphere_radius, blur_size=2 * yx_pixel_size, ) - zyx_phase = ( - sphere - * (index_of_refraction_sample - index_of_refraction_media) - * z_pixel_size - / wavelength_illumination - ) # phase in radians + zyx_phase = sphere * ( + index_of_refraction_sample - index_of_refraction_media + ) # refractive index increment return zyx_phase @@ -120,12 +116,19 @@ def visualize_transfer_function( def apply_transfer_function( - zyx_object, real_potential_transfer_function, z_padding + zyx_object, real_potential_transfer_function, z_padding, brightness ): # This simplified forward model only handles phase, so it resuses the fluorescence forward model # TODO: extend to absorption - return isotropic_fluorescent_thick_3d.apply_transfer_function( - zyx_object, real_potential_transfer_function, z_padding + return ( + isotropic_fluorescent_thick_3d.apply_transfer_function( + zyx_object, + real_potential_transfer_function, + z_padding, + background=0, + ) + * brightness + + brightness ) @@ -134,8 +137,6 @@ def apply_inverse_transfer_function( real_potential_transfer_function: Tensor, imaginary_potential_transfer_function: Tensor, z_padding: int, - z_pixel_size: float, # TODO: MOVE THIS PARAM TO OTF? (leaky param) - wavelength_illumination: float, # TOOD: MOVE THIS PARAM TO OTF? (leaky param) absorption_ratio: float = 0.0, reconstruction_algorithm: Literal["Tikhonov", "TV"] = "Tikhonov", regularization_strength: float = 1e-3, @@ -158,14 +159,6 @@ def apply_inverse_transfer_function( z_padding : int Padding for axial dimension. Use zero for defocus stacks that extend ~3 PSF widths beyond the sample. Pad by ~3 PSF widths otherwise. - z_pixel_size : float - spacing between axial samples in sample space - units must be consistent with wavelength_illumination - TODO: move this leaky parameter to calculate_transfer_function - wavelength_illumination : float, - illumination wavelength - units must be consistent with z_pixel_size - TODO: move this leaky parameter to calculate_transfer_function absorption_ratio : float, optional, Absorption-to-phase ratio in the sample. Use default 0 for purely phase objects. @@ -223,4 +216,4 @@ def apply_inverse_transfer_function( if z_padding != 0: f_real = f_real[z_padding:-z_padding] - return f_real * z_pixel_size / 4 / np.pi * wavelength_illumination + return f_real diff --git a/waveorder/optics.py b/waveorder/optics.py index 8033ab86..cc231c77 100644 --- a/waveorder/optics.py +++ b/waveorder/optics.py @@ -270,7 +270,7 @@ def Source_subsample(Source_cont, NAx_coord, NAy_coord, subsampled_NA=0.1): illu_list.append(i) first_idx = False elif ( - np.product( + np.prod( (NAx_list[i] - NAx_list[illu_list]) ** 2 + (NAy_list[i] - NAy_list[illu_list]) ** 2 >= subsampled_NA**2 @@ -739,19 +739,23 @@ def compute_weak_object_transfer_function_3D( H1 = torch.fft.ifft2(torch.conj(SPHz_hat) * PG_hat, dim=(1, 2)) H1 = H1 * window[:, None, None] - H1 = torch.fft.fft(H1, dim=0) * z_pixel_size + H1 = torch.fft.fft(H1, dim=0) H2 = torch.fft.ifft2(SPHz_hat * torch.conj(PG_hat), dim=(1, 2)) H2 = H2 * window[:, None, None] - H2 = torch.fft.fft(H2, dim=0) * z_pixel_size + H2 = torch.fft.fft(H2, dim=0) - I_norm = torch.sum( + direct_intensity = torch.sum( illumination_pupil_support * detection_pupil * torch.conj(detection_pupil) ) - real_potential_transfer_function = (H1 + H2) / I_norm - imag_potential_transfer_function = 1j * (H1 - H2) / I_norm + real_potential_transfer_function = (H1 + H2) / direct_intensity + imag_potential_transfer_function = 1j * (H1 - H2) / direct_intensity + + # Discretization factor for unitless input and output + real_potential_transfer_function *= z_pixel_size + imag_potential_transfer_function *= z_pixel_size return real_potential_transfer_function, imag_potential_transfer_function From 3573180a1d415e2fc8f069a58540ecd8dbfea224 Mon Sep 17 00:00:00 2001 From: Ziwen Liu <67518483+ziw-liu@users.noreply.github.com> Date: Thu, 11 Jul 2024 17:12:24 -0700 Subject: [PATCH 4/4] Pin numpy <2 (#173) pin numpy <2 https://github.com/pytorch/pytorch/issues/128860 --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 1a57bb79..30dcade1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -45,7 +45,7 @@ classifiers = [ "Operating System :: MacOS", ] dependencies = [ - "numpy>=1.21", + "numpy>=1.21, <2", "matplotlib>=3.1.1", "scipy>=1.3.0", "pywavelets>=1.1.1",