Mercurial > hg > orthanc-stone

comparison Framework/Volumes/VolumeImageGeometry.cpp @ 742:fa5febe0f0c2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
moved OrientedBoundingBox in the Volumes folder
author Sebastien Jodogne <s.jodogne@gmail.com>
date Wed, 22 May 2019 09:41:03 +0200
parents Framework/Toolbox/VolumeImageGeometry.cpp@93a8949a1ef7
children e8fdf29cd0ca
comparison
equal deleted inserted replaced
741:c1d6a566dfd3 742:fa5febe0f0c2
1 /**
2 * Stone of Orthanc
3 * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics
4 * Department, University Hospital of Liege, Belgium
5 * Copyright (C) 2017-2019 Osimis S.A., Belgium
6 *
7 * This program is free software: you can redistribute it and/or
8 * modify it under the terms of the GNU Affero General Public License
9 * as published by the Free Software Foundation, either version 3 of
10 * the License, or (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful, but
13 * WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Affero General Public License for more details.
16 *
17 * You should have received a copy of the GNU Affero General Public License
18 * along with this program. If not, see <http://www.gnu.org/licenses/>.
19 **/
20
21
22 #include "VolumeImageGeometry.h"
23
24 #include "../Toolbox/GeometryToolbox.h"
25
26 #include <Core/OrthancException.h>
27
28
29 namespace OrthancStone
30 {
31 void VolumeImageGeometry::Invalidate()
32 {
33 Vector p = (axialGeometry_.GetOrigin() +
34 static_cast<double>(depth_ - 1) * voxelDimensions_[2] * axialGeometry_.GetNormal());
35
36 coronalGeometry_ = CoordinateSystem3D(p,
37 axialGeometry_.GetAxisX(),
38 -axialGeometry_.GetNormal());
39
40 sagittalGeometry_ = CoordinateSystem3D(p,
41 axialGeometry_.GetAxisY(),
42 axialGeometry_.GetNormal());
43
44 Vector origin = (
45 axialGeometry_.MapSliceToWorldCoordinates(-0.5 * voxelDimensions_[0],
46 -0.5 * voxelDimensions_[1]) -
47 0.5 * voxelDimensions_[2] * axialGeometry_.GetNormal());
48
49 Vector scaling;
50
51 if (width_ == 0 ||
52 height_ == 0 ||
53 depth_ == 0)
54 {
55 LinearAlgebra::AssignVector(scaling, 1, 1, 1);
56 }
57 else
58 {
59 scaling = (
60 axialGeometry_.GetAxisX() * voxelDimensions_[0] * static_cast<double>(width_) +
61 axialGeometry_.GetAxisY() * voxelDimensions_[1] * static_cast<double>(height_) +
62 axialGeometry_.GetNormal() * voxelDimensions_[2] * static_cast<double>(depth_));
63 }
64
65 transform_ = LinearAlgebra::Product(
66 GeometryToolbox::CreateTranslationMatrix(origin[0], origin[1], origin[2]),
67 GeometryToolbox::CreateScalingMatrix(scaling[0], scaling[1], scaling[2]));
68
69 LinearAlgebra::InvertMatrix(transformInverse_, transform_);
70 }
71
72
73 VolumeImageGeometry::VolumeImageGeometry() :
74 width_(0),
75 height_(0),
76 depth_(0)
77 {
78 LinearAlgebra::AssignVector(voxelDimensions_, 1, 1, 1);
79 Invalidate();
80 }
81
82
83 void VolumeImageGeometry::SetSize(unsigned int width,
84 unsigned int height,
85 unsigned int depth)
86 {
87 width_ = width;
88 height_ = height;
89 depth_ = depth;
90 Invalidate();
91 }
92
93
94 void VolumeImageGeometry::SetAxialGeometry(const CoordinateSystem3D& geometry)
95 {
96 axialGeometry_ = geometry;
97 Invalidate();
98 }
99
100
101 void VolumeImageGeometry::SetVoxelDimensions(double x,
102 double y,
103 double z)
104 {
105 if (x <= 0 ||
106 y <= 0 ||
107 z <= 0)
108 {
109 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
110 }
111 else
112 {
113 LinearAlgebra::AssignVector(voxelDimensions_, x, y, z);
114 Invalidate();
115 }
116 }
117
118
119 const CoordinateSystem3D& VolumeImageGeometry::GetProjectionGeometry(VolumeProjection projection) const
120 {
121 switch (projection)
122 {
123 case VolumeProjection_Axial:
124 return axialGeometry_;
125
126 case VolumeProjection_Coronal:
127 return coronalGeometry_;
128
129 case VolumeProjection_Sagittal:
130 return sagittalGeometry_;
131
132 default:
133 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
134 }
135 }
136
137
138 Vector VolumeImageGeometry::GetVoxelDimensions(VolumeProjection projection) const
139 {
140 switch (projection)
141 {
142 case VolumeProjection_Axial:
143 return voxelDimensions_;
144
145 case VolumeProjection_Coronal:
146 return LinearAlgebra::CreateVector(voxelDimensions_[0], voxelDimensions_[2], voxelDimensions_[1]);
147
148 case VolumeProjection_Sagittal:
149 return LinearAlgebra::CreateVector(voxelDimensions_[1], voxelDimensions_[2], voxelDimensions_[0]);
150
151 default:
152 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
153 }
154 }
155
156
157 unsigned int VolumeImageGeometry::GetProjectionWidth(VolumeProjection projection) const
158 {
159 switch (projection)
160 {
161 case VolumeProjection_Axial:
162 return width_;
163
164 case VolumeProjection_Coronal:
165 return width_;
166
167 case VolumeProjection_Sagittal:
168 return height_;
169
170 default:
171 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
172 }
173 }
174
175
176 unsigned int VolumeImageGeometry::GetProjectionHeight(VolumeProjection projection) const
177 {
178 switch (projection)
179 {
180 case VolumeProjection_Axial:
181 return height_;
182
183 case VolumeProjection_Coronal:
184 return depth_;
185
186 case VolumeProjection_Sagittal:
187 return depth_;
188
189 default:
190 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
191 }
192 }
193
194
195 unsigned int VolumeImageGeometry::GetProjectionDepth(VolumeProjection projection) const
196 {
197 switch (projection)
198 {
199 case VolumeProjection_Axial:
200 return depth_;
201
202 case VolumeProjection_Coronal:
203 return height_;
204
205 case VolumeProjection_Sagittal:
206 return width_;
207
208 default:
209 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
210 }
211 }
212
213
214 Vector VolumeImageGeometry::GetCoordinates(float x,
215 float y,
216 float z) const
217 {
218 Vector p = LinearAlgebra::Product(transform_, LinearAlgebra::CreateVector(x, y, z, 1));
219
220 assert(LinearAlgebra::IsNear(p[3], 1)); // Affine transform, no perspective effect
221
222 // Back to non-homogeneous coordinates
223 return LinearAlgebra::CreateVector(p[0], p[1], p[2]);
224 }
225
226
227 bool VolumeImageGeometry::DetectProjection(VolumeProjection& projection,
228 const Vector& planeNormal) const
229 {
230 if (GeometryToolbox::IsParallel(planeNormal, axialGeometry_.GetNormal()))
231 {
232 projection = VolumeProjection_Axial;
233 return true;
234 }
235 else if (GeometryToolbox::IsParallel(planeNormal, coronalGeometry_.GetNormal()))
236 {
237 projection = VolumeProjection_Coronal;
238 return true;
239 }
240 else if (GeometryToolbox::IsParallel(planeNormal, sagittalGeometry_.GetNormal()))
241 {
242 projection = VolumeProjection_Sagittal;
243 return true;
244 }
245 else
246 {
247 return false;
248 }
249 }
250
251
252 bool VolumeImageGeometry::DetectSlice(VolumeProjection& projection,
253 unsigned int& slice,
254 const CoordinateSystem3D& plane) const
255 {
256 if (!DetectProjection(projection, plane.GetNormal()))
257 {
258 return false;
259 }
260
261 // Transforms the coordinates of the origin of the plane, into the
262 // coordinates of the axial geometry
263 const Vector& origin = plane.GetOrigin();
264 Vector p = LinearAlgebra::Product(
265 transformInverse_,
266 LinearAlgebra::CreateVector(origin[0], origin[1], origin[2], 1));
267
268 assert(LinearAlgebra::IsNear(p[3], 1));
269
270 double z;
271
272 switch (projection)
273 {
274 case VolumeProjection_Axial:
275 z = p[2];
276 break;
277
278 case VolumeProjection_Coronal:
279 z = p[1];
280 break;
281
282 case VolumeProjection_Sagittal:
283 z = p[0];
284 break;
285
286 default:
287 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
288 }
289
290 const unsigned int projectionDepth = GetProjectionDepth(projection);
291
292 z *= static_cast<double>(projectionDepth);
293 if (z < 0)
294 {
295 return false;
296 }
297
298 unsigned int d = static_cast<unsigned int>(std::floor(z));
299 if (d >= projectionDepth)
300 {
301 return false;
302 }
303 else
304 {
305 slice = d;
306 return true;
307 }
308 }
309 }