From f09f03b8a098593de0c9a17b728088a393006209 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Wed, 27 Jul 2022 18:08:13 +0100 Subject: [PATCH 1/3] inpainting module to take masks --- .../savu_hpc/savu_installer/savu_installer.sh | 4 +- savu/plugins/filters/inpainting/inpainting.py | 78 ++++++++++++------- .../filters/inpainting/inpainting_tools.py | 47 ++++++++--- 3 files changed, 90 insertions(+), 39 deletions(-) diff --git a/install/savu_hpc/savu_installer/savu_installer.sh b/install/savu_hpc/savu_installer/savu_installer.sh index ac3959e6f..82e35726b 100755 --- a/install/savu_hpc/savu_installer/savu_installer.sh +++ b/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 diff --git a/savu/plugins/filters/inpainting/inpainting.py b/savu/plugins/filters/inpainting/inpainting.py index a15ae8726..1ff902e8a 100644 --- a/savu/plugins/filters/inpainting/inpainting.py +++ b/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_LINCOMB @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) - pars = {'algorithm' : INPAINT_LINCOMB, \ - 'input' : input_temp,\ + if self.parameters['method'] == 'LINEARCOMB': + pars = {'input' : input_temp, 'maskData' : mask, 'number_of_iterations' : self.parameters['iterations'], 'windowsize_half' : self.parameters['windowsize_half'], 'sigma' : np.exp(self.parameters['sigma'])} - - (result, mask_upd) = INPAINT_LINCOMB(pars['input'], + (result, mask_upd) = INPAINT_LINCOMB(pars['input'], pars['maskData'], pars['number_of_iterations'], pars['windowsize_half'], pars['sigma']) + + 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'] + return 1 \ No newline at end of file diff --git a/savu/plugins/filters/inpainting/inpainting_tools.py b/savu/plugins/filters/inpainting/inpainting_tools.py index 251f56925..dc8ec70d3 100644 --- a/savu/plugins/filters/inpainting/inpainting_tools.py +++ b/savu/plugins/filters/inpainting/inpainting_tools.py @@ -9,27 +9,52 @@ 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: [LINEARCOMB, NONLOCAL_MARCH, DIFFUSION] + description: Choose inpainting method + default: LINEARCOMB iterations: visibility: basic - dtype: float - description: controls the smoothing level of the inpainted region. + dtype: int + description: the number of iterations to perform the inpainting (controls the smoothing level) default: 50 windowsize_half: - visibility: basic + visibility: intermediate dtype: int - description: half-size of the smoothing window. + description: half-size of the smoothing window, increase size for more smoothing default: 3 + dependency: + method: [LINEARCOMB] + search_window_increment: + visibility: advanced + dtype: int + description: the increment value with which the searching window is growing in iterations + default: 1 + dependency: + method: [NONLOCAL_MARCH] sigma: visibility: basic dtype: float - description: maximum value for the inpainted region. + description: maximum intensity value for the inpainted region. default: 0.5 - pattern: + dependency: + method: [LINEARCOMB] + regularisation_parameter: visibility: basic - dtype: str - description: Pattern to apply these to. - default: SINOGRAM - + dtype: float + 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] """ \ No newline at end of file From f61dc63a7fc1b1a507ffd3c403078213d519328e Mon Sep 17 00:00:00 2001 From: Yousef Moazzam Date: Fri, 29 Jul 2022 10:50:20 +0100 Subject: [PATCH 2/3] Add plugin to detect stripe artifacts in sinograms --- savu/plugins/ring_removal/detect_stripes.py | 153 ++++++++++++++++++ .../ring_removal/detect_stripes_tools.py | 57 +++++++ 2 files changed, 210 insertions(+) create mode 100644 savu/plugins/ring_removal/detect_stripes.py create mode 100644 savu/plugins/ring_removal/detect_stripes_tools.py diff --git a/savu/plugins/ring_removal/detect_stripes.py b/savu/plugins/ring_removal/detect_stripes.py new file mode 100644 index 000000000..6dc3ba2c7 --- /dev/null +++ b/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 + +""" + +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 diff --git a/savu/plugins/ring_removal/detect_stripes_tools.py b/savu/plugins/ring_removal/detect_stripes_tools.py new file mode 100644 index 000000000..1c56eac3b --- /dev/null +++ b/savu/plugins/ring_removal/detect_stripes_tools.py @@ -0,0 +1,57 @@ +from savu.plugins.plugin_tools import PluginTools + +class DetectStripesTools(PluginTools): + """ + A plugin to detect stripes in a sinogram and return a 2D binary mask that + specifies the position of the stripes. + """ + def define_parameters(self): + """ + size: + visibility: basic + dtype: int + description: Size of the median filter window. Greater is stronger. + default: 51 + snr: + visibility: basic + dtype: float + description: Ratio used to detect locations of stripes. + Greater is less sensitive. + default: 3.0 + + binary_dilation: + visibility: basic + dtype: bool + description: Apply binary dilation to the stripe mask + default: True + """ + + def citation(self): + """ + The code of stripe artifact detection is the implementation of the work + of Nghia T. Vo et al. taken from algorithm 4 in this paper. + bibtex: + @article{vo2018superior, + title = {Superior techniques for eliminating ring artifacts in X-ray micro-tomography}, + author={Vo, Nghia T and Atwood, Robert C and Drakopoulos, Michael}, + journal={Optics express}, + volume={26}, + number={22}, + pages={28396--28412}, + year={2018}, + publisher={Optical Society of America}} + endnote: + %0 Journal Article + %T Superior techniques for eliminating ring artifacts in X-ray micro-tomography + %A Vo, Nghia T + %A Atwood, Robert C + %A Drakopoulos, Michael + %J Optics express + %V 26 + %N 22 + %P 28396-28412 + %@ 1094-4087 + %D 2018 + %I Optical Society of America + doi: "10.1364/OE.26.028396" + """ From 2c12e61e1dd6ca885e97ea1be08b442a8b5fde43 Mon Sep 17 00:00:00 2001 From: dkazanc Date: Mon, 1 Aug 2022 10:56:29 +0100 Subject: [PATCH 3/3] correction to inpainting plugin --- savu/plugins/filters/inpainting/inpainting.py | 16 ++++++++-------- .../filters/inpainting/inpainting_tools.py | 19 ++++++------------- 2 files changed, 14 insertions(+), 21 deletions(-) diff --git a/savu/plugins/filters/inpainting/inpainting.py b/savu/plugins/filters/inpainting/inpainting.py index 1ff902e8a..150aa93e7 100644 --- a/savu/plugins/filters/inpainting/inpainting.py +++ b/savu/plugins/filters/inpainting/inpainting.py @@ -25,7 +25,7 @@ from savu.plugins.utils import register_plugin import numpy as np -from larix.methods.misc import INPAINT_NDF, INPAINT_NM, INPAINT_LINCOMB +from larix.methods.misc import INPAINT_NDF, INPAINT_NM, INPAINT_EUCL_WEIGHTED @register_plugin class Inpainting(Plugin, CpuPlugin): @@ -60,19 +60,19 @@ def process_frames(self, data): 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 - if self.parameters['method'] == 'LINEARCOMB': + 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, @@ -108,4 +108,4 @@ def nInput_datasets(self): return max(len(self.parameters['in_datasets']), 1) def nOutput_datasets(self): - return 1 \ No newline at end of file + return 1 diff --git a/savu/plugins/filters/inpainting/inpainting_tools.py b/savu/plugins/filters/inpainting/inpainting_tools.py index dc8ec70d3..a33bd2d01 100644 --- a/savu/plugins/filters/inpainting/inpainting_tools.py +++ b/savu/plugins/filters/inpainting/inpainting_tools.py @@ -14,21 +14,21 @@ def define_parameters(self): method: visibility: intermediate dtype: str - options: [LINEARCOMB, NONLOCAL_MARCH, DIFFUSION] + options: [INPAINT_EUCL_WEIGHTED, NONLOCAL_MARCH, DIFFUSION] description: Choose inpainting method - default: LINEARCOMB + default: INPAINT_EUCL_WEIGHTED iterations: visibility: basic dtype: int description: the number of iterations to perform the inpainting (controls the smoothing level) - default: 50 + default: 5 windowsize_half: visibility: intermediate dtype: int description: half-size of the smoothing window, increase size for more smoothing - default: 3 + default: 7 dependency: - method: [LINEARCOMB] + method: [INPAINT_EUCL_WEIGHTED] search_window_increment: visibility: advanced dtype: int @@ -36,13 +36,6 @@ def define_parameters(self): default: 1 dependency: method: [NONLOCAL_MARCH] - sigma: - visibility: basic - dtype: float - description: maximum intensity value for the inpainted region. - default: 0.5 - dependency: - method: [LINEARCOMB] regularisation_parameter: visibility: basic dtype: float @@ -57,4 +50,4 @@ def define_parameters(self): default: 0.00001 dependency: method: [DIFFUSION] - """ \ No newline at end of file + """