Mercurial > hg > orthanc-stl

comparison Sources/Plugin.cpp @ 35:ee3bc8f7df5b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
reorganization
author Sebastien Jodogne <s.jodogne@gmail.com>
date Thu, 04 Apr 2024 20:37:08 +0200
parents bee2017f3088
children 13698d34e059
comparison
equal deleted inserted replaced
34:bee2017f3088 35:ee3bc8f7df5b
20 * You should have received a copy of the GNU General Public License 20 * You should have received a copy of the GNU General Public License
21 * along with this program. If not, see <http://www.gnu.org/licenses/>. 21 * along with this program. If not, see <http://www.gnu.org/licenses/>.
22 **/ 22 **/
23 23
24 24
25 #include "StructureSetGeometry.h"
26 #include "StructureSet.h"
25 #include "STLToolbox.h" 27 #include "STLToolbox.h"
26 #include "StructurePolygon.h" 28 #include "StructurePolygon.h"
27 #include "VTKToolbox.h" 29 #include "VTKToolbox.h"
28 #include "Vector3D.h" 30 #include "Vector3D.h"
29 #include "Toolbox.h" 31 #include "Toolbox.h"
152 154
153 155
154 #include <dcmtk/dcmdata/dcfilefo.h> 156 #include <dcmtk/dcmdata/dcfilefo.h>
155 #include <dcmtk/dcmdata/dcsequen.h> 157 #include <dcmtk/dcmdata/dcsequen.h>
156 #include <dcmtk/dcmdata/dcuid.h> 158 #include <dcmtk/dcmdata/dcuid.h>
157
158
159 class StructureSet : public boost::noncopyable
160 {
161 private:
162 std::vector<StructurePolygon*> polygons_;
163 std::string patientId_;
164 std::string studyInstanceUid_;
165 std::string seriesInstanceUid_;
166 std::string sopInstanceUid_;
167 bool hasFrameOfReferenceUid_;
168 std::string frameOfReferenceUid_;
169
170 public:
171 explicit StructureSet(Orthanc::ParsedDicomFile& dicom) :
172 hasFrameOfReferenceUid_(false)
173 {
174 DcmDataset& dataset = *dicom.GetDcmtkObject().getDataset();
175 patientId_ = STLToolbox::GetStringValue(dataset, DCM_PatientID);
176 studyInstanceUid_ = STLToolbox::GetStringValue(dataset, DCM_StudyInstanceUID);
177 seriesInstanceUid_ = STLToolbox::GetStringValue(dataset, DCM_SeriesInstanceUID);
178 sopInstanceUid_ = STLToolbox::GetStringValue(dataset, DCM_SOPInstanceUID);
179
180 DcmSequenceOfItems* frame = NULL;
181 if (!dataset.findAndGetSequence(DCM_ReferencedFrameOfReferenceSequence, frame).good() ||
182 frame == NULL)
183 {
184 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
185 }
186
187 if (frame->card() == 1)
188 {
189 const char* v = NULL;
190 if (frame->getItem(0)->findAndGetString(DCM_FrameOfReferenceUID, v).good() &&
191 v != NULL)
192 {
193 hasFrameOfReferenceUid_ = true;
194 frameOfReferenceUid_.assign(v);
195 }
196 }
197
198 DcmSequenceOfItems* rois = NULL;
199 if (!dataset.findAndGetSequence(DCM_ROIContourSequence, rois).good() ||
200 rois == NULL)
201 {
202 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
203 }
204
205 std::vector<DcmSequenceOfItems*> contours(rois->card());
206 size_t countPolygons = 0;
207
208 for (unsigned long i = 0; i < rois->card(); i++)
209 {
210 DcmSequenceOfItems* contour = NULL;
211 if (!rois->getItem(i)->findAndGetSequence(DCM_ContourSequence, contour).good() ||
212 contour == NULL)
213 {
214 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
215 }
216 else
217 {
218 contours[i] = contour;
219 countPolygons += contour->card();
220 }
221 }
222
223 polygons_.resize(countPolygons);
224
225 size_t pos = 0;
226 for (unsigned long i = 0; i < contours.size(); i++)
227 {
228 for (unsigned long j = 0; j < contours[i]->card(); j++, pos++)
229 {
230 polygons_[pos] = new StructurePolygon(dicom, i, j);
231 }
232 }
233
234 assert(pos == countPolygons);
235 }
236
237 ~StructureSet()
238 {
239 for (size_t i = 0; i < polygons_.size(); i++)
240 {
241 assert(polygons_[i] != NULL);
242 delete polygons_[i];
243 }
244 }
245
246 const std::string& GetPatientId() const
247 {
248 return patientId_;
249 }
250
251 const std::string& GetStudyInstanceUid() const
252 {
253 return studyInstanceUid_;
254 }
255
256 const std::string& GetSeriesInstanceUid() const
257 {
258 return seriesInstanceUid_;
259 }
260
261 const std::string& GetSopInstanceUid() const
262 {
263 return sopInstanceUid_;
264 }
265
266 std::string HashStudy() const
267 {
268 Orthanc::DicomInstanceHasher hasher(patientId_, studyInstanceUid_, seriesInstanceUid_, sopInstanceUid_);
269 return hasher.HashStudy();
270 }
271
272 size_t GetPolygonsCount() const
273 {
274 return polygons_.size();
275 }
276
277 const StructurePolygon& GetPolygon(size_t i) const
278 {
279 if (i >= polygons_.size())
280 {
281 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
282 }
283 else
284 {
285 assert(polygons_[i] != NULL);
286 return *polygons_[i];
287 }
288 }
289
290 bool HasFrameOfReferenceUid() const
291 {
292 return hasFrameOfReferenceUid_;
293 }
294
295 const std::string& GetFrameOfReferenceUid() const
296 {
297 if (hasFrameOfReferenceUid_)
298 {
299 return frameOfReferenceUid_;
300 }
301 else
302 {
303 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
304 }
305 }
306
307 // This static method is faster than constructing the full "StructureSet" object
308 static void ListStructuresNames(std::set<std::string>& target,
309 Orthanc::ParsedDicomFile& source)
310 {
311 target.clear();
312
313 DcmSequenceOfItems* sequence = NULL;
314 if (!source.GetDcmtkObject().getDataset()->findAndGetSequence(DCM_StructureSetROISequence, sequence).good() ||
315 sequence == NULL)
316 {
317 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
318 }
319
320 for (unsigned long i = 0; i < sequence->card(); i++)
321 {
322 DcmItem* item = sequence->getItem(i);
323 if (item == NULL)
324 {
325 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
326 }
327 else
328 {
329 target.insert(STLToolbox::GetStringValue(*item, DCM_ROIName));
330 }
331 }
332 }
333 };
334
335
336 class StructureSetGeometry : public boost::noncopyable
337 {
338 private:
339 bool strict_;
340 Vector3D slicesNormal_;
341 double slicesSpacing_;
342 double minProjectionAlongNormal_;
343 double maxProjectionAlongNormal_;
344 size_t slicesCount_;
345
346 bool LookupProjectionIndex(size_t& index,
347 double z) const
348 {
349 if (slicesCount_ == 0)
350 {
351 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
352 }
353
354 if (z < minProjectionAlongNormal_ ||
355 z > maxProjectionAlongNormal_)
356 {
357 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
358 }
359
360 assert(slicesSpacing_ > 0 &&
361 minProjectionAlongNormal_ < maxProjectionAlongNormal_);
362
363 double d = (z - minProjectionAlongNormal_) / slicesSpacing_;
364
365 if (STLToolbox::IsNear(d, round(d)))
366 {
367 if (d < 0.0 ||
368 d > static_cast<double>(slicesCount_) - 1.0)
369 {
370 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
371 }
372 else
373 {
374 index = static_cast<size_t>(round(d));
375 return true;
376 }
377 }
378 else
379 {
380 return false;
381 }
382 }
383
384 public:
385 StructureSetGeometry(const StructureSet& structures,
386 bool strict) :
387 strict_(strict)
388 {
389 bool isValid = false;
390
391 std::vector<double> projections;
392 projections.reserve(structures.GetPolygonsCount());
393
394 for (size_t i = 0; i < structures.GetPolygonsCount(); i++)
395 {
396 Vector3D normal;
397 if (structures.GetPolygon(i).IsCoplanar(normal))
398 {
399 // Initialize the normal of the whole volume, if need be
400 if (!isValid)
401 {
402 isValid = true;
403 slicesNormal_ = normal;
404 }
405
406 if (Vector3D::AreParallel(normal, slicesNormal_))
407 {
408 // This is a valid slice (it is parallel to the normal)
409 const Vector3D& point = structures.GetPolygon(i).GetPoint(0);
410 projections.push_back(Vector3D::DotProduct(point, slicesNormal_));
411 }
412 else
413 {
414 // RT-STRUCT with non-parallel slices
415 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
416 }
417 }
418 else
419 {
420 // Ignore slices that are not coplanar
421 }
422 }
423
424 if (projections.empty())
425 {
426 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented,
427 "Structure set without a valid geometry");
428 }
429
430 // Only keep unique projections
431
432 std::sort(projections.begin(), projections.end());
433 STLToolbox::RemoveDuplicateValues(projections);
434 assert(!projections.empty());
435
436 if (projections.size() == 1)
437 {
438 // Volume with one single slice
439 minProjectionAlongNormal_ = projections[0];
440 maxProjectionAlongNormal_ = projections[0];
441 slicesSpacing_ = 1; // Arbitrary value
442 slicesCount_ = 1;
443 return;
444 }
445
446
447 // Compute the most probable spacing between the slices
448
449 {
450 std::vector<double> spacings;
451 spacings.resize(projections.size() - 1);
452
453 for (size_t i = 0; i < spacings.size(); i++)
454 {
455 spacings[i] = projections[i + 1] - projections[i];
456 assert(spacings[i] > 0);
457 }
458
459 std::sort(spacings.begin(), spacings.end());
460 STLToolbox::RemoveDuplicateValues(spacings);
461
462 if (spacings.empty())
463 {
464 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
465 }
466
467 slicesSpacing_ = spacings[spacings.size() / 10]; // Take the 90% percentile of smallest spacings
468 assert(slicesSpacing_ > 0);
469 }
470
471
472 // Find the projection along the normal with the largest support
473
474 bool first = true;
475 size_t bestSupport;
476 double bestProjection;
477
478 std::list<size_t> candidates;
479 for (size_t i = 0; i < projections.size(); i++)
480 {
481 candidates.push_back(i);
482 }
483
484 while (!candidates.empty())
485 {
486 std::list<size_t> next;
487
488 size_t countSupport = 0;
489
490 std::list<size_t>::const_iterator it = candidates.begin();
491 size_t reference = *it;
492 it++;
493
494 while (it != candidates.end())
495 {
496 double d = (projections[*it] - projections[reference]) / slicesSpacing_;
497 if (STLToolbox::IsNear(d, round(d)))
498 {
499 countSupport ++;
500 }
501 else
502 {
503 next.push_back(*it);
504 }
505
506 it++;
507 }
508
509 if (first ||
510 countSupport > bestSupport)
511 {
512 first = false;
513 bestSupport = countSupport;
514 bestProjection = projections[reference];
515 }
516
517 if (strict &&
518 !next.empty())
519 {
520 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat,
521 "Structure set with multiple support, which is not allowed in Strict mode");
522 }
523
524 candidates.swap(next);
525 }
526
527
528 // Compute the range of the projections
529
530 minProjectionAlongNormal_ = bestProjection;
531 maxProjectionAlongNormal_ = bestProjection;
532
533 for (size_t i = 0; i < projections.size(); i++)
534 {
535 double d = (projections[i] - bestProjection) / slicesSpacing_;
536 if (STLToolbox::IsNear(d, round(d)))
537 {
538 minProjectionAlongNormal_ = std::min(minProjectionAlongNormal_, projections[i]);
539 maxProjectionAlongNormal_ = std::max(maxProjectionAlongNormal_, projections[i]);
540 }
541 }
542
543 double d = (maxProjectionAlongNormal_ - minProjectionAlongNormal_) / slicesSpacing_;
544 if (STLToolbox::IsNear(d, round(d)))
545 {
546 slicesCount_ = static_cast<size_t>(round(d)) + 1;
547 }
548 else
549 {
550 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
551 }
552
553
554 // Sanity check
555
556 size_t a, b;
557 if (!LookupProjectionIndex(a, minProjectionAlongNormal_) ||
558 !LookupProjectionIndex(b, maxProjectionAlongNormal_) ||
559 a != 0 ||
560 b + 1 != slicesCount_)
561 {
562 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
563 }
564 }
565
566 const Vector3D& GetSlicesNormal() const
567 {
568 return slicesNormal_;
569 }
570
571 double GetSlicesSpacing() const
572 {
573 return slicesSpacing_;
574 }
575
576 double GetMinProjectionAlongNormal() const
577 {
578 return minProjectionAlongNormal_;
579 }
580
581 double GetMaxProjectionAlongNormal() const
582 {
583 return maxProjectionAlongNormal_;
584 }
585
586 bool ProjectAlongNormal(double& z,
587 const StructurePolygon& polygon) const
588 {
589 Vector3D normal;
590 if (polygon.IsCoplanar(normal) &&
591 Vector3D::AreParallel(normal, slicesNormal_))
592 {
593 z = Vector3D::DotProduct(polygon.GetPoint(0), slicesNormal_);
594 return true;
595 }
596 else
597 {
598 return false;
599 }
600 }
601
602 size_t GetSlicesCount() const
603 {
604 return slicesCount_;
605 }
606
607 bool LookupSliceIndex(size_t& slice,
608 const StructurePolygon& polygon) const
609 {
610 double z;
611 return (ProjectAlongNormal(z, polygon) &&
612 LookupProjectionIndex(slice, z));
613 }
614 };
615
616
617 159
618 160
619 class XorFiller : public Orthanc::ImageProcessing::IPolygonFiller 161 class XorFiller : public Orthanc::ImageProcessing::IPolygonFiller
620 { 162 {
621 private: 163 private: