1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
|
# -*- coding: utf-8 -*-
# This work is part of the Core Imaging Library developed by
# Visual Analytics and Imaging System Group of the Science Technology
# Facilities Council, STFC
# Copyright 2018 Jakob Jorgensen, Daniil Kazantsev and Edoardo Pasca
# 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
# 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.
import numpy
from ccpi.optimisation.ops import Operator, PowerMethodNonsquare
from ccpi.framework import ImageData, DataContainer , \
ImageGeometry, AcquisitionGeometry
from ccpi.plugins.processors import CCPiBackwardProjector, \
CCPiForwardProjector , setupCCPiGeometries
class CCPiProjectorSimple(Operator):
"""ASTRA projector modified to use DataSet and geometry."""
def __init__(self, geomv, geomp, default=False):
super(CCPiProjectorSimple, self).__init__()
# Store volume and sinogram geometries.
self.acquisition_geometry = geomp
self.volume_geometry = geomv
if geomp.geom_type == "cone":
raise TypeError('Can only handle parallel beam')
# set-up the geometries if compatible
geoms = setupCCPiGeometries(geomv.voxel_num_x,geomv.voxel_num_y,
geomv.voxel_num_z, geomp.angles, 0)
vg = ImageGeometry(voxel_num_x=geoms['output_volume_x'],
voxel_num_y=geoms['output_volume_y'],
voxel_num_z=geoms['output_volume_z'])
pg = AcquisitionGeometry('parallel',
'3D',
geomp.angles,
geoms['n_h'], geomp.pixel_size_h,
geoms['n_v'], geomp.pixel_size_v #2D in 3D is a slice 1 pixel thick
)
if not default:
# check if geometry is the same (on the voxels)
if not ( vg.voxel_num_x == geomv.voxel_num_x and \
vg.voxel_num_y == geomv.voxel_num_y and \
vg.voxel_num_z == geomv.voxel_num_z ):
msg = 'The required volume geometry will not work\nThe following would\n'
msg += vg.__str__()
raise ValueError(msg)
if not (pg.pixel_num_h == geomp.pixel_num_h and \
pg.pixel_num_v == geomp.pixel_num_v and \
len( pg.angles ) == len( geomp.angles ) ) :
msg = 'The required acquisition geometry will not work\nThe following would\n'
msg += pg.__str__()
raise ValueError(msg)
self.fp = CCPiForwardProjector(image_geometry=vg,
acquisition_geometry=pg,
output_axes_order=['angle','vertical','horizontal'])
self.bp = CCPiBackwardProjector(image_geometry=vg,
acquisition_geometry=pg,
output_axes_order=['horizontal_x','horizontal_y','vertical'])
# Initialise empty for singular value.
self.s1 = None
def direct(self, image_data):
self.fp.set_input(image_data)
out = self.fp.get_output()
return out
def adjoint(self, acquisition_data):
self.bp.set_input(acquisition_data)
out = self.bp.get_output()
return out
#def delete(self):
# astra.data2d.delete(self.proj_id)
def get_max_sing_val(self):
a = PowerMethodNonsquare(self,10)
self.s1 = a[0]
return self.s1
def size(self):
# Only implemented for 3D
return ( (self.acquisition_geometry.angles.size, \
self.acquisition_geometry.pixel_num_v,
self.acquisition_geometry.pixel_num_h), \
(self.volume_geometry.voxel_num_x, \
self.volume_geometry.voxel_num_y,
self.volume_geometry.voxel_num_z) )
def create_image_data(self):
x0 = ImageData(geometry = self.volume_geometry,
dimension_labels=self.bp.output_axes_order)#\
#.subset(['horizontal_x','horizontal_y','vertical'])
x0.fill(numpy.random.randn(*x0.shape))
return x0
|
