Mercurial > hg > orthanc-stone
comparison OrthancStone/Sources/Toolbox/CoordinateSystem3D.cpp @ 1512:244ad1e4e76a
reorganization of folders
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Tue, 07 Jul 2020 16:21:02 +0200 |
parents | Framework/Toolbox/CoordinateSystem3D.cpp@d8af188ab545 |
children | 85e117739eca |
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 "CoordinateSystem3D.h" | |
23 | |
24 #include "LinearAlgebra.h" | |
25 #include "GeometryToolbox.h" | |
26 | |
27 #include <Logging.h> | |
28 #include <Toolbox.h> | |
29 #include <OrthancException.h> | |
30 | |
31 namespace OrthancStone | |
32 { | |
33 void CoordinateSystem3D::CheckAndComputeNormal() | |
34 { | |
35 // DICOM expects normal vectors to define the axes: "The row and | |
36 // column direction cosine vectors shall be normal, i.e., the dot | |
37 // product of each direction cosine vector with itself shall be | |
38 // unity." | |
39 // http://dicom.nema.org/medical/dicom/current/output/chtml/part03/sect_C.7.6.2.html | |
40 if (!LinearAlgebra::IsNear(boost::numeric::ublas::norm_2(axisX_), 1.0) || | |
41 !LinearAlgebra::IsNear(boost::numeric::ublas::norm_2(axisY_), 1.0)) | |
42 { | |
43 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
44 } | |
45 | |
46 // The vectors within "Image Orientation Patient" must be | |
47 // orthogonal, according to the DICOM specification: "The row and | |
48 // column direction cosine vectors shall be orthogonal, i.e., | |
49 // their dot product shall be zero." | |
50 // http://dicom.nema.org/medical/dicom/current/output/chtml/part03/sect_C.7.6.2.html | |
51 if (!LinearAlgebra::IsCloseToZero(boost::numeric::ublas::inner_prod(axisX_, axisY_))) | |
52 { | |
53 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
54 } | |
55 | |
56 LinearAlgebra::CrossProduct(normal_, axisX_, axisY_); | |
57 | |
58 d_ = -(normal_[0] * origin_[0] + normal_[1] * origin_[1] + normal_[2] * origin_[2]); | |
59 | |
60 // Just a sanity check, it should be useless by construction | |
61 assert(LinearAlgebra::IsNear(boost::numeric::ublas::norm_2(normal_), 1.0)); | |
62 } | |
63 | |
64 | |
65 void CoordinateSystem3D::SetupCanonical() | |
66 { | |
67 LinearAlgebra::AssignVector(origin_, 0, 0, 0); | |
68 LinearAlgebra::AssignVector(axisX_, 1, 0, 0); | |
69 LinearAlgebra::AssignVector(axisY_, 0, 1, 0); | |
70 CheckAndComputeNormal(); | |
71 } | |
72 | |
73 | |
74 CoordinateSystem3D::CoordinateSystem3D(const Vector& origin, | |
75 const Vector& axisX, | |
76 const Vector& axisY) : | |
77 origin_(origin), | |
78 axisX_(axisX), | |
79 axisY_(axisY) | |
80 { | |
81 CheckAndComputeNormal(); | |
82 } | |
83 | |
84 | |
85 void CoordinateSystem3D::Setup(const std::string& imagePositionPatient, | |
86 const std::string& imageOrientationPatient) | |
87 { | |
88 std::string tmpPosition = Orthanc::Toolbox::StripSpaces(imagePositionPatient); | |
89 std::string tmpOrientation = Orthanc::Toolbox::StripSpaces(imageOrientationPatient); | |
90 | |
91 Vector orientation; | |
92 if (!LinearAlgebra::ParseVector(origin_, tmpPosition) || | |
93 !LinearAlgebra::ParseVector(orientation, tmpOrientation) || | |
94 origin_.size() != 3 || | |
95 orientation.size() != 6) | |
96 { | |
97 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
98 } | |
99 | |
100 axisX_.resize(3); | |
101 axisX_[0] = orientation[0]; | |
102 axisX_[1] = orientation[1]; | |
103 axisX_[2] = orientation[2]; | |
104 | |
105 axisY_.resize(3); | |
106 axisY_[0] = orientation[3]; | |
107 axisY_[1] = orientation[4]; | |
108 axisY_[2] = orientation[5]; | |
109 | |
110 CheckAndComputeNormal(); | |
111 } | |
112 | |
113 | |
114 CoordinateSystem3D::CoordinateSystem3D(const IDicomDataset& dicom) | |
115 { | |
116 std::string a, b; | |
117 | |
118 if (dicom.GetStringValue(a, Orthanc::DICOM_TAG_IMAGE_POSITION_PATIENT) && | |
119 dicom.GetStringValue(b, Orthanc::DICOM_TAG_IMAGE_ORIENTATION_PATIENT)) | |
120 { | |
121 Setup(a, b); | |
122 } | |
123 else | |
124 { | |
125 SetupCanonical(); | |
126 } | |
127 } | |
128 | |
129 | |
130 CoordinateSystem3D::CoordinateSystem3D(const Orthanc::DicomMap& dicom) | |
131 { | |
132 std::string a, b; | |
133 | |
134 if (dicom.LookupStringValue(a, Orthanc::DICOM_TAG_IMAGE_POSITION_PATIENT, false) && | |
135 dicom.LookupStringValue(b, Orthanc::DICOM_TAG_IMAGE_ORIENTATION_PATIENT, false)) | |
136 { | |
137 Setup(a, b); | |
138 } | |
139 else | |
140 { | |
141 SetupCanonical(); | |
142 } | |
143 } | |
144 | |
145 | |
146 void CoordinateSystem3D::SetOrigin(const Vector& origin) | |
147 { | |
148 if (origin.size() != 3) | |
149 { | |
150 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); | |
151 } | |
152 else | |
153 { | |
154 origin_ = origin; | |
155 } | |
156 } | |
157 | |
158 | |
159 Vector CoordinateSystem3D::MapSliceToWorldCoordinates(double x, | |
160 double y) const | |
161 { | |
162 return origin_ + x * axisX_ + y * axisY_; | |
163 } | |
164 | |
165 | |
166 double CoordinateSystem3D::ProjectAlongNormal(const Vector& point) const | |
167 { | |
168 return boost::numeric::ublas::inner_prod(point, normal_); | |
169 } | |
170 | |
171 void CoordinateSystem3D::ProjectPoint2(double& offsetX, double& offsetY, const Vector& point) const | |
172 { | |
173 // Project the point onto the slice | |
174 double projectionX,projectionY,projectionZ; | |
175 GeometryToolbox::ProjectPointOntoPlane2(projectionX, projectionY, projectionZ, point, normal_, origin_); | |
176 | |
177 // As the axes are orthonormal vectors thanks to | |
178 // CheckAndComputeNormal(), the following dot products give the | |
179 // offset of the origin of the slice wrt. the origin of the | |
180 // reference plane https://en.wikipedia.org/wiki/Vector_projection | |
181 offsetX = axisX_[0] * (projectionX - origin_[0]) + axisX_[1] * (projectionY - origin_[1]) + axisX_[2] * (projectionZ - origin_[2]); | |
182 offsetY = axisY_[0] * (projectionX - origin_[0]) + axisY_[1] * (projectionY - origin_[1]) + axisY_[2] * (projectionZ - origin_[2]); | |
183 } | |
184 | |
185 void CoordinateSystem3D::ProjectPoint(double& offsetX, | |
186 double& offsetY, | |
187 const Vector& point) const | |
188 { | |
189 // Project the point onto the slice | |
190 Vector projection; | |
191 GeometryToolbox::ProjectPointOntoPlane(projection, point, normal_, origin_); | |
192 | |
193 // As the axes are orthonormal vectors thanks to | |
194 // CheckAndComputeNormal(), the following dot products give the | |
195 // offset of the origin of the slice wrt. the origin of the | |
196 // reference plane https://en.wikipedia.org/wiki/Vector_projection | |
197 offsetX = boost::numeric::ublas::inner_prod(axisX_, projection - origin_); | |
198 offsetY = boost::numeric::ublas::inner_prod(axisY_, projection - origin_); | |
199 } | |
200 | |
201 bool CoordinateSystem3D::IntersectSegment(Vector& p, | |
202 const Vector& edgeFrom, | |
203 const Vector& edgeTo) const | |
204 { | |
205 return GeometryToolbox::IntersectPlaneAndSegment(p, normal_, d_, edgeFrom, edgeTo); | |
206 } | |
207 | |
208 | |
209 bool CoordinateSystem3D::IntersectLine(Vector& p, | |
210 const Vector& origin, | |
211 const Vector& direction) const | |
212 { | |
213 return GeometryToolbox::IntersectPlaneAndLine(p, normal_, d_, origin, direction); | |
214 } | |
215 | |
216 | |
217 bool CoordinateSystem3D::ComputeDistance(double& distance, | |
218 const CoordinateSystem3D& a, | |
219 const CoordinateSystem3D& b) | |
220 { | |
221 bool opposite = false; // Ignored | |
222 | |
223 if (OrthancStone::GeometryToolbox::IsParallelOrOpposite( | |
224 opposite, a.GetNormal(), b.GetNormal())) | |
225 { | |
226 distance = std::abs(a.ProjectAlongNormal(a.GetOrigin()) - | |
227 a.ProjectAlongNormal(b.GetOrigin())); | |
228 return true; | |
229 } | |
230 else | |
231 { | |
232 return false; | |
233 } | |
234 } | |
235 | |
236 std::ostream& operator<< (std::ostream& s, const CoordinateSystem3D& that) | |
237 { | |
238 s << "origin: " << that.origin_ << " normal: " << that.normal_ | |
239 << " axisX: " << that.axisX_ << " axisY: " << that.axisY_ | |
240 << " D: " << that.d_; | |
241 return s; | |
242 } | |
243 | |
244 | |
245 CoordinateSystem3D CoordinateSystem3D::NormalizeCuttingPlane(const CoordinateSystem3D& plane) | |
246 { | |
247 double ox, oy; | |
248 plane.ProjectPoint(ox, oy, LinearAlgebra::CreateVector(0, 0, 0)); | |
249 | |
250 CoordinateSystem3D normalized(plane); | |
251 normalized.SetOrigin(plane.MapSliceToWorldCoordinates(ox, oy)); | |
252 return normalized; | |
253 } | |
254 } |