From 1e9d5bc4ad830ff6ddb3cd9fbb8239d7a93b43be Mon Sep 17 00:00:00 2001 From: dkazanc Date: Fri, 13 Dec 2019 17:57:21 +0000 Subject: initial modification of NDF --- Readme.md | 14 ++-- src/Core/regularisers_CPU/Diffusion_core.c | 88 ++++++++++++++++---------- src/Core/regularisers_GPU/NonlDiff_GPU_core.cu | 48 ++++++++++++-- 3 files changed, 104 insertions(+), 46 deletions(-) diff --git a/Readme.md b/Readme.md index 13e182b..077d2f2 100644 --- a/Readme.md +++ b/Readme.md @@ -101,13 +101,13 @@ conda install ccpi-regulariser -c ccpi -c conda-forge #### Python (conda-build) ``` - export CIL_VERSION=19.10 (Unix) / set CIL_VERSION=19.10 (Windows) - conda build recipe/ --numpy 1.15 --python 3.7 - conda install ccpi-regulariser=${CIL_VERSION} --use-local --force-reinstall # doesn't work? - conda install -c file://${CONDA_PREFIX}/conda-bld/ ccpi-regulariser=${CIL_VERSION} --force-reinstall # try this one - cd demos/ - python demo_cpu_regularisers.py # to run CPU demo - python demo_gpu_regularisers.py # to run GPU demo +export CIL_VERSION=19.10.1 (Unix) / set CIL_VERSION=19.10 (Windows) +conda build recipe/ --numpy 1.15 --python 3.7 +conda install ccpi-regulariser=${CIL_VERSION} --use-local --force-reinstall # doesn't work? +conda install -c file://${CONDA_PREFIX}/conda-bld/ ccpi-regulariser=${CIL_VERSION} --force-reinstall # try this one +cd demos/ +python demo_cpu_regularisers.py # to run CPU demo +python demo_gpu_regularisers.py # to run GPU demo ``` #### Python build diff --git a/src/Core/regularisers_CPU/Diffusion_core.c b/src/Core/regularisers_CPU/Diffusion_core.c index f0039c7..1a05c39 100644 --- a/src/Core/regularisers_CPU/Diffusion_core.c +++ b/src/Core/regularisers_CPU/Diffusion_core.c @@ -38,7 +38,7 @@ int signNDFc(float x) { * 3. Edge-preserving parameter (sigma), when sigma equals to zero nonlinear diffusion -> linear diffusion * 4. Number of iterations, for explicit scheme >= 150 is recommended * 5. tau - time-marching step for explicit scheme - * 6. Penalty type: 1 - Huber, 2 - Perona-Malik, 3 - Tukey Biweight + * 6. Penalty type: 1 - Huber, 2 - Perona-Malik, 3 - Tukey Biweight, 4 - Threshold-constrained Linear * 7. eplsilon - tolerance constant * * Output: @@ -60,14 +60,14 @@ float Diffusion_CPU_main(float *Input, float *Output, float *infovector, float l re = 0.0f; re1 = 0.0f; int count = 0; DimTotal = (long)(dimX*dimY*dimZ); - + if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float)); - + /* copy into output */ copyIm(Input, Output, (long)(dimX), (long)(dimY), (long)(dimZ)); - + for(i=0; i < iterationsNumb; i++) { - + if ((epsil != 0.0f) && (i % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ)); if (dimZ == 1) { /* running 2D diffusion iterations */ @@ -93,7 +93,7 @@ float Diffusion_CPU_main(float *Input, float *Output, float *infovector, float l if (count > 3) break; } } - + free(Output_prev); /*adding info into info_vector */ infovector[0] = (float)(i); /*iterations number (if stopped earlier based on tolerance)*/ @@ -109,7 +109,7 @@ float LinearDiff2D(float *Input, float *Output, float lambdaPar, float tau, long { long i,j,i1,i2,j1,j2,index; float e,w,n,s,e1,w1,n1,s1; - + #pragma omp parallel for shared(Input) private(index,i,j,i1,i2,j1,j2,e,w,n,s,e1,w1,n1,s1) for(j=0; j 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; } @@ -195,8 +195,18 @@ float NonLinearDiff2D(float *Input, float *Output, float lambdaPar, float sigmaP if (fabs(s1) <= sigmaPar) s1 = s1*powf((1.0f - powf((s1/sigmaPar),2)), 2); else s1 = 0.0f; } + else if (penaltytype == 4) { + /* Threshold-constrained linear diffusion + This means that the linear diffusion will be performed on pixels with + absolute difference less than the threshold. + */ + if (fabs(e1) > sigmaPar) e1 = 0.0f; + if (fabs(w1) > sigmaPar) w1 = 0.0f; + if (fabs(n1) > sigmaPar) n1 = 0.0f; + if (fabs(s1) > sigmaPar) s1 = 0.0f; + } else { - printf("%s \n", "No penalty function selected! Use 1,2 or 3."); + printf("%s \n", "No penalty function selected! Use 1,2,3 or 4."); break; } Output[index] += tau*(lambdaPar*(e1 + w1 + n1 + s1) - (Output[index] - Input[index])); @@ -211,7 +221,7 @@ float LinearDiff3D(float *Input, float *Output, float lambdaPar, float tau, long { long 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) 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 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; } @@ -322,11 +332,23 @@ float NonLinearDiff3D(float *Input, float *Output, float lambdaPar, float sigmaP if (fabs(d1) <= sigmaPar) d1 = d1*powf((1.0f - powf((d1/sigmaPar),2)), 2); else d1 = 0.0f; } + else if (penaltytype == 4) { + /* Threshold-constrained linear diffusion + This means that the linear diffusion will be performed on pixels with + absolute difference less than the threshold. + */ + if (fabs(e1) > sigmaPar) e1 = 0.0f; + if (fabs(w1) > sigmaPar) w1 = 0.0f; + if (fabs(n1) > sigmaPar) n1 = 0.0f; + if (fabs(s1) > sigmaPar) s1 = 0.0f; + if (fabs(u1) > sigmaPar) u1 = 0.0f; + if (fabs(d1) > sigmaPar) d1 = 0.0f; + } else { - printf("%s \n", "No penalty function selected! Use 1,2 or 3."); + printf("%s \n", "No penalty function selected! Use 1,2,3 or 4."); break; } - + Output[index] += tau*(lambdaPar*(e1 + w1 + n1 + s1 + u1 + d1) - (Output[index] - Input[index])); }}} return *Output; diff --git a/src/Core/regularisers_GPU/NonlDiff_GPU_core.cu b/src/Core/regularisers_GPU/NonlDiff_GPU_core.cu index de9abd4..a8876cd 100644 --- a/src/Core/regularisers_GPU/NonlDiff_GPU_core.cu +++ b/src/Core/regularisers_GPU/NonlDiff_GPU_core.cu @@ -32,7 +32,7 @@ limitations under the License. * 3. Edge-preserving parameter (sigma), when sigma equals to zero nonlinear diffusion -> linear diffusion * 4. Number of iterations, for explicit scheme >= 150 is recommended * 5. tau - time-marching step for explicit scheme - * 6. Penalty type: 1 - Huber, 2 - Perona-Malik, 3 - Tukey Biweight + * 6. Penalty type: 1 - Huber, 2 - Perona-Malik, 3 - Tukey Biweight, 4 - Threshold-constrained Linear * 7. eplsilon: tolerance constant * Output: @@ -108,8 +108,8 @@ __global__ void LinearDiff2D_kernel(float *Input, float *Output, float lambdaPar if ((i >= 0) && (i < N) && (j >= 0) && (j < M)) { /* boundary conditions (Neumann reflections) */ - i1 = i+1; if (i1 == N) i1 = i-1; - i2 = i-1; if (i2 < 0) i2 = i+1; + i1 = i+1; if (i1 == N) i1 = i-1; + i2 = i-1; if (i2 < 0) i2 = i+1; j1 = j+1; if (j1 == M) j1 = j-1; j2 = j-1; if (j2 < 0) j2 = j+1; @@ -155,7 +155,32 @@ __global__ void LinearDiff2D_kernel(float *Input, float *Output, float lambdaPar if (abs(s1) <= sigmaPar) s1 = s1*pow((1.0f - pow((s1/sigmaPar),2)), 2); else s1 = 0.0f; } - else printf("%s \n", "No penalty function selected! Use 1,2 or 3."); + else if (penaltytype == 4) { + /* Threshold-constrained linear diffusion + This means that the linear diffusion will be performed on pixels with + absolute difference less than the threshold. + */ + if (abs(e1) <= 1.5f*sigmaPar) { + if (abs(e1) > sigmaPar) e1 = signNDF(e1); + else e1 = e1/sigmaPar;} + else e1 = 0.0f; + + if (abs(w1) <= 1.5f*sigmaPar) { + if (abs(w1) > sigmaPar) w1 = signNDF(w1); + else w1 = w1/sigmaPar;} + else w1 = 0.0f; + + if (abs(n1) <= 1.5f*sigmaPar) { + if (abs(n1) > sigmaPar) n1 = signNDF(n1); + else n1 = n1/sigmaPar; } + else n1 = 0.0f; + + if (abs(s1) <= 1.5f*sigmaPar) { + if (abs(s1) > sigmaPar) s1 = signNDF(s1); + else s1 = s1/sigmaPar; } + else s1 = 0.0f; + } + else printf("%s \n", "No penalty function selected! Use 1,2,3 or 4."); Output[index] += tau*(lambdaPar*(e1 + w1 + n1 + s1) - (Output[index] - Input[index])); } @@ -281,8 +306,19 @@ __global__ void NonLinearDiff3D_kernel(float *Input, float *Output, float lambda if (abs(d1) <= sigmaPar) d1 = d1*pow((1.0f - pow((d1/sigmaPar),2)), 2); else d1 = 0.0f; } - else printf("%s \n", "No penalty function selected! Use 1,2 or 3."); - + else if (penaltytype == 4) { + /* Threshold-constrained linear diffusion + This means that the linear diffusion will be performed on pixels with + absolute difference less than the threshold. + */ + if (abs(e1) > sigmaPar) e1 = 0.0f; + if (abs(w1) > sigmaPar) w1 = 0.0f; + if (abs(n1) > sigmaPar) n1 = 0.0f; + if (abs(s1) > sigmaPar) s1 = 0.0f; + if (abs(u1) > sigmaPar) u1 = 0.0f; + if (abs(d1) > sigmaPar) d1 = 0.0f; + } + else printf("%s \n", "No penalty function selected! Use 1,2,3 or 4."); Output[index] += tau*(lambdaPar*(e1 + w1 + n1 + s1 + u1 + d1) - (Output[index] - Input[index])); } } -- cgit v1.2.3