diff --git a/nornir_imageregistration/core.py b/nornir_imageregistration/core.py index 3f2f693..afaee6b 100644 --- a/nornir_imageregistration/core.py +++ b/nornir_imageregistration/core.py @@ -3,6 +3,7 @@ """ from collections.abc import Iterable +from logging import exception import math from multiprocessing import shared_memory from multiprocessing.shared_memory import SharedMemory @@ -215,12 +216,15 @@ def remove_duplicate_points(points: NDArray, columns=list[int]) -> tuple[NDArray return sorted_point_pairs -def ReduceImage(image: NDArray, scalar: float) -> NDArray: +def ScaleImage(image: NDArray, scalar: float) -> NDArray: """ - Returns a zoomed array using spline interpolation (CPU/GPU agnostic function) + Returns a scaled array using spline interpolation (CPU/GPU agnostic function) """ xp = cupyx.scipy.get_array_module(image) - return xp.ndimage.zoom(image, scalar) + if nornir_imageregistration.UsingCupy(): + return xp.ndimage.zoom(image, scalar) + else: + return xp.ndimage.zoom(image.astype(np.float32), scalar) def ExtractROI(image: NDArray, center, area) -> NDArray: @@ -666,7 +670,11 @@ def close_shared_memory(input: nornir_imageregistration.Shared_Mem_Metadata | Sh if input.shared_memory is not None: input.shared_memory.close() elif isinstance(input, SharedMemory): - input.close() + try: + input.close() + except Exception as e: + prettyoutput.LogErr(f"Error closing shared memory {input.name}\n{e}") + return # if input.name in __known_shared_memory_allocations: # shared_mem, finalizer = __known_shared_memory_allocations[input.name] @@ -1078,7 +1086,7 @@ def LoadImage(ImageFullPath: str, if not MaxDimension is None: scalar = ScalarForMaxDimension(MaxDimension, image.shape) if scalar < 1.0: - image = ReduceImage(image, scalar) + image = ScaleImage(image, scalar) image_mask = None @@ -1094,7 +1102,7 @@ def LoadImage(ImageFullPath: str, if MaxDimension is not None: scalar = ScalarForMaxDimension(MaxDimension, image_mask.shape) if scalar < 1.0: - image_mask = ReduceImage(image_mask, scalar) + image_mask = ScaleImage(image_mask, scalar) assert (image.shape == image_mask.shape) image = RandomNoiseMask(image, image_mask) @@ -1286,10 +1294,12 @@ def CreateExtremaMask(image: np.ndarray, mask: np.ndarray = None, size_cutoff=0. """ Returns a mask for features above a set size that are at max or min pixel value :param image: + :param mask: Masked regions are excluded from the analysis, the min/max values are calculated from unmasked pixels only :param minima: :param maxima: :param numpy.ndarray mask: Pixels we wish to not include in the analysis :param size_cutoff: Determines how large a continuous region must be before it is masked. If 0 to 1 this is a fraction of total area. If > 1 it is an absolute count of pixels. If None all min/max are masked regardless of size + :returns: Mask of extrema pixels, pixels that are FALSE are extrema to be excluded """ # (minima, maxima, iMin, iMax) = scipy.ndimage.measurements.extrema(image) @@ -1307,6 +1317,9 @@ def CreateExtremaMask(image: np.ndarray, mask: np.ndarray = None, size_cutoff=0. if maxima is None: maxima = image.max() + # Pixels that are TRUE will be excluded, exclude pixels equal to the min or max. + # However, the ndimage.label function finds features that are TRUE. So we start with an + # inverted mask extrema_mask = xp.logical_or(image == maxima, image == minima) if mask is not None: @@ -1319,6 +1332,8 @@ def CreateExtremaMask(image: np.ndarray, mask: np.ndarray = None, size_cutoff=0. if nLabels == 0: # If there are no labels, do not mask anything return xp.ones(image.shape, extrema_mask.dtype) + # Identify the label of non-extrema pixels + label_sums = sp.ndimage.sum_labels( extrema_mask.astype(np.int32) if nornir_imageregistration.UsingCupy() else extrema_mask, extrema_mask_label, xp.array(range(0, nLabels))) @@ -1334,16 +1349,16 @@ def CreateExtremaMask(image: np.ndarray, mask: np.ndarray = None, size_cutoff=0. else: cutoff_value = size_cutoff - labels_to_save = label_sums < cutoff_value - if xp.any(labels_to_save): - cutoff_labels = xp.flatnonzero(labels_to_save) + small_regions = label_sums < cutoff_value + if xp.any(small_regions): + cutoff_labels = xp.flatnonzero(small_regions) extrema_mask_minus_small_features = xp.isin(extrema_mask_label, cutoff_labels) # nornir_imageregistration.ShowGrayscale((image, extrema_mask, extrema_mask_minus_small_features)) return extrema_mask_minus_small_features else: - raise NotImplemented() + return np.ones(image.shape, bool) # No features large enough to exclude, retain the entire image def ReplaceImageExtremaWithNoise(image: np.ndarray, imagemask: np.ndarray = None, diff --git a/nornir_imageregistration/image_permutation_helper.py b/nornir_imageregistration/image_permutation_helper.py index c4174d6..f0a24af 100644 --- a/nornir_imageregistration/image_permutation_helper.py +++ b/nornir_imageregistration/image_permutation_helper.py @@ -91,7 +91,7 @@ def __init__(self, self._extrema_size_cutoff_in_pixels = None if extrema_mask_size_cuttoff is None: - extrema_mask_size_cuttoff = np.array((128, 128)) + extrema_mask_size_cuttoff = 0.01 if isinstance(extrema_mask_size_cuttoff, np.ndarray): self.extrema_size_cutoff_in_pixels = int(np.prod(extrema_mask_size_cuttoff)) diff --git a/nornir_imageregistration/settings/__init__.py b/nornir_imageregistration/settings/__init__.py index 5a2f844..b55adc4 100644 --- a/nornir_imageregistration/settings/__init__.py +++ b/nornir_imageregistration/settings/__init__.py @@ -3,6 +3,8 @@ from .grid_refinement import GridRefinement from .mosaic_tile_offset import LoadMosaicOffsets, SaveMosaicOffsets, TileOffset from .translate import TranslateSettings +from .angle_range import AngleSearchRange +from .stos_brute import StosBruteSettings def GetOrSaveTranslateSettings(settings: TranslateSettings, path: str): diff --git a/nornir_imageregistration/settings/angle_range.py b/nornir_imageregistration/settings/angle_range.py new file mode 100644 index 0000000..83e7504 --- /dev/null +++ b/nornir_imageregistration/settings/angle_range.py @@ -0,0 +1,24 @@ +from numpy.typing import NDArray +import numpy as np +from pydantic import BaseModel +from math import pi + + +class AngleSearchRange(BaseModel): + max_angle: float | None = None # Maximum +/- deflection angle to rotate images when searching for the best control point alignment, if None, a full circle is searched at step-size intervals + angle_step_size: float = 3 # Number of degrees to step between search angles + + @property + def angle_range(self) -> NDArray[float]: + if self.max_angle is None: + angles = np.arange(start=-180, stop=180, step=self.angle_step_size) + else: + angles = np.arange(start=-self.max_angle, + stop=self.max_angle + self.angle_step_size, + step=self.angle_step_size) # numpy.linspace(-7.5, 7.5, 11) + + angles = np.union1d(angles, [0]) + return angles + + def __iter__(self): + return iter(self.angle_range) diff --git a/nornir_imageregistration/settings/grid_refinement.py b/nornir_imageregistration/settings/grid_refinement.py index e8ee2b7..ebe3110 100644 --- a/nornir_imageregistration/settings/grid_refinement.py +++ b/nornir_imageregistration/settings/grid_refinement.py @@ -4,6 +4,7 @@ @author: u0490822 """ from __future__ import annotations +from typing import Iterable, Sequence import numpy as np from numpy.typing import NDArray @@ -11,7 +12,7 @@ import nornir_imageregistration -class GridRefinement(object): +class GridRefinement: """ Settings for grid refinement """ @@ -81,10 +82,10 @@ def __init__(self, target_mask: NDArray[np.bool_] | None = None, source_mask: NDArray[np.bool_] | None = None, num_iterations: int = None, - cell_size=None, - grid_spacing=None, - angles_to_search=None, - final_pass_angles=None, + cell_size: int | NDArray[int] | Iterable[int] | None = None, + grid_spacing: int | NDArray[int] | Iterable[int] | None = None, + angles_to_search: Iterable[float] | NDArray[float] | None = None, + final_pass_angles: Iterable[float] | NDArray[float] | None = None, max_travel_for_finalization: float = None, max_travel_for_finalization_improvement: float = None, min_alignment_overlap: float = None, @@ -112,7 +113,7 @@ def __init__(self, :param bool cupy_processing: True if the refinement will be done on the GPU. When set, arrays are created as cupy arrays instead of NDArrays """ - self._single_thread_processing = single_thread_processing or nornir_imageregistration.UsingCupy() + self._single_thread_processing = single_thread_processing or nornir_imageregistration.UsingCupy() if target_image is None: raise ValueError("target_image must be specified") @@ -177,10 +178,10 @@ def __init__(self, def CreateWithPreprocessedImages(target_img_data: nornir_imageregistration.ImagePermutationHelper, source_img_data: nornir_imageregistration.ImagePermutationHelper, num_iterations: int = None, - cell_size=None, - grid_spacing=None, - angles_to_search=None, - final_pass_angles=None, + cell_size: int | NDArray[int] | Iterable[int] | None = None, + grid_spacing: int | NDArray[int] | Iterable[int] | None = None, + angles_to_search: Iterable[float] | NDArray[float] | None = None, + final_pass_angles: Iterable[float] | NDArray[float] | None = None, max_travel_for_finalization: float = None, max_travel_for_finalization_improvement: float = None, min_alignment_overlap: float = None, diff --git a/nornir_imageregistration/settings/stos_brute.py b/nornir_imageregistration/settings/stos_brute.py new file mode 100644 index 0000000..728caf1 --- /dev/null +++ b/nornir_imageregistration/settings/stos_brute.py @@ -0,0 +1,73 @@ +from numpy.typing import NDArray +import numpy as np +from typing import NamedTuple, Sequence, Iterable +from pydantic import BaseModel +from nornir_imageregistration.settings.angle_range import AngleSearchRange + + +class StosBruteSettings(BaseModel): + """Encodes the settings required or used to invoke StosBrute""" + angles: AngleSearchRange | Sequence[float] | None = None + min_overlap: float = 0.75 # The minimum amount of overlap we require in the images. Higher values reduce false positives but may not register offset images + source_image_scale_factors: tuple[float] | None = None + """Amount to scale the warped image before attempting registration, + this handles cases where multiple scopes are used with slightly differnt magnification values """ + + larget_dimension: int | None = None # The input images should be scaled so the largest image dimension is equal to this value, default is 1024. None means use the actual image size + try_flipped: bool = False # If True the algorithm will test the flipped version of the source image too + + def __init__(self, + angles: AngleSearchRange | Sequence[float] | None = None, + min_overlap: float = 0.75, + source_image_scale_factors: NDArray[float] | None = None, + larget_dimension: int | None = 1024, + try_flipped: bool = False): + """ + + :param angles: Angles to search for the best control point alignment or None if all angles should be searched + :param min_overlap: Minimum amount of overlap we require to consider a registration to be valiud + :param source_image_scale_factors: Amount to scale the warped image before attempting registration, this handles cases where multiple scopes are used with slightly differnt magnification values + :param larget_dimension: The input images should be scaled so the largest image dimension is equal to this value, default is 1024. None means use the actual image size + :param try_flipped: If True the algorithm will test the flipped version of the source image too + """ + super().__init__() + self.angles = angles + self.min_overlap = min_overlap + self.larget_dimension = larget_dimension + self.try_flipped = try_flipped + + self.source_image_scale_factors = source_image_scale_factors + + if self.source_image_scale_factors is None: + pass + elif isinstance(source_image_scale_factors, np.ndarray): + self.source_image_scale_factors = source_image_scale_factors + elif hasattr(source_image_scale_factors, '__iter__'): + self.source_image_scale_factors = np.array(source_image_scale_factors, float) + else: + self.source_image_scale_factors = np.array([source_image_scale_factors, source_image_scale_factors], float) + + def angle_range_defined(self) -> bool: + return self.angles is not None + + @property + def angle_range(self) -> NDArray[float]: + """:return: The range of angles to search for the best control point alignment or None if all angles should be searched""" + if self.angles is None: + return np.array(range(-178, 182, 2), float) + + if isinstance(self.angles, np.ndarray): + return self.angles + elif isinstance(self.angles, AngleSearchRange): + return self.angles.angle_range + elif isinstance(self.angles, Iterable): + return np.array(list(self.angles), float) + + raise ValueError(f"Unexpected type for self.angle_search_settings: {self.angles.__class__}") + + @property + def source_image_scaling_required(self) -> bool: + if self.source_image_scale_factors is None: + return False + + return np.any(self.source_image_scale_factors != 1) diff --git a/nornir_imageregistration/stos_brute.py b/nornir_imageregistration/stos_brute.py index e3057cf..66d3998 100644 --- a/nornir_imageregistration/stos_brute.py +++ b/nornir_imageregistration/stos_brute.py @@ -3,6 +3,7 @@ @author: u0490822 ''' +from typing import NamedTuple import multiprocessing import multiprocessing.sharedctypes from time import sleep @@ -11,6 +12,9 @@ from typing import Sequence import logging +from nornir_imageregistration import AlignmentRecord +from nornir_imageregistration.settings import StosBruteSettings, AngleSearchRange + # Check if cupy is available, and if it is not import thunks that refer to scipy/numpy try: import cupy as cp @@ -60,25 +64,32 @@ def SliceToSliceBruteForce(FixedImageInput: nornir_imageregistration.ImageLike, AngleSearchRange = set(AngleSearchRange) # if isinstance(AngleSearchRange, np.ndarray): - if 0 not in AngleSearchRange: - logger = logging.getLogger(__name__ + '.SliceToSliceBruteForce') - logger.warning("AngleSearchRange should contain 0 degrees to ensure the best match is found") + if 0 not in set(AngleSearchRange): + logger = logging.getLogger(__name__ + '.SliceToSliceBruteForce') + logger.warning("AngleSearchRange should contain 0 degrees to ensure the best match is found") SingleThread = True if use_cp else SingleThread - WarpedImageScalingRequired = False - if WarpedImageScaleFactors is not None: - if hasattr(WarpedImageScaleFactors, '__iter__'): - WarpedImageScaleFactors = nornir_imageregistration.EnsurePointsAre1DNumpyArray(WarpedImageScaleFactors) - WarpedImageScalingRequired = any(WarpedImageScaleFactors != 1) - else: - WarpedImageScalingRequired = WarpedImageScaleFactors != 1 - WarpedImageScaleFactors = nornir_imageregistration.EnsurePointsAre1DNumpyArray( - [WarpedImageScaleFactors, WarpedImageScaleFactors]) - target_image_data = nornir_imageregistration.ImagePermutationHelper(FixedImageInput, FixedImageMaskPath) source_image_data = nornir_imageregistration.ImagePermutationHelper(WarpedImageInput, WarpedImageMaskPath) + settings = StosBruteSettings(angles=AngleSearchRange, + min_overlap=MinOverlap, + source_image_scale_factors=WarpedImageScaleFactors, + larget_dimension=LargestDimension, + try_flipped=TestFlip) + + return SliceToSliceBruteForceWithPreprocessedImages(source_image_data, target_image_data, settings, + SingleThread=SingleThread, Cluster=Cluster) + + +def SliceToSliceBruteForceWithPreprocessedImages(source_image_data: nornir_imageregistration.ImagePermutationHelper, + target_image_data: nornir_imageregistration.ImagePermutationHelper, + settings: StosBruteSettings, + SingleThread: bool = False, + Cluster: bool = False) -> nornir_imageregistration.AlignmentRecord: + use_cp = nornir_imageregistration.GetActiveComputationLib() == nornir_imageregistration.ComputationLib.cupy + target_image = target_image_data.ImageWithMaskAsNoise source_image = source_image_data.ImageWithMaskAsNoise @@ -88,71 +99,72 @@ def SliceToSliceBruteForce(FixedImageInput: nornir_imageregistration.ImageLike, del target_image_data del source_image_data - target_image = cp.asarray(target_image) if use_cp and not isinstance(target_image, cp.ndarray) else target_image + target_image = cp.asarraye(target_image) if use_cp and not isinstance(target_image, cp.ndarray) else target_image source_image = cp.asarray(source_image) if use_cp and not isinstance(source_image, cp.ndarray) else source_image scalar = 1.0 - if LargestDimension is not None: - scalar = nornir_imageregistration.ScalarForMaxDimension(LargestDimension, + if settings.larget_dimension is not None: + scalar = nornir_imageregistration.ScalarForMaxDimension(settings.larget_dimension, [target_image.shape, source_image.shape]) - if scalar < 1.0: - target_image = nornir_imageregistration.ReduceImage(target_image, scalar) - source_image = nornir_imageregistration.ReduceImage(source_image, scalar) + if scalar != 1.0: + target_image = nornir_imageregistration.ScaleImage(target_image, scalar) + source_image = nornir_imageregistration.ScaleImage(source_image, scalar) # Replace extrema with noise - UserDefinedAngleSearchRange = AngleSearchRange is not None - if not UserDefinedAngleSearchRange: - AngleSearchRange = list(range(-178, 182, 2)) - - BestMatch = _find_best_angle(target_image, source_image, - target_stats, source_stats, - AngleSearchRange, MinOverlap=MinOverlap, SingleThread=SingleThread, - use_cluster=Cluster) - - IsFlipped = False - if TestFlip: - # imWarpedFlipped = np.copy(source_image) - imWarpedFlipped = np.flipud(source_image) - - BestMatchFlipped = _find_best_angle(target_image, imWarpedFlipped, - target_stats, source_stats, - AngleSearchRange, MinOverlap=MinOverlap, - SingleThread=SingleThread, use_cluster=Cluster) - BestMatchFlipped.flippedud = True + + best_match = _find_best_angle(source_image=source_image, target_image=target_image, + source_stats=source_stats, target_stats=target_stats, + angle_range=settings.angle_range, + min_overlap=settings.min_overlap, + SingleThread=SingleThread, + use_cluster=Cluster) + + is_flipped = False + if settings.try_flipped: + # source_flipped = np.copy(source_image) + source_flipped = np.flipud(source_image) + + best_match_flipped = _find_best_angle(source_image=source_flipped, target_image=target_image, + source_stats=source_stats, target_stats=target_stats, + angle_range=settings.angle_range, + min_overlap=settings.min_overlap, + SingleThread=SingleThread, use_cluster=Cluster) + best_match_flipped.flippedud = True # Determine if the best match is flipped or not - IsFlipped = BestMatchFlipped.weight > BestMatch.weight + is_flipped = best_match_flipped.weight > best_match.weight - if IsFlipped: - imWarped = imWarpedFlipped - BestMatch = BestMatchFlipped + if is_flipped: + source_image = source_flipped + best_match = best_match_flipped else: - imWarped = source_image + source_image = source_image # Note Clement - the RefinedAngleSearch list below is not centered around the current best angle # Default angle search range every 2 degrees - # Old RefinedAngleSearch list: [(x * 0.1) + BestMatch.angle - 1.9 for x in range(0, 18)] - # New RefinedAngleSearch list (length 39): [(x * 0.1 + BestMatch.angle) for x in range(-19, 20)] - # New optional RefinedAngleSearch list (length 18): [(x * 0.2 + BestMatch.angle) for x in range(-9, 10)] - if not UserDefinedAngleSearchRange: - BestRefinedMatch = _find_best_angle(target_image, imWarped, - target_stats, source_stats, - # [(x * 0.1) + BestMatch.angle - 1.9 for x in range(0, 18)], - [(x * 0.1 + BestMatch.angle) for x in range(-19, 20)], - MinOverlap=MinOverlap, SingleThread=SingleThread) - BestRefinedMatch.flippedud = IsFlipped + # Old RefinedAngleSearch list: [(x * 0.1) + best_match.angle - 1.9 for x in range(0, 18)] + # New RefinedAngleSearch list (length 39): [(x * 0.1 + best_match.angle) for x in range(-19, 20)] + # New optional RefinedAngleSearch list (length 18): [(x * 0.2 + best_match.angle) for x in range(-9, 10)] + if not settings.angle_range_defined(): + best_refined_match = _find_best_angle(source_image=source_image, target_image=target_image, + source_stats=source_stats, target_stats=target_stats, + # [(x * 0.1) + best_match.angle - 1.9 for x in range(0, 18)], + angle_range=[(x * 0.2 + best_match.angle) for x in range(-9, 10)], + min_overlap=settings.min_overlap, SingleThread=SingleThread) + best_refined_match.flippedud = is_flipped else: min_step_size = 0.25 - if len(AngleSearchRange) > 2: - iMatch = AngleSearchRange.index(BestMatch.angle) - iBelow = iMatch - 1 if iMatch - 1 >= 0 else len(AngleSearchRange) - 1 - iAbove = iMatch + 1 if iMatch + 1 < len(AngleSearchRange) else 0 - below = AngleSearchRange[iMatch - 1] if iMatch - 1 >= 0 else AngleSearchRange[0] - np.abs( - AngleSearchRange[1] - AngleSearchRange[0]) - above = AngleSearchRange[iMatch + 1] if iMatch + 1 < len(AngleSearchRange) else AngleSearchRange[ - iMatch] + np.abs( - AngleSearchRange[iMatch] - AngleSearchRange[iMatch - 1]) + if len(settings.angle_range) > 2: + sorted_angles = sorted(settings.angle_range) + iMatch = sorted_angles.index(best_match.angle) + iBelow = iMatch - 1 if iMatch - 1 >= 0 else len(sorted_angles) - 1 + iAbove = iMatch + 1 if iMatch + 1 < len(sorted_angles) else 0 + below = sorted_angles[iMatch - 1] if iMatch - 1 >= 0 else sorted_angles[0] - np.abs( + sorted_angles[1] - sorted_angles[0]) + above = sorted_angles[iMatch + 1] if iMatch + 1 < len(sorted_angles) else sorted_angles[ + iMatch] + np.abs( + sorted_angles[iMatch] - sorted_angles[iMatch - 1]) refine_search_range = above - below nSteps = 20 stepsize = refine_search_range / nSteps @@ -161,162 +173,178 @@ def SliceToSliceBruteForce(FixedImageInput: nornir_imageregistration.ImageLike, nSteps = int(refine_search_range / min_step_size) stepsize = refine_search_range / nSteps - BestRefinedMatch = _find_best_angle(target_image, imWarped, - target_stats, source_stats, - [(x * stepsize) + below for x in range(1, nSteps)], - MinOverlap=MinOverlap, SingleThread=SingleThread) - BestRefinedMatch.flippedud = IsFlipped + refined_angle_search_range = {(x * stepsize) + below for x in range(1, nSteps)} + + # Ensure we include the best match angle + refined_angle_search_range.add(best_match.angle) + + best_refined_match = _find_best_angle(source_image=source_image, target_image=target_image, + source_stats=source_stats, target_stats=target_stats, + angle_range=np.array(list(refined_angle_search_range), float), + min_overlap=settings.min_overlap, SingleThread=SingleThread) + best_refined_match.flippedud = is_flipped else: - BestRefinedMatch = BestMatch - BestRefinedMatch.flippedud = IsFlipped + best_refined_match = best_match + best_refined_match.flippedud = is_flipped if scalar > 1.0: - AdjustedPeak = (BestRefinedMatch.peak[0] * scalar, BestRefinedMatch.peak[1] * scalar) - BestRefinedMatch = nornir_imageregistration.AlignmentRecord(AdjustedPeak, BestRefinedMatch.weight, - BestRefinedMatch.angle, IsFlipped) + AdjustedPeak = (best_refined_match.peak[0] * scalar, best_refined_match.peak[1] * scalar) + best_refined_match = nornir_imageregistration.AlignmentRecord(AdjustedPeak, best_refined_match.weight, + best_refined_match.angle, is_flipped) - if WarpedImageScalingRequired: - # AdjustedPeak = BestRefinedMatch.peak * (1.0 / WarpedImageScaleFactors) - BestRefinedMatch = nornir_imageregistration.AlignmentRecord(BestRefinedMatch.peak, BestRefinedMatch.weight, - BestRefinedMatch.angle, IsFlipped, - WarpedImageScaleFactors) + if settings.source_image_scaling_required: + # AdjustedPeak = best_refined_match.peak * (1.0 / WarpedImageScaleFactors) + best_refined_match = nornir_imageregistration.AlignmentRecord(best_refined_match.peak, + best_refined_match.weight, + best_refined_match.angle, is_flipped, + settings.source_image_scale_factors) - # BestRefinedMatch.CorrectPeakForOriginalImageSize(imFixed.shape, imWarped.shape) + # best_refined_match.CorrectPeakForOriginalImageSize(imFixed.shape, source_image.shape) - return BestRefinedMatch + return best_refined_match -def ScoreOneAngle(imFixed_original: NDArray, imWarped_original: NDArray, - FixedImageShape: tuple[int, int], WarpedImageShape: tuple[int, int], +def ScoreOneAngle(target_original: NDArray, source_original: NDArray, + target_image_shape: tuple[int, int], source_image_shape: tuple[int, int], angle: float, - fixedStats: nornir_imageregistration.ImageStats | None = None, - warpedStats: nornir_imageregistration.ImageStats | None = None, - FixedImagePrePadded: bool = True, MinOverlap: float = 0.75): + target_stats: nornir_imageregistration.ImageStats | None = None, + source_stats: nornir_imageregistration.ImageStats | None = None, + target_image_prepadded: bool = True, min_overlap: float = 0.75): '''Returns an alignment score for a fixed image and an image rotated at a specified angle''' - imFixed = nornir_imageregistration.ImageParamToImageArray(imFixed_original, - dtype=nornir_imageregistration.default_image_dtype()) - imWarped = nornir_imageregistration.ImageParamToImageArray(imWarped_original, - dtype=nornir_imageregistration.default_image_dtype()) - - use_cp = nornir_imageregistration.GetActiveComputationLib() == nornir_imageregistration.ComputationLib.cupy - # Use of cupy or numpy - xp = cp.get_array_module(imFixed_original) - # Use of cupyx.scipy.fft or scipy.fft - xp_scipy = cupyx.scipy.get_array_module(imFixed_original) - rotate = xp_scipy.ndimage.rotate - - # imFixed = cp.asarray(imFixed) if use_cp and not isinstance(imFixed, cp.ndarray) else imFixed - # imWarped = cp.asarray(imWarped) if use_cp and not isinstance(imWarped, cp.ndarray) else imWarped - - # gc.set_debug(gc.DEBUG_LEAK) - if fixedStats is None: - fixedStats = nornir_imageregistration.ImageStats.CalcStats(imFixed) - - if warpedStats is None: - warpedStats = nornir_imageregistration.ImageStats.CalcStats(imWarped) - - OKToDelimWarped = False - if angle != 0: - # This confused me for years, but the implementation of rotate calls affine_transform with - # the rotation matrix. However the docs for affine_transform state it needs to be called - # with the inverse transform. Hence negating the angle here. - try: - if use_cp: - imWarped = rotate(imWarped, axes=(0, 1), angle=-angle, cval=np.nan) - else: - imWarped = rotate(imWarped.astype(np.float32, copy=False), axes=(0, 1), angle=-angle, - cval=np.nan).astype( - imWarped.dtype, copy=False) # Numpy cannot rotate float16 images - except RuntimeWarning as e: - pass - imWarpedEmptyIndicies = xp.isnan(imWarped) - imWarped[imWarpedEmptyIndicies] = warpedStats.GenerateNoise(xp.sum(imWarpedEmptyIndicies), dtype=imWarped.dtype) - OKToDelimWarped = True - - RotatedWarped = nornir_imageregistration.PadImageForPhaseCorrelation(imWarped, ImageMedian=warpedStats.median, - ImageStdDev=warpedStats.std, - MinOverlap=MinOverlap) - - assert (RotatedWarped.shape[0] > 0) - assert (RotatedWarped.shape[1] > 0) - - if not FixedImagePrePadded: - PaddedFixed = nornir_imageregistration.PadImageForPhaseCorrelation(imFixed, ImageMedian=fixedStats.median, - ImageStdDev=fixedStats.std, - MinOverlap=MinOverlap) - else: - PaddedFixed = imFixed - - # print str(PaddedFixed.shape) + ' ' + str(RotatedPaddedWarped.shape) - - TargetHeight = max([PaddedFixed.shape[0], RotatedWarped.shape[0]]) - TargetWidth = max([PaddedFixed.shape[1], RotatedWarped.shape[1]]) - - # Why is MinOverlap hard-coded to 1.0? - # PadImageForPhaseCorrelation will always return a copy, so don't call it unless we need to - if not np.array_equal(imFixed.shape, np.array((TargetHeight, TargetWidth))): - PaddedFixed = nornir_imageregistration.PadImageForPhaseCorrelation(imFixed, NewWidth=TargetWidth, - NewHeight=TargetHeight, - ImageMedian=fixedStats.median, - ImageStdDev=fixedStats.std, MinOverlap=1.0) + try: + im_target = nornir_imageregistration.ImageParamToImageArray(target_original, + dtype=nornir_imageregistration.default_image_dtype()) + im_source = nornir_imageregistration.ImageParamToImageArray(source_original, + dtype=nornir_imageregistration.default_image_dtype()) + + use_cp = nornir_imageregistration.GetActiveComputationLib() == nornir_imageregistration.ComputationLib.cupy + # Use of cupy or numpy + xp = cp.get_array_module(target_original) + # Use of cupyx.scipy.fft or scipy.fft + xp_scipy = cupyx.scipy.get_array_module(target_original) + rotate = xp_scipy.ndimage.rotate + + # im_target = cp.asarray(im_target) if use_cp and not isinstance(im_target, cp.ndarray) else im_target + # im_source = cp.asarray(im_source) if use_cp and not isinstance(im_source, cp.ndarray) else im_source + + # gc.set_debug(gc.DEBUG_LEAK) + if target_stats is None: + target_stats = nornir_imageregistration.ImageStats.CalcStats(im_target) + + if source_stats is None: + source_stats = nornir_imageregistration.ImageStats.CalcStats(im_source) + + OKToDelimWarped = False + if angle != 0: + # This confused me for years, but the implementation of rotate calls affine_transform with + # the rotation matrix. However the docs for affine_transform state it needs to be called + # with the inverse transform. Hence negating the angle here. + try: + if use_cp: + im_source = rotate(im_source, axes=(0, 1), angle=-angle, cval=np.nan) + else: + im_source = rotate(im_source.astype(np.float32, copy=False), axes=(0, 1), angle=-angle, + cval=np.nan).astype( + im_source.dtype, copy=False) # Numpy cannot rotate float16 images + except RuntimeWarning as e: + pass + im_source_empty_entries = xp.isnan(im_source) + im_source[im_source_empty_entries] = source_stats.GenerateNoise(xp.sum(im_source_empty_entries), + dtype=im_source.dtype) + OKToDelimWarped = True + + rotated_source = nornir_imageregistration.PadImageForPhaseCorrelation(im_source, + ImageMedian=source_stats.median, + ImageStdDev=source_stats.std, + MinOverlap=min_overlap) + + assert (rotated_source.shape[0] > 0) + assert (rotated_source.shape[1] > 0) + + if not target_image_prepadded: + padded_target = nornir_imageregistration.PadImageForPhaseCorrelation(im_target, + ImageMedian=target_stats.median, + ImageStdDev=target_stats.std, + MinOverlap=min_overlap) + else: + padded_target = im_target + + # print str(padded_target.shape) + ' ' + str(rotated_padded_source.shape) + + TargetHeight = max([padded_target.shape[0], rotated_source.shape[0]]) + TargetWidth = max([padded_target.shape[1], rotated_source.shape[1]]) + + # Why is MinOverlap hard-coded to 1.0? + # PadImageForPhaseCorrelation will always return a copy, so don't call it unless we need to + if not np.array_equal(im_target.shape, np.array((TargetHeight, TargetWidth))): + padded_target = nornir_imageregistration.PadImageForPhaseCorrelation(im_target, NewWidth=TargetWidth, + NewHeight=TargetHeight, + ImageMedian=target_stats.median, + ImageStdDev=target_stats.std, + MinOverlap=1.0) + print(f"{angle}: Padding target image to {padded_target.shape}") + else: + print(f"{angle}: No additional padding {padded_target.shape}") - if np.array_equal(RotatedWarped.shape, np.array((TargetHeight, TargetWidth))): - RotatedPaddedWarped = RotatedWarped - else: - RotatedPaddedWarped = nornir_imageregistration.PadImageForPhaseCorrelation(RotatedWarped, NewWidth=TargetWidth, - NewHeight=TargetHeight, - ImageMedian=warpedStats.median, - ImageStdDev=warpedStats.std, - MinOverlap=1.0) + if np.array_equal(rotated_source.shape, np.array((TargetHeight, TargetWidth))): + rotated_padded_source = rotated_source + else: + rotated_padded_source = nornir_imageregistration.PadImageForPhaseCorrelation(rotated_source, + NewWidth=TargetWidth, + NewHeight=TargetHeight, + ImageMedian=source_stats.median, + ImageStdDev=source_stats.std, + MinOverlap=1.0) - assert (np.array_equal(PaddedFixed.shape, RotatedPaddedWarped.shape)) + assert (np.array_equal(padded_target.shape, rotated_padded_source.shape)) - # if OKToDelimWarped: - del imWarped - del imFixed + # if OKToDelimWarped: + del im_source + del im_target - del RotatedWarped + del rotated_source - # if use_cp and not isinstance(PaddedFixed, cp.ndarray): - # PaddedFixed = cp.asarray(PaddedFixed) - # - # if use_cp and not isinstance(RotatedPaddedWarped, cp.ndarray): - # RotatedPaddedWarped = cp.asarray(RotatedPaddedWarped) + # if use_cp and not isinstance(padded_target, cp.ndarray): + # padded_target = cp.asarray(padded_target) + # + # if use_cp and not isinstance(rotated_padded_source, cp.ndarray): + # rotated_padded_source = cp.asarray(rotated_padded_source) - CorrelationImage = nornir_imageregistration.ImagePhaseCorrelation(PaddedFixed, RotatedPaddedWarped, fixedStats.mean, - warpedStats.mean) + correlation_image = nornir_imageregistration.ImagePhaseCorrelation(padded_target, rotated_padded_source, + target_stats.mean, + source_stats.mean) - del PaddedFixed - del RotatedPaddedWarped + del padded_target + del rotated_padded_source - CorrelationImage = xp_scipy.fft.fftshift(CorrelationImage) - try: - CorrelationImage -= CorrelationImage.min() - CorrelationImage /= CorrelationImage.max() - except FloatingPointError as e: - print(f"Floating point error: {e} for {CorrelationImage.min()} or {CorrelationImage.max()}") - record = nornir_imageregistration.AlignmentRecord((0, 0), 0, 0) + correlation_image = xp_scipy.fft.fftshift(correlation_image) + try: + correlation_image -= correlation_image.min() + correlation_image /= correlation_image.max() + except FloatingPointError as e: + print(f"Floating point error: {e} for {correlation_image.min()} or {correlation_image.max()}") + record = nornir_imageregistration.AlignmentRecord((0, 0), 0, 0) + return record + + # Timer.Start('Find Peak') + + # Note - Clement: overlap_mask still uses numpy (and not cupy) + overlap_mask = nornir_imageregistration.overlapmasking.GetOverlapMask(target_image_shape, source_image_shape, + correlation_image.shape, min_overlap, + MaxOverlap=1.0) + if use_cp and not isinstance(overlap_mask, cp.ndarray): + overlap_mask = cp.asarray(overlap_mask) + + (peak, weight) = nornir_imageregistration.FindPeak(correlation_image, overlap_mask) + del overlap_mask + del correlation_image + + record = nornir_imageregistration.AlignmentRecord(peak, weight, angle) return record - - # Timer.Start('Find Peak') - - # Note - Clement: OverlapMask still uses numpy (and not cupy) - OverlapMask = nornir_imageregistration.overlapmasking.GetOverlapMask(FixedImageShape, WarpedImageShape, - CorrelationImage.shape, MinOverlap, - MaxOverlap=1.0) - if use_cp and not isinstance(OverlapMask, cp.ndarray): - OverlapMask = cp.asarray(OverlapMask) - - (peak, weight) = nornir_imageregistration.FindPeak(CorrelationImage, OverlapMask) - del OverlapMask - del CorrelationImage - - nornir_imageregistration.close_shared_memory(imFixed_original) - nornir_imageregistration.close_shared_memory(imWarped_original) - - record = nornir_imageregistration.AlignmentRecord(peak, weight, angle) - return record + finally: + nornir_imageregistration.close_shared_memory(target_original) + nornir_imageregistration.close_shared_memory(source_original) def GetFixedAndWarpedImageStats(imFixed, imWarped): @@ -330,142 +358,155 @@ def GetFixedAndWarpedImageStats(imFixed, imWarped): return fixedStats, warpedStats -def _find_best_angle(imFixed: NDArray[np.floating], - imWarped: NDArray[np.floating], - fixed_stats: nornir_imageregistration.ImageStats, - warped_stats: nornir_imageregistration.ImageStats, - AngleList: set[float] | None, - MinOverlap: float = 0.75, +def _find_best_angle(source_image: NDArray[np.floating], + target_image: NDArray[np.floating], + source_stats: nornir_imageregistration.ImageStats, + target_stats: nornir_imageregistration.ImageStats, + angle_range: NDArray[float], + min_overlap: float = 0.75, SingleThread: bool = False, use_cluster: bool = False): '''Find the best angle to align two images. This function can be very memory intensive. Setting SingleThread=True makes debugging easier''' - Debug = False - pool = None - use_cp = nornir_imageregistration.GetActiveComputationLib() == nornir_imageregistration.ComputationLib.cupy + try: + Debug = False + pool = None + use_cp = nornir_imageregistration.GetActiveComputationLib() == nornir_imageregistration.ComputationLib.cupy - # Temporarily disable until we have cluster pool working again. Leaving this on eliminates shared memory which is a big optimization - use_cluster = False + # Temporarily disable until we have cluster pool working again. Leaving this on eliminates shared memory which is a big optimization + use_cluster = False - if len(AngleList) <= 1: - SingleThread = True + if len(angle_range) <= 1: + SingleThread = True - if not SingleThread: - if Debug: - pool = nornir_pools.GetThreadPool(Poolname=None, num_threads=3) - elif use_cluster: - pool = nornir_pools.GetGlobalClusterPool() - else: - pool = nornir_pools.GetGlobalMultithreadingPool() + if not SingleThread: + if Debug: + pool = nornir_pools.GetThreadPool(Poolname=None, num_threads=3) + elif use_cluster: + pool = nornir_pools.GetGlobalClusterPool() + else: + pool = nornir_pools.GetGlobalMultithreadingPool() - AngleMatchValues = list() - taskList = list() + # Preallocate lists to store results of each angle + AngleMatchValues = list() # type: list[AlignmentRecord | None] + taskList = list() # type: list[nornir_pools.Task | None] - # MaxRotatedDimension = max([max(imFixed), max(imWarped)]) * 1.4143 - # MinRotatedDimension = max(min(imFixed), min(imWarped)) - # - # SmallPaddedFixed = PadImageForPhaseCorrelation(imFixed, MaxOffset=0.1) - # LargePaddedFixed = PadImageForPhaseCorrelation(imFixed, MaxOffset=0.1) + # MaxRotatedDimension = max([max(imFixed), max(imWarped)]) * 1.4143 + # MinRotatedDimension = max(min(imFixed), min(imWarped)) + # + # SmallPaddedFixed = PadImageForPhaseCorrelation(imFixed, MaxOffset=0.1) + # LargePaddedFixed = PadImageForPhaseCorrelation(imFixed, MaxOffset=0.1) - PaddedFixed = nornir_imageregistration.PadImageForPhaseCorrelation(imFixed, - MinOverlap=MinOverlap, - ImageMedian=fixed_stats.median, - ImageStdDev=fixed_stats.std) + padded_target = nornir_imageregistration.PadImageForPhaseCorrelation(target_image, + MinOverlap=min_overlap, + ImageMedian=target_stats.median, + ImageStdDev=target_stats.std) - # Create a shared read-only memory map for the Padded fixed image + # Create a shared read-only memory map for the Padded fixed image - if not (use_cluster or SingleThread): - # temp_padded_fixed_memmap = nornir_imageregistration.CreateTemporaryReadonlyMemmapFile(PaddedFixed) - # temp_shared_warp_memmap = nornir_imageregistration.CreateTemporaryReadonlyMemmapFile(imWarped) + if not (use_cluster or SingleThread): + # temp_padded_fixed_memmap = nornir_imageregistration.CreateTemporaryReadonlyMemmapFile(padded_target) + # temp_shared_warp_memmap = nornir_imageregistration.CreateTemporaryReadonlyMemmapFile(imWarped) - # temp_padded_fixed_memmap.mode = 'r' # We do not want functions we pass the memmap modifying the original data - # temp_shared_warp_memmap.mode = 'r' # We do not want functions we pass the memmap modifying the original data + # temp_padded_fixed_memmap.mode = 'r' # We do not want functions we pass the memmap modifying the original data + # temp_shared_warp_memmap.mode = 'r' # We do not want functions we pass the memmap modifying the original data - shared_fixed_metadata, SharedPaddedFixed = nornir_imageregistration.npArrayToSharedArray(PaddedFixed) - shared_warped_metadata, SharedWarped = nornir_imageregistration.npArrayToSharedArray(imWarped) - # SharedPaddedFixed = np.save(PaddedFixed, ) - else: - SharedPaddedFixed = PaddedFixed.astype(nornir_imageregistration.default_image_dtype(), - copy=False) if not use_cp else cp.array(PaddedFixed, - nornir_imageregistration.default_image_dtype()) - SharedWarped = imWarped.astype(nornir_imageregistration.default_image_dtype(), - copy=False) if not use_cp else cp.array(imWarped, - nornir_imageregistration.default_image_dtype()) - - CheckTaskInterval = 16 - - fixed_shape = imFixed.shape - warped_shape = imWarped.shape - max_task_count = multiprocessing.cpu_count() * 1.5 - - for i, theta in enumerate(AngleList): - if SingleThread: - record = ScoreOneAngle(SharedPaddedFixed, SharedWarped, fixed_shape, - warped_shape, theta, fixedStats=fixed_stats, - warpedStats=warped_stats, - MinOverlap=MinOverlap) - AngleMatchValues.append(record) - elif use_cluster: - task = pool.add_task(str(theta), ScoreOneAngle, SharedPaddedFixed, SharedWarped, - fixed_shape, warped_shape, theta, fixedStats=fixed_stats, warpedStats=warped_stats, - MinOverlap=MinOverlap) - taskList.append(task) + shared_target_metadata, shared_padded_target = nornir_imageregistration.npArrayToSharedArray(padded_target) + shared_source_metadata, shared_source = nornir_imageregistration.npArrayToSharedArray(source_image + ) + # shared_padded_target = np.save(padded_target, ) else: - task = pool.add_task(str(theta), ScoreOneAngle, shared_fixed_metadata, shared_warped_metadata, - fixed_shape, warped_shape, - theta, - fixedStats=fixed_stats, warpedStats=warped_stats, MinOverlap=MinOverlap) - taskList.append(task) - - if not i % CheckTaskInterval == 0: - continue - - # I don't like this, but it lets me delete tasks before filling the queue which may save some memory. - # No sense checking unless we've already filled the queue though - if len(taskList) > max_task_count: + shared_target_metadata = None + shared_source_metadata = None + shared_padded_target = padded_target.astype(nornir_imageregistration.default_image_dtype(), + copy=False) if not use_cp else cp.array(padded_target, + nornir_imageregistration.default_image_dtype()) + shared_source = source_image.astype(nornir_imageregistration.default_image_dtype(), + copy=False) if not use_cp else cp.array(source_image, + nornir_imageregistration.default_image_dtype()) + + CheckTaskInterval = 16 + + source_shape = source_image.shape + target_shape = target_image.shape + max_task_count = multiprocessing.cpu_count() * 1.5 + + for i, theta in enumerate(angle_range): + if SingleThread: + record = ScoreOneAngle(target_original=shared_padded_target, source_original=shared_source, + target_image_shape=target_shape, source_image_shape=source_shape, + angle=theta, + target_stats=target_stats, source_stats=source_stats, + min_overlap=min_overlap) + AngleMatchValues.append(record) + elif use_cluster: + task = pool.add_task(str(theta), ScoreOneAngle, + target_original=shared_padded_target, source_original=shared_source, + target_image_shape=target_shape, source_image_shape=source_shape, + angle=theta, + target_stats=target_stats, source_stats=source_stats, + min_overlap=min_overlap) + taskList.append(task) + else: + task = pool.add_task(str(theta), ScoreOneAngle, + target_original=shared_target_metadata, source_original=shared_source_metadata, + target_image_shape=target_shape, source_image_shape=source_shape, + angle=theta, + target_stats=target_stats, source_stats=source_stats, + min_overlap=min_overlap) + taskList.append(task) + + if not i % CheckTaskInterval == 0: + continue + + # I don't like this, but it lets me delete tasks before filling the queue which may save some memory. + # No sense checking unless we've already filled the queue though + if len(taskList) > max_task_count: + for iTask in range(len(taskList) - 1, -1, -1): + if taskList[iTask].iscompleted: + record = taskList[iTask].wait_return() + AngleMatchValues.append(record) + del taskList[iTask] + + # TestOneAngle(shared_padded_target, shared_source, angle, None, MinOverlap) + + # taskList.sort(key=tpool.Task.name) + + while len(taskList) > 0: for iTask in range(len(taskList) - 1, -1, -1): if taskList[iTask].iscompleted: record = taskList[iTask].wait_return() AngleMatchValues.append(record) del taskList[iTask] - # TestOneAngle(SharedPaddedFixed, SharedWarped, angle, None, MinOverlap) - - # taskList.sort(key=tpool.Task.name) - - while len(taskList) > 0: - for iTask in range(len(taskList) - 1, -1, -1): - if taskList[iTask].iscompleted: - record = taskList[iTask].wait_return() - AngleMatchValues.append(record) - del taskList[iTask] - - if len(taskList) > 0: - # Wait a bit before checking the task list - sleep(0.5) + if len(taskList) > 0: + # Wait a bit before checking the task list + sleep(0.5) - # print(str(record.angle) + ' = ' + str(record.peak) + ' weight: ' + str(record.weight) + '\n') + # print(str(record.angle) + ' = ' + str(record.peak) + ' weight: ' + str(record.weight) + '\n') - # ShowGrayscale(NormCorrelationImage) + # ShowGrayscale(NormCorrelationImage) - # print(str(AngleMatchValues)) + # print(str(AngleMatchValues)) - # Delete the pool to ensure extra python threads do not stick around - # if pool is not None: - # pool.shutdown() + # Delete the pool to ensure extra python threads do not stick around + # if pool is not None: + # pool.shutdown() - del PaddedFixed + del padded_target - BestMatch = max(AngleMatchValues, key=nornir_imageregistration.AlignmentRecord.WeightKey) + BestMatch = max(AngleMatchValues, key=nornir_imageregistration.AlignmentRecord.WeightKey) + return BestMatch + finally: - if not (use_cluster or SingleThread): - nornir_imageregistration.unlink_shared_memory(shared_fixed_metadata) - nornir_imageregistration.unlink_shared_memory(shared_warped_metadata) - # os.remove(temp_shared_warp_memmap.path) - # os.remove(temp_padded_fixed_memmap.path) + if shared_target_metadata is not None: + nornir_imageregistration.unlink_shared_memory(shared_target_metadata) + if shared_source_metadata is not None: + nornir_imageregistration.unlink_shared_memory(shared_source_metadata) - return BestMatch + # os.remove(temp_shared_warp_memmap.path) + # os.remove(temp_padded_fixed_memmap.path) def __ExecuteProfiler(): diff --git a/nornir_imageregistration/transforms/rigid.py b/nornir_imageregistration/transforms/rigid.py index d1a1d05..af7b589 100644 --- a/nornir_imageregistration/transforms/rigid.py +++ b/nornir_imageregistration/transforms/rigid.py @@ -342,7 +342,7 @@ def RotateSourcePoints(self, rangle: float, rotation_center: NDArray[np.floating xp = nornir_imageregistration.GetComputationModule() if rotation_center is not None: - self.source_space_center_of_rotation = xp.array(rotation_center) + self._source_space_center_of_rotation = xp.array(rotation_center) self._update_transform_matrix() self.OnTransformChanged() diff --git a/nornir_imageregistration/transforms/triangulation.py b/nornir_imageregistration/transforms/triangulation.py index 7a509d7..824a90f 100644 --- a/nornir_imageregistration/transforms/triangulation.py +++ b/nornir_imageregistration/transforms/triangulation.py @@ -3,7 +3,7 @@ @author: Jamesan ''' - +from collections.abc import Iterable import logging from multiprocessing.managers import Value @@ -244,7 +244,11 @@ def UpdateSourcePointsByPosition(self, old_points: NDArray[np.floating], new_poi return self.UpdateSourcePointsByIndex(index, new_points) def RemovePoint(self, index: int | NDArray[np.integer]): - if self._points.shape[0] - len(index) <= 3: + nToRemove = 1 + if isinstance(index, Iterable): + nToRemove = len(index) + + if self._points.shape[0] - nToRemove < 3: raise ValueError("Cannot remove points, must have at least three points") self._points = np.delete(self._points, index, 0) diff --git a/nornir_imageregistration/transforms/utils.py b/nornir_imageregistration/transforms/utils.py index 5939ab0..1180207 100644 --- a/nornir_imageregistration/transforms/utils.py +++ b/nornir_imageregistration/transforms/utils.py @@ -107,7 +107,7 @@ def ScaleMatrixXY(scale: float | Sequence[float]) -> NDArray[np.floating]: elif isinstance(scale, float): return xp.array([[scale, 0, 0], [0, scale, 0], [0, 0, 1]]) elif isinstance(scale, int): - return xp.array([[scale, 0, 0], [0, scale, 0], [0, 0, 1]], np.floating) + return xp.array([[scale, 0, 0], [0, scale, 0], [0, 0, 1]], float) elif hasattr(scale, "__iter__"): return xp.array([[scale[0], 0, 0], [0, scale[1], 0], [0, 0, 1]]) diff --git a/test/test_SliceToSliceBrute.py b/test/test_SliceToSliceBrute.py index 10f2a66..37b98b6 100644 --- a/test/test_SliceToSliceBrute.py +++ b/test/test_SliceToSliceBrute.py @@ -5,11 +5,15 @@ ''' import os import unittest +import matplotlib -#Check if cupy is available, and if it is not import thunks that refer to scipy/numpy +matplotlib.use('qtAgg') + +# Check if cupy is available, and if it is not import thunks that refer to scipy/numpy try: import cupy as cp import cupyx + init_context = cp.zeros((64, 64)) except ModuleNotFoundError: import nornir_imageregistration.cupy_thunk as cp @@ -17,9 +21,9 @@ except ImportError: import nornir_imageregistration.cupy_thunk as cp import nornir_imageregistration.cupyx_thunk as cupyx - + import nornir_imageregistration -from nornir_imageregistration import alignment_record +from nornir_imageregistration import AlignmentRecord, alignment_record import nornir_imageregistration.core as core import nornir_imageregistration.files import nornir_imageregistration.scripts.nornir_rotate_translate @@ -30,8 +34,8 @@ import setup_imagetest - -def CheckAlignmentRecord(test, arecord, angle, X, Y, flipud=False, adelta=None, sdelta=None): +def CheckAlignmentRecord(test: unittest.TestCase, arecord: alignment_record.AlignmentRecord, angle: float, X: float, + Y: float, flipud: bool = False, adelta: float | None = None, sdelta: float | None = None): '''Verifies that an alignment record is more or less equal to expected values''' angle = float(angle) @@ -53,7 +57,7 @@ def CheckAlignmentRecord(test, arecord, angle, X, Y, flipud=False, adelta=None, class TestStos(setup_imagetest.ImageTestBase): def testStosWrite(self): - InputDir = 'C:\\Buildscript\\Test\\images\\' + InputDir = 'C:\\Buildscriptd\\Test\\images\\' OutputDir = 'C:\\Temp\\' WarpedImagePath = os.path.join(self.ImportedDataPath, @@ -104,7 +108,7 @@ def testStosBrute_Cluster(self): def testStosBrute_GPU(self): if not nornir_imageregistration.HasCupy(): - return + return nornir_imageregistration.SetActiveComputationLib(nornir_imageregistration.ComputationLib.cupy) self.RunBasicBruteAlignment(self.FixedImagePath, self.WarpedImagePath, SingleThread=True, FlipUD=False) @@ -123,8 +127,8 @@ def testStosBruteWithFlip_Cluster(self): def testStosBruteWithFlip_GPU(self): if not nornir_imageregistration.HasCupy(): - return - + return + nornir_imageregistration.SetActiveComputationLib(nornir_imageregistration.ComputationLib.cupy) self.RunBasicBruteAlignment(self.FixedImagePath, self.WarpedImagePathFlipped, SingleThread=True, FlipUD=True) @@ -140,7 +144,7 @@ def RunBasicBruteAlignment(self, FixedImagePath: str, SingleThread = True if nornir_imageregistration.GetActiveComputationLib() == nornir_imageregistration.ComputationLib.cupy else SingleThread - MinOverlap = 0.75 + MinOverlap = 0.5 timer = TaskTimer() @@ -211,12 +215,12 @@ def testStosBruteWithMask_MultiThread(self): savedstosObj = AlignmentRecord.ToStos(self.FixedImagePath, self.WarpedImagePath, self.FixedImageMaskPath, self.WarpedImageMaskPath, PixelSpacing=1) - self.CheckStosObj(savedstosObj,'17-18_brute_WithMask.stos', self.FixedImageMaskPath, self.WarpedImageMaskPath) + self.CheckStosObj(savedstosObj, '17-18_brute_WithMask.stos', self.FixedImageMaskPath, self.WarpedImageMaskPath) def testStosBruteWithMask_GPU(self): if not nornir_imageregistration.HasCupy(): - return - + return + nornir_imageregistration.SetActiveComputationLib(nornir_imageregistration.ComputationLib.cupy) AlignmentRecord = self.RunBasicBruteAlignmentWithMask(self.FixedImagePath, self.WarpedImagePath, self.FixedImageMaskPath, self.WarpedImageMaskPath, @@ -225,7 +229,8 @@ def testStosBruteWithMask_GPU(self): savedstosObj = AlignmentRecord.ToStos(self.FixedImagePath, self.WarpedImagePath, self.FixedImageMaskPath, self.WarpedImageMaskPath, PixelSpacing=1) - self.CheckStosObj(savedstosObj, '17-18_brute_WithMask_GPU.stos', self.FixedImageMaskPath, self.WarpedImageMaskPath) + self.CheckStosObj(savedstosObj, '17-18_brute_WithMask_GPU.stos', self.FixedImageMaskPath, + self.WarpedImageMaskPath) def RunBasicBruteAlignmentWithMask(self, FixedImagePath: str, @@ -250,6 +255,7 @@ def RunBasicBruteAlignmentWithMask(self, WarpedImagePath, FixedImageMaskPath, WarpedImageMaskPath, + LargestDimension=1024, AngleSearchRange=AngleSearchRange, WarpedImageScaleFactors=WarpedImageScaleFactors, SingleThread=SingleThread, @@ -293,7 +299,7 @@ def testStosBruteScaleMismatchWithMask(self): def testStosBruteScaleMismatchWithMask_GPU(self): if not nornir_imageregistration.HasCupy(): - return + return nornir_imageregistration.SetActiveComputationLib(nornir_imageregistration.ComputationLib.cupy) self.runStosBruteScaleMismatchWithMask() @@ -361,14 +367,12 @@ def testStosBruteExecuteWithMask(self): def testStosBruteExecuteWithMask_GPU(self): if not nornir_imageregistration.HasCupy(): - return - + return + nornir_imageregistration.SetActiveComputationLib(nornir_imageregistration.ComputationLib.cupy) self.runStosBruteExecuteWithMask() - - class TestStosBruteToSameImage(setup_imagetest.ImageTestBase): def setUp(self): @@ -431,8 +435,8 @@ def testSameTEMImageFast_MultiThread(self): def testSameTEMImageFast_GPU(self): '''Make sure the same image aligns to itself with peak (0,0) and angle 0''' if not nornir_imageregistration.HasCupy(): - return - + return + nornir_imageregistration.SetActiveComputationLib(nornir_imageregistration.ComputationLib.cupy) self.assertTrue(os.path.exists(self.FixedImagePath), "Missing test input") self.assertTrue(os.path.exists(self.FixedImageMaskPath), "Missing test input") @@ -471,8 +475,8 @@ def testSameTEMImage_MultiThread(self): def testSameTEMImage_GPU(self): '''Make sure the same image aligns to itself with peak (0,0) and angle 0''' if not nornir_imageregistration.HasCupy(): - return - + return + nornir_imageregistration.SetActiveComputationLib(nornir_imageregistration.ComputationLib.cupy) self.assertTrue(os.path.exists(self.FixedImagePath), "Missing test input") self.assertTrue(os.path.exists(self.FixedImageMaskPath), "Missing test input")