0
|
1 /**
|
|
2 * Stone of Orthanc
|
|
3 * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics
|
|
4 * Department, University Hospital of Liege, Belgium
|
|
5 *
|
|
6 * This program is free software: you can redistribute it and/or
|
|
7 * modify it under the terms of the GNU General Public License as
|
|
8 * published by the Free Software Foundation, either version 3 of the
|
|
9 * License, or (at your option) any later version.
|
|
10 *
|
|
11 * In addition, as a special exception, the copyright holders of this
|
|
12 * program give permission to link the code of its release with the
|
|
13 * OpenSSL project's "OpenSSL" library (or with modified versions of it
|
|
14 * that use the same license as the "OpenSSL" library), and distribute
|
|
15 * the linked executables. You must obey the GNU General Public License
|
|
16 * in all respects for all of the code used other than "OpenSSL". If you
|
|
17 * modify file(s) with this exception, you may extend this exception to
|
|
18 * your version of the file(s), but you are not obligated to do so. If
|
|
19 * you do not wish to do so, delete this exception statement from your
|
|
20 * version. If you delete this exception statement from all source files
|
|
21 * in the program, then also delete it here.
|
|
22 *
|
|
23 * This program is distributed in the hope that it will be useful, but
|
|
24 * WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
25 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
26 * General Public License for more details.
|
|
27 *
|
|
28 * You should have received a copy of the GNU General Public License
|
|
29 * along with this program. If not, see <http://www.gnu.org/licenses/>.
|
|
30 **/
|
|
31
|
|
32
|
|
33 #include "SliceGeometry.h"
|
|
34
|
|
35 #include "GeometryToolbox.h"
|
|
36
|
|
37 #include "../Orthanc/Core/Logging.h"
|
|
38 #include "../Orthanc/Core/Toolbox.h"
|
|
39 #include "../Orthanc/Core/OrthancException.h"
|
|
40
|
|
41 namespace OrthancStone
|
|
42 {
|
|
43 void SliceGeometry::CheckAndComputeNormal()
|
|
44 {
|
|
45 // DICOM expects normal vectors to define the axes: "The row and
|
|
46 // column direction cosine vectors shall be normal, i.e., the dot
|
|
47 // product of each direction cosine vector with itself shall be
|
|
48 // unity."
|
|
49 // http://dicom.nema.org/medical/dicom/current/output/chtml/part03/sect_C.7.6.2.html
|
|
50 if (!GeometryToolbox::IsNear(boost::numeric::ublas::norm_2(axisX_), 1.0) ||
|
|
51 !GeometryToolbox::IsNear(boost::numeric::ublas::norm_2(axisY_), 1.0))
|
|
52 {
|
|
53 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
|
|
54 }
|
|
55
|
|
56 // The vectors within "Image Orientation Patient" must be
|
|
57 // orthogonal, according to the DICOM specification: "The row and
|
|
58 // column direction cosine vectors shall be orthogonal, i.e.,
|
|
59 // their dot product shall be zero."
|
|
60 // http://dicom.nema.org/medical/dicom/current/output/chtml/part03/sect_C.7.6.2.html
|
|
61 if (!GeometryToolbox::IsCloseToZero(boost::numeric::ublas::inner_prod(axisX_, axisY_)))
|
|
62 {
|
|
63 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
|
|
64 }
|
|
65
|
|
66 GeometryToolbox::CrossProduct(normal_, axisX_, axisY_);
|
|
67
|
|
68 // Just a sanity check, it should be useless by construction
|
|
69 assert(GeometryToolbox::IsNear(boost::numeric::ublas::norm_2(normal_), 1.0));
|
|
70 }
|
|
71
|
|
72
|
|
73 void SliceGeometry::SetupCanonical()
|
|
74 {
|
|
75 GeometryToolbox::AssignVector(origin_, 0, 0, 0);
|
|
76 GeometryToolbox::AssignVector(axisX_, 1, 0, 0);
|
|
77 GeometryToolbox::AssignVector(axisY_, 0, 1, 0);
|
|
78 CheckAndComputeNormal();
|
|
79 }
|
|
80
|
|
81
|
|
82 SliceGeometry::SliceGeometry(const Vector& origin,
|
|
83 const Vector& axisX,
|
|
84 const Vector& axisY) :
|
|
85 origin_(origin),
|
|
86 axisX_(axisX),
|
|
87 axisY_(axisY)
|
|
88 {
|
|
89 CheckAndComputeNormal();
|
|
90 }
|
|
91
|
|
92
|
|
93 void SliceGeometry::Setup(const std::string& imagePositionPatient,
|
|
94 const std::string& imageOrientationPatient)
|
|
95 {
|
|
96 std::string tmpPosition = Orthanc::Toolbox::StripSpaces(imagePositionPatient);
|
|
97 std::string tmpOrientation = Orthanc::Toolbox::StripSpaces(imageOrientationPatient);
|
|
98
|
|
99 Vector orientation;
|
|
100 if (!GeometryToolbox::ParseVector(origin_, tmpPosition) ||
|
|
101 !GeometryToolbox::ParseVector(orientation, tmpOrientation) ||
|
|
102 origin_.size() != 3 ||
|
|
103 orientation.size() != 6)
|
|
104 {
|
|
105 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
|
|
106 }
|
|
107
|
|
108 axisX_.resize(3);
|
|
109 axisX_[0] = orientation[0];
|
|
110 axisX_[1] = orientation[1];
|
|
111 axisX_[2] = orientation[2];
|
|
112
|
|
113 axisY_.resize(3);
|
|
114 axisY_[0] = orientation[3];
|
|
115 axisY_[1] = orientation[4];
|
|
116 axisY_[2] = orientation[5];
|
|
117
|
|
118 CheckAndComputeNormal();
|
|
119 }
|
|
120
|
|
121
|
|
122 SliceGeometry::SliceGeometry(const DicomDataset& dicom)
|
|
123 {
|
|
124 if (dicom.HasTag(DICOM_TAG_IMAGE_POSITION_PATIENT) &&
|
|
125 dicom.HasTag(DICOM_TAG_IMAGE_ORIENTATION_PATIENT))
|
|
126 {
|
|
127 Setup(dicom.GetStringValue(DICOM_TAG_IMAGE_POSITION_PATIENT),
|
|
128 dicom.GetStringValue(DICOM_TAG_IMAGE_ORIENTATION_PATIENT));
|
|
129 }
|
|
130 else
|
|
131 {
|
|
132 SetupCanonical();
|
|
133 }
|
|
134 }
|
|
135
|
|
136
|
|
137 Vector SliceGeometry::MapSliceToWorldCoordinates(double x,
|
|
138 double y) const
|
|
139 {
|
|
140 return origin_ + x * axisX_ + y * axisY_;
|
|
141 }
|
|
142
|
|
143
|
|
144 double SliceGeometry::ProjectAlongNormal(const Vector& point) const
|
|
145 {
|
|
146 return boost::numeric::ublas::inner_prod(point, normal_);
|
|
147 }
|
|
148
|
|
149
|
|
150 void SliceGeometry::ProjectPoint(double& offsetX,
|
|
151 double& offsetY,
|
|
152 const Vector& point) const
|
|
153 {
|
|
154 // Project the point onto the slice
|
|
155 Vector projection;
|
|
156 GeometryToolbox::ProjectPointOntoPlane(projection, point, normal_, origin_);
|
|
157
|
|
158 // As the axes are orthonormal vectors thanks to
|
|
159 // CheckAndComputeNormal(), the following dot products give the
|
|
160 // offset of the origin of the slice wrt. the origin of the
|
|
161 // reference plane https://en.wikipedia.org/wiki/Vector_projection
|
|
162 offsetX = boost::numeric::ublas::inner_prod(axisX_, projection - origin_);
|
|
163 offsetY = boost::numeric::ublas::inner_prod(axisY_, projection - origin_);
|
|
164 }
|
|
165 }
|