Mercurial > hg > orthanc-stone
annotate Framework/Toolbox/OrientedBoundingBox.cpp @ 146:fb7d602e7025 wasm
OrientedBoundingBox::ComputeExtent
| author | Sebastien Jodogne <s.jodogne@gmail.com> |
|---|---|
| date | Mon, 22 Jan 2018 17:30:34 +0100 |
| parents | 2115530d3703 |
| children | c5044bbfc303 |
| 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 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 } | |
|
146
fb7d602e7025
OrientedBoundingBox::ComputeExtent
Sebastien Jodogne <s.jodogne@gmail.com>
parents:
140
diff
changeset
|
267 |
|
fb7d602e7025
OrientedBoundingBox::ComputeExtent
Sebastien Jodogne <s.jodogne@gmail.com>
parents:
140
diff
changeset
|
268 |
|
fb7d602e7025
OrientedBoundingBox::ComputeExtent
Sebastien Jodogne <s.jodogne@gmail.com>
parents:
140
diff
changeset
|
269 bool OrientedBoundingBox::ComputeExtent(Extent2D& extent, |
|
fb7d602e7025
OrientedBoundingBox::ComputeExtent
Sebastien Jodogne <s.jodogne@gmail.com>
parents:
140
diff
changeset
|
270 const CoordinateSystem3D& plane) const |
|
fb7d602e7025
OrientedBoundingBox::ComputeExtent
Sebastien Jodogne <s.jodogne@gmail.com>
parents:
140
diff
changeset
|
271 { |
|
fb7d602e7025
OrientedBoundingBox::ComputeExtent
Sebastien Jodogne <s.jodogne@gmail.com>
parents:
140
diff
changeset
|
272 extent.Reset(); |
|
fb7d602e7025
OrientedBoundingBox::ComputeExtent
Sebastien Jodogne <s.jodogne@gmail.com>
parents:
140
diff
changeset
|
273 |
|
fb7d602e7025
OrientedBoundingBox::ComputeExtent
Sebastien Jodogne <s.jodogne@gmail.com>
parents:
140
diff
changeset
|
274 std::vector<Vector> points; |
|
fb7d602e7025
OrientedBoundingBox::ComputeExtent
Sebastien Jodogne <s.jodogne@gmail.com>
parents:
140
diff
changeset
|
275 if (HasIntersection(points, plane)) |
|
fb7d602e7025
OrientedBoundingBox::ComputeExtent
Sebastien Jodogne <s.jodogne@gmail.com>
parents:
140
diff
changeset
|
276 { |
|
fb7d602e7025
OrientedBoundingBox::ComputeExtent
Sebastien Jodogne <s.jodogne@gmail.com>
parents:
140
diff
changeset
|
277 for (size_t i = 0; i < points.size(); i++) |
|
fb7d602e7025
OrientedBoundingBox::ComputeExtent
Sebastien Jodogne <s.jodogne@gmail.com>
parents:
140
diff
changeset
|
278 { |
|
fb7d602e7025
OrientedBoundingBox::ComputeExtent
Sebastien Jodogne <s.jodogne@gmail.com>
parents:
140
diff
changeset
|
279 double x, y; |
|
fb7d602e7025
OrientedBoundingBox::ComputeExtent
Sebastien Jodogne <s.jodogne@gmail.com>
parents:
140
diff
changeset
|
280 plane.ProjectPoint(x, y, points[i]); |
|
fb7d602e7025
OrientedBoundingBox::ComputeExtent
Sebastien Jodogne <s.jodogne@gmail.com>
parents:
140
diff
changeset
|
281 extent.AddPoint(x, y); |
|
fb7d602e7025
OrientedBoundingBox::ComputeExtent
Sebastien Jodogne <s.jodogne@gmail.com>
parents:
140
diff
changeset
|
282 } |
|
fb7d602e7025
OrientedBoundingBox::ComputeExtent
Sebastien Jodogne <s.jodogne@gmail.com>
parents:
140
diff
changeset
|
283 |
|
fb7d602e7025
OrientedBoundingBox::ComputeExtent
Sebastien Jodogne <s.jodogne@gmail.com>
parents:
140
diff
changeset
|
284 return true; |
|
fb7d602e7025
OrientedBoundingBox::ComputeExtent
Sebastien Jodogne <s.jodogne@gmail.com>
parents:
140
diff
changeset
|
285 } |
|
fb7d602e7025
OrientedBoundingBox::ComputeExtent
Sebastien Jodogne <s.jodogne@gmail.com>
parents:
140
diff
changeset
|
286 else |
|
fb7d602e7025
OrientedBoundingBox::ComputeExtent
Sebastien Jodogne <s.jodogne@gmail.com>
parents:
140
diff
changeset
|
287 { |
|
fb7d602e7025
OrientedBoundingBox::ComputeExtent
Sebastien Jodogne <s.jodogne@gmail.com>
parents:
140
diff
changeset
|
288 return false; |
|
fb7d602e7025
OrientedBoundingBox::ComputeExtent
Sebastien Jodogne <s.jodogne@gmail.com>
parents:
140
diff
changeset
|
289 } |
|
fb7d602e7025
OrientedBoundingBox::ComputeExtent
Sebastien Jodogne <s.jodogne@gmail.com>
parents:
140
diff
changeset
|
290 } |
| 140 | 291 } |