From 534c848b1dd644ed0823fe819bf299eb64b86d23 Mon Sep 17 00:00:00 2001 From: ashmeigh Date: Fri, 24 May 2024 10:16:32 +0100 Subject: [PATCH 1/2] updated branch --- .../roi_normalisation/roi_normalisation.py | 144 ++++++++---------- 1 file changed, 65 insertions(+), 79 deletions(-) diff --git a/mantidimaging/core/operations/roi_normalisation/roi_normalisation.py b/mantidimaging/core/operations/roi_normalisation/roi_normalisation.py index a03f2d20334..9f2793c0628 100644 --- a/mantidimaging/core/operations/roi_normalisation/roi_normalisation.py +++ b/mantidimaging/core/operations/roi_normalisation/roi_normalisation.py @@ -3,16 +3,12 @@ from __future__ import annotations from functools import partial -from logging import getLogger from typing import TYPE_CHECKING import numpy as np - from mantidimaging import helper as h from mantidimaging.core.operations.base_filter import BaseFilter, FilterGroup from mantidimaging.core.parallel import shared as ps -from mantidimaging.core.parallel import utility as pu -from mantidimaging.core.utility.progress_reporting import Progress from mantidimaging.core.utility.sensible_roi import SensibleROI from mantidimaging.gui.utility import add_property_to_form from mantidimaging.gui.utility.qt_helpers import Type @@ -33,11 +29,8 @@ class RoiNormalisationFilter(BaseFilter): """Normalises the image data by using the average value in a no-sample (air) region of interest (ROI) of the image. This scaling operation allows to account for beam fluctuations and different exposure times of projections. - Intended to be used on: Projections - When: Always, to ensure that any fluctuations in beam intensity are normalised. - Note: Beam fluctuations are visible by horizontal lines in the sinograms, as some projections are brighter/darker than the neighbours. This can be fixed with this operation. """ @@ -51,44 +44,94 @@ def filter_func(images: ImageStack, flat_field: ImageStack | None = None, progress=None): """Normalise by beam intensity. - This does NOT do any checks if the Air Region is out of bounds! If the Air Region is out of bounds, the crop will fail at runtime. If the Air Region is in bounds, but has overlapping coordinates the crop give back a 0 shape of the coordinates that were wrong. - :param images: Sample data which is to be processed. Expected in radiograms - :param region_of_interest: The order is - Left Top Right Bottom. The air region for which grey values are summed up and used for normalisation/scaling. - :param normalisation_mode: Controls what the ROI counts are normalised to. 'Stack Average' : The mean value of the air region across all projections is preserved. 'Flat Field' : The mean value of the air regions in the projections is made equal to the mean value of the air region in the flat field image. - :param flat_field: Flat field to use if 'Flat Field' mode is enabled. - :param progress: Reference to a progress bar object - :returns: Filtered data (stack of images) """ - if normalisation_mode not in modes(): - raise ValueError(f"Unknown normalisation_mode: {normalisation_mode}, should be one of {modes()}") + if not region_of_interest: + raise ValueError('region_of_interest must be provided') - if normalisation_mode == "Flat Field" and flat_field is None: - raise ValueError('flat_field must provided if using normalisation_mode of "Flat Field"') + if isinstance(region_of_interest, list): + region_of_interest = SensibleROI.from_list(region_of_interest) h.check_data_stack(images) - if not region_of_interest: - raise ValueError('region_of_interest must be provided') + global_params = RoiNormalisationFilter.calculate_global(images, region_of_interest, normalisation_mode, + flat_field) + + params = { + 'region_of_interest': region_of_interest, + 'normalisation_mode': normalisation_mode, + 'global_params': global_params + } + + print("Filter Function Params:", params) + + ps.run_compute_func(RoiNormalisationFilter.compute_function, images.data.shape[0], images.shared_array, params) - progress = Progress.ensure_instance(progress, task_name='ROI Normalisation') - _execute(images, region_of_interest, normalisation_mode, flat_field, progress) h.check_data_stack(images) + return images + @staticmethod + def calculate_global(images, region_of_interest, normalisation_mode, flat_field): + global_params = {} + if normalisation_mode == 'Stack Average': + air_means = np.array([ + RoiNormalisationFilter._calc_mean(images.data[i], region_of_interest.left, region_of_interest.top, + region_of_interest.right, region_of_interest.bottom) + for i in range(images.data.shape[0]) + ]) + global_params['global_mean'] = np.mean(air_means) + elif normalisation_mode == 'Flat Field' and flat_field is not None: + flat_field_mean = RoiNormalisationFilter._calc_mean(flat_field.data, region_of_interest.left, + region_of_interest.top, region_of_interest.right, + region_of_interest.bottom) + global_params['flat_field_mean'] = flat_field_mean + else: + raise ValueError(f"Unknown normalisation_mode: {normalisation_mode}") + + print("Global Params:", global_params) + + return global_params + + @staticmethod + def compute_function(i: int, array: np.ndarray, params): + region_of_interest = params['region_of_interest'] + normalisation_mode = params['normalisation_mode'] + global_params = params['global_params'] + + air_mean = RoiNormalisationFilter._calc_mean(array[i], region_of_interest.left, region_of_interest.top, + region_of_interest.right, region_of_interest.bottom) + + print(f"Image {i} Air Mean: {air_mean}") + + if normalisation_mode == 'Stack Average': + normalization_factor = global_params['global_mean'] / air_mean + print(f"Image {i} Stack Average Normalization Factor: {normalization_factor}") + array[i] *= normalization_factor + elif normalisation_mode == 'Flat Field': + normalization_factor = global_params['flat_field_mean'] / air_mean + print(f"Image {i} Flat Field Normalization Factor: {normalization_factor}") + array[i] *= normalization_factor + else: + raise ValueError(f"Unknown normalisation_mode: {normalisation_mode}") + + @staticmethod + def _calc_mean(data, air_left=None, air_top=None, air_right=None, air_bottom=None): + return data[air_top:air_bottom, air_left:air_right].mean() + @staticmethod def register_gui(form, on_change, view): label, roi_field = add_property_to_form("Air Region", @@ -143,63 +186,6 @@ def group_name() -> FilterGroup: return FilterGroup.Basic -def _calc_mean(data, air_left=None, air_top=None, air_right=None, air_bottom=None): - return data[air_top:air_bottom, air_left:air_right].mean() - - -def _divide_by_air(data=None, air_sums=None): - data[:] = np.true_divide(data, air_sums) - - -def _execute(images: ImageStack, - air_region: SensibleROI, - normalisation_mode: str, - flat_field: ImageStack | None, - progress=None): - log = getLogger(__name__) - - with progress: - progress.update(msg="Normalization by air region") - if isinstance(air_region, list): - air_region = SensibleROI.from_list(air_region) - - # initialise same number of air sums - img_num = images.data.shape[0] - air_means = pu.create_array((img_num, ), images.dtype) - - do_calculate_air_means = ps.create_partial(_calc_mean, - ps.return_to_second_at_i, - air_left=air_region.left, - air_top=air_region.top, - air_right=air_region.right, - air_bottom=air_region.bottom) - - arrays = [images.shared_array, air_means] - ps.execute(do_calculate_air_means, arrays, images.data.shape[0], progress) - - if normalisation_mode == 'Stack Average': - air_means.array /= air_means.array.mean() - - elif normalisation_mode == 'Flat Field' and flat_field is not None: - flat_mean = pu.create_array((flat_field.data.shape[0], ), flat_field.dtype) - arrays = [flat_field.shared_array, flat_mean] - ps.execute(do_calculate_air_means, arrays, flat_field.data.shape[0], progress) - air_means.array /= flat_mean.array.mean() - - if np.isnan(air_means.array).any(): - raise ValueError("Air region contains invalid (NaN) pixels") - - do_divide = ps.create_partial(_divide_by_air, fwd_function=ps.inplace2) - arrays = [images.shared_array, air_means] - ps.execute(do_divide, arrays, images.data.shape[0], progress) - - avg = np.average(air_means.array) - max_avg = np.max(air_means.array) / avg - min_avg = np.min(air_means.array) / avg - - log.info(f"Normalization by air region. Average: {avg}, max ratio: {max_avg}, min ratio: {min_avg}.") - - def enable_correct_fields_only(text, flat_file_widget): if text == "Flat Field": flat_file_widget.setEnabled(True) From dd2ca2549b2d4c801becdbfd75c5a3f466d4ce70 Mon Sep 17 00:00:00 2001 From: ashmeigh Date: Tue, 28 May 2024 20:31:08 +0100 Subject: [PATCH 2/2] Modification of calculate_global,compute_function Added Validation Step. Fixed issue with normalization factor --- .../roi_normalisation/roi_normalisation.py | 44 +++++++------------ 1 file changed, 15 insertions(+), 29 deletions(-) diff --git a/mantidimaging/core/operations/roi_normalisation/roi_normalisation.py b/mantidimaging/core/operations/roi_normalisation/roi_normalisation.py index 9f2793c0628..0e834c57e0a 100644 --- a/mantidimaging/core/operations/roi_normalisation/roi_normalisation.py +++ b/mantidimaging/core/operations/roi_normalisation/roi_normalisation.py @@ -87,18 +87,19 @@ def filter_func(images: ImageStack, @staticmethod def calculate_global(images, region_of_interest, normalisation_mode, flat_field): global_params = {} + air_means = np.array([ + RoiNormalisationFilter._calc_mean(images.data[i], region_of_interest.left, region_of_interest.top, + region_of_interest.right, region_of_interest.bottom) + for i in range(images.data.shape[0]) + ]) if normalisation_mode == 'Stack Average': - air_means = np.array([ - RoiNormalisationFilter._calc_mean(images.data[i], region_of_interest.left, region_of_interest.top, - region_of_interest.right, region_of_interest.bottom) - for i in range(images.data.shape[0]) - ]) - global_params['global_mean'] = np.mean(air_means) + normed_air_means = air_means / air_means.mean() + global_params['normalisation_factors'] = normed_air_means elif normalisation_mode == 'Flat Field' and flat_field is not None: - flat_field_mean = RoiNormalisationFilter._calc_mean(flat_field.data, region_of_interest.left, - region_of_interest.top, region_of_interest.right, - region_of_interest.bottom) - global_params['flat_field_mean'] = flat_field_mean + flat_means = flat_field.data[:, region_of_interest.top:region_of_interest.bottom, + region_of_interest.left:region_of_interest.right].mean(axis=(1, 2)) + normed_air_means = air_means / flat_means.mean() + global_params['normalisation_factors'] = normed_air_means else: raise ValueError(f"Unknown normalisation_mode: {normalisation_mode}") @@ -108,25 +109,10 @@ def calculate_global(images, region_of_interest, normalisation_mode, flat_field) @staticmethod def compute_function(i: int, array: np.ndarray, params): - region_of_interest = params['region_of_interest'] - normalisation_mode = params['normalisation_mode'] - global_params = params['global_params'] - - air_mean = RoiNormalisationFilter._calc_mean(array[i], region_of_interest.left, region_of_interest.top, - region_of_interest.right, region_of_interest.bottom) - - print(f"Image {i} Air Mean: {air_mean}") - - if normalisation_mode == 'Stack Average': - normalization_factor = global_params['global_mean'] / air_mean - print(f"Image {i} Stack Average Normalization Factor: {normalization_factor}") - array[i] *= normalization_factor - elif normalisation_mode == 'Flat Field': - normalization_factor = global_params['flat_field_mean'] / air_mean - print(f"Image {i} Flat Field Normalization Factor: {normalization_factor}") - array[i] *= normalization_factor - else: - raise ValueError(f"Unknown normalisation_mode: {normalisation_mode}") + normalisation_factors = params['global_params']['normalisation_factors'] + normalization_factor = normalisation_factors[i] + print(f"Image {i} Normalization Factor: {normalization_factor}") + array[i] /= normalization_factor @staticmethod def _calc_mean(data, air_left=None, air_top=None, air_right=None, air_bottom=None):