Geometry (#31)

* started geometry

* Simple geometry implemented #2

* sino geometry simplified to standard

* Added get min/max methods vol geom

* Basic vol and sino geoms working in astra_op and simple_demo

* Fixed style

* Added center/offset to vol geom.

3 files changed, 168 insertions, 3 deletions

diff --git a/Wrappers/Python/ccpi/reconstruction/astra_ops.py b/Wrappers/Python/ccpi/reconstruction/astra_ops.py
index 0c95b73..bc6f5e1 100644
--- a/Wrappers/Python/ccpi/reconstruction/astra_ops.py
+++ b/Wrappers/Python/ccpi/reconstruction/astra_ops.py
@@ -21,6 +21,69 @@ import numpy
 from ccpi.framework import SinogramData, VolumeData
 from ccpi.reconstruction.ops import PowerMethodNonsquare
 
+class AstraProjectorSimple(Operator):
+    """ASTRA projector modified to use DataSet and geometry."""
+    def __init__(self, geomv, geomp, device):
+        super(AstraProjectorSimple, self).__init__()
+        
+        # ASTRA Volume geometry
+        self.vol_geom = astra.create_vol_geom(geomv.voxel_num_x, \
+                                              geomv.voxel_num_y, \
+                                              geomv.getMinX(), \
+                                              geomv.getMaxX(), \
+                                              geomv.getMinY(), \
+                                              geomv.getMaxY())
+        
+        # ASTRA Projections geometry
+        if geomp.dimension == '2D':
+            if geomp.geom_type == 'parallel':
+                self.proj_geom = astra.create_proj_geom('parallel', \
+                                    geomp.pixel_size_h, \
+                                    geomp.pixel_num_h, \
+                                    geomp.angles)
+            elif geomp.geom_type == 'cone':
+                NotImplemented
+        else:
+            NotImplemented
+        
+        # ASTRA projector
+        if device == 'cpu':
+            # Note that 'line' is only for parallel (2D) and only one option
+            self.proj_id = astra.create_projector('line', self.proj_geom, 
+                                                  self.vol_geom) # for CPU
+        elif device == 'gpu':
+            self.proj_id = astra.create_projector('cuda', self.proj_geom, 
+                                                  self.vol_geom) # for GPU
+        else:
+            NotImplemented
+        
+        self.s1 = None
+    
+    def direct(self, IM):
+        
+        sinogram_id, DATA = astra.create_sino(IM.as_array(), self.proj_id)
+        astra.data2d.delete(sinogram_id)
+        return SinogramData(DATA)
+    
+    def adjoint(self, DATA):
+        rec_id, IM = astra.create_backprojection(DATA.as_array(), self.proj_id)
+        astra.data2d.delete(rec_id)
+        return VolumeData(IM)
+    
+    def delete(self):
+        astra.data2d.delete(self.proj_id)
+    
+    def get_max_sing_val(self):
+        self.s1, sall, svec = PowerMethodNonsquare(self,10)
+        return self.s1
+    
+    def size(self):
+        return ( (self.proj_geom['ProjectionAngles'].size, \
+                  self.proj_geom['DetectorCount']), \
+                 (self.vol_geom['GridColCount'], \
+                  self.vol_geom['GridRowCount']) )
+    
+
 class AstraProjector(Operator):
     """A simple 2D/3D parallel/fan beam projection/backprojection class based 
     on ASTRA toolbox"""
diff --git a/Wrappers/Python/ccpi/reconstruction/geoms.py b/Wrappers/Python/ccpi/reconstruction/geoms.py
new file mode 100644
index 0000000..f7f3329
--- /dev/null
+++ b/Wrappers/Python/ccpi/reconstruction/geoms.py
@@ -0,0 +1,96 @@
+
+class VolumeGeometry:
+    
+    def __init__(self, \
+                 voxel_num_x=None, \
+                 voxel_num_y=None, \
+                 voxel_num_z=None, \
+                 voxel_size_x=None, \
+                 voxel_size_y=None, \
+                 voxel_size_z=None, \
+                 center_x=0, \
+                 center_y=0, \
+                 center_z=0):
+        
+        self.voxel_num_x = voxel_num_x
+        self.voxel_num_y = voxel_num_y
+        self.voxel_num_z = voxel_num_z
+        self.voxel_size_x = voxel_size_x
+        self.voxel_size_y = voxel_size_y
+        self.voxel_size_z = voxel_size_z
+        self.center_x = center_x
+        self.center_y = center_y
+        self.center_z = center_z  
+        
+    def getMinX(self):
+        return self.center_x - 0.5*self.voxel_num_x*self.voxel_size_x
+        
+    def getMaxX(self):
+        return self.center_x + 0.5*self.voxel_num_x*self.voxel_size_x
+        
+    def getMinY(self):
+        return self.center_y - 0.5*self.voxel_num_y*self.voxel_size_y
+        
+    def getMaxY(self):
+        return self.center_y + 0.5*self.voxel_num_y*self.voxel_size_y
+        
+    def getMinZ(self):
+        return self.center_z - 0.5*self.voxel_num_z*self.voxel_size_z
+        
+    def getMaxZ(self):
+        return self.center_z + 0.5*self.voxel_num_z*self.voxel_size_z
+        
+    
+class SinogramGeometry:
+    
+    def __init__(self, \
+                 geom_type, \
+                 dimension, \
+                 angles, \
+                 pixel_num_h=None, \
+                 pixel_size_h=None, \
+                 pixel_num_v=None, \
+                 pixel_size_v=None, \
+                 dist_source_center=None, \
+                 dist_center_detector=None, \
+                 ):
+        """
+        General inputs for standard type projection geometries
+        detectorDomain or detectorpixelSize:
+            If 2D
+                If scalar: Width of detector or single detector pixel
+                If 2-vec: Error
+            If 3D
+                If scalar: Width in both dimensions
+                If 2-vec: Vertical then horizontal size
+        grid
+            If 2D
+                If scalar: number of detectors
+                If 2-vec: error
+            If 3D
+                If scalar: Square grid that size
+                If 2-vec vertical then horizontal size
+        cone or parallel
+        2D or 3D
+        parallel_parameters: ?
+        cone_parameters:
+            source_to_center_dist (if parallel: NaN)
+            center_to_detector_dist (if parallel: NaN)
+        standard or nonstandard (vec) geometry
+        angles
+        angles_format radians or degrees
+        """
+        self.geom_type = geom_type   # 'parallel' or 'cone'
+        self.dimension = dimension # 2D or 3D
+        self.angles = angles
+        
+        self.dist_source_center = dist_source_center
+        self.dist_center_detector = dist_center_detector
+        
+        self.pixel_num_h = pixel_num_h
+        self.pixel_size_h = pixel_size_h
+        self.pixel_num_v = pixel_num_v
+        self.pixel_size_v = pixel_size_v
+
+        
+                
diff --git a/Wrappers/Python/test/simple_demo.py b/Wrappers/Python/test/simple_demo.py
index 4b7e89e..7a6b049 100644
--- a/Wrappers/Python/test/simple_demo.py
+++ b/Wrappers/Python/test/simple_demo.py
@@ -8,6 +8,7 @@ from ccpi.reconstruction.algs import *
 from ccpi.reconstruction.funcs import *
 from ccpi.reconstruction.ops import *
 from ccpi.reconstruction.astra_ops import *
+from ccpi.reconstruction.geoms import *
 
 import numpy as np
 import matplotlib.pyplot as plt
@@ -22,10 +23,10 @@ x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 2.0
 plt.imshow(x)
 plt.show()
 
-#vg = VolumeGeometry(grid=(N,N), domain=((-N/2,N/2),(-N/2,N/2)))
+vg = VolumeGeometry(N,N,None, 1,1,None)
 
-#Phantom = VolumeData(x,geometry=vg)
 Phantom = VolumeData(x)
+#Phantom = VolumeData(x)
 
 # Set up measurement geometry
 angles_num = 20; # angles number
@@ -37,14 +38,19 @@ det_num = N
 SourceOrig = 500
 OrigDetec = 0
 
+# Parallelbeam geometry test
+pg = SinogramGeometry('parallel','2D',angles,det_num,det_w)
+
 # Set up ASTRA projector
 #Aop = AstraProjector(vg, angles, N,'gpu')
 #Aop = AstraProjector(det_w, det_num, angles, projtype='parallel',device='gpu')
 
-Aop = AstraProjector(det_w, det_num, SourceOrig, 
+Aop_old = AstraProjector(det_w, det_num, SourceOrig, 
                      OrigDetec, angles, 
                      N,'fanbeam','gpu') 
 
+Aop = AstraProjectorSimple(vg, pg, 'gpu')
+
 
 # Try forward and backprojection
 b = Aop.direct(Phantom)