Mercurial > hg > orthanc-stone

comparison OrthancStone/Sources/Volumes/OrientedVolumeBoundingBox.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/Volumes/OrientedVolumeBoundingBox.cpp@30deba7bc8e2
children 4fb8fdf03314
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 "OrientedVolumeBoundingBox.h"
23
24 #include "../Toolbox/GeometryToolbox.h"
25 #include "ImageBuffer3D.h"
26
27 #include <OrthancException.h>
28
29 #include <cassert>
30
31 namespace OrthancStone
32 {
33 OrientedVolumeBoundingBox::OrientedVolumeBoundingBox(const VolumeImageGeometry& geometry)
34 {
35 unsigned int n = geometry.GetDepth();
36 if (n < 1)
37 {
38 throw Orthanc::OrthancException(Orthanc::ErrorCode_IncompatibleImageSize);
39 }
40
41 Vector dim = geometry.GetVoxelDimensions(VolumeProjection_Axial);
42
43 u_ = geometry.GetAxialGeometry().GetAxisX();
44 v_ = geometry.GetAxialGeometry().GetAxisY();
45 w_ = geometry.GetAxialGeometry().GetNormal();
46
47 hu_ = static_cast<double>(geometry.GetWidth() * dim[0] / 2.0);
48 hv_ = static_cast<double>(geometry.GetHeight() * dim[1] / 2.0);
49 hw_ = static_cast<double>(geometry.GetDepth() * dim[2] / 2.0);
50
51 c_ = (geometry.GetAxialGeometry().GetOrigin() +
52 (hu_ - dim[0] / 2.0) * u_ +
53 (hv_ - dim[1] / 2.0) * v_ +
54 (hw_ - dim[2] / 2.0) * w_);
55 }
56
57
58 bool OrientedVolumeBoundingBox::HasIntersectionWithPlane(std::vector<Vector>& points,
59 const Vector& normal,
60 double d) const
61 {
62 assert(normal.size() == 3);
63
64 double r = (hu_ * fabs(boost::numeric::ublas::inner_prod(normal, u_)) +
65 hv_ * fabs(boost::numeric::ublas::inner_prod(normal, v_)) +
66 hw_ * fabs(boost::numeric::ublas::inner_prod(normal, w_)));
67
68 double s = boost::numeric::ublas::inner_prod(normal, c_) + d;
69
70 if (fabs(s) >= r)
71 {
72 // No intersection, or intersection is reduced to a single point
73 return false;
74 }
75 else
76 {
77 Vector p;
78
79 // Loop over all the 12 edges (segments) of the oriented
80 // bounding box, and check whether they intersect the plane
81
82 // X-aligned edges
83 if (GeometryToolbox::IntersectPlaneAndSegment
84 (p, normal, d,
85 c_ - u_ * hu_ - v_ * hv_ - w_ * hw_,
86 c_ + u_ * hu_ - v_ * hv_ - w_ * hw_))
87 {
88 points.push_back(p);
89 }
90
91 if (GeometryToolbox::IntersectPlaneAndSegment
92 (p, normal, d,
93 c_ - u_ * hu_ + v_ * hv_ - w_ * hw_,
94 c_ + u_ * hu_ + v_ * hv_ - w_ * hw_))
95 {
96 points.push_back(p);
97 }
98
99 if (GeometryToolbox::IntersectPlaneAndSegment
100 (p, normal, d,
101 c_ - u_ * hu_ - v_ * hv_ + w_ * hw_,
102 c_ + u_ * hu_ - v_ * hv_ + w_ * hw_))
103 {
104 points.push_back(p);
105 }
106
107 if (GeometryToolbox::IntersectPlaneAndSegment
108 (p, normal, d,
109 c_ - u_ * hu_ + v_ * hv_ + w_ * hw_,
110 c_ + u_ * hu_ + v_ * hv_ + w_ * hw_))
111 {
112 points.push_back(p);
113 }
114
115 // Y-aligned edges
116 if (GeometryToolbox::IntersectPlaneAndSegment
117 (p, normal, d,
118 c_ - u_ * hu_ - v_ * hv_ - w_ * hw_,
119 c_ - u_ * hu_ + v_ * hv_ - w_ * hw_))
120 {
121 points.push_back(p);
122 }
123
124 if (GeometryToolbox::IntersectPlaneAndSegment
125 (p, normal, d,
126 c_ + u_ * hu_ - v_ * hv_ - w_ * hw_,
127 c_ + u_ * hu_ + v_ * hv_ - w_ * hw_))
128 {
129 points.push_back(p);
130 }
131
132 if (GeometryToolbox::IntersectPlaneAndSegment
133 (p, normal, d,
134 c_ - u_ * hu_ - v_ * hv_ + w_ * hw_,
135 c_ - u_ * hu_ + v_ * hv_ + w_ * hw_))
136 {
137 points.push_back(p);
138 }
139
140 if (GeometryToolbox::IntersectPlaneAndSegment
141 (p, normal, d,
142 c_ + u_ * hu_ - v_ * hv_ + w_ * hw_,
143 c_ + u_ * hu_ + v_ * hv_ + w_ * hw_))
144 {
145 points.push_back(p);
146 }
147
148 // Z-aligned edges
149 if (GeometryToolbox::IntersectPlaneAndSegment
150 (p, normal, d,
151 c_ - u_ * hu_ - v_ * hv_ - w_ * hw_,
152 c_ - u_ * hu_ - v_ * hv_ + w_ * hw_))
153 {
154 points.push_back(p);
155 }
156
157 if (GeometryToolbox::IntersectPlaneAndSegment
158 (p, normal, d,
159 c_ + u_ * hu_ - v_ * hv_ - w_ * hw_,
160 c_ + u_ * hu_ - v_ * hv_ + w_ * hw_))
161 {
162 points.push_back(p);
163 }
164
165 if (GeometryToolbox::IntersectPlaneAndSegment
166 (p, normal, d,
167 c_ - u_ * hu_ + v_ * hv_ - w_ * hw_,
168 c_ - u_ * hu_ + v_ * hv_ + w_ * hw_))
169 {
170 points.push_back(p);
171 }
172
173 if (GeometryToolbox::IntersectPlaneAndSegment
174 (p, normal, d,
175 c_ + u_ * hu_ + v_ * hv_ - w_ * hw_,
176 c_ + u_ * hu_ + v_ * hv_ + w_ * hw_))
177 {
178 points.push_back(p);
179 }
180
181 return true;
182 }
183 }
184
185
186 bool OrientedVolumeBoundingBox::HasIntersection(std::vector<Vector>& points,
187 const CoordinateSystem3D& plane) const
188 {
189 // From the vector equation of a 3D plane (specified by origin
190 // and normal), to the general equation of a 3D plane (which
191 // looses information about the origin of the coordinate system)
192 const Vector& normal = plane.GetNormal();
193 const Vector& origin = plane.GetOrigin();
194 double d = -(normal[0] * origin[0] + normal[1] * origin[1] + normal[2] * origin[2]);
195
196 return HasIntersectionWithPlane(points, normal, d);
197 }
198
199
200 bool OrientedVolumeBoundingBox::Contains(const Vector& p) const
201 {
202 assert(p.size() == 3);
203
204 const Vector q = p - c_;
205
206 return (fabs(boost::numeric::ublas::inner_prod(q, u_)) <= hu_ &&
207 fabs(boost::numeric::ublas::inner_prod(q, v_)) <= hv_ &&
208 fabs(boost::numeric::ublas::inner_prod(q, w_)) <= hw_);
209 }
210
211
212 void OrientedVolumeBoundingBox::FromInternalCoordinates(Vector& target,
213 double x,
214 double y,
215 double z) const
216 {
217 target = (c_ +
218 u_ * 2.0 * hu_ * (x - 0.5) +
219 v_ * 2.0 * hv_ * (y - 0.5) +
220 w_ * 2.0 * hw_ * (z - 0.5));
221 }
222
223
224 void OrientedVolumeBoundingBox::FromInternalCoordinates(Vector& target,
225 const Vector& source) const
226 {
227 assert(source.size() == 3);
228 FromInternalCoordinates(target, source[0], source[1], source[2]);
229 }
230
231
232 void OrientedVolumeBoundingBox::ToInternalCoordinates(Vector& target,
233 const Vector& source) const
234 {
235 assert(source.size() == 3);
236 const Vector q = source - c_;
237
238 double x = boost::numeric::ublas::inner_prod(q, u_) / (2.0 * hu_) + 0.5;
239 double y = boost::numeric::ublas::inner_prod(q, v_) / (2.0 * hv_) + 0.5;
240 double z = boost::numeric::ublas::inner_prod(q, w_) / (2.0 * hw_) + 0.5;
241
242 LinearAlgebra::AssignVector(target, x, y, z);
243 }
244
245
246 bool OrientedVolumeBoundingBox::ComputeExtent(Extent2D& extent,
247 const CoordinateSystem3D& plane) const
248 {
249 extent.Reset();
250
251 std::vector<Vector> points;
252 if (HasIntersection(points, plane))
253 {
254 for (size_t i = 0; i < points.size(); i++)
255 {
256 double x, y;
257 plane.ProjectPoint(x, y, points[i]);
258 extent.AddPoint(x, y);
259 }
260
261 return true;
262 }
263 else
264 {
265 return false;
266 }
267 }
268 }