From b9fafd363d1d181a4a8b42ea4038924097207913 Mon Sep 17 00:00:00 2001
From: Daniil Kazantsev <dkazanc@hotmail.com>
Date: Mon, 9 Apr 2018 13:41:06 +0100
Subject: major renaming and new 3D demos for Matlab
---
Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m | 44 ++
Wrappers/Matlab/demos/demoMatlab_denoise.m | 8 +-
Wrappers/Matlab/mex_compile/compileCPU_mex.m | 8 +-
Wrappers/Matlab/mex_compile/compileGPU_mex.m | 8 +-
.../Matlab/mex_compile/regularisers_CPU/FGP_TV.c | 96 ++++
.../Matlab/mex_compile/regularisers_CPU/FGP_TV.c~ | 91 ++++
.../Matlab/mex_compile/regularisers_CPU/ROF_TV.c | 73 +++
.../mex_compile/regularisers_GPU/FGP_TV_GPU.cpp | 95 ++++
.../mex_compile/regularisers_GPU/ROF_TV_GPU.cpp | 73 +++
.../Matlab/mex_compile/regularizers_CPU/FGP_TV.c | 96 ----
.../Matlab/mex_compile/regularizers_CPU/ROF_TV.c | 73 ---
.../mex_compile/regularizers_GPU/FGP_TV_GPU.cpp | 95 ----
.../mex_compile/regularizers_GPU/ROF_TV_GPU.cpp | 73 ---
Wrappers/Python/CMakeLists.txt | 20 +-
Wrappers/Python/ccpi/filters/regularisers.py | 44 ++
Wrappers/Python/ccpi/filters/regularizers.py | 44 --
Wrappers/Python/conda-recipe/bld.bat | 4 +-
Wrappers/Python/conda-recipe/build.sh | 4 +-
Wrappers/Python/conda-recipe/meta.yaml | 6 +-
Wrappers/Python/demo/test_cpu_regularisers.py | 323 ++++++++++++++
Wrappers/Python/demo/test_cpu_regularizers.py | 488 ---------------------
Wrappers/Python/fista-recipe/bld.bat | 13 -
Wrappers/Python/fista-recipe/build.sh | 10 -
Wrappers/Python/fista-recipe/meta.yaml | 29 --
Wrappers/Python/setup-fista.py | 27 --
Wrappers/Python/setup-regularisers.py.in | 65 +++
Wrappers/Python/setup-regularizers.py.in | 65 ---
Wrappers/Python/src/cpu_regularisers.pyx | 126 ++++++
Wrappers/Python/src/cpu_regularizers.pyx | 126 ------
Wrappers/Python/src/gpu_regularisers.pyx | 170 +++++++
Wrappers/Python/src/gpu_regularizers.pyx | 170 -------
Wrappers/Python/test/test_cpu_regularisers.py | 442 +++++++++++++++++++
Wrappers/Python/test/test_cpu_regularizers.py | 442 -------------------
.../Python/test/test_cpu_vs_gpu_regularisers.py | 219 +++++++++
.../Python/test/test_cpu_vs_gpu_regularizers.py | 219 ---------
Wrappers/Python/test/test_gpu_regularisers.py | 143 ++++++
Wrappers/Python/test/test_gpu_regularizers.py | 239 ----------
Wrappers/Python/test/test_regularisers_3d.py | 425 ++++++++++++++++++
Wrappers/Python/test/test_regularizers.py | 395 -----------------
Wrappers/Python/test/test_regularizers_3d.py | 425 ------------------
40 files changed, 2458 insertions(+), 3058 deletions(-)
create mode 100644 Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m
create mode 100644 Wrappers/Matlab/mex_compile/regularisers_CPU/FGP_TV.c
create mode 100644 Wrappers/Matlab/mex_compile/regularisers_CPU/FGP_TV.c~
create mode 100644 Wrappers/Matlab/mex_compile/regularisers_CPU/ROF_TV.c
create mode 100644 Wrappers/Matlab/mex_compile/regularisers_GPU/FGP_TV_GPU.cpp
create mode 100644 Wrappers/Matlab/mex_compile/regularisers_GPU/ROF_TV_GPU.cpp
delete mode 100644 Wrappers/Matlab/mex_compile/regularizers_CPU/FGP_TV.c
delete mode 100644 Wrappers/Matlab/mex_compile/regularizers_CPU/ROF_TV.c
delete mode 100644 Wrappers/Matlab/mex_compile/regularizers_GPU/FGP_TV_GPU.cpp
delete mode 100644 Wrappers/Matlab/mex_compile/regularizers_GPU/ROF_TV_GPU.cpp
create mode 100644 Wrappers/Python/ccpi/filters/regularisers.py
delete mode 100644 Wrappers/Python/ccpi/filters/regularizers.py
create mode 100644 Wrappers/Python/demo/test_cpu_regularisers.py
delete mode 100644 Wrappers/Python/demo/test_cpu_regularizers.py
delete mode 100644 Wrappers/Python/fista-recipe/bld.bat
delete mode 100644 Wrappers/Python/fista-recipe/build.sh
delete mode 100644 Wrappers/Python/fista-recipe/meta.yaml
delete mode 100644 Wrappers/Python/setup-fista.py
create mode 100644 Wrappers/Python/setup-regularisers.py.in
delete mode 100644 Wrappers/Python/setup-regularizers.py.in
create mode 100644 Wrappers/Python/src/cpu_regularisers.pyx
delete mode 100644 Wrappers/Python/src/cpu_regularizers.pyx
create mode 100644 Wrappers/Python/src/gpu_regularisers.pyx
delete mode 100644 Wrappers/Python/src/gpu_regularizers.pyx
create mode 100644 Wrappers/Python/test/test_cpu_regularisers.py
delete mode 100644 Wrappers/Python/test/test_cpu_regularizers.py
create mode 100644 Wrappers/Python/test/test_cpu_vs_gpu_regularisers.py
delete mode 100644 Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py
create mode 100644 Wrappers/Python/test/test_gpu_regularisers.py
delete mode 100644 Wrappers/Python/test/test_gpu_regularizers.py
create mode 100644 Wrappers/Python/test/test_regularisers_3d.py
delete mode 100644 Wrappers/Python/test/test_regularizers.py
delete mode 100644 Wrappers/Python/test/test_regularizers_3d.py
(limited to 'Wrappers')
diff --git a/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m b/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m
new file mode 100644
index 0000000..f5c3ad1
--- /dev/null
+++ b/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m
@@ -0,0 +1,44 @@
+% Volume (3D) denoising demo using CCPi-RGL
+
+addpath('../mex_compile/installed');
+addpath('../../../data/');
+
+N = 256;
+slices = 30;
+vol3D = zeros(N,N,slices, 'single');
+Im = double(imread('lena_gray_256.tif'))/255; % loading image
+for i = 1:slices
+vol3D(:,:,i) = Im + .05*randn(size(Im));
+end
+vol3D(vol3D < 0) = 0;
+figure; imshow(vol3D(:,:,15), [0 1]); title('Noisy image');
+
+%%
+fprintf('Denoise using ROF-TV model (CPU) \n');
+lambda_rof = 0.03; % regularisation parameter
+tau_rof = 0.0025; % time-marching constant
+iter_rof = 1000; % number of ROF iterations
+tic; u_rof = ROF_TV(single(vol3D), lambda_rof, iter_rof, tau_rof); toc;
+figure; imshow(u_rof(:,:,15), [0 1]); title('ROF-TV denoised volume (CPU)');
+%%
+% fprintf('Denoise using ROF-TV model (GPU) \n');
+% lambda_rof = 0.03; % regularisation parameter
+% tau_rof = 0.0025; % time-marching constant
+% iter_rof = 1000; % number of ROF iterations
+% tic; u_rofG = ROF_TV_GPU(single(vol3D), lambda_rof, iter_rof, tau_rof); toc;
+% figure; imshow(u_rofG(:,:,15), [0 1]); title('ROF-TV denoised volume (GPU)');
+%%
+fprintf('Denoise using FGP-TV model (CPU) \n');
+lambda_fgp = 0.03; % regularisation parameter
+iter_fgp = 500; % number of FGP iterations
+epsil_tol = 1.0e-05; % tolerance
+tic; u_fgp = FGP_TV(single(vol3D), lambda_fgp, iter_fgp, epsil_tol); toc;
+figure; imshow(u_fgp(:,:,15), [0 1]); title('FGP-TV denoised volume (CPU)');
+%%
+% fprintf('Denoise using FGP-TV model (GPU) \n');
+% lambda_fgp = 0.03; % regularisation parameter
+% iter_fgp = 500; % number of FGP iterations
+% epsil_tol = 1.0e-05; % tolerance
+% tic; u_fgpG = FGP_TV_GPU(single(vol3D), lambda_fgp, iter_fgp, epsil_tol); toc;
+% figure; imshow(u_fgpG(:,:,15), [0 1]); title('FGP-TV denoised volume (GPU)');
+%%
\ No newline at end of file
diff --git a/Wrappers/Matlab/demos/demoMatlab_denoise.m b/Wrappers/Matlab/demos/demoMatlab_denoise.m
index 7258e5e..ab4e95d 100644
--- a/Wrappers/Matlab/demos/demoMatlab_denoise.m
+++ b/Wrappers/Matlab/demos/demoMatlab_denoise.m
@@ -9,28 +9,28 @@ figure; imshow(u0, [0 1]); title('Noisy image');
%%
fprintf('Denoise using ROF-TV model (CPU) \n');
-lambda_rof = 0.03; % regularization parameter
+lambda_rof = 0.03; % regularisation parameter
tau_rof = 0.0025; % time-marching constant
iter_rof = 2000; % number of ROF iterations
tic; u_rof = ROF_TV(single(u0), lambda_rof, iter_rof, tau_rof); toc;
figure; imshow(u_rof, [0 1]); title('ROF-TV denoised image (CPU)');
%%
% fprintf('Denoise using ROF-TV model (GPU) \n');
-% lambda_rof = 0.03; % regularization parameter
+% lambda_rof = 0.03; % regularisation parameter
% tau_rof = 0.0025; % time-marching constant
% iter_rof = 2000; % number of ROF iterations
% tic; u_rofG = ROF_TV_GPU(single(u0), lambda_rof, iter_rof, tau_rof); toc;
% figure; imshow(u_rofG, [0 1]); title('ROF-TV denoised image (GPU)');
%%
fprintf('Denoise using FGP-TV model (CPU) \n');
-lambda_fgp = 0.03; % regularization parameter
+lambda_fgp = 0.03; % regularisation parameter
iter_fgp = 1000; % number of FGP iterations
epsil_tol = 1.0e-05; % tolerance
tic; u_fgp = FGP_TV(single(u0), lambda_fgp, iter_fgp, epsil_tol); toc;
figure; imshow(u_fgp, [0 1]); title('FGP-TV denoised image (CPU)');
%%
% fprintf('Denoise using FGP-TV model (GPU) \n');
-% lambda_fgp = 0.03; % regularization parameter
+% lambda_fgp = 0.03; % regularisation parameter
% iter_fgp = 1000; % number of FGP iterations
% epsil_tol = 1.0e-05; % tolerance
% tic; u_fgpG = FGP_TV_GPU(single(u0), lambda_fgp, iter_fgp, epsil_tol); toc;
diff --git a/Wrappers/Matlab/mex_compile/compileCPU_mex.m b/Wrappers/Matlab/mex_compile/compileCPU_mex.m
index fcee53a..8da81ad 100644
--- a/Wrappers/Matlab/mex_compile/compileCPU_mex.m
+++ b/Wrappers/Matlab/mex_compile/compileCPU_mex.m
@@ -1,10 +1,10 @@
% execute this mex file in Matlab once
-copyfile ../../../Core/regularizers_CPU/ regularizers_CPU/
-copyfile ../../../Core/CCPiDefines.h regularizers_CPU/
+copyfile ../../../Core/regularisers_CPU/ regularisers_CPU/
+copyfile ../../../Core/CCPiDefines.h regularisers_CPU/
-cd regularizers_CPU/
+cd regularisers_CPU/
-fprintf('%s \n', 'Compiling CPU regularizers...');
+fprintf('%s \n', 'Compiling CPU regularisers...');
mex ROF_TV.c ROF_TV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
movefile ROF_TV.mex* ../installed/
diff --git a/Wrappers/Matlab/mex_compile/compileGPU_mex.m b/Wrappers/Matlab/mex_compile/compileGPU_mex.m
index df29a3e..45236fa 100644
--- a/Wrappers/Matlab/mex_compile/compileGPU_mex.m
+++ b/Wrappers/Matlab/mex_compile/compileGPU_mex.m
@@ -9,12 +9,12 @@
% tested on Ubuntu 16.04/MATLAB 2016b
-copyfile ../../../Core/regularizers_GPU/ regularizers_GPU/
-copyfile ../../../Core/CCPiDefines.h regularizers_GPU/
+copyfile ../../../Core/regularisers_GPU/ regularisers_GPU/
+copyfile ../../../Core/CCPiDefines.h regularisers_GPU/
-cd regularizers_GPU/
+cd regularisers_GPU/
-fprintf('%s \n', 'Compiling GPU regularizers (CUDA)...');
+fprintf('%s \n', 'Compiling GPU regularisers (CUDA)...');
!/usr/local/cuda/bin/nvcc -O0 -c TV_ROF_GPU_core.cu -Xcompiler -fPIC -I~/SOFT/MATLAB9/extern/include/
mex -g -I/usr/local/cuda-7.5/include -L/usr/local/cuda-7.5/lib64 -lcudart -lcufft -lmwgpu ROF_TV_GPU.cpp TV_ROF_GPU_core.o
movefile ROF_TV_GPU.mex* ../installed/
diff --git a/Wrappers/Matlab/mex_compile/regularisers_CPU/FGP_TV.c b/Wrappers/Matlab/mex_compile/regularisers_CPU/FGP_TV.c
new file mode 100644
index 0000000..ba06cc7
--- /dev/null
+++ b/Wrappers/Matlab/mex_compile/regularisers_CPU/FGP_TV.c
@@ -0,0 +1,96 @@
+/*
+ * 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 "FGP_TV_core.h"
+
+/* C-OMP implementation of FGP-TV [1] denoising/regularization model (2D/3D case)
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume
+ * 2. lambdaPar - regularization parameter
+ * 3. Number of iterations
+ * 4. eplsilon: tolerance constant
+ * 5. TV-type: methodTV - 'iso' (0) or 'l1' (1)
+ * 6. nonneg: 'nonnegativity (0 is OFF by default)
+ * 7. print information: 0 (off) or 1 (on)
+ *
+ * Output:
+ * [1] Filtered/regularized image
+ *
+ * This function is based on the Matlab's code and paper by
+ * [1] Amir Beck and Marc Teboulle, "Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems"
+ */
+
+
+void mexFunction(
+ int nlhs, mxArray *plhs[],
+ int nrhs, const mxArray *prhs[])
+
+{
+ int number_of_dims, iter, dimX, dimY, dimZ, methTV, printswitch, nonneg;
+ const int *dim_array;
+ float *Input, *Output=NULL, lambda, epsil;
+
+ number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+ dim_array = mxGetDimensions(prhs[0]);
+
+ /*Handling Matlab input data*/
+ if ((nrhs < 2) || (nrhs > 7)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D/3D), Regularization parameter. The full list of parameters: Image(2D/3D), Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1'), nonnegativity switch, print switch");
+
+ Input = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */
+ lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */
+ iter = 300; /* default iterations number */
+ epsil = 0.0001; /* default tolerance constant */
+ methTV = 0; /* default isotropic TV penalty */
+ nonneg = 0; /* default nonnegativity switch, off - 0 */
+ printswitch = 0; /*default print is switched, off - 0 */
+
+ if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+
+ if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7)) iter = (int) mxGetScalar(prhs[2]); /* iterations number */
+ if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7)) epsil = (float) mxGetScalar(prhs[3]); /* tolerance constant */
+ if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7)) {
+ char *penalty_type;
+ penalty_type = mxArrayToString(prhs[4]); /* choosing TV penalty: 'iso' or 'l1', 'iso' is the default */
+ if ((strcmp(penalty_type, "l1") != 0) && (strcmp(penalty_type, "iso") != 0)) mexErrMsgTxt("Choose TV type: 'iso' or 'l1',");
+ if (strcmp(penalty_type, "l1") == 0) methTV = 1; /* enable 'l1' penalty */
+ mxFree(penalty_type);
+ }
+ if ((nrhs == 6) || (nrhs == 7)) {
+ nonneg = (int) mxGetScalar(prhs[5]);
+ if ((nonneg != 0) && (nonneg != 1)) mexErrMsgTxt("Nonnegativity constraint can be enabled by choosing 1 or off - 0");
+ }
+ if (nrhs == 7) {
+ printswitch = (int) mxGetScalar(prhs[6]);
+ if ((printswitch != 0) && (printswitch != 1)) mexErrMsgTxt("Print can be enabled by choosing 1 or off - 0");
+ }
+
+ /*Handling Matlab output data*/
+ dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+
+ if (number_of_dims == 2) {
+ dimZ = 1; /*2D case*/
+ Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+ }
+ if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+
+ /* running the function */
+ TV_FGP_CPU_main(Input, Output, lambda, iter, epsil, methTV, nonneg, printswitch, dimX, dimY, dimZ);
+}
\ No newline at end of file
diff --git a/Wrappers/Matlab/mex_compile/regularisers_CPU/FGP_TV.c~ b/Wrappers/Matlab/mex_compile/regularisers_CPU/FGP_TV.c~
new file mode 100644
index 0000000..30d61cd
--- /dev/null
+++ b/Wrappers/Matlab/mex_compile/regularisers_CPU/FGP_TV.c~
@@ -0,0 +1,91 @@
+/*
+ * 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 "FGP_TV_core.h"
+
+/* C-OMP implementation of FGP-TV [1] denoising/regularization model (2D/3D case)
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume
+ * 2. lambdaPar - regularization parameter
+ * 3. Number of iterations
+ * 4. eplsilon: tolerance constant
+ * 5. TV-type: methodTV - 'iso' (0) or 'l1' (1)
+ * 6. nonneg: 'nonnegativity (0 is OFF by default)
+ * 7. print information: 0 (off) or 1 (on)
+ *
+ * Output:
+ * [1] Filtered/regularized image
+ *
+ * This function is based on the Matlab's code and paper by
+ * [1] Amir Beck and Marc Teboulle, "Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems"
+ */
+
+
+void mexFunction(
+ int nlhs, mxArray *plhs[],
+ int nrhs, const mxArray *prhs[])
+
+{
+ int number_of_dims, iter, dimX, dimY, dimZ, methTV, printswitch;
+ const int *dim_array;
+ float *Input, *Output, lambda, epsil;
+
+ number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+ dim_array = mxGetDimensions(prhs[0]);
+
+ /*Handling Matlab input data*/
+ if ((nrhs < 2) || (nrhs > 6)) mexErrMsgTxt("At least 2 parameters is required: Image(2D/3D), Regularization parameter. The full list of parameters: Image(2D/3D), Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1'), print switch");
+
+ Input = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */
+ lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */
+ iter = 300; /* default iterations number */
+ epsil = 0.0001; /* default tolerance constant */
+ methTV = 0; /* default isotropic TV penalty */
+ printswitch = 0; /*default print is switched off - 0 */
+
+ if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+
+ if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5) || (nrhs == 6)) iter = (int) mxGetScalar(prhs[2]); /* iterations number */
+ if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6)) epsil = (float) mxGetScalar(prhs[3]); /* tolerance constant */
+ if ((nrhs == 5) || (nrhs == 6)) {
+ char *penalty_type;
+ penalty_type = mxArrayToString(prhs[4]); /* choosing TV penalty: 'iso' or 'l1', 'iso' is the default */
+ if ((strcmp(penalty_type, "l1") != 0) && (strcmp(penalty_type, "iso") != 0)) mexErrMsgTxt("Choose TV type: 'iso' or 'l1',");
+ if (strcmp(penalty_type, "l1") == 0) methTV = 1; /* enable 'l1' penalty */
+ mxFree(penalty_type);
+ }
+ if (nrhs == 6) {
+ printswitch = (int) mxGetScalar(prhs[5]);
+ if ((printswitch != 0) || (printswitch != 1)) {mexErrMsgTxt("Print can be enabled by choosing 1 or off - 0"); }
+ }
+
+ /*Handling Matlab output data*/
+ dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+
+ if (number_of_dims == 2) {
+ dimZ = 1; /*2D case*/
+ Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+ }
+ if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+
+
+ TV_FGP_CPU_main(Input, Output, lambda, iter, epsil, methTV, nonneg, printswitch, dimX, dimY, dimZ)
+}
diff --git a/Wrappers/Matlab/mex_compile/regularisers_CPU/ROF_TV.c b/Wrappers/Matlab/mex_compile/regularisers_CPU/ROF_TV.c
new file mode 100644
index 0000000..6b9e1ea
--- /dev/null
+++ b/Wrappers/Matlab/mex_compile/regularisers_CPU/ROF_TV.c
@@ -0,0 +1,73 @@
+/*
+ * 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 "ROF_TV_core.h"
+
+/* ROF-TV denoising/regularization model [1] (2D/3D case)
+ * (MEX wrapper for MATLAB)
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume [REQUIRED]
+ * 2. lambda - regularization parameter [REQUIRED]
+ * 3. Number of iterations, for explicit scheme >= 150 is recommended [REQUIRED]
+ * 4. tau - marching step for explicit scheme, ~1 is recommended [REQUIRED]
+ *
+ * Output:
+ * [1] Regularized image/volume
+ *
+ * This function is based on the paper by
+ * [1] Rudin, Osher, Fatemi, "Nonlinear Total Variation based noise removal algorithms"
+ *
+ * D. Kazantsev, 2016-18
+ */
+
+void mexFunction(
+ int nlhs, mxArray *plhs[],
+ int nrhs, const mxArray *prhs[])
+
+{
+ int number_of_dims, iter_numb, dimX, dimY, dimZ;
+ const int *dim_array;
+ float *Input, *Output=NULL, lambda, tau;
+
+ dim_array = mxGetDimensions(prhs[0]);
+ number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+
+ /*Handling Matlab input data*/
+ Input = (float *) mxGetData(prhs[0]);
+ lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */
+ iter_numb = (int) mxGetScalar(prhs[2]); /* iterations number */
+ tau = (float) mxGetScalar(prhs[3]); /* marching step parameter */
+
+ if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+ if(nrhs != 4) mexErrMsgTxt("Four inputs reqired: Image(2D,3D), regularization parameter, iterations number, marching step constant");
+ /*Handling Matlab output data*/
+ 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 */
+ Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+ }
+ if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+
+ TV_ROF_CPU_main(Input, Output, lambda, iter_numb, tau, dimX, dimY, dimZ);
+}
\ No newline at end of file
diff --git a/Wrappers/Matlab/mex_compile/regularisers_GPU/FGP_TV_GPU.cpp b/Wrappers/Matlab/mex_compile/regularisers_GPU/FGP_TV_GPU.cpp
new file mode 100644
index 0000000..9ed9ae0
--- /dev/null
+++ b/Wrappers/Matlab/mex_compile/regularisers_GPU/FGP_TV_GPU.cpp
@@ -0,0 +1,95 @@
+/*
+ * 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 "TV_FGP_GPU_core.h"
+
+/* GPU (CUDA) implementation of FGP-TV [1] denoising/regularization model (2D/3D case)
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume
+ * 2. lambdaPar - regularization parameter
+ * 3. Number of iterations
+ * 4. eplsilon: tolerance constant
+ * 5. TV-type: methodTV - 'iso' (0) or 'l1' (1)
+ * 6. nonneg: 'nonnegativity (0 is OFF by default)
+ * 7. print information: 0 (off) or 1 (on)
+ *
+ * Output:
+ * [1] Filtered/regularized image
+ *
+ * This function is based on the Matlab's code and paper by
+ * [1] Amir Beck and Marc Teboulle, "Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems"
+ */
+
+void mexFunction(
+ int nlhs, mxArray *plhs[],
+ int nrhs, const mxArray *prhs[])
+
+{
+ int number_of_dims, iter, dimX, dimY, dimZ, methTV, printswitch, nonneg;
+ const int *dim_array;
+ float *Input, *Output=NULL, lambda, epsil;
+
+ number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+ dim_array = mxGetDimensions(prhs[0]);
+
+ /*Handling Matlab input data*/
+ if ((nrhs < 2) || (nrhs > 7)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D/3D), Regularization parameter. The full list of parameters: Image(2D/3D), Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1'), nonnegativity switch, print switch");
+
+ Input = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */
+ lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */
+ iter = 300; /* default iterations number */
+ epsil = 0.0001; /* default tolerance constant */
+ methTV = 0; /* default isotropic TV penalty */
+ nonneg = 0; /* default nonnegativity switch, off - 0 */
+ printswitch = 0; /*default print is switched, off - 0 */
+
+ if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+
+ if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7)) iter = (int) mxGetScalar(prhs[2]); /* iterations number */
+ if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7)) epsil = (float) mxGetScalar(prhs[3]); /* tolerance constant */
+ if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7)) {
+ char *penalty_type;
+ penalty_type = mxArrayToString(prhs[4]); /* choosing TV penalty: 'iso' or 'l1', 'iso' is the default */
+ if ((strcmp(penalty_type, "l1") != 0) && (strcmp(penalty_type, "iso") != 0)) mexErrMsgTxt("Choose TV type: 'iso' or 'l1',");
+ if (strcmp(penalty_type, "l1") == 0) methTV = 1; /* enable 'l1' penalty */
+ mxFree(penalty_type);
+ }
+ if ((nrhs == 6) || (nrhs == 7)) {
+ nonneg = (int) mxGetScalar(prhs[5]);
+ if ((nonneg != 0) && (nonneg != 1)) mexErrMsgTxt("Nonnegativity constraint can be enabled by choosing 1 or off - 0");
+ }
+ if (nrhs == 7) {
+ printswitch = (int) mxGetScalar(prhs[6]);
+ if ((printswitch != 0) && (printswitch != 1)) mexErrMsgTxt("Print can be enabled by choosing 1 or off - 0");
+ }
+
+ /*Handling Matlab output data*/
+ dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+
+ if (number_of_dims == 2) {
+ dimZ = 1; /*2D case*/
+ Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+ }
+ if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+
+ /* running the function */
+ TV_FGP_GPU_main(Input, Output, lambda, iter, epsil, methTV, nonneg, printswitch, dimX, dimY, dimZ);
+}
\ No newline at end of file
diff --git a/Wrappers/Matlab/mex_compile/regularisers_GPU/ROF_TV_GPU.cpp b/Wrappers/Matlab/mex_compile/regularisers_GPU/ROF_TV_GPU.cpp
new file mode 100644
index 0000000..7bbe3af
--- /dev/null
+++ b/Wrappers/Matlab/mex_compile/regularisers_GPU/ROF_TV_GPU.cpp
@@ -0,0 +1,73 @@
+/*
+ * 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 "TV_ROF_GPU_core.h"
+
+/* ROF-TV denoising/regularization model [1] (2D/3D case)
+ * (MEX wrapper for MATLAB)
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume [REQUIRED]
+ * 2. lambda - regularization parameter [REQUIRED]
+ * 3. Number of iterations, for explicit scheme >= 150 is recommended [REQUIRED]
+ * 4. tau - marching step for explicit scheme, ~1 is recommended [REQUIRED]
+ *
+ * Output:
+ * [1] Regularized image/volume
+ *
+ * This function is based on the paper by
+ * [1] Rudin, Osher, Fatemi, "Nonlinear Total Variation based noise removal algorithms"
+ *
+ * D. Kazantsev, 2016-18
+ */
+
+void mexFunction(
+ int nlhs, mxArray *plhs[],
+ int nrhs, const mxArray *prhs[])
+
+{
+ int number_of_dims, iter_numb, dimX, dimY, dimZ;
+ const int *dim_array;
+ float *Input, *Output=NULL, lambda, tau;
+
+ dim_array = mxGetDimensions(prhs[0]);
+ number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+
+ /*Handling Matlab input data*/
+ Input = (float *) mxGetData(prhs[0]);
+ lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */
+ iter_numb = (int) mxGetScalar(prhs[2]); /* iterations number */
+ tau = (float) mxGetScalar(prhs[3]); /* marching step parameter */
+
+ if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+ if(nrhs != 4) mexErrMsgTxt("Four inputs reqired: Image(2D,3D), regularization parameter, iterations number, marching step constant");
+ /*Handling Matlab output data*/
+ 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 */
+ Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+ }
+ if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+
+ TV_ROF_GPU_main(Input, Output, lambda, iter_numb, tau, dimX, dimY, dimZ);
+}
\ No newline at end of file
diff --git a/Wrappers/Matlab/mex_compile/regularizers_CPU/FGP_TV.c b/Wrappers/Matlab/mex_compile/regularizers_CPU/FGP_TV.c
deleted file mode 100644
index ba06cc7..0000000
--- a/Wrappers/Matlab/mex_compile/regularizers_CPU/FGP_TV.c
+++ /dev/null
@@ -1,96 +0,0 @@
-/*
- * 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 "FGP_TV_core.h"
-
-/* C-OMP implementation of FGP-TV [1] denoising/regularization model (2D/3D case)
- *
- * Input Parameters:
- * 1. Noisy image/volume
- * 2. lambdaPar - regularization parameter
- * 3. Number of iterations
- * 4. eplsilon: tolerance constant
- * 5. TV-type: methodTV - 'iso' (0) or 'l1' (1)
- * 6. nonneg: 'nonnegativity (0 is OFF by default)
- * 7. print information: 0 (off) or 1 (on)
- *
- * Output:
- * [1] Filtered/regularized image
- *
- * This function is based on the Matlab's code and paper by
- * [1] Amir Beck and Marc Teboulle, "Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems"
- */
-
-
-void mexFunction(
- int nlhs, mxArray *plhs[],
- int nrhs, const mxArray *prhs[])
-
-{
- int number_of_dims, iter, dimX, dimY, dimZ, methTV, printswitch, nonneg;
- const int *dim_array;
- float *Input, *Output=NULL, lambda, epsil;
-
- number_of_dims = mxGetNumberOfDimensions(prhs[0]);
- dim_array = mxGetDimensions(prhs[0]);
-
- /*Handling Matlab input data*/
- if ((nrhs < 2) || (nrhs > 7)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D/3D), Regularization parameter. The full list of parameters: Image(2D/3D), Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1'), nonnegativity switch, print switch");
-
- Input = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */
- lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */
- iter = 300; /* default iterations number */
- epsil = 0.0001; /* default tolerance constant */
- methTV = 0; /* default isotropic TV penalty */
- nonneg = 0; /* default nonnegativity switch, off - 0 */
- printswitch = 0; /*default print is switched, off - 0 */
-
- if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
-
- if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7)) iter = (int) mxGetScalar(prhs[2]); /* iterations number */
- if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7)) epsil = (float) mxGetScalar(prhs[3]); /* tolerance constant */
- if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7)) {
- char *penalty_type;
- penalty_type = mxArrayToString(prhs[4]); /* choosing TV penalty: 'iso' or 'l1', 'iso' is the default */
- if ((strcmp(penalty_type, "l1") != 0) && (strcmp(penalty_type, "iso") != 0)) mexErrMsgTxt("Choose TV type: 'iso' or 'l1',");
- if (strcmp(penalty_type, "l1") == 0) methTV = 1; /* enable 'l1' penalty */
- mxFree(penalty_type);
- }
- if ((nrhs == 6) || (nrhs == 7)) {
- nonneg = (int) mxGetScalar(prhs[5]);
- if ((nonneg != 0) && (nonneg != 1)) mexErrMsgTxt("Nonnegativity constraint can be enabled by choosing 1 or off - 0");
- }
- if (nrhs == 7) {
- printswitch = (int) mxGetScalar(prhs[6]);
- if ((printswitch != 0) && (printswitch != 1)) mexErrMsgTxt("Print can be enabled by choosing 1 or off - 0");
- }
-
- /*Handling Matlab output data*/
- dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
-
- if (number_of_dims == 2) {
- dimZ = 1; /*2D case*/
- Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
- }
- if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-
- /* running the function */
- TV_FGP_CPU_main(Input, Output, lambda, iter, epsil, methTV, nonneg, printswitch, dimX, dimY, dimZ);
-}
\ No newline at end of file
diff --git a/Wrappers/Matlab/mex_compile/regularizers_CPU/ROF_TV.c b/Wrappers/Matlab/mex_compile/regularizers_CPU/ROF_TV.c
deleted file mode 100644
index 6b9e1ea..0000000
--- a/Wrappers/Matlab/mex_compile/regularizers_CPU/ROF_TV.c
+++ /dev/null
@@ -1,73 +0,0 @@
-/*
- * 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 "ROF_TV_core.h"
-
-/* ROF-TV denoising/regularization model [1] (2D/3D case)
- * (MEX wrapper for MATLAB)
- *
- * Input Parameters:
- * 1. Noisy image/volume [REQUIRED]
- * 2. lambda - regularization parameter [REQUIRED]
- * 3. Number of iterations, for explicit scheme >= 150 is recommended [REQUIRED]
- * 4. tau - marching step for explicit scheme, ~1 is recommended [REQUIRED]
- *
- * Output:
- * [1] Regularized image/volume
- *
- * This function is based on the paper by
- * [1] Rudin, Osher, Fatemi, "Nonlinear Total Variation based noise removal algorithms"
- *
- * D. Kazantsev, 2016-18
- */
-
-void mexFunction(
- int nlhs, mxArray *plhs[],
- int nrhs, const mxArray *prhs[])
-
-{
- int number_of_dims, iter_numb, dimX, dimY, dimZ;
- const int *dim_array;
- float *Input, *Output=NULL, lambda, tau;
-
- dim_array = mxGetDimensions(prhs[0]);
- number_of_dims = mxGetNumberOfDimensions(prhs[0]);
-
- /*Handling Matlab input data*/
- Input = (float *) mxGetData(prhs[0]);
- lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */
- iter_numb = (int) mxGetScalar(prhs[2]); /* iterations number */
- tau = (float) mxGetScalar(prhs[3]); /* marching step parameter */
-
- if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
- if(nrhs != 4) mexErrMsgTxt("Four inputs reqired: Image(2D,3D), regularization parameter, iterations number, marching step constant");
- /*Handling Matlab output data*/
- 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 */
- Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
- }
- if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-
- TV_ROF_CPU_main(Input, Output, lambda, iter_numb, tau, dimX, dimY, dimZ);
-}
\ No newline at end of file
diff --git a/Wrappers/Matlab/mex_compile/regularizers_GPU/FGP_TV_GPU.cpp b/Wrappers/Matlab/mex_compile/regularizers_GPU/FGP_TV_GPU.cpp
deleted file mode 100644
index 9ed9ae0..0000000
--- a/Wrappers/Matlab/mex_compile/regularizers_GPU/FGP_TV_GPU.cpp
+++ /dev/null
@@ -1,95 +0,0 @@
-/*
- * 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 "TV_FGP_GPU_core.h"
-
-/* GPU (CUDA) implementation of FGP-TV [1] denoising/regularization model (2D/3D case)
- *
- * Input Parameters:
- * 1. Noisy image/volume
- * 2. lambdaPar - regularization parameter
- * 3. Number of iterations
- * 4. eplsilon: tolerance constant
- * 5. TV-type: methodTV - 'iso' (0) or 'l1' (1)
- * 6. nonneg: 'nonnegativity (0 is OFF by default)
- * 7. print information: 0 (off) or 1 (on)
- *
- * Output:
- * [1] Filtered/regularized image
- *
- * This function is based on the Matlab's code and paper by
- * [1] Amir Beck and Marc Teboulle, "Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems"
- */
-
-void mexFunction(
- int nlhs, mxArray *plhs[],
- int nrhs, const mxArray *prhs[])
-
-{
- int number_of_dims, iter, dimX, dimY, dimZ, methTV, printswitch, nonneg;
- const int *dim_array;
- float *Input, *Output=NULL, lambda, epsil;
-
- number_of_dims = mxGetNumberOfDimensions(prhs[0]);
- dim_array = mxGetDimensions(prhs[0]);
-
- /*Handling Matlab input data*/
- if ((nrhs < 2) || (nrhs > 7)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D/3D), Regularization parameter. The full list of parameters: Image(2D/3D), Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1'), nonnegativity switch, print switch");
-
- Input = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */
- lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */
- iter = 300; /* default iterations number */
- epsil = 0.0001; /* default tolerance constant */
- methTV = 0; /* default isotropic TV penalty */
- nonneg = 0; /* default nonnegativity switch, off - 0 */
- printswitch = 0; /*default print is switched, off - 0 */
-
- if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
-
- if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7)) iter = (int) mxGetScalar(prhs[2]); /* iterations number */
- if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7)) epsil = (float) mxGetScalar(prhs[3]); /* tolerance constant */
- if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7)) {
- char *penalty_type;
- penalty_type = mxArrayToString(prhs[4]); /* choosing TV penalty: 'iso' or 'l1', 'iso' is the default */
- if ((strcmp(penalty_type, "l1") != 0) && (strcmp(penalty_type, "iso") != 0)) mexErrMsgTxt("Choose TV type: 'iso' or 'l1',");
- if (strcmp(penalty_type, "l1") == 0) methTV = 1; /* enable 'l1' penalty */
- mxFree(penalty_type);
- }
- if ((nrhs == 6) || (nrhs == 7)) {
- nonneg = (int) mxGetScalar(prhs[5]);
- if ((nonneg != 0) && (nonneg != 1)) mexErrMsgTxt("Nonnegativity constraint can be enabled by choosing 1 or off - 0");
- }
- if (nrhs == 7) {
- printswitch = (int) mxGetScalar(prhs[6]);
- if ((printswitch != 0) && (printswitch != 1)) mexErrMsgTxt("Print can be enabled by choosing 1 or off - 0");
- }
-
- /*Handling Matlab output data*/
- dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
-
- if (number_of_dims == 2) {
- dimZ = 1; /*2D case*/
- Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
- }
- if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-
- /* running the function */
- TV_FGP_GPU_main(Input, Output, lambda, iter, epsil, methTV, nonneg, printswitch, dimX, dimY, dimZ);
-}
\ No newline at end of file
diff --git a/Wrappers/Matlab/mex_compile/regularizers_GPU/ROF_TV_GPU.cpp b/Wrappers/Matlab/mex_compile/regularizers_GPU/ROF_TV_GPU.cpp
deleted file mode 100644
index 7bbe3af..0000000
--- a/Wrappers/Matlab/mex_compile/regularizers_GPU/ROF_TV_GPU.cpp
+++ /dev/null
@@ -1,73 +0,0 @@
-/*
- * 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 "TV_ROF_GPU_core.h"
-
-/* ROF-TV denoising/regularization model [1] (2D/3D case)
- * (MEX wrapper for MATLAB)
- *
- * Input Parameters:
- * 1. Noisy image/volume [REQUIRED]
- * 2. lambda - regularization parameter [REQUIRED]
- * 3. Number of iterations, for explicit scheme >= 150 is recommended [REQUIRED]
- * 4. tau - marching step for explicit scheme, ~1 is recommended [REQUIRED]
- *
- * Output:
- * [1] Regularized image/volume
- *
- * This function is based on the paper by
- * [1] Rudin, Osher, Fatemi, "Nonlinear Total Variation based noise removal algorithms"
- *
- * D. Kazantsev, 2016-18
- */
-
-void mexFunction(
- int nlhs, mxArray *plhs[],
- int nrhs, const mxArray *prhs[])
-
-{
- int number_of_dims, iter_numb, dimX, dimY, dimZ;
- const int *dim_array;
- float *Input, *Output=NULL, lambda, tau;
-
- dim_array = mxGetDimensions(prhs[0]);
- number_of_dims = mxGetNumberOfDimensions(prhs[0]);
-
- /*Handling Matlab input data*/
- Input = (float *) mxGetData(prhs[0]);
- lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */
- iter_numb = (int) mxGetScalar(prhs[2]); /* iterations number */
- tau = (float) mxGetScalar(prhs[3]); /* marching step parameter */
-
- if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
- if(nrhs != 4) mexErrMsgTxt("Four inputs reqired: Image(2D,3D), regularization parameter, iterations number, marching step constant");
- /*Handling Matlab output data*/
- 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 */
- Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
- }
- if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-
- TV_ROF_GPU_main(Input, Output, lambda, iter_numb, tau, dimX, dimY, dimZ);
-}
\ No newline at end of file
diff --git a/Wrappers/Python/CMakeLists.txt b/Wrappers/Python/CMakeLists.txt
index eb3bac7..fb00706 100644
--- a/Wrappers/Python/CMakeLists.txt
+++ b/Wrappers/Python/CMakeLists.txt
@@ -1,7 +1,7 @@
# Copyright 2018 Edoardo Pasca
cmake_minimum_required (VERSION 3.0)
-project(RegularizerPython)
+project(regulariserPython)
#https://stackoverflow.com/questions/13298504/using-cmake-with-setup-py
# The version number.
@@ -25,8 +25,8 @@ if (PYTHONINTERP_FOUND)
endif()
-## Build the regularizers package as a library
-message("Creating Regularizers as shared library")
+## Build the regularisers package as a library
+message("Creating Regularisers as shared library")
message("CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS}")
@@ -52,7 +52,7 @@ elseif(UNIX)
)
endif()
-# GPU Regularizers
+# GPU regularisers
find_package(CUDA)
if (CUDA_FOUND)
@@ -60,12 +60,12 @@ if (CUDA_FOUND)
set (SETUP_GPU_WRAPPERS "extra_libraries += ['cilregcuda']\n\
setup( \n\
name='ccpi', \n\
- description='CCPi Core Imaging Library - Image Regularizers GPU',\n\
+ description='CCPi Core Imaging Library - Image regularisers GPU',\n\
version=cil_version,\n\
cmdclass = {'build_ext': build_ext},\n\
- ext_modules = [Extension('ccpi.filters.gpu_regularizers',\n\
+ ext_modules = [Extension('ccpi.filters.gpu_regularisers',\n\
sources=[ \n\
- os.path.join('.' , 'src', 'gpu_regularizers.pyx' ),\n\
+ os.path.join('.' , 'src', 'gpu_regularisers.pyx' ),\n\
],\n\
include_dirs=extra_include_dirs, \n\
library_dirs=extra_library_dirs, \n\
@@ -80,9 +80,9 @@ else()
set(SETUP_GPU_WRAPPERS "#CUDA NOT FOUND")
endif()
-configure_file("setup-regularizers.py.in" "setup-regularizers.py")
+configure_file("setup-regularisers.py.in" "setup-regularisers.py")
-#add_executable(regularizer_test ${CMAKE_CURRENT_SOURCE_DIR}/test/test_regularizer.cpp)
+#add_executable(regulariser_test ${CMAKE_CURRENT_SOURCE_DIR}/test/test_regulariser.cpp)
-#target_link_libraries (regularizer_test LINK_PUBLIC regularizers_lib)
+#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
new file mode 100644
index 0000000..039daab
--- /dev/null
+++ b/Wrappers/Python/ccpi/filters/regularisers.py
@@ -0,0 +1,44 @@
+"""
+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
+from ccpi.filters.gpu_regularisers import TV_ROF_GPU, TV_FGP_GPU
+
+def ROF_TV(inputData, regularisation_parameter, iterations,
+ time_marching_parameter,device='cpu'):
+ if device == 'cpu':
+ return TV_ROF_CPU(inputData,
+ regularisation_parameter,
+ iterations,
+ time_marching_parameter)
+ elif device == 'gpu':
+ return TV_ROF_GPU(inputData,
+ regularisation_parameter,
+ iterations,
+ time_marching_parameter)
+ else:
+ raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
+ .format(device))
+
+def FGP_TV(inputData, regularisation_parameter,iterations,
+ tolerance_param, methodTV, nonneg, printM, device='cpu'):
+ if device == 'cpu':
+ return TV_FGP_CPU(inputData,
+ regularisation_parameter,
+ iterations,
+ tolerance_param,
+ methodTV,
+ nonneg,
+ printM)
+ elif device == 'gpu':
+ return TV_FGP_GPU(inputData,
+ regularisation_parameter,
+ iterations,
+ tolerance_param,
+ methodTV,
+ nonneg,
+ printM)
+ else:
+ raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
+ .format(device))
diff --git a/Wrappers/Python/ccpi/filters/regularizers.py b/Wrappers/Python/ccpi/filters/regularizers.py
deleted file mode 100644
index d6dfa8c..0000000
--- a/Wrappers/Python/ccpi/filters/regularizers.py
+++ /dev/null
@@ -1,44 +0,0 @@
-"""
-script which assigns a proper device core function based on a flag ('cpu' or 'gpu')
-"""
-
-from ccpi.filters.cpu_regularizers_cython import TV_ROF_CPU, TV_FGP_CPU
-from ccpi.filters.gpu_regularizers import TV_ROF_GPU, TV_FGP_GPU
-
-def ROF_TV(inputData, regularization_parameter, iterations,
- time_marching_parameter,device='cpu'):
- if device == 'cpu':
- return TV_ROF_CPU(inputData,
- regularization_parameter,
- iterations,
- time_marching_parameter)
- elif device == 'gpu':
- return TV_ROF_GPU(inputData,
- regularization_parameter,
- iterations,
- time_marching_parameter)
- else:
- raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
- .format(device))
-
-def FGP_TV(inputData, regularization_parameter,iterations,
- tolerance_param, methodTV, nonneg, printM, device='cpu'):
- if device == 'cpu':
- return TV_FGP_CPU(inputData,
- regularization_parameter,
- iterations,
- tolerance_param,
- methodTV,
- nonneg,
- printM)
- elif device == 'gpu':
- return TV_FGP_GPU(inputData,
- regularization_parameter,
- iterations,
- tolerance_param,
- methodTV,
- nonneg,
- printM)
- else:
- raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
- .format(device))
\ No newline at end of file
diff --git a/Wrappers/Python/conda-recipe/bld.bat b/Wrappers/Python/conda-recipe/bld.bat
index 850905c..e47f8d9 100644
--- a/Wrappers/Python/conda-recipe/bld.bat
+++ b/Wrappers/Python/conda-recipe/bld.bat
@@ -11,7 +11,7 @@ cd %SRC_DIR%\ccpi\Python
:: issue cmake to create setup.py
cmake .
-%PYTHON% setup-regularizers.py build_ext
+%PYTHON% setup-regularisers.py build_ext
if errorlevel 1 exit 1
-%PYTHON% setup-regularizers.py install
+%PYTHON% setup-regularisers.py install
if errorlevel 1 exit 1
diff --git a/Wrappers/Python/conda-recipe/build.sh b/Wrappers/Python/conda-recipe/build.sh
index 9ea4161..8b05663 100644
--- a/Wrappers/Python/conda-recipe/build.sh
+++ b/Wrappers/Python/conda-recipe/build.sh
@@ -13,7 +13,7 @@ echo "$SRC_DIR/ccpi/Python"
cmake .
-$PYTHON setup-regularizers.py build_ext
-$PYTHON setup-regularizers.py install
+$PYTHON setup-regularisers.py build_ext
+$PYTHON setup-regularisers.py install
diff --git a/Wrappers/Python/conda-recipe/meta.yaml b/Wrappers/Python/conda-recipe/meta.yaml
index f4cb471..5336d14 100644
--- a/Wrappers/Python/conda-recipe/meta.yaml
+++ b/Wrappers/Python/conda-recipe/meta.yaml
@@ -1,5 +1,5 @@
package:
- name: ccpi-regularizer
+ name: ccpi-regulariser
version: {{ environ['CIL_VERSION'] }}
@@ -17,7 +17,7 @@ requirements:
#- boost ==1.64.0
#- boost-cpp ==1.64.0
- cython
- - cil_regularizer
+ - cil_regulariser
- vc 14 # [win and py35]
- vc 9 # [win and py27]
- cmake
@@ -26,7 +26,7 @@ requirements:
- python
- numpy x.x
#- boost ==1.64
- - cil_regularizer
+ - cil_regulariser
- vc 14 # [win and py35]
- vc 9 # [win and py27]
diff --git a/Wrappers/Python/demo/test_cpu_regularisers.py b/Wrappers/Python/demo/test_cpu_regularisers.py
new file mode 100644
index 0000000..4e4a2dd
--- /dev/null
+++ b/Wrappers/Python/demo/test_cpu_regularisers.py
@@ -0,0 +1,323 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Aug 4 11:10:05 2017
+
+@author: ofn77899
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+from enum import Enum
+import timeit
+from ccpi.filters.regularisers import ROF_TV, FGP_TV
+###############################################################################
+def nrmse(im1, im2):
+ a, b = im1.shape
+ rmse = np.sqrt(np.sum((im2 - im1) ** 2) / float(a * b))
+ max_val = max(np.max(im1), np.max(im2))
+ min_val = min(np.min(im1), np.min(im2))
+ return 1 - (rmse / (max_val - min_val))
+
+def rmse(im1, im2):
+ a, b = im1.shape
+ rmse = np.sqrt(np.sum((im1 - im2) ** 2) / float(a * b))
+ return 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))
+ else:
+ txt += "{0} = {1}".format(key, value)
+ txt += '\n'
+ return txt
+###############################################################################
+#
+# 2D Regularizers
+#
+###############################################################################
+# assumes the script is launched from the test directory
+filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif")
+
+Im = plt.imread(filename)
+Im = np.asarray(Im, dtype='float32')
+
+Im = Im/255
+
+perc = 0.075
+u0 = Im + np.random.normal(loc = Im ,
+ scale = perc * Im ,
+ size = np.shape(Im))
+# map the u0 u0->u0>0
+f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
+u0 = f(u0).astype('float32')
+
+## plot
+fig = plt.figure()
+
+a=fig.add_subplot(2,4,1)
+a.set_title('noise')
+imgplot = plt.imshow(u0,cmap="gray"
+ )
+
+reg_output = []
+##############################################################################
+# Call regularisers
+
+###################### FGP_TV #########################################
+
+start_time = timeit.default_timer()
+pars = {'algorithm' : FGP_TV , \
+ 'input' : u0,\
+ 'regularisation_parameter':0.07, \
+ 'number_of_iterations' :300 ,\
+ 'tolerance_constant':0.00001,\
+ 'methodTV': 0 ,\
+ 'nonneg': 0 ,\
+ 'printingOut': 0
+ }
+
+fgp = FGP_TV(pars['input'],
+ pars['regularisation_parameter'],
+ pars['number_of_iterations'],
+ pars['tolerance_constant'],
+ pars['methodTV'],
+ pars['nonneg'],
+ pars['printingOut'], 'cpu')
+
+rms = rmse(Im, fgp)
+pars['rmse'] = rms
+
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+
+
+a=fig.add_subplot(2,4,3)
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+imgplot = plt.imshow(fgp, \
+ cmap="gray"
+ )
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
+ verticalalignment='top', bbox=props)
+
+# ###################### ROF_TV #########################################
+
+start_time = timeit.default_timer()
+
+pars = {'algorithm': ROF_TV , \
+ 'input' : u0,\
+ 'regularisation_parameter':0.07,\
+ 'marching_step': 0.0025,\
+ 'number_of_iterations': 300
+ }
+rof = ROF_TV(pars['input'],
+ pars['regularisation_parameter'],
+ pars['number_of_iterations'],
+ pars['marching_step'], 'cpu')
+
+rms = rmse(Im, rof)
+pars['rmse'] = rms
+
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(2,4,7)
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
+ verticalalignment='top', bbox=props)
+imgplot = plt.imshow(rof, cmap="gray")
+
+plt.show()
+
+################################################################################
+##
+## 3D Regularizers
+##
+################################################################################
+##Example:
+## figure;
+## Im = double(imread('lena_gray_256.tif'))/255; % loading image
+## u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
+## u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+#
+##filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Reconstruction\python\test\reconstruction_example.mha"
+#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Simpleflex\data\head.mha"
+#
+#reader = vtk.vtkMetaImageReader()
+#reader.SetFileName(os.path.normpath(filename))
+#reader.Update()
+##vtk returns 3D images, let's take just the one slice there is as 2D
+#Im = Converter.vtk2numpy(reader.GetOutput())
+#Im = Im.astype('float32')
+##imgplot = plt.imshow(Im)
+#perc = 0.05
+#u0 = Im + (perc* np.random.normal(size=np.shape(Im)))
+## map the u0 u0->u0>0
+#f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
+#u0 = f(u0).astype('float32')
+#converter = Converter.numpy2vtkImporter(u0, reader.GetOutput().GetSpacing(),
+# reader.GetOutput().GetOrigin())
+#converter.Update()
+#writer = vtk.vtkMetaImageWriter()
+#writer.SetInputData(converter.GetOutput())
+#writer.SetFileName(r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\noisy_head.mha")
+##writer.Write()
+#
+#
+### plot
+#fig3D = plt.figure()
+##a=fig.add_subplot(3,3,1)
+##a.set_title('Original')
+##imgplot = plt.imshow(Im)
+#sliceNo = 32
+#
+#a=fig3D.add_subplot(2,4,1)
+#a.set_title('noise')
+#imgplot = plt.imshow(u0.T[sliceNo])
+#
+#reg_output3d = []
+#
+###############################################################################
+## Call regularizer
+#
+######################## SplitBregman_TV #####################################
+## u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+#
+##reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV)
+#
+##out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30,
+## #tolerance_constant=1e-4,
+## TV_Penalty=Regularizer.TotalVariationPenalty.l1)
+#
+#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30,
+# tolerance_constant=1e-4,
+# TV_Penalty=Regularizer.TotalVariationPenalty.l1)
+#
+#
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,4,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+# verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+#
+####################### FGP_TV #########################################
+## u = FGP_TV(single(u0), 0.05, 100, 1e-04);
+#out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.005,
+# number_of_iterations=200)
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,4,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+# verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+#
+####################### LLT_model #########################################
+## * u0 = Im + .03*randn(size(Im)); % adding noise
+## [Den] = LLT_model(single(u0), 10, 0.1, 1);
+##Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0);
+##input, regularization_parameter , time_step, number_of_iterations,
+## tolerance_constant, restrictive_Z_smoothing=0
+#out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25,
+# time_step=0.0003,
+# tolerance_constant=0.0001,
+# number_of_iterations=300)
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,4,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+# verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+#
+####################### PatchBased_Regul #########################################
+## Quick 2D denoising example in Matlab:
+## Im = double(imread('lena_gray_256.tif'))/255; % loading image
+## u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+## ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05);
+#
+#out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05,
+# searching_window_ratio=3,
+# similarity_window_ratio=1,
+# PB_filtering_parameter=0.08)
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,4,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+# verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+#
+
+###################### TGV_PD #########################################
+# Quick 2D denoising example in Matlab:
+# Im = double(imread('lena_gray_256.tif'))/255; % loading image
+# u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+# u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550);
+
+
+#out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05,
+# first_order_term=1.3,
+# second_order_term=1,
+# number_of_iterations=550)
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,4,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+# verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
diff --git a/Wrappers/Python/demo/test_cpu_regularizers.py b/Wrappers/Python/demo/test_cpu_regularizers.py
deleted file mode 100644
index 76b9de7..0000000
--- a/Wrappers/Python/demo/test_cpu_regularizers.py
+++ /dev/null
@@ -1,488 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Fri Aug 4 11:10:05 2017
-
-@author: ofn77899
-"""
-
-
-import matplotlib.pyplot as plt
-import numpy as np
-import os
-from enum import Enum
-import timeit
-from ccpi.filters.cpu_regularizers_boost import SplitBregman_TV, LLT_model, PatchBased_Regul, TGV_PD
-from ccpi.filters.regularizers import ROF_TV, FGP_TV
-###############################################################################
-#https://stackoverflow.com/questions/13875989/comparing-image-in-url-to-image-in-filesystem-in-python/13884956#13884956
-#NRMSE a normalization of the root of the mean squared error
-#NRMSE is simply 1 - [RMSE / (maxval - minval)]. Where maxval is the maximum
-# intensity from the two images being compared, and respectively the same for
-# minval. RMSE is given by the square root of MSE:
-# sqrt[(sum(A - B) ** 2) / |A|],
-# where |A| means the number of elements in A. By doing this, the maximum value
-# given by RMSE is maxval.
-
-def nrmse(im1, im2):
- a, b = im1.shape
- rmse = np.sqrt(np.sum((im2 - im1) ** 2) / float(a * b))
- max_val = max(np.max(im1), np.max(im2))
- min_val = min(np.min(im1), np.min(im2))
- return 1 - (rmse / (max_val - min_val))
-
-def rmse(im1, im2):
- a, b = im1.shape
- rmse = np.sqrt(np.sum((im1 - im2) ** 2) / float(a * b))
- return 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))
- else:
- txt += "{0} = {1}".format(key, value)
- txt += '\n'
- return txt
-###############################################################################
-#
-# 2D Regularizers
-#
-###############################################################################
-#Example:
-# figure;
-# Im = double(imread('lena_gray_256.tif'))/255; % loading image
-# u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
-# u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
-
-# assumes the script is launched from the test directory
-filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif")
-#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\lena_gray_512.tif"
-#filename = r"/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/data/lena_gray_512.tif"
-#filename = r'/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif'
-
-Im = plt.imread(filename)
-Im = np.asarray(Im, dtype='float32')
-
-Im = Im/255
-
-perc = 0.075
-u0 = Im + np.random.normal(loc = Im ,
- scale = perc * Im ,
- size = np.shape(Im))
-# map the u0 u0->u0>0
-f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
-u0 = f(u0).astype('float32')
-
-## plot
-fig = plt.figure()
-
-a=fig.add_subplot(2,4,1)
-a.set_title('noise')
-imgplot = plt.imshow(u0,cmap="gray"
- )
-
-reg_output = []
-##############################################################################
-# Call regularizer
-
-####################### SplitBregman_TV #####################################
-# u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
-
-start_time = timeit.default_timer()
-pars = {'algorithm' : SplitBregman_TV , \
- 'input' : u0,
- 'regularization_parameter':15. , \
- 'number_of_iterations' :40 ,\
- 'tolerance_constant':0.0001 , \
- 'TV_penalty': 0
-}
-
-out = SplitBregman_TV (pars['input'], pars['regularization_parameter'],
- pars['number_of_iterations'],
- pars['tolerance_constant'],
- pars['TV_penalty'])
-splitbregman = out[0]
-rms = rmse(Im, splitbregman)
-pars['rmse'] = rms
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-
-a=fig.add_subplot(2,4,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow(splitbregman,\
- cmap="gray"
- )
-
-###################### FGP_TV #########################################
-
-start_time = timeit.default_timer()
-pars = {'algorithm' : FGP_TV , \
- 'input' : u0,\
- 'regularization_parameter':0.07, \
- 'number_of_iterations' :300 ,\
- 'tolerance_constant':0.00001,\
- 'methodTV': 0 ,\
- 'nonneg': 0 ,\
- 'printingOut': 0
- }
-
-fgp = FGP_TV(pars['input'],
- pars['regularization_parameter'],
- pars['number_of_iterations'],
- pars['tolerance_constant'],
- pars['methodTV'],
- pars['nonneg'],
- pars['printingOut'], 'cpu')
-
-rms = rmse(Im, fgp)
-pars['rmse'] = rms
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-
-
-a=fig.add_subplot(2,4,3)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-imgplot = plt.imshow(fgp, \
- cmap="gray"
- )
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
- verticalalignment='top', bbox=props)
-
-
-###################### LLT_model #########################################
-start_time = timeit.default_timer()
-
-pars = {'algorithm': LLT_model , \
- 'input' : u0,
- 'regularization_parameter': 5,\
- 'time_step':0.00035, \
- 'number_of_iterations' :350,\
- 'tolerance_constant':0.0001,\
- 'restrictive_Z_smoothing': 0
-}
-out = LLT_model(pars['input'],
- pars['regularization_parameter'],
- pars['time_step'] ,
- pars['number_of_iterations'],
- pars['tolerance_constant'],
- pars['restrictive_Z_smoothing'] )
-
-llt = out[0]
-rms = rmse(Im, out[0])
-pars['rmse'] = rms
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(2,4,4)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow(llt,\
- cmap="gray"
- )
-# ###################### PatchBased_Regul #########################################
-# # Quick 2D denoising example in Matlab:
-# # Im = double(imread('lena_gray_256.tif'))/255; % loading image
-# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
-# # ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05);
-
-start_time = timeit.default_timer()
-
-pars = {'algorithm': PatchBased_Regul , \
- 'input' : u0,
- 'regularization_parameter': 0.05,\
- 'searching_window_ratio':3, \
- 'similarity_window_ratio':1,\
- 'PB_filtering_parameter': 0.06
-}
-out = PatchBased_Regul(pars['input'],
- pars['regularization_parameter'],
- pars['searching_window_ratio'] ,
- pars['similarity_window_ratio'] ,
- pars['PB_filtering_parameter'])
-pbr = out[0]
-rms = rmse(Im, out[0])
-pars['rmse'] = rms
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-
-a=fig.add_subplot(2,4,5)
-
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow(pbr ,cmap="gray")
-
-
-# ###################### TGV_PD #########################################
-# # Quick 2D denoising example in Matlab:
-# # Im = double(imread('lena_gray_256.tif'))/255; % loading image
-# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
-# # u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550);
-
-start_time = timeit.default_timer()
-
-pars = {'algorithm': TGV_PD , \
- 'input' : u0,\
- 'regularization_parameter':0.07,\
- 'first_order_term': 1.3,\
- 'second_order_term': 1, \
- 'number_of_iterations': 550
- }
-out = TGV_PD(pars['input'],
- pars['regularization_parameter'],
- pars['first_order_term'] ,
- pars['second_order_term'] ,
- pars['number_of_iterations'])
-tgv = out[0]
-rms = rmse(Im, out[0])
-pars['rmse'] = rms
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(2,4,6)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow(tgv, cmap="gray")
-
-# ###################### ROF_TV #########################################
-
-start_time = timeit.default_timer()
-
-pars = {'algorithm': ROF_TV , \
- 'input' : u0,\
- 'regularization_parameter':0.07,\
- 'marching_step': 0.0025,\
- 'number_of_iterations': 300
- }
-rof = ROF_TV(pars['input'],
- pars['regularization_parameter'],
- pars['number_of_iterations'],
- pars['marching_step'], 'cpu')
-
-rms = rmse(Im, rof)
-pars['rmse'] = rms
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(2,4,7)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow(rof, cmap="gray")
-
-plt.show()
-
-################################################################################
-##
-## 3D Regularizers
-##
-################################################################################
-##Example:
-## figure;
-## Im = double(imread('lena_gray_256.tif'))/255; % loading image
-## u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
-## u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
-#
-##filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Reconstruction\python\test\reconstruction_example.mha"
-#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Simpleflex\data\head.mha"
-#
-#reader = vtk.vtkMetaImageReader()
-#reader.SetFileName(os.path.normpath(filename))
-#reader.Update()
-##vtk returns 3D images, let's take just the one slice there is as 2D
-#Im = Converter.vtk2numpy(reader.GetOutput())
-#Im = Im.astype('float32')
-##imgplot = plt.imshow(Im)
-#perc = 0.05
-#u0 = Im + (perc* np.random.normal(size=np.shape(Im)))
-## map the u0 u0->u0>0
-#f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
-#u0 = f(u0).astype('float32')
-#converter = Converter.numpy2vtkImporter(u0, reader.GetOutput().GetSpacing(),
-# reader.GetOutput().GetOrigin())
-#converter.Update()
-#writer = vtk.vtkMetaImageWriter()
-#writer.SetInputData(converter.GetOutput())
-#writer.SetFileName(r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\noisy_head.mha")
-##writer.Write()
-#
-#
-### plot
-#fig3D = plt.figure()
-##a=fig.add_subplot(3,3,1)
-##a.set_title('Original')
-##imgplot = plt.imshow(Im)
-#sliceNo = 32
-#
-#a=fig3D.add_subplot(2,4,1)
-#a.set_title('noise')
-#imgplot = plt.imshow(u0.T[sliceNo])
-#
-#reg_output3d = []
-#
-###############################################################################
-## Call regularizer
-#
-######################## SplitBregman_TV #####################################
-## u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
-#
-##reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV)
-#
-##out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30,
-## #tolerance_constant=1e-4,
-## TV_Penalty=Regularizer.TotalVariationPenalty.l1)
-#
-#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30,
-# tolerance_constant=1e-4,
-# TV_Penalty=Regularizer.TotalVariationPenalty.l1)
-#
-#
-#pars = out2[-2]
-#reg_output3d.append(out2)
-#
-#a=fig3D.add_subplot(2,4,2)
-#
-#
-#textstr = out2[-1]
-#
-#
-## these are matplotlib.patch.Patch properties
-#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-## place a text box in upper left in axes coords
-#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
-# verticalalignment='top', bbox=props)
-#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
-#
-####################### FGP_TV #########################################
-## u = FGP_TV(single(u0), 0.05, 100, 1e-04);
-#out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.005,
-# number_of_iterations=200)
-#pars = out2[-2]
-#reg_output3d.append(out2)
-#
-#a=fig3D.add_subplot(2,4,2)
-#
-#
-#textstr = out2[-1]
-#
-#
-## these are matplotlib.patch.Patch properties
-#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-## place a text box in upper left in axes coords
-#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
-# verticalalignment='top', bbox=props)
-#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
-#
-####################### LLT_model #########################################
-## * u0 = Im + .03*randn(size(Im)); % adding noise
-## [Den] = LLT_model(single(u0), 10, 0.1, 1);
-##Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0);
-##input, regularization_parameter , time_step, number_of_iterations,
-## tolerance_constant, restrictive_Z_smoothing=0
-#out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25,
-# time_step=0.0003,
-# tolerance_constant=0.0001,
-# number_of_iterations=300)
-#pars = out2[-2]
-#reg_output3d.append(out2)
-#
-#a=fig3D.add_subplot(2,4,2)
-#
-#
-#textstr = out2[-1]
-#
-#
-## these are matplotlib.patch.Patch properties
-#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-## place a text box in upper left in axes coords
-#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
-# verticalalignment='top', bbox=props)
-#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
-#
-####################### PatchBased_Regul #########################################
-## Quick 2D denoising example in Matlab:
-## Im = double(imread('lena_gray_256.tif'))/255; % loading image
-## u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
-## ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05);
-#
-#out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05,
-# searching_window_ratio=3,
-# similarity_window_ratio=1,
-# PB_filtering_parameter=0.08)
-#pars = out2[-2]
-#reg_output3d.append(out2)
-#
-#a=fig3D.add_subplot(2,4,2)
-#
-#
-#textstr = out2[-1]
-#
-#
-## these are matplotlib.patch.Patch properties
-#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-## place a text box in upper left in axes coords
-#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
-# verticalalignment='top', bbox=props)
-#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
-#
-
-###################### TGV_PD #########################################
-# Quick 2D denoising example in Matlab:
-# Im = double(imread('lena_gray_256.tif'))/255; % loading image
-# u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
-# u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550);
-
-
-#out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05,
-# first_order_term=1.3,
-# second_order_term=1,
-# number_of_iterations=550)
-#pars = out2[-2]
-#reg_output3d.append(out2)
-#
-#a=fig3D.add_subplot(2,4,2)
-#
-#
-#textstr = out2[-1]
-#
-#
-## these are matplotlib.patch.Patch properties
-#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-## place a text box in upper left in axes coords
-#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
-# verticalalignment='top', bbox=props)
-#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
diff --git a/Wrappers/Python/fista-recipe/bld.bat b/Wrappers/Python/fista-recipe/bld.bat
deleted file mode 100644
index 32c1bc6..0000000
--- a/Wrappers/Python/fista-recipe/bld.bat
+++ /dev/null
@@ -1,13 +0,0 @@
-IF NOT DEFINED CIL_VERSION (
-ECHO CIL_VERSION Not Defined.
-exit 1
-)
-
-mkdir "%SRC_DIR%\ccpifista"
-xcopy /e "%RECIPE_DIR%\.." "%SRC_DIR%\ccpifista"
-
-cd "%SRC_DIR%\ccpifista"
-::%PYTHON% setup-fista.py -q bdist_egg
-:: %PYTHON% setup.py install --single-version-externally-managed --record=record.txt
-%PYTHON% setup-fista.py install
-if errorlevel 1 exit 1
diff --git a/Wrappers/Python/fista-recipe/build.sh b/Wrappers/Python/fista-recipe/build.sh
deleted file mode 100644
index e3f3552..0000000
--- a/Wrappers/Python/fista-recipe/build.sh
+++ /dev/null
@@ -1,10 +0,0 @@
-if [ -z "$CIL_VERSION" ]; then
- echo "Need to set CIL_VERSION"
- exit 1
-fi
-mkdir "$SRC_DIR/ccpifista"
-cp -r "$RECIPE_DIR/.." "$SRC_DIR/ccpifista"
-
-cd $SRC_DIR/ccpifista
-
-$PYTHON setup-fista.py install
diff --git a/Wrappers/Python/fista-recipe/meta.yaml b/Wrappers/Python/fista-recipe/meta.yaml
deleted file mode 100644
index 4e5cba6..0000000
--- a/Wrappers/Python/fista-recipe/meta.yaml
+++ /dev/null
@@ -1,29 +0,0 @@
-package:
- name: ccpi-fista
- version: {{ environ['CIL_VERSION'] }}
-
-
-build:
- preserve_egg_dir: False
- script_env:
- - CIL_VERSION
-# number: 0
-
-requirements:
- build:
- - python
- - numpy
- - setuptools
-
- run:
- - python
- - numpy
- #- astra-toolbox
- - ccpi-regularizer
-
-
-
-about:
- home: http://www.ccpi.ac.uk
- license: Apache v.2.0 license
- summary: 'CCPi Core Imaging Library (Viewer)'
diff --git a/Wrappers/Python/setup-fista.py b/Wrappers/Python/setup-fista.py
deleted file mode 100644
index c5c9f4d..0000000
--- a/Wrappers/Python/setup-fista.py
+++ /dev/null
@@ -1,27 +0,0 @@
-from distutils.core import setup
-#from setuptools import setup, find_packages
-import os
-
-cil_version=os.environ['CIL_VERSION']
-if cil_version == '':
- print("Please set the environmental variable CIL_VERSION")
- sys.exit(1)
-
-setup(
- name="ccpi-fista",
- version=cil_version,
- packages=['ccpi','ccpi.reconstruction'],
- install_requires=['numpy'],
-
- zip_safe = False,
-
- # metadata for upload to PyPI
- author="Edoardo Pasca",
- author_email="edo.paskino@gmail.com",
- description='CCPi Core Imaging Library - FISTA Reconstructor module',
- license="Apache v2.0",
- keywords="tomography interative reconstruction",
- url="http://www.ccpi.ac.uk", # project home page, if any
-
- # could also include long_description, download_url, classifiers, etc.
-)
diff --git a/Wrappers/Python/setup-regularisers.py.in b/Wrappers/Python/setup-regularisers.py.in
new file mode 100644
index 0000000..a1c1ab6
--- /dev/null
+++ b/Wrappers/Python/setup-regularisers.py.in
@@ -0,0 +1,65 @@
+#!/usr/bin/env python
+
+import setuptools
+from distutils.core import setup
+from distutils.extension import Extension
+from Cython.Distutils import build_ext
+
+import os
+import sys
+import numpy
+import platform
+
+cil_version=os.environ['CIL_VERSION']
+if cil_version == '':
+ print("Please set the environmental variable CIL_VERSION")
+ sys.exit(1)
+
+library_include_path = ""
+library_lib_path = ""
+try:
+ library_include_path = os.environ['LIBRARY_INC']
+ library_lib_path = os.environ['LIBRARY_LIB']
+except:
+ library_include_path = os.environ['PREFIX']+'/include'
+ pass
+
+extra_include_dirs = [numpy.get_include(), library_include_path]
+#extra_library_dirs = [os.path.join(library_include_path, "..", "lib")]
+extra_compile_args = []
+extra_library_dirs = []
+extra_compile_args = []
+extra_link_args = []
+extra_libraries = ['cilreg']
+
+extra_include_dirs += [os.path.join(".." , ".." , "Core"),
+ os.path.join(".." , ".." , "Core", "regularisers_CPU"),
+ os.path.join(".." , ".." , "Core", "regularisers_GPU" , "TV_FGP" ) ,
+ os.path.join(".." , ".." , "Core", "regularisers_GPU" , "TV_ROF" ) ,
+ "."]
+
+if platform.system() == 'Windows':
+ extra_compile_args[0:] = ['/DWIN32','/EHsc','/DBOOST_ALL_NO_LIB' , '/openmp' ]
+else:
+ extra_compile_args = ['-fopenmp','-O2', '-funsigned-char', '-Wall', '-std=c++0x']
+ extra_libraries += [@EXTRA_OMP_LIB@]
+
+setup(
+ name='ccpi',
+ description='CCPi Core Imaging Library - Image regularisers',
+ version=cil_version,
+ cmdclass = {'build_ext': build_ext},
+ ext_modules = [Extension("ccpi.filters.cpu_regularisers_cython",
+ sources=[os.path.join("." , "src", "cpu_regularisers.pyx" ) ],
+ include_dirs=extra_include_dirs,
+ library_dirs=extra_library_dirs,
+ extra_compile_args=extra_compile_args,
+ libraries=extra_libraries ),
+
+ ],
+ zip_safe = False,
+ packages = {'ccpi','ccpi.filters'},
+)
+
+
+@SETUP_GPU_WRAPPERS@
diff --git a/Wrappers/Python/setup-regularizers.py.in b/Wrappers/Python/setup-regularizers.py.in
deleted file mode 100644
index 0811372..0000000
--- a/Wrappers/Python/setup-regularizers.py.in
+++ /dev/null
@@ -1,65 +0,0 @@
-#!/usr/bin/env python
-
-import setuptools
-from distutils.core import setup
-from distutils.extension import Extension
-from Cython.Distutils import build_ext
-
-import os
-import sys
-import numpy
-import platform
-
-cil_version=os.environ['CIL_VERSION']
-if cil_version == '':
- print("Please set the environmental variable CIL_VERSION")
- sys.exit(1)
-
-library_include_path = ""
-library_lib_path = ""
-try:
- library_include_path = os.environ['LIBRARY_INC']
- library_lib_path = os.environ['LIBRARY_LIB']
-except:
- library_include_path = os.environ['PREFIX']+'/include'
- pass
-
-extra_include_dirs = [numpy.get_include(), library_include_path]
-#extra_library_dirs = [os.path.join(library_include_path, "..", "lib")]
-extra_compile_args = []
-extra_library_dirs = []
-extra_compile_args = []
-extra_link_args = []
-extra_libraries = ['cilreg']
-
-extra_include_dirs += [os.path.join(".." , ".." , "Core"),
- os.path.join(".." , ".." , "Core", "regularizers_CPU"),
- os.path.join(".." , ".." , "Core", "regularizers_GPU" , "TV_FGP" ) ,
- os.path.join(".." , ".." , "Core", "regularizers_GPU" , "TV_ROF" ) ,
- "."]
-
-if platform.system() == 'Windows':
- extra_compile_args[0:] = ['/DWIN32','/EHsc','/DBOOST_ALL_NO_LIB' , '/openmp' ]
-else:
- extra_compile_args = ['-fopenmp','-O2', '-funsigned-char', '-Wall', '-std=c++0x']
- extra_libraries += [@EXTRA_OMP_LIB@]
-
-setup(
- name='ccpi',
- description='CCPi Core Imaging Library - Image Regularizers',
- version=cil_version,
- cmdclass = {'build_ext': build_ext},
- ext_modules = [Extension("ccpi.filters.cpu_regularizers_cython",
- sources=[os.path.join("." , "src", "cpu_regularizers.pyx" ) ],
- include_dirs=extra_include_dirs,
- library_dirs=extra_library_dirs,
- extra_compile_args=extra_compile_args,
- libraries=extra_libraries ),
-
- ],
- zip_safe = False,
- packages = {'ccpi','ccpi.filters'},
-)
-
-
-@SETUP_GPU_WRAPPERS@
diff --git a/Wrappers/Python/src/cpu_regularisers.pyx b/Wrappers/Python/src/cpu_regularisers.pyx
new file mode 100644
index 0000000..248bad1
--- /dev/null
+++ b/Wrappers/Python/src/cpu_regularisers.pyx
@@ -0,0 +1,126 @@
+# distutils: language=c++
+"""
+Copyright 2018 CCPi
+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.
+
+Author: Edoardo Pasca, Daniil Kazantsev
+"""
+
+import cython
+import numpy as np
+cimport numpy as np
+
+cdef extern float TV_ROF_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float tau, int dimX, int dimY, int dimZ);
+cdef extern float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ);
+
+
+#****************************************************************#
+#********************** Total-variation ROF *********************#
+#****************************************************************#
+def TV_ROF_CPU(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter):
+ if inputData.ndim == 2:
+ return TV_ROF_2D(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter)
+ elif inputData.ndim == 3:
+ return TV_ROF_3D(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter)
+
+def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
+ float regularisation_parameter,
+ int iterationsNumb,
+ float marching_step_parameter):
+ 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 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)
+
+ return outputData
+
+def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
+ int iterationsNumb,
+ float regularisation_parameter,
+ float marching_step_parameter):
+ 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 ROF iterations for 3D data
+ TV_ROF_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, iterationsNumb, marching_step_parameter, dims[0], dims[1], dims[2])
+
+ return outputData
+
+#****************************************************************#
+#********************** Total-variation FGP *********************#
+#****************************************************************#
+#******** Total-variation Fast-Gradient-Projection (FGP)*********#
+def TV_FGP_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM):
+ if inputData.ndim == 2:
+ return TV_FGP_2D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM)
+ elif inputData.ndim == 3:
+ return TV_FGP_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM)
+
+def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
+ float regularisation_parameter,
+ int iterationsNumb,
+ float tolerance_param,
+ int methodTV,
+ int nonneg,
+ int printM):
+
+ 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 ROF iterations for 2D data */
+ TV_FGP_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter,
+ iterationsNumb,
+ tolerance_param,
+ methodTV,
+ nonneg,
+ printM,
+ dims[0], dims[1], 1)
+
+ return outputData
+
+def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
+ float regularisation_parameter,
+ int iterationsNumb,
+ float tolerance_param,
+ int methodTV,
+ int nonneg,
+ int printM):
+ 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 ROF iterations for 3D data */
+ TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter,
+ iterationsNumb,
+ tolerance_param,
+ methodTV,
+ nonneg,
+ printM,
+ dims[0], dims[1], dims[2])
+ return outputData
diff --git a/Wrappers/Python/src/cpu_regularizers.pyx b/Wrappers/Python/src/cpu_regularizers.pyx
deleted file mode 100644
index f993e54..0000000
--- a/Wrappers/Python/src/cpu_regularizers.pyx
+++ /dev/null
@@ -1,126 +0,0 @@
-# distutils: language=c++
-"""
-Copyright 2018 CCPi
-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.
-
-Author: Edoardo Pasca, Daniil Kazantsev
-"""
-
-import cython
-import numpy as np
-cimport numpy as np
-
-cdef extern float TV_ROF_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float tau, int dimX, int dimY, int dimZ);
-cdef extern float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ);
-
-
-#****************************************************************#
-#********************** Total-variation ROF *********************#
-#****************************************************************#
-def TV_ROF_CPU(inputData, regularization_parameter, iterationsNumb, marching_step_parameter):
- if inputData.ndim == 2:
- return TV_ROF_2D(inputData, regularization_parameter, iterationsNumb, marching_step_parameter)
- elif inputData.ndim == 3:
- return TV_ROF_3D(inputData, regularization_parameter, iterationsNumb, marching_step_parameter)
-
-def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
- float regularization_parameter,
- int iterationsNumb,
- float marching_step_parameter):
- 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 ROF iterations for 2D data
- TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], regularization_parameter, iterationsNumb, marching_step_parameter, dims[0], dims[1], 1)
-
- return outputData
-
-def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
- int iterationsNumb,
- float regularization_parameter,
- float marching_step_parameter):
- 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 ROF iterations for 3D data
- TV_ROF_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularization_parameter, iterationsNumb, marching_step_parameter, dims[0], dims[1], dims[2])
-
- return outputData
-
-#****************************************************************#
-#********************** Total-variation FGP *********************#
-#****************************************************************#
-#******** Total-variation Fast-Gradient-Projection (FGP)*********#
-def TV_FGP_CPU(inputData, regularization_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM):
- if inputData.ndim == 2:
- return TV_FGP_2D(inputData, regularization_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM)
- elif inputData.ndim == 3:
- return TV_FGP_3D(inputData, regularization_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM)
-
-def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
- float regularization_parameter,
- int iterationsNumb,
- float tolerance_param,
- int methodTV,
- int nonneg,
- int printM):
-
- 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 ROF iterations for 2D data */
- TV_FGP_CPU_main(&inputData[0,0], &outputData[0,0], regularization_parameter,
- iterationsNumb,
- tolerance_param,
- methodTV,
- nonneg,
- printM,
- dims[0], dims[1], 1)
-
- return outputData
-
-def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
- float regularization_parameter,
- int iterationsNumb,
- float tolerance_param,
- int methodTV,
- int nonneg,
- int printM):
- 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 ROF iterations for 3D data */
- TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularization_parameter,
- iterationsNumb,
- tolerance_param,
- methodTV,
- nonneg,
- printM,
- dims[0], dims[1], dims[2])
- return outputData
diff --git a/Wrappers/Python/src/gpu_regularisers.pyx b/Wrappers/Python/src/gpu_regularisers.pyx
new file mode 100644
index 0000000..7ebd011
--- /dev/null
+++ b/Wrappers/Python/src/gpu_regularisers.pyx
@@ -0,0 +1,170 @@
+# distutils: language=c++
+"""
+Copyright 2018 CCPi
+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.
+
+Author: Edoardo Pasca, Daniil Kazantsev
+"""
+
+import cython
+import numpy as np
+cimport numpy as np
+
+cdef extern void TV_ROF_GPU_main(float* Input, float* Output, float lambdaPar, int iter, float tau, int N, int M, int Z);
+cdef extern void TV_FGP_GPU_main(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int N, int M, int Z);
+
+# Total-variation Rudin-Osher-Fatemi (ROF)
+def TV_ROF_GPU(inputData,
+ regularisation_parameter,
+ iterations,
+ time_marching_parameter):
+ if inputData.ndim == 2:
+ return ROFTV2D(inputData,
+ regularisation_parameter,
+ iterations,
+ time_marching_parameter)
+ elif inputData.ndim == 3:
+ return ROFTV3D(inputData,
+ regularisation_parameter,
+ iterations,
+ time_marching_parameter)
+
+# Total-variation Fast-Gradient-Projection (FGP)
+def TV_FGP_GPU(inputData,
+ regularisation_parameter,
+ iterations,
+ tolerance_param,
+ methodTV,
+ nonneg,
+ printM):
+ if inputData.ndim == 2:
+ return FGPTV2D(inputData,
+ regularisation_parameter,
+ iterations,
+ tolerance_param,
+ methodTV,
+ nonneg,
+ printM)
+ elif inputData.ndim == 3:
+ return FGPTV3D(inputData,
+ regularisation_parameter,
+ iterations,
+ tolerance_param,
+ methodTV,
+ nonneg,
+ printM)
+
+#****************************************************************#
+#********************** Total-variation ROF *********************#
+#****************************************************************#
+def ROFTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
+ float regularisation_parameter,
+ int iterations,
+ float time_marching_parameter):
+
+ 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')
+
+ # Running CUDA code here
+ TV_ROF_GPU_main(
+ &inputData[0,0], &outputData[0,0],
+ regularisation_parameter,
+ iterations ,
+ time_marching_parameter,
+ dims[0], dims[1], 1);
+
+ return outputData
+
+def ROFTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
+ float regularisation_parameter,
+ int iterations,
+ float time_marching_parameter):
+
+ 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')
+
+ # Running CUDA code here
+ TV_ROF_GPU_main(
+ &inputData[0,0,0], &outputData[0,0,0],
+ regularisation_parameter,
+ iterations ,
+ time_marching_parameter,
+ dims[0], dims[1], dims[2]);
+
+ return outputData
+#****************************************************************#
+#********************** Total-variation FGP *********************#
+#****************************************************************#
+#******** Total-variation Fast-Gradient-Projection (FGP)*********#
+def FGPTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
+ float regularisation_parameter,
+ int iterations,
+ float tolerance_param,
+ int methodTV,
+ int nonneg,
+ int printM):
+
+ 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')
+
+ # Running CUDA code here
+ TV_FGP_GPU_main(&inputData[0,0], &outputData[0,0],
+ regularisation_parameter,
+ iterations,
+ tolerance_param,
+ methodTV,
+ nonneg,
+ printM,
+ dims[0], dims[1], 1);
+
+ return outputData
+
+def FGPTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
+ float regularisation_parameter,
+ int iterations,
+ float tolerance_param,
+ int methodTV,
+ int nonneg,
+ int printM):
+
+ 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')
+
+ # Running CUDA code here
+ TV_FGP_GPU_main(
+ &inputData[0,0,0], &outputData[0,0,0],
+ regularisation_parameter ,
+ iterations,
+ tolerance_param,
+ methodTV,
+ nonneg,
+ printM,
+ dims[0], dims[1], dims[2]);
+
+ return outputData
diff --git a/Wrappers/Python/src/gpu_regularizers.pyx b/Wrappers/Python/src/gpu_regularizers.pyx
deleted file mode 100644
index a44bd1d..0000000
--- a/Wrappers/Python/src/gpu_regularizers.pyx
+++ /dev/null
@@ -1,170 +0,0 @@
-# distutils: language=c++
-"""
-Copyright 2018 CCPi
-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.
-
-Author: Edoardo Pasca, Daniil Kazantsev
-"""
-
-import cython
-import numpy as np
-cimport numpy as np
-
-cdef extern void TV_ROF_GPU_main(float* Input, float* Output, float lambdaPar, int iter, float tau, int N, int M, int Z);
-cdef extern void TV_FGP_GPU_main(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int N, int M, int Z);
-
-# Total-variation Rudin-Osher-Fatemi (ROF)
-def TV_ROF_GPU(inputData,
- regularization_parameter,
- iterations,
- time_marching_parameter):
- if inputData.ndim == 2:
- return ROFTV2D(inputData,
- regularization_parameter,
- iterations,
- time_marching_parameter)
- elif inputData.ndim == 3:
- return ROFTV3D(inputData,
- regularization_parameter,
- iterations,
- time_marching_parameter)
-
-# Total-variation Fast-Gradient-Projection (FGP)
-def TV_FGP_GPU(inputData,
- regularization_parameter,
- iterations,
- tolerance_param,
- methodTV,
- nonneg,
- printM):
- if inputData.ndim == 2:
- return FGPTV2D(inputData,
- regularization_parameter,
- iterations,
- tolerance_param,
- methodTV,
- nonneg,
- printM)
- elif inputData.ndim == 3:
- return FGPTV3D(inputData,
- regularization_parameter,
- iterations,
- tolerance_param,
- methodTV,
- nonneg,
- printM)
-
-#****************************************************************#
-#********************** Total-variation ROF *********************#
-#****************************************************************#
-def ROFTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
- float regularization_parameter,
- int iterations,
- float time_marching_parameter):
-
- 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')
-
- # Running CUDA code here
- TV_ROF_GPU_main(
- &inputData[0,0], &outputData[0,0],
- regularization_parameter,
- iterations ,
- time_marching_parameter,
- dims[0], dims[1], 1);
-
- return outputData
-
-def ROFTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
- float regularization_parameter,
- int iterations,
- float time_marching_parameter):
-
- 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')
-
- # Running CUDA code here
- TV_ROF_GPU_main(
- &inputData[0,0,0], &outputData[0,0,0],
- regularization_parameter,
- iterations ,
- time_marching_parameter,
- dims[0], dims[1], dims[2]);
-
- return outputData
-#****************************************************************#
-#********************** Total-variation FGP *********************#
-#****************************************************************#
-#******** Total-variation Fast-Gradient-Projection (FGP)*********#
-def FGPTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
- float regularization_parameter,
- int iterations,
- float tolerance_param,
- int methodTV,
- int nonneg,
- int printM):
-
- 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')
-
- # Running CUDA code here
- TV_FGP_GPU_main(&inputData[0,0], &outputData[0,0],
- regularization_parameter,
- iterations,
- tolerance_param,
- methodTV,
- nonneg,
- printM,
- dims[0], dims[1], 1);
-
- return outputData
-
-def FGPTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
- float regularization_parameter,
- int iterations,
- float tolerance_param,
- int methodTV,
- int nonneg,
- int printM):
-
- 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')
-
- # Running CUDA code here
- TV_FGP_GPU_main(
- &inputData[0,0,0], &outputData[0,0,0],
- regularization_parameter ,
- iterations,
- tolerance_param,
- methodTV,
- nonneg,
- printM,
- dims[0], dims[1], dims[2]);
-
- return outputData
diff --git a/Wrappers/Python/test/test_cpu_regularisers.py b/Wrappers/Python/test/test_cpu_regularisers.py
new file mode 100644
index 0000000..9713baa
--- /dev/null
+++ b/Wrappers/Python/test/test_cpu_regularisers.py
@@ -0,0 +1,442 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Aug 4 11:10:05 2017
+
+@author: ofn77899
+"""
+
+
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+from enum import Enum
+import timeit
+from ccpi.filters.cpu_regularizers_boost import SplitBregman_TV , FGP_TV ,\
+ LLT_model, PatchBased_Regul ,\
+ TGV_PD
+
+###############################################################################
+#https://stackoverflow.com/questions/13875989/comparing-image-in-url-to-image-in-filesystem-in-python/13884956#13884956
+#NRMSE a normalization of the root of the mean squared error
+#NRMSE is simply 1 - [RMSE / (maxval - minval)]. Where maxval is the maximum
+# intensity from the two images being compared, and respectively the same for
+# minval. RMSE is given by the square root of MSE:
+# sqrt[(sum(A - B) ** 2) / |A|],
+# where |A| means the number of elements in A. By doing this, the maximum value
+# given by RMSE is maxval.
+
+def nrmse(im1, im2):
+ a, b = im1.shape
+ rmse = np.sqrt(np.sum((im2 - im1) ** 2) / float(a * b))
+ max_val = max(np.max(im1), np.max(im2))
+ min_val = min(np.min(im1), np.min(im2))
+ return 1 - (rmse / (max_val - min_val))
+###############################################################################
+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))
+ else:
+ txt += "{0} = {1}".format(key, value)
+ txt += '\n'
+ return txt
+###############################################################################
+#
+# 2D Regularizers
+#
+###############################################################################
+#Example:
+# figure;
+# Im = double(imread('lena_gray_256.tif'))/255; % loading image
+# u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
+# u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+
+# assumes the script is launched from the test directory
+filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif")
+#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\lena_gray_512.tif"
+#filename = r"/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/data/lena_gray_512.tif"
+#filename = r'/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif'
+
+Im = plt.imread(filename)
+Im = np.asarray(Im, dtype='float32')
+
+perc = 0.15
+u0 = Im + np.random.normal(loc = Im ,
+ scale = perc * Im ,
+ size = np.shape(Im))
+# map the u0 u0->u0>0
+f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
+u0 = f(u0).astype('float32')
+
+## plot
+fig = plt.figure()
+
+a=fig.add_subplot(2,3,1)
+a.set_title('noise')
+imgplot = plt.imshow(u0#,cmap="gray"
+ )
+
+reg_output = []
+##############################################################################
+# Call regularizer
+
+####################### SplitBregman_TV #####################################
+# u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+
+start_time = timeit.default_timer()
+pars = {'algorithm' : SplitBregman_TV , \
+ 'input' : u0,
+ 'regularization_parameter':10. , \
+'number_of_iterations' :35 ,\
+'tolerance_constant':0.0001 , \
+'TV_penalty': 0
+}
+
+out = SplitBregman_TV (pars['input'], pars['regularization_parameter'],
+ pars['number_of_iterations'],
+ pars['tolerance_constant'],
+ pars['TV_penalty'])
+splitbregman = out[0]
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+
+
+a=fig.add_subplot(2,3,2)
+
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
+ verticalalignment='top', bbox=props)
+imgplot = plt.imshow(splitbregman,\
+ #cmap="gray"
+ )
+
+###################### FGP_TV #########################################
+# u = FGP_TV(single(u0), 0.05, 100, 1e-04);
+start_time = timeit.default_timer()
+pars = {'algorithm' : FGP_TV , \
+ 'input' : u0,
+ 'regularization_parameter':5e-4, \
+ 'number_of_iterations' :10 ,\
+ 'tolerance_constant':0.001,\
+ 'TV_penalty': 0
+}
+
+out = FGP_TV (pars['input'],
+ pars['regularization_parameter'],
+ pars['number_of_iterations'],
+ pars['tolerance_constant'],
+ pars['TV_penalty'])
+
+fgp = out[0]
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+
+
+a=fig.add_subplot(2,3,3)
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+imgplot = plt.imshow(fgp, \
+ #cmap="gray"
+ )
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
+ verticalalignment='top', bbox=props)
+
+###################### LLT_model #########################################
+
+start_time = timeit.default_timer()
+
+pars = {'algorithm': LLT_model , \
+ 'input' : u0,
+ 'regularization_parameter': 25,\
+ 'time_step':0.0003, \
+ 'number_of_iterations' :300,\
+ 'tolerance_constant':0.001,\
+ 'restrictive_Z_smoothing': 0
+}
+out = LLT_model(pars['input'],
+ pars['regularization_parameter'],
+ pars['time_step'] ,
+ pars['number_of_iterations'],
+ pars['tolerance_constant'],
+ pars['restrictive_Z_smoothing'] )
+
+llt = out[0]
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(2,3,4)
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
+ verticalalignment='top', bbox=props)
+imgplot = plt.imshow(llt,\
+ #cmap="gray"
+ )
+
+
+# ###################### PatchBased_Regul #########################################
+# # Quick 2D denoising example in Matlab:
+# # Im = double(imread('lena_gray_256.tif'))/255; % loading image
+# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+# # ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05);
+start_time = timeit.default_timer()
+
+pars = {'algorithm': PatchBased_Regul , \
+ 'input' : u0,
+ 'regularization_parameter': 0.05,\
+ 'searching_window_ratio':3, \
+ 'similarity_window_ratio':1,\
+ 'PB_filtering_parameter': 0.08
+}
+out = PatchBased_Regul(pars['input'],
+ pars['regularization_parameter'],
+ pars['searching_window_ratio'] ,
+ pars['similarity_window_ratio'] ,
+ pars['PB_filtering_parameter'])
+pbr = out[0]
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+
+a=fig.add_subplot(2,3,5)
+
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
+ verticalalignment='top', bbox=props)
+imgplot = plt.imshow(pbr #,cmap="gray"
+ )
+
+
+# ###################### TGV_PD #########################################
+# # Quick 2D denoising example in Matlab:
+# # Im = double(imread('lena_gray_256.tif'))/255; % loading image
+# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+# # u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550);
+
+start_time = timeit.default_timer()
+
+pars = {'algorithm': TGV_PD , \
+ 'input' : u0,\
+ 'regularization_parameter':0.05,\
+ 'first_order_term': 1.3,\
+ 'second_order_term': 1, \
+ 'number_of_iterations': 550
+ }
+out = TGV_PD(pars['input'],
+ pars['regularization_parameter'],
+ pars['first_order_term'] ,
+ pars['second_order_term'] ,
+ pars['number_of_iterations'])
+tgv = out[0]
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(2,3,6)
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
+ verticalalignment='top', bbox=props)
+imgplot = plt.imshow(tgv #, cmap="gray")
+ )
+
+
+plt.show()
+
+################################################################################
+##
+## 3D Regularizers
+##
+################################################################################
+##Example:
+## figure;
+## Im = double(imread('lena_gray_256.tif'))/255; % loading image
+## u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
+## u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+#
+##filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Reconstruction\python\test\reconstruction_example.mha"
+#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Simpleflex\data\head.mha"
+#
+#reader = vtk.vtkMetaImageReader()
+#reader.SetFileName(os.path.normpath(filename))
+#reader.Update()
+##vtk returns 3D images, let's take just the one slice there is as 2D
+#Im = Converter.vtk2numpy(reader.GetOutput())
+#Im = Im.astype('float32')
+##imgplot = plt.imshow(Im)
+#perc = 0.05
+#u0 = Im + (perc* np.random.normal(size=np.shape(Im)))
+## map the u0 u0->u0>0
+#f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
+#u0 = f(u0).astype('float32')
+#converter = Converter.numpy2vtkImporter(u0, reader.GetOutput().GetSpacing(),
+# reader.GetOutput().GetOrigin())
+#converter.Update()
+#writer = vtk.vtkMetaImageWriter()
+#writer.SetInputData(converter.GetOutput())
+#writer.SetFileName(r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\noisy_head.mha")
+##writer.Write()
+#
+#
+### plot
+#fig3D = plt.figure()
+##a=fig.add_subplot(3,3,1)
+##a.set_title('Original')
+##imgplot = plt.imshow(Im)
+#sliceNo = 32
+#
+#a=fig3D.add_subplot(2,3,1)
+#a.set_title('noise')
+#imgplot = plt.imshow(u0.T[sliceNo])
+#
+#reg_output3d = []
+#
+###############################################################################
+## Call regularizer
+#
+######################## SplitBregman_TV #####################################
+## u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+#
+##reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV)
+#
+##out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30,
+## #tolerance_constant=1e-4,
+## TV_Penalty=Regularizer.TotalVariationPenalty.l1)
+#
+#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30,
+# tolerance_constant=1e-4,
+# TV_Penalty=Regularizer.TotalVariationPenalty.l1)
+#
+#
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,3,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+# verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+#
+####################### FGP_TV #########################################
+## u = FGP_TV(single(u0), 0.05, 100, 1e-04);
+#out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.005,
+# number_of_iterations=200)
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,3,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+# verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+#
+####################### LLT_model #########################################
+## * u0 = Im + .03*randn(size(Im)); % adding noise
+## [Den] = LLT_model(single(u0), 10, 0.1, 1);
+##Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0);
+##input, regularization_parameter , time_step, number_of_iterations,
+## tolerance_constant, restrictive_Z_smoothing=0
+#out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25,
+# time_step=0.0003,
+# tolerance_constant=0.0001,
+# number_of_iterations=300)
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,3,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+# verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+#
+####################### PatchBased_Regul #########################################
+## Quick 2D denoising example in Matlab:
+## Im = double(imread('lena_gray_256.tif'))/255; % loading image
+## u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+## ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05);
+#
+#out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05,
+# searching_window_ratio=3,
+# similarity_window_ratio=1,
+# PB_filtering_parameter=0.08)
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,3,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+# verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+#
+
+###################### TGV_PD #########################################
+# Quick 2D denoising example in Matlab:
+# Im = double(imread('lena_gray_256.tif'))/255; % loading image
+# u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+# u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550);
+
+
+#out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05,
+# first_order_term=1.3,
+# second_order_term=1,
+# number_of_iterations=550)
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,3,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+# verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
diff --git a/Wrappers/Python/test/test_cpu_regularizers.py b/Wrappers/Python/test/test_cpu_regularizers.py
deleted file mode 100644
index 9713baa..0000000
--- a/Wrappers/Python/test/test_cpu_regularizers.py
+++ /dev/null
@@ -1,442 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Fri Aug 4 11:10:05 2017
-
-@author: ofn77899
-"""
-
-
-import matplotlib.pyplot as plt
-import numpy as np
-import os
-from enum import Enum
-import timeit
-from ccpi.filters.cpu_regularizers_boost import SplitBregman_TV , FGP_TV ,\
- LLT_model, PatchBased_Regul ,\
- TGV_PD
-
-###############################################################################
-#https://stackoverflow.com/questions/13875989/comparing-image-in-url-to-image-in-filesystem-in-python/13884956#13884956
-#NRMSE a normalization of the root of the mean squared error
-#NRMSE is simply 1 - [RMSE / (maxval - minval)]. Where maxval is the maximum
-# intensity from the two images being compared, and respectively the same for
-# minval. RMSE is given by the square root of MSE:
-# sqrt[(sum(A - B) ** 2) / |A|],
-# where |A| means the number of elements in A. By doing this, the maximum value
-# given by RMSE is maxval.
-
-def nrmse(im1, im2):
- a, b = im1.shape
- rmse = np.sqrt(np.sum((im2 - im1) ** 2) / float(a * b))
- max_val = max(np.max(im1), np.max(im2))
- min_val = min(np.min(im1), np.min(im2))
- return 1 - (rmse / (max_val - min_val))
-###############################################################################
-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))
- else:
- txt += "{0} = {1}".format(key, value)
- txt += '\n'
- return txt
-###############################################################################
-#
-# 2D Regularizers
-#
-###############################################################################
-#Example:
-# figure;
-# Im = double(imread('lena_gray_256.tif'))/255; % loading image
-# u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
-# u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
-
-# assumes the script is launched from the test directory
-filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif")
-#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\lena_gray_512.tif"
-#filename = r"/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/data/lena_gray_512.tif"
-#filename = r'/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif'
-
-Im = plt.imread(filename)
-Im = np.asarray(Im, dtype='float32')
-
-perc = 0.15
-u0 = Im + np.random.normal(loc = Im ,
- scale = perc * Im ,
- size = np.shape(Im))
-# map the u0 u0->u0>0
-f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
-u0 = f(u0).astype('float32')
-
-## plot
-fig = plt.figure()
-
-a=fig.add_subplot(2,3,1)
-a.set_title('noise')
-imgplot = plt.imshow(u0#,cmap="gray"
- )
-
-reg_output = []
-##############################################################################
-# Call regularizer
-
-####################### SplitBregman_TV #####################################
-# u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
-
-start_time = timeit.default_timer()
-pars = {'algorithm' : SplitBregman_TV , \
- 'input' : u0,
- 'regularization_parameter':10. , \
-'number_of_iterations' :35 ,\
-'tolerance_constant':0.0001 , \
-'TV_penalty': 0
-}
-
-out = SplitBregman_TV (pars['input'], pars['regularization_parameter'],
- pars['number_of_iterations'],
- pars['tolerance_constant'],
- pars['TV_penalty'])
-splitbregman = out[0]
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-
-
-a=fig.add_subplot(2,3,2)
-
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow(splitbregman,\
- #cmap="gray"
- )
-
-###################### FGP_TV #########################################
-# u = FGP_TV(single(u0), 0.05, 100, 1e-04);
-start_time = timeit.default_timer()
-pars = {'algorithm' : FGP_TV , \
- 'input' : u0,
- 'regularization_parameter':5e-4, \
- 'number_of_iterations' :10 ,\
- 'tolerance_constant':0.001,\
- 'TV_penalty': 0
-}
-
-out = FGP_TV (pars['input'],
- pars['regularization_parameter'],
- pars['number_of_iterations'],
- pars['tolerance_constant'],
- pars['TV_penalty'])
-
-fgp = out[0]
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-
-
-a=fig.add_subplot(2,3,3)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-imgplot = plt.imshow(fgp, \
- #cmap="gray"
- )
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
- verticalalignment='top', bbox=props)
-
-###################### LLT_model #########################################
-
-start_time = timeit.default_timer()
-
-pars = {'algorithm': LLT_model , \
- 'input' : u0,
- 'regularization_parameter': 25,\
- 'time_step':0.0003, \
- 'number_of_iterations' :300,\
- 'tolerance_constant':0.001,\
- 'restrictive_Z_smoothing': 0
-}
-out = LLT_model(pars['input'],
- pars['regularization_parameter'],
- pars['time_step'] ,
- pars['number_of_iterations'],
- pars['tolerance_constant'],
- pars['restrictive_Z_smoothing'] )
-
-llt = out[0]
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(2,3,4)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow(llt,\
- #cmap="gray"
- )
-
-
-# ###################### PatchBased_Regul #########################################
-# # Quick 2D denoising example in Matlab:
-# # Im = double(imread('lena_gray_256.tif'))/255; % loading image
-# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
-# # ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05);
-start_time = timeit.default_timer()
-
-pars = {'algorithm': PatchBased_Regul , \
- 'input' : u0,
- 'regularization_parameter': 0.05,\
- 'searching_window_ratio':3, \
- 'similarity_window_ratio':1,\
- 'PB_filtering_parameter': 0.08
-}
-out = PatchBased_Regul(pars['input'],
- pars['regularization_parameter'],
- pars['searching_window_ratio'] ,
- pars['similarity_window_ratio'] ,
- pars['PB_filtering_parameter'])
-pbr = out[0]
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-
-a=fig.add_subplot(2,3,5)
-
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow(pbr #,cmap="gray"
- )
-
-
-# ###################### TGV_PD #########################################
-# # Quick 2D denoising example in Matlab:
-# # Im = double(imread('lena_gray_256.tif'))/255; % loading image
-# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
-# # u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550);
-
-start_time = timeit.default_timer()
-
-pars = {'algorithm': TGV_PD , \
- 'input' : u0,\
- 'regularization_parameter':0.05,\
- 'first_order_term': 1.3,\
- 'second_order_term': 1, \
- 'number_of_iterations': 550
- }
-out = TGV_PD(pars['input'],
- pars['regularization_parameter'],
- pars['first_order_term'] ,
- pars['second_order_term'] ,
- pars['number_of_iterations'])
-tgv = out[0]
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(2,3,6)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow(tgv #, cmap="gray")
- )
-
-
-plt.show()
-
-################################################################################
-##
-## 3D Regularizers
-##
-################################################################################
-##Example:
-## figure;
-## Im = double(imread('lena_gray_256.tif'))/255; % loading image
-## u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
-## u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
-#
-##filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Reconstruction\python\test\reconstruction_example.mha"
-#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Simpleflex\data\head.mha"
-#
-#reader = vtk.vtkMetaImageReader()
-#reader.SetFileName(os.path.normpath(filename))
-#reader.Update()
-##vtk returns 3D images, let's take just the one slice there is as 2D
-#Im = Converter.vtk2numpy(reader.GetOutput())
-#Im = Im.astype('float32')
-##imgplot = plt.imshow(Im)
-#perc = 0.05
-#u0 = Im + (perc* np.random.normal(size=np.shape(Im)))
-## map the u0 u0->u0>0
-#f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
-#u0 = f(u0).astype('float32')
-#converter = Converter.numpy2vtkImporter(u0, reader.GetOutput().GetSpacing(),
-# reader.GetOutput().GetOrigin())
-#converter.Update()
-#writer = vtk.vtkMetaImageWriter()
-#writer.SetInputData(converter.GetOutput())
-#writer.SetFileName(r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\noisy_head.mha")
-##writer.Write()
-#
-#
-### plot
-#fig3D = plt.figure()
-##a=fig.add_subplot(3,3,1)
-##a.set_title('Original')
-##imgplot = plt.imshow(Im)
-#sliceNo = 32
-#
-#a=fig3D.add_subplot(2,3,1)
-#a.set_title('noise')
-#imgplot = plt.imshow(u0.T[sliceNo])
-#
-#reg_output3d = []
-#
-###############################################################################
-## Call regularizer
-#
-######################## SplitBregman_TV #####################################
-## u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
-#
-##reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV)
-#
-##out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30,
-## #tolerance_constant=1e-4,
-## TV_Penalty=Regularizer.TotalVariationPenalty.l1)
-#
-#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30,
-# tolerance_constant=1e-4,
-# TV_Penalty=Regularizer.TotalVariationPenalty.l1)
-#
-#
-#pars = out2[-2]
-#reg_output3d.append(out2)
-#
-#a=fig3D.add_subplot(2,3,2)
-#
-#
-#textstr = out2[-1]
-#
-#
-## these are matplotlib.patch.Patch properties
-#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-## place a text box in upper left in axes coords
-#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
-# verticalalignment='top', bbox=props)
-#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
-#
-####################### FGP_TV #########################################
-## u = FGP_TV(single(u0), 0.05, 100, 1e-04);
-#out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.005,
-# number_of_iterations=200)
-#pars = out2[-2]
-#reg_output3d.append(out2)
-#
-#a=fig3D.add_subplot(2,3,2)
-#
-#
-#textstr = out2[-1]
-#
-#
-## these are matplotlib.patch.Patch properties
-#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-## place a text box in upper left in axes coords
-#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
-# verticalalignment='top', bbox=props)
-#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
-#
-####################### LLT_model #########################################
-## * u0 = Im + .03*randn(size(Im)); % adding noise
-## [Den] = LLT_model(single(u0), 10, 0.1, 1);
-##Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0);
-##input, regularization_parameter , time_step, number_of_iterations,
-## tolerance_constant, restrictive_Z_smoothing=0
-#out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25,
-# time_step=0.0003,
-# tolerance_constant=0.0001,
-# number_of_iterations=300)
-#pars = out2[-2]
-#reg_output3d.append(out2)
-#
-#a=fig3D.add_subplot(2,3,2)
-#
-#
-#textstr = out2[-1]
-#
-#
-## these are matplotlib.patch.Patch properties
-#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-## place a text box in upper left in axes coords
-#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
-# verticalalignment='top', bbox=props)
-#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
-#
-####################### PatchBased_Regul #########################################
-## Quick 2D denoising example in Matlab:
-## Im = double(imread('lena_gray_256.tif'))/255; % loading image
-## u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
-## ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05);
-#
-#out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05,
-# searching_window_ratio=3,
-# similarity_window_ratio=1,
-# PB_filtering_parameter=0.08)
-#pars = out2[-2]
-#reg_output3d.append(out2)
-#
-#a=fig3D.add_subplot(2,3,2)
-#
-#
-#textstr = out2[-1]
-#
-#
-## these are matplotlib.patch.Patch properties
-#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-## place a text box in upper left in axes coords
-#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
-# verticalalignment='top', bbox=props)
-#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
-#
-
-###################### TGV_PD #########################################
-# Quick 2D denoising example in Matlab:
-# Im = double(imread('lena_gray_256.tif'))/255; % loading image
-# u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
-# u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550);
-
-
-#out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05,
-# first_order_term=1.3,
-# second_order_term=1,
-# number_of_iterations=550)
-#pars = out2[-2]
-#reg_output3d.append(out2)
-#
-#a=fig3D.add_subplot(2,3,2)
-#
-#
-#textstr = out2[-1]
-#
-#
-## these are matplotlib.patch.Patch properties
-#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-## place a text box in upper left in axes coords
-#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
-# verticalalignment='top', bbox=props)
-#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
diff --git a/Wrappers/Python/test/test_cpu_vs_gpu_regularisers.py b/Wrappers/Python/test/test_cpu_vs_gpu_regularisers.py
new file mode 100644
index 0000000..15e9042
--- /dev/null
+++ b/Wrappers/Python/test/test_cpu_vs_gpu_regularisers.py
@@ -0,0 +1,219 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Feb 22 11:39:43 2018
+
+Testing CPU implementation against the GPU one
+
+@author: Daniil Kazantsev
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+import timeit
+from ccpi.filters.regularisers import ROF_TV, FGP_TV
+
+###############################################################################
+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))
+ else:
+ txt += "{0} = {1}".format(key, value)
+ txt += '\n'
+ return txt
+###############################################################################
+def rmse(im1, im2):
+ a, b = im1.shape
+ rmse = np.sqrt(np.sum((im1 - im2) ** 2) / float(a * b))
+ return rmse
+
+filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif")
+
+# read image
+Im = plt.imread(filename)
+Im = np.asarray(Im, dtype='float32')
+
+Im = Im/255
+perc = 0.075
+u0 = Im + np.random.normal(loc = Im ,
+ scale = perc * Im ,
+ size = np.shape(Im))
+# map the u0 u0->u0>0
+f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
+u0 = f(u0).astype('float32')
+
+
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+print ("____________ROF-TV bench___________________")
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+
+## plot
+fig = plt.figure(1)
+plt.suptitle('Comparison of ROF-TV regulariser using CPU and GPU implementations')
+a=fig.add_subplot(1,4,1)
+a.set_title('Noisy Image')
+imgplot = plt.imshow(u0,cmap="gray")
+
+# set parameters
+pars = {'algorithm': ROF_TV, \
+ 'input' : u0,\
+ 'regularisation_parameter':0.04,\
+ 'number_of_iterations': 1200,\
+ 'time_marching_parameter': 0.0025
+ }
+print ("#############ROF TV CPU####################")
+start_time = timeit.default_timer()
+rof_cpu = ROF_TV(pars['input'],
+ pars['regularisation_parameter'],
+ pars['number_of_iterations'],
+ pars['time_marching_parameter'],'cpu')
+rms = rmse(Im, rof_cpu)
+pars['rmse'] = rms
+
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(1,4,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(rof_cpu, cmap="gray")
+plt.title('{}'.format('CPU results'))
+
+
+print ("##############ROF TV GPU##################")
+start_time = timeit.default_timer()
+rof_gpu = ROF_TV(pars['input'],
+ pars['regularisation_parameter'],
+ pars['number_of_iterations'],
+ pars['time_marching_parameter'],'gpu')
+
+rms = rmse(Im, rof_gpu)
+pars['rmse'] = rms
+pars['algorithm'] = ROF_TV
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(1,4,3)
+
+# 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(rof_gpu, cmap="gray")
+plt.title('{}'.format('GPU results'))
+
+
+print ("--------Compare the results--------")
+tolerance = 1e-05
+diff_im = np.zeros(np.shape(rof_cpu))
+diff_im = abs(rof_cpu - rof_gpu)
+diff_im[diff_im > tolerance] = 1
+a=fig.add_subplot(1,4,4)
+imgplot = plt.imshow(diff_im, vmin=0, vmax=1, cmap="gray")
+plt.title('{}'.format('Pixels larger threshold difference'))
+if (diff_im.sum() > 1):
+ print ("Arrays do not match!")
+else:
+ print ("Arrays match")
+
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+print ("____________FGP-TV bench___________________")
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+
+## plot
+fig = plt.figure(2)
+plt.suptitle('Comparison of FGP-TV regulariser using CPU and GPU implementations')
+a=fig.add_subplot(1,4,1)
+a.set_title('Noisy Image')
+imgplot = plt.imshow(u0,cmap="gray")
+
+# set parameters
+pars = {'algorithm' : FGP_TV, \
+ 'input' : u0,\
+ 'regularisation_parameter':0.04, \
+ 'number_of_iterations' :1200 ,\
+ 'tolerance_constant':0.00001,\
+ 'methodTV': 0 ,\
+ 'nonneg': 0 ,\
+ 'printingOut': 0
+ }
+
+print ("#############FGP TV CPU####################")
+start_time = timeit.default_timer()
+fgp_cpu = FGP_TV(pars['input'],
+ pars['regularisation_parameter'],
+ pars['number_of_iterations'],
+ pars['tolerance_constant'],
+ pars['methodTV'],
+ pars['nonneg'],
+ pars['printingOut'],'cpu')
+
+
+rms = rmse(Im, fgp_cpu)
+pars['rmse'] = rms
+
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(1,4,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(fgp_cpu, cmap="gray")
+plt.title('{}'.format('CPU results'))
+
+
+print ("##############FGP TV GPU##################")
+start_time = timeit.default_timer()
+fgp_gpu = FGP_TV(pars['input'],
+ pars['regularisation_parameter'],
+ pars['number_of_iterations'],
+ pars['tolerance_constant'],
+ pars['methodTV'],
+ pars['nonneg'],
+ pars['printingOut'],'gpu')
+
+rms = rmse(Im, fgp_gpu)
+pars['rmse'] = rms
+pars['algorithm'] = FGP_TV
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(1,4,3)
+
+# 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(fgp_gpu, cmap="gray")
+plt.title('{}'.format('GPU results'))
+
+
+print ("--------Compare the results--------")
+tolerance = 1e-05
+diff_im = np.zeros(np.shape(rof_cpu))
+diff_im = abs(fgp_cpu - fgp_gpu)
+diff_im[diff_im > tolerance] = 1
+a=fig.add_subplot(1,4,4)
+imgplot = plt.imshow(diff_im, vmin=0, vmax=1, cmap="gray")
+plt.title('{}'.format('Pixels larger threshold difference'))
+if (diff_im.sum() > 1):
+ print ("Arrays do not match!")
+else:
+ print ("Arrays match")
+
+
diff --git a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py
deleted file mode 100644
index 63be1a0..0000000
--- a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py
+++ /dev/null
@@ -1,219 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Thu Feb 22 11:39:43 2018
-
-Testing CPU implementation against the GPU one
-
-@author: Daniil Kazantsev
-"""
-
-import matplotlib.pyplot as plt
-import numpy as np
-import os
-import timeit
-from ccpi.filters.regularizers import ROF_TV, FGP_TV
-
-###############################################################################
-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))
- else:
- txt += "{0} = {1}".format(key, value)
- txt += '\n'
- return txt
-###############################################################################
-def rmse(im1, im2):
- a, b = im1.shape
- rmse = np.sqrt(np.sum((im1 - im2) ** 2) / float(a * b))
- return rmse
-
-filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif")
-
-# read image
-Im = plt.imread(filename)
-Im = np.asarray(Im, dtype='float32')
-
-Im = Im/255
-perc = 0.075
-u0 = Im + np.random.normal(loc = Im ,
- scale = perc * Im ,
- size = np.shape(Im))
-# map the u0 u0->u0>0
-f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
-u0 = f(u0).astype('float32')
-
-
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("____________ROF-TV bench___________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot
-fig = plt.figure(1)
-plt.suptitle('Comparison of ROF-TV regularizer using CPU and GPU implementations')
-a=fig.add_subplot(1,4,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm': ROF_TV, \
- 'input' : u0,\
- 'regularization_parameter':0.04,\
- 'number_of_iterations': 1200,\
- 'time_marching_parameter': 0.0025
- }
-print ("#############ROF TV CPU####################")
-start_time = timeit.default_timer()
-rof_cpu = ROF_TV(pars['input'],
- pars['regularization_parameter'],
- pars['number_of_iterations'],
- pars['time_marching_parameter'],'cpu')
-rms = rmse(Im, rof_cpu)
-pars['rmse'] = rms
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,4,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(rof_cpu, cmap="gray")
-plt.title('{}'.format('CPU results'))
-
-
-print ("##############ROF TV GPU##################")
-start_time = timeit.default_timer()
-rof_gpu = ROF_TV(pars['input'],
- pars['regularization_parameter'],
- pars['number_of_iterations'],
- pars['time_marching_parameter'],'gpu')
-
-rms = rmse(Im, rof_gpu)
-pars['rmse'] = rms
-pars['algorithm'] = ROF_TV
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,4,3)
-
-# 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(rof_gpu, cmap="gray")
-plt.title('{}'.format('GPU results'))
-
-
-print ("--------Compare the results--------")
-tolerance = 1e-05
-diff_im = np.zeros(np.shape(rof_cpu))
-diff_im = abs(rof_cpu - rof_gpu)
-diff_im[diff_im > tolerance] = 1
-a=fig.add_subplot(1,4,4)
-imgplot = plt.imshow(diff_im, vmin=0, vmax=1, cmap="gray")
-plt.title('{}'.format('Pixels larger threshold difference'))
-if (diff_im.sum() > 1):
- print ("Arrays do not match!")
-else:
- print ("Arrays match")
-
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("____________FGP-TV bench___________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot
-fig = plt.figure(2)
-plt.suptitle('Comparison of FGP-TV regularizer using CPU and GPU implementations')
-a=fig.add_subplot(1,4,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm' : FGP_TV, \
- 'input' : u0,\
- 'regularization_parameter':0.04, \
- 'number_of_iterations' :1200 ,\
- 'tolerance_constant':0.00001,\
- 'methodTV': 0 ,\
- 'nonneg': 0 ,\
- 'printingOut': 0
- }
-
-print ("#############FGP TV CPU####################")
-start_time = timeit.default_timer()
-fgp_cpu = FGP_TV(pars['input'],
- pars['regularization_parameter'],
- pars['number_of_iterations'],
- pars['tolerance_constant'],
- pars['methodTV'],
- pars['nonneg'],
- pars['printingOut'],'cpu')
-
-
-rms = rmse(Im, fgp_cpu)
-pars['rmse'] = rms
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,4,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(fgp_cpu, cmap="gray")
-plt.title('{}'.format('CPU results'))
-
-
-print ("##############FGP TV GPU##################")
-start_time = timeit.default_timer()
-fgp_gpu = FGP_TV(pars['input'],
- pars['regularization_parameter'],
- pars['number_of_iterations'],
- pars['tolerance_constant'],
- pars['methodTV'],
- pars['nonneg'],
- pars['printingOut'],'gpu')
-
-rms = rmse(Im, fgp_gpu)
-pars['rmse'] = rms
-pars['algorithm'] = FGP_TV
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,4,3)
-
-# 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(fgp_gpu, cmap="gray")
-plt.title('{}'.format('GPU results'))
-
-
-print ("--------Compare the results--------")
-tolerance = 1e-05
-diff_im = np.zeros(np.shape(rof_cpu))
-diff_im = abs(fgp_cpu - fgp_gpu)
-diff_im[diff_im > tolerance] = 1
-a=fig.add_subplot(1,4,4)
-imgplot = plt.imshow(diff_im, vmin=0, vmax=1, cmap="gray")
-plt.title('{}'.format('Pixels larger threshold difference'))
-if (diff_im.sum() > 1):
- print ("Arrays do not match!")
-else:
- print ("Arrays match")
-
-
diff --git a/Wrappers/Python/test/test_gpu_regularisers.py b/Wrappers/Python/test/test_gpu_regularisers.py
new file mode 100644
index 0000000..2103c0e
--- /dev/null
+++ b/Wrappers/Python/test/test_gpu_regularisers.py
@@ -0,0 +1,143 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Jan 30 10:24:26 2018
+
+@author: ofn77899
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+from enum import Enum
+import timeit
+from ccpi.filters.regularisers import ROF_TV, FGP_TV
+###############################################################################
+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))
+ else:
+ txt += "{0} = {1}".format(key, value)
+ txt += '\n'
+ return txt
+###############################################################################
+def rmse(im1, im2):
+ a, b = im1.shape
+ rmse = np.sqrt(np.sum((im1 - im2) ** 2) / float(a * b))
+ return rmse
+
+filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif")
+
+Im = plt.imread(filename)
+Im = np.asarray(Im, dtype='float32')
+
+Im = Im/255
+perc = 0.075
+u0 = Im + np.random.normal(loc = Im ,
+ scale = perc * Im ,
+ size = np.shape(Im))
+# map the u0 u0->u0>0
+f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
+u0 = f(u0).astype('float32')
+
+## plot
+fig = plt.figure()
+
+a=fig.add_subplot(2,4,1)
+a.set_title('noise')
+imgplot = plt.imshow(u0,cmap="gray")
+
+
+
+## Rudin-Osher-Fatemi (ROF) TV regularisation
+start_time = timeit.default_timer()
+pars = {
+'algorithm' : ROF_TV , \
+ 'input' : u0,
+ 'regularisation_parameter': 0.04,\
+ 'number_of_iterations':300,\
+ 'time_marching_parameter': 0.0025
+
+ }
+
+rof_tv = TV_ROF_GPU(pars['input'],
+ pars['regularisation_parameter'],
+ pars['number_of_iterations'],
+ pars['time_marching_parameter'],'gpu')
+
+rms = rmse(Im, rof_tv)
+pars['rmse'] = rms
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(2,4,4)
+
+# 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=12,
+ verticalalignment='top', bbox=props)
+imgplot = plt.imshow(rof_tv, cmap="gray")
+
+a=fig.add_subplot(2,4,8)
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, 'rof_tv - u0', transform=a.transAxes, fontsize=14,
+ verticalalignment='top', bbox=props)
+imgplot = plt.imshow((rof_tv - u0)**2, vmin=0, vmax=0.03, cmap="gray")
+plt.colorbar(ticks=[0, 0.03], orientation='vertical')
+plt.show()
+
+## Fast-Gradient Projection TV regularisation
+"""
+start_time = timeit.default_timer()
+
+pars = {'algorithm' : FGP_TV, \
+ 'input' : u0,\
+ 'regularisation_parameter':0.04, \
+ 'number_of_iterations' :1200 ,\
+ 'tolerance_constant':0.00001,\
+ 'methodTV': 0 ,\
+ 'nonneg': 0 ,\
+ 'printingOut': 0
+ }
+
+fgp_gpu = FGP_TV(pars['input'],
+ pars['regularisation_parameter'],
+ pars['number_of_iterations'],
+ pars['tolerance_constant'],
+ pars['methodTV'],
+ pars['nonneg'],
+ pars['printingOut'],'gpu')
+
+rms = rmse(Im, fgp_gpu)
+pars['rmse'] = rms
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(2,4,4)
+
+# 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=12,
+ verticalalignment='top', bbox=props)
+imgplot = plt.imshow(fgp_gpu, cmap="gray")
+
+a=fig.add_subplot(2,4,8)
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, 'fgp_gpu - u0', transform=a.transAxes, fontsize=14,
+ verticalalignment='top', bbox=props)
+imgplot = plt.imshow((fgp_gpu - u0)**2, vmin=0, vmax=0.03, cmap="gray")
+plt.colorbar(ticks=[0, 0.03], orientation='vertical')
+plt.show()
+"""
diff --git a/Wrappers/Python/test/test_gpu_regularizers.py b/Wrappers/Python/test/test_gpu_regularizers.py
deleted file mode 100644
index 640b3f9..0000000
--- a/Wrappers/Python/test/test_gpu_regularizers.py
+++ /dev/null
@@ -1,239 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Tue Jan 30 10:24:26 2018
-
-@author: ofn77899
-"""
-
-import matplotlib.pyplot as plt
-import numpy as np
-import os
-from enum import Enum
-import timeit
-from ccpi.filters.gpu_regularizers import Diff4thHajiaboli, NML
-from ccpi.filters.regularizers import ROF_TV, FGP_TV
-###############################################################################
-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))
- else:
- txt += "{0} = {1}".format(key, value)
- txt += '\n'
- return txt
-###############################################################################
-def rmse(im1, im2):
- a, b = im1.shape
- rmse = np.sqrt(np.sum((im1 - im2) ** 2) / float(a * b))
- return rmse
-
-filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif")
-#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\lena_gray_512.tif"
-#filename = r"/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/data/lena_gray_512.tif"
-#filename = r'/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif'
-
-Im = plt.imread(filename)
-Im = np.asarray(Im, dtype='float32')
-
-Im = Im/255
-perc = 0.075
-u0 = Im + np.random.normal(loc = Im ,
- scale = perc * Im ,
- size = np.shape(Im))
-# map the u0 u0->u0>0
-f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
-u0 = f(u0).astype('float32')
-
-## plot
-fig = plt.figure()
-
-a=fig.add_subplot(2,4,1)
-a.set_title('noise')
-imgplot = plt.imshow(u0,cmap="gray")
-
-
-## Diff4thHajiaboli
-start_time = timeit.default_timer()
-pars = {
-'algorithm' : Diff4thHajiaboli , \
- 'input' : u0,
- 'edge_preserv_parameter':0.1 , \
-'number_of_iterations' :250 ,\
-'time_marching_parameter':0.03 ,\
-'regularization_parameter':0.7
-}
-
-
-d4h = Diff4thHajiaboli(pars['input'],
- pars['edge_preserv_parameter'],
- pars['number_of_iterations'],
- pars['time_marching_parameter'],
- pars['regularization_parameter'])
-rms = rmse(Im, d4h)
-pars['rmse'] = rms
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(2,4,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=12,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow(d4h, cmap="gray")
-
-a=fig.add_subplot(2,4,6)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, 'd4h - u0', transform=a.transAxes, fontsize=12,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow((d4h - u0)**2, vmin=0, vmax=0.03, cmap="gray")
-plt.colorbar(ticks=[0, 0.03], orientation='vertical')
-
-
-## Patch Based Regul NML
-start_time = timeit.default_timer()
-"""
-pars = {'algorithm' : NML , \
- 'input' : u0,
- 'SearchW_real':3 , \
-'SimilW' :1,\
-'h':0.05 ,#
-'lambda' : 0.08
-}
-"""
-pars = {
-'algorithm' : NML , \
- 'input' : u0,
- 'regularization_parameter': 0.01,\
- 'searching_window_ratio':3, \
- 'similarity_window_ratio':1,\
- 'PB_filtering_parameter': 0.2
-}
-
-nml = NML(pars['input'],
- pars['searching_window_ratio'],
- pars['similarity_window_ratio'],
- pars['PB_filtering_parameter'],
- pars['regularization_parameter'])
-rms = rmse(Im, nml)
-pars['rmse'] = rms
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(2,4,3)
-
-# 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=12,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow(nml, cmap="gray")
-
-a=fig.add_subplot(2,4,7)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, 'nml - u0', transform=a.transAxes, fontsize=14,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow((nml - u0)**2, vmin=0, vmax=0.03, cmap="gray")
-plt.colorbar(ticks=[0, 0.03], orientation='vertical')
-
-
-
-## Rudin-Osher-Fatemi (ROF) TV regularization
-start_time = timeit.default_timer()
-pars = {
-'algorithm' : ROF_TV , \
- 'input' : u0,
- 'regularization_parameter': 0.04,\
- 'number_of_iterations':300,\
- 'time_marching_parameter': 0.0025
-
- }
-
-rof_tv = TV_ROF_GPU(pars['input'],
- pars['regularization_parameter'],
- pars['number_of_iterations'],
- pars['time_marching_parameter'],'gpu')
-
-rms = rmse(Im, rof_tv)
-pars['rmse'] = rms
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(2,4,4)
-
-# 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=12,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow(rof_tv, cmap="gray")
-
-a=fig.add_subplot(2,4,8)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, 'rof_tv - u0', transform=a.transAxes, fontsize=14,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow((rof_tv - u0)**2, vmin=0, vmax=0.03, cmap="gray")
-plt.colorbar(ticks=[0, 0.03], orientation='vertical')
-plt.show()
-
-## Fast-Gradient Projection TV regularization
-"""
-start_time = timeit.default_timer()
-
-pars = {'algorithm' : FGP_TV, \
- 'input' : u0,\
- 'regularization_parameter':0.04, \
- 'number_of_iterations' :1200 ,\
- 'tolerance_constant':0.00001,\
- 'methodTV': 0 ,\
- 'nonneg': 0 ,\
- 'printingOut': 0
- }
-
-fgp_gpu = FGP_TV(pars['input'],
- pars['regularization_parameter'],
- pars['number_of_iterations'],
- pars['tolerance_constant'],
- pars['methodTV'],
- pars['nonneg'],
- pars['printingOut'],'gpu')
-
-rms = rmse(Im, fgp_gpu)
-pars['rmse'] = rms
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(2,4,4)
-
-# 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=12,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow(fgp_gpu, cmap="gray")
-
-a=fig.add_subplot(2,4,8)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, 'fgp_gpu - u0', transform=a.transAxes, fontsize=14,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow((fgp_gpu - u0)**2, vmin=0, vmax=0.03, cmap="gray")
-plt.colorbar(ticks=[0, 0.03], orientation='vertical')
-plt.show()
-"""
diff --git a/Wrappers/Python/test/test_regularisers_3d.py b/Wrappers/Python/test/test_regularisers_3d.py
new file mode 100644
index 0000000..2d11a7e
--- /dev/null
+++ b/Wrappers/Python/test/test_regularisers_3d.py
@@ -0,0 +1,425 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Aug 4 11:10:05 2017
+
+@author: ofn77899
+"""
+
+#from ccpi.viewer.CILViewer2D import Converter
+#import vtk
+
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+from enum import Enum
+import timeit
+#from PIL import Image
+#from Regularizer import Regularizer
+from ccpi.imaging.Regularizer import Regularizer
+
+###############################################################################
+#https://stackoverflow.com/questions/13875989/comparing-image-in-url-to-image-in-filesystem-in-python/13884956#13884956
+#NRMSE a normalization of the root of the mean squared error
+#NRMSE is simply 1 - [RMSE / (maxval - minval)]. Where maxval is the maximum
+# intensity from the two images being compared, and respectively the same for
+# minval. RMSE is given by the square root of MSE:
+# sqrt[(sum(A - B) ** 2) / |A|],
+# where |A| means the number of elements in A. By doing this, the maximum value
+# given by RMSE is maxval.
+
+def nrmse(im1, im2):
+ a, b = im1.shape
+ rmse = np.sqrt(np.sum((im2 - im1) ** 2) / float(a * b))
+ max_val = max(np.max(im1), np.max(im2))
+ min_val = min(np.min(im1), np.min(im2))
+ return 1 - (rmse / (max_val - min_val))
+###############################################################################
+
+###############################################################################
+#
+# 2D Regularizers
+#
+###############################################################################
+#Example:
+# figure;
+# Im = double(imread('lena_gray_256.tif'))/255; % loading image
+# u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
+# u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+
+
+#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\lena_gray_512.tif"
+filename = r"/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/data/lena_gray_512.tif"
+#filename = r'/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif'
+
+#reader = vtk.vtkTIFFReader()
+#reader.SetFileName(os.path.normpath(filename))
+#reader.Update()
+Im = plt.imread(filename)
+#Im = Image.open('/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif')/255
+#img.show()
+Im = np.asarray(Im, dtype='float32')
+
+# create a 3D image by stacking N of this images
+
+
+#imgplot = plt.imshow(Im)
+perc = 0.05
+u_n = Im + (perc* np.random.normal(size=np.shape(Im)))
+y,z = np.shape(u_n)
+u_n = np.reshape(u_n , (1,y,z))
+
+u0 = u_n.copy()
+for i in range (19):
+ u_n = Im + (perc* np.random.normal(size=np.shape(Im)))
+ u_n = np.reshape(u_n , (1,y,z))
+
+ u0 = np.vstack ( (u0, u_n) )
+
+# map the u0 u0->u0>0
+f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
+u0 = f(u0).astype('float32')
+
+print ("Passed image shape {0}".format(np.shape(u0)))
+
+## plot
+fig = plt.figure()
+#a=fig.add_subplot(3,3,1)
+#a.set_title('Original')
+#imgplot = plt.imshow(Im)
+sliceno = 10
+
+a=fig.add_subplot(2,3,1)
+a.set_title('noise')
+imgplot = plt.imshow(u0[sliceno],cmap="gray")
+
+reg_output = []
+##############################################################################
+# Call regularizer
+
+####################### SplitBregman_TV #####################################
+# u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+
+use_object = True
+if use_object:
+ reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV)
+ print (reg.pars)
+ reg.setParameter(input=u0)
+ reg.setParameter(regularization_parameter=10.)
+ # or
+ # reg.setParameter(input=u0, regularization_parameter=10., #number_of_iterations=30,
+ #tolerance_constant=1e-4,
+ #TV_Penalty=Regularizer.TotalVariationPenalty.l1)
+ plotme = reg() [0]
+ pars = reg.pars
+ textstr = reg.printParametersToString()
+
+ #out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30,
+ #tolerance_constant=1e-4,
+ # TV_Penalty=Regularizer.TotalVariationPenalty.l1)
+
+#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30,
+# tolerance_constant=1e-4,
+# TV_Penalty=Regularizer.TotalVariationPenalty.l1)
+
+else:
+ out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10. )
+ pars = out2[2]
+ reg_output.append(out2)
+ plotme = reg_output[-1][0]
+ textstr = out2[-1]
+
+a=fig.add_subplot(2,3,2)
+
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+ verticalalignment='top', bbox=props)
+imgplot = plt.imshow(plotme[sliceno],cmap="gray")
+
+###################### FGP_TV #########################################
+# u = FGP_TV(single(u0), 0.05, 100, 1e-04);
+out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.0005,
+ number_of_iterations=50)
+pars = out2[-2]
+
+reg_output.append(out2)
+
+a=fig.add_subplot(2,3,3)
+
+textstr = out2[-1]
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+ verticalalignment='top', bbox=props)
+imgplot = plt.imshow(reg_output[-1][0][sliceno])
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+ verticalalignment='top', bbox=props)
+imgplot = plt.imshow(reg_output[-1][0][sliceno],cmap="gray")
+
+###################### LLT_model #########################################
+# * u0 = Im + .03*randn(size(Im)); % adding noise
+# [Den] = LLT_model(single(u0), 10, 0.1, 1);
+#Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0);
+#input, regularization_parameter , time_step, number_of_iterations,
+# tolerance_constant, restrictive_Z_smoothing=0
+out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25,
+ time_step=0.0003,
+ tolerance_constant=0.0001,
+ number_of_iterations=300)
+pars = out2[-2]
+
+reg_output.append(out2)
+
+a=fig.add_subplot(2,3,4)
+
+textstr = out2[-1]
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+ verticalalignment='top', bbox=props)
+imgplot = plt.imshow(reg_output[-1][0][sliceno],cmap="gray")
+
+
+# ###################### PatchBased_Regul #########################################
+# # Quick 2D denoising example in Matlab:
+# # Im = double(imread('lena_gray_256.tif'))/255; % loading image
+# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+# # ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05);
+
+out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05,
+ searching_window_ratio=3,
+ similarity_window_ratio=1,
+ PB_filtering_parameter=0.08)
+pars = out2[-2]
+reg_output.append(out2)
+
+a=fig.add_subplot(2,3,5)
+
+
+textstr = out2[-1]
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+ verticalalignment='top', bbox=props)
+imgplot = plt.imshow(reg_output[-1][0][sliceno],cmap="gray")
+
+
+# ###################### TGV_PD #########################################
+# # Quick 2D denoising example in Matlab:
+# # Im = double(imread('lena_gray_256.tif'))/255; % loading image
+# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+# # u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550);
+
+
+out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05,
+ first_order_term=1.3,
+ second_order_term=1,
+ number_of_iterations=550)
+pars = out2[-2]
+reg_output.append(out2)
+
+a=fig.add_subplot(2,3,6)
+
+
+textstr = out2[-1]
+
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+ verticalalignment='top', bbox=props)
+imgplot = plt.imshow(reg_output[-1][0][sliceno],cmap="gray")
+
+
+plt.show()
+
+################################################################################
+##
+## 3D Regularizers
+##
+################################################################################
+##Example:
+## figure;
+## Im = double(imread('lena_gray_256.tif'))/255; % loading image
+## u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
+## u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+#
+##filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Reconstruction\python\test\reconstruction_example.mha"
+#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Simpleflex\data\head.mha"
+#
+#reader = vtk.vtkMetaImageReader()
+#reader.SetFileName(os.path.normpath(filename))
+#reader.Update()
+##vtk returns 3D images, let's take just the one slice there is as 2D
+#Im = Converter.vtk2numpy(reader.GetOutput())
+#Im = Im.astype('float32')
+##imgplot = plt.imshow(Im)
+#perc = 0.05
+#u0 = Im + (perc* np.random.normal(size=np.shape(Im)))
+## map the u0 u0->u0>0
+#f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
+#u0 = f(u0).astype('float32')
+#converter = Converter.numpy2vtkImporter(u0, reader.GetOutput().GetSpacing(),
+# reader.GetOutput().GetOrigin())
+#converter.Update()
+#writer = vtk.vtkMetaImageWriter()
+#writer.SetInputData(converter.GetOutput())
+#writer.SetFileName(r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\noisy_head.mha")
+##writer.Write()
+#
+#
+### plot
+#fig3D = plt.figure()
+##a=fig.add_subplot(3,3,1)
+##a.set_title('Original')
+##imgplot = plt.imshow(Im)
+#sliceNo = 32
+#
+#a=fig3D.add_subplot(2,3,1)
+#a.set_title('noise')
+#imgplot = plt.imshow(u0.T[sliceNo])
+#
+#reg_output3d = []
+#
+###############################################################################
+## Call regularizer
+#
+######################## SplitBregman_TV #####################################
+## u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+#
+##reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV)
+#
+##out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30,
+## #tolerance_constant=1e-4,
+## TV_Penalty=Regularizer.TotalVariationPenalty.l1)
+#
+#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30,
+# tolerance_constant=1e-4,
+# TV_Penalty=Regularizer.TotalVariationPenalty.l1)
+#
+#
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,3,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+# verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+#
+####################### FGP_TV #########################################
+## u = FGP_TV(single(u0), 0.05, 100, 1e-04);
+#out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.005,
+# number_of_iterations=200)
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,3,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+# verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+#
+####################### LLT_model #########################################
+## * u0 = Im + .03*randn(size(Im)); % adding noise
+## [Den] = LLT_model(single(u0), 10, 0.1, 1);
+##Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0);
+##input, regularization_parameter , time_step, number_of_iterations,
+## tolerance_constant, restrictive_Z_smoothing=0
+#out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25,
+# time_step=0.0003,
+# tolerance_constant=0.0001,
+# number_of_iterations=300)
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,3,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+# verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+#
+####################### PatchBased_Regul #########################################
+## Quick 2D denoising example in Matlab:
+## Im = double(imread('lena_gray_256.tif'))/255; % loading image
+## u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+## ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05);
+#
+#out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05,
+# searching_window_ratio=3,
+# similarity_window_ratio=1,
+# PB_filtering_parameter=0.08)
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,3,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+# verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+#
+
+###################### TGV_PD #########################################
+# Quick 2D denoising example in Matlab:
+# Im = double(imread('lena_gray_256.tif'))/255; % loading image
+# u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+# u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550);
+
+
+#out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05,
+# first_order_term=1.3,
+# second_order_term=1,
+# number_of_iterations=550)
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,3,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+# verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
diff --git a/Wrappers/Python/test/test_regularizers.py b/Wrappers/Python/test/test_regularizers.py
deleted file mode 100644
index cf5da2b..0000000
--- a/Wrappers/Python/test/test_regularizers.py
+++ /dev/null
@@ -1,395 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Fri Aug 4 11:10:05 2017
-
-@author: ofn77899
-"""
-
-#from ccpi.viewer.CILViewer2D import Converter
-#import vtk
-
-import matplotlib.pyplot as plt
-import numpy as np
-import os
-from enum import Enum
-import timeit
-#from PIL import Image
-#from Regularizer import Regularizer
-from ccpi.filters.Regularizer import Regularizer
-
-###############################################################################
-#https://stackoverflow.com/questions/13875989/comparing-image-in-url-to-image-in-filesystem-in-python/13884956#13884956
-#NRMSE a normalization of the root of the mean squared error
-#NRMSE is simply 1 - [RMSE / (maxval - minval)]. Where maxval is the maximum
-# intensity from the two images being compared, and respectively the same for
-# minval. RMSE is given by the square root of MSE:
-# sqrt[(sum(A - B) ** 2) / |A|],
-# where |A| means the number of elements in A. By doing this, the maximum value
-# given by RMSE is maxval.
-
-def nrmse(im1, im2):
- a, b = im1.shape
- rmse = np.sqrt(np.sum((im2 - im1) ** 2) / float(a * b))
- max_val = max(np.max(im1), np.max(im2))
- min_val = min(np.min(im1), np.min(im2))
- return 1 - (rmse / (max_val - min_val))
-###############################################################################
-
-###############################################################################
-#
-# 2D Regularizers
-#
-###############################################################################
-#Example:
-# figure;
-# Im = double(imread('lena_gray_256.tif'))/255; % loading image
-# u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
-# u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
-
-
-filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\lena_gray_512.tif"
-#filename = r"/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/data/lena_gray_512.tif"
-#filename = r'/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif'
-
-#reader = vtk.vtkTIFFReader()
-#reader.SetFileName(os.path.normpath(filename))
-#reader.Update()
-Im = plt.imread(filename)
-#Im = Image.open('/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif')/255
-#img.show()
-Im = np.asarray(Im, dtype='float32')
-
-
-
-
-#imgplot = plt.imshow(Im)
-perc = 0.05
-u0 = Im + (perc* np.random.normal(size=np.shape(Im)))
-# map the u0 u0->u0>0
-f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
-u0 = f(u0).astype('float32')
-
-## plot
-fig = plt.figure()
-#a=fig.add_subplot(3,3,1)
-#a.set_title('Original')
-#imgplot = plt.imshow(Im)
-
-a=fig.add_subplot(2,3,1)
-a.set_title('noise')
-imgplot = plt.imshow(u0,cmap="gray")
-
-reg_output = []
-##############################################################################
-# Call regularizer
-
-####################### SplitBregman_TV #####################################
-# u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
-
-start_time = timeit.default_timer()
-reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV)
-print (reg.pars)
-reg.setParameter(input=u0)
-reg.setParameter(regularization_parameter=10.)
-# or
-# reg.setParameter(input=u0, regularization_parameter=10., #number_of_iterations=30,
- #tolerance_constant=1e-4,
- #TV_Penalty=Regularizer.TotalVariationPenalty.l1)
-plotme = reg(output_all=True) [0]
-pars = reg.pars
-txtstr = reg.printParametersToString()
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-
-
-a=fig.add_subplot(2,3,2)
-
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow(plotme,cmap="gray")
-
-###################### FGP_TV #########################################
-# u = FGP_TV(single(u0), 0.05, 100, 1e-04);
-start_time = timeit.default_timer()
-reg = Regularizer(Regularizer.Algorithm.FGP_TV)
-out2 = reg(input=u0, regularization_parameter=5e-4,
- number_of_iterations=10)
-txtstr = reg.printParametersToString()
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-
-
-a=fig.add_subplot(2,3,3)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-imgplot = plt.imshow(out2,cmap="gray")
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
- verticalalignment='top', bbox=props)
-
-###################### LLT_model #########################################
-# * u0 = Im + .03*randn(size(Im)); % adding noise
-# [Den] = LLT_model(single(u0), 10, 0.1, 1);
-#Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0);
-#input, regularization_parameter , time_step, number_of_iterations,
-# tolerance_constant, restrictive_Z_smoothing=0
-
-del out2
-start_time = timeit.default_timer()
-reg = Regularizer(Regularizer.Algorithm.LLT_model)
-out2 = reg(input=u0, regularization_parameter=25,
- time_step=0.0003,
- tolerance_constant=0.001,
- number_of_iterations=300)
-txtstr = reg.printParametersToString()
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(2,3,4)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow(out2,cmap="gray")
-
-
-# ###################### PatchBased_Regul #########################################
-# # Quick 2D denoising example in Matlab:
-# # Im = double(imread('lena_gray_256.tif'))/255; % loading image
-# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
-# # ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05);
-start_time = timeit.default_timer()
-reg = Regularizer(Regularizer.Algorithm.PatchBased_Regul)
-out2 = reg(input=u0, regularization_parameter=0.05,
- searching_window_ratio=3,
- similarity_window_ratio=1,
- PB_filtering_parameter=0.08)
-txtstr = reg.printParametersToString()
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-
-a=fig.add_subplot(2,3,5)
-
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow(out2,cmap="gray")
-
-
-# ###################### TGV_PD #########################################
-# # Quick 2D denoising example in Matlab:
-# # Im = double(imread('lena_gray_256.tif'))/255; % loading image
-# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
-# # u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550);
-
-start_time = timeit.default_timer()
-reg = Regularizer(Regularizer.Algorithm.TGV_PD)
-out2 = reg(input=u0, regularization_parameter=0.05,
- first_order_term=1.3,
- second_order_term=1,
- number_of_iterations=550)
-txtstr = reg.printParametersToString()
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(2,3,6)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow(out2,cmap="gray")
-
-
-plt.show()
-
-################################################################################
-##
-## 3D Regularizers
-##
-################################################################################
-##Example:
-## figure;
-## Im = double(imread('lena_gray_256.tif'))/255; % loading image
-## u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
-## u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
-#
-##filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Reconstruction\python\test\reconstruction_example.mha"
-#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Simpleflex\data\head.mha"
-#
-#reader = vtk.vtkMetaImageReader()
-#reader.SetFileName(os.path.normpath(filename))
-#reader.Update()
-##vtk returns 3D images, let's take just the one slice there is as 2D
-#Im = Converter.vtk2numpy(reader.GetOutput())
-#Im = Im.astype('float32')
-##imgplot = plt.imshow(Im)
-#perc = 0.05
-#u0 = Im + (perc* np.random.normal(size=np.shape(Im)))
-## map the u0 u0->u0>0
-#f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
-#u0 = f(u0).astype('float32')
-#converter = Converter.numpy2vtkImporter(u0, reader.GetOutput().GetSpacing(),
-# reader.GetOutput().GetOrigin())
-#converter.Update()
-#writer = vtk.vtkMetaImageWriter()
-#writer.SetInputData(converter.GetOutput())
-#writer.SetFileName(r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\noisy_head.mha")
-##writer.Write()
-#
-#
-### plot
-#fig3D = plt.figure()
-##a=fig.add_subplot(3,3,1)
-##a.set_title('Original')
-##imgplot = plt.imshow(Im)
-#sliceNo = 32
-#
-#a=fig3D.add_subplot(2,3,1)
-#a.set_title('noise')
-#imgplot = plt.imshow(u0.T[sliceNo])
-#
-#reg_output3d = []
-#
-###############################################################################
-## Call regularizer
-#
-######################## SplitBregman_TV #####################################
-## u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
-#
-##reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV)
-#
-##out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30,
-## #tolerance_constant=1e-4,
-## TV_Penalty=Regularizer.TotalVariationPenalty.l1)
-#
-#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30,
-# tolerance_constant=1e-4,
-# TV_Penalty=Regularizer.TotalVariationPenalty.l1)
-#
-#
-#pars = out2[-2]
-#reg_output3d.append(out2)
-#
-#a=fig3D.add_subplot(2,3,2)
-#
-#
-#textstr = out2[-1]
-#
-#
-## these are matplotlib.patch.Patch properties
-#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-## place a text box in upper left in axes coords
-#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
-# verticalalignment='top', bbox=props)
-#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
-#
-####################### FGP_TV #########################################
-## u = FGP_TV(single(u0), 0.05, 100, 1e-04);
-#out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.005,
-# number_of_iterations=200)
-#pars = out2[-2]
-#reg_output3d.append(out2)
-#
-#a=fig3D.add_subplot(2,3,2)
-#
-#
-#textstr = out2[-1]
-#
-#
-## these are matplotlib.patch.Patch properties
-#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-## place a text box in upper left in axes coords
-#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
-# verticalalignment='top', bbox=props)
-#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
-#
-####################### LLT_model #########################################
-## * u0 = Im + .03*randn(size(Im)); % adding noise
-## [Den] = LLT_model(single(u0), 10, 0.1, 1);
-##Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0);
-##input, regularization_parameter , time_step, number_of_iterations,
-## tolerance_constant, restrictive_Z_smoothing=0
-#out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25,
-# time_step=0.0003,
-# tolerance_constant=0.0001,
-# number_of_iterations=300)
-#pars = out2[-2]
-#reg_output3d.append(out2)
-#
-#a=fig3D.add_subplot(2,3,2)
-#
-#
-#textstr = out2[-1]
-#
-#
-## these are matplotlib.patch.Patch properties
-#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-## place a text box in upper left in axes coords
-#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
-# verticalalignment='top', bbox=props)
-#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
-#
-####################### PatchBased_Regul #########################################
-## Quick 2D denoising example in Matlab:
-## Im = double(imread('lena_gray_256.tif'))/255; % loading image
-## u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
-## ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05);
-#
-#out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05,
-# searching_window_ratio=3,
-# similarity_window_ratio=1,
-# PB_filtering_parameter=0.08)
-#pars = out2[-2]
-#reg_output3d.append(out2)
-#
-#a=fig3D.add_subplot(2,3,2)
-#
-#
-#textstr = out2[-1]
-#
-#
-## these are matplotlib.patch.Patch properties
-#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-## place a text box in upper left in axes coords
-#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
-# verticalalignment='top', bbox=props)
-#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
-#
-
-###################### TGV_PD #########################################
-# Quick 2D denoising example in Matlab:
-# Im = double(imread('lena_gray_256.tif'))/255; % loading image
-# u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
-# u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550);
-
-
-#out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05,
-# first_order_term=1.3,
-# second_order_term=1,
-# number_of_iterations=550)
-#pars = out2[-2]
-#reg_output3d.append(out2)
-#
-#a=fig3D.add_subplot(2,3,2)
-#
-#
-#textstr = out2[-1]
-#
-#
-## these are matplotlib.patch.Patch properties
-#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-## place a text box in upper left in axes coords
-#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
-# verticalalignment='top', bbox=props)
-#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
diff --git a/Wrappers/Python/test/test_regularizers_3d.py b/Wrappers/Python/test/test_regularizers_3d.py
deleted file mode 100644
index 2d11a7e..0000000
--- a/Wrappers/Python/test/test_regularizers_3d.py
+++ /dev/null
@@ -1,425 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Fri Aug 4 11:10:05 2017
-
-@author: ofn77899
-"""
-
-#from ccpi.viewer.CILViewer2D import Converter
-#import vtk
-
-import matplotlib.pyplot as plt
-import numpy as np
-import os
-from enum import Enum
-import timeit
-#from PIL import Image
-#from Regularizer import Regularizer
-from ccpi.imaging.Regularizer import Regularizer
-
-###############################################################################
-#https://stackoverflow.com/questions/13875989/comparing-image-in-url-to-image-in-filesystem-in-python/13884956#13884956
-#NRMSE a normalization of the root of the mean squared error
-#NRMSE is simply 1 - [RMSE / (maxval - minval)]. Where maxval is the maximum
-# intensity from the two images being compared, and respectively the same for
-# minval. RMSE is given by the square root of MSE:
-# sqrt[(sum(A - B) ** 2) / |A|],
-# where |A| means the number of elements in A. By doing this, the maximum value
-# given by RMSE is maxval.
-
-def nrmse(im1, im2):
- a, b = im1.shape
- rmse = np.sqrt(np.sum((im2 - im1) ** 2) / float(a * b))
- max_val = max(np.max(im1), np.max(im2))
- min_val = min(np.min(im1), np.min(im2))
- return 1 - (rmse / (max_val - min_val))
-###############################################################################
-
-###############################################################################
-#
-# 2D Regularizers
-#
-###############################################################################
-#Example:
-# figure;
-# Im = double(imread('lena_gray_256.tif'))/255; % loading image
-# u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
-# u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
-
-
-#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\lena_gray_512.tif"
-filename = r"/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/data/lena_gray_512.tif"
-#filename = r'/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif'
-
-#reader = vtk.vtkTIFFReader()
-#reader.SetFileName(os.path.normpath(filename))
-#reader.Update()
-Im = plt.imread(filename)
-#Im = Image.open('/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif')/255
-#img.show()
-Im = np.asarray(Im, dtype='float32')
-
-# create a 3D image by stacking N of this images
-
-
-#imgplot = plt.imshow(Im)
-perc = 0.05
-u_n = Im + (perc* np.random.normal(size=np.shape(Im)))
-y,z = np.shape(u_n)
-u_n = np.reshape(u_n , (1,y,z))
-
-u0 = u_n.copy()
-for i in range (19):
- u_n = Im + (perc* np.random.normal(size=np.shape(Im)))
- u_n = np.reshape(u_n , (1,y,z))
-
- u0 = np.vstack ( (u0, u_n) )
-
-# map the u0 u0->u0>0
-f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
-u0 = f(u0).astype('float32')
-
-print ("Passed image shape {0}".format(np.shape(u0)))
-
-## plot
-fig = plt.figure()
-#a=fig.add_subplot(3,3,1)
-#a.set_title('Original')
-#imgplot = plt.imshow(Im)
-sliceno = 10
-
-a=fig.add_subplot(2,3,1)
-a.set_title('noise')
-imgplot = plt.imshow(u0[sliceno],cmap="gray")
-
-reg_output = []
-##############################################################################
-# Call regularizer
-
-####################### SplitBregman_TV #####################################
-# u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
-
-use_object = True
-if use_object:
- reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV)
- print (reg.pars)
- reg.setParameter(input=u0)
- reg.setParameter(regularization_parameter=10.)
- # or
- # reg.setParameter(input=u0, regularization_parameter=10., #number_of_iterations=30,
- #tolerance_constant=1e-4,
- #TV_Penalty=Regularizer.TotalVariationPenalty.l1)
- plotme = reg() [0]
- pars = reg.pars
- textstr = reg.printParametersToString()
-
- #out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30,
- #tolerance_constant=1e-4,
- # TV_Penalty=Regularizer.TotalVariationPenalty.l1)
-
-#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30,
-# tolerance_constant=1e-4,
-# TV_Penalty=Regularizer.TotalVariationPenalty.l1)
-
-else:
- out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10. )
- pars = out2[2]
- reg_output.append(out2)
- plotme = reg_output[-1][0]
- textstr = out2[-1]
-
-a=fig.add_subplot(2,3,2)
-
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow(plotme[sliceno],cmap="gray")
-
-###################### FGP_TV #########################################
-# u = FGP_TV(single(u0), 0.05, 100, 1e-04);
-out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.0005,
- number_of_iterations=50)
-pars = out2[-2]
-
-reg_output.append(out2)
-
-a=fig.add_subplot(2,3,3)
-
-textstr = out2[-1]
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow(reg_output[-1][0][sliceno])
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow(reg_output[-1][0][sliceno],cmap="gray")
-
-###################### LLT_model #########################################
-# * u0 = Im + .03*randn(size(Im)); % adding noise
-# [Den] = LLT_model(single(u0), 10, 0.1, 1);
-#Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0);
-#input, regularization_parameter , time_step, number_of_iterations,
-# tolerance_constant, restrictive_Z_smoothing=0
-out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25,
- time_step=0.0003,
- tolerance_constant=0.0001,
- number_of_iterations=300)
-pars = out2[-2]
-
-reg_output.append(out2)
-
-a=fig.add_subplot(2,3,4)
-
-textstr = out2[-1]
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow(reg_output[-1][0][sliceno],cmap="gray")
-
-
-# ###################### PatchBased_Regul #########################################
-# # Quick 2D denoising example in Matlab:
-# # Im = double(imread('lena_gray_256.tif'))/255; % loading image
-# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
-# # ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05);
-
-out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05,
- searching_window_ratio=3,
- similarity_window_ratio=1,
- PB_filtering_parameter=0.08)
-pars = out2[-2]
-reg_output.append(out2)
-
-a=fig.add_subplot(2,3,5)
-
-
-textstr = out2[-1]
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow(reg_output[-1][0][sliceno],cmap="gray")
-
-
-# ###################### TGV_PD #########################################
-# # Quick 2D denoising example in Matlab:
-# # Im = double(imread('lena_gray_256.tif'))/255; % loading image
-# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
-# # u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550);
-
-
-out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05,
- first_order_term=1.3,
- second_order_term=1,
- number_of_iterations=550)
-pars = out2[-2]
-reg_output.append(out2)
-
-a=fig.add_subplot(2,3,6)
-
-
-textstr = out2[-1]
-
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
- verticalalignment='top', bbox=props)
-imgplot = plt.imshow(reg_output[-1][0][sliceno],cmap="gray")
-
-
-plt.show()
-
-################################################################################
-##
-## 3D Regularizers
-##
-################################################################################
-##Example:
-## figure;
-## Im = double(imread('lena_gray_256.tif'))/255; % loading image
-## u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
-## u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
-#
-##filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Reconstruction\python\test\reconstruction_example.mha"
-#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Simpleflex\data\head.mha"
-#
-#reader = vtk.vtkMetaImageReader()
-#reader.SetFileName(os.path.normpath(filename))
-#reader.Update()
-##vtk returns 3D images, let's take just the one slice there is as 2D
-#Im = Converter.vtk2numpy(reader.GetOutput())
-#Im = Im.astype('float32')
-##imgplot = plt.imshow(Im)
-#perc = 0.05
-#u0 = Im + (perc* np.random.normal(size=np.shape(Im)))
-## map the u0 u0->u0>0
-#f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
-#u0 = f(u0).astype('float32')
-#converter = Converter.numpy2vtkImporter(u0, reader.GetOutput().GetSpacing(),
-# reader.GetOutput().GetOrigin())
-#converter.Update()
-#writer = vtk.vtkMetaImageWriter()
-#writer.SetInputData(converter.GetOutput())
-#writer.SetFileName(r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\noisy_head.mha")
-##writer.Write()
-#
-#
-### plot
-#fig3D = plt.figure()
-##a=fig.add_subplot(3,3,1)
-##a.set_title('Original')
-##imgplot = plt.imshow(Im)
-#sliceNo = 32
-#
-#a=fig3D.add_subplot(2,3,1)
-#a.set_title('noise')
-#imgplot = plt.imshow(u0.T[sliceNo])
-#
-#reg_output3d = []
-#
-###############################################################################
-## Call regularizer
-#
-######################## SplitBregman_TV #####################################
-## u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
-#
-##reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV)
-#
-##out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30,
-## #tolerance_constant=1e-4,
-## TV_Penalty=Regularizer.TotalVariationPenalty.l1)
-#
-#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30,
-# tolerance_constant=1e-4,
-# TV_Penalty=Regularizer.TotalVariationPenalty.l1)
-#
-#
-#pars = out2[-2]
-#reg_output3d.append(out2)
-#
-#a=fig3D.add_subplot(2,3,2)
-#
-#
-#textstr = out2[-1]
-#
-#
-## these are matplotlib.patch.Patch properties
-#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-## place a text box in upper left in axes coords
-#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
-# verticalalignment='top', bbox=props)
-#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
-#
-####################### FGP_TV #########################################
-## u = FGP_TV(single(u0), 0.05, 100, 1e-04);
-#out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.005,
-# number_of_iterations=200)
-#pars = out2[-2]
-#reg_output3d.append(out2)
-#
-#a=fig3D.add_subplot(2,3,2)
-#
-#
-#textstr = out2[-1]
-#
-#
-## these are matplotlib.patch.Patch properties
-#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-## place a text box in upper left in axes coords
-#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
-# verticalalignment='top', bbox=props)
-#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
-#
-####################### LLT_model #########################################
-## * u0 = Im + .03*randn(size(Im)); % adding noise
-## [Den] = LLT_model(single(u0), 10, 0.1, 1);
-##Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0);
-##input, regularization_parameter , time_step, number_of_iterations,
-## tolerance_constant, restrictive_Z_smoothing=0
-#out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25,
-# time_step=0.0003,
-# tolerance_constant=0.0001,
-# number_of_iterations=300)
-#pars = out2[-2]
-#reg_output3d.append(out2)
-#
-#a=fig3D.add_subplot(2,3,2)
-#
-#
-#textstr = out2[-1]
-#
-#
-## these are matplotlib.patch.Patch properties
-#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-## place a text box in upper left in axes coords
-#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
-# verticalalignment='top', bbox=props)
-#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
-#
-####################### PatchBased_Regul #########################################
-## Quick 2D denoising example in Matlab:
-## Im = double(imread('lena_gray_256.tif'))/255; % loading image
-## u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
-## ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05);
-#
-#out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05,
-# searching_window_ratio=3,
-# similarity_window_ratio=1,
-# PB_filtering_parameter=0.08)
-#pars = out2[-2]
-#reg_output3d.append(out2)
-#
-#a=fig3D.add_subplot(2,3,2)
-#
-#
-#textstr = out2[-1]
-#
-#
-## these are matplotlib.patch.Patch properties
-#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-## place a text box in upper left in axes coords
-#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
-# verticalalignment='top', bbox=props)
-#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
-#
-
-###################### TGV_PD #########################################
-# Quick 2D denoising example in Matlab:
-# Im = double(imread('lena_gray_256.tif'))/255; % loading image
-# u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
-# u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550);
-
-
-#out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05,
-# first_order_term=1.3,
-# second_order_term=1,
-# number_of_iterations=550)
-#pars = out2[-2]
-#reg_output3d.append(out2)
-#
-#a=fig3D.add_subplot(2,3,2)
-#
-#
-#textstr = out2[-1]
-#
-#
-## these are matplotlib.patch.Patch properties
-#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-## place a text box in upper left in axes coords
-#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
-# verticalalignment='top', bbox=props)
-#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
--
cgit v1.2.3