From 78310f8915feb34a53ba06126e58945296c7d955 Mon Sep 17 00:00:00 2001 From: Roy Hoitink Date: Mon, 19 Feb 2024 13:13:21 +0100 Subject: [PATCH] Fix bug in noise generation --- pyproject.toml | 1 + src/simulatedmicroscopy/image.py | 97 +++++++++++++++++++++----------- tests/test_image.py | 38 ++++++++++--- 3 files changed, 96 insertions(+), 40 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 0932d25..4e0a681 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -28,6 +28,7 @@ select = [ ignore = [ "E501", "E203", + "UP007", ] exclude = [ diff --git a/src/simulatedmicroscopy/image.py b/src/simulatedmicroscopy/image.py index 2bf4c96..e145725 100644 --- a/src/simulatedmicroscopy/image.py +++ b/src/simulatedmicroscopy/image.py @@ -1,13 +1,16 @@ from __future__ import annotations -import numpy as np -from typing import Optional, Union, Tuple, List + +import warnings from pathlib import Path +from typing import Optional + import h5py +import numpy as np import scipy.signal +import skimage.measure + from .input import Coordinates from .particle import BaseParticle -import warnings -import skimage.measure class Image: @@ -78,7 +81,7 @@ def _number_of_dimensions(self) -> int: return len(self.image.shape) def get_pixel_sizes( - self, dimensions: Optional[list[str]] = list("zyx"), unit: str = "m" + self, dimensions: Optional[list[str]] = None, unit: str = "m" ) -> list[float]: """Get pixel sizes for given dimensions in given unit @@ -94,6 +97,9 @@ def get_pixel_sizes( list[float] List of pixel sizes in given unit, in order of requested dimensions. """ + if dimensions is None: + dimensions = ["zyx"] + if self._number_of_dimensions() < 3 and "z" in dimensions: raise ValueError("Cannot retrieve z dimension for non 3D image") @@ -114,7 +120,9 @@ def get_pixel_sizes( return ps - def save_h5file(self, filename: str, description: Optional[str] = None) -> bool: + def save_h5file( + self, filename: str, description: Optional[str] = None + ) -> bool: """Save image as h5 file (custom format, not compatible with Huygens) Parameters @@ -138,9 +146,9 @@ def save_h5file(self, filename: str, description: Optional[str] = None) -> bool: f.create_group("Metadata") for dim in self.DIMENSIONS_ORDER: - f[f"Metadata/DimensionScale{dim.upper()}"] = self.get_pixel_sizes( - [dim] - )[0] + f[ + f"Metadata/DimensionScale{dim.upper()}" + ] = self.get_pixel_sizes([dim])[0] if self.metadata is not None: for k, v in self.metadata.items(): @@ -173,7 +181,7 @@ def load_h5file(cls, filename: str) -> type[Image]: float(f[f"Metadata/DimensionScale{dim.upper()}"][()]) for dim in list("zyx") ] - if "PixelCoordinates" in f["Metadata"].keys(): + if "PixelCoordinates" in f["Metadata"]: pixel_coordinates = f["Metadata/PixelCoordinates"][()] else: pixel_coordinates = None @@ -231,7 +239,11 @@ def _get_point_image_array( @classmethod def create_point_image( - cls, coordinates: type[Coordinates], pixel_sizes: list[float], *args, **kwargs + cls, + coordinates: type[Coordinates], + pixel_sizes: list[float], + *args, + **kwargs, ) -> type[Image]: """Create point source image in which every point from the set of coordinates is represented by a single white pixel @@ -247,9 +259,11 @@ def create_point_image( type[Image] Genereated image """ - (zs, ys, xs), image = cls._get_point_image_array(coordinates, pixel_sizes) + (zs, ys, xs), image = cls._get_point_image_array( + coordinates, pixel_sizes + ) - im = cls(image=image, pixel_sizes=pixel_sizes, *args, **kwargs) + im = cls(image, pixel_sizes, *args, **kwargs) im.pixel_coordinates = np.transpose([zs, ys, xs]) return im @@ -258,7 +272,7 @@ def create_particle_image( cls, coordinates: type[Coordinates], particle: type[BaseParticle], - edge_pixel_margin: Union[int, List[int]] = 0, + edge_pixel_margin: int | list[int] = 0, *args, **kwargs, ) -> type[Image]: @@ -270,7 +284,7 @@ def create_particle_image( Set of coordinates particle : list[BaseParticle] Particle to use for image, will also use its pixel size for the final image - edge_pixel_margin : Union[int,List[int]] + edge_pixel_margin : Union[int,list[int]] Number of extra empty pixels to apply as a margin around the edges of the image. Either single int for all directions or one int for every dimension (in z,y,x order). @@ -281,7 +295,7 @@ def create_particle_image( """ # convert pixel sizes to micrometers for calculatation pixel_sizes_um = np.array(particle.pixel_sizes) * 1e6 - if isinstance(edge_pixel_margin, int) or isinstance(edge_pixel_margin, float): + if isinstance(edge_pixel_margin, (int, float)): # not iterable, repeat for all dims edge_pixel_margin = np.array( [edge_pixel_margin] * len(pixel_sizes_um), dtype=int @@ -304,7 +318,8 @@ def create_particle_image( # round to integer to create point at certain pixel xs, ys, zs = ( - np.round(scaled_coords).astype(int) + edge_pixel_margin[::-1, np.newaxis] + np.round(scaled_coords).astype(int) + + edge_pixel_margin[::-1, np.newaxis] ) # edge_pixel_margin to x,y,z order by reversing @@ -318,9 +333,11 @@ def create_particle_image( image = np.zeros(shape=image_size) for x, y, z in zip(xs, ys, zs): - image[z : z + prs[0], y : y + prs[1], x : x + prs[2]] += particle_response + image[ + z : z + prs[0], y : y + prs[1], x : x + prs[2] + ] += particle_response - im = cls(image=image, pixel_sizes=particle.pixel_sizes, *args, **kwargs) + im = cls(image, particle.pixel_sizes, *args, **kwargs) im.pixel_coordinates = np.transpose([zs, ys, xs]) + particle_offset.T return im @@ -347,7 +364,9 @@ def downsample(self, downsample_factor: list[int]) -> type[Image]: ds_fac = np.round(downsample_factor).astype(int) # reduce image size by taking the mean of a AxBxC (defined by downsample_factor) sized block and use that as new pixel value - self.image = skimage.measure.block_reduce(self.image, tuple(ds_fac), np.mean) + self.image = skimage.measure.block_reduce( + self.image, tuple(ds_fac), np.mean + ) # adapt pixel sizes self.pixel_sizes = self.pixel_sizes * ds_fac @@ -375,9 +394,13 @@ def convolve(self, other: type[Image]) -> type[Image]: The convolved image """ if not np.isclose(self.pixel_sizes, other.pixel_sizes).all(): - raise ValueError("Cannot convolve images with different pixel sizes") + raise ValueError( + "Cannot convolve images with different pixel sizes" + ) - self.image = scipy.signal.convolve(self.image, other.image, mode="same") + self.image = scipy.signal.convolve( + self.image, other.image, mode="same" + ) self.is_convolved = True self.metadata["is_convolved"] = True @@ -417,10 +440,13 @@ def add_shot_noise(self, SNR: float = 30.0) -> type[Image]: The image with shot noise """ if self.has_shot_noise: - warnings.warn("This image already has shot noise, proceeding anyway") + warnings.warn( + "This image already has shot noise, proceeding anyway", + stacklevel=2, + ) # Poisson cannot deal with negative values, if these are present, shift whole image - if np.any(self.image < 0.0): + if np.any(self.image < 0.0): self.image += np.abs(self.image.min()) scaling_factor = ( @@ -428,7 +454,7 @@ def add_shot_noise(self, SNR: float = 30.0) -> type[Image]: ) # scaling factor to get a proper Poisson sampling self.image = ( - np.random.poisson(scaling_factor * self.image.mean()) / scaling_factor + np.random.poisson(scaling_factor * self.image) / scaling_factor ) # division by `scaling_factor` at the end to revert pixel values back to their (roughly) original range @@ -436,7 +462,9 @@ def add_shot_noise(self, SNR: float = 30.0) -> type[Image]: return self - def add_read_noise(self, SNR: float = 50.0, background: float = 0.0) -> type[Image]: + def add_read_noise( + self, SNR: float = 50.0, background: float = 0.0 + ) -> type[Image]: """Add read noise to the image, realised by adding noise sampled from a normal distribution. The standard deviation of the distribution is constant across the entire image and determined by the SNR parameter, it is therefore signal-independent. @@ -454,7 +482,10 @@ def add_read_noise(self, SNR: float = 50.0, background: float = 0.0) -> type[Ima The image with read noise """ if self.has_read_noise: - warnings.warn("This image already has read noise, proceeding anyway") + warnings.warn( + "This image already has read noise, proceeding anyway", + stacklevel=2, + ) peak = np.percentile( self.image, 99.9 @@ -464,7 +495,7 @@ def add_read_noise(self, SNR: float = 50.0, background: float = 0.0) -> type[Ima sigma = peak / SNR self.image += np.random.normal( - loc=background, scale=sigma + loc=background, scale=sigma, size=self.image.shape ) # centred around value of `background` self.has_read_noise = True # set read noise flag @@ -482,7 +513,7 @@ def get_pixel_coordinates(self) -> np.ndarray: return self.pixel_coordinates def get_particle_image( - self, particle_id: int = 0, box_size: Tuple[int] = (5, 5, 5) + self, particle_id: int = 0, box_size: tuple[int] = (5, 5, 5) ) -> np.ndarray: """Get 3D box around the coordinates of particle with index `particle_id` @@ -490,7 +521,7 @@ def get_particle_image( ---------- particle_id : int, optional Index of the particle, by default 0 - box_size : Tuple[int], optional + box_size : tuple[int], optional Size of the surrounding box in pixels, by default (5, 5, 5) Returns @@ -503,7 +534,9 @@ def get_particle_image( "Not a image genereted from particle, cannot extract single particles" ) - if int(particle_id) < 0 or int(particle_id) >= len(self.pixel_coordinates): + if int(particle_id) < 0 or int(particle_id) >= len( + self.pixel_coordinates + ): raise IndexError("Invalid particle index") coords = self.pixel_coordinates[int(particle_id)].astype(int) @@ -518,7 +551,7 @@ def get_particle_image( class HuygensImage(Image): - def __init__(self, filename: Union[str, Path]) -> None: + def __init__(self, filename: str | Path) -> None: """Wrapper for Huygens-generated .h5 images Extends the `Image` class. diff --git a/tests/test_image.py b/tests/test_image.py index 40d4c78..6c429be 100644 --- a/tests/test_image.py +++ b/tests/test_image.py @@ -1,6 +1,6 @@ -from simulatedmicroscopy import Image, HuygensImage import numpy as np import pytest +from simulatedmicroscopy import HuygensImage, Image def create_demo_image(): @@ -56,7 +56,9 @@ def test_huygens_notexisting(tmp_path): def test_pixel_size_conversion(unit, conversionfactor): pixel_sizes = np.array([1e-9, 1e-8, 1e-7]) image = Image(np.zeros(shape=(5, 5, 5)), pixel_sizes) - assert (image.get_pixel_sizes(unit=unit) == pixel_sizes * conversionfactor).all() + assert ( + image.get_pixel_sizes(unit=unit) == pixel_sizes * conversionfactor + ).all() def test_image_equality_image(): @@ -113,7 +115,8 @@ def test_downsample(downsample_factor): # check if the shape of the image has decreased by `downsample_factor` assert ( - downsampled.image.shape == np.array(im_array.shape) // downsample_factor + downsampled.image.shape + == np.array(im_array.shape) // downsample_factor ).all() # check if the pixel_size of the image has increased by `downsample_factor` @@ -159,11 +162,12 @@ def test_noise_deprecated(): im = create_demo_image() im.noisify() + def test_shot_noise(): im = create_demo_image() # check if returns correct image - assert im.add_shot_noise(SNR = 10.0) != create_demo_image() + assert im.add_shot_noise(SNR=10.0) != create_demo_image() # check if original image was also changed assert im != create_demo_image() @@ -171,11 +175,12 @@ def test_shot_noise(): # should be set to True assert im.has_shot_noise + def test_read_noise(): im = create_demo_image() # check if returns correct image - assert im.add_read_noise(SNR = 10.0, background = 1e-3) != create_demo_image() + assert im.add_read_noise(SNR=10.0, background=1e-3) != create_demo_image() # check if original image was also changed assert im != create_demo_image() @@ -183,6 +188,17 @@ def test_read_noise(): # should be set to True assert im.has_read_noise + +def test_noise_apply_both(): + im = create_demo_image() + + assert im.add_read_noise().add_shot_noise() != create_demo_image() + + assert im.has_read_noise + + assert im.has_shot_noise + + def test_point_image(tmp_path): from simulatedmicroscopy import Coordinates @@ -234,7 +250,9 @@ def test_pixel_coordinates_after_downscale(downsample_factor): im.downsample([downsample_factor] * 3) - assert (pixel_coords_before == im.get_pixel_coordinates() * downsample_factor).all() + assert ( + pixel_coords_before == im.get_pixel_coordinates() * downsample_factor + ).all() def test_pixel_coordinates_after_downscale_onlyz(): @@ -263,7 +281,9 @@ def test_pixel_coordinates_after_downscale_onlyz(): def test_image_metadata_wrongtype(): with pytest.raises(ValueError): - Image(np.zeros(shape=(5, 5, 5)), [1e-6, 1e-6, 1e-6], metadata=[1, 2, 3]) + Image( + np.zeros(shape=(5, 5, 5)), [1e-6, 1e-6, 1e-6], metadata=[1, 2, 3] + ) def test_image_metadata_is_set(): @@ -319,7 +339,9 @@ def test_particle_image_desired_shape(): ] ) - particle_im = Image.create_particle_image(cs, Sphere([1e-6, 1e-6, 1e-6], 1e-6)) + particle_im = Image.create_particle_image( + cs, Sphere([1e-6, 1e-6, 1e-6], 1e-6) + ) box_shape = (2, 2, 2)