From 093c87ea2bca2689f192bdc72c166bfe3d31e3f5 Mon Sep 17 00:00:00 2001 From: alisterburt Date: Wed, 6 Nov 2024 16:26:38 -0800 Subject: [PATCH] updates including tests (#6) * updates including tests * add scipy test dep and make 3.10+ * fix docstring typos * 3.10+ in CI * fix tests --- .github/workflows/ci.yml | 2 +- pyproject.toml | 12 +- src/torch_fourier_slice/backproject.py | 2 +- src/torch_fourier_slice/grids/__init__.py | 2 +- ..._grid.py => central_slice_fftfreq_grid.py} | 4 +- src/torch_fourier_slice/grids/fftfreq_grid.py | 158 ------------------ src/torch_fourier_slice/project.py | 4 +- .../_extract_central_slices_rfft_3d.py | 2 +- .../_insert_central_slices_rfft_3d.py | 2 +- tests/conftest.py | 10 ++ tests/test_torch_fourier_slice.py | 61 ++++++- 11 files changed, 85 insertions(+), 174 deletions(-) rename src/torch_fourier_slice/grids/{central_slice_grid.py => central_slice_fftfreq_grid.py} (92%) delete mode 100644 src/torch_fourier_slice/grids/fftfreq_grid.py create mode 100644 tests/conftest.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 883220a..5a67bc4 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -32,7 +32,7 @@ jobs: strategy: fail-fast: false matrix: - python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"] + python-version: ["3.10", "3.11", "3.12"] platform: [ubuntu-latest, ] # macos-latest, windows-latest] steps: diff --git a/pyproject.toml b/pyproject.toml index fe287df..596114d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -19,7 +19,7 @@ name = "torch-fourier-slice" dynamic = ["version"] description = "Fourier slice extraction/insertion in PyTorch." readme = "README.md" -requires-python = ">=3.8" +requires-python = ">=3.10" license = { text = "BSD-3-Clause" } authors = [{ name = "Alister Burt", email = "alisterburt@gmail.com" }] # https://pypi.org/classifiers/ @@ -27,8 +27,6 @@ classifiers = [ "Development Status :: 3 - Alpha", "License :: OSI Approved :: BSD License", "Programming Language :: Python :: 3", - "Programming Language :: Python :: 3.8", - "Programming Language :: Python :: 3.9", "Programming Language :: Python :: 3.10", "Programming Language :: Python :: 3.11", "Programming Language :: Python :: 3.12", @@ -40,13 +38,19 @@ dependencies = [ "numpy", "einops", "torch_image_lerp", + "torch_grid_utils", ] # https://peps.python.org/pep-0621/#dependencies-optional-dependencies # "extras" (e.g. for `pip install .[test]`) [project.optional-dependencies] # add dependencies used for testing here -test = ["pytest", "pytest-cov"] +test = [ + "pytest", + "pytest-cov", + "torch-fourier-shell-correlation", + "scipy" +] # add anything else you like to have in your dev environment here dev = [ "ipython", diff --git a/src/torch_fourier_slice/backproject.py b/src/torch_fourier_slice/backproject.py index efed8c7..fbc8ae9 100644 --- a/src/torch_fourier_slice/backproject.py +++ b/src/torch_fourier_slice/backproject.py @@ -1,7 +1,7 @@ import torch import torch.nn.functional as F +from torch_grid_utils import fftfreq_grid -from .grids import fftfreq_grid from .slice_insertion import insert_central_slices_rfft_3d diff --git a/src/torch_fourier_slice/grids/__init__.py b/src/torch_fourier_slice/grids/__init__.py index 70a89bb..74d3492 100644 --- a/src/torch_fourier_slice/grids/__init__.py +++ b/src/torch_fourier_slice/grids/__init__.py @@ -1 +1 @@ -from .fftfreq_grid import fftfreq_grid +from .central_slice_fftfreq_grid import central_slice_fftfreq_grid \ No newline at end of file diff --git a/src/torch_fourier_slice/grids/central_slice_grid.py b/src/torch_fourier_slice/grids/central_slice_fftfreq_grid.py similarity index 92% rename from src/torch_fourier_slice/grids/central_slice_grid.py rename to src/torch_fourier_slice/grids/central_slice_fftfreq_grid.py index 3b6f2c1..b09d327 100644 --- a/src/torch_fourier_slice/grids/central_slice_grid.py +++ b/src/torch_fourier_slice/grids/central_slice_fftfreq_grid.py @@ -1,7 +1,7 @@ import einops import torch -from .fftfreq_grid import _construct_fftfreq_grid_2d +from torch_grid_utils import fftfreq_grid from ..dft_utils import rfft_shape, fftshift_2d @@ -13,7 +13,7 @@ def central_slice_fftfreq_grid( ) -> torch.Tensor: # generate 2d grid of DFT sample frequencies, shape (h, w, 2) h, w = volume_shape[-2:] - grid = _construct_fftfreq_grid_2d( + grid = fftfreq_grid( image_shape=(h, w), rfft=rfft, device=device diff --git a/src/torch_fourier_slice/grids/fftfreq_grid.py b/src/torch_fourier_slice/grids/fftfreq_grid.py deleted file mode 100644 index cf55da4..0000000 --- a/src/torch_fourier_slice/grids/fftfreq_grid.py +++ /dev/null @@ -1,158 +0,0 @@ -import functools -from typing import Sequence, Tuple - -import einops -import torch - -from ..dft_utils import rfft_shape, fftshift_2d, fftshift_3d - - -@functools.lru_cache(maxsize=1) -def fftfreq_grid( - image_shape: tuple[int, int] | tuple[int, int, int], - rfft: bool, - fftshift: bool = False, - spacing: float | tuple[float, float] | tuple[float, float, float] = 1, - norm: bool = False, - device: torch.device | None = None, -): - """Construct a 2D or 3D grid of DFT sample frequencies. - - For a 2D image with shape `(h, w)` and `rfft=False` this function will produce - a `(h, w, 2)` array of DFT sample frequencies in the `h` and `w` dimensions. - If `norm` is True the Euclidean norm will be calculated over the last dimension - leaving a `(h, w)` grid. - - Parameters - ---------- - image_shape: tuple[int, int] | tuple[int, int, int] - Shape of the 2D or 3D image before computing the DFT. - rfft: bool - Whether the output should contain frequencies for a real-valued DFT. - fftshift: bool - Whether to fftshift the output grid. - spacing: float | tuple[float, float] | tuple[float, float, float] - Spacing between samples in each dimension. Sampling is considered to be - isotropic if a single value is passed. - norm: bool - Whether to compute the Euclidean norm over the last dimension. - device: torch.device | None - PyTorch device on which the returned grid will be stored. - - Returns - ------- - frequency_grid: torch.Tensor - `(*image_shape, ndim)` array of DFT sample frequencies in each - image dimension if `norm` is `False` else `(*image_shape, )`. - """ - if len(image_shape) == 2: - frequency_grid = _construct_fftfreq_grid_2d( - image_shape=image_shape, - rfft=rfft, - spacing=spacing, - device=device, - ) - if fftshift is True: - frequency_grid = einops.rearrange(frequency_grid, '... freq -> freq ...') - frequency_grid = fftshift_2d(frequency_grid, rfft=rfft) - frequency_grid = einops.rearrange(frequency_grid, 'freq ... -> ... freq') - elif len(image_shape) == 3: - frequency_grid = _construct_fftfreq_grid_3d( - image_shape=image_shape, - rfft=rfft, - spacing=spacing, - device=device, - ) - if fftshift is True: - frequency_grid = einops.rearrange(frequency_grid, '... freq -> freq ...') - frequency_grid = fftshift_3d(frequency_grid, rfft=rfft) - frequency_grid = einops.rearrange(frequency_grid, 'freq ... -> ... freq') - else: - raise NotImplementedError( - "Construction of fftfreq grids is currently only supported for " - "2D and 3D images." - ) - if norm is True: - frequency_grid = einops.reduce( - frequency_grid ** 2, '... d -> ...', reduction='sum' - ) ** 0.5 - return frequency_grid - - -def _construct_fftfreq_grid_2d( - image_shape: Tuple[int, int], - rfft: bool, - spacing: float | tuple[float, float] = 1, - device: torch.device = None -) -> torch.Tensor: - """Construct a grid of DFT sample freqs for a 2D image. - - Parameters - ---------- - image_shape: Sequence[int] - A 2D shape `(h, w)` of the input image for which a grid of DFT sample freqs - should be calculated. - rfft: bool - Whether the frequency grid is for a real fft (rfft). - spacing: float | Tuple[float, float] - Sample spacing in `h` and `w` dimensions of the grid. - device: torch.device - Torch device for the resulting grid. - - Returns - ------- - frequency_grid: torch.Tensor - `(h, w, 2)` array of DFT sample freqs. - Order of freqs in the last dimension corresponds to the order of - the two dimensions of the grid. - """ - dh, dw = spacing if isinstance(spacing, Sequence) else [spacing] * 2 - last_axis_fftfreq = torch.fft.rfftfreq if rfft is True else torch.fft.fftfreq - h, w = image_shape - freq_y = torch.fft.fftfreq(h, d=dh, device=device) - freq_x = last_axis_fftfreq(w, d=dw, device=device) - h, w = rfft_shape(image_shape) if rfft is True else image_shape - freq_yy = einops.repeat(freq_y, 'h -> h w', w=w) - freq_xx = einops.repeat(freq_x, 'w -> h w', h=h) - return einops.rearrange([freq_yy, freq_xx], 'freq h w -> h w freq') - - -def _construct_fftfreq_grid_3d( - image_shape: Sequence[int], - rfft: bool, - spacing: float | Tuple[float, float, float] = 1, - device: torch.device = None -) -> torch.Tensor: - """Construct a grid of DFT sample freqs for a 3D image. - - Parameters - ---------- - image_shape: Sequence[int] - A 3D shape `(d, h, w)` of the input image for which a grid of DFT sample freqs - should be calculated. - rfft: bool - Controls Whether the frequency grid is for a real fft (rfft). - spacing: float | Tuple[float, float, float] - Sample spacing in `d`, `h` and `w` dimensions of the grid. - device: torch.device - Torch device for the resulting grid. - - Returns - ------- - frequency_grid: torch.Tensor - `(h, w, 3)` array of DFT sample freqs. - Order of freqs in the last dimension corresponds to the order of dimensions - of the grid. - """ - dd, dh, dw = spacing if isinstance(spacing, Sequence) else [spacing] * 3 - last_axis_frequency_func = torch.fft.rfftfreq if rfft is True else torch.fft.fftfreq - d, h, w = image_shape - freq_z = torch.fft.fftfreq(d, d=dd, device=device) - freq_y = torch.fft.fftfreq(h, d=dh, device=device) - freq_x = last_axis_frequency_func(w, d=dw, device=device) - d, h, w = rfft_shape(image_shape) if rfft is True else image_shape - freq_zz = einops.repeat(freq_z, 'd -> d h w', h=h, w=w) - freq_yy = einops.repeat(freq_y, 'h -> d h w', d=d, w=w) - freq_xx = einops.repeat(freq_x, 'w -> d h w', d=d, h=h) - return einops.rearrange([freq_zz, freq_yy, freq_xx], 'freq ... -> ... freq') - diff --git a/src/torch_fourier_slice/project.py b/src/torch_fourier_slice/project.py index 1072234..448ca4b 100644 --- a/src/torch_fourier_slice/project.py +++ b/src/torch_fourier_slice/project.py @@ -1,7 +1,7 @@ import torch import torch.nn.functional as F +from torch_grid_utils import fftfreq_grid -from .grids import fftfreq_grid from .slice_extraction import extract_central_slices_rfft_3d @@ -18,7 +18,7 @@ def project_3d_to_2d( volume: torch.Tensor `(d, d, d)` volume. rotation_matrices: torch.Tensor - `(..., 3, 3)` array of rotation matrices for insert of `images`. + `(..., 3, 3)` array of rotation matrices for insertion of `images`. Rotation matrices left-multiply column vectors containing xyz coordinates. pad: bool Whether to pad the volume 2x with zeros to increase sampling rate in the DFT. diff --git a/src/torch_fourier_slice/slice_extraction/_extract_central_slices_rfft_3d.py b/src/torch_fourier_slice/slice_extraction/_extract_central_slices_rfft_3d.py index d216aac..8987443 100644 --- a/src/torch_fourier_slice/slice_extraction/_extract_central_slices_rfft_3d.py +++ b/src/torch_fourier_slice/slice_extraction/_extract_central_slices_rfft_3d.py @@ -3,7 +3,7 @@ from torch_image_lerp import sample_image_3d from ..dft_utils import fftfreq_to_dft_coordinates -from ..grids.central_slice_grid import central_slice_fftfreq_grid +from ..grids.central_slice_fftfreq_grid import central_slice_fftfreq_grid def extract_central_slices_rfft_3d( diff --git a/src/torch_fourier_slice/slice_insertion/_insert_central_slices_rfft_3d.py b/src/torch_fourier_slice/slice_insertion/_insert_central_slices_rfft_3d.py index 7b09840..cf02b21 100644 --- a/src/torch_fourier_slice/slice_insertion/_insert_central_slices_rfft_3d.py +++ b/src/torch_fourier_slice/slice_insertion/_insert_central_slices_rfft_3d.py @@ -3,7 +3,7 @@ from torch_image_lerp import insert_into_image_3d from ..dft_utils import fftfreq_to_dft_coordinates, rfft_shape -from ..grids.central_slice_grid import central_slice_fftfreq_grid +from ..grids.central_slice_fftfreq_grid import central_slice_fftfreq_grid def insert_central_slices_rfft_3d( diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..6ecef69 --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,10 @@ +import torch +from pytest import fixture + + +@fixture +def cube() -> torch.Tensor: + volume = torch.zeros((32, 32, 32)) + volume[8:24, 8:24, 8:24] = 1 + volume[16, 16, 16] = 32 + return volume diff --git a/tests/test_torch_fourier_slice.py b/tests/test_torch_fourier_slice.py index 4c6e45f..bc3ef33 100644 --- a/tests/test_torch_fourier_slice.py +++ b/tests/test_torch_fourier_slice.py @@ -1,4 +1,59 @@ -# temporary stub +import pytest +import torch -def test_something(): - pass +from torch_fourier_slice import project_3d_to_2d, backproject_2d_to_3d +from torch_fourier_shell_correlation import fsc +from scipy.stats import special_ortho_group + + +def test_project_3d_to_2d_rotation_center(): + # rotation center should be at position of DC in DFT + volume = torch.zeros((32, 32, 32)) + volume[16, 16, 16] = 1 + + # make projections + rotation_matrices = torch.tensor(special_ortho_group.rvs(dim=3, size=100)).float() + projections = project_3d_to_2d( + volume=volume, + rotation_matrices=rotation_matrices, + ) + + # check max is always at (16, 16), implying point (16, 16) never moves + for image in projections: + max = torch.argmax(image) + i, j = divmod(max.item(), 32) + assert (i, j) == (16, 16) + + +def test_3d_2d_projection_backprojection_cycle(cube): + # make projections + rotation_matrices = torch.tensor(special_ortho_group.rvs(dim=3, size=1500)).float() + projections = project_3d_to_2d( + volume=cube, + rotation_matrices=rotation_matrices, + ) + + # reconstruct + volume = backproject_2d_to_3d( + images=projections, + rotation_matrices=rotation_matrices, + ) + + # calculate FSC between the projections and the reconstructions + _fsc = fsc(cube, volume.float()) + + assert torch.all(_fsc[-5:] > 0.99) # few low res shells at 0.98... + + +@pytest.mark.parametrize( + "images, rotation_matrices", + [ + ( + torch.rand((10, 28, 28)).float(), + torch.tensor(special_ortho_group.rvs(dim=3, size=10)).float() + ), + ] +) +def test_dtypes_slice_insertion(images, rotation_matrices): + result = backproject_2d_to_3d(images, rotation_matrices) + assert result.dtype == torch.float64