Mercurial > hg > orthanc-stone

comparison Framework/Toolbox/DicomStructureSet.cpp @ 132:35c2b85836ce wasm

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