work for precond pdhg

5 files changed, 93 insertions, 5 deletions

diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
index ee8f609..19da3d4 100755
--- a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
@@ -219,5 +219,55 @@ class BlockOperator(Operator):
         shape = (self.shape[1], 1)
         return BlockGeometry(*[el.range_geometry() for el in self.operators],
                     shape=shape)
+        
+    def sum_abs_row(self):
+        
+        res = []
+        for row in range(self.shape[0]):
+            for col in range(self.shape[1]):
+                if col == 0:
+                    prod = self.get_item(row,col).sum_abs_row()
+                else:
+                    prod += self.get_item(row,col).sum_abs_row()
+            res.append(prod)
+            
+        if self.shape[1]==1:
+            tmp = sum(res)
+            return ImageData(tmp)
+        else:
+            return BlockDataContainer(*res)
+                 
+        
 if __name__ == '__main__':
-    pass
+    
+    from ccpi.framework import ImageGeometry
+    from ccpi.optimisation.operators import Gradient, Identity, SparseFiniteDiff
+
+    from ccpi.framework import AcquisitionData, ImageData, BlockDataContainer
+    from ccpi.optimisation.operators import Operator, LinearOperator
+    
+    
+    M, N= 4, 3
+    ig = ImageGeometry(M, N)
+    arr = ig.allocate('random_int')    
+    G = Gradient(ig)
+    Id = Identity(ig)
+    
+    B = BlockOperator(G, Id)
+    
+    print(B.sum_abs_row().as_array())
+    
+    Gx = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
+    Gy = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
+    
+    d1 = abs(Gx.matrix()).toarray().sum(axis=0)
+    d2 = abs(Gy.matrix()).toarray().sum(axis=0)
+    d3 = abs(Id.matrix()).toarray().sum(axis=0)
+    
+    d_res = numpy.reshape(d1 + d2 + d3, ig.shape, 'F')
+    
+    print(d_res)
+    
+    
+    
+    
diff --git a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py
index 0f50e82..df6c076 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py
@@ -7,6 +7,9 @@ Created on Wed Mar  6 19:30:51 2019
 """
 
 from ccpi.optimisation.operators import LinearOperator
+import scipy.sparse as sp
+import numpy as np
+from ccpi.framework import ImageData
 
 
 class Identity(LinearOperator):
@@ -39,4 +42,38 @@ class Identity(LinearOperator):
         return self.gm_domain
         
     def range_geometry(self):
-        return self.gm_range
\ No newline at end of file
+        return self.gm_range
+    
+    def matrix(self):
+        
+        return sp.eye(np.prod(self.gm_domain.shape))
+    
+    def sum_abs_row(self):
+        
+        return ImageData(np.array(np.reshape(abs(self.matrix()).sum(axis=0), self.gm_domain.shape, 'F')))
+ 
+    def sum_abs_col(self):
+        
+        return ImageData(np.array(np.reshape(abs(self.matrix()).sum(axis=1), self.gm_domain.shape, 'F')))
+            
+    
+if __name__ == '__main__':
+    
+    from ccpi.framework import ImageGeometry
+
+    M, N= 2, 3
+    ig = ImageGeometry(M, N)
+    arr = ig.allocate('random_int')
+    
+    Id = Identity(ig)
+    d = Id.matrix()
+    print(d.toarray())
+    
+    d1 = Id.sum_abs_col()
+    print(d1.as_array())
+    
+    
+
+            
+    
+    
\ No newline at end of file
diff --git a/Wrappers/Python/ccpi/optimisation/operators/__init__.py b/Wrappers/Python/ccpi/optimisation/operators/__init__.py
index 57d89ad..63c1320 100755
--- a/Wrappers/Python/ccpi/optimisation/operators/__init__.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/__init__.py
@@ -11,6 +11,7 @@ from .ScaledOperator import ScaledOperator
 from .BlockOperator import BlockOperator

 from .BlockScaledOperator import BlockScaledOperator

 

+from .SparseFiniteDiff import SparseFiniteDiff

 

 from .FiniteDifferenceOperator import FiniteDiff

 from .GradientOperator import Gradient

@@ -18,4 +19,4 @@ from .SymmetrizedGradientOperator import SymmetrizedGradient
 from .IdentityOperator import Identity

 from .ZeroOperator import ZeroOp

 

-from .SparseFiniteDiff import SparseFiniteDiff

+

diff --git a/Wrappers/Python/wip/pdhg_TV_tomography2D.py b/Wrappers/Python/wip/pdhg_TV_tomography2D.py
index 52b7922..8bd63c2 100644
--- a/Wrappers/Python/wip/pdhg_TV_tomography2D.py
+++ b/Wrappers/Python/wip/pdhg_TV_tomography2D.py
@@ -55,7 +55,7 @@ detectors = 150
 angles = np.linspace(0,np.pi,100)
 
 ag = AcquisitionGeometry('parallel','2D',angles, detectors)
-Aop = AstraProjectorSimple(ig, ag, 'cpu')
+Aop = AstraProjectorSimple(ig, ag, 'gpu')
 sin = Aop.direct(data)
 
 plt.imshow(sin.as_array())
diff --git a/Wrappers/Python/wip/pdhg_TV_tomography2D_time.py b/Wrappers/Python/wip/pdhg_TV_tomography2D_time.py
index 7ac1566..e9a85cc 100644
--- a/Wrappers/Python/wip/pdhg_TV_tomography2D_time.py
+++ b/Wrappers/Python/wip/pdhg_TV_tomography2D_time.py
@@ -66,7 +66,7 @@ detectors = 150
 angles = np.linspace(0,np.pi,100)
 
 ag = AcquisitionGeometry('parallel','2D',angles, detectors, channels = np.shape(phantom_2Dt)[0])
-Aop = AstraProjectorMC(ig, ag, 'cpu')
+Aop = AstraProjectorMC(ig, ag, 'gpu')
 sin = Aop.direct(data)
 
 plt.imshow(sin.as_array()[10])