Mercurial > hg > orthanc-stone
comparison Framework/Toolbox/DicomStructureSet.cpp @ 1013:53cc787bd7bc toa2019092301
- Added an optimized ProjectPoint2 to CoordinateSystem3D. It has *not* replaced
the ProjectPoint method because more tests need to be written.
- ProjectPointOntoPlane2 is a faster version of
GeometryToolbox::ProjectPointOntoPlane. Same remark as above.
- DicomStructureSet.cpp now uses this optimized call.
author | Benjamin Golinvaux <bgo@osimis.io> |
---|---|
date | Mon, 23 Sep 2019 15:18:33 +0200 |
parents | 4f28d9459e31 |
children | 709aa65aca17 2d8ab34c8c91 |
comparison
equal
deleted
inserted
replaced
1012:050f01d7951b | 1013:53cc787bd7bc |
---|---|
258 for (Points::const_iterator it = points_.begin(); it != points_.end(); ++it) | 258 for (Points::const_iterator it = points_.begin(); it != points_.end(); ++it) |
259 { | 259 { |
260 if (IsPointOnSliceIfAny(*it)) | 260 if (IsPointOnSliceIfAny(*it)) |
261 { | 261 { |
262 double x, y; | 262 double x, y; |
263 geometry.ProjectPoint(x, y, *it); | 263 geometry.ProjectPoint2(x, y, *it); |
264 extent_.AddPoint(x, y); | 264 extent_.AddPoint(x, y); |
265 } | 265 } |
266 } | 266 } |
267 return true; | 267 return true; |
268 } | 268 } |
299 { | 299 { |
300 return false; | 300 return false; |
301 } | 301 } |
302 | 302 |
303 double x, y; | 303 double x, y; |
304 geometry_.ProjectPoint(x, y, slice.GetOrigin()); | 304 geometry_.ProjectPoint2(x, y, slice.GetOrigin()); |
305 | 305 |
306 bool isOpposite; | 306 bool isOpposite; |
307 if (GeometryToolbox::IsParallelOrOpposite | 307 if (GeometryToolbox::IsParallelOrOpposite |
308 (isOpposite, slice.GetNormal(), geometry_.GetAxisY())) | 308 (isOpposite, slice.GetNormal(), geometry_.GetAxisY())) |
309 { | 309 { |
319 bool isFirst = true; | 319 bool isFirst = true; |
320 double xmin = std::numeric_limits<double>::infinity(); | 320 double xmin = std::numeric_limits<double>::infinity(); |
321 double xmax = -std::numeric_limits<double>::infinity(); | 321 double xmax = -std::numeric_limits<double>::infinity(); |
322 | 322 |
323 double prevX, prevY; | 323 double prevX, prevY; |
324 geometry_.ProjectPoint(prevX, prevY, points_[points_.size() - 1]); | 324 geometry_.ProjectPoint2(prevX, prevY, points_[points_.size() - 1]); |
325 | 325 |
326 for (size_t i = 0; i < points_.size(); i++) | 326 for (size_t i = 0; i < points_.size(); i++) |
327 { | 327 { |
328 // Reference: ../../Resources/Computations/IntersectSegmentAndHorizontalLine.py | 328 // Reference: ../../Resources/Computations/IntersectSegmentAndHorizontalLine.py |
329 double curX, curY; | 329 double curX, curY; |
330 geometry_.ProjectPoint(curX, curY, points_[i]); | 330 geometry_.ProjectPoint2(curX, curY, points_[i]); |
331 | 331 |
332 // if prev* and cur* are on opposite sides of y, this means that the | 332 // if prev* and cur* are on opposite sides of y, this means that the |
333 // segment intersects the plane. | 333 // segment intersects the plane. |
334 if ((prevY < y && curY > y) || | 334 if ((prevY < y && curY > y) || |
335 (prevY > y && curY < y)) | 335 (prevY > y && curY < y)) |
363 sliceThickness_ / 2.0 * geometry_.GetNormal()); | 363 sliceThickness_ / 2.0 * geometry_.GetNormal()); |
364 Vector p2 = (geometry_.MapSliceToWorldCoordinates(xmax, y) - | 364 Vector p2 = (geometry_.MapSliceToWorldCoordinates(xmax, y) - |
365 sliceThickness_ / 2.0 * geometry_.GetNormal()); | 365 sliceThickness_ / 2.0 * geometry_.GetNormal()); |
366 | 366 |
367 // then to the cutting plane geometry... | 367 // then to the cutting plane geometry... |
368 slice.ProjectPoint(x1, y1, p1); | 368 slice.ProjectPoint2(x1, y1, p1); |
369 slice.ProjectPoint(x2, y2, p2); | 369 slice.ProjectPoint2(x2, y2, p2); |
370 return true; | 370 return true; |
371 } | 371 } |
372 } | 372 } |
373 else if (GeometryToolbox::IsParallelOrOpposite | 373 else if (GeometryToolbox::IsParallelOrOpposite |
374 (isOpposite, slice.GetNormal(), geometry_.GetAxisX())) | 374 (isOpposite, slice.GetNormal(), geometry_.GetAxisX())) |
390 bool isFirst = true; | 390 bool isFirst = true; |
391 double ymin = std::numeric_limits<double>::infinity(); | 391 double ymin = std::numeric_limits<double>::infinity(); |
392 double ymax = -std::numeric_limits<double>::infinity(); | 392 double ymax = -std::numeric_limits<double>::infinity(); |
393 | 393 |
394 double prevX, prevY; | 394 double prevX, prevY; |
395 geometry_.ProjectPoint(prevX, prevY, points_[points_.size() - 1]); | 395 geometry_.ProjectPoint2(prevX, prevY, points_[points_.size() - 1]); |
396 | 396 |
397 for (size_t i = 0; i < points_.size(); i++) | 397 for (size_t i = 0; i < points_.size(); i++) |
398 { | 398 { |
399 // Reference: ../../Resources/Computations/IntersectSegmentAndVerticalLine.py | 399 // Reference: ../../Resources/Computations/IntersectSegmentAndVerticalLine.py |
400 double curX, curY; | 400 double curX, curY; |
401 geometry_.ProjectPoint(curX, curY, points_[i]); | 401 geometry_.ProjectPoint2(curX, curY, points_[i]); |
402 | 402 |
403 if ((prevX < x && curX > x) || | 403 if ((prevX < x && curX > x) || |
404 (prevX > x && curX < x)) | 404 (prevX > x && curX < x)) |
405 { | 405 { |
406 double p = (curX * prevY - curY * prevX + x * (curY - prevY)) / (curX - prevX); | 406 double p = (curX * prevY - curY * prevX + x * (curY - prevY)) / (curX - prevX); |
422 Vector p1 = (geometry_.MapSliceToWorldCoordinates(x, ymin) + | 422 Vector p1 = (geometry_.MapSliceToWorldCoordinates(x, ymin) + |
423 sliceThickness_ / 2.0 * geometry_.GetNormal()); | 423 sliceThickness_ / 2.0 * geometry_.GetNormal()); |
424 Vector p2 = (geometry_.MapSliceToWorldCoordinates(x, ymax) - | 424 Vector p2 = (geometry_.MapSliceToWorldCoordinates(x, ymax) - |
425 sliceThickness_ / 2.0 * geometry_.GetNormal()); | 425 sliceThickness_ / 2.0 * geometry_.GetNormal()); |
426 | 426 |
427 slice.ProjectPoint(x1, y1, p1); | 427 slice.ProjectPoint2(x1, y1, p1); |
428 slice.ProjectPoint(x2, y2, p2); | 428 slice.ProjectPoint2(x2, y2, p2); |
429 | 429 |
430 // TODO WHY THIS??? | 430 // TODO WHY THIS??? |
431 y1 = -y1; | 431 y1 = -y1; |
432 y2 = -y2; | 432 y2 = -y2; |
433 | 433 |
836 | 836 |
837 for (Points::const_iterator p = polygon->GetPoints().begin(); | 837 for (Points::const_iterator p = polygon->GetPoints().begin(); |
838 p != polygon->GetPoints().end(); ++p) | 838 p != polygon->GetPoints().end(); ++p) |
839 { | 839 { |
840 double x, y; | 840 double x, y; |
841 slice.ProjectPoint(x, y, *p); | 841 slice.ProjectPoint2(x, y, *p); |
842 polygons.back().push_back(Point2D(x, y)); | 842 polygons.back().push_back(Point2D(x, y)); |
843 } | 843 } |
844 #else | 844 #else |
845 // we need to add all the segments corresponding to this polygon | 845 // we need to add all the segments corresponding to this polygon |
846 const std::vector<Vector>& points3D = polygon->GetPoints(); | 846 const std::vector<Vector>& points3D = polygon->GetPoints(); |
848 { | 848 { |
849 Point2D prev2D; | 849 Point2D prev2D; |
850 { | 850 { |
851 Vector prev = points3D[0]; | 851 Vector prev = points3D[0]; |
852 double prevX, prevY; | 852 double prevX, prevY; |
853 slice.ProjectPoint(prevX, prevY, prev); | 853 slice.ProjectPoint2(prevX, prevY, prev); |
854 prev2D = Point2D(prevX, prevY); | 854 prev2D = Point2D(prevX, prevY); |
855 } | 855 } |
856 | 856 |
857 size_t pointCount = points3D.size(); | 857 size_t pointCount = points3D.size(); |
858 for (size_t ipt = 1; ipt < pointCount; ++ipt) | 858 for (size_t ipt = 1; ipt < pointCount; ++ipt) |
859 { | 859 { |
860 Vector next = points3D[ipt]; | 860 Vector next = points3D[ipt]; |
861 double nextX, nextY; | 861 double nextX, nextY; |
862 slice.ProjectPoint(nextX, nextY, next); | 862 slice.ProjectPoint2(nextX, nextY, next); |
863 Point2D next2D(nextX, nextY); | 863 Point2D next2D(nextX, nextY); |
864 segments.push_back(std::pair<Point2D, Point2D>(prev2D, next2D)); | 864 segments.push_back(std::pair<Point2D, Point2D>(prev2D, next2D)); |
865 prev2D = next2D; | 865 prev2D = next2D; |
866 } | 866 } |
867 } | 867 } |