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