Mercurial > hg > orthanc-stone
comparison OrthancStone/Sources/Toolbox/DicomStructureSet.cpp @ 2161:e65fe2e50fde dicom-sr tip
integration mainline->dicom-sr
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Fri, 27 Sep 2024 22:34:17 +0200 |
parents | 917e40af6b45 |
children |
comparison
equal
deleted
inserted
replaced
2152:f68f9a8d0d63 | 2161:e65fe2e50fde |
---|---|
119 #endif | 119 #endif |
120 | 120 |
121 | 121 |
122 namespace OrthancStone | 122 namespace OrthancStone |
123 { | 123 { |
124 static const Orthanc::DicomTag DICOM_TAG_CONTOUR_DATA(0x3006, 0x0050); | |
124 static const Orthanc::DicomTag DICOM_TAG_CONTOUR_GEOMETRIC_TYPE(0x3006, 0x0042); | 125 static const Orthanc::DicomTag DICOM_TAG_CONTOUR_GEOMETRIC_TYPE(0x3006, 0x0042); |
125 static const Orthanc::DicomTag DICOM_TAG_CONTOUR_IMAGE_SEQUENCE(0x3006, 0x0016); | 126 static const Orthanc::DicomTag DICOM_TAG_CONTOUR_IMAGE_SEQUENCE(0x3006, 0x0016); |
126 static const Orthanc::DicomTag DICOM_TAG_CONTOUR_SEQUENCE(0x3006, 0x0040); | 127 static const Orthanc::DicomTag DICOM_TAG_CONTOUR_SEQUENCE(0x3006, 0x0040); |
127 static const Orthanc::DicomTag DICOM_TAG_CONTOUR_DATA(0x3006, 0x0050); | |
128 static const Orthanc::DicomTag DICOM_TAG_NUMBER_OF_CONTOUR_POINTS(0x3006, 0x0046); | 128 static const Orthanc::DicomTag DICOM_TAG_NUMBER_OF_CONTOUR_POINTS(0x3006, 0x0046); |
129 static const Orthanc::DicomTag DICOM_TAG_REFERENCED_ROI_NUMBER(0x3006, 0x0084); | |
129 static const Orthanc::DicomTag DICOM_TAG_REFERENCED_SOP_INSTANCE_UID(0x0008, 0x1155); | 130 static const Orthanc::DicomTag DICOM_TAG_REFERENCED_SOP_INSTANCE_UID(0x0008, 0x1155); |
130 static const Orthanc::DicomTag DICOM_TAG_ROI_CONTOUR_SEQUENCE(0x3006, 0x0039); | 131 static const Orthanc::DicomTag DICOM_TAG_ROI_CONTOUR_SEQUENCE(0x3006, 0x0039); |
131 static const Orthanc::DicomTag DICOM_TAG_ROI_DISPLAY_COLOR(0x3006, 0x002a); | 132 static const Orthanc::DicomTag DICOM_TAG_ROI_DISPLAY_COLOR(0x3006, 0x002a); |
132 static const Orthanc::DicomTag DICOM_TAG_ROI_NAME(0x3006, 0x0026); | 133 static const Orthanc::DicomTag DICOM_TAG_ROI_NAME(0x3006, 0x0026); |
134 static const Orthanc::DicomTag DICOM_TAG_ROI_NUMBER(0x3006, 0x0022); | |
133 static const Orthanc::DicomTag DICOM_TAG_RT_ROI_INTERPRETED_TYPE(0x3006, 0x00a4); | 135 static const Orthanc::DicomTag DICOM_TAG_RT_ROI_INTERPRETED_TYPE(0x3006, 0x00a4); |
134 static const Orthanc::DicomTag DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE(0x3006, 0x0080); | 136 static const Orthanc::DicomTag DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE(0x3006, 0x0080); |
135 static const Orthanc::DicomTag DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE(0x3006, 0x0020); | 137 static const Orthanc::DicomTag DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE(0x3006, 0x0020); |
136 | 138 |
137 | 139 |
279 } | 281 } |
280 } | 282 } |
281 } | 283 } |
282 | 284 |
283 | 285 |
284 static void AddSegmentIfIntersection(Extent2D& extent, | 286 void DicomStructureSet::Polygon::Project(std::list<Extent2D>& target, |
285 const CoordinateSystem3D& cuttingPlane, | |
286 const Vector& p1, | |
287 const Vector& p2, | |
288 double originDistance) | |
289 { | |
290 // Does this segment intersects the cutting plane? | |
291 double d1 = cuttingPlane.ProjectAlongNormal(p1); | |
292 double d2 = cuttingPlane.ProjectAlongNormal(p2); | |
293 | |
294 if ((d1 < originDistance && d2 > originDistance) || | |
295 (d1 > originDistance && d2 < originDistance)) | |
296 { | |
297 // This is an intersection: Add the segment | |
298 double x, y; | |
299 cuttingPlane.ProjectPoint2(x, y, p1); | |
300 extent.AddPoint(x, y); | |
301 cuttingPlane.ProjectPoint2(x, y, p2); | |
302 extent.AddPoint(x, y); | |
303 } | |
304 } | |
305 | |
306 bool DicomStructureSet::Polygon::Project(double& x1, | |
307 double& y1, | |
308 double& x2, | |
309 double& y2, | |
310 const CoordinateSystem3D& cuttingPlane, | 287 const CoordinateSystem3D& cuttingPlane, |
311 const Vector& estimatedNormal, | 288 const Vector& estimatedNormal, |
312 double estimatedSliceThickness) const | 289 double estimatedSliceThickness) const |
313 { | 290 { |
314 if (points_.size() <= 1) | 291 CoordinateSystem3D geometry; |
315 { | |
316 return false; | |
317 } | |
318 | |
319 Vector normal = estimatedNormal; | |
320 double thickness = estimatedSliceThickness; | 292 double thickness = estimatedSliceThickness; |
293 | |
294 /** | |
295 * 1. Estimate the 3D plane associated with this polygon. | |
296 **/ | |
297 | |
321 if (hasSlice_) | 298 if (hasSlice_) |
322 { | 299 { |
323 normal = geometry_.GetNormal(); | 300 // The exact geometry is known for this slice |
301 geometry = geometry_; | |
324 thickness = sliceThickness_; | 302 thickness = sliceThickness_; |
325 } | 303 } |
304 else if (points_.size() < 2) | |
305 { | |
306 return; | |
307 } | |
308 else | |
309 { | |
310 // Estimate the geometry | |
311 Vector origin = points_[0]; | |
312 | |
313 Vector axisX; | |
314 bool found = false; | |
315 for (size_t i = 1; i < points_.size(); i++) | |
316 { | |
317 axisX = points_[1] - origin; | |
318 if (boost::numeric::ublas::norm_2(axisX) > 10.0 * std::numeric_limits<double>::epsilon()) | |
319 { | |
320 found = true; | |
321 break; | |
322 } | |
323 } | |
324 | |
325 if (!found) | |
326 { | |
327 return; // The polygon is too small to extract a reliable geometry out of it | |
328 } | |
329 | |
330 LinearAlgebra::NormalizeVector(axisX); | |
331 | |
332 Vector axisY; | |
333 LinearAlgebra::CrossProduct(axisY, axisX, estimatedNormal); | |
334 geometry = CoordinateSystem3D(origin, axisX, axisY); | |
335 } | |
336 | |
337 | |
338 /** | |
339 * 2. Project the 3D cutting plane as a 2D line onto the polygon plane. | |
340 **/ | |
341 | |
342 double cuttingX1, cuttingY1, cuttingX2, cuttingY2; | |
343 geometry.ProjectPoint(cuttingX1, cuttingY1, cuttingPlane.GetOrigin()); | |
326 | 344 |
327 bool isOpposite; | 345 bool isOpposite; |
328 if (!GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, cuttingPlane.GetAxisX()) && | 346 if (GeometryToolbox::IsParallelOrOpposite(isOpposite, geometry.GetNormal(), cuttingPlane.GetAxisX())) |
329 !GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, cuttingPlane.GetAxisY())) | 347 { |
330 { | 348 geometry.ProjectPoint(cuttingX2, cuttingY2, cuttingPlane.GetOrigin() + cuttingPlane.GetAxisY()); |
331 return false; | 349 } |
332 } | 350 else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, geometry.GetNormal(), cuttingPlane.GetAxisY())) |
333 | 351 { |
334 const double d = cuttingPlane.ProjectAlongNormal(cuttingPlane.GetOrigin()); | 352 geometry.ProjectPoint(cuttingX2, cuttingY2, cuttingPlane.GetOrigin() + cuttingPlane.GetAxisX()); |
335 | |
336 Extent2D extent; | |
337 AddSegmentIfIntersection(extent, cuttingPlane, points_[points_.size() - 1], points_[0], d); | |
338 for (size_t i = 1; i < points_.size(); i++) | |
339 { | |
340 AddSegmentIfIntersection(extent, cuttingPlane, points_[i - 1], points_[i], d); | |
341 } | |
342 | |
343 if (extent.GetWidth() > 0 || | |
344 extent.GetHeight() > 0) | |
345 { | |
346 // Let's convert them to 3D world geometry to add the slice thickness | |
347 Vector p1 = (cuttingPlane.MapSliceToWorldCoordinates(extent.GetX1(), extent.GetY1()) + | |
348 thickness / 2.0 * normal); | |
349 Vector p2 = (cuttingPlane.MapSliceToWorldCoordinates(extent.GetX2(), extent.GetY2()) - | |
350 thickness / 2.0 * normal); | |
351 | |
352 // Then back to the cutting plane geometry | |
353 cuttingPlane.ProjectPoint2(x1, y1, p1); | |
354 cuttingPlane.ProjectPoint2(x2, y2, p2); | |
355 return true; | |
356 } | 353 } |
357 else | 354 else |
358 { | 355 { |
359 return false; | 356 return; |
357 } | |
358 | |
359 | |
360 /** | |
361 * 3. Compute the intersection of the 2D cutting line with the polygon. | |
362 **/ | |
363 | |
364 // Initialize the projection of a point onto a line: | |
365 // https://stackoverflow.com/a/64330724 | |
366 const double abx = cuttingX2 - cuttingX1; | |
367 const double aby = cuttingY2 - cuttingY1; | |
368 const double denominator = abx * abx + aby * aby; | |
369 if (LinearAlgebra::IsCloseToZero(denominator)) | |
370 { | |
371 return; // Should never happen | |
372 } | |
373 | |
374 std::vector<double> intersections; | |
375 intersections.reserve(points_.size()); | |
376 | |
377 for (size_t i = 0; i < points_.size(); i++) | |
378 { | |
379 double segmentX1, segmentY1, segmentX2, segmentY2; | |
380 geometry.ProjectPoint(segmentX1, segmentY1, points_[i]); | |
381 geometry.ProjectPoint(segmentX2, segmentY2, points_[(i + 1) % points_.size()]); | |
382 | |
383 double x, y; | |
384 if (GeometryToolbox::IntersectLineAndSegment(x, y, cuttingX1, cuttingY1, cuttingX2, cuttingY2, | |
385 segmentX1, segmentY1, segmentX2, segmentY2)) | |
386 { | |
387 // For each polygon segment that intersect the cutting line, | |
388 // register its offset over the cutting line | |
389 const double acx = x - cuttingX1; | |
390 const double acy = y - cuttingY1; | |
391 intersections.push_back((abx * acx + aby * acy) / denominator); | |
392 } | |
393 } | |
394 | |
395 | |
396 /** | |
397 * 4. Sort the intersection offsets, then generates one 2D rectangle on the | |
398 * cutting plane from each pair of successive intersections. | |
399 **/ | |
400 | |
401 std::sort(intersections.begin(), intersections.end()); | |
402 | |
403 if (intersections.size() % 2 == 1) | |
404 { | |
405 return; // Should never happen | |
406 } | |
407 | |
408 for (size_t i = 0; i < intersections.size(); i += 2) | |
409 { | |
410 Vector p1, p2; | |
411 | |
412 { | |
413 // Let's convert them to 3D world geometry to add the slice thickness | |
414 const double x1 = cuttingX1 + intersections[i] * abx; | |
415 const double y1 = cuttingY1 + intersections[i] * aby; | |
416 const double x2 = cuttingX1 + intersections[i + 1] * abx; | |
417 const double y2 = cuttingY1 + intersections[i + 1] * aby; | |
418 | |
419 p1 = (geometry.MapSliceToWorldCoordinates(x1, y1) + thickness / 2.0 * geometry.GetNormal()); | |
420 p2 = (geometry.MapSliceToWorldCoordinates(x2, y2) - thickness / 2.0 * geometry.GetNormal()); | |
421 } | |
422 | |
423 { | |
424 // Then back to the cutting plane geometry | |
425 double x1, y1, x2, y2; | |
426 cuttingPlane.ProjectPoint2(x1, y1, p1); | |
427 cuttingPlane.ProjectPoint2(x2, y2, p2); | |
428 | |
429 target.push_back(Extent2D(x1, y1, x2, y2)); | |
430 } | |
360 } | 431 } |
361 } | 432 } |
362 | 433 |
363 | 434 |
364 const DicomStructureSet::Structure& DicomStructureSet::GetStructure(size_t index) const | 435 const DicomStructureSet::Structure& DicomStructureSet::GetStructure(size_t index) const |
386 { | 457 { |
387 #if STONE_TIME_BLOCKING_OPS | 458 #if STONE_TIME_BLOCKING_OPS |
388 boost::posix_time::ptime timerStart = boost::posix_time::microsec_clock::universal_time(); | 459 boost::posix_time::ptime timerStart = boost::posix_time::microsec_clock::universal_time(); |
389 #endif | 460 #endif |
390 | 461 |
462 std::map<int, size_t> roiNumbersIndex; | |
463 | |
391 DicomDatasetReader reader(tags); | 464 DicomDatasetReader reader(tags); |
392 | 465 |
393 size_t count, tmp; | 466 |
394 if (!tags.GetSequenceSize(count, Orthanc::DicomPath(DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE)) || | 467 /** |
395 !tags.GetSequenceSize(tmp, Orthanc::DicomPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE)) || | 468 * 1. Read all the available ROIs. |
396 tmp != count || | 469 **/ |
397 !tags.GetSequenceSize(tmp, Orthanc::DicomPath(DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE)) || | 470 |
398 tmp != count) | 471 { |
399 { | 472 size_t count; |
400 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | 473 if (!tags.GetSequenceSize(count, Orthanc::DicomPath(DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE))) |
401 } | 474 { |
402 | 475 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); |
403 structures_.resize(count); | 476 } |
404 structureNamesIndex_.clear(); | 477 |
405 | 478 structures_.resize(count); |
406 for (size_t i = 0; i < count; i++) | 479 structureNamesIndex_.clear(); |
407 { | 480 |
408 structures_[i].interpretation_ = reader.GetStringValue | 481 for (size_t i = 0; i < count; i++) |
409 (Orthanc::DicomPath(DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE, i, | 482 { |
410 DICOM_TAG_RT_ROI_INTERPRETED_TYPE), | 483 int roiNumber; |
411 "No interpretation"); | 484 if (!reader.GetIntegerValue |
412 | 485 (roiNumber, Orthanc::DicomPath(DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE, i, DICOM_TAG_ROI_NUMBER))) |
413 structures_[i].name_ = reader.GetStringValue | 486 { |
414 (Orthanc::DicomPath(DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE, i, | 487 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); |
415 DICOM_TAG_ROI_NAME), | 488 } |
416 "No name"); | 489 |
417 | 490 if (roiNumbersIndex.find(roiNumber) != roiNumbersIndex.end()) |
418 if (structureNamesIndex_.find(structures_[i].name_) == structureNamesIndex_.end()) | 491 { |
419 { | 492 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat, |
420 structureNamesIndex_[structures_[i].name_] = i; | 493 "Twice the same ROI number: " + boost::lexical_cast<std::string>(roiNumber)); |
421 } | 494 } |
422 else | 495 |
423 { | 496 roiNumbersIndex[roiNumber] = i; |
424 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat, | 497 |
425 "RT-STRUCT with twice the same name for a structure: " + structures_[i].name_); | 498 structures_[i].name_ = reader.GetStringValue |
426 } | 499 (Orthanc::DicomPath(DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE, i, DICOM_TAG_ROI_NAME), "No name"); |
427 | 500 structures_[i].interpretation_ = "No interpretation"; |
428 Vector color; | 501 |
429 if (FastParseVector(color, tags, Orthanc::DicomPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | 502 if (structureNamesIndex_.find(structures_[i].name_) == structureNamesIndex_.end()) |
430 DICOM_TAG_ROI_DISPLAY_COLOR)) && | 503 { |
431 color.size() == 3) | 504 structureNamesIndex_[structures_[i].name_] = i; |
432 { | 505 } |
433 structures_[i].red_ = ConvertColor(color[0]); | 506 else |
434 structures_[i].green_ = ConvertColor(color[1]); | 507 { |
435 structures_[i].blue_ = ConvertColor(color[2]); | 508 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat, |
436 } | 509 "RT-STRUCT with twice the same name for a structure: " + structures_[i].name_); |
437 else | 510 } |
438 { | 511 } |
439 structures_[i].red_ = 255; | 512 } |
440 structures_[i].green_ = 0; | 513 |
441 structures_[i].blue_ = 0; | 514 |
442 } | 515 /** |
443 | 516 * 2. Read the interpretation of the ROIs (if available). |
444 size_t countSlices; | 517 **/ |
445 if (!tags.GetSequenceSize(countSlices, Orthanc::DicomPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | 518 |
446 DICOM_TAG_CONTOUR_SEQUENCE))) | 519 { |
447 { | 520 size_t count; |
448 countSlices = 0; | 521 if (!tags.GetSequenceSize(count, Orthanc::DicomPath(DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE))) |
449 } | 522 { |
450 | 523 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); |
451 LOG(INFO) << "New RT structure: \"" << structures_[i].name_ | 524 } |
452 << "\" with interpretation \"" << structures_[i].interpretation_ | 525 |
453 << "\" containing " << countSlices << " slices (color: " | 526 for (size_t i = 0; i < count; i++) |
454 << static_cast<int>(structures_[i].red_) << "," | 527 { |
455 << static_cast<int>(structures_[i].green_) << "," | 528 std::string interpretation; |
456 << static_cast<int>(structures_[i].blue_) << ")"; | 529 if (reader.GetDataset().GetStringValue(interpretation, |
457 | 530 Orthanc::DicomPath(DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE, i, |
458 /** | 531 DICOM_TAG_RT_ROI_INTERPRETED_TYPE))) |
459 * These temporary variables avoid allocating many vectors in | 532 { |
460 * the loop below (indeed, "Orthanc::DicomPath" handles a | 533 int roiNumber; |
461 * "std::vector<PrefixItem>") | 534 if (!reader.GetIntegerValue(roiNumber, |
462 **/ | 535 Orthanc::DicomPath(DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE, i, |
463 Orthanc::DicomPath countPointsPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | 536 DICOM_TAG_REFERENCED_ROI_NUMBER))) |
464 DICOM_TAG_CONTOUR_SEQUENCE, 0, | 537 { |
465 DICOM_TAG_NUMBER_OF_CONTOUR_POINTS); | 538 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); |
466 | 539 } |
467 Orthanc::DicomPath geometricTypePath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | 540 |
541 std::map<int, size_t>::const_iterator found = roiNumbersIndex.find(roiNumber); | |
542 if (found == roiNumbersIndex.end()) | |
543 { | |
544 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
545 } | |
546 | |
547 structures_[found->second].interpretation_ = interpretation; | |
548 } | |
549 } | |
550 } | |
551 | |
552 | |
553 /** | |
554 * 3. Read the contours. | |
555 **/ | |
556 | |
557 { | |
558 size_t count; | |
559 if (!tags.GetSequenceSize(count, Orthanc::DicomPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE))) | |
560 { | |
561 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
562 } | |
563 | |
564 for (size_t i = 0; i < count; i++) | |
565 { | |
566 int roiNumber; | |
567 if (!reader.GetIntegerValue(roiNumber, | |
568 Orthanc::DicomPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | |
569 DICOM_TAG_REFERENCED_ROI_NUMBER))) | |
570 { | |
571 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
572 } | |
573 | |
574 std::map<int, size_t>::const_iterator found = roiNumbersIndex.find(roiNumber); | |
575 if (found == roiNumbersIndex.end()) | |
576 { | |
577 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
578 } | |
579 | |
580 Structure& target = structures_[found->second]; | |
581 | |
582 Vector color; | |
583 if (FastParseVector(color, tags, Orthanc::DicomPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | |
584 DICOM_TAG_ROI_DISPLAY_COLOR)) && | |
585 color.size() == 3) | |
586 { | |
587 target.red_ = ConvertColor(color[0]); | |
588 target.green_ = ConvertColor(color[1]); | |
589 target.blue_ = ConvertColor(color[2]); | |
590 } | |
591 else | |
592 { | |
593 target.red_ = 255; | |
594 target.green_ = 0; | |
595 target.blue_ = 0; | |
596 } | |
597 | |
598 size_t countSlices; | |
599 if (!tags.GetSequenceSize(countSlices, Orthanc::DicomPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | |
600 DICOM_TAG_CONTOUR_SEQUENCE))) | |
601 { | |
602 countSlices = 0; | |
603 } | |
604 | |
605 LOG(INFO) << "New RT structure: \"" << target.name_ | |
606 << "\" with interpretation \"" << target.interpretation_ | |
607 << "\" containing " << countSlices << " slices (color: " | |
608 << static_cast<int>(target.red_) << "," | |
609 << static_cast<int>(target.green_) << "," | |
610 << static_cast<int>(target.blue_) << ")"; | |
611 | |
612 /** | |
613 * These temporary variables avoid allocating many vectors in | |
614 * the loop below (indeed, "Orthanc::DicomPath" handles a | |
615 * "std::vector<PrefixItem>") | |
616 **/ | |
617 Orthanc::DicomPath countPointsPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | |
468 DICOM_TAG_CONTOUR_SEQUENCE, 0, | 618 DICOM_TAG_CONTOUR_SEQUENCE, 0, |
469 DICOM_TAG_CONTOUR_GEOMETRIC_TYPE); | 619 DICOM_TAG_NUMBER_OF_CONTOUR_POINTS); |
470 | 620 |
471 Orthanc::DicomPath imageSequencePath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | 621 Orthanc::DicomPath geometricTypePath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, |
622 DICOM_TAG_CONTOUR_SEQUENCE, 0, | |
623 DICOM_TAG_CONTOUR_GEOMETRIC_TYPE); | |
624 | |
625 Orthanc::DicomPath imageSequencePath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | |
626 DICOM_TAG_CONTOUR_SEQUENCE, 0, | |
627 DICOM_TAG_CONTOUR_IMAGE_SEQUENCE); | |
628 | |
629 // (3006,0039)[i] / (0x3006, 0x0040)[0] / (0x3006, 0x0016)[0] / (0x0008, 0x1155) | |
630 Orthanc::DicomPath referencedInstancePath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | |
631 DICOM_TAG_CONTOUR_SEQUENCE, 0, | |
632 DICOM_TAG_CONTOUR_IMAGE_SEQUENCE, 0, | |
633 DICOM_TAG_REFERENCED_SOP_INSTANCE_UID); | |
634 | |
635 Orthanc::DicomPath contourDataPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | |
472 DICOM_TAG_CONTOUR_SEQUENCE, 0, | 636 DICOM_TAG_CONTOUR_SEQUENCE, 0, |
473 DICOM_TAG_CONTOUR_IMAGE_SEQUENCE); | 637 DICOM_TAG_CONTOUR_DATA); |
474 | 638 |
475 // (3006,0039)[i] / (0x3006, 0x0040)[0] / (0x3006, 0x0016)[0] / (0x0008, 0x1155) | 639 for (size_t j = 0; j < countSlices; j++) |
476 Orthanc::DicomPath referencedInstancePath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | 640 { |
477 DICOM_TAG_CONTOUR_SEQUENCE, 0, | 641 unsigned int countPoints; |
478 DICOM_TAG_CONTOUR_IMAGE_SEQUENCE, 0, | 642 |
479 DICOM_TAG_REFERENCED_SOP_INSTANCE_UID); | 643 countPointsPath.SetPrefixIndex(1, j); |
480 | 644 if (!reader.GetUnsignedIntegerValue(countPoints, countPointsPath)) |
481 Orthanc::DicomPath contourDataPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, | 645 { |
482 DICOM_TAG_CONTOUR_SEQUENCE, 0, | 646 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); |
483 DICOM_TAG_CONTOUR_DATA); | 647 } |
484 | 648 |
485 for (size_t j = 0; j < countSlices; j++) | 649 //LOG(INFO) << "Parsing slice containing " << countPoints << " vertices"; |
486 { | 650 |
487 unsigned int countPoints; | 651 geometricTypePath.SetPrefixIndex(1, j); |
488 | 652 std::string type = reader.GetMandatoryStringValue(geometricTypePath); |
489 countPointsPath.SetPrefixIndex(1, j); | 653 if (type != "CLOSED_PLANAR") |
490 if (!reader.GetUnsignedIntegerValue(countPoints, countPointsPath)) | 654 { |
491 { | 655 LOG(WARNING) << "Ignoring contour with geometry type: " << type; |
492 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | 656 continue; |
493 } | 657 } |
494 | 658 |
495 //LOG(INFO) << "Parsing slice containing " << countPoints << " vertices"; | 659 size_t size; |
496 | 660 |
497 geometricTypePath.SetPrefixIndex(1, j); | 661 imageSequencePath.SetPrefixIndex(1, j); |
498 std::string type = reader.GetMandatoryStringValue(geometricTypePath); | 662 if (!tags.GetSequenceSize(size, imageSequencePath) || size != 1) |
499 if (type != "CLOSED_PLANAR") | 663 { |
500 { | 664 LOG(ERROR) << "The ContourImageSequence sequence (tag 3006,0016) must be present and contain one entry."; |
501 LOG(WARNING) << "Ignoring contour with geometry type: " << type; | 665 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); |
502 continue; | 666 } |
503 } | 667 |
504 | 668 referencedInstancePath.SetPrefixIndex(1, j); |
505 size_t size; | 669 std::string sopInstanceUid = reader.GetMandatoryStringValue(referencedInstancePath); |
506 | 670 |
507 imageSequencePath.SetPrefixIndex(1, j); | 671 contourDataPath.SetPrefixIndex(1, j); |
508 if (!tags.GetSequenceSize(size, imageSequencePath) || size != 1) | 672 std::string slicesData = reader.GetMandatoryStringValue(contourDataPath); |
509 { | 673 |
510 LOG(ERROR) << "The ContourImageSequence sequence (tag 3006,0016) must be present and contain one entry."; | 674 Vector points; |
511 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); | 675 |
512 } | 676 if (!GenericToolbox::FastParseVector(points, slicesData) || |
513 | 677 points.size() != 3 * countPoints) |
514 referencedInstancePath.SetPrefixIndex(1, j); | 678 { |
515 std::string sopInstanceUid = reader.GetMandatoryStringValue(referencedInstancePath); | 679 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); |
516 | 680 } |
517 contourDataPath.SetPrefixIndex(1, j); | 681 |
518 std::string slicesData = reader.GetMandatoryStringValue(contourDataPath); | 682 // seen in real world |
519 | 683 if(Orthanc::Toolbox::StripSpaces(sopInstanceUid) == "") |
520 Vector points; | 684 { |
521 | 685 LOG(ERROR) << "WARNING. The following Dicom tag (Referenced SOP Instance UID) contains an empty value : // (3006,0039)[" << i << "] / (0x3006, 0x0040)[0] / (0x3006, 0x0016)[0] / (0x0008, 0x1155)"; |
522 if (!GenericToolbox::FastParseVector(points, slicesData) || | 686 } |
523 points.size() != 3 * countPoints) | 687 |
524 { | 688 Polygon polygon(sopInstanceUid); |
525 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | 689 polygon.Reserve(countPoints); |
526 } | 690 |
527 | 691 for (size_t k = 0; k < countPoints; k++) |
528 // seen in real world | 692 { |
529 if(Orthanc::Toolbox::StripSpaces(sopInstanceUid) == "") | 693 Vector v(3); |
530 { | 694 v[0] = points[3 * k]; |
531 LOG(ERROR) << "WARNING. The following Dicom tag (Referenced SOP Instance UID) contains an empty value : // (3006,0039)[" << i << "] / (0x3006, 0x0040)[0] / (0x3006, 0x0016)[0] / (0x0008, 0x1155)"; | 695 v[1] = points[3 * k + 1]; |
532 } | 696 v[2] = points[3 * k + 2]; |
533 | 697 polygon.AddPoint(v); |
534 Polygon polygon(sopInstanceUid); | 698 } |
535 polygon.Reserve(countPoints); | 699 |
536 | 700 target.polygons_.push_back(polygon); |
537 for (size_t k = 0; k < countPoints; k++) | 701 } |
538 { | |
539 Vector v(3); | |
540 v[0] = points[3 * k]; | |
541 v[1] = points[3 * k + 1]; | |
542 v[2] = points[3 * k + 2]; | |
543 polygon.AddPoint(v); | |
544 } | |
545 | |
546 structures_[i].polygons_.push_back(polygon); | |
547 } | 702 } |
548 } | 703 } |
549 | 704 |
550 EstimateGeometry(); | 705 EstimateGeometry(); |
551 | 706 |
802 std::vector<BoostPolygon> projected; | 957 std::vector<BoostPolygon> projected; |
803 | 958 |
804 for (Polygons::const_iterator polygon = structure.polygons_.begin(); | 959 for (Polygons::const_iterator polygon = structure.polygons_.begin(); |
805 polygon != structure.polygons_.end(); ++polygon) | 960 polygon != structure.polygons_.end(); ++polygon) |
806 { | 961 { |
807 double x1, y1, x2, y2; | 962 std::list<Extent2D> rectangles; |
808 | 963 polygon->Project(rectangles, cutting, GetEstimatedNormal(), GetEstimatedSliceThickness()); |
809 if (polygon->Project(x1, y1, x2, y2, cutting, GetEstimatedNormal(), GetEstimatedSliceThickness())) | 964 |
810 { | 965 for (std::list<Extent2D>::const_iterator it = rectangles.begin(); it != rectangles.end(); ++it) |
811 projected.push_back(CreateRectangle(x1, y1, x2, y2)); | 966 { |
967 projected.push_back(CreateRectangle(it->GetX1(), it->GetY1(), it->GetX2(), it->GetY2())); | |
812 } | 968 } |
813 } | 969 } |
814 | 970 |
815 BoostMultiPolygon merged; | 971 BoostMultiPolygon merged; |
816 Union(merged, projected); | 972 Union(merged, projected); |
832 std::list<Extent2D> rectangles; | 988 std::list<Extent2D> rectangles; |
833 | 989 |
834 for (Polygons::const_iterator polygon = structure.polygons_.begin(); | 990 for (Polygons::const_iterator polygon = structure.polygons_.begin(); |
835 polygon != structure.polygons_.end(); ++polygon) | 991 polygon != structure.polygons_.end(); ++polygon) |
836 { | 992 { |
837 double x1, y1, x2, y2; | 993 polygon->Project(rectangles, cutting, GetEstimatedNormal(), GetEstimatedSliceThickness()); |
838 | |
839 if (polygon->Project(x1, y1, x2, y2, cutting, GetEstimatedNormal(), GetEstimatedSliceThickness())) | |
840 { | |
841 rectangles.push_back(Extent2D(x1, y1, x2, y2)); | |
842 } | |
843 } | 994 } |
844 | 995 |
845 typedef std::list< std::vector<ScenePoint2D> > Contours; | 996 typedef std::list< std::vector<ScenePoint2D> > Contours; |
846 | 997 |
847 Contours contours; | 998 Contours contours; |