diff options
| -rw-r--r-- | src/Core/regularisers_CPU/DiffusionMASK_core.c | 106 | ||||
| -rw-r--r-- | src/Core/regularisers_CPU/DiffusionMASK_core.h | 2 | ||||
| -rw-r--r-- | src/Python/ccpi/filters/regularisers.py | 6 | ||||
| -rw-r--r-- | src/Python/src/cpu_regularisers.pyx | 15 |
4 files changed, 86 insertions, 43 deletions
diff --git a/src/Core/regularisers_CPU/DiffusionMASK_core.c b/src/Core/regularisers_CPU/DiffusionMASK_core.c index 2af3cb6..9d128bc 100644 --- a/src/Core/regularisers_CPU/DiffusionMASK_core.c +++ b/src/Core/regularisers_CPU/DiffusionMASK_core.c @@ -53,7 +53,14 @@ 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, 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) +void swapVAL(unsigned char *xp, unsigned char *yp) +{ + unsigned char temp = *xp; + *xp = *yp; + *yp = temp; +} + +float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *MASK_upd, unsigned char *SelClassesList, int SelClassesList_length, float *Output, float *infovector, int classesNumb, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ) { long i,j,k; int counterG; @@ -64,9 +71,43 @@ float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *M float re, re1; re = 0.0f; re1 = 0.0f; int count = 0; - unsigned char *MASK_temp; + unsigned char *MASK_temp, *ClassesList, CurrClass, temp; DimTotal = (long)(dimX*dimY*dimZ); + /* defines the list for all classes in the mask */ + ClassesList = (unsigned char*) calloc (classesNumb,sizeof(unsigned char)); + /* indentify all classes in the MASK */ + /* find the smaller class value first */ + /* + unsigned char smallestClass; + smallestClass = MASK[0]; + for(i=0; i<DimTotal; i++) { + if (MASK[i] < smallestClass) smallestClass = MASK[i]; + } + */ + MASK_temp = (unsigned char*) calloc (DimTotal,sizeof(unsigned char)); + copyIm_unchar(MASK, MASK_temp, (long)(dimX), (long)(dimY), (long)(dimZ)); + + int x, y; + for(x=0; x<classesNumb; x++) { + for(y=0; y<classesNumb-1; y++) { + if(MASK_temp[y] > MASK_temp[y+1]) { + temp = MASK_temp[y+1]; + MASK_temp[y+1] = MASK_temp[y]; + MASK_temp[y] = temp; + }}} + + //printf("[%u]\n", MASK_temp[0]); + + CurrClass = MASK_temp[0]; ClassesList[0]= MASK_temp[0]; counterG = 0; + for(i=0; i<DimTotal; i++) { + if (MASK_temp[i] > CurrClass) { + CurrClass = MASK_temp[i]; + ClassesList[counterG] = MASK_temp[i]; + printf("[%u]\n", ClassesList[counterG]); + counterG++; } + } + if (dimZ == 1) { DiffusWindow_tot = (2*DiffusWindow + 1)*(2*DiffusWindow + 1); /* generate a 2D Gaussian kernel for NLM procedure */ @@ -92,7 +133,7 @@ float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *M if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float)); - MASK_temp = (unsigned char*) calloc (DimTotal,sizeof(unsigned char)); + /* copy input into output */ copyIm(Input, Output, (long)(dimX), (long)(dimY), (long)(dimZ)); @@ -110,25 +151,27 @@ float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *M /* copy the updated MASK (clean of outliers) */ copyIm_unchar(MASK_upd, MASK_temp, (long)(dimX), (long)(dimY), (long)(dimZ)); -// #pragma omp parallel for shared(MASK_temp,MASK_upd) private(i,j) -// for(i=0; i<dimX; i++) { -// for(j=0; j<dimY; j++) { -// if (MASK_temp[j*dimX+i] == MASK[j*dimX+i]) { - /*the class of the central pixel has not changed, i.e. the central pixel is not an outlier -> continue */ - i = 162; j = 258; + #pragma omp parallel for shared(MASK_temp,MASK_upd) private(i,j) + for(i=0; i<dimX; i++) { + for(j=0; j<dimY; j++) { + if (MASK_temp[j*dimX+i] == MASK[j*dimX+i]) { + /* !One needs to work with a specific class to avoid overlaps! + hence it is crucial to establish relevant classes */ + if (MASK_temp[j*dimX+i] == 149) { + /* The class of the central pixel has not changed, i.e. the central pixel is not an outlier -> continue */ + /* i = 258; j = 165; */ Mask_update2D(MASK_temp, MASK_upd, i, j, DiffusWindow, (long)(dimX), (long)(dimY)); -// } - //}} + }} + }} } - /* start diffusivity iterations */ + /* The mask has been processed, 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_upd, Output, Eucl_Vec, DiffusWindow, lambdaPar, tau, (long)(dimX), (long)(dimY)); /* constrained linear diffusion */ + else NonLinearDiff_MASK2D(Input, MASK_upd, Output, Eucl_Vec, DiffusWindow, lambdaPar, sigmaPar2, tau, penaltytype, (long)(dimX), (long)(dimY)); /* constrained nonlinear diffusion */ } else { /* running 3D diffusion iterations */ @@ -207,11 +250,10 @@ float Mask_update2D(unsigned char *MASK_temp, unsigned char *MASK_upd, long i, l 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*/ + /* 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)); + 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)); } } }} @@ -301,20 +343,20 @@ float NonLinearDiff_MASK2D(float *Input, unsigned char *MASK, float *Output, flo /********************************************************************/ -int bresenham2D(int i, int j, int i1, int j1, unsigned char *MASK, unsigned char *MASK_upd, long dimX, long dimY) { - +int bresenham2D(int i, int j, int i1, int j1, unsigned char *MASK, unsigned char *MASK_upd, long dimX, long dimY) +{ + int n; 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 + // swaping int a, b; a = x[0]; @@ -349,18 +391,15 @@ int bresenham2D(int i, int j, int i1, int j1, unsigned char *MASK, unsigned char if (y[0] < y[1]) {ystep = 1;} else {ystep = -1;} - for(int n = 0; n<delx+1; n++) { + for(n = 0; n<delx+1; n++) { if (steep == 1) { /* X_new[n] = x_n; Y_new[n] = y_n; */ /*printf("[%i][%i][%u]\n", x_n, y_n, MASK[y_n*dimX+x_n]);*/ - - MASK_upd[y_n*dimX+x_n] = 10; - //if (MASK[j*dimX+i] != MASK[y_n*dimX+x_n]) { - // MASK_upd[y_n*dimX+x_n] = MASK[j*dimX+i]; - //} + // MASK_upd[x_n*dimX+y_n] = 10; + if (MASK[j*dimX+i] != MASK[x_n*dimX+y_n]) MASK_upd[x_n*dimX+y_n] = MASK[j*dimX+i]; } else { /* @@ -368,11 +407,8 @@ int bresenham2D(int i, int j, int i1, int j1, unsigned char *MASK, unsigned char Y_new[n] = x_n; */ // printf("[%i][%i][%u]\n", y_n, x_n, MASK[x_n*dimX+y_n]); - // MASK_upd[x_n*dimX+y_n] = 1; - //if (MASK[j*dimX+i] != MASK[x_n*dimX+y_n]) { - // MASK_upd[x_n*dimX+y_n] = MASK[j*dimX+i]; - //} - MASK_upd[x_n*dimX+y_n] = 20; + // MASK_upd[y_n*dimX+x_n] = 20; + if (MASK[j*dimX+i] != MASK[y_n*dimX+x_n]) MASK_upd[y_n*dimX+x_n] = MASK[j*dimX+i]; } x_n = x_n + 1; error = error + dely; diff --git a/src/Core/regularisers_CPU/DiffusionMASK_core.h b/src/Core/regularisers_CPU/DiffusionMASK_core.h index 48142af..a032905 100644 --- a/src/Core/regularisers_CPU/DiffusionMASK_core.h +++ b/src/Core/regularisers_CPU/DiffusionMASK_core.h @@ -54,7 +54,7 @@ limitations under the License. #ifdef __cplusplus extern "C" { #endif -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 DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *MASK_upd, unsigned char *SelClassesList, int SelClassesList_length, float *Output, float *infovector, int classesNumb, 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); diff --git a/src/Python/ccpi/filters/regularisers.py b/src/Python/ccpi/filters/regularisers.py index 00afcc2..610907d 100644 --- a/src/Python/ccpi/filters/regularisers.py +++ b/src/Python/ccpi/filters/regularisers.py @@ -127,11 +127,13 @@ 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, maskdata, diffuswindow, regularisation_parameter, edge_parameter, iterations, +def NDF_MASK(inputData, maskdata, select_classes, total_classesNum, diffuswindow, regularisation_parameter, edge_parameter, iterations, time_marching_parameter, penalty_type, tolerance_param, device='cpu'): if device == 'cpu': return NDF_MASK_CPU(inputData, maskdata, + select_classes, + total_classesNum, diffuswindow, regularisation_parameter, edge_parameter, @@ -142,6 +144,8 @@ def NDF_MASK(inputData, maskdata, diffuswindow, regularisation_parameter, edge_p elif device == 'gpu' and gpu_enabled: return NDF_MASK_CPU(inputData, maskdata, + select_classes, + total_classesNum, diffuswindow, regularisation_parameter, edge_parameter, diff --git a/src/Python/src/cpu_regularisers.pyx b/src/Python/src/cpu_regularisers.pyx index ca402c1..78c46e2 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, 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 DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *MASK_upd, unsigned char *SelClassesList, int SelClassesList_length, float *Output, float *infovector, int classesNumb, 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); @@ -384,14 +384,16 @@ def NDF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, #****************************************************************# #********Constrained Nonlinear(Isotropic) Diffusion**************# #****************************************************************# -def NDF_MASK_CPU(inputData, maskData, diffuswindow, regularisation_parameter, edge_parameter, iterationsNumb,time_marching_parameter, penalty_type, tolerance_param): +def NDF_MASK_CPU(inputData, maskData, select_classes, total_classesNum, diffuswindow, regularisation_parameter, edge_parameter, iterationsNumb,time_marching_parameter, penalty_type, tolerance_param): if inputData.ndim == 2: - return NDF_MASK_2D(inputData, maskData, diffuswindow, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, tolerance_param) + return NDF_MASK_2D(inputData, maskData, select_classes, total_classesNum, diffuswindow, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, tolerance_param) elif inputData.ndim == 3: return 0 def NDF_MASK_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, np.ndarray[np.uint8_t, ndim=2, mode="c"] maskData, + np.ndarray[np.uint8_t, ndim=1, mode="c"] select_classes, + int total_classesNum, int diffuswindow, float regularisation_parameter, float edge_parameter, @@ -403,7 +405,8 @@ def NDF_MASK_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] - + select_classes_length = select_classes.shape[0] + 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 = \ @@ -413,8 +416,8 @@ def NDF_MASK_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, # Run constrained nonlinear diffusion iterations for 2D data - DiffusionMASK_CPU_main(&inputData[0,0], &maskData[0,0], &mask_upd[0,0], &outputData[0,0], &infovec[0], - diffuswindow, regularisation_parameter, edge_parameter, iterationsNumb, + DiffusionMASK_CPU_main(&inputData[0,0], &maskData[0,0], &mask_upd[0,0], &select_classes[0], select_classes_length, &outputData[0,0], &infovec[0], + total_classesNum, diffuswindow, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, tolerance_param, dims[1], dims[0], 1) |