Mercurial > hg > orthanc-stone

comparison OrthancStone/Sources/Toolbox/CoordinateSystem3D.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/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 }