Mercurial > hg > orthanc-stone

comparison OrthancStone/Sources/Volumes/VolumeImageGeometry.cpp @ 1512:244ad1e4e76a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
reorganization of folders
author Sebastien Jodogne <s.jodogne@gmail.com>
date Tue, 07 Jul 2020 16:21:02 +0200
parents Framework/Volumes/VolumeImageGeometry.cpp@30deba7bc8e2
children 8563ea5d8ae4
comparison
equal deleted inserted replaced
1511:9dfeee74c1e6 1512:244ad1e4e76a
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-2020 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 <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 LOG(TRACE) << "VolumeImageGeometry::Invalidate() origin = " << origin(0) << "," << origin(1) << "," << origin(2) << " | width_ = " << width_ << " | height_ = " << height_ << " | depth_ = " << depth_;
50
51 Vector scaling;
52
53 if (width_ == 0 ||
54 height_ == 0 ||
55 depth_ == 0)
56 {
57 LinearAlgebra::AssignVector(scaling, 1, 1, 1);
58 }
59 else
60 {
61 scaling = (
62 axialGeometry_.GetAxisX() * voxelDimensions_[0] * static_cast<double>(width_) +
63 axialGeometry_.GetAxisY() * voxelDimensions_[1] * static_cast<double>(height_) +
64 axialGeometry_.GetNormal() * voxelDimensions_[2] * static_cast<double>(depth_));
65 }
66
67 transform_ = LinearAlgebra::Product(
68 GeometryToolbox::CreateTranslationMatrix(origin[0], origin[1], origin[2]),
69 GeometryToolbox::CreateScalingMatrix(scaling[0], scaling[1], scaling[2]));
70
71 LinearAlgebra::InvertMatrix(transformInverse_, transform_);
72 }
73
74
75 VolumeImageGeometry::VolumeImageGeometry() :
76 width_(0),
77 height_(0),
78 depth_(0)
79 {
80 LinearAlgebra::AssignVector(voxelDimensions_, 1, 1, 1);
81 Invalidate();
82 }
83
84
85 void VolumeImageGeometry::SetSizeInVoxels(unsigned int width,
86 unsigned int height,
87 unsigned int depth)
88 {
89 width_ = width;
90 height_ = height;
91 depth_ = depth;
92 Invalidate();
93 }
94
95
96 void VolumeImageGeometry::SetAxialGeometry(const CoordinateSystem3D& geometry)
97 {
98 axialGeometry_ = geometry;
99 Invalidate();
100 }
101
102
103 void VolumeImageGeometry::SetVoxelDimensions(double x,
104 double y,
105 double z)
106 {
107 if (x <= 0 ||
108 y <= 0 ||
109 z <= 0)
110 {
111 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
112 }
113 else
114 {
115 LinearAlgebra::AssignVector(voxelDimensions_, x, y, z);
116 Invalidate();
117 }
118 }
119
120
121 const CoordinateSystem3D& VolumeImageGeometry::GetProjectionGeometry(VolumeProjection projection) const
122 {
123 switch (projection)
124 {
125 case VolumeProjection_Axial:
126 return axialGeometry_;
127
128 case VolumeProjection_Coronal:
129 return coronalGeometry_;
130
131 case VolumeProjection_Sagittal:
132 return sagittalGeometry_;
133
134 default:
135 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
136 }
137 }
138
139
140 Vector VolumeImageGeometry::GetVoxelDimensions(VolumeProjection projection) const
141 {
142 switch (projection)
143 {
144 case VolumeProjection_Axial:
145 return voxelDimensions_;
146
147 case VolumeProjection_Coronal:
148 return LinearAlgebra::CreateVector(voxelDimensions_[0], voxelDimensions_[2], voxelDimensions_[1]);
149
150 case VolumeProjection_Sagittal:
151 return LinearAlgebra::CreateVector(voxelDimensions_[1], voxelDimensions_[2], voxelDimensions_[0]);
152
153 default:
154 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
155 }
156 }
157
158
159 unsigned int VolumeImageGeometry::GetProjectionWidth(VolumeProjection projection) const
160 {
161 switch (projection)
162 {
163 case VolumeProjection_Axial:
164 return width_;
165
166 case VolumeProjection_Coronal:
167 return width_;
168
169 case VolumeProjection_Sagittal:
170 return height_;
171
172 default:
173 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
174 }
175 }
176
177
178 unsigned int VolumeImageGeometry::GetProjectionHeight(VolumeProjection projection) const
179 {
180 switch (projection)
181 {
182 case VolumeProjection_Axial:
183 return height_;
184
185 case VolumeProjection_Coronal:
186 return depth_;
187
188 case VolumeProjection_Sagittal:
189 return depth_;
190
191 default:
192 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
193 }
194 }
195
196
197 unsigned int VolumeImageGeometry::GetProjectionDepth(VolumeProjection projection) const
198 {
199 switch (projection)
200 {
201 case VolumeProjection_Axial:
202 return depth_;
203
204 case VolumeProjection_Coronal:
205 return height_;
206
207 case VolumeProjection_Sagittal:
208 return width_;
209
210 default:
211 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
212 }
213 }
214
215
216 Vector VolumeImageGeometry::GetCoordinates(float x,
217 float y,
218 float z) const
219 {
220 Vector p = LinearAlgebra::Product(transform_, LinearAlgebra::CreateVector(x, y, z, 1));
221
222 assert(LinearAlgebra::IsNear(p[3], 1)); // Affine transform, no perspective effect
223
224 // Back to non-homogeneous coordinates
225 return LinearAlgebra::CreateVector(p[0], p[1], p[2]);
226 }
227
228
229 bool VolumeImageGeometry::DetectProjection(VolumeProjection& projection,
230 const Vector& planeNormal) const
231 {
232 if (GeometryToolbox::IsParallel(planeNormal, axialGeometry_.GetNormal()))
233 {
234 projection = VolumeProjection_Axial;
235 return true;
236 }
237 else if (GeometryToolbox::IsParallel(planeNormal, coronalGeometry_.GetNormal()))
238 {
239 projection = VolumeProjection_Coronal;
240 return true;
241 }
242 else if (GeometryToolbox::IsParallel(planeNormal, sagittalGeometry_.GetNormal()))
243 {
244 projection = VolumeProjection_Sagittal;
245 return true;
246 }
247 else
248 {
249 return false;
250 }
251 }
252
253
254 bool VolumeImageGeometry::DetectSlice(VolumeProjection& projection,
255 unsigned int& slice,
256 const CoordinateSystem3D& plane) const
257 {
258 if (!DetectProjection(projection, plane.GetNormal()))
259 {
260 return false;
261 }
262
263 // Transforms the coordinates of the origin of the plane, into the
264 // coordinates of the axial geometry
265 const Vector& origin = plane.GetOrigin();
266 Vector p = LinearAlgebra::Product(
267 transformInverse_,
268 LinearAlgebra::CreateVector(origin[0], origin[1], origin[2], 1));
269
270 assert(LinearAlgebra::IsNear(p[3], 1));
271
272 double z;
273
274 switch (projection)
275 {
276 case VolumeProjection_Axial:
277 z = p[2];
278 break;
279
280 case VolumeProjection_Coronal:
281 z = p[1];
282 break;
283
284 case VolumeProjection_Sagittal:
285 z = p[0];
286 break;
287
288 default:
289 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
290 }
291
292 const unsigned int projectionDepth = GetProjectionDepth(projection);
293
294 z *= static_cast<double>(projectionDepth);
295 if (z < 0)
296 {
297 return false;
298 }
299 else
300 {
301 unsigned int d = static_cast<unsigned int>(std::floor(z));
302 if (d >= projectionDepth)
303 {
304 return false;
305 }
306 else
307 {
308 slice = d;
309 return true;
310 }
311 }
312 }
313
314
315 CoordinateSystem3D VolumeImageGeometry::GetProjectionSlice(VolumeProjection projection,
316 unsigned int z) const
317 {
318 if (z >= GetProjectionDepth(projection))
319 {
320 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
321 }
322
323 Vector dim = GetVoxelDimensions(projection);
324 CoordinateSystem3D plane = GetProjectionGeometry(projection);
325
326 Vector normal = plane.GetNormal();
327 if (projection == VolumeProjection_Sagittal)
328 {
329 /**
330 * WARNING: In sagittal geometry, the normal points to REDUCING
331 * X-axis in the 3D world. This is necessary to keep the
332 * right-hand coordinate system. Hence the negation.
333 **/
334 normal = -normal;
335 }
336
337 plane.SetOrigin(plane.GetOrigin() + static_cast<double>(z) * dim[2] * normal);
338
339 return plane;
340 }
341
342 std::ostream& operator<<(std::ostream& s, const VolumeImageGeometry& v)
343 {
344 s << "width: " << v.width_ << " height: " << v.height_
345 << " depth: " << v.depth_
346 << " axialGeometry: " << v.axialGeometry_
347 << " coronalGeometry: " << v.coronalGeometry_
348 << " sagittalGeometry: " << v.sagittalGeometry_
349 << " voxelDimensions_: " << v.voxelDimensions_
350 << " height: " << v.height_
351 << " transform: " << v.transform_
352 << " transformInverse: " << v.transformInverse_;
353 return s;
354 }
355
356 }