Mercurial > hg > orthanc-stone

comparison Framework/Toolbox/VolumeImageGeometry.cpp @ 684:7719eb852dd5

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
new class: VolumeImageGeometry
author Sebastien Jodogne <s.jodogne@gmail.com>
date Thu, 16 May 2019 16:47:46 +0200
parents
children 93a8949a1ef7
comparison
equal deleted inserted replaced
683:dbc1d8bfc68a 684:7719eb852dd5
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 Vector VolumeImageGeometry::GetVoxelDimensions(VolumeProjection projection) const
120 {
121 switch (projection)
122 {
123 case VolumeProjection_Axial:
124 return voxelDimensions_;
125
126 case VolumeProjection_Coronal:
127 return LinearAlgebra::CreateVector(voxelDimensions_[0], voxelDimensions_[2], voxelDimensions_[1]);
128
129 case VolumeProjection_Sagittal:
130 return LinearAlgebra::CreateVector(voxelDimensions_[1], voxelDimensions_[2], voxelDimensions_[0]);
131
132 default:
133 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
134 }
135 }
136
137
138 void VolumeImageGeometry::GetSliceSize(unsigned int& width,
139 unsigned int& height,
140 VolumeProjection projection)
141 {
142 switch (projection)
143 {
144 case VolumeProjection_Axial:
145 width = width_;
146 height = height_;
147 break;
148
149 case VolumeProjection_Coronal:
150 width = width_;
151 height = depth_;
152 break;
153
154 case VolumeProjection_Sagittal:
155 width = height_;
156 height = depth_;
157 break;
158
159 default:
160 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
161 }
162 }
163
164
165
166 Vector VolumeImageGeometry::GetCoordinates(float x,
167 float y,
168 float z) const
169 {
170 Vector p = LinearAlgebra::Product(transform_, LinearAlgebra::CreateVector(x, y, z, 1));
171
172 assert(LinearAlgebra::IsNear(p[3], 1)); // Affine transform, no perspective effect
173
174 // Back to non-homogeneous coordinates
175 return LinearAlgebra::CreateVector(p[0], p[1], p[2]);
176 }
177
178
179 bool VolumeImageGeometry::DetectProjection(VolumeProjection& projection,
180 const CoordinateSystem3D& plane) const
181 {
182 if (GeometryToolbox::IsParallel(plane.GetNormal(),
183 axialGeometry_.GetNormal()))
184 {
185 projection = VolumeProjection_Axial;
186 return true;
187 }
188 else if (GeometryToolbox::IsParallel(plane.GetNormal(),
189 coronalGeometry_.GetNormal()))
190 {
191 projection = VolumeProjection_Coronal;
192 return true;
193 }
194 else if (GeometryToolbox::IsParallel(plane.GetNormal(),
195 sagittalGeometry_.GetNormal()))
196 {
197 projection = VolumeProjection_Sagittal;
198 return true;
199 }
200 else
201 {
202 return false;
203 }
204 }
205 }