progress with bresenham

4 files changed, 207 insertions, 24 deletions

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*****************#