Per-projection geometry cone3D_souv

Wrappers/Python/cil/framework/labels.py

@@ -241,6 +241,7 @@ class AcquisitionType(Flag):
     CONE = auto()
     DIM2 = auto()
     DIM3 = auto()
+    CONE_SOUV = auto()
 
     def validate(self):
         """
@@ -262,7 +263,7 @@ def geometry(self):
         """
         Returns the label for the geometry type
         """
-        return self & (self.PARALLEL | self.CONE)
+        return self & (self.PARALLEL | self.CONE | self.CONE_SOUV)
 
     @classmethod
     def _missing_(cls, value):

Wrappers/Python/cil/plugins/astra/processors/FDK_Flexible.py

@@ -68,7 +68,7 @@ def check_input(self, dataset):
             raise ValueError("Expected input data to be single channel, got {0}"\
                  .format(self.sinogram_geometry.channels))
 
-        if self.sinogram_geometry.geom_type != 'cone':
+        if self.sinogram_geometry.geom_type != 'cone' and self.sinogram_geometry.geom_type != 'cone_souv':
             raise ValueError("Expected input data to be cone beam geometry , got {0}"\
                  .format(self.sinogram_geometry.geom_type))
 

Wrappers/Python/cil/plugins/astra/utilities/convert_geometry_to_astra.py

@@ -39,15 +39,24 @@ def convert_geometry_to_astra(volume_geometry, sinogram_geometry):
         The ASTRA vol_geom and proj_geom
 
     """
+    if sinogram_geometry.geom_type == AcquisitionType.CONE_SOUV:
+        raise ValueError('Cone-SOUV geometry is not supported by this function, use convert_geometry_to_astra_vec_3D instead')
+
     # determine if the geometry is 2D or 3D
     dimension = AcquisitionType.DIM3 if sinogram_geometry.pixel_num_v > 1 else AcquisitionType.DIM2
 
     #get units
 
     if sinogram_geometry.config.angles.angle_unit == AngleUnit.DEGREE:
-        angles_rad = sinogram_geometry.config.angles.angle_data * np.pi / 180.0
+        angles_rad = -sinogram_geometry.config.angles.angle_data * np.pi / 180.0
     else:
-        angles_rad = sinogram_geometry.config.angles.angle_data
+        angles_rad = -sinogram_geometry.config.angles.angle_data
+
+    if 'right' in sinogram_geometry.config.panel.origin:
+        angles_rad += np.pi
+
+    if AcquisitionType.DIM3 & dimension and 'top' in sinogram_geometry.config.panel.origin:
+        raise ValueError('Top origin is not supported for ASTRA 3D simple geometries, either flip the data or use convert_geometry_to_astra_vec_3D')
 
     if AcquisitionType.DIM2 & dimension:
         vol_geom = astra.create_vol_geom(volume_geometry.voxel_num_y,
@@ -61,12 +70,12 @@ def convert_geometry_to_astra(volume_geometry, sinogram_geometry):
             proj_geom = astra.create_proj_geom('parallel',
                                                sinogram_geometry.pixel_size_h,
                                                sinogram_geometry.pixel_num_h,
-                                               -angles_rad)
+                                               angles_rad)
         elif sinogram_geometry.geom_type == 'cone':
             proj_geom = astra.create_proj_geom('fanflat',
                                                sinogram_geometry.pixel_size_h,
                                                sinogram_geometry.pixel_num_h,
-                                               -angles_rad,
+                                               angles_rad,
                                                np.abs(sinogram_geometry.dist_source_center),
                                                np.abs(sinogram_geometry.dist_center_detector))
         else:
@@ -89,14 +98,14 @@ def convert_geometry_to_astra(volume_geometry, sinogram_geometry):
                                                sinogram_geometry.pixel_size_v,
                                                sinogram_geometry.pixel_num_v,
                                                sinogram_geometry.pixel_num_h,
-                                               -angles_rad)
+                                               angles_rad)
         elif sinogram_geometry.geom_type == 'cone':
             proj_geom = astra.create_proj_geom('cone',
                                                sinogram_geometry.pixel_size_h,
                                                sinogram_geometry.pixel_size_v,
                                                sinogram_geometry.pixel_num_v,
                                                sinogram_geometry.pixel_num_h,
-                                               -angles_rad,
+                                               angles_rad,
                                                np.abs(sinogram_geometry.dist_source_center),
                                                np.abs(sinogram_geometry.dist_center_detector))
         else:

Wrappers/Python/cil/plugins/astra/utilities/convert_geometry_to_astra_vec_2D.py

@@ -20,6 +20,7 @@
 import astra
 import numpy as np
 from cil.framework.labels import AngleUnit
+from cil.framework.labels import AcquisitionType, AcquisitionType
 
 def convert_geometry_to_astra_vec_2D(volume_geometry, sinogram_geometry_in):
 
@@ -42,6 +43,9 @@ def convert_geometry_to_astra_vec_2D(volume_geometry, sinogram_geometry_in):
     """
     sinogram_geometry = sinogram_geometry_in.copy()
 
+    if sinogram_geometry.geom_type == AcquisitionType.CONE_SOUV:
+        raise ValueError('Cone-SOUV geometry is not supported by this function, use convert_geometry_to_astra_vec_3D instead')
+
     #this catches behaviour modified after CIL 21.3.1
     try:
         sinogram_geometry.config.system.align_reference_frame('cil')

Wrappers/Python/cil/plugins/astra/utilities/convert_geometry_to_astra_vec_3D.py

@@ -43,21 +43,29 @@ def convert_geometry_to_astra_vec_3D(volume_geometry, sinogram_geometry_in):
 
     sinogram_geometry = sinogram_geometry_in.copy()
 
-    #this catches behaviour modified after CIL 21.3.1
-    try:
-        sinogram_geometry.config.system.align_reference_frame('cil')
-    except:
-        sinogram_geometry.config.system.update_reference_frame()
+    # Only needed when using the traditional geometry set with rotation angles 
+    if sinogram_geometry.geom_type != 'cone_souv':
+
+        #this catches behaviour modified after CIL 21.3.1    
+        try:
+            sinogram_geometry.config.system.align_reference_frame('cil')
+        except:
+            sinogram_geometry.config.system.update_reference_frame()
 
 
-    angles = sinogram_geometry.config.angles
+        angles = sinogram_geometry.config.angles
+
+        #get units
+        degrees = angles.angle_unit == AngleUnit.DEGREE
+
     system = sinogram_geometry.config.system
     panel = sinogram_geometry.config.panel
+
+    # Translation vector that will modify the centre of the reconstructed volume
+    # (by defalut, no translation)
+    translation = [0.0, 0.0, 0.0]
 
-    #get units
-    degrees = angles.angle_unit == AngleUnit.DEGREE
-
-    if AcquisitionType.DIM2 & sinogram_geometry.dimension:
+    if sinogram_geometry.dimension == '2D':
         #create a 3D astra geom from 2D CIL geometry
         volume_geometry_temp = volume_geometry.copy()
         volume_geometry_temp.voxel_num_z = 1
@@ -89,7 +97,9 @@ def convert_geometry_to_astra_vec_3D(volume_geometry, sinogram_geometry_in):
             src[1] = system.source.position[1]
             projector = 'cone_vec'
 
-    else:
+    # Only needed when using the traditional geometry set with rotation angles 
+    elif sinogram_geometry.geom_type != 'cone_souv':
+
         volume_geometry_temp = volume_geometry.copy()
 
         row = panel.pixel_size[0] * system.detector.direction_x.reshape(3,1)
@@ -104,33 +114,65 @@ def convert_geometry_to_astra_vec_3D(volume_geometry, sinogram_geometry_in):
         if sinogram_geometry.geom_type == 'parallel':
             src = system.ray.direction.reshape(3,1)
             projector = 'parallel3d_vec'
-        else:
+        elif sinogram_geometry.geom_type != 'cone_souv':
             src = system.source.position.reshape(3,1)
             projector = 'cone_vec'
+    # Use the per-projection geometry
+    else:
+        volume_geometry_temp = volume_geometry.copy()
+
+        # roi 
+        current_centre = [volume_geometry_temp.center_x, volume_geometry_temp.center_y, volume_geometry_temp.center_z]
+
+        # Compute a translation vector that will modify the centre of the reconstructed volume
+        translation = np.array(system.volume_centre.position)
+
+        projector = 'cone_vec'
 
     #Build for astra 3D only
-    vectors = np.zeros((angles.num_positions, 12))
+    # Use the traditional geometry set with rotation angles
+    if sinogram_geometry.geom_type != 'cone_souv':
+        vectors = np.zeros((angles.num_positions, 12))
+
+        for i, theta in enumerate(angles.angle_data):
+            ang = - angles.initial_angle - theta
 
-    for i, theta in enumerate(angles.angle_data):
-        ang = - angles.initial_angle - theta
+            rotation_matrix = rotation_matrix_z_from_euler(ang, degrees=degrees)
+
+            vectors[i, :3]  = rotation_matrix.dot(src).reshape(3)
+            vectors[i, 3:6] = rotation_matrix.dot(det).reshape(3)
+            vectors[i, 6:9] = rotation_matrix.dot(row).reshape(3)
+            vectors[i, 9:]  = rotation_matrix.dot(col).reshape(3)
+    # Use the per-projection geometry
+    else:
+        vectors = np.zeros((system.num_positions, 12))
+
+        sign_h = 1
+        sign_v = 1
+
+        if 'right' in panel.origin:
+            sign_h = -1
+        if 'top' in panel.origin:
+            sign_v = -1
 
-        rotation_matrix = rotation_matrix_z_from_euler(ang, degrees=degrees)
+        #for i, (src, det, row, col) in enumerate(zip(system.source.position_set, system.detector.position_set, system.detector.direction_x_set, system.detector.direction_y_set)):
+        for i in range(system.num_positions):
+            vectors[i, :3] = system.source[i].position
+            vectors[i, 3:6]  = system.detector[i].position
+            vectors[i, 6:9] = sign_h * system.detector[i].direction_x * panel.pixel_size[0]
+            vectors[i, 9:]  = sign_v * system.detector[i].direction_y * panel.pixel_size[1]
 
-        vectors[i, :3]  = rotation_matrix.dot(src).reshape(3)
-        vectors[i, 3:6] = rotation_matrix.dot(det).reshape(3)
-        vectors[i, 6:9] = rotation_matrix.dot(row).reshape(3)
-        vectors[i, 9:]  = rotation_matrix.dot(col).reshape(3)
 
     proj_geom = astra.creators.create_proj_geom(projector, panel.num_pixels[1], panel.num_pixels[0], vectors)
     vol_geom = astra.create_vol_geom(volume_geometry_temp.voxel_num_y,
                                     volume_geometry_temp.voxel_num_x,
                                     volume_geometry_temp.voxel_num_z,
-                                    volume_geometry_temp.get_min_x(),
-                                    volume_geometry_temp.get_max_x(),
-                                    volume_geometry_temp.get_min_y(),
-                                    volume_geometry_temp.get_max_y(),
-                                    volume_geometry_temp.get_min_z(),
-                                    volume_geometry_temp.get_max_z()
+                                    volume_geometry_temp.get_min_x() + translation[0],
+                                    volume_geometry_temp.get_max_x() + translation[0],
+                                    volume_geometry_temp.get_min_y() + translation[1],
+                                    volume_geometry_temp.get_max_y() + translation[1],
+                                    volume_geometry_temp.get_min_z() + translation[2],
+                                    volume_geometry_temp.get_max_z() + translation[2]
                                     )
 
 

Wrappers/Python/cil/utilities/display.py

@@ -1075,3 +1075,97 @@ def __init__(self,acquisition_geometry, image_geometry=None, elevation=20, azimu
 
         self.display = _ShowGeometry(acquisition_geometry, image_geometry)
         self.figure = self.display.draw(elev=elevation, azim=azimuthal, view_distance=view_distance, grid=grid, figsize=figsize, fontsize=fontsize)
+
+
+
+class show_SOUV_geometry_vectors(show_base):
+    '''
+    Displays four plots to show i) the source position, 
+    ii) the imager centre, iii) the imager x-direction, and 
+    iv) the imager y-direction for each projection.
+
+
+    Parameters
+    ----------
+    acquisition_geometry: AcquisitionGeometry
+        CIL acquisition geometry
+    figsize: tuple (x, y)
+        Set figure size (inches), default (10,10)
+    fontsize: int
+        Set fontsize, default 10
+
+    Returns
+    -------
+    matplotlib.figure.Figure
+        returns a matplotlib.pyplot figure object
+    '''
+
+
+    def __init__(self, acquisition_geometry:AcquisitionGeometry, figsize=(10,10), fontsize=10):
+
+        # Only applicable for AcquisitionGeometry
+        if not isinstance(acquisition_geometry, AcquisitionGeometry):
+            raise ValueError("The data type of `acquisition_geometry` must be \"<class 'cil.framework.framework.AcquisitionGeometry'>\". It is \"%s\", which is not currently supported by this function." % type(acquisition_geometry))
+
+        # Only applicable for cone_souv geometry type
+        if acquisition_geometry.geom_type != "cone_souv":
+            raise ValueError("The geometry type of `acquisition_geometry` must be \"cone_souv\". It is \"%s\", which is not currently supported by this function." % acquisition_geometry.geom_type)
+
+        self.figure = self._draw(acquisition_geometry, figsize, fontsize);
+
+
+    def _draw(self, acquisition_geometry, figsize, fontsize):
+
+        # Plot the data
+        self.fig, self.axs = plt.subplots(2, 2, figsize=figsize);
+
+        x_axis_values = np.arange(acquisition_geometry.num_projections);
+        i = 0; j = 0;
+        x = 0; y = 1; z = 2;
+        self.axs[j,i].set_title("Source position");
+        self.axs[j,i].plot(x_axis_values, acquisition_geometry.config.system.source.position_set[:,0], label="X axis");
+        self.axs[j,i].plot(x_axis_values, acquisition_geometry.config.system.source.position_set[:,1], label="Y axis");
+        self.axs[j,i].plot(x_axis_values, acquisition_geometry.config.system.source.position_set[:,2], label="Z axis");
+        self.axs[j,i].legend(fontsize=fontsize);
+        # self.axs[j,i].set_xlabel("Projection #");
+        self.axs[j,i].set_ylabel("Position");
+
+        i = 1; j = 0;
+        x += 3; y += 3; z += 3;
+        self.axs[j,i].set_title("Imager Center");
+        self.axs[j,i].plot(x_axis_values, acquisition_geometry.config.system.detector.position_set[:,0], label="X axis");
+        self.axs[j,i].plot(x_axis_values, acquisition_geometry.config.system.detector.position_set[:,1], label="Y axis");
+        self.axs[j,i].plot(x_axis_values, acquisition_geometry.config.system.detector.position_set[:,2], label="Z axis");
+        self.axs[j,i].legend(fontsize=fontsize);
+        # axs[j,i].set_xlabel("Projection #");
+        # axs[j,i].set_ylabel("Position in (cm)");
+
+        i = 0; j = 1;
+        x += 3; y += 3; z += 3;
+        self.axs[j,i].set_title("Imager X-direction");
+        self.axs[j,i].plot(x_axis_values, acquisition_geometry.config.system.detector.direction_x_set[:,0], label="X axis");
+        self.axs[j,i].plot(x_axis_values, acquisition_geometry.config.system.detector.direction_x_set[:,1], label="Y axis");
+        self.axs[j,i].plot(x_axis_values, acquisition_geometry.config.system.detector.direction_x_set[:,2], label="Z axis");
+        self.axs[j,i].legend(fontsize=fontsize);
+        self.axs[j,i].set_xlabel("Projection #");
+        self.axs[j,i].set_ylabel("Position");
+
+        i = 1; j = 1;
+        x += 3; y += 3; z += 3;
+        self.axs[j,i].set_title("Imager Y-direction");
+        self.axs[j,i].plot(x_axis_values, acquisition_geometry.config.system.detector.direction_y_set[:,0], label="X axis");
+        self.axs[j,i].plot(x_axis_values, acquisition_geometry.config.system.detector.direction_y_set[:,1], label="Y axis");
+        self.axs[j,i].plot(x_axis_values, acquisition_geometry.config.system.detector.direction_y_set[:,2], label="Z axis");
+        self.axs[j,i].legend(fontsize=fontsize);
+        self.axs[j,i].set_xlabel("Projection #");
+        # axs[j,i].set_ylabel("Position in (mm)");
+
+        # Resize the text
+        for ax in self.axs.flatten():
+            for item in ([ax.title, ax.xaxis.label, ax.yaxis.label] + ax.get_xticklabels() + ax.get_yticklabels()):
+                item.set_fontsize(fontsize)
+
+        plt.tight_layout()
+        fig2 = plt.gcf()
+        return fig2
+