Mercurial > hg > orthanc-stone
comparison Framework/Volumes/OrientedVolumeBoundingBox.cpp @ 742:fa5febe0f0c2
moved OrientedBoundingBox in the Volumes folder
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Wed, 22 May 2019 09:41:03 +0200 |
parents | Framework/Toolbox/OrientedBoundingBox.cpp@c3bbb130abc4 |
children | 2d8ab34c8c91 |
comparison
equal
deleted
inserted
replaced
741:c1d6a566dfd3 | 742:fa5febe0f0c2 |
---|---|
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-2019 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 <Core/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 } |