Mercurial > hg > orthanc-stone

comparison Framework/Toolbox/CoordinateSystem3D.cpp @ 110:53025eecbc95 wasm

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
renamed SliceGeometry as CoordinateSystem3D
author Sebastien Jodogne <s.jodogne@gmail.com>
date Wed, 14 Jun 2017 15:50:38 +0200
parents Framework/Toolbox/SliceGeometry.cpp@f244018a4e4b
children 2eca030792aa
comparison
equal deleted inserted replaced
109:53bd9277b025 110:53025eecbc95
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 Osimis, 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 "GeometryToolbox.h"
25
26 #include "../../Resources/Orthanc/Core/Logging.h"
27 #include "../../Resources/Orthanc/Core/Toolbox.h"
28 #include "../../Resources/Orthanc/Core/OrthancException.h"
29
30 namespace OrthancStone
31 {
32 void CoordinateSystem3D::CheckAndComputeNormal()
33 {
34 // DICOM expects normal vectors to define the axes: "The row and
35 // column direction cosine vectors shall be normal, i.e., the dot
36 // product of each direction cosine vector with itself shall be
37 // unity."
38 // http://dicom.nema.org/medical/dicom/current/output/chtml/part03/sect_C.7.6.2.html
39 if (!GeometryToolbox::IsNear(boost::numeric::ublas::norm_2(axisX_), 1.0) ||
40 !GeometryToolbox::IsNear(boost::numeric::ublas::norm_2(axisY_), 1.0))
41 {
42 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
43 }
44
45 // The vectors within "Image Orientation Patient" must be
46 // orthogonal, according to the DICOM specification: "The row and
47 // column direction cosine vectors shall be orthogonal, i.e.,
48 // their dot product shall be zero."
49 // http://dicom.nema.org/medical/dicom/current/output/chtml/part03/sect_C.7.6.2.html
50 if (!GeometryToolbox::IsCloseToZero(boost::numeric::ublas::inner_prod(axisX_, axisY_)))
51 {
52 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
53 }
54
55 GeometryToolbox::CrossProduct(normal_, axisX_, axisY_);
56
57 // Just a sanity check, it should be useless by construction
58 assert(GeometryToolbox::IsNear(boost::numeric::ublas::norm_2(normal_), 1.0));
59 }
60
61
62 void CoordinateSystem3D::SetupCanonical()
63 {
64 GeometryToolbox::AssignVector(origin_, 0, 0, 0);
65 GeometryToolbox::AssignVector(axisX_, 1, 0, 0);
66 GeometryToolbox::AssignVector(axisY_, 0, 1, 0);
67 CheckAndComputeNormal();
68 }
69
70
71 CoordinateSystem3D::CoordinateSystem3D(const Vector& origin,
72 const Vector& axisX,
73 const Vector& axisY) :
74 origin_(origin),
75 axisX_(axisX),
76 axisY_(axisY)
77 {
78 CheckAndComputeNormal();
79 }
80
81
82 void CoordinateSystem3D::Setup(const std::string& imagePositionPatient,
83 const std::string& imageOrientationPatient)
84 {
85 std::string tmpPosition = Orthanc::Toolbox::StripSpaces(imagePositionPatient);
86 std::string tmpOrientation = Orthanc::Toolbox::StripSpaces(imageOrientationPatient);
87
88 Vector orientation;
89 if (!GeometryToolbox::ParseVector(origin_, tmpPosition) ||
90 !GeometryToolbox::ParseVector(orientation, tmpOrientation) ||
91 origin_.size() != 3 ||
92 orientation.size() != 6)
93 {
94 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
95 }
96
97 axisX_.resize(3);
98 axisX_[0] = orientation[0];
99 axisX_[1] = orientation[1];
100 axisX_[2] = orientation[2];
101
102 axisY_.resize(3);
103 axisY_[0] = orientation[3];
104 axisY_[1] = orientation[4];
105 axisY_[2] = orientation[5];
106
107 CheckAndComputeNormal();
108 }
109
110
111 CoordinateSystem3D::CoordinateSystem3D(const OrthancPlugins::IDicomDataset& dicom)
112 {
113 std::string a, b;
114
115 if (dicom.GetStringValue(a, OrthancPlugins::DICOM_TAG_IMAGE_POSITION_PATIENT) &&
116 dicom.GetStringValue(b, OrthancPlugins::DICOM_TAG_IMAGE_ORIENTATION_PATIENT))
117 {
118 Setup(a, b);
119 }
120 else
121 {
122 SetupCanonical();
123 }
124 }
125
126
127 Vector CoordinateSystem3D::MapSliceToWorldCoordinates(double x,
128 double y) const
129 {
130 return origin_ + x * axisX_ + y * axisY_;
131 }
132
133
134 double CoordinateSystem3D::ProjectAlongNormal(const Vector& point) const
135 {
136 return boost::numeric::ublas::inner_prod(point, normal_);
137 }
138
139
140 void CoordinateSystem3D::ProjectPoint(double& offsetX,
141 double& offsetY,
142 const Vector& point) const
143 {
144 // Project the point onto the slice
145 Vector projection;
146 GeometryToolbox::ProjectPointOntoPlane(projection, point, normal_, origin_);
147
148 // As the axes are orthonormal vectors thanks to
149 // CheckAndComputeNormal(), the following dot products give the
150 // offset of the origin of the slice wrt. the origin of the
151 // reference plane https://en.wikipedia.org/wiki/Vector_projection
152 offsetX = boost::numeric::ublas::inner_prod(axisX_, projection - origin_);
153 offsetY = boost::numeric::ublas::inner_prod(axisY_, projection - origin_);
154 }
155 }