diff --git a/docs/source/optical_and_acoustic_simulation_US_heterogeneity.rst b/docs/source/optical_and_acoustic_simulation_US_heterogeneity.rst new file mode 100644 index 00000000..f9828d40 --- /dev/null +++ b/docs/source/optical_and_acoustic_simulation_US_heterogeneity.rst @@ -0,0 +1,7 @@ +optical_and_acoustic_simulation_US_heterogeneity +========================================= + +.. literalinclude:: ../../simpa_examples/optical_and_acoustic_simulation_US_heterogeneity.py + :language: python + :lines: 1- + diff --git a/docs/source/simpa_examples.rst b/docs/source/simpa_examples.rst index 9d229be6..b93836d8 100644 --- a/docs/source/simpa_examples.rst +++ b/docs/source/simpa_examples.rst @@ -12,6 +12,7 @@ simpa\_examples minimal_optical_simulation_uniform_cube msot_invision_simulation optical_and_acoustic_simulation + optical_and_acoustic_simulation_US_heterogeneity perform_image_reconstruction perform_iterative_qPAI_reconstruction segmentation_loader diff --git a/simpa/utils/libraries/heterogeneity_generator.py b/simpa/utils/libraries/heterogeneity_generator.py index 8d01b3f1..71d3b57e 100644 --- a/simpa/utils/libraries/heterogeneity_generator.py +++ b/simpa/utils/libraries/heterogeneity_generator.py @@ -1,14 +1,18 @@ # SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT +from abc import abstractmethod import numpy as np from sklearn.datasets import make_blobs from scipy.ndimage.filters import gaussian_filter from skimage import transform from simpa.utils import Tags -from typing import Union, Optional +from typing import Union from simpa.log import Logger +import os +import requests +import zipfile class HeterogeneityGeneratorBase(object): @@ -42,9 +46,13 @@ def __init__(self, xdim, ydim, zdim, spacing_mm, target_mean=None, self.map = np.ones((self._xdim, self._ydim, self._zdim), dtype=float) + @abstractmethod def get_map(self): - self.normalise_map() - return self.map.astype(float) + """ + A method to return the 3D heterogeneity map. In some cases, this will mean changing from 2D to 3D at this step + :return: 3D heterogeneity map + """ + pass def normalise_map(self): """ @@ -110,6 +118,10 @@ def __init__(self, xdim, ydim, zdim, spacing_mm, gaussian_blur_size_mm=None, tar _gaussian_blur_size_voxels = gaussian_blur_size_mm / spacing_mm self.map = gaussian_filter(self.map, _gaussian_blur_size_voxels) + def get_map(self): + self.normalise_map() + return self.map.astype(float) + class BlobHeterogeneity(HeterogeneityGeneratorBase): """ @@ -150,22 +162,33 @@ def __init__(self, xdim, ydim, zdim, spacing_mm, num_centers=None, cluster_std=N np.percentile(x[:, 2], 95))))[0] self.map = gaussian_filter(self.map, 5) + def get_map(self): + self.normalise_map() + return self.map.astype(float) + class ImageHeterogeneity(HeterogeneityGeneratorBase): """ This heterogeneity generator takes a pre-specified 2D image, currently only supporting numpy arrays, and uses them - as a map for heterogeneity within the tissue. + as a map for heterogeneity within the tissue. By default, it will use download and use beef ultrasound images taken + by our team. + + ########## + This class will (if not previously downloaded in the directory of the simulation) download a folder with beef + ultrasound images + ########## Attributes: map: the np array of the heterogeneity map transformed and augments to fill the area """ - def __init__(self, xdim: int, ydim: int, zdim: int, heterogeneity_image: np.ndarray, - spacing_mm: Union[int, float] = None, image_pixel_spacing_mm: Union[int, float] = None, - scaling_type: [None, str] = None, constant: Union[int, float] = 0, + def __init__(self, xdim: int, ydim: int, zdim: int, spacing_mm: Union[int, float], + heterogeneity_image: np.ndarray = None, image_pixel_spacing_mm: Union[int, float] = None, + scaling_type: str = Tags.IMAGE_SCALING_SYMMETRIC, constant: Union[int, float] = 0, crop_placement=('centre', 'centre'), target_mean: Union[int, float] = None, target_std: Union[int, float] = None, target_min: Union[int, float] = None, - target_max: Union[int, float] = None): + target_max: Union[int, float] = None, beef_ultrasound_database_path: str = None, + ultrasound_image_type: str = Tags.MEAT_ULTRASOUND_CROPPED, scan_number: int = None): """ :param xdim: the x dimension of the volume in voxels :param ydim: the y dimension of the volume in voxels @@ -175,11 +198,11 @@ def __init__(self, xdim: int, ydim: int, zdim: int, heterogeneity_image: np.ndar :param image_pixel_spacing_mm: the pixel spacing of the image in mm (pixel spacing) :param scaling_type: the scaling type of the heterogeneity map, with default being that no scaling occurs OPTIONS: - TAGS.IMAGE_SCALING_SYMMETRIC: symmetric reflections of the image to span the area - TAGS.IMAGE_SCALING_STRETCH: stretch the image to span the area - TAGS.IMAGE_SCALING_WRAP: multiply the image to span the area - TAGS.IMAGE_SCALING_EDGE: continue the values at the edge of the area to fill the shape - TAGS.IMAGE_SCALING_CONSTANT: span the left-over area with a constant + Tags.IMAGE_SCALING_SYMMETRIC: symmetric reflections of the image to span the area + Tags.IMAGE_SCALING_STRETCH: stretch the image to span the area + Tags.IMAGE_SCALING_WRAP: multiply the image to span the area + Tags.IMAGE_SCALING_EDGE: continue the values at the edge of the area to fill the shape + Tags.IMAGE_SCALING_CONSTANT: span the left-over area with a constant :param constant: the scaling constant of the heterogeneity map, used only for scaling type 'constant' WARNING: scaling constant must be in reference to the values in the heterogeneity_image :param crop_placement: the placement of where the heterogeneity map is cropped @@ -190,44 +213,45 @@ def __init__(self, xdim: int, ydim: int, zdim: int, heterogeneity_image: np.ndar """ super().__init__(xdim, ydim, zdim, spacing_mm, target_mean, target_std, target_min, target_max) self.logger = Logger() + self.heterogeneity_image = heterogeneity_image + + if self.heterogeneity_image is None: + self.heterogeneity_image = get_ultrasound_image(beef_ultrasound_database_path=beef_ultrasound_database_path, + image_type=ultrasound_image_type, + scan_number=scan_number) + image_pixel_spacing_mm = 0.2 if image_pixel_spacing_mm is None: image_pixel_spacing_mm = spacing_mm - (image_width_pixels, image_height_pixels) = heterogeneity_image.shape - [image_width_mm, image_height_mm] = np.array([image_width_pixels, image_height_pixels]) * image_pixel_spacing_mm - (xdim_mm, ydim_mm, zdim_mm) = np.array([xdim, ydim, zdim]) * spacing_mm + if scaling_type == Tags.IMAGE_SCALING_STRETCH: + self.heterogeneity_image = transform.resize(self.heterogeneity_image, output_shape=(xdim, zdim), + mode='symmetric') + else: + (image_width_pixels, image_height_pixels) = self.heterogeneity_image.shape + [image_width_mm, image_height_mm] = np.array( + [image_width_pixels, image_height_pixels]) * image_pixel_spacing_mm + (xdim_mm, ydim_mm, zdim_mm) = np.array([xdim, ydim, zdim]) * spacing_mm - wider = image_width_mm > xdim_mm - taller = image_height_mm > zdim_mm + wider = image_width_mm > xdim_mm + taller = image_height_mm > zdim_mm - if taller or wider: - pixel_scaled_image = self.change_resolution(heterogeneity_image, spacing_mm=spacing_mm, - image_pixel_spacing_mm=image_pixel_spacing_mm) - cropped_image = self.crop_image(xdim, zdim, pixel_scaled_image, crop_placement) + if taller or wider: + self.change_resolution(spacing_mm=spacing_mm, image_pixel_spacing_mm=image_pixel_spacing_mm) + self.crop_image(xdim, zdim, crop_placement) + if not taller and not wider: + self.upsize_to_fill_area(xdim, zdim, scaling_type, constant) - if taller and wider: - area_fill_image = cropped_image else: - area_fill_image = self.upsize_to_fill_area(xdim, zdim, cropped_image, scaling_type, constant) - - else: - pixel_scaled_image = self.change_resolution(heterogeneity_image, spacing_mm=spacing_mm, - image_pixel_spacing_mm=image_pixel_spacing_mm) - area_fill_image = self.upsize_to_fill_area(xdim, zdim, pixel_scaled_image, scaling_type, constant) - - if scaling_type == Tags.IMAGE_SCALING_STRETCH: - area_fill_image = transform.resize(heterogeneity_image, output_shape=(xdim, zdim), mode='symmetric') - - self.map = np.repeat(area_fill_image[:, np.newaxis, :], ydim, axis=1) + self.change_resolution(spacing_mm=spacing_mm, image_pixel_spacing_mm=image_pixel_spacing_mm) + self.upsize_to_fill_area(xdim, zdim, scaling_type, constant) - def upsize_to_fill_area(self, xdim: int, zdim: int, image: np.ndarray, scaling_type: Optional[str] = None, - constant: Union[int, float] = 0) -> np.ndarray: + def upsize_to_fill_area(self, xdim: int, zdim: int, scaling_type: str = Tags.IMAGE_SCALING_SYMMETRIC, + constant: Union[int, float] = 0): """ Fills an area with an image through various methods of expansion :param xdim: the x dimension of the area to be filled in voxels :param zdim: the z dimension of the area to be filled in voxels - :param image: the input image :param scaling_type: the scaling type of the heterogeneity map, with default being that no scaling occurs OPTIONS: TAGS.IMAGE_SCALING_SYMMETRIC: symmetric reflections of the image to span the area @@ -236,38 +260,36 @@ def upsize_to_fill_area(self, xdim: int, zdim: int, image: np.ndarray, scaling_t TAGS.IMAGE_SCALING_EDGE: continue the values at the edge of the area to fill the shape TAGS.IMAGE_SCALING_CONSTANT: span the left-over area with a constant :param constant: the scaling constant of the heterogeneity map, used only for scaling type 'constant' - :return:A numpy array of size (xdim, zdim) containing the filled image expanded to fill the shape """ - if scaling_type is None or scaling_type == Tags.IMAGE_SCALING_STRETCH: - scaled_image = image + if scaling_type == Tags.IMAGE_SCALING_STRETCH: + pass elif scaling_type == Tags.IMAGE_SCALING_CONSTANT: - pad_left = int((xdim - len(image)) / 2) - pad_height = int(zdim - len(image[0])) - pad_right = xdim - pad_left - len(image) - scaled_image = np.pad(image, ((pad_left, pad_right), (0, pad_height)), - mode=scaling_type, constant_values=constant) + pad_left = int((xdim - len(self.heterogeneity_image)) / 2) + pad_height = int(zdim - len(self.heterogeneity_image[0])) + pad_right = xdim - pad_left - len(self.heterogeneity_image) + self.heterogeneity_image = np.pad(array=self.heterogeneity_image, + pad_width=((pad_left, pad_right), (0, pad_height)), mode=scaling_type, + constant_values=constant) else: - pad_left = int((xdim - len(image)) / 2) - pad_height = int(zdim - len(image[0])) - pad_right = xdim - pad_left - len(image) - scaled_image = np.pad(image, ((pad_left, pad_right), (0, pad_height)), - mode=scaling_type) + pad_left = int((xdim - len(self.heterogeneity_image)) / 2) + pad_height = int(zdim - len(self.heterogeneity_image[0])) + pad_right = xdim - pad_left - len(self.heterogeneity_image) + self.heterogeneity_image = np.pad(array=self.heterogeneity_image, + pad_width=((pad_left, pad_right), (0, pad_height)), mode=scaling_type) self.logger.warning("The input image has filled the area by using {} scaling type".format(scaling_type)) - return scaled_image - def crop_image(self, xdim: int, zdim: int, image: np.ndarray, - crop_placement: Union[str, tuple] = Tags.CROP_POSITION_CENTRE) -> np.ndarray: + def crop_image(self, xdim: int, zdim: int, + crop_placement: Union[str, tuple] = Tags.CROP_POSITION_CENTRE): """ Crop the image to fit specified dimensions xdim and zdim :param xdim: the x dimension of the area to be filled in voxels :param zdim: the z dimension of the area to be filled in voxels - :param image: the input image :param crop_placement: the placement of where the heterogeneity map is cropped OPTIONS: TAGS.CROP_PLACEMENT_[TOP,BOTTOM,LEFT,RIGHT,CENTRE,RANDOM] or position of left hand corner on image - :return: cropped image + :raises: ValueError for invalid placements """ - (image_width_pixels, image_height_pixels) = image.shape + (image_width_pixels, image_height_pixels) = self.heterogeneity_image.shape crop_width = min(xdim, image_width_pixels) crop_height = min(zdim, image_height_pixels) @@ -280,6 +302,9 @@ def crop_image(self, xdim: int, zdim: int, image: np.ndarray, crop_horizontal = round((image_width_pixels - crop_width) / 2) elif isinstance(crop_placement[0], int): crop_horizontal = crop_placement[0] + else: + raise ValueError(f"Invalid crop placement {crop_placement[0]}. Please check Tags.CROP_POSITION_... for" + f"valid string arguments and that numbers are of type int") if crop_placement[1] == Tags.CROP_POSITION_TOP: crop_vertical = 0 @@ -289,6 +314,9 @@ def crop_image(self, xdim: int, zdim: int, image: np.ndarray, crop_vertical = round((image_height_pixels - crop_height) / 2) elif isinstance(crop_placement[1], int): crop_vertical = crop_placement[1] + else: + raise ValueError(f"Invalid crop placement {crop_placement[1]}. Please check Tags.CROP_POSITION_... for" + f"valid string arguments and that numbers are of type int") elif isinstance(crop_placement, str): if crop_placement == Tags.CROP_POSITION_CENTRE: @@ -301,27 +329,120 @@ def crop_image(self, xdim: int, zdim: int, image: np.ndarray, crop_vertical = image_height_pixels - crop_height if crop_vertical != 0: crop_vertical = np.random.randint(0, crop_vertical) + else: + raise ValueError(f"Invalid crop placement {crop_placement}. Please check Tags.CROP_POSITION_... for" + f"valid arguments") - cropped_image = image[crop_horizontal: crop_horizontal + crop_width, crop_vertical: crop_vertical + crop_height] + else: + raise ValueError("Crop placement must be tuple or str") + + self.heterogeneity_image = self.heterogeneity_image[crop_horizontal: crop_horizontal + + crop_width, crop_vertical: crop_vertical + crop_height] self.logger.warning( "The input image has been cropped to the dimensions of the simulation volume ({} {})".format(xdim, zdim)) - return cropped_image - def change_resolution(self, image: np.ndarray, spacing_mm: Union[int, float], - image_pixel_spacing_mm: Union[int, float]) -> np.ndarray: + def change_resolution(self, spacing_mm: Union[int, float], + image_pixel_spacing_mm: Union[int, float]): """ Method to change the resolution of an image - :param image: input image :param image_pixel_spacing_mm: original image pixel spacing mm :param spacing_mm: target pixel spacing mm - :return: image with new pixel spacing """ - (image_width_pixels, image_height_pixels) = image.shape + (image_width_pixels, image_height_pixels) = self.heterogeneity_image.shape [image_width_mm, image_height_mm] = np.array([image_width_pixels, image_height_pixels]) * image_pixel_spacing_mm new_image_pixel_width = round(image_width_mm / spacing_mm) new_image_pixel_height = round(image_height_mm / spacing_mm) self.logger.warning( "The input image has changed pixel spacing to {} to match the simulation volume".format(spacing_mm)) - return transform.resize(image, (new_image_pixel_width, new_image_pixel_height)) + self.heterogeneity_image = transform.resize(self.heterogeneity_image, (new_image_pixel_width, + new_image_pixel_height)) + + def exponential(self, factor: Union[int, float] = 6): + """ + Method to put an exponential weighting on the image. This is implemented as we believe the images created the + MSOT Acuity Echo device might be exponential, and hence this method will reverse this. + :param factor: The exponential factor + """ + self.heterogeneity_image = np.exp(factor * self.heterogeneity_image / np.max(self.heterogeneity_image)) + + def invert_image(self): + """ + Method to invert the image + """ + self.heterogeneity_image = np.max(self.heterogeneity_image) - self.heterogeneity_image + + def get_map(self): + self.map = np.repeat(self.heterogeneity_image[:, np.newaxis, :], self._ydim, axis=1) + self.normalise_map() + return self.map.astype(float) + + +def download_ultrasound_images(save_dir: str): + """ + Downloads the latest beef ultrasound images from nextcloud. The metadata about these images can be found in the + folder + + :param save_dir: directory to save the images to + :return: None + """ + logger = Logger() + # nextcloud url with the reference images + nextcloud_url = "https://hub.dkfz.de/s/g8fLZiY5D6ZC3Hx" # shared "beef_ultrasound_database" folder on nextcloud + # Specify the local directory to save the files + zip_filepath = os.path.join(save_dir, "downloaded.zip") + # Construct the download URL based on the public share link + download_url = nextcloud_url.replace('/s/', '/index.php/s/') + '/download' + # Send a GET request to download the file + logger.debug(f'Download folder with ultrasound figures from nextcloud...') + response = requests.get(download_url) + if response.status_code == 200: + # Save the file + with open(zip_filepath, 'wb') as f: + f.write(response.content) + logger.debug(f'File downloaded successfully and stored at {zip_filepath}.') + else: + logger.critical(f'Failed to download file. Status code: {response.status_code}') + raise requests.exceptions.HTTPError(f'Failed to download file. Status code: {response.status_code}') + + # Open the zip file + with zipfile.ZipFile(zip_filepath, 'r') as zip_ref: + # Extract all the contents into the specified directory + zip_ref.extractall(save_dir) + + logger.debug(f'Files extracted to {save_dir}') + + # Remove the zip file after extraction + os.remove(zip_filepath) + logger.debug(f'{zip_filepath} removed successfully.') + + +def get_ultrasound_image(beef_ultrasound_database_path: str = None, image_type: str = Tags.MEAT_ULTRASOUND_CROPPED, + scan_number: int = None): + """ + A method to retrieve an ultrasound image from the beef ultrasound database. If the ultrasound database has not + been downloaded in the current working directory and beef_ultrasound_database_path is not given, the database + will be downloaded in the current working directory. To avoíd accidentally downloading twice, after the first + occasion of using this method, set the beef_ultrasound_database_path parameter to point to the downloaded + folder, which will ensure you won't accidentally download again when working in other directories. + + :param beef_ultrasound_database_path: the path to the beef ultrasound database + :param image_type: whether you would like to use the regular or cropped images + :param scan_number: the scan number of the beef ultrasound database you wish to use. If not set, a random scan will be chosen + :return: the ultrasound image + """ + logger = Logger() + if not beef_ultrasound_database_path: + current_dir = os.getcwd() + beef_ultrasound_database_path = os.path.join(current_dir, "beef_ultrasound_database") + if not os.path.exists(beef_ultrasound_database_path): + download_ultrasound_images(current_dir) + # if there is no chosen scan number, let it be random + if scan_number is None: + rng = np.random.default_rng() + scan_number = rng.integers(low=2, high=63) + logger.debug("Scan number {} was used for this simulation".format(scan_number)) + ultrasound_image = np.load(beef_ultrasound_database_path + "/" + image_type + "/Scan_" + str(scan_number) + + ".npy") + return ultrasound_image diff --git a/simpa/utils/tags.py b/simpa/utils/tags.py index 1da32334..3add28b1 100644 --- a/simpa/utils/tags.py +++ b/simpa/utils/tags.py @@ -1585,3 +1585,15 @@ class Tags: """ Identifier for the volume fraction for the simulation """ + + MEAT_ULTRASOUND_CROPPED = "cropped_images" + """ + Flag to use the cropped beef ultrasound images + Usage: simpa.utils.libraries.structure_library.heterogeneity_generator + """ + + MEAT_ULTRASOUND_FULL = "full_images" + """ + Flag to use the full beef ultrasound images + Usage: simpa.utils.libraries.structure_library.heterogeneity_generator + """ diff --git a/simpa_examples/optical_and_acoustic_simulation_US_heterogeneity.py b/simpa_examples/optical_and_acoustic_simulation_US_heterogeneity.py new file mode 100644 index 00000000..d8baaf04 --- /dev/null +++ b/simpa_examples/optical_and_acoustic_simulation_US_heterogeneity.py @@ -0,0 +1,219 @@ +# SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ +# SPDX-FileCopyrightText: 2021 Janek Groehl +# SPDX-License-Identifier: MIT + +import simpa.utils.libraries.heterogeneity_generator +from simpa import Tags +import simpa as sp +import numpy as np +from simpa.utils.profiling import profile +from argparse import ArgumentParser + +# FIXME temporary workaround for newest Intel architectures +import os +os.environ["KMP_DUPLICATE_LIB_OK"] = "TRUE" + +# TODO: Please make sure that a valid path_config.env file is located in your home directory, or that you +# point to the correct file in the PathManager(). + + +@profile +def run_optical_and_acoustic_simulation(spacing: float | int = 0.2, path_manager=None, + visualise: bool = True): + """ + ########## + This example will (if not previously downloaded) download a folder with beef ultrasound images + ########## + + An example of the full phptoacoustic pipeline and reconstruction with a heterogeneous muscle blood volume fraction + and the MSOT AcuityEcho. + The Novelty of this example comes in the origin of its heterogeneous background. Here the heterogeneity come from an + ultrasound image taken of a piece of beef. For a full description of how the data was obtained, please refeer to the + md file in the downloaded folder. + + :param spacing: The simulation spacing between voxels + :param path_manager: the path manager to be used, typically sp.PathManager + :param visualise: If VISUALIZE is set to True, the reconstruction result will be plotted + :return: a run through of the example + """ + if path_manager is None: + path_manager = sp.PathManager() + VOLUME_TRANSDUCER_DIM_IN_MM = 60 + VOLUME_PLANAR_DIM_IN_MM = 40 + VOLUME_HEIGHT_IN_MM = 34 + RANDOM_SEED = 4711 + + # If VISUALIZE is set to True, the simulation result will be plotted + VISUALIZE = True + + def create_example_tissue(): + """ + This is a very simple example script of how to create a tissue definition. + It contains a muscular background, an epidermis layer on top of the muscles + and a blood vessel. + """ + dim_x, dim_y, dim_z = settings.get_volume_dimensions_voxels() + background_dictionary = sp.Settings() + background_dictionary[Tags.MOLECULE_COMPOSITION] = sp.TISSUE_LIBRARY.constant(1e-10, 1e-10, 1.0) + background_dictionary[Tags.STRUCTURE_TYPE] = Tags.BACKGROUND + + tissue_dict = sp.Settings() + tissue_dict[Tags.BACKGROUND] = background_dictionary + + us_heterogeneity = sp.ImageHeterogeneity(xdim=dim_x, ydim=dim_y, zdim=dim_z, + spacing_mm=spacing, target_min=0, target_max=0.05, + scaling_type=Tags.IMAGE_SCALING_SYMMETRIC) + us_heterogeneity.exponential(2) + us_heterogeneity.invert_image() + bvf = us_heterogeneity.get_map() + + muscle_dictionary = sp.Settings() + muscle_dictionary[Tags.PRIORITY] = 1 + muscle_dictionary[Tags.STRUCTURE_START_MM] = [0, 0, 0] + muscle_dictionary[Tags.STRUCTURE_END_MM] = [0, 0, 34] + muscle_dictionary[Tags.MOLECULE_COMPOSITION] = sp.TISSUE_LIBRARY.muscle(oxygenation=0.6, + blood_volume_fraction=bvf) + muscle_dictionary[Tags.CONSIDER_PARTIAL_VOLUME] = True + muscle_dictionary[Tags.ADHERE_TO_DEFORMATION] = True + muscle_dictionary[Tags.STRUCTURE_TYPE] = Tags.HORIZONTAL_LAYER_STRUCTURE + + tissue_dict["muscle"] = muscle_dictionary + return tissue_dict + + # Seed the numpy random configuration prior to creating the global_settings file in + # order to ensure that the same volume + # is generated with the same random seed every time. + + np.random.seed(RANDOM_SEED) + VOLUME_NAME = "CompletePipelineTestMSOT_"+str(RANDOM_SEED) + + general_settings = { + # These parameters set the general properties of the simulated volume + Tags.RANDOM_SEED: RANDOM_SEED, + Tags.VOLUME_NAME: "CompletePipelineExample_" + str(RANDOM_SEED), + Tags.SIMULATION_PATH: path_manager.get_hdf5_file_save_path(), + Tags.SPACING_MM: spacing, + Tags.DIM_VOLUME_Z_MM: VOLUME_HEIGHT_IN_MM, + Tags.DIM_VOLUME_X_MM: VOLUME_TRANSDUCER_DIM_IN_MM, + Tags.DIM_VOLUME_Y_MM: VOLUME_PLANAR_DIM_IN_MM, + Tags.VOLUME_CREATOR: Tags.VOLUME_CREATOR_VERSATILE, + Tags.GPU: True, + Tags.WAVELENGTHS: [700, 800], + Tags.DO_FILE_COMPRESSION: True, + Tags.DO_IPASC_EXPORT: True + } + settings = sp.Settings(general_settings) + np.random.seed(RANDOM_SEED) + + settings.set_volume_creation_settings({ + Tags.STRUCTURES: create_example_tissue(), + Tags.SIMULATE_DEFORMED_LAYERS: True + }) + + settings.set_optical_settings({ + Tags.OPTICAL_MODEL_NUMBER_PHOTONS: 1e7, + Tags.OPTICAL_MODEL_BINARY_PATH: path_manager.get_mcx_binary_path(), + Tags.ILLUMINATION_TYPE: Tags.ILLUMINATION_TYPE_MSOT_ACUITY_ECHO, + Tags.LASER_PULSE_ENERGY_IN_MILLIJOULE: 50, + Tags.MCX_ASSUMED_ANISOTROPY: 0.9, + Tags.ADDITIONAL_FLAGS: ['--printgpu'] # to print MCX GPU information + }) + + settings.set_acoustic_settings({ + Tags.ACOUSTIC_SIMULATION_3D: False, + Tags.ACOUSTIC_MODEL_BINARY_PATH: path_manager.get_matlab_binary_path(), + Tags.KWAVE_PROPERTY_ALPHA_POWER: 0.00, + Tags.KWAVE_PROPERTY_SENSOR_RECORD: "p", + Tags.KWAVE_PROPERTY_PMLInside: False, + Tags.KWAVE_PROPERTY_PMLSize: [31, 32], + Tags.KWAVE_PROPERTY_PMLAlpha: 1.5, + Tags.KWAVE_PROPERTY_PlotPML: False, + Tags.RECORDMOVIE: False, + Tags.MOVIENAME: "visualization_log", + Tags.ACOUSTIC_LOG_SCALE: True + }) + + settings.set_reconstruction_settings({ + Tags.RECONSTRUCTION_PERFORM_BANDPASS_FILTERING: False, + Tags.ACOUSTIC_MODEL_BINARY_PATH: path_manager.get_matlab_binary_path(), + Tags.ACOUSTIC_SIMULATION_3D: False, + Tags.KWAVE_PROPERTY_ALPHA_POWER: 0.00, + Tags.TUKEY_WINDOW_ALPHA: 0.5, + Tags.BANDPASS_CUTOFF_LOWPASS_IN_HZ: int(8e6), + Tags.BANDPASS_CUTOFF_HIGHPASS_IN_HZ: int(0.1e4), + Tags.RECONSTRUCTION_BMODE_AFTER_RECONSTRUCTION: False, + Tags.RECONSTRUCTION_BMODE_METHOD: Tags.RECONSTRUCTION_BMODE_METHOD_HILBERT_TRANSFORM, + Tags.RECONSTRUCTION_APODIZATION_METHOD: Tags.RECONSTRUCTION_APODIZATION_BOX, + Tags.RECONSTRUCTION_MODE: Tags.RECONSTRUCTION_MODE_PRESSURE, + Tags.KWAVE_PROPERTY_SENSOR_RECORD: "p", + Tags.KWAVE_PROPERTY_PMLInside: False, + Tags.KWAVE_PROPERTY_PMLSize: [31, 32], + Tags.KWAVE_PROPERTY_PMLAlpha: 1.5, + Tags.KWAVE_PROPERTY_PlotPML: False, + Tags.RECORDMOVIE: False, + Tags.MOVIENAME: "visualization_log", + Tags.ACOUSTIC_LOG_SCALE: True, + Tags.DATA_FIELD_SPEED_OF_SOUND: 1540, + Tags.DATA_FIELD_ALPHA_COEFF: 0.01, + Tags.DATA_FIELD_DENSITY: 1000, + Tags.SPACING_MM: spacing + }) + + settings["noise_initial_pressure"] = { + Tags.NOISE_MEAN: 1, + Tags.NOISE_STD: 0.01, + Tags.NOISE_MODE: Tags.NOISE_MODE_MULTIPLICATIVE, + Tags.DATA_FIELD: Tags.DATA_FIELD_INITIAL_PRESSURE, + Tags.NOISE_NON_NEGATIVITY_CONSTRAINT: True + } + + settings["noise_time_series"] = { + Tags.NOISE_STD: 1, + Tags.NOISE_MODE: Tags.NOISE_MODE_ADDITIVE, + Tags.DATA_FIELD: Tags.DATA_FIELD_TIME_SERIES_DATA + } + + # TODO: For the device choice, uncomment the undesired device + + device = sp.MSOTAcuityEcho(device_position_mm=np.array([VOLUME_TRANSDUCER_DIM_IN_MM/2, + VOLUME_PLANAR_DIM_IN_MM/2, + 0])) + device.update_settings_for_use_of_model_based_volume_creator(settings) + + SIMULATION_PIPELINE = [ + sp.ModelBasedAdapter(settings), + sp.MCXAdapter(settings), + sp.GaussianNoise(settings, "noise_initial_pressure"), + sp.KWaveAdapter(settings), + sp.GaussianNoise(settings, "noise_time_series"), + sp.DelayAndSumAdapter(settings), + sp.FieldOfViewCropping(settings) + ] + + sp.simulate(SIMULATION_PIPELINE, settings, device) + + if Tags.WAVELENGTH in settings: + WAVELENGTH = settings[Tags.WAVELENGTH] + else: + WAVELENGTH = 700 + + if visualise: + sp.visualise_data(path_to_hdf5_file=settings[Tags.SIMPA_OUTPUT_PATH], + wavelength=WAVELENGTH, + show_time_series_data=True, + show_initial_pressure=True, + show_reconstructed_data=True, + log_scale=False, + show_xz_only=False, + show_blood_volume_fraction=True) + + +if __name__ == "__main__": + parser = ArgumentParser(description='Run the optical and acoustic simulation example') + parser.add_argument("--spacing", default=0.2, type=float, help='the voxel spacing in mm') + parser.add_argument("--path_manager", default=None, help='the path manager, None uses sp.PathManager') + parser.add_argument("--visualise", default=True, type=bool, help='whether to visualise the result') + config = parser.parse_args() + + run_optical_and_acoustic_simulation(spacing=config.spacing, path_manager=config.path_manager, + visualise=config.visualise) diff --git a/simpa_examples/segmentation_loader.py b/simpa_examples/segmentation_loader.py index 4d6a383e..f80c3bba 100644 --- a/simpa_examples/segmentation_loader.py +++ b/simpa_examples/segmentation_loader.py @@ -8,6 +8,8 @@ from skimage.data import shepp_logan_phantom from scipy.ndimage import zoom from skimage.transform import resize +import nrrd +import random # FIXME temporary workaround for newest Intel architectures import os @@ -33,17 +35,65 @@ def run_segmentation_loader(spacing: float | int = 1.0, input_spacing: float | i if path_manager is None: path_manager = sp.PathManager() - label_mask = shepp_logan_phantom() + def us_heterogeneity(): + """ + + """ + scans = [2, 3, 14, 19, 32, 39, 50, 51] + scan = random.choice(scans) + + blood_volume_fraction = sp.ImageHeterogeneity(xdim=400, ydim=200, zdim=400, spacing_mm=spacing, target_min=0, + target_max=0.05, ultrasound_image_type=Tags.MEAT_ULTRASOUND_FULL, + scan_number=scan) + blood_volume_fraction.exponential(6) + blood_volume_fraction.invert_image() + + segmentation_mask, header = nrrd.read("./beef_ultrasound_database/segmentations/Scan_"+str(scan)+"_labels.nrrd") + segmentation_mask = np.repeat(np.swapaxes(segmentation_mask, 1, 2), 200, axis=1).astype(np.int32) + labels = header['org.mitk.multilabel.segmentation.labelgroups'] + + # Here, we extract the labels from the nrrd file created by the MITK Workbench used to segment our data. + lab_no = 1 + label_dict = {} + while True: + label, number, labels = labels.rpartition('"value":'+str(lab_no)+',') + actual_label = label.rpartition('"name":"', )[2].rpartition('","opacity":')[0] + if not actual_label: + break + label_dict[actual_label] = lab_no + lab_no += 1 + + # A fix for the fact only one has a fat layer on top. + try: + fat = label_dict['5 Fat'] + except KeyError: + fat = 7 + + background_dictionary = sp.Settings() + background_dictionary[Tags.MOLECULE_COMPOSITION] = sp.TISSUE_LIBRARY.constant(1e-10, 1e-10, 1.0) + background_dictionary[Tags.STRUCTURE_TYPE] = Tags.BACKGROUND + + tissue_dict = dict() + tissue_dict[Tags.BACKGROUND] = background_dictionary + tissue_dict[label_dict['1 Heavy Water']] = sp.TISSUE_LIBRARY.heavy_water() + tissue_dict[label_dict['2 Mediprene']] = sp.TISSUE_LIBRARY.mediprene() + tissue_dict[label_dict['3 US Gel']] = sp.TISSUE_LIBRARY.ultrasound_gel() + tissue_dict[label_dict['4 Muscle']] = sp.TISSUE_LIBRARY.muscle(oxygenation=0.4, + blood_volume_fraction=blood_volume_fraction.get_map()) + tissue_dict[fat] = sp.TISSUE_LIBRARY.subcutaneous_fat() + return tissue_dict, segmentation_mask + + def shepp_logan_phantom(): + label_mask = shepp_logan_phantom() + + label_mask = np.digitize(label_mask, bins=np.linspace(0.0, 1.0, 11), right=True) + label_mask = label_mask[100:300, 100:300] + label_mask = np.reshape(label_mask, (label_mask.shape[0], 1, label_mask.shape[1])) + + segmentation_volume_tiled = np.tile(label_mask, (1, 128, 1)) + segmentation_volume_mask = np.round(zoom(segmentation_volume_tiled, input_spacing / spacing, + order=0)).astype(int) - label_mask = np.digitize(label_mask, bins=np.linspace(0.0, 1.0, 11), right=True) - label_mask = label_mask[100:300, 100:300] - label_mask = np.reshape(label_mask, (label_mask.shape[0], 1, label_mask.shape[1])) - - segmentation_volume_tiled = np.tile(label_mask, (1, 128, 1)) - segmentation_volume_mask = np.round(zoom(segmentation_volume_tiled, input_spacing/spacing, - order=0)).astype(int) - - def segmentation_class_mapping(): ret_dict = dict() ret_dict[0] = sp.TISSUE_LIBRARY.heavy_water() ret_dict[1] = sp.TISSUE_LIBRARY.blood() @@ -60,7 +110,14 @@ def segmentation_class_mapping(): ret_dict[8] = sp.TISSUE_LIBRARY.heavy_water() ret_dict[9] = sp.TISSUE_LIBRARY.heavy_water() ret_dict[10] = sp.TISSUE_LIBRARY.heavy_water() - return ret_dict + + return ret_dict, segmentation_volume_mask + + ############################################### + # TODO: uncomment which example you wish to run + volume_settings = us_heterogeneity() + # volume_settings = shepp_logan_phantom() + ############################################### settings = sp.Settings() settings[Tags.SIMULATION_PATH] = path_manager.get_hdf5_file_save_path() @@ -68,13 +125,13 @@ def segmentation_class_mapping(): settings[Tags.RANDOM_SEED] = 1234 settings[Tags.WAVELENGTHS] = [700, 800] settings[Tags.SPACING_MM] = spacing - settings[Tags.DIM_VOLUME_X_MM] = segmentation_volume_mask.shape[0] * spacing - settings[Tags.DIM_VOLUME_Y_MM] = segmentation_volume_mask.shape[1] * spacing - settings[Tags.DIM_VOLUME_Z_MM] = segmentation_volume_mask.shape[2] * spacing + settings[Tags.DIM_VOLUME_X_MM] = volume_settings[1].shape[0] * spacing + settings[Tags.DIM_VOLUME_Y_MM] = volume_settings[1].shape[1] * spacing + settings[Tags.DIM_VOLUME_Z_MM] = volume_settings[1].shape[2] * spacing settings.set_volume_creation_settings({ - Tags.INPUT_SEGMENTATION_VOLUME: segmentation_volume_mask, - Tags.SEGMENTATION_CLASS_MAPPING: segmentation_class_mapping(), + Tags.INPUT_SEGMENTATION_VOLUME: volume_settings[1], + Tags.SEGMENTATION_CLASS_MAPPING: volume_settings[0], }) @@ -90,7 +147,16 @@ def segmentation_class_mapping(): sp.MCXAdapter(settings) ] - sp.simulate(pipeline, settings, sp.RSOMExplorerP50(element_spacing_mm=1.0)) + # TODO: For the device choice, uncomment the undesired device + device = sp.RSOMExplorerP50(element_spacing_mm=1.0) + # device = sp.MSOTAcuityEcho(device_position_mm=np.array([settings[Tags.DIM_VOLUME_X_MM] / 2, + # settings[Tags.DIM_VOLUME_Y_MM] / 2, + # 0])) + # device.update_settings_for_use_of_segmentation_based_volume_creator(settings, add_layers=[Tags.ADD_HEAVY_WATER], + # heavy_water_tag=2, + # current_heavy_water_depth=100 * spacing) + + sp.simulate(pipeline, settings, device) if Tags.WAVELENGTH in settings: WAVELENGTH = settings[Tags.WAVELENGTH]