summaryrefslogtreecommitdiffstats
path: root/Wrappers
diff options
context:
space:
mode:
Diffstat (limited to 'Wrappers')
-rw-r--r--Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m9
-rw-r--r--Wrappers/Matlab/demos/demoMatlab_denoise.m9
-rw-r--r--Wrappers/Matlab/demos/demoMatlab_inpaint.m35
-rw-r--r--Wrappers/Matlab/mex_compile/compileCPU_mex.m16
-rw-r--r--Wrappers/Matlab/mex_compile/regularisers_CPU/NonlDiff_Inp.c101
-rw-r--r--Wrappers/Matlab/mex_compile/regularisers_CPU/NonlocalMarching_Inpaint.c82
-rw-r--r--Wrappers/Python/CMakeLists.txt5
-rw-r--r--Wrappers/Python/ccpi/filters/regularisers.py9
-rw-r--r--Wrappers/Python/demos/demo_cpu_inpainters.py192
-rw-r--r--Wrappers/Python/demos/demo_cpu_regularisers.py43
-rw-r--r--Wrappers/Python/demos/demo_gpu_regularisers.py18
-rw-r--r--Wrappers/Python/setup-regularisers.py.in3
-rw-r--r--Wrappers/Python/src/cpu_regularisers.pyx90
-rw-r--r--Wrappers/Python/src/gpu_regularisers.pyx10
14 files changed, 578 insertions, 44 deletions
diff --git a/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m b/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m
index 5a54d18..c087433 100644
--- a/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m
+++ b/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m
@@ -1,8 +1,9 @@
% Volume (3D) denoising demo using CCPi-RGL
-clear
-close all
-addpath('../mex_compile/installed');
-addpath('../../../data/');
+clear; close all
+Path1 = sprintf(['..' filesep 'mex_compile' filesep 'installed'], 1i);
+Path2 = sprintf(['..' filesep '..' filesep '..' filesep 'data' filesep], 1i);
+addpath(Path1);
+addpath(Path2);
N = 512;
slices = 30;
diff --git a/Wrappers/Matlab/demos/demoMatlab_denoise.m b/Wrappers/Matlab/demos/demoMatlab_denoise.m
index 151a604..d93f477 100644
--- a/Wrappers/Matlab/demos/demoMatlab_denoise.m
+++ b/Wrappers/Matlab/demos/demoMatlab_denoise.m
@@ -1,8 +1,9 @@
% Image (2D) denoising demo using CCPi-RGL
-clear
-close all
-addpath('../mex_compile/installed');
-addpath('../../../data/');
+clear; close all
+Path1 = sprintf(['..' filesep 'mex_compile' filesep 'installed'], 1i);
+Path2 = sprintf(['..' filesep '..' filesep '..' filesep 'data' filesep], 1i);
+addpath(Path1);
+addpath(Path2);
Im = double(imread('lena_gray_512.tif'))/255; % loading image
u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
diff --git a/Wrappers/Matlab/demos/demoMatlab_inpaint.m b/Wrappers/Matlab/demos/demoMatlab_inpaint.m
new file mode 100644
index 0000000..66f9c15
--- /dev/null
+++ b/Wrappers/Matlab/demos/demoMatlab_inpaint.m
@@ -0,0 +1,35 @@
+% Image (2D) inpainting demo using CCPi-RGL
+clear; close all
+Path1 = sprintf(['..' filesep 'mex_compile' filesep 'installed'], 1i);
+Path2 = sprintf(['..' filesep '..' filesep '..' filesep 'data' filesep], 1i);
+addpath(Path1);
+addpath(Path2);
+
+load('SinoInpaint.mat');
+Sinogram = Sinogram./max(Sinogram(:));
+Sino_mask = Sinogram.*(1-single(Mask));
+figure;
+subplot(1,2,1); imshow(Sino_mask, [0 1]); title('Missing data sinogram');
+subplot(1,2,2); imshow(Mask, [0 1]); title('Mask');
+%%
+fprintf('Inpaint using Linear-Diffusion model (CPU) \n');
+iter_diff = 5000; % number of diffusion iterations
+lambda_regDiff = 6000; % regularisation for the diffusivity
+sigmaPar = 0.0; % edge-preserving parameter
+tau_param = 0.000075; % time-marching constant
+tic; u_diff = NonlDiff_Inp(single(Sino_mask), Mask, lambda_regDiff, sigmaPar, iter_diff, tau_param); toc;
+figure; imshow(u_diff, [0 1]); title('Linear-Diffusion inpainted sinogram (CPU)');
+%%
+fprintf('Inpaint using Nonlinear-Diffusion model (CPU) \n');
+iter_diff = 1500; % number of diffusion iterations
+lambda_regDiff = 80; % regularisation for the diffusivity
+sigmaPar = 0.00009; % edge-preserving parameter
+tau_param = 0.000008; % time-marching constant
+tic; u_diff = NonlDiff_Inp(single(Sino_mask), Mask, lambda_regDiff, sigmaPar, iter_diff, tau_param, 'Huber'); toc;
+figure; imshow(u_diff, [0 1]); title('Non-Linear Diffusion inpainted sinogram (CPU)');
+%%
+fprintf('Inpaint using Nonlocal Vertical Marching model (CPU) \n');
+Increment = 1; % linear increment for the searching window
+tic; [u_nom,maskupd] = NonlocalMarching_Inpaint(single(Sino_mask), Mask, Increment); toc;
+figure; imshow(u_nom, [0 1]); title('NVM inpainted sinogram (CPU)');
+%% \ No newline at end of file
diff --git a/Wrappers/Matlab/mex_compile/compileCPU_mex.m b/Wrappers/Matlab/mex_compile/compileCPU_mex.m
index ee0d99e..b232f33 100644
--- a/Wrappers/Matlab/mex_compile/compileCPU_mex.m
+++ b/Wrappers/Matlab/mex_compile/compileCPU_mex.m
@@ -2,9 +2,11 @@
pathcopyFrom = sprintf(['..' filesep '..' filesep '..' filesep 'Core' filesep 'regularisers_CPU'], 1i);
pathcopyFrom1 = sprintf(['..' filesep '..' filesep '..' filesep 'Core' filesep 'CCPiDefines.h'], 1i);
+pathcopyFrom2 = sprintf(['..' filesep '..' filesep '..' filesep 'Core' filesep 'inpainters_CPU'], 1i);
copyfile(pathcopyFrom, 'regularisers_CPU');
copyfile(pathcopyFrom1, 'regularisers_CPU');
+copyfile(pathcopyFrom2, 'regularisers_CPU');
cd regularisers_CPU
@@ -32,9 +34,17 @@ movefile('NonlDiff.mex*',Pathmove);
mex TV_energy.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
movefile('TV_energy.mex*',Pathmove);
+%############Inpainters##############%
+mex NonlDiff_Inp.c Diffusion_Inpaint_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+movefile('NonlDiff_Inp.mex*',Pathmove);
+
+mex NonlocalMarching_Inpaint.c NonlocalMarching_Inpaint_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+movefile('NonlocalMarching_Inpaint.mex*',Pathmove);
+
delete SB_TV_core* ROF_TV_core* FGP_TV_core* FGP_dTV_core* TNV_core* utils* Diffusion_core* CCPiDefines.h
-fprintf('%s \n', 'All successfully compiled!');
+delete Diffusion_Inpaint_core* NonlocalMarching_Inpaint_core*
+fprintf('%s \n', 'Regularisers successfully compiled!');
-pathA = sprintf(['..' filesep '..' filesep], 1i);
-cd(pathA);
+pathA2 = sprintf(['..' filesep '..' filesep], 1i);
+cd(pathA2);
cd demos
diff --git a/Wrappers/Matlab/mex_compile/regularisers_CPU/NonlDiff_Inp.c b/Wrappers/Matlab/mex_compile/regularisers_CPU/NonlDiff_Inp.c
new file mode 100644
index 0000000..eaab4a7
--- /dev/null
+++ b/Wrappers/Matlab/mex_compile/regularisers_CPU/NonlDiff_Inp.c
@@ -0,0 +1,101 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * 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.
+ */
+#include "matrix.h"
+#include "mex.h"
+#include "Diffusion_Inpaint_core.h"
+
+/* C-OMP implementation of linear and nonlinear diffusion [1,2] for inpainting task (2D/3D case)
+ * The minimisation is performed using explicit scheme.
+ *
+ * Input Parameters:
+ * 1. Image/volume to inpaint
+ * 2. Inpainting Mask of the same size as (1) in 'unsigned char' format (ones mark the region to inpaint, zeros belong to the data)
+ * 3. lambda - regularization parameter
+ * 4. Edge-preserving parameter (sigma), when sigma equals to zero nonlinear diffusion -> linear diffusion
+ * 5. Number of iterations, for explicit scheme >= 150 is recommended
+ * 6. tau - time-marching step for explicit scheme
+ * 7. Penalty type: 1 - Huber, 2 - Perona-Malik, 3 - Tukey Biweight
+ *
+ * Output:
+ * [1] Inpainted image/volume
+ *
+ * This function is based on the paper by
+ * [1] Perona, P. and Malik, J., 1990. Scale-space and edge detection using anisotropic diffusion. IEEE Transactions on pattern analysis and machine intelligence, 12(7), pp.629-639.
+ * [2] Black, M.J., Sapiro, G., Marimont, D.H. and Heeger, D., 1998. Robust anisotropic diffusion. IEEE Transactions on image processing, 7(3), pp.421-432.
+ */
+
+void mexFunction(
+ int nlhs, mxArray *plhs[],
+ int nrhs, const mxArray *prhs[])
+
+{
+ int number_of_dims, iter_numb, dimX, dimY, dimZ, penaltytype, i, inpaint_elements;
+ const int *dim_array;
+ const int *dim_array2;
+ float *Input, *Output=NULL, lambda, tau, sigma;
+ unsigned char *Mask;
+
+ dim_array = mxGetDimensions(prhs[0]);
+ dim_array2 = mxGetDimensions(prhs[1]);
+ number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+
+ /*Handling Matlab input data*/
+ Input = (float *) mxGetData(prhs[0]);
+ Mask = (unsigned char *) mxGetData(prhs[1]); /* MASK */
+ lambda = (float) mxGetScalar(prhs[2]); /* regularization parameter */
+ sigma = (float) mxGetScalar(prhs[3]); /* Edge-preserving parameter */
+ iter_numb = 300; /* iterations number */
+ tau = 0.025; /* marching step parameter */
+ penaltytype = 1; /* Huber penalty by default */
+
+ if ((nrhs < 4) || (nrhs > 7)) mexErrMsgTxt("At least 4 parameters is required, all parameters are: Image(2D/3D), Mask(2D/3D), Regularisation parameter, Edge-preserving parameter, iterations number, time-marching constant, penalty type - Huber, PM or Tukey");
+ if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7)) iter_numb = (int) mxGetScalar(prhs[4]); /* iterations number */
+ if ((nrhs == 6) || (nrhs == 7)) tau = (float) mxGetScalar(prhs[5]); /* marching step parameter */
+ if (nrhs == 7) {
+ char *penalty_type;
+ penalty_type = mxArrayToString(prhs[6]); /* Huber, PM or Tukey 'Huber' is the default */
+ if ((strcmp(penalty_type, "Huber") != 0) && (strcmp(penalty_type, "PM") != 0) && (strcmp(penalty_type, "Tukey") != 0)) mexErrMsgTxt("Choose penalty: 'Huber', 'PM' or 'Tukey',");
+ if (strcmp(penalty_type, "Huber") == 0) penaltytype = 1; /* enable 'Huber' penalty */
+ if (strcmp(penalty_type, "PM") == 0) penaltytype = 2; /* enable Perona-Malik penalty */
+ if (strcmp(penalty_type, "Tukey") == 0) penaltytype = 3; /* enable Tikey Biweight penalty */
+ mxFree(penalty_type);
+ }
+
+ if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+ if (mxGetClassID(prhs[1]) != mxUINT8_CLASS) {mexErrMsgTxt("The mask must be in uint8 precision");}
+
+ dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+
+ /* output arrays*/
+ if (number_of_dims == 2) {
+ dimZ = 1; /*2D case*/
+ /* output image/volume */
+ if ((dimX != dim_array2[0]) || (dimY != dim_array2[1])) mexErrMsgTxt("Input image and the provided mask are of different dimensions!");
+ Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+ }
+ if (number_of_dims == 3) {
+ if ((dimX != dim_array2[0]) || (dimY != dim_array2[1]) || (dimZ != dim_array2[2])) mexErrMsgTxt("Input image and the provided mask are of different dimensions!");
+ Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+ }
+
+ inpaint_elements = 0;
+ for (i=0; i<dimY*dimX*dimZ; i++) if (Mask[i] == 1) inpaint_elements++;
+ if (inpaint_elements == 0) mexErrMsgTxt("The mask is full of zeros, nothing to inpaint");
+ Diffusion_Inpaint_CPU_main(Input, Mask, Output, lambda, sigma, iter_numb, tau, penaltytype, dimX, dimY, dimZ);
+} \ No newline at end of file
diff --git a/Wrappers/Matlab/mex_compile/regularisers_CPU/NonlocalMarching_Inpaint.c b/Wrappers/Matlab/mex_compile/regularisers_CPU/NonlocalMarching_Inpaint.c
new file mode 100644
index 0000000..36cf05c
--- /dev/null
+++ b/Wrappers/Matlab/mex_compile/regularisers_CPU/NonlocalMarching_Inpaint.c
@@ -0,0 +1,82 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * 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.
+ */
+#include "matrix.h"
+#include "mex.h"
+#include "NonlocalMarching_Inpaint_core.h"
+
+/* C-OMP implementation of Nonlocal Vertical Marching inpainting method (2D case)
+ * The method is heuristic but computationally efficent (especially for larger images).
+ * It developed specifically to smoothly inpaint horizontal or inclined missing data regions in sinograms
+ * The method WILL not work satisfactory if you have lengthy vertical stripes of missing data
+ *
+ * Input:
+ * 1. 2D image or sinogram [REQUIRED]
+ * 2. Mask of the same size as A in 'unsigned char' format (ones mark the region to inpaint, zeros belong to the data) [REQUIRED]
+ * 3. Linear increment to increase searching window size in iterations, values from 1-3 is a good choice [OPTIONAL, default 1]
+ * 4. Number of iterations [OPTIONAL, default - calculate based on the mask]
+ *
+ * Output:
+ * 1. Inpainted sinogram
+ * 2. updated mask
+ * Reference: TBA
+ */
+
+void mexFunction(
+ int nlhs, mxArray *plhs[],
+ int nrhs, const mxArray *prhs[])
+
+{
+ int number_of_dims, dimX, dimY, dimZ, iterations, SW_increment;
+ const int *dim_array;
+ const int *dim_array2;
+ float *Input, *Output=NULL;
+ unsigned char *Mask, *Mask_upd=NULL;
+
+ dim_array = mxGetDimensions(prhs[0]);
+ dim_array2 = mxGetDimensions(prhs[1]);
+ number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+
+ /*Handling Matlab input data*/
+ Input = (float *) mxGetData(prhs[0]);
+ Mask = (unsigned char *) mxGetData(prhs[1]); /* MASK */
+ SW_increment = 1;
+ iterations = 0;
+
+ if ((nrhs < 2) || (nrhs > 4)) mexErrMsgTxt("At least 4 parameters is required, all parameters are: Image(2D/3D), Mask(2D/3D), Linear increment, Iterations number");
+ if ((nrhs == 3) || (nrhs == 4)) SW_increment = (int) mxGetScalar(prhs[2]); /* linear increment */
+ if ((nrhs == 4)) iterations = (int) mxGetScalar(prhs[3]); /* iterations number */
+
+ if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+ if (mxGetClassID(prhs[1]) != mxUINT8_CLASS) {mexErrMsgTxt("The mask must be in uint8 precision");}
+
+ dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+
+ /* output arrays*/
+ if (number_of_dims == 2) {
+ dimZ = 1; /*2D case*/
+ /* output image/volume */
+ if ((dimX != dim_array2[0]) || (dimY != dim_array2[1])) mexErrMsgTxt("Input image and the provided mask are of different dimensions!");
+ Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+ Mask_upd = (unsigned char*)mxGetPr(plhs[1] = mxCreateNumericArray(2, dim_array, mxUINT8_CLASS, mxREAL));
+ }
+ if (number_of_dims == 3) {
+ mexErrMsgTxt("Currently 2D supported only");
+ }
+ NonlocalMarching_Inpaint_main(Input, Mask, Output, Mask_upd, SW_increment, iterations, 0, dimX, dimY, dimZ);
+} \ No newline at end of file
diff --git a/Wrappers/Python/CMakeLists.txt b/Wrappers/Python/CMakeLists.txt
index fb00706..7833b54 100644
--- a/Wrappers/Python/CMakeLists.txt
+++ b/Wrappers/Python/CMakeLists.txt
@@ -81,8 +81,3 @@ else()
endif()
configure_file("setup-regularisers.py.in" "setup-regularisers.py")
-
-
-#add_executable(regulariser_test ${CMAKE_CURRENT_SOURCE_DIR}/test/test_regulariser.cpp)
-
-#target_link_libraries (regulariser_test LINK_PUBLIC regularisers_lib)
diff --git a/Wrappers/Python/ccpi/filters/regularisers.py b/Wrappers/Python/ccpi/filters/regularisers.py
index eec8c4d..a07b39a 100644
--- a/Wrappers/Python/ccpi/filters/regularisers.py
+++ b/Wrappers/Python/ccpi/filters/regularisers.py
@@ -2,7 +2,7 @@
script which assigns a proper device core function based on a flag ('cpu' or 'gpu')
"""
-from ccpi.filters.cpu_regularisers_cython import TV_ROF_CPU, TV_FGP_CPU, TV_SB_CPU, dTV_FGP_CPU, TNV_CPU, NDF_CPU
+from ccpi.filters.cpu_regularisers import TV_ROF_CPU, TV_FGP_CPU, TV_SB_CPU, dTV_FGP_CPU, TNV_CPU, NDF_CPU, NDF_INPAINT_CPU, NVM_INPAINT_CPU
from ccpi.filters.gpu_regularisers import TV_ROF_GPU, TV_FGP_GPU, TV_SB_GPU, dTV_FGP_GPU, NDF_GPU
def ROF_TV(inputData, regularisation_parameter, iterations,
@@ -110,3 +110,10 @@ def NDF(inputData, regularisation_parameter, edge_parameter, iterations,
else:
raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
.format(device))
+def NDF_INP(inputData, maskData, regularisation_parameter, edge_parameter, iterations,
+ time_marching_parameter, penalty_type):
+ return NDF_INPAINT_CPU(inputData, maskData, regularisation_parameter,
+ edge_parameter, iterations, time_marching_parameter, penalty_type)
+
+def NVM_INP(inputData, maskData, SW_increment, iterations):
+ return NVM_INPAINT_CPU(inputData, maskData, SW_increment, iterations)
diff --git a/Wrappers/Python/demos/demo_cpu_inpainters.py b/Wrappers/Python/demos/demo_cpu_inpainters.py
new file mode 100644
index 0000000..3b4191b
--- /dev/null
+++ b/Wrappers/Python/demos/demo_cpu_inpainters.py
@@ -0,0 +1,192 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Demonstration of CPU inpainters
+@authors: Daniil Kazantsev, Edoardo Pasca
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+import timeit
+from scipy import io
+from ccpi.filters.regularisers import NDF_INP, NVM_INP
+from qualitymetrics import rmse
+###############################################################################
+def printParametersToString(pars):
+ txt = r''
+ for key, value in pars.items():
+ if key== 'algorithm' :
+ txt += "{0} = {1}".format(key, value.__name__)
+ elif key == 'input':
+ txt += "{0} = {1}".format(key, np.shape(value))
+ elif key == 'maskData':
+ txt += "{0} = {1}".format(key, np.shape(value))
+ else:
+ txt += "{0} = {1}".format(key, value)
+ txt += '\n'
+ return txt
+###############################################################################
+
+# read sinogram and the mask
+filename = os.path.join(".." , ".." , ".." , "data" ,"SinoInpaint.mat")
+sino = io.loadmat(filename)
+sino_full = sino.get('Sinogram')
+Mask = sino.get('Mask')
+[angles_dim,detectors_dim] = sino_full.shape
+sino_full = sino_full/np.max(sino_full)
+#apply mask to sinogram
+sino_cut = sino_full*(1-Mask)
+#sino_cut_new = np.zeros((angles_dim,detectors_dim),'float32')
+#sino_cut_new = sino_cut.copy(order='c')
+#sino_cut_new[:] = sino_cut[:]
+sino_cut_new = np.ascontiguousarray(sino_cut, dtype=np.float32);
+#mask = np.zeros((angles_dim,detectors_dim),'uint8')
+#mask =Mask.copy(order='c')
+#mask[:] = Mask[:]
+mask = np.ascontiguousarray(Mask, dtype=np.uint8);
+
+plt.figure(1)
+plt.subplot(121)
+plt.imshow(sino_cut_new,vmin=0.0, vmax=1)
+plt.title('Missing Data sinogram')
+plt.subplot(122)
+plt.imshow(mask)
+plt.title('Mask')
+plt.show()
+#%%
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+print ("___Inpainting using linear diffusion (2D)__")
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+
+## plot
+fig = plt.figure(2)
+plt.suptitle('Performance of linear inpainting using the CPU')
+a=fig.add_subplot(1,2,1)
+a.set_title('Missing data sinogram')
+imgplot = plt.imshow(sino_cut_new,cmap="gray")
+
+# set parameters
+pars = {'algorithm' : NDF_INP, \
+ 'input' : sino_cut_new,\
+ 'maskData' : mask,\
+ 'regularisation_parameter':5000,\
+ 'edge_parameter':0,\
+ 'number_of_iterations' :5000 ,\
+ 'time_marching_parameter':0.000075,\
+ 'penalty_type':0
+ }
+
+start_time = timeit.default_timer()
+ndf_inp_linear = NDF_INP(pars['input'],
+ pars['maskData'],
+ pars['regularisation_parameter'],
+ pars['edge_parameter'],
+ pars['number_of_iterations'],
+ pars['time_marching_parameter'],
+ pars['penalty_type'])
+
+rms = rmse(sino_full, ndf_inp_linear)
+pars['rmse'] = rms
+
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(1,2,2)
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
+# place a text box in upper left in axes coords
+a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
+ verticalalignment='top', bbox=props)
+imgplot = plt.imshow(ndf_inp_linear, cmap="gray")
+plt.title('{}'.format('Linear diffusion inpainting results'))
+#%%
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+print ("_Inpainting using nonlinear diffusion (2D)_")
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+
+## plot
+fig = plt.figure(3)
+plt.suptitle('Performance of nonlinear diffusion inpainting using the CPU')
+a=fig.add_subplot(1,2,1)
+a.set_title('Missing data sinogram')
+imgplot = plt.imshow(sino_cut_new,cmap="gray")
+
+# set parameters
+pars = {'algorithm' : NDF_INP, \
+ 'input' : sino_cut_new,\
+ 'maskData' : mask,\
+ 'regularisation_parameter':80,\
+ 'edge_parameter':0.00009,\
+ 'number_of_iterations' :1500 ,\
+ 'time_marching_parameter':0.000008,\
+ 'penalty_type':1
+ }
+
+start_time = timeit.default_timer()
+ndf_inp_nonlinear = NDF_INP(pars['input'],
+ pars['maskData'],
+ pars['regularisation_parameter'],
+ pars['edge_parameter'],
+ pars['number_of_iterations'],
+ pars['time_marching_parameter'],
+ pars['penalty_type'])
+
+rms = rmse(sino_full, ndf_inp_nonlinear)
+pars['rmse'] = rms
+
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(1,2,2)
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
+# place a text box in upper left in axes coords
+a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
+ verticalalignment='top', bbox=props)
+imgplot = plt.imshow(ndf_inp_nonlinear, cmap="gray")
+plt.title('{}'.format('Nonlinear diffusion inpainting results'))
+#%%
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+print ("Inpainting using nonlocal vertical marching")
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+
+## plot
+fig = plt.figure(4)
+plt.suptitle('Performance of NVM inpainting using the CPU')
+a=fig.add_subplot(1,2,1)
+a.set_title('Missing data sinogram')
+imgplot = plt.imshow(sino_cut,cmap="gray")
+
+# set parameters
+pars = {'algorithm' : NVM_INP, \
+ 'input' : sino_cut_new,\
+ 'maskData' : mask,\
+ 'SW_increment': 1,\
+ 'number_of_iterations' : 150
+ }
+
+start_time = timeit.default_timer()
+(nvm_inp, mask_upd) = NVM_INP(pars['input'],
+ pars['maskData'],
+ pars['SW_increment'],
+ pars['number_of_iterations'])
+
+rms = rmse(sino_full, nvm_inp)
+pars['rmse'] = rms
+
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(1,2,2)
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
+# place a text box in upper left in axes coords
+a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
+ verticalalignment='top', bbox=props)
+imgplot = plt.imshow(nvm_inp, cmap="gray")
+plt.title('{}'.format('Nonlocal Vertical Marching inpainting results'))
+#%%
diff --git a/Wrappers/Python/demos/demo_cpu_regularisers.py b/Wrappers/Python/demos/demo_cpu_regularisers.py
index 3567f91..986e3e9 100644
--- a/Wrappers/Python/demos/demo_cpu_regularisers.py
+++ b/Wrappers/Python/demos/demo_cpu_regularisers.py
@@ -44,13 +44,31 @@ u0 = Im + np.random.normal(loc = 0 ,
u_ref = Im + np.random.normal(loc = 0 ,
scale = 0.01 * Im ,
size = np.shape(Im))
-
+(N,M) = np.shape(u0)
# map the u0 u0->u0>0
# f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
u0 = u0.astype('float32')
u_ref = u_ref.astype('float32')
-
+# change dims to check that modules work with non-squared images
+"""
+M = M-100
+u_ref2 = np.zeros([N,M],dtype='float32')
+u_ref2[:,0:M] = u_ref[:,0:M]
+u_ref = u_ref2
+del u_ref2
+
+u02 = np.zeros([N,M],dtype='float32')
+u02[:,0:M] = u0[:,0:M]
+u0 = u02
+del u02
+
+Im2 = np.zeros([N,M],dtype='float32')
+Im2[:,0:M] = Im[:,0:M]
+Im = Im2
+del Im2
+"""
+#%%
print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
print ("_______________ROF-TV (2D)_________________")
print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
@@ -288,7 +306,6 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
imgplot = plt.imshow(fgp_dtv_cpu, cmap="gray")
plt.title('{}'.format('CPU results'))
-
print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
print ("__________Total nuclear Variation__________")
print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
@@ -301,9 +318,8 @@ a.set_title('Noisy Image')
imgplot = plt.imshow(u0,cmap="gray")
channelsNo = 5
-N = 512
-noisyVol = np.zeros((channelsNo,N,N),dtype='float32')
-idealVol = np.zeros((channelsNo,N,N),dtype='float32')
+noisyVol = np.zeros((channelsNo,N,M),dtype='float32')
+idealVol = np.zeros((channelsNo,N,M),dtype='float32')
for i in range (channelsNo):
noisyVol[i,:,:] = Im + np.random.normal(loc = 0 , scale = perc * Im , size = np.shape(Im))
@@ -344,25 +360,19 @@ plt.title('{}'.format('CPU results'))
# Uncomment to test 3D regularisation performance
#%%
"""
-N = 512
slices = 20
-
-filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif")
-Im = plt.imread(filename)
-Im = np.asarray(Im, dtype='float32')
-
-Im = Im/255
perc = 0.05
-noisyVol = np.zeros((slices,N,N),dtype='float32')
-noisyRef = np.zeros((slices,N,N),dtype='float32')
-idealVol = np.zeros((slices,N,N),dtype='float32')
+noisyVol = np.zeros((slices,N,M),dtype='float32')
+noisyRef = np.zeros((slices,N,M),dtype='float32')
+idealVol = np.zeros((slices,N,M),dtype='float32')
for i in range (slices):
noisyVol[i,:,:] = Im + np.random.normal(loc = 0 , scale = perc * Im , size = np.shape(Im))
noisyRef[i,:,:] = Im + np.random.normal(loc = 0 , scale = 0.01 * Im , size = np.shape(Im))
idealVol[i,:,:] = Im
+
print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
print ("_______________ROF-TV (3D)_________________")
print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
@@ -403,6 +413,7 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
imgplot = plt.imshow(rof_cpu3D[10,:,:], cmap="gray")
plt.title('{}'.format('Recovered volume on the CPU using ROF-TV'))
+
print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
print ("_______________FGP-TV (3D)__________________")
print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
diff --git a/Wrappers/Python/demos/demo_gpu_regularisers.py b/Wrappers/Python/demos/demo_gpu_regularisers.py
index b873700..f3ed50c 100644
--- a/Wrappers/Python/demos/demo_gpu_regularisers.py
+++ b/Wrappers/Python/demos/demo_gpu_regularisers.py
@@ -44,10 +44,28 @@ u0 = Im + np.random.normal(loc = 0 ,
u_ref = Im + np.random.normal(loc = 0 ,
scale = 0.01 * Im ,
size = np.shape(Im))
+(N,M) = np.shape(u0)
# map the u0 u0->u0>0
# f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
u0 = u0.astype('float32')
u_ref = u_ref.astype('float32')
+"""
+M = M-100
+u_ref2 = np.zeros([N,M],dtype='float32')
+u_ref2[:,0:M] = u_ref[:,0:M]
+u_ref = u_ref2
+del u_ref2
+
+u02 = np.zeros([N,M],dtype='float32')
+u02[:,0:M] = u0[:,0:M]
+u0 = u02
+del u02
+
+Im2 = np.zeros([N,M],dtype='float32')
+Im2[:,0:M] = Im[:,0:M]
+Im = Im2
+del Im2
+"""
print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
print ("____________ROF-TV regulariser_____________")
diff --git a/Wrappers/Python/setup-regularisers.py.in b/Wrappers/Python/setup-regularisers.py.in
index b900efe..f55c6fe 100644
--- a/Wrappers/Python/setup-regularisers.py.in
+++ b/Wrappers/Python/setup-regularisers.py.in
@@ -34,6 +34,7 @@ extra_libraries = ['cilreg']
extra_include_dirs += [os.path.join(".." , ".." , "Core"),
os.path.join(".." , ".." , "Core", "regularisers_CPU"),
+ os.path.join(".." , ".." , "Core", "inpainters_CPU"),
os.path.join(".." , ".." , "Core", "regularisers_GPU" , "TV_FGP" ) ,
os.path.join(".." , ".." , "Core", "regularisers_GPU" , "TV_ROF" ) ,
os.path.join(".." , ".." , "Core", "regularisers_GPU" , "TV_SB" ) ,
@@ -52,7 +53,7 @@ setup(
description='CCPi Core Imaging Library - Image regularisers',
version=cil_version,
cmdclass = {'build_ext': build_ext},
- ext_modules = [Extension("ccpi.filters.cpu_regularisers_cython",
+ ext_modules = [Extension("ccpi.filters.cpu_regularisers",
sources=[os.path.join("." , "src", "cpu_regularisers.pyx" ) ],
include_dirs=extra_include_dirs,
library_dirs=extra_library_dirs,
diff --git a/Wrappers/Python/src/cpu_regularisers.pyx b/Wrappers/Python/src/cpu_regularisers.pyx
index 7ed8fa1..c934f1d 100644
--- a/Wrappers/Python/src/cpu_regularisers.pyx
+++ b/Wrappers/Python/src/cpu_regularisers.pyx
@@ -25,6 +25,8 @@ cdef extern float Diffusion_CPU_main(float *Input, float *Output, float lambdaPa
cdef extern float TNV_CPU_main(float *Input, float *u, float lambdaPar, int maxIter, float tol, int dimX, int dimY, int dimZ);
cdef extern float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float lambdaPar, int iterationsNumb, float epsil, float eta, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ);
+cdef extern float Diffusion_Inpaint_CPU_main(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, int dimX, int dimY, int dimZ);
+cdef extern float NonlocalMarching_Inpaint_main(float *Input, unsigned char *M, float *Output, unsigned char *M_upd, int SW_increment, int iterationsNumb, int trigger, int dimX, int dimY, int dimZ);
#****************************************************************#
#********************** Total-variation ROF *********************#
#****************************************************************#
@@ -46,7 +48,7 @@ def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
np.zeros([dims[0],dims[1]], dtype='float32')
# Run ROF iterations for 2D data
- TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, iterationsNumb, marching_step_parameter, dims[0], dims[1], 1)
+ TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, iterationsNumb, marching_step_parameter, dims[1], dims[0], 1)
return outputData
@@ -99,7 +101,7 @@ def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
methodTV,
nonneg,
printM,
- dims[0], dims[1], 1)
+ dims[1],dims[0],1)
return outputData
@@ -158,7 +160,7 @@ def TV_SB_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
tolerance_param,
methodTV,
printM,
- dims[0], dims[1], 1)
+ dims[1],dims[0],1)
return outputData
@@ -219,7 +221,7 @@ def dTV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
methodTV,
nonneg,
printM,
- dims[0], dims[1], 1)
+ dims[1], dims[0], 1)
return outputData
@@ -298,7 +300,7 @@ def NDF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
np.zeros([dims[0],dims[1]], dtype='float32')
# Run Nonlinear Diffusion iterations for 2D data
- Diffusion_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[0], dims[1], 1)
+ Diffusion_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[1], dims[0], 1)
return outputData
def NDF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
@@ -319,3 +321,81 @@ def NDF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
Diffusion_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[2], dims[1], dims[0])
return outputData
+
+#*********************Inpainting WITH****************************#
+#***************Nonlinear (Isotropic) Diffusion******************#
+#****************************************************************#
+def NDF_INPAINT_CPU(inputData, maskData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type):
+ if inputData.ndim == 2:
+ return NDF_INP_2D(inputData, maskData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type)
+ elif inputData.ndim == 3:
+ return NDF_INP_3D(inputData, maskData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type)
+
+def NDF_INP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
+ np.ndarray[np.uint8_t, ndim=2, mode="c"] maskData,
+ float regularisation_parameter,
+ float edge_parameter,
+ int iterationsNumb,
+ float time_marching_parameter,
+ int penalty_type):
+
+ cdef long dims[2]
+ dims[0] = inputData.shape[0]
+ dims[1] = inputData.shape[1]
+
+
+ cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
+ np.zeros([dims[0],dims[1]], dtype='float32')
+
+ # Run Inpaiting by Diffusion iterations for 2D data
+ Diffusion_Inpaint_CPU_main(&inputData[0,0], &maskData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[1], dims[0], 1)
+ return outputData
+
+def NDF_INP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
+ np.ndarray[np.uint8_t, ndim=3, mode="c"] maskData,
+ float regularisation_parameter,
+ float edge_parameter,
+ int iterationsNumb,
+ float time_marching_parameter,
+ int penalty_type):
+ cdef long dims[3]
+ dims[0] = inputData.shape[0]
+ dims[1] = inputData.shape[1]
+ dims[2] = inputData.shape[2]
+
+ cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
+ np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
+
+ # Run Inpaiting by Diffusion iterations for 3D data
+ Diffusion_Inpaint_CPU_main(&inputData[0,0,0], &maskData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[2], dims[1], dims[0])
+
+ return outputData
+#*********************Inpainting WITH****************************#
+#***************Nonlocal Vertical Marching method****************#
+#****************************************************************#
+def NVM_INPAINT_CPU(inputData, maskData, SW_increment, iterationsNumb):
+ if inputData.ndim == 2:
+ return NVM_INP_2D(inputData, maskData, SW_increment, iterationsNumb)
+ elif inputData.ndim == 3:
+ return
+
+def NVM_INP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
+ np.ndarray[np.uint8_t, ndim=2, mode="c"] maskData,
+ int SW_increment,
+ int iterationsNumb):
+ cdef long dims[2]
+ dims[0] = inputData.shape[0]
+ dims[1] = inputData.shape[1]
+
+ cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
+ np.zeros([dims[0],dims[1]], dtype='float32')
+
+ cdef np.ndarray[np.uint8_t, ndim=2, mode="c"] maskData_upd = \
+ np.zeros([dims[0],dims[1]], dtype='uint8')
+
+ # Run Inpaiting by Nonlocal vertical marching method for 2D data
+ NonlocalMarching_Inpaint_main(&inputData[0,0], &maskData[0,0], &outputData[0,0],
+ &maskData_upd[0,0],
+ SW_increment, iterationsNumb, 1, dims[1], dims[0], 1)
+
+ return (outputData, maskData_upd)
diff --git a/Wrappers/Python/src/gpu_regularisers.pyx b/Wrappers/Python/src/gpu_regularisers.pyx
index b0775054..7eab5d5 100644
--- a/Wrappers/Python/src/gpu_regularisers.pyx
+++ b/Wrappers/Python/src/gpu_regularisers.pyx
@@ -157,7 +157,7 @@ def ROFTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
regularisation_parameter,
iterations ,
time_marching_parameter,
- dims[0], dims[1], 1);
+ dims[1], dims[0], 1);
return outputData
@@ -210,7 +210,7 @@ def FGPTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
methodTV,
nonneg,
printM,
- dims[0], dims[1], 1);
+ dims[1], dims[0], 1);
return outputData
@@ -266,7 +266,7 @@ def SBTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
tolerance_param,
methodTV,
printM,
- dims[0], dims[1], 1);
+ dims[1], dims[0], 1);
return outputData
@@ -325,7 +325,7 @@ def FGPdTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
methodTV,
nonneg,
printM,
- dims[0], dims[1], 1);
+ dims[1], dims[0], 1);
return outputData
@@ -381,7 +381,7 @@ def NDF_GPU_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
# Run Nonlinear Diffusion iterations for 2D data
# Running CUDA code here
- NonlDiff_GPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[0], dims[1], 1)
+ NonlDiff_GPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[1], dims[0], 1)
return outputData
def NDF_GPU_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
generated by cgit v1.2.3 (git 2.25.1) at 2025-01-09 19:30:36 +0000