summaryrefslogtreecommitdiffstats
path: root/Wrappers/Matlab
diff options
context:
space:
mode:
Diffstat (limited to 'Wrappers/Matlab')
-rw-r--r--Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m9
-rw-r--r--Wrappers/Matlab/demos/demoMatlab_denoise.m9
-rw-r--r--Wrappers/Matlab/demos/demoMatlab_inpaint.m35
-rw-r--r--Wrappers/Matlab/mex_compile/compileCPU_mex.m16
-rw-r--r--Wrappers/Matlab/mex_compile/regularisers_CPU/NonlDiff_Inp.c101
-rw-r--r--Wrappers/Matlab/mex_compile/regularisers_CPU/NonlocalMarching_Inpaint.c82
6 files changed, 241 insertions, 11 deletions
diff --git a/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m b/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m
index 5a54d18..c087433 100644
--- a/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m
+++ b/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m
@@ -1,8 +1,9 @@
% Volume (3D) denoising demo using CCPi-RGL
-clear
-close all
-addpath('../mex_compile/installed');
-addpath('../../../data/');
+clear; close all
+Path1 = sprintf(['..' filesep 'mex_compile' filesep 'installed'], 1i);
+Path2 = sprintf(['..' filesep '..' filesep '..' filesep 'data' filesep], 1i);
+addpath(Path1);
+addpath(Path2);
N = 512;
slices = 30;
diff --git a/Wrappers/Matlab/demos/demoMatlab_denoise.m b/Wrappers/Matlab/demos/demoMatlab_denoise.m
index 151a604..d93f477 100644
--- a/Wrappers/Matlab/demos/demoMatlab_denoise.m
+++ b/Wrappers/Matlab/demos/demoMatlab_denoise.m
@@ -1,8 +1,9 @@
% Image (2D) denoising demo using CCPi-RGL
-clear
-close all
-addpath('../mex_compile/installed');
-addpath('../../../data/');
+clear; close all
+Path1 = sprintf(['..' filesep 'mex_compile' filesep 'installed'], 1i);
+Path2 = sprintf(['..' filesep '..' filesep '..' filesep 'data' filesep], 1i);
+addpath(Path1);
+addpath(Path2);
Im = double(imread('lena_gray_512.tif'))/255; % loading image
u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
diff --git a/Wrappers/Matlab/demos/demoMatlab_inpaint.m b/Wrappers/Matlab/demos/demoMatlab_inpaint.m
new file mode 100644
index 0000000..66f9c15
--- /dev/null
+++ b/Wrappers/Matlab/demos/demoMatlab_inpaint.m
@@ -0,0 +1,35 @@
+% Image (2D) inpainting demo using CCPi-RGL
+clear; close all
+Path1 = sprintf(['..' filesep 'mex_compile' filesep 'installed'], 1i);
+Path2 = sprintf(['..' filesep '..' filesep '..' filesep 'data' filesep], 1i);
+addpath(Path1);
+addpath(Path2);
+
+load('SinoInpaint.mat');
+Sinogram = Sinogram./max(Sinogram(:));
+Sino_mask = Sinogram.*(1-single(Mask));
+figure;
+subplot(1,2,1); imshow(Sino_mask, [0 1]); title('Missing data sinogram');
+subplot(1,2,2); imshow(Mask, [0 1]); title('Mask');
+%%
+fprintf('Inpaint using Linear-Diffusion model (CPU) \n');
+iter_diff = 5000; % number of diffusion iterations
+lambda_regDiff = 6000; % regularisation for the diffusivity
+sigmaPar = 0.0; % edge-preserving parameter
+tau_param = 0.000075; % time-marching constant
+tic; u_diff = NonlDiff_Inp(single(Sino_mask), Mask, lambda_regDiff, sigmaPar, iter_diff, tau_param); toc;
+figure; imshow(u_diff, [0 1]); title('Linear-Diffusion inpainted sinogram (CPU)');
+%%
+fprintf('Inpaint using Nonlinear-Diffusion model (CPU) \n');
+iter_diff = 1500; % number of diffusion iterations
+lambda_regDiff = 80; % regularisation for the diffusivity
+sigmaPar = 0.00009; % edge-preserving parameter
+tau_param = 0.000008; % time-marching constant
+tic; u_diff = NonlDiff_Inp(single(Sino_mask), Mask, lambda_regDiff, sigmaPar, iter_diff, tau_param, 'Huber'); toc;
+figure; imshow(u_diff, [0 1]); title('Non-Linear Diffusion inpainted sinogram (CPU)');
+%%
+fprintf('Inpaint using Nonlocal Vertical Marching model (CPU) \n');
+Increment = 1; % linear increment for the searching window
+tic; [u_nom,maskupd] = NonlocalMarching_Inpaint(single(Sino_mask), Mask, Increment); toc;
+figure; imshow(u_nom, [0 1]); title('NVM inpainted sinogram (CPU)');
+%% \ No newline at end of file
diff --git a/Wrappers/Matlab/mex_compile/compileCPU_mex.m b/Wrappers/Matlab/mex_compile/compileCPU_mex.m
index ee0d99e..b232f33 100644
--- a/Wrappers/Matlab/mex_compile/compileCPU_mex.m
+++ b/Wrappers/Matlab/mex_compile/compileCPU_mex.m
@@ -2,9 +2,11 @@
pathcopyFrom = sprintf(['..' filesep '..' filesep '..' filesep 'Core' filesep 'regularisers_CPU'], 1i);
pathcopyFrom1 = sprintf(['..' filesep '..' filesep '..' filesep 'Core' filesep 'CCPiDefines.h'], 1i);
+pathcopyFrom2 = sprintf(['..' filesep '..' filesep '..' filesep 'Core' filesep 'inpainters_CPU'], 1i);
copyfile(pathcopyFrom, 'regularisers_CPU');
copyfile(pathcopyFrom1, 'regularisers_CPU');
+copyfile(pathcopyFrom2, 'regularisers_CPU');
cd regularisers_CPU
@@ -32,9 +34,17 @@ movefile('NonlDiff.mex*',Pathmove);
mex TV_energy.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
movefile('TV_energy.mex*',Pathmove);
+%############Inpainters##############%
+mex NonlDiff_Inp.c Diffusion_Inpaint_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+movefile('NonlDiff_Inp.mex*',Pathmove);
+
+mex NonlocalMarching_Inpaint.c NonlocalMarching_Inpaint_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+movefile('NonlocalMarching_Inpaint.mex*',Pathmove);
+
delete SB_TV_core* ROF_TV_core* FGP_TV_core* FGP_dTV_core* TNV_core* utils* Diffusion_core* CCPiDefines.h
-fprintf('%s \n', 'All successfully compiled!');
+delete Diffusion_Inpaint_core* NonlocalMarching_Inpaint_core*
+fprintf('%s \n', 'Regularisers successfully compiled!');
-pathA = sprintf(['..' filesep '..' filesep], 1i);
-cd(pathA);
+pathA2 = sprintf(['..' filesep '..' filesep], 1i);
+cd(pathA2);
cd demos
diff --git a/Wrappers/Matlab/mex_compile/regularisers_CPU/NonlDiff_Inp.c b/Wrappers/Matlab/mex_compile/regularisers_CPU/NonlDiff_Inp.c
new file mode 100644
index 0000000..eaab4a7
--- /dev/null
+++ b/Wrappers/Matlab/mex_compile/regularisers_CPU/NonlDiff_Inp.c
@@ -0,0 +1,101 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "matrix.h"
+#include "mex.h"
+#include "Diffusion_Inpaint_core.h"
+
+/* C-OMP implementation of linear and nonlinear diffusion [1,2] for inpainting task (2D/3D case)
+ * The minimisation is performed using explicit scheme.
+ *
+ * Input Parameters:
+ * 1. Image/volume to inpaint
+ * 2. Inpainting Mask of the same size as (1) in 'unsigned char' format (ones mark the region to inpaint, zeros belong to the data)
+ * 3. lambda - regularization parameter
+ * 4. Edge-preserving parameter (sigma), when sigma equals to zero nonlinear diffusion -> linear diffusion
+ * 5. Number of iterations, for explicit scheme >= 150 is recommended
+ * 6. tau - time-marching step for explicit scheme
+ * 7. Penalty type: 1 - Huber, 2 - Perona-Malik, 3 - Tukey Biweight
+ *
+ * Output:
+ * [1] Inpainted image/volume
+ *
+ * This function is based on the paper by
+ * [1] Perona, P. and Malik, J., 1990. Scale-space and edge detection using anisotropic diffusion. IEEE Transactions on pattern analysis and machine intelligence, 12(7), pp.629-639.
+ * [2] Black, M.J., Sapiro, G., Marimont, D.H. and Heeger, D., 1998. Robust anisotropic diffusion. IEEE Transactions on image processing, 7(3), pp.421-432.
+ */
+
+void mexFunction(
+ int nlhs, mxArray *plhs[],
+ int nrhs, const mxArray *prhs[])
+
+{
+ int number_of_dims, iter_numb, dimX, dimY, dimZ, penaltytype, i, inpaint_elements;
+ const int *dim_array;
+ const int *dim_array2;
+ float *Input, *Output=NULL, lambda, tau, sigma;
+ unsigned char *Mask;
+
+ dim_array = mxGetDimensions(prhs[0]);
+ dim_array2 = mxGetDimensions(prhs[1]);
+ number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+
+ /*Handling Matlab input data*/
+ Input = (float *) mxGetData(prhs[0]);
+ Mask = (unsigned char *) mxGetData(prhs[1]); /* MASK */
+ lambda = (float) mxGetScalar(prhs[2]); /* regularization parameter */
+ sigma = (float) mxGetScalar(prhs[3]); /* Edge-preserving parameter */
+ iter_numb = 300; /* iterations number */
+ tau = 0.025; /* marching step parameter */
+ penaltytype = 1; /* Huber penalty by default */
+
+ if ((nrhs < 4) || (nrhs > 7)) mexErrMsgTxt("At least 4 parameters is required, all parameters are: Image(2D/3D), Mask(2D/3D), Regularisation parameter, Edge-preserving parameter, iterations number, time-marching constant, penalty type - Huber, PM or Tukey");
+ if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7)) iter_numb = (int) mxGetScalar(prhs[4]); /* iterations number */
+ if ((nrhs == 6) || (nrhs == 7)) tau = (float) mxGetScalar(prhs[5]); /* marching step parameter */
+ if (nrhs == 7) {
+ char *penalty_type;
+ penalty_type = mxArrayToString(prhs[6]); /* Huber, PM or Tukey 'Huber' is the default */
+ if ((strcmp(penalty_type, "Huber") != 0) && (strcmp(penalty_type, "PM") != 0) && (strcmp(penalty_type, "Tukey") != 0)) mexErrMsgTxt("Choose penalty: 'Huber', 'PM' or 'Tukey',");
+ if (strcmp(penalty_type, "Huber") == 0) penaltytype = 1; /* enable 'Huber' penalty */
+ if (strcmp(penalty_type, "PM") == 0) penaltytype = 2; /* enable Perona-Malik penalty */
+ if (strcmp(penalty_type, "Tukey") == 0) penaltytype = 3; /* enable Tikey Biweight penalty */
+ mxFree(penalty_type);
+ }
+
+ if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+ if (mxGetClassID(prhs[1]) != mxUINT8_CLASS) {mexErrMsgTxt("The mask must be in uint8 precision");}
+
+ dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+
+ /* output arrays*/
+ if (number_of_dims == 2) {
+ dimZ = 1; /*2D case*/
+ /* output image/volume */
+ if ((dimX != dim_array2[0]) || (dimY != dim_array2[1])) mexErrMsgTxt("Input image and the provided mask are of different dimensions!");
+ Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+ }
+ if (number_of_dims == 3) {
+ if ((dimX != dim_array2[0]) || (dimY != dim_array2[1]) || (dimZ != dim_array2[2])) mexErrMsgTxt("Input image and the provided mask are of different dimensions!");
+ Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+ }
+
+ inpaint_elements = 0;
+ for (i=0; i<dimY*dimX*dimZ; i++) if (Mask[i] == 1) inpaint_elements++;
+ if (inpaint_elements == 0) mexErrMsgTxt("The mask is full of zeros, nothing to inpaint");
+ Diffusion_Inpaint_CPU_main(Input, Mask, Output, lambda, sigma, iter_numb, tau, penaltytype, dimX, dimY, dimZ);
+} \ No newline at end of file
diff --git a/Wrappers/Matlab/mex_compile/regularisers_CPU/NonlocalMarching_Inpaint.c b/Wrappers/Matlab/mex_compile/regularisers_CPU/NonlocalMarching_Inpaint.c
new file mode 100644
index 0000000..5e4ab1f
--- /dev/null
+++ b/Wrappers/Matlab/mex_compile/regularisers_CPU/NonlocalMarching_Inpaint.c
@@ -0,0 +1,82 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "matrix.h"
+#include "mex.h"
+#include "NonlocalMarching_Inpaint_core.h"
+
+/* C-OMP implementation of Nonlocal Vertical Marching inpainting method (2D case)
+ * The method is heuristic but computationally efficent (especially for larger images).
+ * It developed specifically to smoothly inpaint horizontal or inclined missing data regions in sinograms
+ * The method WILL not work satisfactory if you have lengthy vertical stripes of missing data
+ *
+ * Input:
+ * 1. 2D image or sinogram [REQUIRED]
+ * 2. Mask of the same size as A in 'unsigned char' format (ones mark the region to inpaint, zeros belong to the data) [REQUIRED]
+ * 3. Linear increment to increase searching window size in iterations, values from 1-3 is a good choice [OPTIONAL, default 1]
+ * 4. Number of iterations [OPTIONAL, default - calculate based on the mask]
+ *
+ * Output:
+ * 1. Inpainted sinogram
+ * 2. updated mask
+ * Reference: TBA
+ */
+
+void mexFunction(
+ int nlhs, mxArray *plhs[],
+ int nrhs, const mxArray *prhs[])
+
+{
+ int number_of_dims, dimX, dimY, dimZ, iterations, SW_increment;
+ const int *dim_array;
+ const int *dim_array2;
+ float *Input, *Output=NULL;
+ unsigned char *Mask, *Mask_upd=NULL;
+
+ dim_array = mxGetDimensions(prhs[0]);
+ dim_array2 = mxGetDimensions(prhs[1]);
+ number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+
+ /*Handling Matlab input data*/
+ Input = (float *) mxGetData(prhs[0]);
+ Mask = (unsigned char *) mxGetData(prhs[1]); /* MASK */
+ SW_increment = 1;
+ iterations = 0;
+
+ if ((nrhs < 2) || (nrhs > 4)) mexErrMsgTxt("At least 4 parameters is required, all parameters are: Image(2D/3D), Mask(2D/3D), Linear increment, Iterations number");
+ if ((nrhs == 3) || (nrhs == 4)) SW_increment = (int) mxGetScalar(prhs[2]); /* linear increment */
+ if ((nrhs == 4)) iterations = (int) mxGetScalar(prhs[3]); /* iterations number */
+
+ if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+ if (mxGetClassID(prhs[1]) != mxUINT8_CLASS) {mexErrMsgTxt("The mask must be in uint8 precision");}
+
+ dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+
+ /* output arrays*/
+ if (number_of_dims == 2) {
+ dimZ = 1; /*2D case*/
+ /* output image/volume */
+ if ((dimX != dim_array2[0]) || (dimY != dim_array2[1])) mexErrMsgTxt("Input image and the provided mask are of different dimensions!");
+ Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+ Mask_upd = (unsigned char*)mxGetPr(plhs[1] = mxCreateNumericArray(2, dim_array, mxUINT8_CLASS, mxREAL));
+ }
+ if (number_of_dims == 3) {
+ mexErrMsgTxt("Currently 2D supported only");
+ }
+ NonlocalMarching_Inpaint_main(Input, Mask, Output, Mask_upd, SW_increment, iterations, dimX, dimY, dimZ);
+} \ No newline at end of file
generated by cgit v1.2.3 (git 2.25.1) at 2024-12-03 20:40:25 +0000