summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rwxr-xr-xbuild/run.sh2
-rw-r--r--src/Core/regularisers_CPU/DiffusionMASK_core.c210
-rw-r--r--src/Core/regularisers_CPU/DiffusionMASK_core.h7
-rw-r--r--src/Python/ccpi/filters/regularisers.py4
-rw-r--r--src/Python/src/cpu_regularisers.pyx10
5 files changed, 208 insertions, 25 deletions
diff --git a/build/run.sh b/build/run.sh
index b40d222..0a153aa 100755
--- a/build/run.sh
+++ b/build/run.sh
@@ -21,7 +21,7 @@ make install
cp install/lib/libcilreg.so install/python/ccpi/filters
cd install/python
export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:../lib
-# spyder
+spyder
############### Matlab(linux)###############
# export LD_PRELOAD=/home/algol/anaconda3/lib/libstdc++.so.6 # if there is libstdc error in matlab
# PATH="/path/to/mex/:$PATH" LD_LIBRARY_PATH="/path/to/library:$LD_LIBRARY_PATH" matlab
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<dimX; i++) {
+ for(j=0; j<dimY; j++) {
+ /* STEP1: in a smaller neighbourhood check that the current pixel is NOT an outlier */
+ OutiersRemoval2D(MASK, MASK_upd, i, j, (long)(dimX), (long)(dimY));
+ }}
+ /* 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;
+ 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<dimX; i++) {
- for(j=0; j<dimY; j++) {
+ for(j=0; j<dimY; j++) {
index = j*dimX+i; /* current pixel index */
counter = 0; diffVal = 0.0f;
for(i_m=-DiffusWindow; i_m<=DiffusWindow; i_m++) {
- for(j_m=-DiffusWindow; j_m<=DiffusWindow; j_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 (((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<dimX; i++) {
- for(j=0; j<dimY; j++) {
+ for(j=0; j<dimY; j++) {
index = j*dimX+i; /* current pixel index */
counter = 0; diffVal = 0.0f; funcVal = 0.0f;
for(i_m=-DiffusWindow; i_m<=DiffusWindow; i_m++) {
- for(j_m=-DiffusWindow; j_m<=DiffusWindow; j_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 (((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+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];
+ //}
+ }
+ else {
+ /*
+ X_new[n] = y_n;
+ 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;
+ }
+ x_n = x_n + 1;
+ error = error + dely;
+
+ if (2*error >= 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<delx+1; n++)
+
+ return 0;
+}
diff --git a/src/Core/regularisers_CPU/DiffusionMASK_core.h b/src/Core/regularisers_CPU/DiffusionMASK_core.h
index 8890c73..48142af 100644
--- a/src/Core/regularisers_CPU/DiffusionMASK_core.h
+++ b/src/Core/regularisers_CPU/DiffusionMASK_core.h
@@ -33,7 +33,7 @@ limitations under the License.
* Input Parameters:
* 1. Noisy image/volume
* 2. MASK (in unsigned short 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
@@ -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*****************#
generated by cgit v1.2.3 (git 2.39.1) at 2025-11-04 06:18:35 +0000