Mercurial > hg > orthanc-stone
comparison Framework/Toolbox/OrientedBoundingBox.cpp @ 140:2115530d3703 wasm
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 } |