From 09eb48ffbb4ad699e2eefd25678e10dc59d6a177 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Tue, 1 May 2018 09:44:07 +0100 Subject: new inpainters --- Core/CMakeLists.txt | 4 +- Core/inpainters_CPU/Diffusion_Inpaint_core.c | 317 +++++++++++++++++++++ Core/inpainters_CPU/Diffusion_Inpaint_core.h | 60 ++++ .../inpainters_CPU/NonlocalMarching_Inpaint_core.c | 158 ++++++++++ .../inpainters_CPU/NonlocalMarching_Inpaint_core.h | 54 ++++ Core/regularisers_CPU/utils.c | 33 ++- Core/regularisers_CPU/utils.h | 5 +- Readme.md | 27 +- Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m | 9 +- Wrappers/Matlab/demos/demoMatlab_denoise.m | 9 +- Wrappers/Matlab/demos/demoMatlab_inpaint.m | 35 +++ Wrappers/Matlab/mex_compile/compileCPU_mex.m | 16 +- .../mex_compile/regularisers_CPU/NonlDiff_Inp.c | 101 +++++++ .../regularisers_CPU/NonlocalMarching_Inpaint.c | 82 ++++++ Wrappers/Python/ccpi/filters/regularisers.py | 7 +- Wrappers/Python/setup-regularisers.py.in | 3 +- Wrappers/Python/src/cpu_regularisers.pyx | 50 ++++ data/SinoInpaint.mat | Bin 0 -> 3303762 bytes 18 files changed, 937 insertions(+), 33 deletions(-) create mode 100644 Core/inpainters_CPU/Diffusion_Inpaint_core.c create mode 100644 Core/inpainters_CPU/Diffusion_Inpaint_core.h create mode 100644 Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c create mode 100644 Core/inpainters_CPU/NonlocalMarching_Inpaint_core.h create mode 100644 Wrappers/Matlab/demos/demoMatlab_inpaint.m create mode 100644 Wrappers/Matlab/mex_compile/regularisers_CPU/NonlDiff_Inp.c create mode 100644 Wrappers/Matlab/mex_compile/regularisers_CPU/NonlocalMarching_Inpaint.c create mode 100644 data/SinoInpaint.mat 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 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 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 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 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 +#include +#include +#include +#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= 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 +#include +#include +#include +#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 -/* 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 -//#include #include #include #include "CCPiDefines.h" -//#include #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 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 Binary files /dev/null and b/data/SinoInpaint.mat differ -- cgit v1.2.3