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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
|
/************************************************************************
*
* Copyright (C) 2025 IRCAD France
*
* This file is part of Sight.
*
* Sight is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Sight 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 Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with Sight. If not, see <https://www.gnu.org/licenses/>.
*
***********************************************************************/
#pragma once
#include <data/helper/medical_image.hpp>
#include <data/image.hpp>
#include <geometry/__/matrix.hpp>
#include <glm/ext/matrix_transform.hpp>
#include <glm/glm.hpp>
#include <glm/gtc/type_ptr.hpp>
#define GLM_ENABLE_EXPERIMENTAL
#include <glm/gtx/vec_swizzle.hpp>
#undef GLM_ENABLE_EXPERIMENTAL
namespace sight::geometry::data
{
/**
* @brief Return the transform matrix that convert image voxels in world coordinates, taking into account origin,
* orientation and spacing.
*
* @param _image reference image
* @return transform matrix
*/
template<typename M = glm::dmat4>
requires std::is_same_v<M, glm::mat4>|| std::is_same_v<M, glm::dmat4>
inline M image_to_world_transform(const sight::data::image& _image)
{
const auto origin = glm::make_vec3(_image.origin().data());
const auto spacing = glm::make_vec3(_image.spacing().data());
const auto orientation = _image.orientation();
const M transform {
orientation[0], orientation[3], orientation[6], 0.,
orientation[1], orientation[4], orientation[7], 0.,
orientation[2], orientation[5], orientation[8], 0.,
origin[0], origin[1], origin[2], 1
};
return glm::scale(transform, glm::vec<3, typename M::value_type>(spacing));
}
/**
* @brief Return the transform matrix that convert world coordinates into image voxels coordinates, taking into account
* origin, orientation and spacing.
*
* @param _image reference image
* @return transform matrix
*/
template<typename M = glm::dmat4>
requires std::is_same_v<M, glm::mat4>|| std::is_same_v<M, glm::dmat4>
inline M world_to_image_transform(const sight::data::image& _image)
{
const auto origin = _image.origin();
const auto orientation = _image.orientation();
const M transform {
orientation[0], orientation[3], orientation[6], 0.,
orientation[1], orientation[4], orientation[7], 0.,
orientation[2], orientation[5], orientation[8], 0.,
origin[0], origin[1], origin[2], 1
};
const auto spacing = glm::make_vec3(_image.spacing().data());
SIGHT_ASSERT("Spacing is null", spacing[0] != 0. && spacing[1] != 0. && spacing[2] != 0.);
const auto scale = glm::scale(glm::identity<M>(), glm::vec<3, typename M::value_type>(1.0 / spacing));
return scale * geometry::inverse_translation_rotation(transform);
}
/**
* @brief Convert to image voxel coordinates
*
* @tparam W input container type (default glm::dvec3), I output container type (default glm::ivec3)
* @param _image reference image
* @param _world image coordinates
* @param _round if true, round final computed values using std::round so 0.9877.. will be 1.
* This mitigated floating point precision issues
* @param _clamp if true, returns the closest valid image coordinates, if the world coordinates are out of bounds
* @return image coordinates
*/
template<typename I = glm::ivec3, typename W = glm::dvec3>
I world_to_image(const sight::data::image& _image, const W& _world, bool _round = false, bool _clamp = false)
{
const glm::dvec4 world_vector {
glm::dvec4::value_type(_world[0]),
glm::dvec4::value_type(_world[1]),
glm::dvec4::value_type(_world[2]),
glm::dvec4::value_type(1)
};
const glm::dmat4 transform = world_to_image_transform(_image);
const auto image_vector = transform * world_vector;
auto voxel = glm::xyz(_round ? glm::round(image_vector) : image_vector);
if(_clamp)
{
const glm::dvec3 size = glm::make_vec3(_image.size().data());
voxel = glm::clamp(voxel, glm::dvec3(0.0), size - 1.0);
}
const I res = {
typename I::value_type(voxel[0]),
typename I::value_type(voxel[1]),
typename I::value_type(voxel[2])
};
return res;
}
/**
* @brief Convert to world coordinates
*
* @tparam I input container type (default glm::ivec3), W output container type (default glm::dvec3)
* @param _image reference image
* @param _voxel image coordinates
* @return world coordinates
*/
template<typename W = glm::dvec3, typename I = glm::ivec3>
W image_to_world(const sight::data::image& _image, const I& _voxel)
{
const glm::dvec4 image_vector = {
glm::dvec4::value_type(_voxel[0]),
glm::dvec4::value_type(_voxel[1]),
glm::dvec4::value_type(_voxel[2]),
glm::dvec4::value_type(1)
};
const glm::dmat4 transform = image_to_world_transform(_image);
const auto world_vector = transform * image_vector;
return {world_vector.x, world_vector.y, world_vector.z};
}
/**
* @brief Helper function to calculate the slice index of a given fiducial point in a specified orientation within a
* medical image.
*
* @param _image : The input image as a constant reference to data::image.
* @param _point : The coordinates of the fiducial point as a std::array of three doubles.
* @param _axis : The orientation (axial, sagittal, or frontal) to calculate the slice index.
* @return std::optional<std::int64_t> : The calculated slice index as an integer.
*/
template<typename V>
std::optional<std::int64_t> get_fiducial_slice_index(
const sight::data::image& _image,
const V& _point,
sight::data::helper::medical_image::axis_t _axis
)
{
const auto point_in_image = world_to_image(_image, _point, true);
const auto slice_index = point_in_image[_axis];
if(slice_index < 0 || slice_index >= std::int64_t(_image.size()[_axis]))
{
return std::nullopt;
}
return slice_index;
}
} // namespace sight::geometry::data
|
