From 2d436f37be6029d57d7f876d4c7c378ee712a11e Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Wed, 10 Apr 2019 22:39:11 +0100 Subject: progress with bresenham --- src/Core/regularisers_CPU/DiffusionMASK_core.c | 210 ++++++++++++++++++++++--- src/Core/regularisers_CPU/DiffusionMASK_core.h | 7 +- src/Python/ccpi/filters/regularisers.py | 4 +- src/Python/src/cpu_regularisers.pyx | 10 +- 4 files changed, 207 insertions(+), 24 deletions(-) (limited to 'src') diff --git a/src/Core/regularisers_CPU/DiffusionMASK_core.c b/src/Core/regularisers_CPU/DiffusionMASK_core.c index eef173d..2af3cb6 100644 --- a/src/Core/regularisers_CPU/DiffusionMASK_core.c +++ b/src/Core/regularisers_CPU/DiffusionMASK_core.c @@ -36,7 +36,7 @@ int signNDF_m(float x) { * Input Parameters: * 1. Noisy image/volume * 2. MASK (in unsigned char format) - * 3. Diffusivity window (half-size of the searching window, e.g. 3) + * 3. Diffusivity window (half-size of the searching window, e.g. 3) * 4. lambda - regularization parameter * 5. Edge-preserving parameter (sigma), when sigma equals to zero nonlinear diffusion -> linear diffusion * 6. Number of iterations, for explicit scheme >= 150 is recommended @@ -53,9 +53,10 @@ int signNDF_m(float x) { * [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 DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, float *Output, float *infovector, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ) +float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *MASK_upd, float *Output, float *infovector, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ) { - int i,j,k,counterG; + long i,j,k; + int counterG; float sigmaPar2, *Output_prev=NULL, *Eucl_Vec; int DiffusWindow_tot; sigmaPar2 = sigmaPar/sqrt(2.0f); @@ -63,6 +64,7 @@ float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, float *Output, f float re, re1; re = 0.0f; re1 = 0.0f; int count = 0; + unsigned char *MASK_temp; DimTotal = (long)(dimX*dimY*dimZ); if (dimZ == 1) { @@ -85,21 +87,48 @@ float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, float *Output, f for(k=-DiffusWindow; k<=DiffusWindow; k++) { Eucl_Vec[counterG] = (float)expf(-(powf(((float) i), 2) + powf(((float) j), 2) + powf(((float) k), 2))/(2*DiffusWindow*DiffusWindow*DiffusWindow)); counterG++; - }}} /*main neighb loop */ + }}} /*main neighb loop */ } if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float)); - /* copy into output */ + MASK_temp = (unsigned char*) calloc (DimTotal,sizeof(unsigned char)); + + /* copy input into output */ copyIm(Input, Output, (long)(dimX), (long)(dimY), (long)(dimZ)); + /* copy given MASK */ + copyIm_unchar(MASK, MASK_upd, (long)(dimX), (long)(dimY), (long)(dimZ)); + /********************** PERFORM MASK PROCESSING ************************/ + if (dimZ == 1) { + #pragma omp parallel for shared(MASK,MASK_upd) private(i,j) + for(i=0; i continue */ + i = 162; j = 258; + Mask_update2D(MASK_temp, MASK_upd, i, j, DiffusWindow, (long)(dimX), (long)(dimY)); +// } + //}} + } + + /* start diffusivity iterations */ 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 */ - if (sigmaPar == 0.0f) LinearDiff_MASK2D(Input, MASK, Output, Eucl_Vec, DiffusWindow, lambdaPar, tau, (long)(dimX), (long)(dimY)); /* constrained linear diffusion */ - else NonLinearDiff_MASK2D(Input, MASK, Output, Eucl_Vec, DiffusWindow, lambdaPar, sigmaPar2, tau, penaltytype, (long)(dimX), (long)(dimY)); /* constrained nonlinear diffusion */ + /* running 2D diffusion iterations */ + //if (sigmaPar == 0.0f) LinearDiff_MASK2D(Input, MASK, Output, Eucl_Vec, DiffusWindow, lambdaPar, tau, (long)(dimX), (long)(dimY)); /* constrained linear diffusion */ + //else NonLinearDiff_MASK2D(Input, MASK, Output, Eucl_Vec, DiffusWindow, lambdaPar, sigmaPar2, tau, penaltytype, (long)(dimX), (long)(dimY)); /* constrained nonlinear diffusion */ } else { /* running 3D diffusion iterations */ @@ -122,6 +151,8 @@ float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, float *Output, f } free(Output_prev); + free(Eucl_Vec); + free(MASK_temp); /*adding info into info_vector */ infovector[0] = (float)(i); /*iterations number (if stopped earlier based on tolerance)*/ infovector[1] = re; /* reached tolerance */ @@ -132,6 +163,62 @@ float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, float *Output, f /********************************************************************/ /***************************2D Functions*****************************/ /********************************************************************/ +float OutiersRemoval2D(unsigned char *MASK, unsigned char *MASK_upd, long i, long j, long dimX, long dimY) +{ + /*if the ROI pixel does not belong to the surrondings, turn it into the surronding*/ + long i_m, j_m, i1, j1, counter; + counter = 0; + for(i_m=-1; i_m<=1; i_m++) { + for(j_m=-1; j_m<=1; j_m++) { + i1 = i+i_m; + j1 = j+j_m; + if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) { + if (MASK[j*dimX+i] != MASK[j1*dimX+i1]) counter++; + } + }} + if (counter >= 8) MASK_upd[j*dimX+i] = MASK[j1*dimX+i1]; + return *MASK_upd; +} + +float Mask_update2D(unsigned char *MASK_temp, unsigned char *MASK_upd, long i, long j, int DiffusWindow, long dimX, long dimY) +{ + long i_m, j_m, i1, j1, CounterOtherClass; + + /* STEP2: in a larger neighbourhood check that the other class is present */ + CounterOtherClass = 0; + for(i_m=-DiffusWindow; i_m<=DiffusWindow; i_m++) { + for(j_m=-DiffusWindow; j_m<=DiffusWindow; j_m++) { + i1 = i+i_m; + j1 = j+j_m; + if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) { + if (MASK_temp[j*dimX+i] != MASK_temp[j1*dimX+i1]) CounterOtherClass++; + } + }} + if (CounterOtherClass > 0) { + /* the other class is present in the vicinity of DiffusWindow, continue to STEP 3 */ + /* + STEP 3: Loop through all neighbours in DiffusWindow and check the spatial connection. + Meaning that we're instrested if there are any classes between points A and B that + does not belong to A and B (A,B \in C) + */ + for(i_m=-DiffusWindow; i_m<=DiffusWindow; i_m++) { + for(j_m=-DiffusWindow; j_m<=DiffusWindow; j_m++) { + i1 = i+i_m; + j1 = j+j_m; + if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) { + if (MASK_temp[j*dimX+i] == MASK_temp[j1*dimX+i1]) { + /*A and B points belong to the same class, do STEP 4*/ + /* STEP 4: Run Bresenham line algorithm between A and B points + and convert all points on the way to the class of A/B. + */ + bresenham2D(i, j, i1, j1, MASK_temp, MASK_upd, (long)(dimX), (long)(dimY)); + } + } + }} + } + return *MASK_upd; +} + /* MASKED-constrained 2D linear diffusion (PDE heat equation) */ float LinearDiff_MASK2D(float *Input, unsigned char *MASK, float *Output, float *Eucl_Vec, int DiffusWindow, float lambdaPar, float tau, long dimX, long dimY) { @@ -142,18 +229,18 @@ float diffVal; #pragma omp parallel for shared(Input) private(index,i,j,i1,j1,i_m,j_m,counter,diffVal,indexneighb,class_c,class_n) for(i=0; i= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) { + if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) { indexneighb = j1*dimX+i1; /* neighbour pixel index */ class_c = MASK[index]; /* current class value */ class_n = MASK[indexneighb]; /* neighbour class value */ - + /* perform diffusion only within the same class (given by MASK) */ if (class_n == class_c) diffVal += Output[indexneighb] - Output[index]; } @@ -170,21 +257,21 @@ float NonLinearDiff_MASK2D(float *Input, unsigned char *MASK, float *Output, flo long i,j,i1,j1,i_m,j_m,index,indexneighb,counter; unsigned char class_c, class_n; float diffVal, funcVal; - + #pragma omp parallel for shared(Input) private(index,i,j,i1,j1,i_m,j_m,counter,diffVal,funcVal,indexneighb,class_c,class_n) for(i=0; i= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) { + if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) { indexneighb = j1*dimX+i1; /* neighbour pixel index */ class_c = MASK[index]; /* current class value */ class_n = MASK[indexneighb]; /* neighbour class value */ - + /* perform diffusion only within the same class (given by MASK) */ if (class_n == class_c) { diffVal = Output[indexneighb] - Output[index]; @@ -200,7 +287,7 @@ float NonLinearDiff_MASK2D(float *Input, unsigned char *MASK, float *Output, flo if (fabs(diffVal) <= sigmaPar) funcVal += diffVal*powf((1.0f - powf((diffVal/sigmaPar),2)), 2); } else { printf("%s \n", "No penalty function selected! Use Huber,2 or 3."); - break; } + break; } } } counter++; @@ -212,3 +299,90 @@ float NonLinearDiff_MASK2D(float *Input, unsigned char *MASK, float *Output, flo /********************************************************************/ /***************************3D Functions*****************************/ /********************************************************************/ + + +int bresenham2D(int i, int j, int i1, int j1, unsigned char *MASK, unsigned char *MASK_upd, long dimX, long dimY) { + + int x[] = {i, i1}; + int y[] = {j, j1}; + int steep = (fabs(y[1]-y[0]) > fabs(x[1]-x[0])); + int ystep = 0; + + //printf("[%i][%i][%i][%i]\n", x[1], y[1], steep, kk) ; + + //if (steep == 1) {swap(x[0],y[0]); swap(x[1],y[1]);} + + if (steep == 1) { + + ///swaping + int a, b; + + a = x[0]; + b = y[0]; + x[0] = b; + y[0] = a; + + a = x[1]; + b = y[1]; + x[1] = b; + y[1] = a; + } + + if (x[0] > x[1]) { + int a, b; + a = x[0]; + b = x[1]; + x[0] = b; + x[1] = a; + + a = y[0]; + b = y[1]; + y[0] = b; + y[1] = a; + } //(x[0] > x[1]) + + int delx = x[1]-x[0]; + int dely = fabs(y[1]-y[0]); + int error = 0; + int x_n = x[0]; + int y_n = y[0]; + if (y[0] < y[1]) {ystep = 1;} + else {ystep = -1;} + + for(int n = 0; n= delx) { + y_n = y_n + ystep; + error = error - delx; + } // (2*error >= delx) + //printf("[%i][%i][%i]\n", X_new[n], Y_new[n], n) ; + } // for(int n = 0; n linear diffusion * 6. Number of iterations, for explicit scheme >= 150 is recommended @@ -54,9 +54,12 @@ limitations under the License. #ifdef __cplusplus extern "C" { #endif -CCPI_EXPORT float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, float *Output, float *infovector, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ); +CCPI_EXPORT float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *MASK_upd, float *Output, float *infovector, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ); CCPI_EXPORT float LinearDiff_MASK2D(float *Input, unsigned char *MASK, float *Output, float *Eucl_Vec, int DiffusWindow, float lambdaPar, float tau, long dimX, long dimY); CCPI_EXPORT float NonLinearDiff_MASK2D(float *Input, unsigned char *MASK, float *Output, float *Eucl_Vec, int DiffusWindow, float lambdaPar, float sigmaPar, float tau, int penaltytype, long dimX, long dimY); +CCPI_EXPORT float OutiersRemoval2D(unsigned char *MASK, unsigned char *MASK_upd, long i, long j, long dimX, long dimY); +CCPI_EXPORT float Mask_update2D(unsigned char *MASK_temp, unsigned char *MASK_upd, long i, long j, int DiffusWindow, long dimX, long dimY); +CCPI_EXPORT int bresenham2D(int i, int j, int i1, int j1, unsigned char *MASK, unsigned char *MASK_upd, long dimX, long dimY); #ifdef __cplusplus } #endif diff --git a/src/Python/ccpi/filters/regularisers.py b/src/Python/ccpi/filters/regularisers.py index 1e427bf..00afcc2 100644 --- a/src/Python/ccpi/filters/regularisers.py +++ b/src/Python/ccpi/filters/regularisers.py @@ -127,10 +127,11 @@ def NDF(inputData, regularisation_parameter, edge_parameter, iterations, raise ValueError ('GPU is not available') raise ValueError('Unknown device {0}. Expecting gpu or cpu'\ .format(device)) -def NDF_MASK(inputData, diffuswindow, regularisation_parameter, edge_parameter, iterations, +def NDF_MASK(inputData, maskdata, diffuswindow, regularisation_parameter, edge_parameter, iterations, time_marching_parameter, penalty_type, tolerance_param, device='cpu'): if device == 'cpu': return NDF_MASK_CPU(inputData, + maskdata, diffuswindow, regularisation_parameter, edge_parameter, @@ -140,6 +141,7 @@ def NDF_MASK(inputData, diffuswindow, regularisation_parameter, edge_parameter, tolerance_param) elif device == 'gpu' and gpu_enabled: return NDF_MASK_CPU(inputData, + maskdata, diffuswindow, regularisation_parameter, edge_parameter, diff --git a/src/Python/src/cpu_regularisers.pyx b/src/Python/src/cpu_regularisers.pyx index 305ee1f..ca402c1 100644 --- a/src/Python/src/cpu_regularisers.pyx +++ b/src/Python/src/cpu_regularisers.pyx @@ -24,7 +24,7 @@ cdef extern float SB_TV_CPU_main(float *Input, float *Output, float *infovector, cdef extern float LLT_ROF_CPU_main(float *Input, float *Output, float *infovector, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ); cdef extern float TGV_main(float *Input, float *Output, float *infovector, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, float epsil, int dimX, int dimY, int dimZ); cdef extern float Diffusion_CPU_main(float *Input, float *Output, float *infovector, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ); -cdef extern float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, float *Output, float *infovector, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ); +cdef extern float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *MASK_upd, float *Output, float *infovector, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ); cdef extern float Diffus4th_CPU_main(float *Input, float *Output, float *infovector, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ); cdef extern float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float *infovector, float lambdaPar, int iterationsNumb, float epsil, float eta, int methodTV, int nonneg, 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); @@ -403,18 +403,22 @@ def NDF_MASK_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] + + cdef np.ndarray[np.uint8_t, ndim=2, mode="c"] mask_upd = \ + np.zeros([dims[0],dims[1]], dtype='uint8') cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \ np.zeros([dims[0],dims[1]], dtype='float32') cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \ np.zeros([2], dtype='float32') + # Run constrained nonlinear diffusion iterations for 2D data - DiffusionMASK_CPU_main(&inputData[0,0], &maskData[0,0], &outputData[0,0], &infovec[0], + DiffusionMASK_CPU_main(&inputData[0,0], &maskData[0,0], &mask_upd[0,0], &outputData[0,0], &infovec[0], diffuswindow, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, tolerance_param, dims[1], dims[0], 1) - return (outputData,infovec) + return (mask_upd,outputData,infovec) #****************************************************************# #*************Anisotropic Fourth-Order diffusion*****************# -- cgit v1.2.3