Merge pull request #214 from DiamondLightSource/datavalidator

Data validation checker

Author: dkazanc

httomolibgpu/cuda_kernels/remove_nan_inf.cu

@@ -0,0 +1,19 @@
+template <typename Type>
+__global__ void remove_nan_inf(Type *data, int Z, int M, int N, int *result) {
+  const long i = blockDim.x * blockIdx.x + threadIdx.x;
+  const long j = blockDim.y * blockIdx.y + threadIdx.y;
+  const long k = blockDim.z * blockIdx.z + threadIdx.z;
+
+  if (i >= N || j >= M || k >= Z)
+    return;
+
+  long long index = static_cast<long long>(i) + N * static_cast<long long>(j) + N * M * static_cast<long long>(k);
+
+  float val = float(data[index]); /*needs a cast to float for isnan isinf functions to work*/
+  Type zero = 0;
+  if (isnan(val) || isinf(val)) {
+    result[0] = 1;
+    data[index] = zero;
+  }
+
+}

httomolibgpu/misc/corr.py

@@ -18,9 +18,7 @@
 # Created By  : Tomography Team at DLS <[email protected]>
 # Created Date: 21/October/2022
 # ---------------------------------------------------------------------------
-""" Module for data correction. For more detailed information see :ref:`data_correction_module`.
-
-"""
+"""Module for data correction. For more detailed information see :ref:`data_correction_module`."""
 
 import numpy as np
 from typing import Union
@@ -38,6 +36,7 @@
 else:
     load_cuda_module = Mock()
 
+from httomolibgpu.misc.supp_func import data_checker
 
 __all__ = [
     "median_filter",
@@ -74,7 +73,6 @@ def median_filter(
         If the input array is not three dimensional.
     """
     input_type = data.dtype
-
     if input_type not in ["float32", "uint16"]:
         raise ValueError("The input data should be either float32 or uint16 data type")
 
@@ -84,6 +82,10 @@ def median_filter(
     else:
         raise ValueError("The input array must be a 3D array")
 
+    data = data_checker(
+        data, verbosity=True, method_name="median_filter_or_remove_outlier"
+    )
+
     if kernel_size not in [3, 5, 7, 9, 11, 13]:
         raise ValueError("Please select a correct kernel size: 3, 5, 7, 9, 11, 13")
 

httomolibgpu/misc/denoise.py

@@ -18,8 +18,7 @@
 # Created By  : Tomography Team at DLS <[email protected]>
 # Created Date: 18/December/2024
 # ---------------------------------------------------------------------------
-""" Module for data denoising. For more detailed information see :ref:`data_denoising_module`.
-"""
+"""Module for data denoising. For more detailed information see :ref:`data_denoising_module`."""
 
 import numpy as np
 from typing import Union, Optional
@@ -29,9 +28,10 @@
 cp = cupywrapper.cp
 cupy_run = cupywrapper.cupy_run
 
-from numpy import float32
 from unittest.mock import Mock
 
+from httomolibgpu.misc.supp_func import data_checker
+
 if cupy_run:
     from ccpi.filters.regularisersCuPy import ROF_TV, PD_TV
 else:
@@ -82,6 +82,8 @@ def total_variation_ROF(
         If the input array is not float32 data type.
     """
 
+    data = data_checker(data, verbosity=True, method_name="total_variation_ROF")
+
     return ROF_TV(
         data, regularisation_parameter, iterations, time_marching_parameter, gpu_id
     )
@@ -127,6 +129,8 @@ def total_variation_PD(
         If the input array is not float32 data type.
     """
 
+    data_checker(data, verbosity=True, method_name="total_variation_PD")
+
     methodTV = 0
     if not isotropic:
         methodTV = 1

httomolibgpu/misc/morph.py

@@ -35,6 +35,8 @@
 
 from typing import Literal
 
+from httomolibgpu.misc.supp_func import data_checker
+
 __all__ = [
     "sino_360_to_180",
     "data_resampler",
@@ -66,6 +68,8 @@ def sino_360_to_180(
     if data.ndim != 3:
         raise ValueError("only 3D data is supported")
 
+    data = data_checker(data, verbosity=True, method_name="sino_360_to_180")
+
     dx, dy, dz = data.shape
 
     overlap = int(np.round(overlap))
@@ -136,6 +140,8 @@ def data_resampler(
         data = cp.expand_dims(data, 1)
         axis = 1
 
+    data = data_checker(data, verbosity=True, method_name="data_resampler")
+
     N, M, Z = cp.shape(data)
 
     if axis == 0:

httomolibgpu/misc/rescale.py

@@ -18,9 +18,7 @@
 # Created By  : Tomography Team at DLS <[email protected]>
 # Created Date: 1 March 2024
 # ---------------------------------------------------------------------------
-""" Module for data rescaling. For more detailed information see :ref:`data_rescale_module`.
-
-"""
+"""Module for data rescaling. For more detailed information see :ref:`data_rescale_module`."""
 
 import numpy as np
 from httomolibgpu import cupywrapper
@@ -30,6 +28,8 @@
 
 from typing import Literal, Optional, Tuple, Union
 
+from httomolibgpu.misc.supp_func import data_checker
+
 __all__ = [
     "rescale_to_int",
 ]
@@ -80,6 +80,8 @@ def rescale_to_int(
     else:
         output_dtype = np.uint32
 
+    data = data_checker(data, verbosity=True, method_name="rescale_to_int")
+
     if cupy_run:
         xp = cp.get_array_module(data)
     else:
@@ -109,7 +111,6 @@ def rescale_to_int(
     if xp.__name__ == "numpy":
         if input_max == pow(2, 32):
             input_max -= 1
-        data[np.logical_not(np.isfinite(data))] = 0
         res = np.copy(data.astype(float))
         res[data.astype(float) < input_min] = int(input_min)
         res[data.astype(float) > input_max] = int(input_max)

httomolibgpu/misc/supp_func.py

@@ -0,0 +1,191 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+# ---------------------------------------------------------------------------
+# Copyright 2022 Diamond Light Source Ltd.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+# ---------------------------------------------------------------------------
+# Created By  : Tomography Team at DLS <[email protected]>
+# Created Date: 02/June/2025
+# ---------------------------------------------------------------------------
+"""This is a collection of supplementary functions (utils) to perform various data checks"""
+
+from httomolibgpu import cupywrapper
+from typing import Optional
+
+cp = cupywrapper.cp
+cupy_run = cupywrapper.cupy_run
+
+import numpy as np
+
+from unittest.mock import Mock
+
+if cupy_run:
+    from httomolibgpu.cuda_kernels import load_cuda_module
+else:
+    load_cuda_module = Mock()
+
+
+def _naninfs_check(
+    data: cp.ndarray,
+    verbosity: bool = True,
+    method_name: Optional[str] = None,
+) -> cp.ndarray:
+    """
+    This function finds NaN's, +-Inf's in the input data and then prints the warnings and correct the data if correction is enabled.
+
+    Parameters
+    ----------
+    data : cp.ndarray
+        Input CuPy or Numpy array either float32 or uint16 data type.
+    verbosity : bool
+        If enabled, then the printing of the warning happens when data contains infs or nans
+    method_name : str, optional.
+        Method's name for which input data is tested.
+
+    Returns
+    -------
+    ndarray
+        Uncorrected or corrected (nans and infs converted to zeros) input array.
+    """
+    present_nans_infs_b = False
+
+    if cupy_run:
+        xp = cp.get_array_module(data)
+    else:
+        import numpy as xp
+
+    if xp.__name__ == "cupy":
+        input_type = data.dtype
+        if len(data.shape) == 2:
+            dy, dx = data.shape
+            dz = 1
+        else:
+            dz, dy, dx = data.shape
+
+        present_nans_infs = cp.zeros(shape=(1)).astype(cp.uint8)
+
+        block_x = 128
+        # setting grid/block parameters
+        block_dims = (block_x, 1, 1)
+        grid_x = (dx + block_x - 1) // block_x
+        grid_y = dy
+        grid_z = dz
+        grid_dims = (grid_x, grid_y, grid_z)
+        params = (data, dz, dy, dx, present_nans_infs)
+
+        kernel_args = "remove_nan_inf<{0}>".format(
+            "float" if input_type == "float32" else "unsigned short"
+        )
+
+        module = load_cuda_module("remove_nan_inf", name_expressions=[kernel_args])
+        remove_nan_inf_kernel = module.get_function(kernel_args)
+        remove_nan_inf_kernel(grid_dims, block_dims, params)
+
+        if present_nans_infs[0].get() == 1:
+            present_nans_infs_b = True
+    else:
+        if not np.all(np.isfinite(data)):
+            present_nans_infs_b = True
+            np.nan_to_num(data, copy=False, nan=0.0, posinf=0.0, neginf=0.0)
+
+    if present_nans_infs_b:
+        if verbosity:
+            print(
+                f"Warning!!! Input data to method: {method_name} contains Inf's or/and NaN's. This will be corrected but it is recommended to check the validity of input to the method."
+            )
+
+    return data
+
+
+def _zeros_check(
+    data: cp.ndarray,
+    verbosity: bool = True,
+    percentage_threshold: float = 50,
+    method_name: Optional[str] = None,
+) -> bool:
+    """
+    This function finds all zeros present in the data. If the amount of zeros is larger than percentage_threshold it prints the warning.
+
+    Parameters
+    ----------
+    data : cp.ndarray
+        Input CuPy or Numpy array.
+    verbosity : bool
+        If enabled, then the printing of the warning happens when data contains infs or nans.
+    percentage_threshold: float:
+        If the number of zeros in input data is more than the percentage of all data points, then print the data warning
+    method_name : str, optional.
+        Method's name for which input data is tested.
+
+    Returns
+    -------
+    bool
+        True if the data contains too many zeros
+    """
+    if cupy_run:
+        xp = cp.get_array_module(data)
+    else:
+        import numpy as xp
+
+    nonzero_elements_total = 1
+    for tot_elements_mult in data.shape:
+        nonzero_elements_total *= tot_elements_mult
+
+    warning_zeros = False
+    zero_elements_total = nonzero_elements_total - int(xp.count_nonzero(data))
+
+    if (zero_elements_total / nonzero_elements_total) * 100 >= percentage_threshold:
+        warning_zeros = True
+        if verbosity:
+            print(
+                f"Warning!!! Input data to method: {method_name} contains more than {percentage_threshold} percent of zeros."
+            )
+
+    return warning_zeros
+
+
+def data_checker(
+    data: cp.ndarray,
+    verbosity: bool = True,
+    method_name: Optional[str] = None,
+) -> bool:
+    """
+    Function that performs the variety of checks on input data, in some cases also correct the data and prints warnings.
+    Currently it checks for: the presence of infs and nans in data; the number of zero elements.
+
+    Parameters
+    ----------
+    data : xp.ndarray
+        Input CuPy or Numpy array either float32 or uint16 data type.
+    verbosity : bool
+        If enabled, then the printing of the warning happens when data contains infs or nans.
+    method_name : str, optional.
+        Method's name for which input data is tested.
+
+    Returns
+    -------
+    cp.ndarray
+        Returns corrected or not data array.
+    """
+
+    data = _naninfs_check(data, verbosity=verbosity, method_name=method_name)
+
+    _zeros_check(
+        data,
+        verbosity=verbosity,
+        percentage_threshold=50,
+        method_name=method_name,
+    )
+
+    return data

httomolibgpu/prep/alignment.py

@@ -35,6 +35,8 @@
 
 from typing import Dict, List, Tuple
 
+from httomolibgpu.misc.supp_func import data_checker
+
 __all__ = [
     "distortion_correction_proj_discorpy",
 ]
@@ -86,6 +88,10 @@ def distortion_correction_proj_discorpy(
     if len(data.shape) == 2:
         data = cp.expand_dims(data, axis=0)
 
+    data = data_checker(
+        data, verbosity=True, method_name="distortion_correction_proj_discorpy"
+    )
+
     # Get info from metadata txt file
     xcenter, ycenter, list_fact = _load_metadata_txt(metadata_path)