summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Core/CMakeLists.txt4
-rw-r--r--Core/inpainters_CPU/Diffusion_Inpaint_core.c317
-rw-r--r--Core/inpainters_CPU/Diffusion_Inpaint_core.h60
-rw-r--r--Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c158
-rw-r--r--Core/inpainters_CPU/NonlocalMarching_Inpaint_core.h54
-rw-r--r--Core/regularisers_CPU/utils.c33
-rw-r--r--Core/regularisers_CPU/utils.h5
-rw-r--r--Readme.md27
-rw-r--r--Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m9
-rw-r--r--Wrappers/Matlab/demos/demoMatlab_denoise.m9
-rw-r--r--Wrappers/Matlab/demos/demoMatlab_inpaint.m35
-rw-r--r--Wrappers/Matlab/mex_compile/compileCPU_mex.m16
-rw-r--r--Wrappers/Matlab/mex_compile/regularisers_CPU/NonlDiff_Inp.c101
-rw-r--r--Wrappers/Matlab/mex_compile/regularisers_CPU/NonlocalMarching_Inpaint.c82
-rw-r--r--Wrappers/Python/ccpi/filters/regularisers.py7
-rw-r--r--Wrappers/Python/setup-regularisers.py.in3
-rw-r--r--Wrappers/Python/src/cpu_regularisers.pyx50
-rw-r--r--data/SinoInpaint.matbin0 -> 3303762 bytes
18 files changed, 937 insertions, 33 deletions
diff --git a/Core/CMakeLists.txt b/Core/CMakeLists.txt
index 61986dc..322732b 100644
--- a/Core/CMakeLists.txt
+++ b/Core/CMakeLists.txt
@@ -85,7 +85,7 @@ message("Adding regularisers as a shared library")
add_library(cilreg SHARED
${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/FGP_TV_core.c
${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/SB_TV_core.c
- ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/Diffusion_core.c
+ ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/Diffusion_core.c
#${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/LLT_model_core.c
#${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/PatchBased_Regul_core.c
#${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/TGV_PD_core.c
@@ -93,6 +93,8 @@ add_library(cilreg SHARED
${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/FGP_dTV_core.c
${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/TNV_core.c
${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/utils.c
+ ${CMAKE_CURRENT_SOURCE_DIR}/inpainters_CPU/Diffusion_Inpaint_core.c
+ ${CMAKE_CURRENT_SOURCE_DIR}/inpainters_CPU/NonlocalMarching_Inpaint_core.c
)
target_link_libraries(cilreg ${EXTRA_LIBRARIES} )
include_directories(cilreg PUBLIC
diff --git a/Core/inpainters_CPU/Diffusion_Inpaint_core.c b/Core/inpainters_CPU/Diffusion_Inpaint_core.c
new file mode 100644
index 0000000..16e87de
--- /dev/null
+++ b/Core/inpainters_CPU/Diffusion_Inpaint_core.c
@@ -0,0 +1,317 @@
+/*
+ * 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 "Diffusion_Inpaint_core.h"
+#include "utils.h"
+
+/*sign function*/
+int signNDFc(float x) {
+ return (x > 0) - (x < 0);
+}
+
+/* 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. 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.
+ */
+
+float Diffusion_Inpaint_CPU_main(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, int dimX, int dimY, int dimZ)
+{
+ int i;
+ float sigmaPar2;
+ sigmaPar2 = sigmaPar/sqrt(2.0f);
+
+ /* copy into output */
+ copyIm(Input, Output, dimX, dimY, dimZ);
+
+ if (dimZ == 1) {
+ /* running 2D diffusion iterations */
+ for(i=0; i < iterationsNumb; i++) {
+ if (sigmaPar == 0.0f) LinearDiff_Inp_2D(Input, Mask, Output, lambdaPar, tau, dimX, dimY); /* linear diffusion (heat equation) */
+ else NonLinearDiff_Inp_2D(Input, Mask, Output, lambdaPar, sigmaPar2, tau, penaltytype, dimX, dimY); /* nonlinear diffusion */
+ }
+ }
+ else {
+ /* running 3D diffusion iterations */
+ for(i=0; i < iterationsNumb; i++) {
+ if (sigmaPar == 0.0f) LinearDiff_Inp_3D(Input, Mask, Output, lambdaPar, tau, dimX, dimY, dimZ);
+ else NonLinearDiff_Inp_3D(Input, Mask, Output, lambdaPar, sigmaPar2, tau, penaltytype, dimX, dimY, dimZ);
+ }
+ }
+ return *Output;
+}
+
+
+/********************************************************************/
+/***************************2D Functions*****************************/
+/********************************************************************/
+/* linear diffusion (heat equation) */
+float LinearDiff_Inp_2D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float tau, int dimX, int dimY)
+{
+ int i,j,i1,i2,j1,j2,index;
+ float e,w,n,s,e1,w1,n1,s1;
+
+#pragma omp parallel for shared(Input,Mask) private(index,i,j,i1,i2,j1,j2,e,w,n,s,e1,w1,n1,s1)
+ for(i=0; i<dimX; i++) {
+ /* symmetric boundary conditions (Neuman) */
+ i1 = i+1; if (i1 == dimX) i1 = i-1;
+ i2 = i-1; if (i2 < 0) i2 = i+1;
+ for(j=0; j<dimY; j++) {
+ /* symmetric boundary conditions (Neuman) */
+ j1 = j+1; if (j1 == dimY) j1 = j-1;
+ j2 = j-1; if (j2 < 0) j2 = j+1;
+ index = j*dimX+i;
+
+ if (Mask[index] > 0) {
+ /*inpainting process*/
+ e = Output[j*dimX+i1];
+ w = Output[j*dimX+i2];
+ n = Output[j1*dimX+i];
+ s = Output[j2*dimX+i];
+
+ e1 = e - Output[index];
+ w1 = w - Output[index];
+ n1 = n - Output[index];
+ s1 = s - Output[index];
+
+ Output[index] += tau*(lambdaPar*(e1 + w1 + n1 + s1) - (Output[index] - Input[index]));
+ }
+ }}
+ return *Output;
+}
+
+/* nonlinear diffusion */
+float NonLinearDiff_Inp_2D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, int dimX, int dimY)
+{
+ int i,j,i1,i2,j1,j2,index;
+ float e,w,n,s,e1,w1,n1,s1;
+
+#pragma omp parallel for shared(Input,Mask) private(index,i,j,i1,i2,j1,j2,e,w,n,s,e1,w1,n1,s1)
+ for(i=0; i<dimX; i++) {
+ /* symmetric boundary conditions (Neuman) */
+ i1 = i+1; if (i1 == dimX) i1 = i-1;
+ i2 = i-1; if (i2 < 0) i2 = i+1;
+ for(j=0; j<dimY; j++) {
+ /* symmetric boundary conditions (Neuman) */
+ j1 = j+1; if (j1 == dimY) j1 = j-1;
+ j2 = j-1; if (j2 < 0) j2 = j+1;
+ index = j*dimX+i;
+
+ if (Mask[index] > 0) {
+ /*inpainting process*/
+ e = Output[j*dimX+i1];
+ w = Output[j*dimX+i2];
+ n = Output[j1*dimX+i];
+ s = Output[j2*dimX+i];
+
+ e1 = e - Output[index];
+ w1 = w - Output[index];
+ n1 = n - Output[index];
+ s1 = s - Output[index];
+
+ if (penaltytype == 1){
+ /* Huber penalty */
+ if (fabs(e1) > sigmaPar) e1 = signNDFc(e1);
+ else e1 = e1/sigmaPar;
+
+ if (fabs(w1) > sigmaPar) w1 = signNDFc(w1);
+ else w1 = w1/sigmaPar;
+
+ if (fabs(n1) > sigmaPar) n1 = signNDFc(n1);
+ else n1 = n1/sigmaPar;
+
+ if (fabs(s1) > sigmaPar) s1 = signNDFc(s1);
+ else s1 = s1/sigmaPar;
+ }
+ else if (penaltytype == 2) {
+ /* Perona-Malik */
+ e1 = (e1)/(1.0f + powf((e1/sigmaPar),2));
+ w1 = (w1)/(1.0f + powf((w1/sigmaPar),2));
+ n1 = (n1)/(1.0f + powf((n1/sigmaPar),2));
+ s1 = (s1)/(1.0f + powf((s1/sigmaPar),2));
+ }
+ else if (penaltytype == 3) {
+ /* Tukey Biweight */
+ if (fabs(e1) <= sigmaPar) e1 = e1*powf((1.0f - powf((e1/sigmaPar),2)), 2);
+ else e1 = 0.0f;
+ if (fabs(w1) <= sigmaPar) w1 = w1*powf((1.0f - powf((w1/sigmaPar),2)), 2);
+ else w1 = 0.0f;
+ if (fabs(n1) <= sigmaPar) n1 = n1*powf((1.0f - powf((n1/sigmaPar),2)), 2);
+ else n1 = 0.0f;
+ if (fabs(s1) <= sigmaPar) s1 = s1*powf((1.0f - powf((s1/sigmaPar),2)), 2);
+ else s1 = 0.0f;
+ }
+ else {
+ printf("%s \n", "No penalty function selected! Use 1,2 or 3.");
+ break;
+ }
+ Output[index] += tau*(lambdaPar*(e1 + w1 + n1 + s1) - (Output[index] - Input[index]));
+ }
+ }}
+ return *Output;
+}
+/********************************************************************/
+/***************************3D Functions*****************************/
+/********************************************************************/
+/* linear diffusion (heat equation) */
+float LinearDiff_Inp_3D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float tau, int dimX, int dimY, int dimZ)
+{
+ int i,j,k,i1,i2,j1,j2,k1,k2,index;
+ float e,w,n,s,u,d,e1,w1,n1,s1,u1,d1;
+
+#pragma omp parallel for shared(Input,Mask) private(index,i,j,i1,i2,j1,j2,e,w,n,s,e1,w1,n1,s1,k,k1,k2,u1,d1,u,d)
+for(k=0; k<dimZ; k++) {
+ k1 = k+1; if (k1 == dimZ) k1 = k-1;
+ k2 = k-1; if (k2 < 0) k2 = k+1;
+ for(i=0; i<dimX; i++) {
+ /* symmetric boundary conditions (Neuman) */
+ i1 = i+1; if (i1 == dimX) i1 = i-1;
+ i2 = i-1; if (i2 < 0) i2 = i+1;
+ for(j=0; j<dimY; j++) {
+ /* symmetric boundary conditions (Neuman) */
+ j1 = j+1; if (j1 == dimY) j1 = j-1;
+ j2 = j-1; if (j2 < 0) j2 = j+1;
+ index = (dimX*dimY)*k + j*dimX+i;
+
+ if (Mask[index] > 0) {
+ /*inpainting process*/
+
+ e = Output[(dimX*dimY)*k + j*dimX+i1];
+ w = Output[(dimX*dimY)*k + j*dimX+i2];
+ n = Output[(dimX*dimY)*k + j1*dimX+i];
+ s = Output[(dimX*dimY)*k + j2*dimX+i];
+ u = Output[(dimX*dimY)*k1 + j*dimX+i];
+ d = Output[(dimX*dimY)*k2 + j*dimX+i];
+
+ e1 = e - Output[index];
+ w1 = w - Output[index];
+ n1 = n - Output[index];
+ s1 = s - Output[index];
+ u1 = u - Output[index];
+ d1 = d - Output[index];
+
+ Output[index] += tau*(lambdaPar*(e1 + w1 + n1 + s1 + u1 + d1) - (Output[index] - Input[index]));
+ }
+ }}}
+ return *Output;
+}
+
+float NonLinearDiff_Inp_3D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, int dimX, int dimY, int dimZ)
+{
+ int i,j,k,i1,i2,j1,j2,k1,k2,index;
+ float e,w,n,s,u,d,e1,w1,n1,s1,u1,d1;
+
+#pragma omp parallel for shared(Input,Mask) private(index,i,j,i1,i2,j1,j2,e,w,n,s,e1,w1,n1,s1,k,k1,k2,u1,d1,u,d)
+for(k=0; k<dimZ; k++) {
+ k1 = k+1; if (k1 == dimZ) k1 = k-1;
+ k2 = k-1; if (k2 < 0) k2 = k+1;
+ for(i=0; i<dimX; i++) {
+ /* symmetric boundary conditions (Neuman) */
+ i1 = i+1; if (i1 == dimX) i1 = i-1;
+ i2 = i-1; if (i2 < 0) i2 = i+1;
+ for(j=0; j<dimY; j++) {
+ /* symmetric boundary conditions (Neuman) */
+ j1 = j+1; if (j1 == dimY) j1 = j-1;
+ j2 = j-1; if (j2 < 0) j2 = j+1;
+ index = (dimX*dimY)*k + j*dimX+i;
+
+ if (Mask[index] > 0) {
+ /*inpainting process*/
+ e = Output[(dimX*dimY)*k + j*dimX+i1];
+ w = Output[(dimX*dimY)*k + j*dimX+i2];
+ n = Output[(dimX*dimY)*k + j1*dimX+i];
+ s = Output[(dimX*dimY)*k + j2*dimX+i];
+ u = Output[(dimX*dimY)*k1 + j*dimX+i];
+ d = Output[(dimX*dimY)*k2 + j*dimX+i];
+
+ e1 = e - Output[index];
+ w1 = w - Output[index];
+ n1 = n - Output[index];
+ s1 = s - Output[index];
+ u1 = u - Output[index];
+ d1 = d - Output[index];
+
+ if (penaltytype == 1){
+ /* Huber penalty */
+ if (fabs(e1) > sigmaPar) e1 = signNDFc(e1);
+ else e1 = e1/sigmaPar;
+
+ if (fabs(w1) > sigmaPar) w1 = signNDFc(w1);
+ else w1 = w1/sigmaPar;
+
+ if (fabs(n1) > sigmaPar) n1 = signNDFc(n1);
+ else n1 = n1/sigmaPar;
+
+ if (fabs(s1) > sigmaPar) s1 = signNDFc(s1);
+ else s1 = s1/sigmaPar;
+
+ if (fabs(u1) > sigmaPar) u1 = signNDFc(u1);
+ else u1 = u1/sigmaPar;
+
+ if (fabs(d1) > sigmaPar) d1 = signNDFc(d1);
+ else d1 = d1/sigmaPar;
+ }
+ else if (penaltytype == 2) {
+ /* Perona-Malik */
+ e1 = (e1)/(1.0f + powf((e1/sigmaPar),2));
+ w1 = (w1)/(1.0f + powf((w1/sigmaPar),2));
+ n1 = (n1)/(1.0f + powf((n1/sigmaPar),2));
+ s1 = (s1)/(1.0f + powf((s1/sigmaPar),2));
+ u1 = (u1)/(1.0f + powf((u1/sigmaPar),2));
+ d1 = (d1)/(1.0f + powf((d1/sigmaPar),2));
+ }
+ else if (penaltytype == 3) {
+ /* Tukey Biweight */
+ if (fabs(e1) <= sigmaPar) e1 = e1*powf((1.0f - powf((e1/sigmaPar),2)), 2);
+ else e1 = 0.0f;
+ if (fabs(w1) <= sigmaPar) w1 = w1*powf((1.0f - powf((w1/sigmaPar),2)), 2);
+ else w1 = 0.0f;
+ if (fabs(n1) <= sigmaPar) n1 = n1*powf((1.0f - powf((n1/sigmaPar),2)), 2);
+ else n1 = 0.0f;
+ if (fabs(s1) <= sigmaPar) s1 = s1*powf((1.0f - powf((s1/sigmaPar),2)), 2);
+ else s1 = 0.0f;
+ if (fabs(u1) <= sigmaPar) u1 = u1*powf((1.0f - powf((u1/sigmaPar),2)), 2);
+ else u1 = 0.0f;
+ if (fabs(d1) <= sigmaPar) d1 = d1*powf((1.0f - powf((d1/sigmaPar),2)), 2);
+ else d1 = 0.0f;
+ }
+ else {
+ printf("%s \n", "No penalty function selected! Use 1,2 or 3.");
+ break;
+ }
+
+ Output[index] += tau*(lambdaPar*(e1 + w1 + n1 + s1 + u1 + d1) - (Output[index] - Input[index]));
+ }
+ }}}
+ return *Output;
+}
diff --git a/Core/inpainters_CPU/Diffusion_Inpaint_core.h b/Core/inpainters_CPU/Diffusion_Inpaint_core.h
new file mode 100644
index 0000000..0fc2608
--- /dev/null
+++ b/Core/inpainters_CPU/Diffusion_Inpaint_core.h
@@ -0,0 +1,60 @@
+/*
+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 <math.h>
+#include <stdlib.h>
+#include <memory.h>
+#include <stdio.h>
+#include "omp.h"
+#include "utils.h"
+#include "CCPiDefines.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. 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.
+ */
+
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+CCPI_EXPORT float Diffusion_Inpaint_CPU_main(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, int dimX, int dimY, int dimZ);
+CCPI_EXPORT float LinearDiff_Inp_2D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float tau, int dimX, int dimY);
+CCPI_EXPORT float NonLinearDiff_Inp_2D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, int dimX, int dimY);
+CCPI_EXPORT float LinearDiff_Inp_3D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float tau, int dimX, int dimY, int dimZ);
+CCPI_EXPORT float NonLinearDiff_Inp_3D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, int dimX, int dimY, int dimZ);
+#ifdef __cplusplus
+}
+#endif
diff --git a/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c b/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c
new file mode 100644
index 0000000..d6cdf93
--- /dev/null
+++ b/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c
@@ -0,0 +1,158 @@
+/*
+ * 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 "NonlocalMarching_Inpaint_core.h"
+#include "utils.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 with horizontal or inclined regions of missing data
+ * 2. Mask of the same size as A in 'unsigned char' format (ones mark the region to inpaint, zeros belong to the data)
+ * 3. Linear increment to increase searching window size in iterations, values from 1-3 is a good choice
+
+ * Output:
+ * 1. Inpainted image or a sinogram
+ * 2. updated mask
+ *
+ * Reference: TBA
+ */
+
+float NonlocalMarching_Inpaint_main(float *Input, unsigned char *M, float *Output, unsigned char *M_upd, int SW_increment, int iterationsNumb, int dimX, int dimY, int dimZ)
+{
+ int i, j, i_m, j_m, counter, iter, iterations_number, W_fullsize, switchmask, switchcurr, counterElements;
+ float *Gauss_weights;
+
+ /* copying M to M_upd */
+ copyIm_unchar(M, M_upd, dimX, dimY, 1);
+
+ /* Copying the image */
+ copyIm(Input, Output, dimX, dimY, 1);
+
+ /* Find how many inpainting iterations (equal to the number of ones) required based on a mask */
+ if (iterationsNumb == 0) {
+ iterations_number = 0;
+ for (i=0; i<dimY*dimX; i++) {
+ if (M[i] == 1) iterations_number++;
+ }
+ }
+ else iterations_number = iterationsNumb;
+ if (iterations_number == 0) printf("%s \n", "Nothing to inpaint, zero mask!");
+ else {
+
+ printf("%s %i \n", "Max iteration number equals to:", iterations_number);
+
+ /* Inpainting iterations run here*/
+ int W_halfsize = 1;
+ for(iter=0; iter < iterations_number; iter++) {
+
+ //if (mod (iter, 2) == 0) {W_halfsize += 1;}
+ // printf("%i \n", W_halfsize);
+
+ /* pre-calculation of Gaussian distance weights */
+ W_fullsize = (int)(2*W_halfsize + 1); /*full size of similarity window */
+ Gauss_weights = (float*)calloc(W_fullsize*W_fullsize,sizeof(float ));
+ counter = 0;
+ for(i_m=-W_halfsize; i_m<=W_halfsize; i_m++) {
+ for(j_m=-W_halfsize; j_m<=W_halfsize; j_m++) {
+ Gauss_weights[counter] = exp(-(pow((i_m), 2) + pow((j_m), 2))/(2*W_halfsize*W_halfsize));
+ counter++;
+ }
+ }
+ /* find a point in the mask to inpaint */
+#pragma omp parallel for shared(Output, M_upd, Gauss_weights) private(i, j, switchmask, switchcurr)
+ for(j=0; j<dimY; j++) {
+ switchmask = 0;
+ for(i=0; i<dimX; i++) {
+ switchcurr = 0;
+ if ((M_upd[j*dimX + i] == 1) && (switchmask == 0)) {
+ /* perform inpainting of the current pixel */
+ inpaint_func(Output, M_upd, Gauss_weights, i, j, dimX, dimY, W_halfsize, W_fullsize);
+ /* add value to the mask*/
+ M_upd[j*dimX + i] = 0;
+ switchmask = 1; switchcurr = 1;
+ }
+ if ((M_upd[j*dimX + i] == 0) && (switchmask == 1) && (switchcurr == 0)) {
+ /* perform inpainting of the previous (j-1) pixel */
+ inpaint_func(Output, M_upd, Gauss_weights, i-1, j, dimX, dimY, W_halfsize, W_fullsize);
+ /* add value to the mask*/
+ M_upd[(j)*dimX + i-1] = 0;
+ switchmask = 0;
+ }
+ }
+ }
+ free(Gauss_weights);
+
+ /* check if possible to terminate iterations earlier */
+ counterElements = 0;
+ for(i=0; i<dimX*dimY; i++) if (M_upd[i] == 0) counterElements++;
+
+ if (counterElements == dimX*dimY) {
+ printf("%s \n", "Padding completed!");
+ break;
+ }
+ W_halfsize += SW_increment;
+ }
+ printf("%s %i \n", "Iterations stopped at:", iter);
+ }
+ return *Output;
+}
+
+float inpaint_func(float *U, unsigned char *M_upd, float *Gauss_weights, int i, int j, int dimX, int dimY, int W_halfsize, int W_fullsize)
+{
+ int i1, j1, i_m, j_m, counter;
+ float sum_val, sumweight;
+
+ /*method 1: inpainting based on Euclidian weights */
+ sumweight = 0.0f;
+ counter = 0; sum_val = 0.0f;
+ for(i_m=-W_halfsize; i_m<=W_halfsize; i_m++) {
+ i1 = i+i_m;
+ for(j_m=-W_halfsize; j_m<=W_halfsize; j_m++) {
+ j1 = j+j_m;
+ if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) {
+ if (M_upd[j1*dimX + i1] == 0) {
+ sumweight += Gauss_weights[counter];
+ }
+ }
+ counter++;
+ }
+ }
+ counter = 0; sum_val = 0.0f;
+ for(i_m=-W_halfsize; i_m<=W_halfsize; i_m++) {
+ i1 = i+i_m;
+ for(j_m=-W_halfsize; j_m<=W_halfsize; j_m++) {
+ j1 = j+j_m;
+ if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) {
+ if ((M_upd[j1*dimX + i1] == 0) && (sumweight != 0.0f)) {
+ /* we have data so add it with Euc weight */
+ sum_val += (Gauss_weights[counter]/sumweight)*U[j1*dimX + i1];
+ }
+ }
+ counter++;
+ }
+ }
+ U[j*dimX + i] = sum_val;
+ return *U;
+}
+
diff --git a/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.h b/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.h
new file mode 100644
index 0000000..199820c
--- /dev/null
+++ b/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.h
@@ -0,0 +1,54 @@
+/*
+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 <math.h>
+#include <stdlib.h>
+#include <memory.h>
+#include <stdio.h>
+#include "omp.h"
+#include "utils.h"
+#include "CCPiDefines.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
+ *
+ * Inputs:
+ * 1. 2D image or sinogram with horizontal or inclined regions of missing data
+ * 2. Mask of the same size as A in 'unsigned char' format (ones mark the region to inpaint, zeros belong to the data)
+ * 3. Linear increment to increase searching window size in iterations, values from 1-3 is a good choice
+
+ * Output:
+ * 1. Inpainted image or a sinogram
+ * 2. updated mask
+ *
+ * Reference: TBA
+ */
+
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+CCPI_EXPORT float NonlocalMarching_Inpaint_main(float *Input, unsigned char *M, float *Output, unsigned char *M_upd, int SW_increment, int iterationsNumb, int dimX, int dimY, int dimZ);
+CCPI_EXPORT float inpaint_func(float *U, unsigned char *M_upd, float *Gauss_weights, int i, int j, int dimX, int dimY, int W_halfsize, int W_fullsize);
+#ifdef __cplusplus
+}
+#endif
diff --git a/Core/regularisers_CPU/utils.c b/Core/regularisers_CPU/utils.c
index ca5c59a..3eebdaa 100644
--- a/Core/regularisers_CPU/utils.c
+++ b/Core/regularisers_CPU/utils.c
@@ -20,7 +20,7 @@ limitations under the License.
#include "utils.h"
#include <math.h>
-/* Copy Image */
+/* Copy Image (float) */
float copyIm(float *A, float *U, int dimX, int dimY, int dimZ)
{
int j;
@@ -29,12 +29,33 @@ float copyIm(float *A, float *U, int dimX, int dimY, int dimZ)
return *U;
}
-/*static inline int8_t SIGN(int val) {
- if (val < 0) return -1;
- if (val==0) return 0;
- return 1;
+/* Copy Image */
+unsigned char copyIm_unchar(unsigned char *A, unsigned char *U, int dimX, int dimY, int dimZ)
+{
+ int j;
+#pragma omp parallel for shared(A, U) private(j)
+ for (j = 0; j<dimX*dimY*dimZ; j++) U[j] = A[j];
+ return *U;
+}
+
+/*Roll image symmetrically from top to bottom*/
+float copyIm_roll(float *A, float *U, int dimX, int dimY, int roll_value, int switcher)
+{
+ int i, j;
+#pragma omp parallel for shared(U, A) private(i,j)
+ for (i=0; i<dimX; i++) {
+ for (j=0; j<dimY; j++) {
+ if (switcher == 0) {
+ if (j < (dimY - roll_value)) U[j*dimX + i] = A[(j+roll_value)*dimX + i];
+ else U[j*dimX + i] = A[(j - (dimY - roll_value))*dimX + i];
+ }
+ else {
+ if (j < roll_value) U[j*dimX + i] = A[(j+(dimY - roll_value))*dimX + i];
+ else U[j*dimX + i] = A[(j - roll_value)*dimX + i];
+ }
+ }}
+ return *U;
}
-*/
/* function that calculates TV energy
* type - 1: 2*lambda*min||\nabla u|| + ||u -u0||^2
diff --git a/Core/regularisers_CPU/utils.h b/Core/regularisers_CPU/utils.h
index 866bc01..7ad83ab 100644
--- a/Core/regularisers_CPU/utils.h
+++ b/Core/regularisers_CPU/utils.h
@@ -17,17 +17,16 @@ See the License for the specific language governing permissions and
limitations under the License.
*/
-//#include <matrix.h>
-//#include <math.h>
#include <stdlib.h>
#include <memory.h>
#include "CCPiDefines.h"
-//#include <stdio.h>
#include "omp.h"
#ifdef __cplusplus
extern "C" {
#endif
CCPI_EXPORT float copyIm(float *A, float *U, int dimX, int dimY, int dimZ);
+CCPI_EXPORT unsigned char copyIm_unchar(unsigned char *A, unsigned char *U, int dimX, int dimY, int dimZ);
+CCPI_EXPORT float copyIm_roll(float *A, float *U, int dimX, int dimY, int roll_value, int switcher);
CCPI_EXPORT float TV_energy2D(float *U, float *U0, float *E_val, float lambda, int type, int dimX, int dimY);
CCPI_EXPORT float TV_energy3D(float *U, float *U0, float *E_val, float lambda, int type, int dimX, int dimY, int dimZ);
#ifdef __cplusplus
diff --git a/Readme.md b/Readme.md
index e73b4fb..8a67c60 100644
--- a/Readme.md
+++ b/Readme.md
@@ -16,17 +16,22 @@ the toolkit can be used independently to solve image denoising problems. The cor
* C compilers
* nvcc (CUDA SDK) compilers
-## Package modules (regularisers):
+## Package modules:
-### Single-channel
+### Single-channel (denoising):
1. Rudin-Osher-Fatemi (ROF) Total Variation (explicit PDE minimisation scheme) **2D/3D CPU/GPU** (Ref. *1*)
2. Fast-Gradient-Projection (FGP) Total Variation **2D/3D CPU/GPU** (Ref. *2*)
-3. Split-Bregman (SB) Total Variation **2D/3D CPU/GPU** (Ref. *4*)
-4. Linear and nonlinear diffusion (explicit PDE minimisation scheme) **2D/3D CPU/GPU** (Ref. *6*)
+3. Split-Bregman (SB) Total Variation **2D/3D CPU/GPU** (Ref. *5*)
+4. Linear and nonlinear diffusion (explicit PDE minimisation scheme) **2D/3D CPU/GPU** (Ref. *7*)
+
+### Multi-channel (denoising):
+1. Fast-Gradient-Projection (FGP) Directional Total Variation **2D/3D CPU/GPU** (Ref. *3,4,2*)
+2. Total Nuclear Variation (TNV) penalty **2D+channels CPU** (Ref. *6*)
+
+### Inpainting:
+1. Linear and nonlinear diffusion (explicit PDE minimisation scheme) **2D/3D CPU** (Ref. *7*)
+2. Iterative nonlocal vertical marching method **2D CPU**
-### Multi-channel
-1. Fast-Gradient-Projection (FGP) Directional Total Variation **2D/3D CPU/GPU** (Ref. *3,2*)
-2. Total Nuclear Variation (TNV) penalty **2D+channels CPU** (Ref. *5*)
## Installation:
@@ -56,11 +61,13 @@ the toolkit can be used independently to solve image denoising problems. The cor
*3. Ehrhardt, M.J. and Betcke, M.M., 2016. Multicontrast MRI reconstruction with structure-guided total variation. SIAM Journal on Imaging Sciences, 9(3), pp.1084-1106.*
-*4. Goldstein, T. and Osher, S., 2009. The split Bregman method for L1-regularized problems. SIAM journal on imaging sciences, 2(2), pp.323-343.*
+*4. Kazantsev, D., Jørgensen, J.S., Andersen, M., Lionheart, W.R., Lee, P.D. and Withers, P.J., 2018. Joint image reconstruction method with correlative multi-channel prior for X-ray spectral computed tomography. Inverse Problems*
+
+*5. Goldstein, T. and Osher, S., 2009. The split Bregman method for L1-regularized problems. SIAM journal on imaging sciences, 2(2), pp.323-343.*
-*5. Duran, J., Moeller, M., Sbert, C. and Cremers, D., 2016. Collaborative total variation: a general framework for vectorial TV models. SIAM Journal on Imaging Sciences, 9(1), pp.116-151.*
+*6. Duran, J., Moeller, M., Sbert, C. and Cremers, D., 2016. Collaborative total variation: a general framework for vectorial TV models. SIAM Journal on Imaging Sciences, 9(1), pp.116-151.*
-*6. 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.*
+*7. 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.*
### License:
[Apache License, Version 2.0](http://www.apache.org/licenses/LICENSE-2.0)
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
diff --git a/Wrappers/Python/ccpi/filters/regularisers.py b/Wrappers/Python/ccpi/filters/regularisers.py
index eec8c4d..e62c020 100644
--- a/Wrappers/Python/ccpi/filters/regularisers.py
+++ b/Wrappers/Python/ccpi/filters/regularisers.py
@@ -2,7 +2,7 @@
script which assigns a proper device core function based on a flag ('cpu' or 'gpu')
"""
-from ccpi.filters.cpu_regularisers_cython import TV_ROF_CPU, TV_FGP_CPU, TV_SB_CPU, dTV_FGP_CPU, TNV_CPU, NDF_CPU
+from ccpi.filters.cpu_regularisers import TV_ROF_CPU, TV_FGP_CPU, TV_SB_CPU, dTV_FGP_CPU, TNV_CPU, NDF_CPU, NDF_INPAINT_CPU
from ccpi.filters.gpu_regularisers import TV_ROF_GPU, TV_FGP_GPU, TV_SB_GPU, dTV_FGP_GPU, NDF_GPU
def ROF_TV(inputData, regularisation_parameter, iterations,
@@ -110,3 +110,8 @@ def NDF(inputData, regularisation_parameter, edge_parameter, iterations,
else:
raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
.format(device))
+def NDF_INP(inputData, maskData, regularisation_parameter, edge_parameter, iterations,
+ time_marching_parameter, penalty_type):
+ return NDF_INPAINT_CPU(inputData, maskData, regularisation_parameter,
+ edge_parameter, iterationsNumb,
+ time_marching_parameter, penalty_type)
diff --git a/Wrappers/Python/setup-regularisers.py.in b/Wrappers/Python/setup-regularisers.py.in
index b900efe..f55c6fe 100644
--- a/Wrappers/Python/setup-regularisers.py.in
+++ b/Wrappers/Python/setup-regularisers.py.in
@@ -34,6 +34,7 @@ extra_libraries = ['cilreg']
extra_include_dirs += [os.path.join(".." , ".." , "Core"),
os.path.join(".." , ".." , "Core", "regularisers_CPU"),
+ os.path.join(".." , ".." , "Core", "inpainters_CPU"),
os.path.join(".." , ".." , "Core", "regularisers_GPU" , "TV_FGP" ) ,
os.path.join(".." , ".." , "Core", "regularisers_GPU" , "TV_ROF" ) ,
os.path.join(".." , ".." , "Core", "regularisers_GPU" , "TV_SB" ) ,
@@ -52,7 +53,7 @@ setup(
description='CCPi Core Imaging Library - Image regularisers',
version=cil_version,
cmdclass = {'build_ext': build_ext},
- ext_modules = [Extension("ccpi.filters.cpu_regularisers_cython",
+ ext_modules = [Extension("ccpi.filters.cpu_regularisers",
sources=[os.path.join("." , "src", "cpu_regularisers.pyx" ) ],
include_dirs=extra_include_dirs,
library_dirs=extra_library_dirs,
diff --git a/Wrappers/Python/src/cpu_regularisers.pyx b/Wrappers/Python/src/cpu_regularisers.pyx
index 7ed8fa1..3625106 100644
--- a/Wrappers/Python/src/cpu_regularisers.pyx
+++ b/Wrappers/Python/src/cpu_regularisers.pyx
@@ -25,6 +25,9 @@ cdef extern float Diffusion_CPU_main(float *Input, float *Output, float lambdaPa
cdef extern float TNV_CPU_main(float *Input, float *u, float lambdaPar, int maxIter, float tol, int dimX, int dimY, int dimZ);
cdef extern float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float lambdaPar, int iterationsNumb, float epsil, float eta, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ);
+cdef extern float Diffusion_Inpaint_CPU_main(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, int dimX, int dimY, int dimZ);
+#cdef extern float NonlocalMarching_Inpaint_main(float *Input, unsigned char *M, float *Output, unsigned char *M_upd, int SW_increment, int iterationsNumb, int dimX, int dimY, int dimZ);
+
#****************************************************************#
#********************** Total-variation ROF *********************#
#****************************************************************#
@@ -319,3 +322,50 @@ def NDF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
Diffusion_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[2], dims[1], dims[0])
return outputData
+
+#*********************Inpainting WITH****************************#
+#***************Nonlinear (Isotropic) Diffusion******************#
+#****************************************************************#
+def NDF_INPAINT_CPU(inputData, maskData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type):
+ if inputData.ndim == 2:
+ return NDF_INP_2D(inputData, maskData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type)
+ elif inputData.ndim == 3:
+ return NDF_INP_3D(inputData, maskData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type)
+
+def NDF_INP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
+ np.ndarray[np.uint8_t, ndim=2, mode="c"] maskData,
+ float regularisation_parameter,
+ float edge_parameter,
+ int iterationsNumb,
+ float time_marching_parameter,
+ int penalty_type):
+ cdef long dims[2]
+ dims[0] = inputData.shape[0]
+ dims[1] = inputData.shape[1]
+
+ cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
+ np.zeros([dims[0],dims[1]], dtype='float32')
+
+ # Run Inpaiting by Diffusion iterations for 2D data
+ Diffusion_Inpaint_CPU_main(&inputData[0,0], &maskData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[0], dims[1], 1)
+ return outputData
+
+def NDF_INP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
+ np.ndarray[np.uint8_t, ndim=3, mode="c"] maskData,
+ float regularisation_parameter,
+ float edge_parameter,
+ int iterationsNumb,
+ float time_marching_parameter,
+ int penalty_type):
+ cdef long dims[3]
+ dims[0] = inputData.shape[0]
+ dims[1] = inputData.shape[1]
+ dims[2] = inputData.shape[2]
+
+ cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
+ np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
+
+ # Run Inpaiting by Diffusion iterations for 3D data
+ Diffusion_Inpaint_CPU_main(&inputData[0,0,0], &maskData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[2], dims[1], dims[0])
+
+ return outputData
diff --git a/data/SinoInpaint.mat b/data/SinoInpaint.mat
new file mode 100644
index 0000000..d2fbe12
--- /dev/null
+++ b/data/SinoInpaint.mat
Binary files differ
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-15 15:49:13 +0000