SLIX
Scattered Light Imaging Toolbox (SLIX) – an open-source Python package that allows a fully automated evaluation of SLI measurements and the generation of different parameter maps
View Source
""" Scattered Light Imaging Toolbox (SLIX) – an open-source Python package that allows a fully automated evaluation of SLI measurements and the generation of different parameter maps """ __version__ = '2.4.0' __all__ = ['toolbox', 'io', 'visualization', 'preparation', 'attributemanager', 'classification'] from . import toolbox, io, visualization, preparation, attributemanager, classification
View Source
from SLIX.CPU import toolbox as cpu_toolbox from SLIX._logging import get_logger import numpy try: try: import cupy from numba import cuda cupy.empty(0, dtype=float) from SLIX.GPU import toolbox as gpu_toolbox gpu_available = True except cupy.cuda.runtime.CUDARuntimeError: print('[WARNING] CuPy is installed but an error was thrown by the ' 'runtime. SLIX will fall back to the CPU variant.') gpu_available = False except (cuda.cudadrv.driver.CudaAPIError, cuda.cudadrv.driver.LinkerError): get_logger("SLIX").info("Numba CUDA couldn't be initialized. " "Please check if there are problems with your CUDA / Numba " "version. SLIX will fall back to the CPU variant.") gpu_available = False except (ModuleNotFoundError, NameError): gpu_available = False get_logger("SLIX").info('CuPy is not installed. The toolbox will use the CPU ' 'variant instead. If you want to use the GPU variant, please run ' '`pip install cupy`.') __all__ = ['background_mask', 'centroid_correction', 'direction', 'unit_vectors', 'num_peaks', 'mean_peak_prominence', 'peaks', 'peak_prominence', 'peak_width', 'mean_peak_distance', 'peak_distance', 'mean_peak_width', 'significant_peaks'] def background_mask(image, use_gpu=gpu_available, return_numpy=True): """ Creates a background mask by setting all image pixels with low scattering signals to zero. As all background pixels are near zero for all images in the SLI image stack, this method should remove most of the background allowing for better approximations using the available features. It is advised to use this function. Args: image: Complete SLI measurement image stack as a 2D/3D NumPy array use_gpu: If available use the GPU for calculation return_numpy: Necessary if using `use_gpu`. Specifies if a CuPy or NumPy array will be returned. Returns: numpy.array: 1D/2D-image which masks the background as True and foreground as False """ if use_gpu: return gpu_toolbox.background_mask(image, return_numpy) else: return cpu_toolbox.background_mask(image) def peaks(image, use_gpu=gpu_available, return_numpy=True): """ Detect all peaks from a full SLI measurement. Peaks will not be filtered in any way. To detect only significant peaks, filter the peaks by using the prominence as a threshold. Args: image: Complete SLI measurement image stack as a 2D/3D NumPy array use_gpu: If available use the GPU for calculation return_numpy: Necessary if using `use_gpu`. Specifies if a CuPy or NumPy array will be returned. Returns: 2D/3D boolean image containing masking the peaks with `True` """ if use_gpu: return gpu_toolbox.peaks(image, return_numpy) else: return cpu_toolbox.peaks(image) def significant_peaks(image, low_prominence=cpu_toolbox.TARGET_PROMINENCE, high_prominence=numpy.inf, use_gpu=gpu_available, return_numpy=True): """ Detect all peaks from a full SLI measurement and filter them by passing thresholds. Args: image: Complete SLI measurement image stack as a 2D/3D NumPy array low_prominence: Minimum prominence needed by peak to count as a peak. Peaks below this threshold will not be considered as a peak. high_prominence: Maximum prominence needed by peak to count as a peak. Peaks below this threshold will not be considered as a peak. use_gpu: If available use the GPU for calculation return_numpy: Necessary if using `use_gpu`. Specifies if a CuPy or NumPy array will be returned. Returns: 2D/3D boolean image containing masking the peaks with `True` """ if use_gpu: peaks = gpu_toolbox.peaks(image, return_numpy=return_numpy) prominences = gpu_toolbox.peak_prominence(image, peaks, return_numpy=return_numpy) peaks[prominences < low_prominence] = False peaks[prominences > high_prominence] = False else: peaks = cpu_toolbox.peaks(image) prominences = cpu_toolbox.peak_prominence(image, peaks) peaks[prominences < low_prominence] = False peaks[prominences > high_prominence] = False return peaks def num_peaks(image, low_prominence=cpu_toolbox.TARGET_PROMINENCE, high_prominence=numpy.inf, use_gpu=gpu_available, return_numpy=True): """ Calculate the number of peaks from each line profile in an SLI image series by detecting all peaks and applying thresholds to remove unwanted peaks. Args: image: Full SLI measurement (series of images) which is prepared for the pipeline using the SLIX toolbox methods. low_prominence: Lower prominence bound for detecting a peak. high_prominence: Higher prominence bound for detecting a peak. use_gpu: If available use the GPU for calculation return_numpy: Necessary if using `use_gpu`. Specifies if a CuPy or NumPy array will be returned. Returns: Array where each entry corresponds to the number of detected peaks within the first dimension of the SLI image series. """ if use_gpu: peaks = significant_peaks(image, low_prominence, high_prominence, return_numpy=False) return gpu_toolbox.num_peaks(peak_image=peaks, return_numpy=return_numpy) else: peaks = significant_peaks(image, low_prominence, high_prominence, use_gpu=False) return cpu_toolbox.num_peaks(peak_image=peaks) def direction(peak_image, centroids, correction_angle=0, number_of_directions=3, use_gpu=gpu_available, return_numpy=True): """ Calculate up to `number_of_directions` direction angles based on the given peak positions. If more than `number_of_directions*2` peaks are present, no direction angle will be calculated to avoid errors. This will result in a direction angle of BACKGROUND_COLOR. The peak positions are determined by the position of the corresponding peak pairs (i.e. 6 peaks: 1+4, 2+5, 3+6). If two peaks are too far away or too near (outside of 180°±35°), the direction angle will be considered as invalid, resulting in a direction angle of BACKGROUND_COLOR. Args: correction_angle: Correct the resulting direction angle by the value. This is useful when the stack or camera was rotated. peak_image: Boolean NumPy array specifying the peak positions in the full SLI stack centroids: Centroids resulting from `centroid_correction` for more accurate results number_of_directions: Number of directions which shall be generated. use_gpu: If available use the GPU for calculation return_numpy: Necessary if using `use_gpu`. Specifies if a CuPy or NumPy array will be returned. Returns: NumPy array with the shape (x, y, `number_of_directions`) containing up to `number_of_directions` direction angles. x equals the number of pixels of the SLI image series. If a direction angle is invalid or missing, the array entry will be BACKGROUND_COLOR instead. """ if use_gpu: return gpu_toolbox.direction(peak_image, centroids, correction_angle, number_of_directions, return_numpy) else: return cpu_toolbox.direction(peak_image, centroids, correction_angle, number_of_directions) def peak_distance(peak_image, centroids, use_gpu=gpu_available, return_numpy=True): """ Calculate the mean peak distance in degrees between two corresponding peaks for each line profile in an SLI image series. Args: peak_image: Boolean NumPy array specifying the peak positions in the full SLI stack centroids: Use centroid calculation to better determine the peak position regardless of the number of measurements / illumination angles used. use_gpu: If available use the GPU for calculation return_numpy: Necessary if using `use_gpu`. Specifies if a CuPy or NumPy array will be returned. Returns: NumPy array of floating point values containing the peak distance of the line profiles in degrees in their respective peak position. The first peak of each peak pair will show the distance between peak_1 and peak_2 while the second peak will show 360 - (peak_2 - peak_1). """ if use_gpu: return gpu_toolbox.peak_distance(peak_image, centroids, return_numpy) else: return cpu_toolbox.peak_distance(peak_image, centroids) def mean_peak_distance(peak_image, centroids, use_gpu=gpu_available, return_numpy=True): """ Calculate the mean peak distance in degrees between two corresponding peaks for each line profile in an SLI image series. Args: peak_image: Boolean NumPy array specifying the peak positions in the full SLI stack centroids: Use centroid calculation to better determine the peak position regardless of the number of measurements / illumination angles used. use_gpu: If available use the GPU for calculation return_numpy: Necessary if using `use_gpu`. Specifies if a CuPy or NumPy array will be returned. Returns: NumPy array of floating point values containing the mean peak distance of the line profiles in degrees. """ if use_gpu: return gpu_toolbox.mean_peak_distance(peak_image, centroids, return_numpy) else: return cpu_toolbox.mean_peak_distance(peak_image, centroids) def peak_prominence(image, peak_image=None, kind_of_normalization=1, use_gpu=gpu_available, return_numpy=True): """ Calculate the peak prominence of all given peak positions within a line profile. The line profile will be normalized by dividing the line profile through its mean value. Therefore, values above 1 are possible. Args: image: Original line profile used to detect all peaks. This array will be further analyzed to better determine the peak positions. peak_image: Boolean NumPy array specifying the peak positions in the full SLI stack kind_of_normalization: Normalize given line profile by using a normalization technique based on the kind_of_normalization parameter. 0 : Scale line profile to be between 0 and 1 1 : Divide line profile through its mean value use_gpu: If available use the GPU for calculation return_numpy: Necessary if using `use_gpu`. Specifies if a CuPy or NumPy array will be returned. Returns: Floating point value containing the mean peak prominence of the line profile in degrees. """ if use_gpu: return gpu_toolbox.peak_prominence(image, peak_image, kind_of_normalization, return_numpy) else: return cpu_toolbox.peak_prominence(image, peak_image, kind_of_normalization) def mean_peak_prominence(image, peak_image=None, kind_of_normalization=1, use_gpu=gpu_available, return_numpy=True): """ Calculate the mean peak prominence of all given peak positions within a line profile. The line profile will be normalized by dividing the line profile through its mean value. Therefore, values above 1 are possible. Args: image: Original line profile used to detect all peaks. This array will be further analyzed to better determine the peak positions. peak_image: Boolean NumPy array specifying the peak positions in the full SLI stack kind_of_normalization: Normalize given line profile by using a normalization technique based on the kind_of_normalization parameter. 0 : Scale line profile to be between 0 and 1 1 : Divide line profile through its mean value use_gpu: If available use the GPU for calculation return_numpy: Necessary if using `use_gpu`. Specifies if a CuPy or NumPy array will be returned. Returns: Floating point value containing the mean peak prominence of the line profile in degrees. """ if use_gpu: return gpu_toolbox.mean_peak_prominence(image, peak_image, kind_of_normalization, return_numpy) else: return cpu_toolbox.mean_peak_prominence(image, peak_image, kind_of_normalization) def peak_width(image, peak_image=None, target_height=0.5, use_gpu=gpu_available, return_numpy=True): """ Calculate the peak width of all given peak positions within a line profile. Args: image: Original line profile used to detect all peaks. This array will be further analyzed to better determine the peak positions. peak_image: Boolean NumPy array specifying the peak positions in the full SLI stack target_height: Relative peak height in relation to the prominence of the given peak. use_gpu: If available use the GPU for calculation return_numpy: Necessary if using `use_gpu`. Specifies if a CuPy or NumPy array will be returned. Returns: NumPy array where each entry corresponds to the peak width of the line profile. The values are in degree. """ if use_gpu: return gpu_toolbox.peak_width(image, peak_image, target_height, return_numpy=return_numpy) else: return cpu_toolbox.peak_width(image, peak_image, target_height) def mean_peak_width(image, peak_image=None, target_height=0.5, use_gpu=gpu_available, return_numpy=True): """ Calculate the mean peak width of all given peak positions within a line profile. Args: image: Original line profile used to detect all peaks. This array will be further analyzed to better determine the peak positions. peak_image: Boolean NumPy array specifying the peak positions in the full SLI stack target_height: Relative peak height in relation to the prominence of the given peak. use_gpu: If available use the GPU for calculation return_numpy: Necessary if using `use_gpu`. Specifies if a CuPy or NumPy array will be returned. Returns: NumPy array where each entry corresponds to the mean peak width of the line profile. The values are in degree. """ if use_gpu: return gpu_toolbox.mean_peak_width(image, peak_image, target_height, return_numpy=return_numpy) else: return cpu_toolbox.mean_peak_width(image, peak_image, target_height) def centroid_correction(image, peak_image, low_prominence=cpu_toolbox.TARGET_PROMINENCE, high_prominence=numpy.inf, use_gpu=gpu_available, return_numpy=True): """ Correct peak positions from a line profile by looking at only the peak with a given threshold using a centroid calculation. If a minimum is found in the considered interval, this minimum will be used as the limit instead. The range for the peak correction is limited by MAX_DISTANCE_FOR_CENTROID_ESTIMATION. Args: image: Original line profile used to detect all peaks. This array will be further analyzed to better determine the peak positions. peak_image: Boolean NumPy array specifying the peak positions in the full SLI stack low_prominence: Lower prominence bound for detecting a peak. high_prominence: Higher prominence bound for detecting a peak. use_gpu: If available use the GPU for calculation return_numpy: Necessary if using `use_gpu`. Specifies if a CuPy or NumPy array will be returned. Returns: _numpy array with the positions of all detected peak positions corrected with the centroid calculation. """ if use_gpu: return gpu_toolbox.centroid_correction(image, peak_image, low_prominence, high_prominence, return_numpy) else: return cpu_toolbox.centroid_correction(image, peak_image, low_prominence, high_prominence) def unit_vectors(direction, use_gpu=gpu_available, return_numpy=True): """ Calculate the unit vectors (UnitX, UnitY) from a given direction angle. Args: direction: 3D NumPy array - direction angles in degrees use_gpu: If available use the GPU for calculation return_numpy: Necessary if using `use_gpu`. Specifies if a CuPy or NumPy array will be returned. Returns: UnitX, UnitY: 3D NumPy array, 3D NumPy array x- and y-vector component in arrays """ if use_gpu: return gpu_toolbox.unit_vectors(direction, return_numpy=return_numpy) else: return cpu_toolbox.unit_vectors(direction) def unit_vectors_3d(direction, inclination, use_gpu=gpu_available, return_numpy=True): """ Calculate the unit vectors (UnitX, UnitY) from a given direction angle. Args: direction: 3D NumPy array - direction angles in degrees inclination: 3D NumPy array - inclination angles in degrees use_gpu: If available use the GPU for calculation return_numpy: Necessary if using `use_gpu`. Specifies if a CuPy or NumPy array will be returned. Returns: UnitX, UnitY: 3D NumPy array, 3D NumPy array x- and y-vector component in arrays """ if use_gpu: return gpu_toolbox.unit_vectors_3d(direction, inclination, return_numpy=return_numpy) else: return cpu_toolbox.unit_vectors_3d(direction, inclination) def inclination_sign(peak_image, centroids, correction_angle=0, use_gpu=gpu_available, return_numpy=True): """ Calculate the inclination sign from the peak positions. The inclination sign is based on the peak distance between two peaks. Explanation of the results: -1: The minimal peak distance is behind the first peak (wrapping around) 0: This pixel / line profile has more than two peaks 1: The minimal peak distance is in front of the first peak. Args: peaks: 3D NumPy array - peak positions use_gpu: If available use the GPU for calculation return_numpy: Necessary if using `use_gpu`. Specifies if a CuPy or NumPy array will be returned. Returns: inclination_sign: 3D NumPy array inclination sign Raises: ValueError: If neither peaks nor peak_distance are given. ValueError: If peaks or peak_distance are given, but are a 2D array. """ if use_gpu: return gpu_toolbox.inclination_sign(peak_image, centroids, correction_angle, return_numpy=return_numpy) else: return cpu_toolbox.inclination_sign(peak_image, centroids, correction_angle)
View Source
from PIL import Image import numpy import tifffile import nibabel import h5py import sys import re import os import glob import datetime import SLIX import nibabel.openers from .attributemanager import AttributeHandler from ._logging import get_logger __all__ = ['H5FileReader', 'H5FileWriter', 'imread', 'imwrite', 'imwrite_rgb'] nibabel.openers.Opener.default_compresslevel = 9 _fileregex = r'.*_+p[0-9]+_?.*\.(tif{1,2}|jpe*g|nii|h5|png)' class H5FileReader: """ This class allows to read HDF5 files from your file system. It supports reading datasets but not reading attributes. """ def __init__(self): self.path = None self.file = None self.content = None def open(self, path): """ Args: path: Path on the filesystem to the HDF5 file which will be read Returns: None """ if not path == self.path: self.close() self.path = path self.file = h5py.File(path, 'r') def close(self): """ Close the currently opened file, if any is open. Returns: None """ if self.file is not None: self.file.close() self.file = None self.path = None self.content = None def read(self, dataset): """ Read a dataset from the currently opened file, if any is open. The content of the dataset will be stored for future use. Args: dataset: Path to the dataset within the HDF5 Returns: The opened dataset. """ if self.content is None: self.content = {} if dataset not in self.content.keys(): self.content[dataset] = numpy.squeeze(self.file[dataset][:]) return self.content[dataset] class H5FileWriter: """ This class allows to write HDF5 files to your file system. It supports writing datasets and writing attributes. """ def __init__(self): self.path = None self.file = None def add_symlink(self, dataset, symlink_path): """ Adds a symbolic link from one dataset to another path. Args: dataset: Dataset path within the HDF5 symlink_path: Path to the symlink. Returns: None """ if self.file is None: return self.file[symlink_path] = self.file[dataset] self.file.flush() def add_plim_attributes(self, stack_path, dataset='/Image'): """ PLIM is a package used in the 3D-PLI group to read and write multiple attributes from/to a HDF5 file. The basic functionality is added in attributemanager.py. Calling this method will write many attributes to the HDF5 file at the given dataset. This includes: A unique ID, the creator, software parameters, creation time, software_revision, image_modality and all attributes from the original stack, if it was a HDF5 Args: stack_path: Path of the stack that was used to calculate the content which will be written to the HDF5 file. dataset: Dataset where the attributes shall be written to. Returns: None """ if self.path is None or self.file is None: return if dataset not in self.file: self.file.create_dataset(dataset) output_handler = AttributeHandler(self.file[dataset]) if stack_path[:-3] == ".h5": original_file = h5py.File(stack_path, 'r') original_dataset = original_file[dataset] original_handler = AttributeHandler(original_dataset) original_handler.copy_all_attributes_to(output_handler) original_file.close() output_handler.add_creator() output_handler.set_attribute('software', sys.argv[0]) output_handler.set_attribute('software_revision', SLIX.__version__) output_handler.set_attribute('creation_time', datetime.datetime.now() .strftime('%Y-%m-%d %H:%M:%S')) output_handler.set_attribute('software_parameters', ' '.join(sys.argv[1:])) output_handler.set_attribute('image_modality', "Placeholder") output_handler.add_id() self.file.flush() def write_attribute(self, dataset, attr_name, value): """ Write a single attribute to a dataset. Args: dataset: Path to the dataset within the HDF5 attr_name: Name of the attribute which shall be written. value: Value of the attribute that shall be written. Returns: None """ if self.file is None: return if dataset not in self.file: self.file.create_dataset(dataset) output_handler = AttributeHandler(self.file[dataset]) output_handler.set_attribute(attr_name, value) self.file.flush() def write_dataset(self, dataset, content): """ Write a dataset to the currently opened HDF5 file, if any is open. If no HDF5 file is open, this method will return without writing. Args: dataset: Path to the dataset within the HDF5 file. content: Content which shall be written. Returns: None """ if self.file is None: return if dataset not in self.file: # Change compression algorithm for large files as it can take # very long for the compression to finish if len(content.shape) == 3: self.file.create_dataset(dataset, content.shape, dtype=content.dtype, data=content, compression='gzip', compression_opts=5, shuffle=True) else: self.file.create_dataset(dataset, content.shape, dtype=content.dtype, data=content, compression='gzip', compression_opts=9, shuffle=True) else: self.file[dataset] = content self.file.flush() def close(self): """ Close the currently opened file. Returns: None """ if self.file is None: return self.file.flush() self.file.close() self.path = None self.file = None def open(self, path): """ Open a new HDF5 file with the given path. If another file was opened, it will be closed first. Args: path: Path to the HDF5 file. Returns: None """ if self.path != path: self.close() self.path = path self.file = h5py.File(path, mode='w') def read_folder(filepath): """ Reads multiple image files from a folder and returns the resulting stack. To find the images in the right order, a regex is used which will search for files with the following pattern: [prefix]_p[Nr][suffix]. The start number doesn't need to be 0. The files are sorted with a natural sort, meaning that files like 0002, 1, 004, 3 will be sorted as 1, 0002, 3, 004. The follwing regex is used to find the measurements: ".*_+p[0-9]+_?.*.(tif{1,2}|jpe*g|nii|h5|png)" Supported file formats for the image file equal the supported formats of SLIX.imread. Args: filepath: Path to folder Returns: numpy.array: Image with shape [x, y, z] where [x, y] is the size of a single image and z specifies the number of measurements """ files_in_folder = glob.glob(filepath + '/*') matching_files = [] for file in files_in_folder: if re.match(_fileregex, file) is not None: matching_files.append(file) matching_files.sort(key=__natural_sort_filenames_key) image = None # Check if files contain the needed regex for our measurements for file in matching_files: measurement_image = imread(file) if image is None: image = measurement_image elif len(image.shape) == 2: image = numpy.stack((image, measurement_image), axis=-1) else: image = numpy.concatenate((image, measurement_image [:, :, numpy.newaxis]), axis=-1) return image def __natural_sort_filenames_key(string, regex=re.compile('([0-9]+)')): return [int(text) if text.isdigit() else text.lower() for text in regex.split(string)] def imread(filepath, dataset="/Image"): """ Reads image file and returns it. Supported file formats: HDF5, NIfTI, Tiff. Args: filepath: Path to image dataset: When reading a HDF5 file, a dataset is required. Default: '/Image' Returns: numpy.array: Image with shape [x, y, z] where [x, y] is the size of a single image and z specifies the number of measurements """ if os.path.isdir(filepath): data = read_folder(filepath) # Load NIfTI dataset elif filepath.endswith('.nii') or filepath.endswith('.nii.gz'): data = nibabel.load(filepath).get_fdata() data = numpy.squeeze(numpy.swapaxes(data, 0, 1)) elif filepath.endswith('.tiff') or filepath.endswith('.tif'): data = tifffile.imread(filepath) if len(data.shape) == 3: data = numpy.squeeze(numpy.moveaxis(data, 0, -1)) elif filepath.endswith('.h5'): reader = H5FileReader() reader.open(filepath) data = reader.read(dataset) if len(data.shape) == 3: data = numpy.squeeze(numpy.moveaxis(data, 0, -1)) reader.close() return data else: data = numpy.array(Image.open(filepath)) return data def imwrite(filepath, data, dataset='/Image', original_stack_path=""): """ Write generated image to given filepath. Supported file formats: HDF5, NIfTI, Tiff. Other file formats are only indirectly supported and might result in errors. Args: filepath: Path to image data: Data which will be written to the disk dataset: When writing a HDF5 file, a dataset is required. Default: '/Image' original_stack_path: Path to the original image stack used to create this content. Only required when a HDF5 file is written. Returns: None """ save_data = data.copy() if save_data.dtype == bool: save_data = save_data.astype(numpy.uint8) elif save_data.dtype == numpy.float64: save_data = save_data.astype(numpy.float32) elif save_data.dtype == numpy.int64: save_data = save_data.astype(numpy.int32) elif save_data.dtype == numpy.uint64: save_data = save_data.astype(numpy.uint32) if filepath.endswith('.nii') or filepath.endswith('.nii.gz'): save_data = numpy.swapaxes(save_data, 0, 1) nibabel.save(nibabel.Nifti1Image(save_data, numpy.eye(4)), filepath) elif filepath.endswith('.tiff') or filepath.endswith('.tif'): if len(save_data.shape) == 3: save_data = numpy.moveaxis(save_data, -1, 0) tifffile_version_date = datetime.datetime.strptime( tifffile.__version__, '%Y.%m.%d') tifffile_comparison_date = datetime.datetime.strptime( '2020.10.02', '%Y.%m.%d') if tifffile_version_date > tifffile_comparison_date: tifffile.imwrite(filepath, save_data, compression=8) else: tifffile.imwrite(filepath, save_data, compress=9) elif filepath.endswith('.h5'): if len(save_data.shape) == 3: save_data = numpy.moveaxis(save_data, -1, 0) writer = H5FileWriter() writer.open(filepath) writer.write_dataset(dataset, save_data) writer.add_plim_attributes(original_stack_path, dataset) writer.add_symlink(dataset, '/pyramid/00') writer.close() else: Image.fromarray(save_data).save(filepath) def imwrite_rgb(filepath, data, dataset='/Image', original_stack_path=""): """ Write generated RGB image to given filepath. Supported file formats: HDF5, Tiff. Other file formats are only indirectly supported and might result in errors. Args: filepath: Path to image data: Data which will be written to the disk dataset: When reading a HDF5 file, a dataset is required. Default: '/Image' original_stack_path: Path to the original image stack used to create this content. Only required when a HDF5 file is written. Returns: None """ logger = get_logger("imwrite_rgb") save_data = data.copy() axis = numpy.argwhere(numpy.array(save_data.shape) == 3).flatten() if len(axis) == 0: logger.error('Cannot create RGB image as no dimension has a depth of 3.') return if filepath.endswith('.tiff') or filepath.endswith('.tif'): save_data = numpy.moveaxis(save_data, axis[0], 0) tifffile.imwrite(filepath, save_data, photometric='rgb', compression=8) elif filepath.endswith('.h5'): writer = H5FileWriter() writer.open(filepath) writer.write_dataset(dataset, save_data) writer.add_plim_attributes(original_stack_path, dataset) writer.add_symlink(dataset, '/pyramid/00') writer.close() else: logger.error("File type is not supported. " "Supported file types are .h5, .tif(f)") def check_output_dir(folder: str) -> bool: """ Check if the output directory exists. If not, create it. Also checks for read/write permissions. Args: folder: Path to the output directory Returns: True if the directory exists and is writable, False otherwise """ logger = get_logger("check_output_dir") if not os.path.exists(folder): os.makedirs(folder, exist_ok=True) # Check if the output path is writable if not os.access(folder, os.W_OK): logger.error('Output path is not writable. Please choose a valid output ' 'path!') return False return True
View Source
import copy import numpy import tqdm from PIL import Image from matplotlib import pyplot as plt from matplotlib.colors import hsv_to_rgb, rgb_to_hsv from SLIX._visualization import _downsample, _plot_axes_unit_vectors, \ _visualize_multiple_direction, \ _visualize_one_direction __all__ = ['parameter_map', 'unit_vectors', 'unit_vector_distribution', 'direction', 'Colormap'] class Colormap: @staticmethod def prepare(direction: numpy.ndarray, inclination: numpy.ndarray) -> (numpy.ndarray, numpy.ndarray): if direction.max(axis=None) > numpy.pi and not numpy.isclose(direction.max(axis=None), numpy.pi): direction = numpy.deg2rad(direction) if inclination.max(axis=None) > numpy.pi and not numpy.isclose(inclination.max(axis=None), numpy.pi): inclination = numpy.deg2rad(inclination) # If inclination is only 2D and direction is 3D, we need to make sure that the # inclination matches the shape of the direction. if inclination.ndim == 2 and direction.ndim == 3: inclination = inclination[..., numpy.newaxis] if inclination.ndim == 3 and inclination.shape[-1] != direction.shape[-1]: inclination = numpy.repeat(inclination, direction.shape[-1], axis=-1) return direction, inclination @staticmethod def hsv_white(direction: numpy.ndarray, inclination: numpy.ndarray) -> numpy.ndarray: direction, inclination = Colormap.prepare(direction, inclination) hsv_stack = numpy.stack((direction / numpy.pi, 1.0 - (2 * inclination / numpy.pi), numpy.ones(direction.shape))) hsv_stack = numpy.moveaxis(hsv_stack, 0, -1) return numpy.clip(hsv_to_rgb(hsv_stack), 0, 1) @staticmethod def hsv_black(direction: numpy.ndarray, inclination: numpy.ndarray) -> numpy.ndarray: direction, inclination = Colormap.prepare(direction, inclination) hsv_stack = numpy.stack((direction / numpy.pi, numpy.ones(direction.shape), 1.0 - (2 * inclination / numpy.pi))) hsv_stack = numpy.moveaxis(hsv_stack, 0, -1) return numpy.clip(hsv_to_rgb(hsv_stack), 0, 1) @staticmethod def rgb(direction: numpy.ndarray, inclination: numpy.ndarray) -> numpy.ndarray: direction, inclination = Colormap.prepare(direction, inclination) direction[direction > numpy.pi / 2] = numpy.pi - direction[direction > numpy.pi / 2] rgb_stack = numpy.stack(( numpy.cos(inclination) * numpy.cos(direction), numpy.cos(inclination) * numpy.sin(direction), numpy.sin(inclination) )) rgb_stack = numpy.moveaxis(rgb_stack, 0, -1) return numpy.clip(rgb_stack, 0, 1) @staticmethod def hsv_black_reverse(direction: numpy.ndarray, inclination: numpy.ndarray) -> numpy.ndarray: direction, inclination = Colormap.prepare(direction, inclination) direction = numpy.clip(numpy.abs(-numpy.pi + direction), 0, numpy.pi) return Colormap.hsv_black(direction, inclination) @staticmethod def hsv_white_reverse(direction: numpy.ndarray, inclination: numpy.ndarray) -> numpy.ndarray: direction, inclination = Colormap.prepare(direction, inclination) direction = numpy.clip(numpy.abs(-numpy.pi + direction), 0, numpy.pi) return Colormap.hsv_white(direction, inclination) @staticmethod def rgb_reverse(direction: numpy.ndarray, inclination: numpy.ndarray) -> numpy.ndarray: direction, inclination = Colormap.prepare(direction, inclination) direction = numpy.clip(numpy.abs(-numpy.pi + direction), 0, numpy.pi) return Colormap.rgb(direction, inclination) def parameter_map(parameter_map, fig=None, ax=None, alpha=1, cmap='viridis', vmin=0, vmax=None, colorbar=True): """ This method will create a Matplotlib plot based on imshow to display the given parameter map in different colors. The parameter map is plotted to the current axis and figure. If neither is given, the method will create a new subfigure. To show the results, please use pyplot.show(). Args: parameter_map: 2D parameter map calculated with SLIX.toolbox. fig: Matplotlib figure. If None, a new subfigure will be created for fig and ax. ax: Matplotlib axis. If None, a new subfigure will be created for fig and ax. alpha: Apply alpha to Matplotlib plots to overlay them with some other image like the averaged transmitted light intensity. cmap: Matplotlib color map which is used for displaying the image. vmin: Minimum value in the resulting plot. If any value is below vmin, it will be displayed in black. vmax: Maximum value in the resulting plot. If any value is above vmax, it will be displayed in white. colorbar: Boolean value controlling if a color bar will be displayed in the current subplot. Returns: The current Matplotlib figure and axis. The image can be shown with pyplot.show(). """ if fig is None or ax is None: fig, ax = plt.subplots(1, 1) cmap_mod = copy.copy(plt.get_cmap(cmap)) im = ax.imshow(parameter_map, interpolation='nearest', cmap=cmap_mod, alpha=alpha) im.cmap.set_under(color='k') # Color for values less than vmin im.cmap.set_over(color='w') # Color for values more than vmax im.set_clim(vmin, vmax) ax.axis('off') if colorbar: fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04) return fig, ax def color_bubble(colormap: Colormap, shape=(1000, 1000, 3)) -> numpy.ndarray: """ Based on the chosen colormap in methods like unit_vectors or direction, the user might want to see the actual color bubble to understand the shown orientations. This method creates an empty numpy array and fills it with values based on the circular orientation from the middle point. The color can be directed from the colormap argument Args: colormap: Colormap which will be used to create the color bubble shape: Shape of the resulting color bubble. Returns: NumPy array containing the color bubble """ # create a meshgrid of the shape with the position of each pixel x, y = numpy.meshgrid(numpy.arange(shape[0]), numpy.arange(shape[1])) # center of our color_bubble center = numpy.array([shape[0]/2, shape[1]/2]) # radius where a full circle is still visible radius = numpy.minimum(numpy.minimum(center[0], center[1]), numpy.minimum(shape[0] - center[0], shape[1] - center[1])) # calculate the direction as the angle between the center and the pixel direction = numpy.pi - numpy.arctan2(y - center[0], x - center[1]) % numpy.pi # calculate the inclination as the distance between the center and the pixel inclination = numpy.sqrt((y - center[0])**2 + (x - center[1])**2) # normalize the inclination to a range of 0 to 90 degrees where 0 degree is at a distance of radius # and 90 degree is at a distance of 0 inclination = 90 - inclination / radius * 90 # create the color bubble color_bubble = colormap(direction, inclination) color_bubble[inclination < 0] = 0 return (255.0 * color_bubble).astype('uint8') def unit_vectors(unit_x, unit_y, ax=None, thinout=20, scale=-1, vector_width=1, alpha=0.8, background_threshold=0.5, background_value=0, colormap=Colormap.hsv_black, weighting=None): """ This method will create a Matplotlib plot based on quiver to represent the given unit vectors as colored lines (vector map). Parameters like thinout can be used to reduce the computing load. If thinout = 1, the resulting vectors might not be visible without zooming in significantly. Here, the vectors will only be plotted to the current axis. To show the results, please use pyplot.show(). Args: unit_x: Unit vector components along the x-axis (3D NumPy array). unit_y: Unit vector components along the y-axis (3D NumPy array). thinout: Downscaling parameter N (defines how many vectors N x N are replaced by one vector). Unit vectors will be thinned out using downscaling and thinning in combination. This will increase the vector size in the resulting image but will also reduce the information density. Please use with caution. scale: Increase the vector length by the given scale. Vectors will be longer and might overlap if the scale is too high. ax: Matplotlib axis. If None, the current context axis will be used. vector_width: When choosing a high scale, the vectors might appear quite thin which results in hard to read images. This option allows to increase the vector thickness to improve visibility. alpha: Apply alpha to Matplotlib plots to overlay them with some other other image like the averaged transmitted light intensity. background_threshold: If the fraction of background pixels (number of pixels without vector within N x N pixels) is below this threshold, the downscaled pixel will not show a vector. background_value: Background value of the parameter map. This is generally 0 in both axes for unit vector maps but can differ if another threshold was set. background_value: Fraction of background pixels in the considered (N x N) area for which the image pixels are set to background_value. If the fraction of background pixels lies above this defined threshold, background pixels will not be considered for computing the median. colormap: The colormap to use. Default is HSV black. The available color maps can be found in the colormap class. weighting: Weighting of the vectors. If None, the vectors will be weighted by a value of one, resulting in normal unit vectors. Returns: The current Matplotlib axis. The image can be shown with pyplot.show(). """ if ax is None: ax = plt.gca() while len(unit_x.shape) < 3: unit_x = unit_x[..., numpy.newaxis] while len(unit_y.shape) < 3: unit_y = unit_y[..., numpy.newaxis] # The default scale is below zero to allow the user to define his own scale # A scale below zero isn't valid for visualization. If the user # defines no scale, we suspect that the user wants an image # where each vector has a scale of one. Therefor we set the scale to # the same as our thinout when we draw the image. if scale < 0: scale = thinout if thinout > 1: downscaled_unit_x = _downsample(unit_x, thinout, background_threshold, background_value) downscaled_unit_y = _downsample(unit_y, thinout, background_threshold, background_value) while len(downscaled_unit_x.shape) < 3: downscaled_unit_x = downscaled_unit_x[..., numpy.newaxis] while len(downscaled_unit_y.shape) < 3: downscaled_unit_y = downscaled_unit_y[..., numpy.newaxis] # Rescale images to original dimensions for i in range(unit_x.shape[2]): unit_x[:, :, i] = numpy.array( Image.fromarray(downscaled_unit_x[:, :, i]) .resize(unit_x.shape[:2][::-1], Image.NEAREST) ) unit_y[:, :, i] = numpy.array( Image.fromarray(downscaled_unit_y[:, :, i]) .resize(unit_y.shape[:2][::-1], Image.NEAREST) ) del downscaled_unit_y del downscaled_unit_x if weighting is not None: while len(weighting.shape) < 3: weighting = weighting[..., numpy.newaxis] downscaled_weighting = _downsample(weighting, thinout, 0, 0) weighting = numpy.array( Image.fromarray(downscaled_weighting) .resize(weighting.shape[:2][::-1], Image.NEAREST) ) weighting = weighting[thinout // 2::thinout, thinout // 2::thinout] weighting = weighting.flatten() for i in range(unit_x.shape[2]): mesh_x, mesh_y = numpy.meshgrid(numpy.arange(thinout // 2, unit_x.shape[1], thinout), numpy.arange(thinout // 2, unit_x.shape[0], thinout)) mesh_u = unit_x[thinout // 2::thinout, thinout // 2::thinout, i] mesh_v = unit_y[thinout // 2::thinout, thinout // 2::thinout, i] _plot_axes_unit_vectors(ax, mesh_x.flatten(), mesh_y.flatten(), mesh_u.flatten(), mesh_v.flatten(), scale, alpha, vector_width, weighting, colormap) return ax def unit_vector_distribution(unit_x, unit_y, ax=None, thinout=20, scale=-1, vector_width=1, alpha=0.01, colormap=Colormap.hsv_black, weighting=None): """ This method will create a Matplotlib plot based on quiver to represent the given unit vectors as colored lines (vector map). Instead of showing a single vector like in unit_vector, here each vector will be shown in the resulting image. The thinout parameter will determine how many vectors will be overlapping. It is recommended to use a very small alpha value to see which directions in the resulting plot are dominant. Here, the vectors will only be plotted to the current axis. To show the results, please use pyplot.show(). The result might need some time to show depending on the input image size. Args: unit_x: Unit vector components along the x-axis (3D NumPy array). unit_y: Unit vector components along the y-axis (3D NumPy array). thinout: Downscaling parameter N (defines how many vectors N x N are replaced by one vector). Unit vectors will be thinned out using downscaling and thinning in combination. This will increase the vector size in the resulting image but will also reduce the information density. Please use with caution. scale: Increase the vector length by the given scale. Vectors will be longer and might overlap if the scale is too high. ax: Matplotlib axis. If None, the current context axis will be used. vector_width: When choosing a high scale, the vectors might appear quite thin which results in hard to read images. This option allows to increase the vector thickness to improve visibility. alpha: Apply alpha to Matplotlib plots to overlay them with some other other image like the averaged transmitted light intensity. colormap: The colormap to use. Default is HSV black. The available color maps can be found in the colormap class. weighting: Weighting of the vectors. If None, the vectors will be weighted by a value of one, resulting in normal unit vectors. Returns: The current Matplotlib axis. The image can be shown with pyplot.show(). """ if ax is None: ax = plt.gca() while len(unit_x.shape) < 3: unit_x = unit_x[..., numpy.newaxis] while len(unit_y.shape) < 3: unit_y = unit_y[..., numpy.newaxis] # The default scale is below zero to allow the user to define his own scale # A scale below zero isn't valid for visualization. If the user # defines no scale, we suspect that the user wants an image # where each vector has a scale of one. Therefore we set the scale to # the same as our thinout when we draw the image. if scale < 0: scale = thinout mesh_x = numpy.empty(unit_x.size) mesh_y = numpy.empty(unit_x.size) mesh_u = numpy.empty(unit_x.size) mesh_v = numpy.empty(unit_x.size) mesh_weighting = numpy.empty(unit_x.size) idx = 0 progress_bar = tqdm.tqdm(total=thinout * thinout, desc='Creating unit vectors.') for offset_x in range(thinout): for offset_y in range(thinout): progress_bar.update(1) for i in range(unit_x.shape[2]): mesh_x_it, mesh_y_it = numpy.meshgrid( numpy.arange(0, unit_x.shape[1] - offset_x, thinout), numpy.arange(0, unit_x.shape[0] - offset_y, thinout) ) mesh_x_it = mesh_x_it.flatten() mesh_y_it = mesh_y_it.flatten() mesh_u_it = unit_x[offset_y::thinout, offset_x::thinout, i] \ .flatten() mesh_v_it = unit_y[offset_y::thinout, offset_x::thinout, i] \ .flatten() if weighting is not None: mesh_weighting_it = weighting[offset_y::thinout, offset_x::thinout] \ .flatten() else: mesh_weighting_it = numpy.ones(mesh_u_it.size) mesh_x[idx:idx + len(mesh_x_it)] = mesh_x_it mesh_y[idx:idx + len(mesh_y_it)] = mesh_y_it mesh_u[idx:idx + len(mesh_u_it)] = mesh_u_it mesh_v[idx:idx + len(mesh_v_it)] = mesh_v_it mesh_weighting[idx:idx + len(mesh_weighting_it)] = mesh_weighting_it idx = idx + len(mesh_x_it) progress_bar.set_description('Finished. Plotting unit vectors.') _plot_axes_unit_vectors(ax, mesh_x, mesh_y, mesh_u, mesh_v, scale, alpha, vector_width, mesh_weighting, colormap) progress_bar.set_description('Done') progress_bar.close() return ax def direction(direction, inclination=None, saturation=None, value=None, colormap=Colormap.hsv_black): """ Generate a 2D colorized direction image in the HSV color space based on the original direction. Value and saturation of the color will always be one. The hue is determined by the direction. If the direction parameter is only a 2D numpy array, the result will be a simple orientation map where each pixel contains the HSV value corresponding to the direction angle. When a 3D stack with max. three directions is used, the result will be different. The resulting image will have two times the width and height. Each 2x2 square will show the direction angle of up to three directions. Depending on the number of directions, the following pattern is used to show the different direction angles. 1 direction: 1 1 1 1 2 directions: 1 2 2 1 3 directions: 1 2 3 0 Args: direction: 2D or 3D Numpy array containing the direction of the image stack inclination: Optional inclination of the image in degrees. If none is set, an inclination of 0° is assumed. saturation: Weight image by using the saturation value. Use either a 2D image or a 3D image with the same shape as the direction. If no image is used, the saturation for all image pixels will be set to 1 value: Weight image by using the value. Use either a 2D image or a 3D image with the same shape as the direction. If no image is used, the value for all image pixels will be set to 1 colormap: The colormap to use. Default is HSV black. The available color maps can be found in the colormap class. Returns: numpy.ndarray: 2D image containing the resulting HSV orientation map """ direction = numpy.array(direction) direction_shape = direction.shape if inclination is None: inclination = numpy.zeros_like(direction) colors = colormap(direction, inclination) hsv_colors = rgb_to_hsv(colors) # If no saturation is given, create an "empty" saturation image that will be used if saturation is None: saturation = numpy.ones(direction.shape) # Normalize saturation image saturation = saturation / saturation.max(axis=None) # If we have a saturation image, check if the shape matches (3D) and correct accordingly while len(saturation.shape) < len(direction.shape): saturation = saturation[..., numpy.newaxis] if not saturation.shape[-1] == direction_shape[-1]: saturation = numpy.repeat(saturation, direction_shape[-1], axis=-1) # If no value is given, create an "empty" value image that will be used if value is None: value = numpy.ones(direction.shape) # Normalize value image value = value / value.max(axis=None) # If we have a value image, check if the shape matches (3D) and correct accordingly while len(value.shape) < len(direction.shape): value = value[..., numpy.newaxis] if not value.shape[-1] == direction_shape[-1]: value = numpy.repeat(value, direction_shape[-1], axis=-1) hsv_colors[..., 1] *= saturation hsv_colors[..., 2] *= value colors = hsv_to_rgb(hsv_colors) if len(direction_shape) > 2: return (255.0 * _visualize_multiple_direction(direction, colors)).astype(numpy.uint8) else: return (255.0 * _visualize_one_direction(direction, colors)).astype(numpy.uint8)
View Source
from functools import partial import numpy from multiprocessing import Pool, current_process from multiprocessing.sharedctypes import RawArray import os import scipy.signal as signal from SLIX._preparation import _thin_out_median, _thin_out_plain, \ _thin_out_average, _init_worker_fourier_smoothing, \ _worker_function_fourier_smoothing, _fourier_smoothing __all__ = ['thin_out', 'savitzky_golay_smoothing', 'low_pass_fourier_smoothing'] def low_pass_fourier_smoothing(image, threshold=0.2, smoothing_factor=0.025): """ Applies Low Pass fourier filter to given line profiles / image and returns the smoothened measurement. Args: image: Complete SLI measurement image stack as a 2D/3D Numpy array threshold: Threshold percentage of low frequencies which will completely pass smoothing_factor: Apply a smoothing factor which will smooth out the applied multiplication factor for the low pass filter. A higher value will result in more smoothing of the curve. Values between 1e-15 and 1 are accepted. Other values might result in an error. Returns: Complete SLI measurement image with applied Low Pass fourier filter and the same shape as the original image. """ if not current_process().name == 'MainProcess': x_shape = image.shape x = RawArray('d', x_shape[0] * x_shape[1] * x_shape[2]) x_np = numpy.frombuffer(x).reshape(x_shape) numpy.copyto(x_np, image) partial_worker_function = partial(_worker_function_fourier_smoothing, threshold=threshold, window=smoothing_factor) with Pool(processes=os.cpu_count(), initializer=_init_worker_fourier_smoothing, initargs=(x, x_shape)) as pool: pool.map(partial_worker_function, range(x_shape[0] * x_shape[1])) return x_np.astype(image.dtype) image = _fourier_smoothing(image, threshold, smoothing_factor) return image def savitzky_golay_smoothing(image, window_length=45, polyorder=2): """ Applies Savitzky-Golay filter to given line profiles / image and returns the smoothened measurement. Args: image: Complete SLI measurement image stack as a 2D/3D Numpy array window_length: Used window length for the Savitzky-Golay filter polyorder: Used polynomial order for the Savitzky-Golay filter Returns: Complete SLI measurement image with applied Savitzky-Golay filter and the same shape as the original image. """ conc_image = numpy.concatenate((image[:, :, -window_length:], image, image[:, :, :window_length]), axis=-1) conc_image = signal.savgol_filter(conc_image, window_length, polyorder, axis=-1) return conc_image[:, :, window_length:-window_length] def thin_out(image, factor=2, strategy='plain'): """ Thin out the image stack used for SLIX. This can be useful when the image stack is quite large and should be processed quickly. This can also prove useful if there is a lot of noise that could be filtered by using a lower resolution image. Args: image: Image that should be thinned out. factor: Factor which will be used for thinning the image. A factor of N means that every N-th pixel will be kept. strategy: Strategy used for thinning out the image. Available methods: 'plain' (keep the pixel), 'average' (calculate average of area), 'median' (calculate the median of area) Returns: numpy.ndarray with the thinned out image """ strategy = strategy.lower() if strategy == 'plain': return _thin_out_plain(image, factor) elif strategy == 'average': return _thin_out_average(image, factor) elif strategy == 'median': return _thin_out_median(image, factor) else: raise ValueError('Strategy not implemented. Known strategies are:' ' plain, average, median.')
View Source
import hashlib import typing import h5py import numpy import secrets import getpass from SLIX._logging import get_logger __all__ = ['AttributeHandler'] class AttributeHandler: def __init__(self, dataset: h5py.Dataset): """ Initialize the AttributeHandler with a already opened HDF5 dataset. This dataset will be used for all operations of this class. Args: dataset: h5py dataset """ self.dataset: h5py.Dataset = dataset self.attrs: h5py.AttributeManager = dataset.attrs self.logger = get_logger(__name__) def does_attribute_exist(self, attribute_name: str) -> bool: """ Check if the attribute already exists in the HDF5 dataset. This has to be done before doing any operations because writing to an HDF5 attribute without properly deleting it first can result in errors. Args: attribute_name: Name of the attribute you want to check in the dataset. Returns: None """ attribute_names: typing.AbstractSet[str] = self.attrs.keys() return attribute_name in attribute_names def delete_attribute(self, attribute_name: str) -> None: """ Delete an attribute from a HDF5 dataset. Args: attribute_name: Name of the attribute you want to delete in the dataset. Returns: None """ if self.does_attribute_exist(attribute_name): del self.attrs[attribute_name] else: self.logger.warning("Attribute %s was not deleted as it does not " "exist" % {attribute_name}) def get_attribute(self, attribute_name: str) -> \ typing.Union[str, float, int, bool, numpy.array]: """ Get an attribute from the HDF5 dataset. Args: attribute_name: Name of the attribute you want to get from the dataset. Returns: Value from the dataset (string, float, int, bool or numpy.array) if the attribute es present. Otherwise None will be returned. """ if not self.does_attribute_exist(attribute_name): self.logger.error('Attribute %s does not exist!' % {attribute_name}) return None return self.attrs[attribute_name] def set_attribute(self, attribute_name: str, value: typing.Union[str, float, int, bool, numpy.array]) \ -> None: """ Set an attribute in the HDF5 dataset. Args: attribute_name: Name of the attribute you want to get from the dataset. value: String, Float, Integer, Boolean or numpy array you want to set in the HDF5 attribute. Returns: None """ if self.does_attribute_exist(attribute_name): self.delete_attribute(attribute_name) self.attrs.create(attribute_name, value) def set_reference_modality_to(self, reference: "AttributeHandler") -> None: """ When SLIX generates an image based on a SLI measurement, the original HDF5 file can be saves as a reference for the future. This method adds the reference file and dataset to the output. However, the input HDF5 and dataset must contain an id attribute. Args: reference: Reference AttributeHandler containing the dataset of the input file. Returns: None """ self.set_reference_modalities_to([reference]) def set_reference_modalities_to(self, references: typing.List["AttributeHandler"]) \ -> None: """ When SLIX generates an image based on a SLI measurement, the original HDF5 file can be saves as a reference for the future. This method adds the reference file and dataset to the output. However, the input HDF5 and dataset must contain an id attribute. Args: references: Reference list of AttributeHandlers containing the dataset of the input file. Returns: None """ ref_id: typing.List[str] = [] for reference in references: if reference.does_attribute_exist('id'): ref_id.append(reference.get_attribute('id')) else: self.logger.error('Could not set reference_images because id is not ' 'present in at least one dataset') return self.set_attribute('reference_images', ref_id) def add_creator(self) -> None: """ Adds the creator of the HDF5 file to the dataset. Returns: None """ creator: str = getpass.getuser() self.set_attribute('created_by', creator) def add_id(self) -> None: """ Computes a unique ID that will be added to the dataset of the HDF5 file. The ID is a sha256 hash containing some of the attributes as well as a 50 digit randomized prefix. Returns: None """ hashstr: str = secrets.token_hex(50) attributes: list = ['brain_id', 'brain_part_id', 'section_id', 'image_modality', 'creation_time', 'created_by', 'software', 'software_revision', 'software_parameters'] for attribute in attributes: hashstr: str = hashstr + attribute hash_obj: hashlib.sha256 = hashlib.sha256() hash_obj.update(hashstr.encode('ascii')) self.set_attribute('id', hash_obj.hexdigest()) def copy_all_attributes_to(self, dest: "AttributeHandler", exceptions: typing.List[str] = None) -> None: """ Copies all attributes from one AttributeHandler to another. Exceptions can be given in a list. In general, the following attributes will not be copied: "created_by", "creation_time", "id", "image_modality", "reference_images", "software", "software_revision", "software_parameters", "filename", "path", "scale" Args: dest: Destination where the attributes of this handler should be copied to. exceptions: Exceptions in form of a list with strings. Those attributes will not be copied when calling the method. Returns: None """ if exceptions is None: exceptions = [] fixed_exceptions: typing.List[str] = [ # Attributes that MUST be manually set "created_by", "creation_time", "dashboard_id", "id", "image_modality", "reference_images", "software", "software_revision", "software_parameters", # Attributes set by DB software "filename", "path", "scale" ] attribute_names: typing.AbstractSet[str] = self.attrs.keys() for attribute_name in attribute_names: if attribute_name not in fixed_exceptions \ and attribute_name not in exceptions: self.copy_attribute_to(dest, attribute_name) def copy_attributes_to(self, dest: "AttributeHandler", attributes: typing.List[str] = None) -> None: """ Copies given attributes from one AttributeHandler to another. Args: dest: Destination where the attributes of this handler should be copied to. attributes: Attributes as a list of strings which will be copied. Returns: None """ if attributes is None: return for attribute in attributes: self.copy_attribute_to(dest, attribute) def copy_attribute_to(self, dest: "AttributeHandler", attribute_name: str) -> None: """ Copy a single attribute from one AttributeHandler to another. Args: dest: Destination where the attributes of this handler should be copied to. attribute_name: Attribute name which will be copied. Returns: None """ dest.set_attribute(attribute_name, self.get_attribute(attribute_name))
View Source
import numpy def full_mask(high_prominence_peaks, low_prominence_peaks, peakdistance, max_image): """ This method classifies the resulting parameter maps of an SLI measurement to identify the areas of the image. Based on the other methods in this module, a mask combining the results of the other methods is created. The resulting mask will therefore contain values to separate flat, crossing and inclined fibers. ---------- The resulting mask is a binary mask with the following values: 0: The area neither flat, crossing or inclined. 1: The area is a flat fiber. 2: The area is contains two crossing fibers. 3: The area is contains three crossing fibers. 4: The area is lightly inclined. 5: The area is inclined. 6: The area is strongly inclined. Args: high_prominence_peaks: numpy.ndarray containing the number of peaks with a high prominence. low_prominence_peaks: numpy.ndarray containing the number of peaks with a low prominence. peakdistance: numpy.ndarray containing the distance between the peaks. max_image: numpy.ndarray containing the maximum signal of the image. Returns: numpy.ndarray containing the binary mask. """ crossing = crossing_mask(high_prominence_peaks, max_image) flat = flat_mask(high_prominence_peaks, low_prominence_peaks, peakdistance) inclined = inclinated_mask(high_prominence_peaks, peakdistance, max_image, flat) return_mask = flat.astype(numpy.uint8) return_mask[crossing == 1] = 2 return_mask[crossing == 2] = 3 return_mask[inclined == 2] = 4 return_mask[inclined == 3] = 5 return_mask[inclined == 4] = 6 return return_mask def crossing_mask(high_prominence_peaks, max_image): """ This method classifies the resulting parameter maps of an SLI measurement to identify the areas of the image where the underlying fiber structure is probably a crossing one. To do so, the method uses the following criteria: 1. The maximum signal during the measurement is above the mean signal of the maximum image. 2. The number of peaks is either four (two crossing fibers) or six (three crossing fibers). -------------- The resulting mask is a binary mask with the following values: 0: The area is not a crossing one. 1: The area is a crossing one with two crossing fibers. 2: The area is a crossing one with three crossing fibers. Args: high_prominence_peaks: numpy.ndarray containing the number of peaks with a high prominence. max_image: numpy.ndarray containing the maximum signal of the image. Returns: numpy.ndarray containing the binary mask. """ crossing = (high_prominence_peaks == 4) | (high_prominence_peaks == 6) crossing = crossing.astype(numpy.uint8) mean_flat_signal = numpy.average(max_image) crossing[(max_image < mean_flat_signal)] = 0 crossing[high_prominence_peaks == 4] = 1 crossing[high_prominence_peaks == 6] = 2 return crossing def inclinated_mask(high_prominence_peaks, peakdistance, max_image, flat_mask): """ This method classifies the resulting parameter maps of an SLI measurement to identify the areas of the image where the underlying fiber structure is probably an inclined one. To do so, the method uses the following criteria: Flat fibers: 1. The maximum signal during the measurement is above the mean signal of the maximum image. 2. The number of peaks is two (one flat fiber) Inclined fibers: Three different scenarios are possible: 1. Two peaks are present and the peak distance is between 120° and 150° (lightly inclined fiber) 2. Two peaks are present and the peak distance is below 120° (inclined fiber) 3. One single peak is present (steep fiber) ---------- The resulting mask is a binary mask with the following values: 0: The area is neither a flat nor an inclined one. 1: The area is a flat one. 2: The area is a lightly inclined one. 3: The area is an inclined one. 4: The area is a steep one. Args: high_prominence_peaks: numpy.ndarray containing the number of peaks with a high prominence. peakdistance: Mean distance between the prominent peaks in high_prominence_peaks. max_image: numpy.ndarray containing the maximum signal of the image. flat_mask: numpy.ndarray containing the binary mask of the flat areas. Returns: numpy.ndarray containing the binary mask. """ inclinated_areas = numpy.zeros(high_prominence_peaks.shape, dtype=numpy.uint8) mean_flat_signal = numpy.average(max_image[flat_mask > 0]) two_peak_mask = (high_prominence_peaks == 2) & \ (max_image > mean_flat_signal) inclinated_areas[flat_mask > 0] = 1 inclinated_areas[two_peak_mask & (peakdistance > 120) & (peakdistance < 150)] = 2 inclinated_areas[two_peak_mask & (peakdistance < 120)] = 3 inclinated_areas[high_prominence_peaks == 1] = 4 return inclinated_areas def flat_mask(high_prominence_peaks, low_prominence_peaks, peakdistance): """ This method classifies the resulting parameter maps of an SLI measurement to identify the areas of the image where the underlying fiber structure is probably a flat one. To do so, the method uses the following criteria: 1. Two prominent peaks are present. 2. The peak distance is between 145° and 215°. A peak distance of 180° is expected for a completely flat fiber, but small deviations for example through the sampling steps of the measurement are possible. 3. No more than two low prominent peaks are present. Completely flat fibers generally have a very stable signal and therefore a low amount of low prominent peaks. ---------- The resulting mask is a binary mask with the following values: 0: The area is not a flat fiber. 1: The area is a flat fiber. Args: high_prominence_peaks: numpy.ndarray containing the number of peaks low_prominence_peaks: numpy.ndarray containing the number of peaks which are below the threshold of prominence needed to be considered as a prominent peak. peakdistance: Mean distance between the prominent peaks in high_prominence_peaks. Returns: numpy.ndarray containing the binary mask. """ return (high_prominence_peaks == 2) & (low_prominence_peaks < 2) & \ (peakdistance > 145) & (peakdistance < 215)