Mercurial > hg > orthanc-stone

comparison Framework/Toolbox/DicomStructureSet.cpp @ 1013:53cc787bd7bc toa2019092301

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
- 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 }