diff --git a/src/pystencils_reco/projection.py b/src/pystencils_reco/projection.py index 83a785d78c23cdecee6814d4f85751f0f21200fa..3d5f8da6be03a1d62733cd248c86e1b6b5bb866e 100644 --- a/src/pystencils_reco/projection.py +++ b/src/pystencils_reco/projection.py @@ -24,7 +24,8 @@ def forward_projection(volume: pystencils.Field, projection_matrix, step_size=1, cubic_bspline_interpolation=False, - add_to_projector=False): + add_to_projector=False, + central_ray_point=None): # is_projection_stack = projection.spatial_dimensions == volume.spatial_dimensions interpolation_mode = 'cubic_spline' if cubic_bspline_interpolation else 'linear' @@ -59,9 +60,10 @@ def forward_projection(volume: pystencils.Field, conditions = pystencils_reco._geometry.coordinate_in_field_conditions( volume, ray_equations) - central_ray = sympy.Matrix( - projection_matrix.nullspace()[0][:volume.spatial_dimensions]) - central_ray /= central_ray.norm() + if not central_ray_point: + central_ray_point = [0] * projection.spatial_dimensions + central_ray = projection_vector.subs({i: j for i, j in zip( + pystencils.x_vector(projection.spatial_dimensions), central_ray_point)}) intersection_candidates = [] for i in range(ndim): @@ -104,10 +106,10 @@ def forward_projection(volume: pystencils.Field, # tex_coord = ray_equations.subs({t: min_t_tmp + i * step}) tex_coord = ray_equations.subs({t: min_t_tmp}) + projection_vector * i - try: - intensity_weighting_sym = (volume.coordinate_transform(projection_vector)).dot(central_ray) ** 2 - except Exception: - intensity_weighting_sym = (volume.coordinate_transform @ projection_vector).dot(central_ray) ** 2 + if callable(volume.coordinate_transform): + intensity_weighting_sym = projection_vector.dot(central_ray) ** 2 + else: + intensity_weighting_sym = projection_vector.dot(central_ray) ** 2 assignments = { min_t_tmp: min_t, @@ -117,7 +119,7 @@ def forward_projection(volume: pystencils.Field, (i, 0, num_steps)), intensity_weighting: intensity_weighting_sym, projection.center(): (line_integral * step_size * intensity_weighting) + - projection.center() if add_to_projector else 0 + (projection.center() if add_to_projector else 0) # projection.center(): (max_t_tmp - min_t_tmp) / step # Uncomment to get path length } diff --git a/tests/test_projection.py b/tests/test_projection.py index 1238d559a34cf2de694c3161790ee5c277f0afbf..522454877b28bc9bd5f803e2d5b4291b1b37d231 100644 --- a/tests/test_projection.py +++ b/tests/test_projection.py @@ -17,6 +17,7 @@ import sympy import pystencils import pystencils_reco +from pystencils.data_types import create_type from pystencils_reco.projection import forward_projection try: @@ -46,13 +47,12 @@ def test_genereric_projection(): projection_matrix = pystencils_reco.matrix_symbols('T', pystencils.data_types.create_type('float32'), 3, 4) assignments = forward_projection(volume, projections, projection_matrix) - print(assignments) kernel = assignments.compile('gpu') pystencils.show_code(kernel) def test_projection_cpu(): - volume = pystencils.fields('volume: float32[100,200,300]') + volume = pystencils.fields('volume: float32[3d]') projections = pystencils.fields('projections: float32[2D]') projection_matrix = sympy.Matrix([[-289.0098977737411, -1205.2274801832275, 0.0, 186000.0],