projection.py 6.53 KB
Newer Older
Stephan Seitz's avatar
Stephan Seitz committed
1
2
3
4
5
6
7
8
9
10
# -*- coding: utf-8 -*-
#
# Copyright © 2019 seitz_local <seitz_local@lmeXX>
#
# Distributed under terms of the GPLv3 license.

"""
Implements a generic forward and backprojection projections
"""

Stephan Seitz's avatar
Stephan Seitz committed
11
12
import sympy

Stephan Seitz's avatar
Stephan Seitz committed
13
import pystencils
Stephan Seitz's avatar
Stephan Seitz committed
14
import pystencils.astnodes
Stephan Seitz's avatar
Stephan Seitz committed
15
import pystencils.autodiff
Stephan Seitz's avatar
Stephan Seitz committed
16
import pystencils.interpolation_astnodes
Stephan Seitz's avatar
Stephan Seitz committed
17
import pystencils_reco._geometry
18
from pystencils_reco import crazy
Stephan Seitz's avatar
Stephan Seitz committed
19
20


21
@crazy
22
23
def forward_projection(volume: pystencils.Field,
                       projection: pystencils.Field,
Stephan Seitz's avatar
Use    
Stephan Seitz committed
24
25
                       projection_matrix,
                       step_size=1,
26
27
28
                       cubic_bspline_interpolation=False,
                       add_to_projector=False):
    # is_projection_stack = projection.spatial_dimensions == volume.spatial_dimensions
Stephan Seitz's avatar
Use    
Stephan Seitz committed
29

Stephan Seitz's avatar
Stephan Seitz committed
30
    interpolation_mode = 'cubic_spline' if cubic_bspline_interpolation else 'linear'
31
    volume_texture = pystencils.interpolation_astnodes.Interpolator(volume,
Stephan Seitz's avatar
Stephan Seitz committed
32
                                                                    interpolation_mode)
33
    ndim = volume.spatial_dimensions
Stephan Seitz's avatar
Stephan Seitz committed
34
35
    projection_matrix = pystencils_reco.ProjectiveMatrix(projection_matrix)

36
37
    t = pystencils_reco.typed_symbols('_parametrization', 'float32')
    texture_coordinates = sympy.Matrix(pystencils_reco.typed_symbols(f'_t:{ndim}', 'float32'))
38
39
    u = projection.physical_coordinates_staggered
    x = volume.index_to_physical(texture_coordinates)
Stephan Seitz's avatar
Stephan Seitz committed
40

41
42
43
    is_perspective = projection_matrix.matrix.cols == ndim + 1

    if is_perspective:
44
        eqn = projection_matrix @ sympy.Matrix([*x, 1]) - sympy.Matrix([*(t * u), t])
45
    else:
Stephan Seitz's avatar
Use    
Stephan Seitz committed
46
47
        # this also works for perspective/cone beam projection (but may lead to instable parametrization)
        eqn = projection_matrix @ x - u
Stephan Seitz's avatar
Stephan Seitz committed
48
    ray_equations = sympy.solve(eqn, texture_coordinates, rational=False)
Stephan Seitz's avatar
Stephan Seitz committed
49

50
51
52
    if not is_perspective:
        t = [t for t in texture_coordinates if t not in ray_equations.keys()][0]
        assert len(ray_equations.keys()) == ndim - 1, "projection_matrix does not appear to define a projection"
53
    ray_equations = sympy.Matrix([ray_equations[s] if s != t else t for s in texture_coordinates])
Stephan Seitz's avatar
Stephan Seitz committed
54

Stephan Seitz's avatar
Stephan Seitz committed
55
    projection_vector = sympy.diff(ray_equations, t)
56
57
    projection_vector_norm = projection_vector.norm()
    projection_vector /= projection_vector_norm
Stephan Seitz's avatar
Stephan Seitz committed
58

Stephan Seitz's avatar
Stephan Seitz committed
59
    conditions = pystencils_reco._geometry.coordinate_in_field_conditions(
60
        volume, ray_equations)
Stephan Seitz's avatar
Stephan Seitz committed
61

Stephan Seitz's avatar
Stephan Seitz committed
62
    central_ray = sympy.Matrix(
63
        projection_matrix.nullspace()[0][:volume.spatial_dimensions])
Stephan Seitz's avatar
Stephan Seitz committed
64
    central_ray /= central_ray.norm()
Stephan Seitz's avatar
Stephan Seitz committed
65

Stephan Seitz's avatar
Stephan Seitz committed
66
    intersection_candidates = []
Stephan Seitz's avatar
Stephan Seitz committed
67
68
    for i in range(ndim):
        solution_min = sympy.solve(ray_equations[i], t, rational=False)
69
        solution_max = sympy.solve(ray_equations[i] - volume.spatial_shape[i],
Stephan Seitz's avatar
Stephan Seitz committed
70
71
72
                                   t,
                                   rational=False)
        intersection_candidates.extend(solution_min + solution_max)
Stephan Seitz's avatar
Stephan Seitz committed
73

74
    intersection_point1 = sympy.Piecewise(
75
        *[(f, sympy.And(*conditions).subs({t: f})) for f in intersection_candidates], (-0, True))
76
    intersection_point2 = sympy.Piecewise(*[(f, sympy.And(*conditions).subs({t: f}))
77
                                            for f in reversed(intersection_candidates)], (-0, True))
78
    assert intersection_point1 != intersection_point2, \
Stephan Seitz's avatar
Use    
Stephan Seitz committed
79
        "The intersections are unconditionally equal, reconstruction volume is not in detector FOV!"
80

Stephan Seitz's avatar
Stephan Seitz committed
81
82
83
84
85
86
87
    # perform a integer set analysis here?
    # space = isl.Space.create_from_names(isl.DEFAULT_CONTEXT, set=[str(t) for t in texture_coordinates])
    # ray_set = isl.BasicSet.universe(space)
    # for i, t in enumerate(texture_coordinates):
    #    # dafaq?
    #    ray_set.add_constraint(isl.Constraint.ineq_from_names(space, {str(texture_coordinates): 1}))
    #    ray_set.add_constraint(isl.Constraint.ineq_from_names(space,
88
    #                                                        # {1: -volume.shape[i],
Stephan Seitz's avatar
Lint    
Stephan Seitz committed
89
    # str(texture_coordinates): -1}))
Stephan Seitz's avatar
Stephan Seitz committed
90
91
    #    ray_set.add_constraint(isl.Constraint.eq_from_name(space, ray_equations[i].subs({ #TODO

92
93
    min_t = sympy.Min(intersection_point1, intersection_point2)
    max_t = sympy.Max(intersection_point1, intersection_point2)
94
95
    # parametrization_dim = list(ray_equations).index(t)
    # min_t = 0
96
    # max_t = volume.spatial_shape[parametrization_dim]
Stephan Seitz's avatar
Stephan Seitz committed
97

Stephan Seitz's avatar
Stephan Seitz committed
98
99
    line_integral, num_steps, min_t_tmp, max_t_tmp, intensity_weighting, step = pystencils.data_types.typed_symbols(
        'line_integral, num_steps, min_t_tmp, max_t_tmp, intensity_weighting, step', 'float32')
Stephan Seitz's avatar
Stephan Seitz committed
100
    i = pystencils.data_types.TypedSymbol('i', 'int32')
101
102
103
104
105
    num_steps = pystencils.data_types.TypedSymbol('num_steps', 'int32')

    # step = step_size / projection_vector_norm
    # tex_coord = ray_equations.subs({t: min_t_tmp + i * step})
    tex_coord = ray_equations.subs({t: min_t_tmp}) + projection_vector * i
Stephan Seitz's avatar
Stephan Seitz committed
106

107
    try:
108
        intensity_weighting_sym = (volume.coordinate_transform(projection_vector)).dot(central_ray) ** 2
109
    except Exception:
110
        intensity_weighting_sym = (volume.coordinate_transform @ projection_vector).dot(central_ray) ** 2
111

Stephan Seitz's avatar
Stephan Seitz committed
112
    assignments = {
Stephan Seitz's avatar
Stephan Seitz committed
113
114
        min_t_tmp: min_t,
        max_t_tmp: max_t,
115
        num_steps: sympy.ceiling((max_t_tmp - min_t_tmp) / (step_size / projection_vector_norm)),
116
        line_integral: sympy.Sum(volume_texture.at(tex_coord),
Stephan Seitz's avatar
Stephan Seitz committed
117
                                 (i, 0, num_steps)),
118
        intensity_weighting: intensity_weighting_sym,
119
120
121
        projection.center(): (line_integral * step_size * intensity_weighting) +
        projection.center() if add_to_projector else 0
        # projection.center(): (max_t_tmp - min_t_tmp) / step # Uncomment to get path length
Stephan Seitz's avatar
Stephan Seitz committed
122
    }
123

Stephan Seitz's avatar
Stephan Seitz committed
124
    # def create_autodiff(self, constant_fields=None):
125
126
    # backward_assignments = backward_projection(AdjointField(projections),
    # AdjointField(volume),
Stephan Seitz's avatar
Stephan Seitz committed
127
128
129
130
    # projection_matrix,
    # 1)
    # self._autodiff = pystencils.autodiff.AutoDiffOp(
    # assignments, "op", constant_fields=constant_fields, backward_assignments=backward_assignments)
131

Stephan Seitz's avatar
Stephan Seitz committed
132
    # assignments._create_autodiff = types.MethodType(create_autodiff, assignments)
133

134
    return assignments
135
136


137
@crazy
138
139
140
def backward_projection(input_projection, output_volume, projection_matrix, normalization):
    projection_matrix = pystencils_reco.ProjectiveMatrix(projection_matrix)
    assignments = pystencils_reco.resampling.generic_spatial_matrix_transform(
141
142
143
144
        input_projection,
        output_volume,
        None,
        inverse_matrix=projection_matrix)
145

146
147
148
149
    for a in assignments.all_assignments:
        a = pystencils.Assignment(a.lhs, a.rhs / normalization)

    return assignments
Stephan Seitz's avatar
Stephan Seitz committed
150
151
152
    a = pystencils.Assignment(a.lhs, a.rhs / normalization)

    return assignments