summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorVaggelis Papoutsellis <22398586+epapoutsellis@users.noreply.github.com>2019-10-15 15:04:41 +0100
committerEdoardo Pasca <edo.paskino@gmail.com>2019-10-15 15:04:41 +0100
commitcd0b0bf9e9be2f5aa2135d589ca7d07d8a73a761 (patch)
tree5ad8790514674d2eff4e9bcf32d7debc98914eee
parentb0c8d3c31fc5d58e74da4745169b2426d2206975 (diff)
downloadastra-wrapper-cd0b0bf9e9be2f5aa2135d589ca7d07d8a73a761.tar.gz
astra-wrapper-cd0b0bf9e9be2f5aa2135d589ca7d07d8a73a761.tar.bz2
astra-wrapper-cd0b0bf9e9be2f5aa2135d589ca7d07d8a73a761.tar.xz
astra-wrapper-cd0b0bf9e9be2f5aa2135d589ca7d07d8a73a761.zip
Add fbp (#36)
* fix norm for MC operator * fix norm for MC operator * fix demo_astra_simple * fix norm * add 3DMC * add tests * fix 3D case * fix labels 2DMC * test for all AstraOps * fix 2DMC * minor fix * FBP 2D * check 3D FBP * add FBP FDK * merge from fix_Astra_ops * Delete Untitled-checkpoint.ipynb * Delete Untitled1-checkpoint.ipynb * fix tests * fix tests * fix tests * fix tests * new FBP with tests * add FBM for 2D/3D MC * remove demos * remove demos * updata scaling 3D FDK cone * add warning form FanFlat and GPU case * remove comments and old FBP from operators * remove old demo * remove old fbp * remove old imports
-rw-r--r--Wrappers/Python/ccpi/astra/operators/AstraProjector3DMC.py116
-rw-r--r--Wrappers/Python/ccpi/astra/operators/AstraProjector3DSimple.py49
-rw-r--r--Wrappers/Python/ccpi/astra/operators/AstraProjectorMC.py84
-rw-r--r--Wrappers/Python/ccpi/astra/operators/AstraProjectorSimple.py66
-rw-r--r--Wrappers/Python/ccpi/astra/operators/__init__.py1
-rw-r--r--Wrappers/Python/ccpi/astra/processors/AstraBackProjector.py2
-rw-r--r--Wrappers/Python/ccpi/astra/processors/AstraBackProjector3D.py2
-rw-r--r--Wrappers/Python/ccpi/astra/processors/AstraBackProjectorMC.py4
-rw-r--r--Wrappers/Python/ccpi/astra/processors/AstraForwardProjector.py2
-rw-r--r--Wrappers/Python/ccpi/astra/processors/AstraForwardProjector3D.py4
-rw-r--r--Wrappers/Python/ccpi/astra/processors/AstraForwardProjectorMC.py3
-rw-r--r--Wrappers/Python/ccpi/astra/processors/FBP.py349
-rw-r--r--Wrappers/Python/ccpi/astra/processors/__init__.py1
13 files changed, 617 insertions, 66 deletions
diff --git a/Wrappers/Python/ccpi/astra/operators/AstraProjector3DMC.py b/Wrappers/Python/ccpi/astra/operators/AstraProjector3DMC.py
new file mode 100644
index 0000000..b7b52ef
--- /dev/null
+++ b/Wrappers/Python/ccpi/astra/operators/AstraProjector3DMC.py
@@ -0,0 +1,116 @@
+# -*- coding: utf-8 -*-
+# This work is independent part of the Core Imaging Library developed by
+# Visual Analytics and Imaging System Group of the Science Technology
+# Facilities Council, STFC
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+from ccpi.optimisation.operators import LinearOperator
+from ccpi.astra.operators import AstraProjector3DSimple
+import numpy as np
+
+
+class AstraProjector3DMC(LinearOperator):
+
+ """ASTRA projector modified to use DataSet and geometry."""
+ def __init__(self, geomv, geomp):
+ super(AstraProjector3DMC, self).__init__()
+
+ # Store volume and sinogram geometries.
+ self.sinogram_geometry = geomp
+ self.volume_geometry = geomv
+
+ self.igtmp3D = self.volume_geometry.clone()
+ self.igtmp3D.channels = 1
+ self.igtmp3D.shape = self.volume_geometry.shape[1:]
+ self.igtmp3D.dimension_labels = ['vertical', 'horizontal_y', 'horizontal_x']
+
+ self.agtmp3D = self.sinogram_geometry.clone()
+ self.agtmp3D.channels = 1
+ self.agtmp3D.shape = self.sinogram_geometry.shape[1:]
+ self.agtmp3D.dimension_labels = ['vertical', 'angle', 'horizontal']
+
+ self.A3D = AstraProjector3DSimple(self.igtmp3D, self.agtmp3D)
+
+ self.s1 = None
+ self.channels = self.volume_geometry.channels
+
+ def direct(self, x, out = None):
+
+ if out is None:
+
+ tmp = self.sinogram_geometry.allocate()
+ for k in range(self.channels):
+ t = self.A3D.direct(x.subset(channel=k))
+ np.copyto(tmp.array[k], t.array)
+ return tmp
+
+ else:
+
+ for k in range(self.channels):
+ tmp = self.A3D.direct(x.subset(channel=k))
+ np.copyto(out.array[k], tmp.array)
+
+ def adjoint(self, x, out = None):
+
+ if out is None:
+
+ tmp = self.volume_geometry.allocate()
+ for k in range(self.channels):
+ t = self.A3D.adjoint(x.subset(channel=k))
+ np.copyto(tmp.array[k], t.array)
+ return tmp
+
+ else:
+
+ for k in range(self.channels):
+ tmp = self.A3D.adjoint(x.subset(channel=k))
+ np.copyto(out.array[k], tmp.array)
+
+ def domain_geometry(self):
+ return self.volume_geometry
+
+ def range_geometry(self):
+ return self.sinogram_geometry
+
+ def calculate_norm(self):
+
+ return self.A3D.norm()
+
+
+
+if __name__ == '__main__':
+
+ from ccpi.framework import ImageGeometry, AcquisitionGeometry
+
+ N = 30
+ channels = 5
+ angles = np.linspace(0, np.pi, 180)
+ ig = ImageGeometry(N, N, N, channels = channels)
+ ag = AcquisitionGeometry('parallel','3D', angles, pixel_num_h = N, pixel_num_v=5, channels = channels)
+
+ A3DMC = AstraProjector3DMC(ig, ag)
+ z = A3DMC.norm()
+
+# igtmp3D = A3DMC.volume_geometry.clone()
+# igtmp3D.channels = 1
+# igtmp3D.shape = A3DMC.volume_geometry.shape[1:]
+# igtmp3D.dimension_labels = ['vertical', 'horizontal_y', 'horizontal_x']
+#
+# agtmp3D = A3DMC.sinogram_geometry.clone()
+# agtmp3D.channels = 1
+# agtmp3D.shape = A3DMC.sinogram_geometry.shape[1:]
+# agtmp3D.dimension_labels = ['vertical', 'angle', 'horizontal']
+#
+# A3D = AstraProjector3DSimple(igtmp3D, agtmp3D)
+# z = A3D.norm() \ No newline at end of file
diff --git a/Wrappers/Python/ccpi/astra/operators/AstraProjector3DSimple.py b/Wrappers/Python/ccpi/astra/operators/AstraProjector3DSimple.py
index fabfe33..a139f64 100644
--- a/Wrappers/Python/ccpi/astra/operators/AstraProjector3DSimple.py
+++ b/Wrappers/Python/ccpi/astra/operators/AstraProjector3DSimple.py
@@ -15,19 +15,27 @@
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
-from ccpi.optimisation.operators import Operator, LinearOperator
-from ccpi.framework import AcquisitionData, ImageData, DataContainer
-from ccpi.astra.processors import AstraForwardProjector, AstraBackProjector, \
- AstraForwardProjector3D, AstraBackProjector3D
+from ccpi.optimisation.operators import LinearOperator
+from ccpi.framework import ImageGeometry, AcquisitionGeometry
+from ccpi.astra.processors import AstraForwardProjector3D, AstraBackProjector3D
+import numpy as np
class AstraProjector3DSimple(LinearOperator):
+
"""ASTRA projector modified to use DataSet and geometry."""
def __init__(self, geomv, geomp):
+
super(AstraProjector3DSimple, self).__init__()
# Store volume and sinogram geometries.
- self.sinogram_geometry = geomp
- self.volume_geometry = geomv
+ # The order of the ouput sinogram is not the default one from the acquistion geometry
+ # The order of the backprojection is the default one from the image geometry
+
+ geomp.dimension_labels = ['vertical','angle','horizontal']
+ geomp.shape = (geomp.pixel_num_v, len(geomp.angles), geomp.pixel_num_h)
+
+ self.sinogram_geometry = geomp
+ self.volume_geometry = geomv
self.fp = AstraForwardProjector3D(volume_geometry=geomv,
sinogram_geometry=geomp,
@@ -36,7 +44,7 @@ class AstraProjector3DSimple(LinearOperator):
self.bp = AstraBackProjector3D(volume_geometry=geomv,
sinogram_geometry=geomp,
output_axes_order=['vertical','horizontal_y','horizontal_x'])
-
+
# Initialise empty for singular value.
self.s1 = None
@@ -60,9 +68,30 @@ class AstraProjector3DSimple(LinearOperator):
return self.volume_geometry
def range_geometry(self):
- return self.sinogram_geometry
-
+ return self.sinogram_geometry
+
def norm(self):
x0 = self.volume_geometry.allocate('random')
self.s1, sall, svec = LinearOperator.PowerMethod(self, 50, x0)
- return self.s1 \ No newline at end of file
+ return self.s1
+
+
+if __name__ == '__main__':
+
+ N = 30
+ angles = np.linspace(0, np.pi, 180)
+ ig = ImageGeometry(N, N, N)
+ ag = AcquisitionGeometry('parallel','3D', angles, pixel_num_h = N, pixel_num_v=5)
+
+ A = AstraProjector3DSimple(ig, ag)
+ print(A.norm())
+
+ x = ig.allocate('random_int')
+ sin = A.direct(x)
+
+ y = ag.allocate('random_int')
+ im = A.adjoint(y)
+
+
+
+ \ No newline at end of file
diff --git a/Wrappers/Python/ccpi/astra/operators/AstraProjectorMC.py b/Wrappers/Python/ccpi/astra/operators/AstraProjectorMC.py
index 9a982a5..fb0375f 100644
--- a/Wrappers/Python/ccpi/astra/operators/AstraProjectorMC.py
+++ b/Wrappers/Python/ccpi/astra/operators/AstraProjectorMC.py
@@ -1,25 +1,29 @@
-# -*- coding: utf-8 -*-
-# This work is independent part of the Core Imaging Library developed by
-# Visual Analytics and Imaging System Group of the Science Technology
-# Facilities Council, STFC
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
+#========================================================================
+# Copyright 2019 Science Technology Facilities Council
+# Copyright 2019 University of Manchester
#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU General Public License for more details.
+# This work is part of the Core Imaging Library developed by Science Technology
+# Facilities Council and University of Manchester
#
-# You should have received a copy of the GNU General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0.txt
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+#=========================================================================
+
+
+from ccpi.optimisation.operators import LinearOperator
+from ccpi.astra.processors import AstraForwardProjectorMC, AstraBackProjectorMC
+from ccpi.astra.operators import AstraProjectorSimple
-from ccpi.optimisation.operators import Operator, LinearOperator
-from ccpi.framework import AcquisitionData, ImageData, DataContainer
-from ccpi.astra.processors import AstraForwardProjector, AstraBackProjector, \
- AstraForwardProjectorMC, AstraBackProjectorMC, AstraForwardProjector3D, \
- AstraBackProjector3D
class AstraProjectorMC(LinearOperator):
"""ASTRA Multichannel projector"""
@@ -42,6 +46,7 @@ class AstraProjectorMC(LinearOperator):
# Initialise empty for singular value.
self.s1 = None
+ self.device = device
def direct(self, IM, out=None):
self.fp.set_input(IM)
@@ -63,9 +68,40 @@ class AstraProjectorMC(LinearOperator):
return self.volume_geometry
def range_geometry(self):
- return self.sinogram_geometry
+ return self.sinogram_geometry
- def norm(self):
- x0 = self.volume_geometry.allocate('random')
- self.s1, sall, svec = LinearOperator.PowerMethod(self, 50, x0)
- return self.s1
+ def calculate_norm(self):
+
+ igtmp = self.volume_geometry.clone()
+ igtmp.shape = self.volume_geometry.shape[1:]
+ igtmp.dimension_labels = ['horizontal_y', 'horizontal_x']
+ igtmp.channels = 1
+
+ agtmp = self.sinogram_geometry.clone()
+ agtmp.shape = self.sinogram_geometry.shape[1:]
+ agtmp.dimension_labels = ['angle', 'horizontal']
+ agtmp.channels = 1
+
+ Atmp = AstraProjectorSimple(igtmp, agtmp, device = self.device)
+
+ return Atmp.norm()
+
+
+
+if __name__ == '__main__':
+
+ from ccpi.framework import ImageGeometry, AcquisitionGeometry
+ import numpy as np
+
+ N = 30
+ angles = np.linspace(0, np.pi, 180)
+ ig = ImageGeometry(N, N, channels = 5)
+ ag = AcquisitionGeometry('parallel','2D', angles, pixel_num_h = N, channels = 5)
+
+ A = AstraProjectorMC(ig, ag, 'gpu')
+ print(A.norm())
+
+
+
+
+
diff --git a/Wrappers/Python/ccpi/astra/operators/AstraProjectorSimple.py b/Wrappers/Python/ccpi/astra/operators/AstraProjectorSimple.py
index 618ec96..a1d6ef6 100644
--- a/Wrappers/Python/ccpi/astra/operators/AstraProjectorSimple.py
+++ b/Wrappers/Python/ccpi/astra/operators/AstraProjectorSimple.py
@@ -1,24 +1,29 @@
-# -*- coding: utf-8 -*-
-# This work is independent part of the Core Imaging Library developed by
-# Visual Analytics and Imaging System Group of the Science Technology
-# Facilities Council, STFC
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
+#========================================================================
+# Copyright 2019 Science Technology Facilities Council
+# Copyright 2019 University of Manchester
#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU General Public License for more details.
+# This work is part of the Core Imaging Library developed by Science Technology
+# Facilities Council and University of Manchester
#
-# You should have received a copy of the GNU General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0.txt
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+#=========================================================================
-from ccpi.optimisation.operators import Operator, LinearOperator
-from ccpi.framework import AcquisitionData, ImageData, DataContainer
+from ccpi.optimisation.operators import LinearOperator
from ccpi.astra.processors import AstraForwardProjector, AstraBackProjector
+
+
class AstraProjectorSimple(LinearOperator):
"""ASTRA projector modified to use DataSet and geometry."""
def __init__(self, geomv, geomp, device):
@@ -30,14 +35,14 @@ class AstraProjectorSimple(LinearOperator):
self.fp = AstraForwardProjector(volume_geometry=geomv,
sinogram_geometry=geomp,
- proj_id=None,
+ proj_id = None,
device=device)
- self.bp = AstraBackProjector(volume_geometry=geomv,
- sinogram_geometry=geomp,
- proj_id=None,
- device=device)
-
+ self.bp = AstraBackProjector(volume_geometry = geomv,
+ sinogram_geometry = geomp,
+ proj_id = None,
+ device = device)
+
# Initialise empty for singular value.
self.s1 = None
@@ -62,8 +67,23 @@ class AstraProjectorSimple(LinearOperator):
def range_geometry(self):
return self.sinogram_geometry
-
+
def norm(self):
x0 = self.volume_geometry.allocate('random')
self.s1, sall, svec = LinearOperator.PowerMethod(self, 50, x0)
return self.s1
+
+
+if __name__ == '__main__':
+
+ from ccpi.framework import ImageGeometry, AcquisitionGeometry
+ import numpy as np
+
+ N = 30
+ angles = np.linspace(0, np.pi, 180)
+ ig = ImageGeometry(N, N)
+ ag = AcquisitionGeometry('parallel','2D', angles, pixel_num_h = N)
+ A = AstraProjectorSimple(ig, ag, 'cpu')
+ print(A.norm())
+
+
diff --git a/Wrappers/Python/ccpi/astra/operators/__init__.py b/Wrappers/Python/ccpi/astra/operators/__init__.py
index 6848776..58028a7 100644
--- a/Wrappers/Python/ccpi/astra/operators/__init__.py
+++ b/Wrappers/Python/ccpi/astra/operators/__init__.py
@@ -1,4 +1,5 @@
from .AstraProjectorSimple import AstraProjectorSimple
from .AstraProjector3DSimple import AstraProjector3DSimple
+from .AstraProjector3DMC import AstraProjector3DMC
from .AstraProjectorMC import AstraProjectorMC \ No newline at end of file
diff --git a/Wrappers/Python/ccpi/astra/processors/AstraBackProjector.py b/Wrappers/Python/ccpi/astra/processors/AstraBackProjector.py
index d094e79..14ece80 100644
--- a/Wrappers/Python/ccpi/astra/processors/AstraBackProjector.py
+++ b/Wrappers/Python/ccpi/astra/processors/AstraBackProjector.py
@@ -1,4 +1,4 @@
-from ccpi.framework import DataProcessor, ImageData, AcquisitionData
+from ccpi.framework import DataProcessor, ImageData
from ccpi.astra.utils import convert_geometry_to_astra
import astra
diff --git a/Wrappers/Python/ccpi/astra/processors/AstraBackProjector3D.py b/Wrappers/Python/ccpi/astra/processors/AstraBackProjector3D.py
index 26a549a..2e03af3 100644
--- a/Wrappers/Python/ccpi/astra/processors/AstraBackProjector3D.py
+++ b/Wrappers/Python/ccpi/astra/processors/AstraBackProjector3D.py
@@ -1,4 +1,4 @@
-from ccpi.framework import DataProcessor, ImageData, AcquisitionData
+from ccpi.framework import DataProcessor, ImageData
from ccpi.astra.utils import convert_geometry_to_astra
import astra
diff --git a/Wrappers/Python/ccpi/astra/processors/AstraBackProjectorMC.py b/Wrappers/Python/ccpi/astra/processors/AstraBackProjectorMC.py
index ef8fe5d..ef63fa5 100644
--- a/Wrappers/Python/ccpi/astra/processors/AstraBackProjectorMC.py
+++ b/Wrappers/Python/ccpi/astra/processors/AstraBackProjectorMC.py
@@ -1,5 +1,5 @@
-from ccpi.framework import DataProcessor, ImageData, AcquisitionData
-from ccpi.astra.utils import convert_geometry_to_astra
+from ccpi.framework import ImageData
+
from ccpi.astra.processors import AstraBackProjector
diff --git a/Wrappers/Python/ccpi/astra/processors/AstraForwardProjector.py b/Wrappers/Python/ccpi/astra/processors/AstraForwardProjector.py
index 2cf5716..8c8ea83 100644
--- a/Wrappers/Python/ccpi/astra/processors/AstraForwardProjector.py
+++ b/Wrappers/Python/ccpi/astra/processors/AstraForwardProjector.py
@@ -1,4 +1,4 @@
-from ccpi.framework import DataProcessor, ImageData, AcquisitionData
+from ccpi.framework import DataProcessor, AcquisitionData
from ccpi.astra.utils import convert_geometry_to_astra
import astra
diff --git a/Wrappers/Python/ccpi/astra/processors/AstraForwardProjector3D.py b/Wrappers/Python/ccpi/astra/processors/AstraForwardProjector3D.py
index 18b4078..77a5858 100644
--- a/Wrappers/Python/ccpi/astra/processors/AstraForwardProjector3D.py
+++ b/Wrappers/Python/ccpi/astra/processors/AstraForwardProjector3D.py
@@ -1,4 +1,4 @@
-from ccpi.framework import DataProcessor, ImageData, AcquisitionData
+from ccpi.framework import DataProcessor, AcquisitionData
from ccpi.astra.utils import convert_geometry_to_astra
import astra
@@ -27,7 +27,6 @@ class AstraForwardProjector3D(DataProcessor):
'output_axes_order' : output_axes_order
}
- #DataProcessor.__init__(self, **kwargs)
super(AstraForwardProjector3D, self).__init__(**kwargs)
self.set_ImageGeometry(volume_geometry)
@@ -58,6 +57,7 @@ class AstraForwardProjector3D(DataProcessor):
self.vol_geom = vol_geom
def process(self, out=None):
+
IM = self.get_input()
DATA = AcquisitionData(geometry=self.sinogram_geometry,
dimension_labels=self.output_axes_order)
diff --git a/Wrappers/Python/ccpi/astra/processors/AstraForwardProjectorMC.py b/Wrappers/Python/ccpi/astra/processors/AstraForwardProjectorMC.py
index 34f4222..55ef4cd 100644
--- a/Wrappers/Python/ccpi/astra/processors/AstraForwardProjectorMC.py
+++ b/Wrappers/Python/ccpi/astra/processors/AstraForwardProjectorMC.py
@@ -1,5 +1,4 @@
-from ccpi.framework import DataProcessor, ImageData, AcquisitionData
-from ccpi.astra.utils import convert_geometry_to_astra
+from ccpi.framework import AcquisitionData
from ccpi.astra.processors import AstraForwardProjector
diff --git a/Wrappers/Python/ccpi/astra/processors/FBP.py b/Wrappers/Python/ccpi/astra/processors/FBP.py
new file mode 100644
index 0000000..4405af6
--- /dev/null
+++ b/Wrappers/Python/ccpi/astra/processors/FBP.py
@@ -0,0 +1,349 @@
+from ccpi.framework import DataProcessor, ImageGeometry, AcquisitionGeometry
+from ccpi.astra.utils import convert_geometry_to_astra
+import astra
+import numpy as np
+import warnings
+
+class FBP(DataProcessor):
+
+ '''Filtered Back Projection is a reconstructor
+
+ Input: Volume Geometry
+ Sinogram Geometry
+ Filter_type
+ Device = cpu/gpu. For 2D cases we have the option cpu/gpu
+ For 3D cases we have the option of gpu
+
+ Cases: 1) Parallel 3D FBP: using 2D FBP per slice ( CPU/GPU )
+
+
+ Example: FBP(ig, ag, 'ram-lak', 'cpu')
+ FBP.set_input(sinogram)
+ reconstruction = FBP.get_ouput()
+
+ Output: ImageData
+
+
+ '''
+
+ def __init__(self, volume_geometry,
+ sinogram_geometry,
+ device = 'cpu',
+ filter_type = 'ram-lak',
+ **kwargs):
+
+
+ if sinogram_geometry.dimension == '3D':
+
+
+ # For 3D FBP and parallel goemetry, we perform 2D FBP per slice
+ # Therofore, we need vol_geom2D, porj_geom2D attributes to setup in the DataProcessor.
+ # There are both CPU/GPU options for this case
+
+ if sinogram_geometry.geom_type=='parallel':
+
+ # By default, we select GPU device
+ super(FBP, self).__init__(volume_geometry = volume_geometry,
+ sinogram_geometry = sinogram_geometry,
+ device = 'gpu', proj_id = None,
+ vol_geom2D = None, proj_geom2D = None,
+ filter_type = filter_type)
+
+ # we need a 2D image and acquisition geometry
+ ig2D = ImageGeometry(voxel_num_x = self.volume_geometry.voxel_num_x,
+ voxel_num_y = self.volume_geometry.voxel_num_y,
+ voxel_size_x = self.volume_geometry.voxel_size_x,
+ voxel_size_y = self.volume_geometry.voxel_size_y)
+
+ ag2D = AcquisitionGeometry(geom_type = self.sinogram_geometry.geom_type,
+ dimension = '2D',
+ angles = self.sinogram_geometry.angles,
+ pixel_num_h=self.sinogram_geometry.pixel_num_h,
+ pixel_size_h=self.sinogram_geometry.pixel_size_h)
+
+ # Convert to Astra Geometries
+ self.vol_geom2D, self.proj_geom2D = convert_geometry_to_astra(ig2D,ag2D)
+
+ if self.device == 'gpu':
+ self.set_projector(astra.create_projector('cuda', self.proj_geom2D, self.vol_geom2D) )
+ else:
+ self.set_projector(astra.create_projector('line', self.proj_geom2D, self.vol_geom2D) )
+
+
+ elif sinogram_geometry.geom_type=='cone':
+
+ # For 3D cone we do not need ASTRA proj_id
+ super(FBP, self).__init__(volume_geometry = volume_geometry,
+ sinogram_geometry = sinogram_geometry,
+ device = 'gpu',
+ filter_type = 'ram-lak')
+
+ #Raise an error if the user select a filter other than ram-lak
+ if self.filter_type !='ram-lak':
+ raise NotImplementedError('Currently in astra, FDK has only ram-lak available')
+ else:
+
+
+ # Setup for 2D case
+ super(FBP, self).__init__(volume_geometry = volume_geometry,
+ sinogram_geometry = sinogram_geometry,
+ device = device, proj_id = None,
+ filter_type = filter_type)
+
+ # Raise an error if the user select a filter other than ram-lak, for 2D case
+ if self.filter_type !='ram-lak':
+ raise NotImplementedError('Currently in astra, 2D FBP is using only the ram-lak filter, switch to gpu for other filters')
+
+ if (self.sinogram_geometry.geom_type == 'cone' and self.device == 'gpu') and self.sinogram_geometry.channels>=1:
+ warnings.warn("For FanFlat geometry there is a mismatch between CPU & GPU filtered back projection. Currently, only CPU case provides a reasonable result.")
+
+ self.set_ImageGeometry(volume_geometry)
+ self.set_AcquisitionGeometry(sinogram_geometry)
+
+ # Set up ASTRA Volume and projection geometry, not to be stored in self
+ vol_geom, proj_geom = convert_geometry_to_astra(self.volume_geometry,
+ self.sinogram_geometry)
+
+ # Here we define proj_id that will be used for ASTRA
+ # ASTRA projector, to be stored
+ if self.device == 'cpu':
+
+ # Note that 'line' only one option
+ if self.sinogram_geometry.geom_type == 'parallel':
+ self.set_projector(astra.create_projector('line', proj_geom, vol_geom) )
+
+ elif self.sinogram_geometry.geom_type == 'cone':
+ self.set_projector(astra.create_projector('line_fanflat', proj_geom, vol_geom) )
+
+ else:
+ NotImplemented
+
+ elif self.device == 'gpu':
+
+ self.set_projector(astra.create_projector('cuda', proj_geom, vol_geom) )
+
+ else:
+ NotImplemented
+
+
+
+ def check_input(self, dataset):
+
+ if self.sinogram_geometry.dimension == '2D':
+ if dataset.number_of_dimensions == 2 or dataset.geometry.channels>1:
+ return True
+ else:
+ raise ValueError("Expected input dimensions is 2, got {0}"\
+ .format(dataset.number_of_dimensions))
+
+ elif self.sinogram_geometry.dimension=='3D':
+ if dataset.number_of_dimensions == 3 or dataset.geometry.channels>1:
+ return True
+ else:
+ raise ValueError("Expected input dimensions is 3, got {0}"\
+ .format(dataset.number_of_dimensions))
+
+ def set_projector(self, proj_id):
+ self.proj_id = proj_id
+
+ def set_ImageGeometry(self, volume_geometry):
+ self.volume_geometry = volume_geometry
+
+ def set_AcquisitionGeometry(self, sinogram_geometry):
+ self.sinogram_geometry = sinogram_geometry
+
+ def set_filter(self, filter_type):
+ self.filter_type = filter_type
+
+ def process(self):
+
+ # Get DATA
+ DATA = self.get_input()
+
+ # Allocate FBP reconstruction
+ IM = self.volume_geometry.allocate()
+
+ # Set up ASTRA Volume and projection geometry, not to be stored in self
+ vol_geom, proj_geom = convert_geometry_to_astra(self.volume_geometry,
+ self.sinogram_geometry)
+ # Get number of channels
+ num_chan = DATA.geometry.channels
+
+ #######################################################################
+ #######################################################################
+ #######################################################################
+
+ # Run FBP for one channel data
+ if num_chan == 1:
+
+ if self.sinogram_geometry.dimension == '2D':
+
+ # Create a data object for the reconstruction and to hold sinogram for 2D
+ rec_id = astra.data2d.create( '-vol', vol_geom)
+ sinogram_id = astra.data2d.create('-sino', proj_geom, DATA.as_array())
+
+ # ASTRA configuration for reconstruction algorithm
+ if self.device == 'cpu':
+ cfg = astra.astra_dict('FBP')
+ cfg['ProjectorId'] = self.proj_id
+ else:
+ cfg = astra.astra_dict('FBP_CUDA')
+ cfg['ReconstructionDataId'] = rec_id
+ cfg['ProjectionDataId'] = sinogram_id
+ cfg['FilterType'] = self.filter_type
+ alg_id = astra.algorithm.create(cfg)
+ astra.algorithm.run(alg_id)
+
+ # Get the result
+ IM.array = astra.data2d.get(rec_id)
+ astra.data2d.delete(rec_id)
+ astra.data2d.delete(sinogram_id)
+ astra.algorithm.delete(alg_id)
+
+ if self.sinogram_geometry.geom_type == 'parallel':
+ if self.device == 'cpu':
+ return IM / (self.volume_geometry.voxel_size_x**2)
+ else:
+ scaling = self.volume_geometry.voxel_size_x
+ return scaling * IM
+ else:
+ if self.device == 'cpu':
+ return IM / (self.volume_geometry.voxel_size_x**2)
+ else:
+ return IM
+
+
+
+ elif self.sinogram_geometry.dimension == '3D':
+
+ # Here we use 2D geometries to run FBP for parallel geometry for each slice
+ if self.sinogram_geometry.geom_type == 'parallel':
+
+ rec_id = astra.data2d.create( '-vol', self.vol_geom2D)
+
+ for i in range(DATA.shape[0]):
+
+ sinogram_id = astra.data2d.create('-sino', self.proj_geom2D, DATA.as_array()[i])
+
+ cfg = astra.astra_dict('FBP_CUDA')
+ cfg['ReconstructionDataId'] = rec_id
+ cfg['ProjectionDataId'] = sinogram_id
+ cfg['FilterType'] = self.filter_type
+ alg_id = astra.algorithm.create(cfg)
+ astra.algorithm.run(alg_id)
+
+ # since we run gpu we need to rescale and flip it
+ IM.array[i] = np.flip(astra.data2d.get(rec_id) * self.volume_geometry.voxel_size_x,0)
+
+ astra.algorithm.delete(alg_id)
+ astra.data2d.delete(rec_id)
+ astra.data2d.delete(sinogram_id)
+
+ return IM
+
+ elif self.sinogram_geometry.geom_type == 'cone':
+
+ rec_id = astra.data3d.create('-vol', vol_geom)
+ sinogram_id = astra.data3d.create('-sino', proj_geom, DATA.as_array())
+
+ cfg = astra.astra_dict('FDK_CUDA')
+ cfg['ReconstructionDataId'] = rec_id
+ cfg['ProjectionDataId'] = sinogram_id
+ alg_id = astra.algorithm.create(cfg)
+ astra.algorithm.run(alg_id)
+ IM.array = astra.data3d.get(rec_id)
+ astra.data3d.delete(rec_id)
+ astra.data3d.delete(sinogram_id)
+ astra.algorithm.delete(alg_id)
+ return IM/(self.volume_geometry.voxel_size_x**4)
+
+ #######################################################################
+ #######################################################################
+ #######################################################################
+
+ # Here we preform FBP for each channel for the 2D/3D cases
+ else:
+
+ if self.sinogram_geometry.dimension == '2D':
+
+ if self.device == 'cpu':
+ cfg = astra.astra_dict('FBP')
+ cfg['ProjectorId'] = self.proj_id
+ else:
+ cfg = astra.astra_dict('FBP_CUDA')
+ cfg['FilterType'] = self.filter_type
+
+ for i in range(num_chan):
+
+ sinogram_id = astra.data2d.create('-sino', proj_geom, DATA.as_array()[i])
+ rec_id = astra.data2d.create( '-vol', vol_geom)
+
+ cfg['ReconstructionDataId'] = rec_id
+ cfg['ProjectionDataId'] = sinogram_id
+
+ alg_id = astra.algorithm.create(cfg)
+ astra.algorithm.run(alg_id)
+
+ # Get the result
+ IM.array[i] = astra.data2d.get(rec_id)
+
+
+ if self.device == 'cpu':
+ IM.array[i] /= (self.volume_geometry.voxel_size_x**2)
+ else:
+ IM.array[i] *= self.volume_geometry.voxel_size_x
+ return IM
+
+ elif self.sinogram_geometry.dimension == '3D':
+
+ if self.sinogram_geometry.geom_type == 'parallel':
+
+ cfg = astra.astra_dict('FBP_CUDA')
+ cfg['FilterType'] = self.filter_type
+
+ for i in range(num_chan):
+
+ for j in range(DATA.shape[1]):
+
+ sinogram_id = astra.data2d.create('-sino', self.proj_geom2D, DATA.as_array()[i,j])
+ rec_id = astra.data2d.create( '-vol', self.vol_geom2D)
+
+ cfg['ReconstructionDataId'] = rec_id
+ cfg['ProjectionDataId'] = sinogram_id
+
+ alg_id = astra.algorithm.create(cfg)
+ astra.algorithm.run(alg_id)
+
+ # Get the result
+ IM.array[i,j] = np.flip(astra.data2d.get(rec_id) * self.volume_geometry.voxel_size_x,0)
+
+ astra.algorithm.delete(alg_id)
+ astra.data2d.delete(rec_id)
+ astra.data2d.delete(sinogram_id)
+
+ return IM
+
+ elif self.sinogram_geometry.geom_type == 'cone':
+
+ for i in range(num_chan):
+
+ rec_id = astra.data3d.create('-vol', vol_geom)
+ sinogram_id = astra.data3d.create('-sino', proj_geom, DATA.as_array()[i])
+
+ cfg = astra.astra_dict('FDK_CUDA')
+ cfg['ReconstructionDataId'] = rec_id
+ cfg['ProjectionDataId'] = sinogram_id
+ alg_id = astra.algorithm.create(cfg)
+ astra.algorithm.run(alg_id)
+ IM.array[i] = astra.data3d.get(rec_id) / self.volume_geometry.voxel_size_x ** 4
+
+ astra.data3d.delete(rec_id)
+ astra.data3d.delete(sinogram_id)
+ astra.algorithm.delete(alg_id)
+
+ return IM
+
+
+
+
+
diff --git a/Wrappers/Python/ccpi/astra/processors/__init__.py b/Wrappers/Python/ccpi/astra/processors/__init__.py
index 6bfd7f5..d4cc90b 100644
--- a/Wrappers/Python/ccpi/astra/processors/__init__.py
+++ b/Wrappers/Python/ccpi/astra/processors/__init__.py
@@ -5,3 +5,4 @@ from .AstraForwardProjectorMC import AstraForwardProjectorMC
from .AstraBackProjectorMC import AstraBackProjectorMC
from .AstraForwardProjector3D import AstraForwardProjector3D
from .AstraBackProjector3D import AstraBackProjector3D
+from .FBP import FBP \ No newline at end of file