Mercurial > hg > orthanc-stone
comparison Framework/Toolbox/DicomStructureSet.cpp @ 132:35c2b85836ce wasm
fix rtstruct projections
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Fri, 17 Nov 2017 18:01:31 +0100 |
parents | 3e6163a53b16 |
children | e2fe9352f240 |
comparison
equal
deleted
inserted
replaced
131:3e6163a53b16 | 132:35c2b85836ce |
---|---|
211 double& y1, | 211 double& y1, |
212 double& x2, | 212 double& x2, |
213 double& y2, | 213 double& y2, |
214 const CoordinateSystem3D& slice) const | 214 const CoordinateSystem3D& slice) const |
215 { | 215 { |
216 // TODO Optimize this method using a sweep-line algorithm for polygons | |
217 | |
216 if (!hasSlice_ || | 218 if (!hasSlice_ || |
217 points_.empty()) | 219 points_.size() <= 1) |
218 { | 220 { |
219 return false; | 221 return false; |
220 } | 222 } |
221 | 223 |
222 double x, y; | 224 double x, y; |
224 | 226 |
225 bool isOpposite; | 227 bool isOpposite; |
226 if (GeometryToolbox::IsParallelOrOpposite | 228 if (GeometryToolbox::IsParallelOrOpposite |
227 (isOpposite, slice.GetNormal(), geometry_.GetAxisY())) | 229 (isOpposite, slice.GetNormal(), geometry_.GetAxisY())) |
228 { | 230 { |
229 if (y >= extent_.GetY1() && | 231 if (y < extent_.GetY1() || |
230 y <= extent_.GetY2()) | 232 y > extent_.GetY2()) |
231 { | 233 { |
232 Vector p1 = (geometry_.MapSliceToWorldCoordinates(extent_.GetX1(), y) + | 234 // The polygon does not intersect the input slice |
235 return false; | |
236 } | |
237 | |
238 bool isFirst = true; | |
239 double xmin = std::numeric_limits<double>::infinity(); | |
240 double xmax = -std::numeric_limits<double>::infinity(); | |
241 | |
242 double prevX, prevY; | |
243 geometry_.ProjectPoint(prevX, prevY, points_[points_.size() - 1]); | |
244 | |
245 for (size_t i = 0; i < points_.size(); i++) | |
246 { | |
247 // Reference: ../../Resources/Computations/IntersectSegmentAndHorizontalLine.py | |
248 double curX, curY; | |
249 geometry_.ProjectPoint(curX, curY, points_[i]); | |
250 | |
251 if ((prevY < y && curY > y) || | |
252 (prevY > y && curY < y)) | |
253 { | |
254 double p = (curX * prevY - curY * prevX + y * (prevX - curX)) / (prevY - curY); | |
255 xmin = std::min(xmin, p); | |
256 xmax = std::max(xmax, p); | |
257 isFirst = false; | |
258 } | |
259 | |
260 prevX = curX; | |
261 prevY = curY; | |
262 } | |
263 | |
264 if (isFirst) | |
265 { | |
266 return false; | |
267 } | |
268 else | |
269 { | |
270 Vector p1 = (geometry_.MapSliceToWorldCoordinates(xmin, y) + | |
233 sliceThickness_ / 2.0 * geometry_.GetNormal()); | 271 sliceThickness_ / 2.0 * geometry_.GetNormal()); |
234 Vector p2 = (geometry_.MapSliceToWorldCoordinates(extent_.GetX2(), y) - | 272 Vector p2 = (geometry_.MapSliceToWorldCoordinates(xmax, y) - |
235 sliceThickness_ / 2.0 * geometry_.GetNormal()); | 273 sliceThickness_ / 2.0 * geometry_.GetNormal()); |
236 | 274 |
237 slice.ProjectPoint(x1, y1, p1); | 275 slice.ProjectPoint(x1, y1, p1); |
238 slice.ProjectPoint(x2, y2, p2); | 276 slice.ProjectPoint(x2, y2, p2); |
239 return true; | 277 return true; |
240 } | 278 } |
241 else | |
242 { | |
243 return false; | |
244 } | |
245 } | 279 } |
246 else if (GeometryToolbox::IsParallelOrOpposite | 280 else if (GeometryToolbox::IsParallelOrOpposite |
247 (isOpposite, slice.GetNormal(), geometry_.GetAxisX())) | 281 (isOpposite, slice.GetNormal(), geometry_.GetAxisX())) |
248 { | 282 { |
249 if (x >= extent_.GetX1() && | 283 if (x < extent_.GetX1() || |
250 x <= extent_.GetX2()) | 284 x > extent_.GetX2()) |
251 { | 285 { |
252 Vector p1 = (geometry_.MapSliceToWorldCoordinates(x, extent_.GetY1()) + | 286 return false; |
287 } | |
288 | |
289 bool isFirst = true; | |
290 double ymin = std::numeric_limits<double>::infinity(); | |
291 double ymax = -std::numeric_limits<double>::infinity(); | |
292 | |
293 double prevX, prevY; | |
294 geometry_.ProjectPoint(prevX, prevY, points_[points_.size() - 1]); | |
295 | |
296 for (size_t i = 0; i < points_.size(); i++) | |
297 { | |
298 // Reference: ../../Resources/Computations/IntersectSegmentAndVerticalLine.py | |
299 double curX, curY; | |
300 geometry_.ProjectPoint(curX, curY, points_[i]); | |
301 | |
302 if ((prevX < x && curX > x) || | |
303 (prevX > x && curX < x)) | |
304 { | |
305 double p = (curX * prevY - curY * prevX + x * (curY - prevY)) / (curX - prevX); | |
306 ymin = std::min(ymin, p); | |
307 ymax = std::max(ymax, p); | |
308 isFirst = false; | |
309 } | |
310 | |
311 prevX = curX; | |
312 prevY = curY; | |
313 } | |
314 | |
315 if (isFirst) | |
316 { | |
317 return false; | |
318 } | |
319 else | |
320 { | |
321 Vector p1 = (geometry_.MapSliceToWorldCoordinates(x, ymin) + | |
253 sliceThickness_ / 2.0 * geometry_.GetNormal()); | 322 sliceThickness_ / 2.0 * geometry_.GetNormal()); |
254 Vector p2 = (geometry_.MapSliceToWorldCoordinates(x, extent_.GetY2()) - | 323 Vector p2 = (geometry_.MapSliceToWorldCoordinates(x, ymax) - |
255 sliceThickness_ / 2.0 * geometry_.GetNormal()); | 324 sliceThickness_ / 2.0 * geometry_.GetNormal()); |
256 | 325 |
257 slice.ProjectPoint(x1, y1, p1); | 326 slice.ProjectPoint(x1, y1, p1); |
258 slice.ProjectPoint(x2, y2, p2); | 327 slice.ProjectPoint(x2, y2, p2); |
328 | |
329 // TODO WHY THIS??? | |
330 y1 = -y1; | |
331 y2 = -y2; | |
332 | |
259 return true; | 333 return true; |
260 } | |
261 else | |
262 { | |
263 return false; | |
264 } | 334 } |
265 } | 335 } |
266 else | 336 else |
267 { | 337 { |
268 // Should not happen | 338 // Should not happen |
402 { | 472 { |
403 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | 473 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); |
404 } | 474 } |
405 | 475 |
406 Polygon polygon(sopInstanceUid); | 476 Polygon polygon(sopInstanceUid); |
477 polygon.Reserve(countPoints); | |
407 | 478 |
408 for (size_t k = 0; k < countPoints; k++) | 479 for (size_t k = 0; k < countPoints; k++) |
409 { | 480 { |
410 Vector v(3); | 481 Vector v(3); |
411 v[0] = points[3 * k]; | 482 v[0] = points[3 * k]; |
669 return true; | 740 return true; |
670 } | 741 } |
671 else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisX()) || | 742 else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisX()) || |
672 GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisY())) | 743 GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisY())) |
673 { | 744 { |
745 #if 1 | |
674 // Sagittal or coronal projection | 746 // Sagittal or coronal projection |
675 std::vector<BoostPolygon> projected; | 747 std::vector<BoostPolygon> projected; |
676 | 748 |
677 for (Polygons::iterator polygon = structure.polygons_.begin(); | 749 for (Polygons::iterator polygon = structure.polygons_.begin(); |
678 polygon != structure.polygons_.end(); ++polygon) | 750 polygon != structure.polygons_.end(); ++polygon) |
696 for (size_t j = 0; j < outer.size(); j++) | 768 for (size_t j = 0; j < outer.size(); j++) |
697 { | 769 { |
698 polygons[i][j] = std::make_pair(outer[j].x(), outer[j].y()); | 770 polygons[i][j] = std::make_pair(outer[j].x(), outer[j].y()); |
699 } | 771 } |
700 } | 772 } |
773 #else | |
774 for (Polygons::iterator polygon = structure.polygons_.begin(); | |
775 polygon != structure.polygons_.end(); ++polygon) | |
776 { | |
777 double x1, y1, x2, y2; | |
778 if (polygon->Project(x1, y1, x2, y2, slice)) | |
779 { | |
780 std::vector<PolygonPoint> p(4); | |
781 p[0] = std::make_pair(x1, y1); | |
782 p[1] = std::make_pair(x2, y1); | |
783 p[2] = std::make_pair(x2, y2); | |
784 p[3] = std::make_pair(x1, y2); | |
785 polygons.push_back(p); | |
786 } | |
787 } | |
788 #endif | |
701 | 789 |
702 return true; | 790 return true; |
703 } | 791 } |
704 else | 792 else |
705 { | 793 { |