Mercurial > hg > orthanc-stone

comparison Framework/Toolbox/OrientedBoundingBox.cpp @ 140:2115530d3703 wasm

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