3D version work

Demos/Python/DemoFISTA_artifacts2D.py

@@ -210,7 +210,7 @@
                     device='gpu')
 
 RecFISTA_Huber_reg = RectoolsIR.FISTA(noisy_zing_stripe, 
-                                   huber_data_threshold=3.5,
+                                   huber_data_threshold=4.0,
                                    iterationsFISTA = 20, 
                                    regularisation = 'ROF_TV', 
                                    regularisation_parameter = 0.01,
@@ -250,7 +250,7 @@
                     device='gpu')
 
 RecFISTA_Huber_RingModel = RectoolsIR.FISTA(noisy_zing_stripe, 
-                                   huber_data_threshold=3.5,
+                                   huber_data_threshold=4.0,
                                    ring_model_horiz_size=9,
                                    iterationsFISTA = 20, 
                                    regularisation = 'ROF_TV', 

Demos/Python/DemoFISTA_artifacts3D.py

@@ -26,8 +26,8 @@
 
 print ("Building 3D phantom using TomoPhantom software")
 tic=timeit.default_timer()
-model = 9 # select a model number from the library
-N_size = 256 # Define phantom dimensions using a scalar value (cubic phantom)
+model = 13 # select a model number from the library
+N_size = 150 # Define phantom dimensions using a scalar value (cubic phantom)
 path = os.path.dirname(tomophantom.__file__)
 path_library3D = os.path.join(path, "Phantom3DLibrary.dat")
 #This will generate a N_size x N_size x N_size phantom (3D)
@@ -63,7 +63,7 @@
 projData3D_analyt= TomoP3D.ModelSino(model, N_size, Horiz_det, Vert_det, angles, path_library3D)
 
 intens_max = 70
-sliceSel = 150
+sliceSel = 100
 plt.figure() 
 plt.subplot(131)
 plt.imshow(projData3D_analyt[:,sliceSel,:],vmin=0, vmax=intens_max)
@@ -105,15 +105,41 @@
 #%%
 # initialise tomobar DIRECT reconstruction class ONCE
 from tomobar.methodsDIR import RecToolsDIR
-RectoolsDIR = RecToolsDIR(DetectorsDimH = Horiz_det,  # DetectorsDimH # detector dimension (horizontal)
+Rectools = RecToolsDIR(DetectorsDimH = Horiz_det,  # DetectorsDimH # detector dimension (horizontal)
                     DetectorsDimV = Vert_det,  # DetectorsDimV # detector dimension (vertical) for 3D case only
-                    CenterRotOffset = None, # Center of Rotation (CoR) scalar (for 3D case only)
+                    CenterRotOffset = 0.0, # Center of Rotation (CoR) scalar (for 3D case only)
                     AnglesVec = angles_rad, # array of angles in radians
                     ObjSize = N_size, # a scalar to define reconstructed object dimensions
                     device = 'gpu')
 #%%
+# testing RING
+Forwproj = Rectools.FORWPROJ(phantom_tm)
+residual = Forwproj - projData3D_norm
+#residual = np.swapaxes(residual,0,1)
+#%%
+from tomobar.supp.addmodules import RING_WEIGHTS
+ring_model_horiz_size = 13
+ring_model_vert_size=0
+ring_model_slices_size=0
+
+offset_rings = RING_WEIGHTS(residual, ring_model_horiz_size, ring_model_vert_size, ring_model_slices_size)
+
+intens_max = 70
+sliceSel = 120
+plt.figure() 
+plt.subplot(131)
+plt.imshow(offset_rings[:,sliceSel,:])
+plt.title('2D Projection (erroneous)')
+plt.subplot(132)
+plt.imshow(offset_rings[sliceSel,:,:])
+plt.title('Sinogram view')
+plt.subplot(133)
+plt.imshow(offset_rings[:,:,sliceSel])
+plt.title('Tangentogram view')
+plt.show()
+#%%
 print ("Reconstruction using FBP from tomobar")
-recNumerical= RectoolsDIR.FBP(projData3D_norm) # FBP reconstruction
+recNumerical= Rectools.FBP(projData3D_norm) # FBP reconstruction
 
 sliceSel = int(0.5*N_size)
 max_val = 1
@@ -150,7 +176,7 @@
                     datafidelity='LS',# data fidelity, choose LS, PWLS (wip), GH (wip), Student (wip)
                     nonnegativity='ENABLE', # enable nonnegativity constraint (set to 'ENABLE')
                     OS_number = 10, # the number of subsets, NONE/(or > 1) ~ classical / ordered subsets
-                    tolerance = 1e-08, # tolerance to stop outer iterations earlier
+                    tolerance = 1e-10, # tolerance to stop outer iterations earlier
                     device='gpu')
 lc = RectoolsIR.powermethod() # calculate Lipschitz constant
 #%%
@@ -165,6 +191,11 @@
                                 regularisation_iterations = 300,
                                 lipschitz_const = lc)
 
+
+Qtools = QualityTools(phantom_tm, RecFISTA_reg)
+RMSE_FISTA_TV = Qtools.rmse()
+print("RMSE for FISTA-OS-TV is {}".format(RMSE_FISTA_TV))
+
 sliceSel = int(0.5*N_size)
 max_val = 1
 plt.figure() 
@@ -190,6 +221,10 @@
                                 regularisation_iterations = 300,
                                 lipschitz_const = lc)
 
+Qtools = QualityTools(phantom_tm, RecFISTA_Huber_TV)
+RMSE_FISTA_HUBER_TV = Qtools.rmse()
+print("RMSE for FISTA-OS-Huber-TV is {}".format(RMSE_FISTA_HUBER_TV))
+
 sliceSel = int(0.5*N_size)
 max_val = 1
 plt.figure() 
@@ -206,18 +241,20 @@
 plt.title('3D Huber Rec, sagittal')
 plt.show()
 #%%
-#%%
 # Run FISTA reconstrucion algorithm with 3D regularisation
 RecFISTA_HuberRING_TV = RectoolsIR.FISTA(projData3D_norm, 
                                 iterationsFISTA = 20, 
                                 huber_data_threshold = 1.0,
                                 ring_model_horiz_size= 7,
-                                ring_model_vert_size = 1,
                                 regularisation = 'FGP_TV', 
                                 regularisation_parameter = 0.00015, 
                                 regularisation_iterations = 300,
                                 lipschitz_const = lc)
 
+Qtools = QualityTools(phantom_tm, RecFISTA_HuberRING_TV)
+RMSE_FISTA_HUBER_RING_TV = Qtools.rmse()
+print("RMSE for FISTA-OS-Huber-Ring-TV is {}".format(RMSE_FISTA_HUBER_RING_TV))
+
 sliceSel = int(0.5*N_size)
 max_val = 1
 plt.figure() 

src/Core/functions_CPU/RingWeights_core.c

@@ -46,54 +46,59 @@ void swap(float *xp, float *yp)
 }
 
 
-float RingWeights_main(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, long anglesDim, long detectorsDim, long slices)
+float RingWeights_main(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, int slices_window_halfsize, long anglesDim, long detectorsDim, long slices)
 {
     long i, j, k;
     int full_window;
+    full_window = (int)(2*slices_window_halfsize+1)*(2*vert_window_halfsize+1)*(2*horiz_window_halfsize+1);
 
     if (slices == 1) {
     /****************2D INPUT ***************/
-    full_window = (2*horiz_window_halfsize+1);
     #pragma omp parallel for shared(residual, weights) private(j, i)
     for(j=0; j<detectorsDim; j++) {
         for(i=0; i<anglesDim; i++) {
-          RingWeights2D(residual, weights, horiz_window_halfsize, full_window, anglesDim, detectorsDim, i, j);
+          RingWeights2D(residual, weights, horiz_window_halfsize, vert_window_halfsize, full_window, anglesDim, detectorsDim, i, j);
         }}
     }
     else {
     /****************3D INPUT ***************/
-    full_window = (2*horiz_window_halfsize+1)*(2*vert_window_halfsize+1);
+    printf("%i %i %i \n", slices, anglesDim, detectorsDim);
     #pragma omp parallel for shared(residual, weights) private(k, j, i)
-    for(k=0; k<slices; k++) {
-      for(j=0; j<detectorsDim; j++) {
-        for(i=0; i<anglesDim; i++) {
-          RingWeights3D(residual, weights, horiz_window_halfsize, vert_window_halfsize, full_window, anglesDim, detectorsDim, slices, i, j, k);
+    for(i = 0; i<anglesDim; i++) {
+      for(j = 0; j<detectorsDim; j++) {
+        for(k = 0; k<slices; k++) {
+          RingWeights3D(residual, weights, horiz_window_halfsize, vert_window_halfsize, slices_window_halfsize, full_window, anglesDim, detectorsDim, slices, i, j, k);
         }}}
     }
   return *weights;
 }
 /********************************************************************/
 /***************************2D Functions*****************************/
 /********************************************************************/
-float RingWeights2D(float *residual, float *weights, int horiz_window_halfsize, int full_window, long anglesDim, long detectorsDim, long i, long j)
+float RingWeights2D(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, int full_window, long anglesDim, long detectorsDim, long i, long j)
 {
             float *Values_Vec;
             long k, j1, index;
-            int counter, x, y;
+            int counter, x, y, midval;
 
             index = i*detectorsDim+j;
             Values_Vec = (float*) calloc (full_window, sizeof(float));
+            midval = (int)(0.5f*full_window) - 1;
 
             /* intiate the estimation of the backround using strictly horizontal values */
             counter = 0;
+            long l, i1;
             for(k=-horiz_window_halfsize; k < horiz_window_halfsize; k++) {
                 j1 = j + k;
-                if ((j1 >= 0) && (j1 < detectorsDim)) {
-                  Values_Vec[counter] = residual[i*detectorsDim+j1];
+                for(l=-vert_window_halfsize; l < vert_window_halfsize; l++) {
+                  i1 = i + l;
+                if ((j1 >= 0) && (j1 < detectorsDim) && (i1 >= 0) && (i1 < anglesDim)) {
+                  Values_Vec[counter] = residual[i1*detectorsDim+j1];
                 }
                 else Values_Vec[counter] = residual[index];
                 counter ++;
-            }
+            }}
+
             /* perform sorting of the vector array */
             for (x = 0; x < counter-1; x++)  {
                 for (y = 0; y < counter-x-1; y++)  {
@@ -102,59 +107,39 @@ float RingWeights2D(float *residual, float *weights, int horiz_window_halfsize,
                     }
                 }
             }
-            /* include diagonal values in the estimation of the backround */
-            int i1;
-            float *Values_Vec_diag1;
-            Values_Vec_diag1 = (float*) calloc (full_window, sizeof(float));
-            counter = 0;
-            for(k=-horiz_window_halfsize; k < horiz_window_halfsize; k++) {
-                j1 = j + k;
-                i1 = i + k;
-                if ((j1 >= 0) && (j1 < detectorsDim) && (i1 >= 0) && (i1 < anglesDim)) {
-                  Values_Vec_diag1[counter] = residual[i1*detectorsDim+j1];
-                }
-                else Values_Vec_diag1[counter] = residual[index];
-                counter ++;
-            }
-            /* perform sorting of the vector array */
-            for (x = 0; x < counter-1; x++)  {
-                for (y = 0; y < counter-x-1; y++)  {
-                    if (Values_Vec_diag1[y] > Values_Vec_diag1[y+1]) {
-                        swap(&Values_Vec_diag1[y], &Values_Vec_diag1[y+1]);
-                    }
-                }
-            }
 
-            weights[index] = residual[index] - 0.5f*(Values_Vec[horiz_window_halfsize] + Values_Vec_diag1[horiz_window_halfsize]);
+            weights[index] = residual[index] - Values_Vec[midval];
 
       free(Values_Vec);
-      free(Values_Vec_diag1);
       return *weights;
 }
 
-float RingWeights3D(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, int full_window, long anglesDim, long detectorsDim, long slices, long i, long j, long k)
+float RingWeights3D(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, int slices_window_halfsize, int full_window, long anglesDim, long detectorsDim, long slices, long i, long j, long k)
 {
   float *Values_Vec;
-  long k1, j1, l, v, index;
+  long k1, j1, i1, l, v, m, index;
   int counter, x, y, midval;
 
   midval = (int)(0.5f*full_window) - 1;
   index = (detectorsDim*anglesDim*k) + i*detectorsDim+j;
   Values_Vec = (float*) calloc (full_window, sizeof(float));
-
-  /* intiate the estimation of the backround using strictly horizontal values */
+/*
   counter = 0;
   for(l=-horiz_window_halfsize; l < horiz_window_halfsize; l++) {
-      for(v=-vert_window_halfsize; v < vert_window_halfsize; v++) {
-      j1 = j + l;
-      k1 = k + v;
-      if ((j1 >= 0) && (j1 < detectorsDim) && (k1 >= 0) && (k1 < slices)) {
-        Values_Vec[counter] = residual[(detectorsDim*anglesDim*k1) + i*detectorsDim+j1];
+    j1 = j + l;
+    for(m=-vert_window_halfsize; m < vert_window_halfsize; m++) {
+      i1 = i + m;
+      for(v=-slices_window_halfsize; v < slices_window_halfsize; v++) {
+        k1 = k + v;
+      if ((j1 >= 0) && (j1 < detectorsDim) && (k1 >= 0) && (k1 < slices) &&  (i1 >= 0) && (i1 < anglesDim)) {
+        Values_Vec[counter] = residual[(detectorsDim*anglesDim*k1) + i1*detectorsDim+j1];
       }
       else Values_Vec[counter] = residual[index];
       counter ++;
-  }}
+  }}}
+*/
   /* perform sorting of the vector array */
+  /*
   for (x = 0; x < counter-1; x++)  {
       for (y = 0; y < counter-x-1; y++)  {
           if (Values_Vec[y] > Values_Vec[y+1]) {
@@ -163,6 +148,9 @@ float RingWeights3D(float *residual, float *weights, int horiz_window_halfsize,
       }
   }
   weights[index] = residual[index] - Values_Vec[midval];
+  */
+  j1 = j + 15;
+  if ((j1 > 0) && (j1 < detectorsDim)) weights[index] = residual[(detectorsDim*anglesDim*k) + i*detectorsDim+j1];
 
   free(Values_Vec);
   return *weights;

src/Core/functions_CPU/RingWeights_core.h

@@ -39,11 +39,11 @@
 #ifdef __cplusplus
 extern "C" {
 #endif
-float RingWeights_main(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, long anglesDim, long detectorsDim, long slices);
+float RingWeights_main(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, int slices_window_halfsize, long anglesDim, long detectorsDim, long slices);
 /************2D functions ***********/
-float RingWeights2D(float *residual, float *weights, int horiz_window_halfsize, int full_window, long anglesDim, long detectorsDim, long i, long j);
+float RingWeights2D(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, int full_window, long anglesDim, long detectorsDim, long i, long j);
 /************3D functions ***********/
-float RingWeights3D(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, int full_window, long anglesDim, long detectorsDim, long slices, long i, long j, long k);
+float RingWeights3D(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, int slices_window_halfsize, int full_window, long anglesDim, long detectorsDim, long slices, long i, long j, long k);
 #ifdef __cplusplus
 }
 #endif

src/Python/src/cpu_wrappers.pyx

@@ -15,17 +15,18 @@ import cython
 import numpy as np
 cimport numpy as np
 
-cdef extern float RingWeights_main(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, long anglesDim, long detectorsDim, long slices);
+cdef extern float RingWeights_main(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, int slices_window_halfsize, long anglesDim, long detectorsDim, long slices);
 
 ##############################################################################
-def RING_WEIGHTS(residual, horiz_window_halfsize, vert_window_halfsize):
+def RING_WEIGHTS(residual, horiz_window_halfsize, vert_window_halfsize, slices_window_halfsize):
     if residual.ndim == 2:
-        return RING_WEIGHTS_2D(residual, horiz_window_halfsize)
+        return RING_WEIGHTS_2D(residual, horiz_window_halfsize, vert_window_halfsize)
     elif residual.ndim == 3:
-        return RING_WEIGHTS_3D(residual, horiz_window_halfsize, vert_window_halfsize)
+        return RING_WEIGHTS_3D(residual, horiz_window_halfsize, vert_window_halfsize, slices_window_halfsize)
 
 def RING_WEIGHTS_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] residual,
-                     int horiz_window_halfsize):
+                     int horiz_window_halfsize,
+                     int vert_window_halfsize):
 
     cdef long dims[2]
     dims[0] = residual.shape[0]
@@ -34,13 +35,13 @@ def RING_WEIGHTS_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] residual,
     cdef np.ndarray[np.float32_t, ndim=2, mode="c"] weights = \
             np.zeros([dims[0],dims[1]], dtype='float32')
 
-    RingWeights_main(&residual[0,0], &weights[0,0], horiz_window_halfsize, 0,  dims[1], dims[0], 1);
-
+    RingWeights_main(&residual[0,0], &weights[0,0], horiz_window_halfsize, vert_window_halfsize, 0, dims[1], dims[0], 1);
     return weights
 
 def RING_WEIGHTS_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] residual,
                      int horiz_window_halfsize,
-                     int vert_window_halfsize):
+                     int vert_window_halfsize,
+                     int slices_window_halfsize):
 
     cdef long dims[3]
     dims[0] = residual.shape[0]
@@ -50,6 +51,5 @@ def RING_WEIGHTS_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] residual,
     cdef np.ndarray[np.float32_t, ndim=3, mode="c"] weights = \
             np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
 
-    RingWeights_main(&residual[0,0,0], &weights[0,0,0], horiz_window_halfsize, vert_window_halfsize, dims[2], dims[1], dims[0]);
-
+    RingWeights_main(&residual[0,0,0], &weights[0,0,0], horiz_window_halfsize, vert_window_halfsize, slices_window_halfsize, dims[1], dims[2], dims[0]);
     return weights

src/Python/tomobar/methodsDIR.py

@@ -77,7 +77,7 @@ def FORWPROJ(self, image):
             sinogram = Atools.forwproj(image)
         if (self.geom == '3D'):
             from tomobar.supp.astraOP import AstraTools3D
-            Atools = AstraTools3D(self.DetectorsDimH, self.DetectorsDimV, self.AnglesVec, self.ObjSize) # initiate 3D ASTRA class object
+            Atools = AstraTools3D(self.DetectorsDimH, self.DetectorsDimV, self.AnglesVec, self.CenterRotOffset, self.ObjSize) # initiate 3D ASTRA class object
             sinogram = Atools.forwproj(image)
         return sinogram
     def BACKPROJ(self, sinogram):

src/Python/tomobar/methodsIR.py

@@ -172,7 +172,8 @@ def FISTA(self,
               lambdaR_L1 = 0.0, # regularisation parameter for GH data model
               alpha_ring = 50, # GH data model convergence accelerator (use carefully -> can generate artifacts)
               ring_model_horiz_size = None, # enable a better model to supress ring artifacts, size of the window defines a possible thickness of ring artifacts
-              ring_model_vert_size = 1, # for 3D case only define a vertical window
+              ring_model_vert_size = 0, # define a vertical window
+              ring_model_slices_size = 0, # 3D case only, define a slices window
               huber_data_threshold = 0.0, # threshold parameter for __Huber__ data fidelity 
               student_data_threshold = 0.0, # threshold parameter for __Students't__ data fidelity 
               initialise = None, # initialise reconstruction with an array
@@ -321,7 +322,7 @@ def FISTA(self,
                         multHuber[(np.where(np.abs(res) > huber_data_threshold))] = np.divide(huber_data_threshold, np.abs(res[(np.where(np.abs(res) > huber_data_threshold))]))
                     else:
                         # Apply smarter Huber model to supress ring artifacts
-                        offset_rings = RING_WEIGHTS(res, ring_model_horiz_size, ring_model_vert_size)
+                        offset_rings = RING_WEIGHTS(res, ring_model_horiz_size, ring_model_vert_size, ring_model_slices_size)
                         tempRes = res+offset_rings
                         multHuber[(np.where(np.abs(tempRes) > huber_data_threshold))] = np.divide(huber_data_threshold, np.abs(tempRes[(np.where(np.abs(tempRes) > huber_data_threshold))]))
                     if (self.OS_number > 1):
@@ -337,7 +338,7 @@ def FISTA(self,
                         multStudent = np.divide(2.0, student_data_threshold**2 + res**2)
                     else:
                         # Apply smarter Students't model to supress ring artifacts
-                        offset_rings = RING_WEIGHTS(res, ring_model_horiz_size, ring_model_vert_size)
+                        offset_rings = RING_WEIGHTS(res, ring_model_horiz_size, ring_model_vert_size, ring_model_slices_size)
                         tempRes = res+offset_rings
                         multStudent = np.divide(2.0, student_data_threshold**2 + tempRes**2)
                     if (self.OS_number > 1):