diff options
author | Daniil Kazantsev <dkazanc@hotmail.com> | 2018-12-02 19:01:42 +0000 |
---|---|---|
committer | Daniil Kazantsev <dkazanc@hotmail.com> | 2018-12-02 19:01:42 +0000 |
commit | a48c9e69e941ec4046aca9d5d6ec453b9e9debdc (patch) | |
tree | f62cbc2b1d51aff9aaff14e1675f932f1922dde8 | |
parent | d252fcf6889855bb276cf6f9bf516e61910c064f (diff) | |
download | regularization-a48c9e69e941ec4046aca9d5d6ec453b9e9debdc.tar.gz regularization-a48c9e69e941ec4046aca9d5d6ec453b9e9debdc.tar.bz2 regularization-a48c9e69e941ec4046aca9d5d6ec453b9e9debdc.tar.xz regularization-a48c9e69e941ec4046aca9d5d6ec453b9e9debdc.zip |
cythonised nltv and updated demo, readme, bash run added
18 files changed, 327 insertions, 618 deletions
diff --git a/Core/regularisers_CPU/Nonlocal_TV_core.c b/Core/regularisers_CPU/Nonlocal_TV_core.c index d327dd5..c4c9118 100644 --- a/Core/regularisers_CPU/Nonlocal_TV_core.c +++ b/Core/regularisers_CPU/Nonlocal_TV_core.c @@ -42,36 +42,37 @@ */ /*****************************************************************************/ -float Nonlocal_TV_CPU_main(float *A_orig, float *Output, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, int dimX, int dimY, int dimZ, int NumNeighb, float lambda, int IterNumb) +float Nonlocal_TV_CPU_main(float *A_orig, float *Output, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, int dimX, int dimY, int dimZ, int NumNeighb, float lambdaReg, int IterNumb) { long i, j, k; int iter; + lambdaReg = 1.0f/lambdaReg; /*****2D INPUT *****/ if (dimZ == 0) { - copyIm(A_orig, Output, (long)(dimX), (long)(dimY), 1l); + copyIm(A_orig, Output, (long)(dimX), (long)(dimY), 1l); /* for each pixel store indeces of the most similar neighbours (patches) */ for(iter=0; iter<IterNumb; iter++) { #pragma omp parallel for shared (A_orig, Output, Weights, H_i, H_j, iter) private(i,j) for(i=0; i<(long)(dimX); i++) { for(j=0; j<(long)(dimY); j++) { - /* NLM_H1_2D(Output, A_orig, H_i, H_j, Weights, i, j, dimX, dimY, NumNeighb, lambda); */ /* NLM - H1 penalty */ - NLM_TV_2D(Output, A_orig, H_i, H_j, Weights, i, j, (long)(dimX), (long)(dimY), NumNeighb, lambda); /* NLM - TV penalty */ - }} + /*NLM_H1_2D(Output, A_orig, H_i, H_j, Weights, i, j, (long)(dimX), (long)(dimY), NumNeighb, lambdaReg);*/ /* NLM - H1 penalty */ + NLM_TV_2D(Output, A_orig, H_i, H_j, Weights, i, j, (long)(dimX), (long)(dimY), NumNeighb, lambdaReg); /* NLM - TV penalty */ + }} } } else { /*****3D INPUT *****/ - copyIm(A_orig, Output, (long)(dimX), (long)(dimY), (long)(dimZ)); + copyIm(A_orig, Output, (long)(dimX), (long)(dimY), (long)(dimZ)); /* for each pixel store indeces of the most similar neighbours (patches) */ for(iter=0; iter<IterNumb; iter++) { #pragma omp parallel for shared (A_orig, Output, Weights, H_i, H_j, H_k, iter) private(i,j,k) for(i=0; i<(long)(dimX); i++) { for(j=0; j<(long)(dimY); j++) { for(k=0; k<(long)(dimZ); k++) { - /* NLM_H1_3D(Output, A_orig, H_i, H_j, H_k, Weights, i, j, k, dimX, dimY, dimZ, NumNeighb, lambda); */ /* NLM - H1 penalty */ - NLM_TV_3D(Output, A_orig, H_i, H_j, H_k, Weights, i, j, k, (long)(dimX), (long)(dimY), (long)(dimZ), NumNeighb, lambda); /* NLM - TV penalty */ + /* NLM_H1_3D(Output, A_orig, H_i, H_j, H_k, Weights, i, j, k, dimX, dimY, dimZ, NumNeighb, lambdaReg); */ /* NLM - H1 penalty */ + NLM_TV_3D(Output, A_orig, H_i, H_j, H_k, Weights, i, j, k, (long)(dimX), (long)(dimY), (long)(dimZ), NumNeighb, lambdaReg); /* NLM - TV penalty */ }}} } } @@ -79,23 +80,24 @@ float Nonlocal_TV_CPU_main(float *A_orig, float *Output, unsigned short *H_i, un } /***********<<<<Main Function for NLM - H1 penalty>>>>**********/ -float NLM_H1_2D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimX, long dimY, int NumNeighb, float lambda) +float NLM_H1_2D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimX, long dimY, int NumNeighb, float lambdaReg) { - long x, i1, j1, index; + long x, i1, j1, index, index_m; float value = 0.0f, normweight = 0.0f; + index_m = j*dimX+i; for(x=0; x < NumNeighb; x++) { index = (dimX*dimY*x) + j*dimX+i; i1 = H_i[index]; j1 = H_j[index]; value += A[j1*dimX+i1]*Weights[index]; normweight += Weights[index]; - } - A[j*dimX+i] = (lambda*A_orig[j*dimX+i] + value)/(lambda + normweight); + } + A[index_m] = (lambdaReg*A_orig[index_m] + value)/(lambdaReg + normweight); return *A; } /*3D version*/ -float NLM_H1_3D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long i, long j, long k, long dimX, long dimY, long dimZ, int NumNeighb, float lambda) +float NLM_H1_3D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long i, long j, long k, long dimX, long dimY, long dimZ, int NumNeighb, float lambdaReg) { long x, i1, j1, k1, index; float value = 0.0f, normweight = 0.0f; @@ -108,39 +110,41 @@ float NLM_H1_3D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_ value += A[(dimX*dimY*k1) + j1*dimX+i1]*Weights[index]; normweight += Weights[index]; } - A[(dimX*dimY*k) + j*dimX+i] = (lambda*A_orig[(dimX*dimY*k) + j*dimX+i] + value)/(lambda + normweight); + A[(dimX*dimY*k) + j*dimX+i] = (lambdaReg*A_orig[(dimX*dimY*k) + j*dimX+i] + value)/(lambdaReg + normweight); return *A; } /***********<<<<Main Function for NLM - TV penalty>>>>**********/ -float NLM_TV_2D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimX, long dimY, int NumNeighb, float lambda) +float NLM_TV_2D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimX, long dimY, int NumNeighb, float lambdaReg) { - long x, i1, j1, index; + long x, i1, j1, index, index_m; float value = 0.0f, normweight = 0.0f, NLgrad_magn = 0.0f, NLCoeff; + index_m = j*dimX+i; + for(x=0; x < NumNeighb; x++) { - index = (dimX*dimY*x) + j*dimX+i; + index = (dimX*dimY*x) + j*dimX+i; /*c*/ i1 = H_i[index]; j1 = H_j[index]; - NLgrad_magn += powf((A[j1*dimX+i1] - A[j*dimX+i]),2)*Weights[index]; + NLgrad_magn += powf((A[j1*dimX+i1] - A[index_m]),2)*Weights[index]; } NLgrad_magn = sqrtf(NLgrad_magn); /*Non Local Gradients Magnitude */ NLCoeff = 2.0f*(1.0f/(NLgrad_magn + EPS)); for(x=0; x < NumNeighb; x++) { - index = (dimX*dimY*x) + j*dimX+i; + index = (dimX*dimY*x) + j*dimX+i; /*c*/ i1 = H_i[index]; j1 = H_j[index]; value += A[j1*dimX+i1]*NLCoeff*Weights[index]; normweight += Weights[index]*NLCoeff; } - A[j*dimX+i] = (lambda*A_orig[j*dimX+i] + value)/(lambda + normweight); + A[index_m] = (lambdaReg*A_orig[index_m] + value)/(lambdaReg + normweight); return *A; } /*3D version*/ -float NLM_TV_3D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long i, long j, long k, long dimX, long dimY, long dimZ, int NumNeighb, float lambda) +float NLM_TV_3D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long i, long j, long k, long dimX, long dimY, long dimZ, int NumNeighb, float lambdaReg) { long x, i1, j1, k1, index; float value = 0.0f, normweight = 0.0f, NLgrad_magn = 0.0f, NLCoeff; @@ -164,6 +168,6 @@ float NLM_TV_3D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_ value += A[(dimX*dimY*k1) + j1*dimX+i1]*NLCoeff*Weights[index]; normweight += Weights[index]*NLCoeff; } - A[(dimX*dimY*k) + j*dimX+i] = (lambda*A_orig[(dimX*dimY*k) + j*dimX+i] + value)/(lambda + normweight); + A[(dimX*dimY*k) + j*dimX+i] = (lambdaReg*A_orig[(dimX*dimY*k) + j*dimX+i] + value)/(lambdaReg + normweight); return *A; } diff --git a/Core/regularisers_CPU/Nonlocal_TV_core.h b/Core/regularisers_CPU/Nonlocal_TV_core.h index 5b6963e..6d55101 100644 --- a/Core/regularisers_CPU/Nonlocal_TV_core.h +++ b/Core/regularisers_CPU/Nonlocal_TV_core.h @@ -51,11 +51,11 @@ #ifdef __cplusplus extern "C" { #endif -CCPI_EXPORT float Nonlocal_TV_CPU_main(float *A_orig, float *Output, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, int dimX, int dimY, int dimZ, int NumNeighb, float lambda, int IterNumb); -CCPI_EXPORT float NLM_H1_2D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimX, long dimY, int NumNeighb, float lambda); -CCPI_EXPORT float NLM_TV_2D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimX, long dimY, int NumNeighb, float lambda); -CCPI_EXPORT float NLM_H1_3D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long i, long j, long k, long dimX, long dimY, long dimZ, int NumNeighb, float lambda); -CCPI_EXPORT float NLM_TV_3D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long i, long j, long k, long dimX, long dimY, long dimZ, int NumNeighb, float lambda); +CCPI_EXPORT float Nonlocal_TV_CPU_main(float *A_orig, float *Output, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, int dimX, int dimY, int dimZ, int NumNeighb, float lambdaReg, int IterNumb); +CCPI_EXPORT float NLM_H1_2D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimX, long dimY, int NumNeighb, float lambdaReg); +CCPI_EXPORT float NLM_TV_2D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimX, long dimY, int NumNeighb, float lambdaReg); +CCPI_EXPORT float NLM_H1_3D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long i, long j, long k, long dimX, long dimY, long dimZ, int NumNeighb, float lambdaReg); +CCPI_EXPORT float NLM_TV_3D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long i, long j, long k, long dimX, long dimY, long dimZ, int NumNeighb, float lambdaReg); #ifdef __cplusplus } #endif diff --git a/Core/regularisers_CPU/PatchSelect_core.c b/Core/regularisers_CPU/PatchSelect_core.c index efc5498..45cf8c8 100644 --- a/Core/regularisers_CPU/PatchSelect_core.c +++ b/Core/regularisers_CPU/PatchSelect_core.c @@ -46,15 +46,15 @@ /**************************************************/ -float PatchSelect_CPU_main(float *A, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long dimX, long dimY, long dimZ, int SearchWindow, int SimilarWin, int NumNeighb, float h) +float PatchSelect_CPU_main(float *A, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, int dimX, int dimY, int dimZ, int SearchWindow, int SimilarWin, int NumNeighb, float h, int switchM) { int counterG; long i, j, k; float *Eucl_Vec, h2; - h2 = h*h; + h2 = h*h; /****************2D INPUT ***************/ - if (dimZ == 0) { + if (dimZ == 0) { /* generate a 2D Gaussian kernel for NLM procedure */ Eucl_Vec = (float*) calloc ((2*SimilarWin+1)*(2*SimilarWin+1),sizeof(float)); counterG = 0; @@ -63,13 +63,21 @@ float PatchSelect_CPU_main(float *A, unsigned short *H_i, unsigned short *H_j, u Eucl_Vec[counterG] = (float)exp(-(pow(((float) i), 2) + pow(((float) j), 2))/(2*SimilarWin*SimilarWin)); counterG++; }} /*main neighb loop */ - /* for each pixel store indeces of the most similar neighbours (patches) */ + if (switchM == 1) { #pragma omp parallel for shared (A, Weights, H_i, H_j) private(i,j) - for(i=0; i<dimX; i++) { - for(j=0; j<dimY; j++) { - Indeces2D(A, H_i, H_j, Weights, (i), (j), (dimX), (dimY), Eucl_Vec, NumNeighb, SearchWindow, SimilarWin, h2); + for(i=0; i<(long)(dimX); i++) { + for(j=0; j<(long)(dimY); j++) { + Indeces2D_p(A, H_i, H_j, Weights, i, j, (long)(dimX), (long)(dimY), Eucl_Vec, NumNeighb, SearchWindow, SimilarWin, h2); + }} + } + else { +#pragma omp parallel for shared (A, Weights, H_i, H_j) private(i,j) + for(i=0; i<(long)(dimX); i++) { + for(j=0; j<(long)(dimY); j++) { + Indeces2D(A, H_i, H_j, Weights, i, j, (long)(dimX), (long)(dimY), Eucl_Vec, NumNeighb, SearchWindow, SimilarWin, h2); }} + } } else { /****************3D INPUT ***************/ @@ -84,18 +92,28 @@ float PatchSelect_CPU_main(float *A, unsigned short *H_i, unsigned short *H_j, u }}} /*main neighb loop */ /* for each voxel store indeces of the most similar neighbours (patches) */ + if (switchM == 1) { +#pragma omp parallel for shared (A, Weights, H_i, H_j, H_k) private(i,j,k) + for(i=0; i<dimX; i++) { + for(j=0; j<dimY; j++) { + for(k=0; k<dimZ; k++) { + Indeces3D(A, H_i, H_j, H_k, Weights, j, i, (k), (dimX), (dimY), (dimZ), Eucl_Vec, NumNeighb, SearchWindow, SimilarWin, h2); + }}} + } + else { #pragma omp parallel for shared (A, Weights, H_i, H_j, H_k) private(i,j,k) for(i=0; i<dimX; i++) { for(j=0; j<dimY; j++) { for(k=0; k<dimZ; k++) { Indeces3D(A, H_i, H_j, H_k, Weights, (i), (j), (k), (dimX), (dimY), (dimZ), Eucl_Vec, NumNeighb, SearchWindow, SimilarWin, h2); }}} + } } free(Eucl_Vec); return 1; } -float Indeces2D(float *Aorig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimY, long dimX, float *Eucl_Vec, int NumNeighb, int SearchWindow, int SimilarWin, float h2) +float Indeces2D(float *Aorig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimX, long dimY, float *Eucl_Vec, int NumNeighb, int SearchWindow, int SimilarWin, float h2) { long i1, j1, i_m, j_m, i_c, j_c, i2, j2, i3, j3, counter, x, y, index, sizeWin_tot, counterG; float *Weight_Vec, normsum, temp; @@ -167,6 +185,81 @@ float Indeces2D(float *Aorig, unsigned short *H_i, unsigned short *H_j, float *W return 1; } + +float Indeces2D_p(float *Aorig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimX, long dimY, float *Eucl_Vec, int NumNeighb, int SearchWindow, int SimilarWin, float h2) +{ + long i1, j1, i_m, j_m, i_c, j_c, i2, j2, i3, j3, counter, x, y, index, sizeWin_tot, counterG; + float *Weight_Vec, normsum, temp; + unsigned short *ind_i, *ind_j, temp_i, temp_j; + + sizeWin_tot = (2*SearchWindow + 1)*(2*SearchWindow + 1); + + Weight_Vec = (float*) calloc(sizeWin_tot, sizeof(float)); + ind_i = (unsigned short*) calloc(sizeWin_tot, sizeof(unsigned short)); + ind_j = (unsigned short*) calloc(sizeWin_tot, sizeof(unsigned short)); + + counter = 0; + for(i_m=-SearchWindow; i_m<=SearchWindow; i_m++) { + for(j_m=-SearchWindow; j_m<=SearchWindow; j_m++) { + i1 = i+i_m; + j1 = j+j_m; + if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) { + normsum = 0.0f; counterG = 0; + for(i_c=-SimilarWin; i_c<=SimilarWin; i_c++) { + for(j_c=-SimilarWin; j_c<=SimilarWin; j_c++) { + i2 = i1 + i_c; + j2 = j1 + j_c; + i3 = i + i_c; + j3 = j + j_c; + if (((i2 >= 0) && (i2 < dimX)) && ((j2 >= 0) && (j2 < dimY))) { + if (((i3 >= 0) && (i3 < dimX)) && ((j3 >= 0) && (j3 < dimY))) { + //normsum += Eucl_Vec[counterG]*pow(Aorig[j3*dimX + (i3)] - Aorig[j2*dimX + (i2)], 2); + normsum += Eucl_Vec[counterG]*pow(Aorig[i3*dimY + (j3)] - Aorig[i2*dimY + (j2)], 2); + counterG++; + }} + + }} + /* writing temporarily into vectors */ + if (normsum > EPS) { + Weight_Vec[counter] = expf(-normsum/h2); + ind_i[counter] = i1; + ind_j[counter] = j1; + counter++; + } + } + }} + /* do sorting to choose the most prominent weights [HIGH to LOW] */ + /* and re-arrange indeces accordingly */ + for (x = 0; x < counter; x++) { + for (y = 0; y < counter; y++) { + if (Weight_Vec[y] < Weight_Vec[x]) { + temp = Weight_Vec[y+1]; + temp_i = ind_i[y+1]; + temp_j = ind_j[y+1]; + Weight_Vec[y+1] = Weight_Vec[y]; + Weight_Vec[y] = temp; + ind_i[y+1] = ind_i[y]; + ind_i[y] = temp_i; + ind_j[y+1] = ind_j[y]; + ind_j[y] = temp_j; + }}} + /*sorting loop finished*/ + + /*now select the NumNeighb more prominent weights and store into arrays */ + for(x=0; x < NumNeighb; x++) { + //index = (dimX*dimY*x) + j*dimX+i; + index = (dimX*dimY*x) + i*dimY+j; + H_i[index] = ind_i[x]; + H_j[index] = ind_j[x]; + Weights[index] = Weight_Vec[x]; + } + + free(ind_i); + free(ind_j); + free(Weight_Vec); + return 1; +} + float Indeces3D(float *Aorig, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long i, long j, long k, long dimY, long dimX, long dimZ, float *Eucl_Vec, int NumNeighb, int SearchWindow, int SimilarWin, float h2) { long i1, j1, k1, i_m, j_m, k_m, i_c, j_c, k_c, i2, j2, k2, i3, j3, k3, counter, x, y, index, sizeWin_tot, counterG; diff --git a/Core/regularisers_CPU/PatchSelect_core.h b/Core/regularisers_CPU/PatchSelect_core.h index 43fce87..ddaa428 100644 --- a/Core/regularisers_CPU/PatchSelect_core.h +++ b/Core/regularisers_CPU/PatchSelect_core.h @@ -54,8 +54,9 @@ #ifdef __cplusplus extern "C" { #endif -CCPI_EXPORT float PatchSelect_CPU_main(float *A, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long dimX, long dimY, long dimZ, int SearchWindow, int SimilarWin, int NumNeighb, float h); -CCPI_EXPORT float Indeces2D(float *Aorig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimY, long dimX, float *Eucl_Vec, int NumNeighb, int SearchWindow, int SimilarWin, float h2); +CCPI_EXPORT float PatchSelect_CPU_main(float *A, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, int dimX, int dimY, int dimZ, int SearchWindow, int SimilarWin, int NumNeighb, float h, int switchM); +CCPI_EXPORT float Indeces2D(float *Aorig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimX, long dimY, float *Eucl_Vec, int NumNeighb, int SearchWindow, int SimilarWin, float h2); +CCPI_EXPORT float Indeces2D_p(float *Aorig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimX, long dimY, float *Eucl_Vec, int NumNeighb, int SearchWindow, int SimilarWin, float h2); CCPI_EXPORT float Indeces3D(float *Aorig, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long i, long j, long k, long dimY, long dimX, long dimZ, float *Eucl_Vec, int NumNeighb, int SearchWindow, int SimilarWin, float h2); #ifdef __cplusplus } @@ -83,7 +83,7 @@ conda install ccpi-regulariser -c ccpi -c conda-forge #### Python (conda-build) ``` - export CIL_VERSION=0.10.1 + export CIL_VERSION=0.10.2 conda build Wrappers/Python/conda-recipe --numpy 1.12 --python 3.5 conda install ccpi-regulariser=${CIL_VERSION} --use-local --force cd demos/ diff --git a/Wrappers/Matlab/demos/demoMatlab_denoise.m b/Wrappers/Matlab/demos/demoMatlab_denoise.m index 2cbdb56..54b8bac 100644 --- a/Wrappers/Matlab/demos/demoMatlab_denoise.m +++ b/Wrappers/Matlab/demos/demoMatlab_denoise.m @@ -136,21 +136,20 @@ figure; imshow(u_diff4, [0 1]); title('Diffusion 4thO denoised image (CPU)'); % figure; imshow(u_diff4_g, [0 1]); title('Diffusion 4thO denoised image (GPU)'); %% fprintf('Weights pre-calculation for Non-local TV (takes time on CPU) \n'); -SearchingWindow = 9; -PatchWindow = 3; +SearchingWindow = 7; +PatchWindow = 2; NeighboursNumber = 15; % the number of neibours to include -h = 0.25; % edge related parameter for NLM +h = 0.23; % edge related parameter for NLM [H_i, H_j, Weights] = PatchSelect(single(u0), SearchingWindow, PatchWindow, NeighboursNumber, h); %% fprintf('Denoise using Non-local Total Variation (CPU) \n'); -iter_nltv = 3; % number of nltv iterations +iter_nltv = 2; % number of nltv iterations lambda_nltv = 0.085; % regularisation parameter for nltv tic; u_nltv = Nonlocal_TV(single(u0), H_i, H_j, 0, Weights, lambda_nltv, iter_nltv); toc; rmse_nltv = (RMSE(u_nltv(:),Im(:))); fprintf('%s %f \n', 'RMSE error for Non-local Total Variation is:', rmse_nltv); figure; imshow(u_nltv, [0 1]); title('Non-local Total Variation denoised image (CPU)'); %% - %>>>>>>>>>>>>>> MULTI-CHANNEL priors <<<<<<<<<<<<<<< % fprintf('Denoise using the FGP-dTV model (CPU) \n'); diff --git a/Wrappers/Matlab/mex_compile/compileCPU_mex_Linux.m b/Wrappers/Matlab/mex_compile/compileCPU_mex_Linux.m index 49b5dfd..72a828e 100644 --- a/Wrappers/Matlab/mex_compile/compileCPU_mex_Linux.m +++ b/Wrappers/Matlab/mex_compile/compileCPU_mex_Linux.m @@ -72,6 +72,7 @@ mex NonlocalMarching_Inpaint.c NonlocalMarching_Inpaint_core.c utils.c CFLAGS="\ movefile('NonlocalMarching_Inpaint.mex*',Pathmove); delete SB_TV_core* ROF_TV_core* FGP_TV_core* FGP_dTV_core* TNV_core* utils* Diffusion_core* Diffus4th_order_core* TGV_core* LLT_ROF_core* CCPiDefines.h +delete PatchSelect_core* Nonlocal_TV_core* delete Diffusion_Inpaint_core* NonlocalMarching_Inpaint_core* fprintf('%s \n', '<<<<<<< Regularisers successfully compiled! >>>>>>>'); diff --git a/Wrappers/Matlab/mex_compile/compileCPU_mex_WINDOWS.m b/Wrappers/Matlab/mex_compile/compileCPU_mex_WINDOWS.m index 1b59dc2..6f7541c 100644 --- a/Wrappers/Matlab/mex_compile/compileCPU_mex_WINDOWS.m +++ b/Wrappers/Matlab/mex_compile/compileCPU_mex_WINDOWS.m @@ -60,6 +60,12 @@ fprintf('%s \n', 'Compiling ROF-LLT...'); mex LLT_ROF.c LLT_ROF_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99" movefile('LLT_ROF.mex*',Pathmove); +fprintf('%s \n', 'Compiling NonLocal-TV...'); +mex PatchSelect.c PatchSelect_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99" +mex Nonlocal_TV.c Nonlocal_TV_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99" +movefile('Nonlocal_TV.mex*',Pathmove); +movefile('PatchSelect.mex*',Pathmove); + fprintf('%s \n', 'Compiling additional tools...'); mex TV_energy.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99" movefile('TV_energy.mex*',Pathmove); @@ -73,9 +79,7 @@ fprintf('%s \n', 'Compiling Nonlocal marching method for inpaiting...'); mex NonlocalMarching_Inpaint.c NonlocalMarching_Inpaint_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99" movefile('NonlocalMarching_Inpaint.mex*',Pathmove); -delete SB_TV_core* ROF_TV_core* FGP_TV_core* FGP_dTV_core* TNV_core* utils* Diffusion_core* Diffus4th_order_core* TGV_core* CCPiDefines.h -delete Diffusion_Inpaint_core* NonlocalMarching_Inpaint_core* -fprintf('%s \n', 'Regularisers successfully compiled!'); + %% %%% The second approach to compile using TDM-GCC which follows this %%% discussion: @@ -105,15 +109,24 @@ fprintf('%s \n', 'Regularisers successfully compiled!'); % movefile('TGV.mex*',Pathmove); % mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" LLT_ROF.c LLT_ROF_core.c utils.c % movefile('LLT_ROF.mex*',Pathmove); +% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" PatchSelect.c PatchSelect_core.c utils.c +% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" Nonlocal_TV.c Nonlocal_TV_core.c utils.c +% movefile('Nonlocal_TV.mex*',Pathmove); +% movefile('PatchSelect.mex*',Pathmove); % mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" TV_energy.c utils.c % movefile('TV_energy.mex*',Pathmove); % mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" NonlDiff_Inp.c Diffusion_Inpaint_core.c utils.c % movefile('NonlDiff_Inp.mex*',Pathmove); % mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" NonlocalMarching_Inpaint.c NonlocalMarching_Inpaint_core.c utils.c % movefile('NonlocalMarching_Inpaint.mex*',Pathmove); -% delete SB_TV_core* ROF_TV_core* FGP_TV_core* FGP_dTV_core* TNV_core* utils* Diffusion_core* Diffus4th_order_core* TGV_core* CCPiDefines.h -% delete Diffusion_Inpaint_core* NonlocalMarching_Inpaint_core* -% fprintf('%s \n', 'Regularisers successfully compiled!'); + + +delete SB_TV_core* ROF_TV_core* FGP_TV_core* FGP_dTV_core* TNV_core* utils* Diffusion_core* Diffus4th_order_core* TGV_core* CCPiDefines.h +delete PatchSelect_core* Nonlocal_TV_core* +delete Diffusion_Inpaint_core* NonlocalMarching_Inpaint_core* +fprintf('%s \n', 'Regularisers successfully compiled!'); + + %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/Wrappers/Matlab/mex_compile/regularisers_CPU/Nonlocal_TV.c b/Wrappers/Matlab/mex_compile/regularisers_CPU/Nonlocal_TV.c index dea343c..014c0a0 100644 --- a/Wrappers/Matlab/mex_compile/regularisers_CPU/Nonlocal_TV.c +++ b/Wrappers/Matlab/mex_compile/regularisers_CPU/Nonlocal_TV.c @@ -68,7 +68,6 @@ void mexFunction( lambda = (float) mxGetScalar(prhs[5]); /* regularisation parameter */ IterNumb = (int) mxGetScalar(prhs[6]); /* the number of iterations */ - lambda = 1.0f/lambda; dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; /*****2D INPUT *****/ @@ -81,9 +80,9 @@ void mexFunction( /****************************************************/ if (number_of_dims == 3) { NumNeighb = dim_array2[3]; - Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); } /* run the main function here */ - Nonlocal_TV_CPU_main(A_orig, Output, H_i, H_j, H_k, Weights, dimX, dimY, dimZ, NumNeighb, lambda, IterNumb); + Nonlocal_TV_CPU_main(A_orig, Output, H_i, H_j, H_k, Weights, dimX, dimY, dimZ, NumNeighb, lambda, IterNumb); } diff --git a/Wrappers/Matlab/mex_compile/regularisers_CPU/Nonlocal_TV_core.c b/Wrappers/Matlab/mex_compile/regularisers_CPU/Nonlocal_TV_core.c deleted file mode 100644 index d327dd5..0000000 --- a/Wrappers/Matlab/mex_compile/regularisers_CPU/Nonlocal_TV_core.c +++ /dev/null @@ -1,169 +0,0 @@ -/* - * This work is part of the Core Imaging Library developed by - * Visual Analytics and Imaging System Group of the Science Technology - * Facilities Council, STFC and Diamond Light Source Ltd. - * - * Copyright 2017 Daniil Kazantsev - * Copyright 2017 Srikanth Nagella, Edoardo Pasca - * Copyright 2018 Diamond Light Source Ltd. - * - * 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 "Nonlocal_TV_core.h" - -/* C-OMP implementation of non-local regulariser - * Weights and associated indices must be given as an input. - * Gauss-Seidel fixed point iteration requires ~ 3 iterations, so the main effort - * goes in pre-calculation of weights and selection of patches - * - * - * Input Parameters: - * 1. 2D/3D grayscale image/volume - * 2. AR_i - indeces of i neighbours - * 3. AR_j - indeces of j neighbours - * 4. AR_k - indeces of k neighbours (0 - for 2D case) - * 5. Weights_ij(k) - associated weights - * 6. regularisation parameter - * 7. iterations number - - * Output: - * 1. denoised image/volume - * Elmoataz, Abderrahim, Olivier Lezoray, and Sébastien Bougleux. "Nonlocal discrete regularization on weighted graphs: a framework for image and manifold processing." IEEE Trans. Image Processing 17, no. 7 (2008): 1047-1060. - - */ -/*****************************************************************************/ - -float Nonlocal_TV_CPU_main(float *A_orig, float *Output, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, int dimX, int dimY, int dimZ, int NumNeighb, float lambda, int IterNumb) -{ - - long i, j, k; - int iter; - - /*****2D INPUT *****/ - if (dimZ == 0) { - copyIm(A_orig, Output, (long)(dimX), (long)(dimY), 1l); - /* for each pixel store indeces of the most similar neighbours (patches) */ - for(iter=0; iter<IterNumb; iter++) { -#pragma omp parallel for shared (A_orig, Output, Weights, H_i, H_j, iter) private(i,j) - for(i=0; i<(long)(dimX); i++) { - for(j=0; j<(long)(dimY); j++) { - /* NLM_H1_2D(Output, A_orig, H_i, H_j, Weights, i, j, dimX, dimY, NumNeighb, lambda); */ /* NLM - H1 penalty */ - NLM_TV_2D(Output, A_orig, H_i, H_j, Weights, i, j, (long)(dimX), (long)(dimY), NumNeighb, lambda); /* NLM - TV penalty */ - }} - } - } - else { - /*****3D INPUT *****/ - copyIm(A_orig, Output, (long)(dimX), (long)(dimY), (long)(dimZ)); - /* for each pixel store indeces of the most similar neighbours (patches) */ - for(iter=0; iter<IterNumb; iter++) { -#pragma omp parallel for shared (A_orig, Output, Weights, H_i, H_j, H_k, iter) private(i,j,k) - for(i=0; i<(long)(dimX); i++) { - for(j=0; j<(long)(dimY); j++) { - for(k=0; k<(long)(dimZ); k++) { - /* NLM_H1_3D(Output, A_orig, H_i, H_j, H_k, Weights, i, j, k, dimX, dimY, dimZ, NumNeighb, lambda); */ /* NLM - H1 penalty */ - NLM_TV_3D(Output, A_orig, H_i, H_j, H_k, Weights, i, j, k, (long)(dimX), (long)(dimY), (long)(dimZ), NumNeighb, lambda); /* NLM - TV penalty */ - }}} - } - } - return *Output; -} - -/***********<<<<Main Function for NLM - H1 penalty>>>>**********/ -float NLM_H1_2D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimX, long dimY, int NumNeighb, float lambda) -{ - long x, i1, j1, index; - float value = 0.0f, normweight = 0.0f; - - for(x=0; x < NumNeighb; x++) { - index = (dimX*dimY*x) + j*dimX+i; - i1 = H_i[index]; - j1 = H_j[index]; - value += A[j1*dimX+i1]*Weights[index]; - normweight += Weights[index]; - } - A[j*dimX+i] = (lambda*A_orig[j*dimX+i] + value)/(lambda + normweight); - return *A; -} -/*3D version*/ -float NLM_H1_3D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long i, long j, long k, long dimX, long dimY, long dimZ, int NumNeighb, float lambda) -{ - long x, i1, j1, k1, index; - float value = 0.0f, normweight = 0.0f; - - for(x=0; x < NumNeighb; x++) { - index = dimX*dimY*dimZ*x + (dimX*dimY*k) + j*dimX+i; - i1 = H_i[index]; - j1 = H_j[index]; - k1 = H_k[index]; - value += A[(dimX*dimY*k1) + j1*dimX+i1]*Weights[index]; - normweight += Weights[index]; - } - A[(dimX*dimY*k) + j*dimX+i] = (lambda*A_orig[(dimX*dimY*k) + j*dimX+i] + value)/(lambda + normweight); - return *A; -} - - -/***********<<<<Main Function for NLM - TV penalty>>>>**********/ -float NLM_TV_2D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimX, long dimY, int NumNeighb, float lambda) -{ - long x, i1, j1, index; - float value = 0.0f, normweight = 0.0f, NLgrad_magn = 0.0f, NLCoeff; - - for(x=0; x < NumNeighb; x++) { - index = (dimX*dimY*x) + j*dimX+i; - i1 = H_i[index]; - j1 = H_j[index]; - NLgrad_magn += powf((A[j1*dimX+i1] - A[j*dimX+i]),2)*Weights[index]; - } - - NLgrad_magn = sqrtf(NLgrad_magn); /*Non Local Gradients Magnitude */ - NLCoeff = 2.0f*(1.0f/(NLgrad_magn + EPS)); - - for(x=0; x < NumNeighb; x++) { - index = (dimX*dimY*x) + j*dimX+i; - i1 = H_i[index]; - j1 = H_j[index]; - value += A[j1*dimX+i1]*NLCoeff*Weights[index]; - normweight += Weights[index]*NLCoeff; - } - A[j*dimX+i] = (lambda*A_orig[j*dimX+i] + value)/(lambda + normweight); - return *A; -} -/*3D version*/ -float NLM_TV_3D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long i, long j, long k, long dimX, long dimY, long dimZ, int NumNeighb, float lambda) -{ - long x, i1, j1, k1, index; - float value = 0.0f, normweight = 0.0f, NLgrad_magn = 0.0f, NLCoeff; - - for(x=0; x < NumNeighb; x++) { - index = dimX*dimY*dimZ*x + (dimX*dimY*k) + j*dimX+i; - i1 = H_i[index]; - j1 = H_j[index]; - k1 = H_k[index]; - NLgrad_magn += powf((A[(dimX*dimY*k1) + j1*dimX+i1] - A[(dimX*dimY*k1) + j*dimX+i]),2)*Weights[index]; - } - - NLgrad_magn = sqrtf(NLgrad_magn); /*Non Local Gradients Magnitude */ - NLCoeff = 2.0f*(1.0f/(NLgrad_magn + EPS)); - - for(x=0; x < NumNeighb; x++) { - index = dimX*dimY*dimZ*x + (dimX*dimY*k) + j*dimX+i; - i1 = H_i[index]; - j1 = H_j[index]; - k1 = H_k[index]; - value += A[(dimX*dimY*k1) + j1*dimX+i1]*NLCoeff*Weights[index]; - normweight += Weights[index]*NLCoeff; - } - A[(dimX*dimY*k) + j*dimX+i] = (lambda*A_orig[(dimX*dimY*k) + j*dimX+i] + value)/(lambda + normweight); - return *A; -} diff --git a/Wrappers/Matlab/mex_compile/regularisers_CPU/Nonlocal_TV_core.h b/Wrappers/Matlab/mex_compile/regularisers_CPU/Nonlocal_TV_core.h deleted file mode 100644 index 5b6963e..0000000 --- a/Wrappers/Matlab/mex_compile/regularisers_CPU/Nonlocal_TV_core.h +++ /dev/null @@ -1,61 +0,0 @@ -/* - * This work is part of the Core Imaging Library developed by - * Visual Analytics and Imaging System Group of the Science Technology - * Facilities Council, STFC and Diamond Light Source Ltd. - * - * Copyright 2017 Daniil Kazantsev - * Copyright 2017 Srikanth Nagella, Edoardo Pasca - * Copyright 2018 Diamond Light Source Ltd. - * - * 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" - -#define EPS 1.0000e-9 - -/* C-OMP implementation of non-local regulariser - * Weights and associated indices must be given as an input. - * Gauss-Seidel fixed point iteration requires ~ 3 iterations, so the main effort - * goes in pre-calculation of weights and selection of patches - * - * - * Input Parameters: - * 1. 2D/3D grayscale image/volume - * 2. AR_i - indeces of i neighbours - * 3. AR_j - indeces of j neighbours - * 4. AR_k - indeces of k neighbours (0 - for 2D case) - * 5. Weights_ij(k) - associated weights - * 6. regularisation parameter - * 7. iterations number - - * Output: - * 1. denoised image/volume - * Elmoataz, Abderrahim, Olivier Lezoray, and Sébastien Bougleux. "Nonlocal discrete regularization on weighted graphs: a framework for image and manifold processing." IEEE Trans. Image Processing 17, no. 7 (2008): 1047-1060. - */ - -#ifdef __cplusplus -extern "C" { -#endif -CCPI_EXPORT float Nonlocal_TV_CPU_main(float *A_orig, float *Output, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, int dimX, int dimY, int dimZ, int NumNeighb, float lambda, int IterNumb); -CCPI_EXPORT float NLM_H1_2D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimX, long dimY, int NumNeighb, float lambda); -CCPI_EXPORT float NLM_TV_2D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimX, long dimY, int NumNeighb, float lambda); -CCPI_EXPORT float NLM_H1_3D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long i, long j, long k, long dimX, long dimY, long dimZ, int NumNeighb, float lambda); -CCPI_EXPORT float NLM_TV_3D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long i, long j, long k, long dimX, long dimY, long dimZ, int NumNeighb, float lambda); -#ifdef __cplusplus -} -#endif diff --git a/Wrappers/Matlab/mex_compile/regularisers_CPU/PatchSelect.c b/Wrappers/Matlab/mex_compile/regularisers_CPU/PatchSelect.c index fdd9a97..f942539 100644 --- a/Wrappers/Matlab/mex_compile/regularisers_CPU/PatchSelect.c +++ b/Wrappers/Matlab/mex_compile/regularisers_CPU/PatchSelect.c @@ -87,6 +87,6 @@ void mexFunction( Weights = (float*)mxGetPr(plhs[3] = mxCreateNumericArray(4, dim_array3, mxSINGLE_CLASS, mxREAL)); } - PatchSelect_CPU_main(A, H_i, H_j, H_k, Weights, (long)(dimX), (long)(dimY), (long)(dimZ), SearchWindow, SimilarWin, NumNeighb, h); + PatchSelect_CPU_main(A, H_i, H_j, H_k, Weights, (long)(dimX), (long)(dimY), (long)(dimZ), SearchWindow, SimilarWin, NumNeighb, h, 0); } diff --git a/Wrappers/Matlab/mex_compile/regularisers_CPU/PatchSelect_core.c b/Wrappers/Matlab/mex_compile/regularisers_CPU/PatchSelect_core.c deleted file mode 100644 index efc5498..0000000 --- a/Wrappers/Matlab/mex_compile/regularisers_CPU/PatchSelect_core.c +++ /dev/null @@ -1,254 +0,0 @@ -/* - * This work is part of the Core Imaging Library developed by - * Visual Analytics and Imaging System Group of the Science Technology - * Facilities Council, STFC and Diamond Light Source Ltd. - * - * Copyright 2017 Daniil Kazantsev - * Copyright 2017 Srikanth Nagella, Edoardo Pasca - * Copyright 2018 Diamond Light Source Ltd. - * - * 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 "PatchSelect_core.h" - -/* C-OMP implementation of non-local weight pre-calculation for non-local priors - * Weights and associated indices are stored into pre-allocated arrays and passed - * to the regulariser - * - * - * Input Parameters: - * 1. 2D/3D grayscale image/volume - * 2. Searching window (half-size of the main bigger searching window, e.g. 11) - * 3. Similarity window (half-size of the patch window, e.g. 2) - * 4. The number of neighbours to take (the most prominent after sorting neighbours will be taken) - * 5. noise-related parameter to calculate non-local weights - * - * Output [2D]: - * 1. AR_i - indeces of i neighbours - * 2. AR_j - indeces of j neighbours - * 3. Weights_ij - associated weights - * - * Output [3D]: - * 1. AR_i - indeces of i neighbours - * 2. AR_j - indeces of j neighbours - * 3. AR_k - indeces of j neighbours - * 4. Weights_ijk - associated weights - */ - -/**************************************************/ - -float PatchSelect_CPU_main(float *A, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long dimX, long dimY, long dimZ, int SearchWindow, int SimilarWin, int NumNeighb, float h) -{ - int counterG; - long i, j, k; - float *Eucl_Vec, h2; - h2 = h*h; - - /****************2D INPUT ***************/ - if (dimZ == 0) { - /* generate a 2D Gaussian kernel for NLM procedure */ - Eucl_Vec = (float*) calloc ((2*SimilarWin+1)*(2*SimilarWin+1),sizeof(float)); - counterG = 0; - for(i=-SimilarWin; i<=SimilarWin; i++) { - for(j=-SimilarWin; j<=SimilarWin; j++) { - Eucl_Vec[counterG] = (float)exp(-(pow(((float) i), 2) + pow(((float) j), 2))/(2*SimilarWin*SimilarWin)); - counterG++; - }} /*main neighb loop */ - - /* for each pixel store indeces of the most similar neighbours (patches) */ -#pragma omp parallel for shared (A, Weights, H_i, H_j) private(i,j) - for(i=0; i<dimX; i++) { - for(j=0; j<dimY; j++) { - Indeces2D(A, H_i, H_j, Weights, (i), (j), (dimX), (dimY), Eucl_Vec, NumNeighb, SearchWindow, SimilarWin, h2); - }} - } - else { - /****************3D INPUT ***************/ - /* generate a 3D Gaussian kernel for NLM procedure */ - Eucl_Vec = (float*) calloc ((2*SimilarWin+1)*(2*SimilarWin+1)*(2*SimilarWin+1),sizeof(float)); - counterG = 0; - for(i=-SimilarWin; i<=SimilarWin; i++) { - for(j=-SimilarWin; j<=SimilarWin; j++) { - for(k=-SimilarWin; k<=SimilarWin; k++) { - Eucl_Vec[counterG] = (float)exp(-(pow(((float) i), 2) + pow(((float) j), 2) + pow(((float) k), 2))/(2*SimilarWin*SimilarWin*SimilarWin)); - counterG++; - }}} /*main neighb loop */ - - /* for each voxel store indeces of the most similar neighbours (patches) */ -#pragma omp parallel for shared (A, Weights, H_i, H_j, H_k) private(i,j,k) - for(i=0; i<dimX; i++) { - for(j=0; j<dimY; j++) { - for(k=0; k<dimZ; k++) { - Indeces3D(A, H_i, H_j, H_k, Weights, (i), (j), (k), (dimX), (dimY), (dimZ), Eucl_Vec, NumNeighb, SearchWindow, SimilarWin, h2); - }}} - } - free(Eucl_Vec); - return 1; -} - -float Indeces2D(float *Aorig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimY, long dimX, float *Eucl_Vec, int NumNeighb, int SearchWindow, int SimilarWin, float h2) -{ - long i1, j1, i_m, j_m, i_c, j_c, i2, j2, i3, j3, counter, x, y, index, sizeWin_tot, counterG; - float *Weight_Vec, normsum, temp; - unsigned short *ind_i, *ind_j, temp_i, temp_j; - - sizeWin_tot = (2*SearchWindow + 1)*(2*SearchWindow + 1); - - Weight_Vec = (float*) calloc(sizeWin_tot, sizeof(float)); - ind_i = (unsigned short*) calloc(sizeWin_tot, sizeof(unsigned short)); - ind_j = (unsigned short*) calloc(sizeWin_tot, sizeof(unsigned short)); - - counter = 0; - for(i_m=-SearchWindow; i_m<=SearchWindow; i_m++) { - for(j_m=-SearchWindow; j_m<=SearchWindow; j_m++) { - i1 = i+i_m; - j1 = j+j_m; - if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) { - normsum = 0.0f; counterG = 0; - for(i_c=-SimilarWin; i_c<=SimilarWin; i_c++) { - for(j_c=-SimilarWin; j_c<=SimilarWin; j_c++) { - i2 = i1 + i_c; - j2 = j1 + j_c; - i3 = i + i_c; - j3 = j + j_c; - if (((i2 >= 0) && (i2 < dimX)) && ((j2 >= 0) && (j2 < dimY))) { - if (((i3 >= 0) && (i3 < dimX)) && ((j3 >= 0) && (j3 < dimY))) { - normsum += Eucl_Vec[counterG]*pow(Aorig[j3*dimX + (i3)] - Aorig[j2*dimX + (i2)], 2); - counterG++; - }} - - }} - /* writing temporarily into vectors */ - if (normsum > EPS) { - Weight_Vec[counter] = expf(-normsum/h2); - ind_i[counter] = i1; - ind_j[counter] = j1; - counter++; - } - } - }} - /* do sorting to choose the most prominent weights [HIGH to LOW] */ - /* and re-arrange indeces accordingly */ - for (x = 0; x < counter; x++) { - for (y = 0; y < counter; y++) { - if (Weight_Vec[y] < Weight_Vec[x]) { - temp = Weight_Vec[y+1]; - temp_i = ind_i[y+1]; - temp_j = ind_j[y+1]; - Weight_Vec[y+1] = Weight_Vec[y]; - Weight_Vec[y] = temp; - ind_i[y+1] = ind_i[y]; - ind_i[y] = temp_i; - ind_j[y+1] = ind_j[y]; - ind_j[y] = temp_j; - }}} - /*sorting loop finished*/ - - /*now select the NumNeighb more prominent weights and store into arrays */ - for(x=0; x < NumNeighb; x++) { - index = (dimX*dimY*x) + j*dimX+i; - H_i[index] = ind_i[x]; - H_j[index] = ind_j[x]; - Weights[index] = Weight_Vec[x]; - } - - free(ind_i); - free(ind_j); - free(Weight_Vec); - return 1; -} - -float Indeces3D(float *Aorig, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long i, long j, long k, long dimY, long dimX, long dimZ, float *Eucl_Vec, int NumNeighb, int SearchWindow, int SimilarWin, float h2) -{ - long i1, j1, k1, i_m, j_m, k_m, i_c, j_c, k_c, i2, j2, k2, i3, j3, k3, counter, x, y, index, sizeWin_tot, counterG; - float *Weight_Vec, normsum, temp; - unsigned short *ind_i, *ind_j, *ind_k, temp_i, temp_j, temp_k; - - sizeWin_tot = (2*SearchWindow + 1)*(2*SearchWindow + 1)*(2*SearchWindow + 1); - - Weight_Vec = (float*) calloc(sizeWin_tot, sizeof(float)); - ind_i = (unsigned short*) calloc(sizeWin_tot, sizeof(unsigned short)); - ind_j = (unsigned short*) calloc(sizeWin_tot, sizeof(unsigned short)); - ind_k = (unsigned short*) calloc(sizeWin_tot, sizeof(unsigned short)); - - counter = 0l; - for(i_m=-SearchWindow; i_m<=SearchWindow; i_m++) { - for(j_m=-SearchWindow; j_m<=SearchWindow; j_m++) { - for(k_m=-SearchWindow; k_m<=SearchWindow; k_m++) { - k1 = k+k_m; - i1 = i+i_m; - j1 = j+j_m; - if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY)) && ((k1 >= 0) && (k1 < dimZ))) { - normsum = 0.0f; counterG = 0l; - for(i_c=-SimilarWin; i_c<=SimilarWin; i_c++) { - for(j_c=-SimilarWin; j_c<=SimilarWin; j_c++) { - for(k_c=-SimilarWin; k_c<=SimilarWin; k_c++) { - i2 = i1 + i_c; - j2 = j1 + j_c; - k2 = k1 + k_c; - i3 = i + i_c; - j3 = j + j_c; - k3 = k + k_c; - if (((i2 >= 0) && (i2 < dimX)) && ((j2 >= 0) && (j2 < dimY)) && ((k2 >= 0) && (k2 < dimZ))) { - if (((i3 >= 0) && (i3 < dimX)) && ((j3 >= 0) && (j3 < dimY)) && ((k3 >= 0) && (k3 < dimZ))) { - normsum += Eucl_Vec[counterG]*pow(Aorig[(dimX*dimY*k3) + j3*dimX + (i3)] - Aorig[(dimX*dimY*k2) + j2*dimX + (i2)], 2); - counterG++; - }} - }}} - /* writing temporarily into vectors */ - if (normsum > EPS) { - Weight_Vec[counter] = expf(-normsum/h2); - ind_i[counter] = i1; - ind_j[counter] = j1; - ind_k[counter] = k1; - counter ++; - } - } - }}} - /* do sorting to choose the most prominent weights [HIGH to LOW] */ - /* and re-arrange indeces accordingly */ - for (x = 0; x < counter; x++) { - for (y = 0; y < counter; y++) { - if (Weight_Vec[y] < Weight_Vec[x]) { - temp = Weight_Vec[y+1]; - temp_i = ind_i[y+1]; - temp_j = ind_j[y+1]; - temp_k = ind_k[y+1]; - Weight_Vec[y+1] = Weight_Vec[y]; - Weight_Vec[y] = temp; - ind_i[y+1] = ind_i[y]; - ind_i[y] = temp_i; - ind_j[y+1] = ind_j[y]; - ind_j[y] = temp_j; - ind_k[y+1] = ind_k[y]; - ind_k[y] = temp_k; - }}} - /*sorting loop finished*/ - - /*now select the NumNeighb more prominent weights and store into arrays */ - for(x=0; x < NumNeighb; x++) { - index = dimX*dimY*dimZ*x + (dimX*dimY*k) + j*dimX+i; - - H_i[index] = ind_i[x]; - H_j[index] = ind_j[x]; - H_k[index] = ind_k[x]; - - Weights[index] = Weight_Vec[x]; - } - - free(ind_i); - free(ind_j); - free(ind_k); - free(Weight_Vec); - return 1; -} - diff --git a/Wrappers/Matlab/mex_compile/regularisers_CPU/PatchSelect_core.h b/Wrappers/Matlab/mex_compile/regularisers_CPU/PatchSelect_core.h deleted file mode 100644 index 43fce87..0000000 --- a/Wrappers/Matlab/mex_compile/regularisers_CPU/PatchSelect_core.h +++ /dev/null @@ -1,62 +0,0 @@ -/* - * This work is part of the Core Imaging Library developed by - * Visual Analytics and Imaging System Group of the Science Technology - * Facilities Council, STFC and Diamond Light Source Ltd. - * - * Copyright 2017 Daniil Kazantsev - * Copyright 2017 Srikanth Nagella, Edoardo Pasca - * Copyright 2018 Diamond Light Source Ltd. - * - * 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" -#define EPS 1.0000e-12 - -/* C-OMP implementation of non-local weight pre-calculation for non-local priors - * Weights and associated indices are stored into pre-allocated arrays and passed - * to the regulariser - * - * - * Input Parameters: - * 1. 2D/3D grayscale image/volume - * 2. Searching window (half-size of the main bigger searching window, e.g. 11) - * 3. Similarity window (half-size of the patch window, e.g. 2) - * 4. The number of neighbours to take (the most prominent after sorting neighbours will be taken) - * 5. noise-related parameter to calculate non-local weights - * - * Output [2D]: - * 1. AR_i - indeces of i neighbours - * 2. AR_j - indeces of j neighbours - * 3. Weights_ij - associated weights - * - * Output [3D]: - * 1. AR_i - indeces of i neighbours - * 2. AR_j - indeces of j neighbours - * 3. AR_k - indeces of j neighbours - * 4. Weights_ijk - associated weights - */ -/*****************************************************************************/ -#ifdef __cplusplus -extern "C" { -#endif -CCPI_EXPORT float PatchSelect_CPU_main(float *A, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long dimX, long dimY, long dimZ, int SearchWindow, int SimilarWin, int NumNeighb, float h); -CCPI_EXPORT float Indeces2D(float *Aorig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimY, long dimX, float *Eucl_Vec, int NumNeighb, int SearchWindow, int SimilarWin, float h2); -CCPI_EXPORT float Indeces3D(float *Aorig, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long i, long j, long k, long dimY, long dimX, long dimZ, float *Eucl_Vec, int NumNeighb, int SearchWindow, int SimilarWin, float h2); -#ifdef __cplusplus -} -#endif diff --git a/Wrappers/Python/ccpi/filters/regularisers.py b/Wrappers/Python/ccpi/filters/regularisers.py index c3c3c7e..bf7e23c 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 import TV_ROF_CPU, TV_FGP_CPU, TV_SB_CPU, dTV_FGP_CPU, TNV_CPU, NDF_CPU, Diff4th_CPU, TGV_CPU, LLT_ROF_CPU, PATCHSEL_CPU +from ccpi.filters.cpu_regularisers import TV_ROF_CPU, TV_FGP_CPU, TV_SB_CPU, dTV_FGP_CPU, TNV_CPU, NDF_CPU, Diff4th_CPU, TGV_CPU, LLT_ROF_CPU, PATCHSEL_CPU, NLTV_CPU try: from ccpi.filters.gpu_regularisers import TV_ROF_GPU, TV_FGP_GPU, TV_SB_GPU, dTV_FGP_GPU, NDF_GPU, Diff4th_GPU, TGV_GPU, LLT_ROF_GPU gpu_enabled = True @@ -145,7 +145,7 @@ def DIFF4th(inputData, regularisation_parameter, edge_parameter, iterations, raise ValueError('Unknown device {0}. Expecting gpu or cpu'\ .format(device)) -def PatchSelect_CPU(inputData, searchwindow, patchwindow, neighbours, edge_parameter, device='cpu'): +def PatchSelect(inputData, searchwindow, patchwindow, neighbours, edge_parameter, device='cpu'): if device == 'cpu': return PATCHSEL_CPU(inputData, searchwindow, @@ -159,7 +159,16 @@ def PatchSelect_CPU(inputData, searchwindow, patchwindow, neighbours, edge_param raise ValueError ('GPU is not available') raise ValueError('Unknown device {0}. Expecting gpu or cpu'\ .format(device)) - + +def NLTV(inputData, H_i, H_j, H_k, Weights, regularisation_parameter, iterations): + return NLTV_CPU(inputData, + H_i, + H_j, + H_k, + Weights, + regularisation_parameter, + iterations) + def TGV(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst, device='cpu'): if device == 'cpu': diff --git a/Wrappers/Python/demos/demo_cpu_regularisers.py b/Wrappers/Python/demos/demo_cpu_regularisers.py index e99b271..31e4cad 100644 --- a/Wrappers/Python/demos/demo_cpu_regularisers.py +++ b/Wrappers/Python/demos/demo_cpu_regularisers.py @@ -13,6 +13,7 @@ import numpy as np import os import timeit from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, TNV, NDF, DIFF4th +from ccpi.filters.regularisers import PatchSelect, NLTV from qualitymetrics import rmse ############################################################################### def printParametersToString(pars): @@ -350,7 +351,7 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(ndf_cpu, cmap="gray") plt.title('{}'.format('CPU results')) - +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("___Anisotropic Diffusion 4th Order (2D)____") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") @@ -395,7 +396,71 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(diff4_cpu, cmap="gray") plt.title('{}'.format('CPU results')) +#%% +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("___Nonlocal patches pre-calculation____") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +# set parameters +pars = {'algorithm' : PatchSelect, \ + 'input' : u0,\ + 'searchwindow': 7, \ + 'patchwindow': 2,\ + 'neighbours' : 15 ,\ + 'edge_parameter':0.23} + +H_i, H_j, Weights = PatchSelect(pars['input'], + pars['searchwindow'], + pars['patchwindow'], + pars['neighbours'], + pars['edge_parameter'],'cpu') + +#%% +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("___Nonlocal Total Variation penalty____") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +## plot +fig = plt.figure() +plt.suptitle('Performance of NLTV regulariser using the CPU') +a=fig.add_subplot(1,2,1) +a.set_title('Noisy Image') +imgplot = plt.imshow(u0,cmap="gray") + +pars2 = {'algorithm' : NLTV, \ + 'input' : u0,\ + 'H_i': H_i, \ + 'H_j': H_j,\ + 'H_k' : 0,\ + 'Weights' : Weights,\ + 'regularisation_parameter': 0.085,\ + 'iterations': 2 + } +#%% +start_time = timeit.default_timer() +nltv_cpu = NLTV(pars2['input'], + pars2['H_i'], + pars2['H_j'], + pars2['H_k'], + pars2['Weights'], + pars2['regularisation_parameter'], + pars2['iterations']) + +rms = rmse(Im, nltv_cpu) +pars['rmse'] = rms + +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +a=fig.add_subplot(1,2,2) +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) +# place a text box in upper left in axes coords +a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(nltv_cpu, cmap="gray") +plt.title('{}'.format('CPU results')) + +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("_____________FGP-dTV (2D)__________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") @@ -447,7 +512,7 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, verticalalignment='top', bbox=props) imgplot = plt.imshow(fgp_dtv_cpu, cmap="gray") plt.title('{}'.format('CPU results')) - +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("__________Total nuclear Variation__________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") diff --git a/Wrappers/Python/src/cpu_regularisers.pyx b/Wrappers/Python/src/cpu_regularisers.pyx index b056bba..e51e6d8 100644 --- a/Wrappers/Python/src/cpu_regularisers.pyx +++ b/Wrappers/Python/src/cpu_regularisers.pyx @@ -27,10 +27,9 @@ cdef extern float Diffusion_CPU_main(float *Input, float *Output, float lambdaPa cdef extern float Diffus4th_CPU_main(float *Input, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int dimX, int dimY, int dimZ); 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 PatchSelect_CPU_main(float *Input, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long dimX, long dimY, long dimZ, int SearchWindow, int SimilarWin, int NumNeighb, float h); +cdef extern float PatchSelect_CPU_main(float *Input, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, int dimX, int dimY, int dimZ, int SearchWindow, int SimilarWin, int NumNeighb, float h, int switchM); cdef extern float Nonlocal_TV_CPU_main(float *A_orig, float *Output, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, int dimX, int dimY, int dimZ, int NumNeighb, float lambdaReg, int IterNumb); - 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 trigger, int dimX, int dimY, int dimZ); cdef extern float TV_energy2D(float *U, float *U0, float *E_val, float lambdaPar, int type, int dimX, int dimY); @@ -450,7 +449,6 @@ def Diff4th_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, return outputData - #****************************************************************# #***************Patch-based weights calculation******************# #****************************************************************# @@ -458,31 +456,85 @@ def PATCHSEL_CPU(inputData, searchwindow, patchwindow, neighbours, edge_paramete if inputData.ndim == 2: return PatchSel_2D(inputData, searchwindow, patchwindow, neighbours, edge_parameter) elif inputData.ndim == 3: - return 1 -# PatchSel_3D(inputData, searchwindow, patchwindow, neighbours, edge_parameter) -def PatchSel_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, + return PatchSel_3D(inputData, searchwindow, patchwindow, neighbours, edge_parameter) +def PatchSel_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, int searchwindow, int patchwindow, int neighbours, float edge_parameter): cdef long dims[3] - dims[0] = inputData.shape[0] - dims[1] = inputData.shape[1] - dims[2] = neighbours + dims[0] = neighbours + dims[1] = inputData.shape[0] + dims[2] = inputData.shape[1] + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] Weights = \ - np.zeros([dims[0],dims[1],dims[2]], dtype='float32') + np.zeros([dims[0], dims[1],dims[2]], dtype='float32') cdef np.ndarray[np.uint16_t, ndim=3, mode="c"] H_i = \ - np.zeros([dims[0],dims[1],dims[2]], dtype='uint16 ') + np.zeros([dims[0], dims[1],dims[2]], dtype='uint16') cdef np.ndarray[np.uint16_t, ndim=3, mode="c"] H_j = \ - np.zeros([dims[0],dims[1],dims[2]], dtype='uint16 ') + np.zeros([dims[0], dims[1],dims[2]], dtype='uint16') # Run patch-based weight selection function - PatchSelect_CPU_main(&inputData[0,0], &H_i[0,0,0], &H_j[0,0,0], &H_i[0,0,0], &Weights[0,0,0], dims[1], dims[0], 1, searchwindow, patchwindow, neighbours, edge_parameter) + PatchSelect_CPU_main(&inputData[0,0], &H_j[0,0,0], &H_i[0,0,0], &H_i[0,0,0], &Weights[0,0,0], dims[2], dims[1], 0, searchwindow, patchwindow, neighbours, edge_parameter, 1) return H_i, H_j, Weights - + +def PatchSel_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + int searchwindow, + int patchwindow, + int neighbours, + float edge_parameter): + cdef long dims[4] + dims[0] = inputData.shape[0] + dims[1] = inputData.shape[1] + dims[2] = inputData.shape[2] + dims[3] = neighbours + + cdef np.ndarray[np.float32_t, ndim=4, mode="c"] Weights = \ + np.zeros([dims[3],dims[0],dims[1],dims[2]], dtype='float32') + + cdef np.ndarray[np.uint16_t, ndim=4, mode="c"] H_i = \ + np.zeros([dims[3],dims[0],dims[1],dims[2]], dtype='uint16') + + cdef np.ndarray[np.uint16_t, ndim=4, mode="c"] H_j = \ + np.zeros([dims[3],dims[0],dims[1],dims[2]], dtype='uint16') + + cdef np.ndarray[np.uint16_t, ndim=4, mode="c"] H_k = \ + np.zeros([dims[3],dims[0],dims[1],dims[2]], dtype='uint16') + + # Run patch-based weight selection function + PatchSelect_CPU_main(&inputData[0,0,0], &H_i[0,0,0,0], &H_j[0,0,0,0], &H_k[0,0,0,0], &Weights[0,0,0,0], dims[2], dims[1], dims[0], searchwindow, patchwindow, neighbours, edge_parameter, 1) + return H_i, H_j, H_k, Weights + + +#****************************************************************# +#***************Non-local Total Variation******************# +#****************************************************************# +def NLTV_CPU(inputData, H_i, H_j, H_k, Weights, regularisation_parameter, iterations): + if inputData.ndim == 2: + return NLTV_2D(inputData, H_i, H_j, Weights, regularisation_parameter, iterations) + elif inputData.ndim == 3: + return 1 +def NLTV_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, + np.ndarray[np.uint16_t, ndim=3, mode="c"] H_i, + np.ndarray[np.uint16_t, ndim=3, mode="c"] H_j, + np.ndarray[np.float32_t, ndim=3, mode="c"] Weights, + float regularisation_parameter, + int iterations): + + cdef long dims[2] + dims[0] = inputData.shape[0] + dims[1] = inputData.shape[1] + neighbours = H_i.shape[0] + + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \ + np.zeros([dims[0],dims[1]], dtype='float32') + + # Run nonlocal TV regularisation + Nonlocal_TV_CPU_main(&inputData[0,0], &outputData[0,0], &H_i[0,0,0], &H_j[0,0,0], &H_i[0,0,0], &Weights[0,0,0], dims[1], dims[0], 0, neighbours, regularisation_parameter, iterations) + return outputData #*********************Inpainting WITH****************************# #***************Nonlinear (Isotropic) Diffusion******************# @@ -0,0 +1,19 @@ +#!/bin/bash +echo "Building CCPi-regularisation Toolkit using CMake" +# rm -r build +# Requires Cython, install it first: +# pip install cython +mkdir build +cd build/ +make clean +# install Python modules only without CUDA +cmake ../ -DBUILD_PYTHON_WRAPPER=ON -DBUILD_MATLAB_WRAPPER=OFF -DBUILD_CUDA=OFF -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install +# install Python modules only with CUDA +# cmake ../ -DBUILD_PYTHON_WRAPPER=ON -DBUILD_MATLAB_WRAPPER=OFF -DBUILD_CUDA=ON -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install +make install +# cp install/lib/libcilreg.so install/python/ccpi/filters +cd install/python +export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:../lib +# spyder +# one can also run Matlab in Linux as: +# PATH="/path/to/mex/:$PATH" LD_LIBRARY_PATH="/path/to/library:$LD_LIBRARY_PATH" matlab |