some WIP on impainting, merging to finalise master

install/savu_hpc/savu_installer/savu_installer.sh

@@ -281,8 +281,8 @@ if [ ! $test_flag ]; then
 
     if [ $EXPLICIT_FILE = false ]; then
     echo "Installing mpi4py from savu-dep conda channel"
-    export VERSION_BUILD_MPI4PI=$mpi4py_version"_openmpi_"$openmpi_version
-    conda install --yes -c savu-dep mpi4py=$VERSION_BUILD_MPI4PI --no-deps
+    export VERSION_BUILD_MPI4PY=$mpi4py_version"_openmpi_"$openmpi_version
+    conda install --yes -c savu-dep mpi4py=$VERSION_BUILD_MPI4PY --no-deps
 
     echo "Installing hdf5 from savu-dep conda channel"
     export VERSION_BUILD_HDF5=$hdf5_version"_openmpi_"$openmpi_version

savu/plugins/filters/inpainting/inpainting.py

@@ -25,19 +25,13 @@
 from savu.plugins.utils import register_plugin
 
 import numpy as np
-from larix.methods.misc import INPAINT_LINCOMB
+from larix.methods.misc import INPAINT_NDF, INPAINT_NM, INPAINT_EUCL_WEIGHTED
 
 @register_plugin
 class Inpainting(Plugin, CpuPlugin):
     """
     A plugin to apply 2D/3D inpainting technique to data. If there is a chunk of
-    data missing or one needs to inpaint some data features.
-
-    :u*param mask_range: mask for inpainting is set as a threhsold in a range. Default: [1.0,100].
-    :u*param iterations: controls the smoothing level of the inpainted region. Default: 50.
-    :u*param windowsize_half: half-size of the smoothing window. Default: 3.
-    :u*param sigma: maximum value for the inpainted region. Default: 0.5.
-    :u*param pattern: pattern to apply this to. Default: "SINOGRAM".
+    data missing or one needs to inpaint some missing data features.
     """
 
     def __init__(self):
@@ -47,39 +41,71 @@ def setup(self):
         in_dataset, out_dataset = self.get_datasets()
         out_dataset[0].create_dataset(in_dataset[0])
         in_pData, out_pData = self.get_plugin_datasets()
-        in_pData[0].plugin_data_setup(self.parameters['pattern'], 'single')
-        out_pData[0].plugin_data_setup(self.parameters['pattern'], 'single')
+        pattern = list(in_dataset[0].get_data_patterns().keys())[0]
+        in_pData[0].plugin_data_setup(pattern, 'single')
+        out_pData[0].plugin_data_setup(pattern, 'single')
+        if self.nInput_datasets() > 1:
+            in_pData[1].plugin_data_setup(pattern, 'single')
         self.mask_range = self.parameters['mask_range']
 
     def process_frames(self, data):
-        input_temp = np.float32(data[0])
-        mask = np.zeros(np.shape(input_temp))
-        indices = np.where(np.isnan(input_temp))
-        input_temp[indices] = 0.0
-        mask[(input_temp >= self.mask_range[0]) & (input_temp < self.mask_range[1])] = 1.0
-        input_temp = np.ascontiguousarray(input_temp, dtype=np.float32);
-        mask = np.ascontiguousarray(mask, dtype=np.uint8);
+        input_temp = np.float32(data[0]) # get the data in for inpaintng
+        if self.nInput_datasets() > 1:
+            mask = np.uint8(data[1]) # get the mask if given
+        else:
+            # build the mask if not given based on threshold values
+            mask = np.uint8(np.zeros(np.shape(input_temp)))
+            indices = np.where(np.isnan(input_temp))
+            input_temp[indices] = 0
+            mask[(input_temp >= self.mask_range[0]) & (input_temp < self.mask_range[1])] = 1
+        input_temp = np.ascontiguousarray(input_temp, dtype=np.float32)
+        mask = np.ascontiguousarray(mask, dtype=np.uint8)
+        # modify input to crop out masked values
+        input_temp[mask == 1] = 0.0
 
-        pars = {'algorithm' : INPAINT_LINCOMB, \
-                'input' : input_temp,\
+        if self.parameters['method'] == 'INPAINT_EUCL_WEIGHTED':
+            pars = {'input' : input_temp,
                 'maskData' : mask,
                 'number_of_iterations' : self.parameters['iterations'],
-                'windowsize_half' : self.parameters['windowsize_half'],
-                'sigma' : np.exp(self.parameters['sigma'])}
-
+                'windowsize_half' : self.parameters['windowsize_half']}
 
-        (result, mask_upd) = INPAINT_LINCOMB(pars['input'],
+            (result, mask_upd) = INPAINT_EUCL_WEIGHTED(pars['input'],
                               pars['maskData'],
                               pars['number_of_iterations'],
-                              pars['windowsize_half'],
-                              pars['sigma'])
+                              pars['windowsize_half'])
+
+        elif self.parameters['method'] == 'NONLOCAL_MARCH':
+            pars = {'input' : input_temp,
+                'maskData' : mask,
+                'SW_increment' : self.parameters['search_window_increment'],
+                'number_of_iterations' : self.parameters['iterations']}
+
+            (result, mask_upd) = INPAINT_NM(pars['input'],
+                                             pars['maskData'],
+                                             pars['SW_increment'],
+                                             pars['number_of_iterations'])
+        elif self.parameters['method'] == 'DIFFUSION':
+            pars = {'input': input_temp,
+                    'maskData': mask,
+                    'regularisation_parameter': self.parameters['regularisation_parameter'], \
+                    'edge_parameter': 0,
+                    'number_of_iterations': self.parameters['iterations'],
+                    'time_marching_parameter': self.parameters['time_marching_parameter'],
+                    'penalty_type': 1
+                    }
+            result = INPAINT_NDF(pars['input'],
+                                         pars['maskData'],
+                                         pars['regularisation_parameter'],
+                                         pars['edge_parameter'],
+                                         pars['number_of_iterations'],
+                                         pars['time_marching_parameter'],
+                                         pars['penalty_type'])
+        else:
+            print("The inpainting method is not selected!")
         return result
 
     def nInput_datasets(self):
-        return 1
+        return max(len(self.parameters['in_datasets']), 1)
 
     def nOutput_datasets(self):
         return 1
-
-    def get_plugin_pattern(self):
-        return self.parameters['pattern']

savu/plugins/filters/inpainting/inpainting_tools.py

@@ -9,27 +9,45 @@ def define_parameters(self):
         mask_range:
               visibility: basic
               dtype: list[float,float]
-              description: mask for inpainting is set as a threshold in a range.
+              description: generate mask as a threshold of the input image given the intensity range (ignored if the mask is provided).
               default: [1.0,100]
+        method:
+              visibility: intermediate
+              dtype: str
+              options: [INPAINT_EUCL_WEIGHTED, NONLOCAL_MARCH, DIFFUSION]
+              description: Choose inpainting method
+              default: INPAINT_EUCL_WEIGHTED
         iterations:
               visibility: basic
-              dtype: float
-              description: controls the smoothing level of the inpainted region.
-              default: 50
+              dtype: int
+              description: the number of iterations to perform the inpainting (controls the smoothing level)
+              default: 5
         windowsize_half:
-              visibility: basic
+              visibility: intermediate
+              dtype: int
+              description: half-size of the smoothing window, increase size for more smoothing
+              default: 7
+              dependency:
+                method: [INPAINT_EUCL_WEIGHTED]
+        search_window_increment:
+              visibility: advanced
               dtype: int
-              description: half-size of the smoothing window.
-              default: 3
-        sigma:
+              description: the increment value with which the searching window is growing in iterations
+              default: 1
+              dependency:
+                method: [NONLOCAL_MARCH]
+        regularisation_parameter:
               visibility: basic
               dtype: float
-              description: maximum value for the inpainted region.
-              default: 0.5
-        pattern:
-              visibility: basic
-              dtype: str
-              description: Pattern to apply these to.
-              default: SINOGRAM
-
-        """
+              description: regularisation parameter for diffusion
+              default: 1000
+              dependency:
+                method: [DIFFUSION]
+        time_marching_parameter:
+              visibility: intermediate
+              dtype: float
+              description: ensuring convergence of the diffusion scheme
+              default: 0.00001
+              dependency:
+                method: [DIFFUSION]
+        """

savu/plugins/ring_removal/detect_stripes.py

@@ -0,0 +1,153 @@
+# 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.
+
+"""
+.. module:: detect_stripes
+   :platform: Unix
+   :synopsis: Method working in the sinogram space to detect stripe artifacts.
+.. moduleauthor:: Nghia Vo, Yousef Moazzam <scientificsoftware@diamond.ac.uk>
+
+"""
+
+from savu.plugins.plugin import Plugin
+from savu.plugins.driver.cpu_plugin import CpuPlugin
+from savu.plugins.utils import register_plugin
+
+import numpy as np
+from scipy.ndimage import median_filter
+from scipy.ndimage import binary_dilation
+from scipy.ndimage import uniform_filter1d
+
+
+@register_plugin
+class DetectStripes(Plugin, CpuPlugin):
+    """
+    A plugin to detect stripes in a sinogram and return a 2D binary mask that
+    specifies the position of the stripes.
+    """
+    def __init__(self):
+        super(DetectStripes, self).__init__('DetectStripes')
+
+    def setup(self):
+        in_dataset, out_dataset = self.get_datasets()
+        in_pData, out_pData = self.get_plugin_datasets()
+        in_pData[0].plugin_data_setup('SINOGRAM',self.get_max_frames())
+        out_dataset[0].create_dataset(in_dataset[0])
+        out_pData[0].plugin_data_setup('SINOGRAM', self.get_max_frames())
+
+    def pre_process(self):
+        in_pData = self.get_plugin_in_datasets()
+        width_dim = \
+            in_pData[0].get_data_dimension_by_axis_label('detector_x')
+        height_dim = \
+            in_pData[0].get_data_dimension_by_axis_label('rotation_angle')
+        sino_shape = list(in_pData[0].get_shape())
+        self.width1 = sino_shape[width_dim]
+        self.height1 = sino_shape[height_dim]
+        listindex = np.arange(0.0, self.height1, 1.0)
+        self.matindex = np.tile(listindex, (self.width1, 1))
+        self.size = np.clip(np.int16(self.parameters['size']), 1,
+                            self.width1 - 1)
+        self.snr = np.clip(np.float32(self.parameters['snr']), 1.0, None)
+
+    def process_frames(self, data):
+        sinogram = np.copy(data[0])
+        # begin pre-processing sinogram as preparation for the stripe detection
+        # algorithm
+        #
+        # NOTE: the pre-processing steps here are specifically for unresponsive
+        # and fluctuating stripe artifacts
+        sinosmoothed = np.apply_along_axis(uniform_filter1d, 0, sinogram, 10)
+        listdiff = np.sum(np.abs(sinogram - sinosmoothed), axis=0)
+        nmean = np.mean(listdiff)
+        listdiffbck = median_filter(listdiff, self.size)
+        listdiffbck[listdiffbck == 0.0] = nmean
+        listfact = listdiff / listdiffbck
+        # finish pre-processing sinogram
+
+        listmask = self.detect_stripe(listfact, self.snr)
+        if self.parameters['binary_dilation']:
+            listmask = binary_dilation(listmask, iterations=1).astype(
+                listmask.dtype)
+            listmask[0:2] = 0.0
+            listmask[-2:] = 0.0
+
+        stripe_locations = np.squeeze(np.argwhere(listmask==1))
+        mask_2d = np.zeros(sinogram.shape)
+        # separate detected columns in `listmask` into groups that correspond to
+        # the stripes that have been detected
+        stripes = []
+        for i in range(len(stripe_locations)):
+            if i == 0:
+                stripe = [stripe_locations[i]]
+                continue
+            elif i == len(stripe_locations) - 1:
+                stripe.append(stripe_locations[i])
+                stripes.append(stripe)
+                continue
+            if stripe_locations[i] != stripe_locations[i-1] + 1:
+                stripe.append(stripe_locations[i-1])
+                stripes.append(stripe)
+                stripe = [stripe_locations[i]]
+
+        # apply 1's in 2D mask at stripe locations
+        for stripe in stripes:
+            mask_2d[:, stripe[0]:stripe[1]] = 1
+
+        return mask_2d
+
+    def detect_stripe(self, listdata, snr):
+        """Algorithm 4 in the paper. To locate stripe positions.
+
+        Parameters
+        ----------
+        listdata : 1D normalized array.
+        snr : Ratio (>1.0) used to detect stripe locations.
+
+        Returns
+        -------
+        listmask : 1D binary mask.
+        """
+        numdata = len(listdata)
+        listsorted = np.sort(listdata)[::-1]
+        xlist = np.arange(0, numdata, 1.0)
+        ndrop = np.int16(0.25 * numdata)
+        (_slope, _intercept) = np.polyfit(
+            xlist[ndrop:-ndrop - 1], listsorted[ndrop:-ndrop - 1], 1)
+        numt1 = _intercept + _slope * xlist[-1]
+        noiselevel = np.abs(numt1 - _intercept)
+        if noiselevel == 0.0:
+            raise ValueError(
+                "The method doesn't work on noise-free data. If you " \
+                "apply the method on simulated data, please add" \
+                " noise!")
+        val1 = np.abs(listsorted[0] - _intercept) / noiselevel
+        val2 = np.abs(listsorted[-1] - numt1) / noiselevel
+        listmask = np.zeros_like(listdata)
+        if val1 >= snr:
+            upper_thresh = _intercept + noiselevel * snr * 0.5
+            listmask[listdata > upper_thresh] = 1.0
+        if val2 >= snr:
+            lower_thresh = numt1 - noiselevel * snr * 0.5
+            listmask[listdata <= lower_thresh] = 1.0
+        return listmask
+
+    def get_max_frames(self):
+        return 'single'
+
+    def nInput_datasets(self):
+        return 1
+
+    def nOutput_datasets(self):
+        return 1