Mercurial > hg > orthanc-stl

comparison Sources/Plugin.cpp @ 33:2460b376d3f7

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 18:50:11 +0200
parents 976da5476810
children bee2017f3088
comparison
equal deleted inserted replaced
32:976da5476810 33:2460b376d3f7
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 "StructurePolygon.h"
25 #include "VTKToolbox.h" 26 #include "VTKToolbox.h"
26 #include "Vector3D.h" 27 #include "Vector3D.h"
27 #include "Toolbox.h" 28 #include "Toolbox.h"
28 #include "Extent2D.h" 29 #include "Extent2D.h"
29 30
38 #include <Logging.h> 39 #include <Logging.h>
39 #include <OrthancFramework.h> 40 #include <OrthancFramework.h>
40 #include <SerializationToolbox.h> 41 #include <SerializationToolbox.h>
41 #include <SystemToolbox.h> 42 #include <SystemToolbox.h>
42 43
44 #include <vtkNew.h>
45
43 #include <boost/thread/shared_mutex.hpp> 46 #include <boost/thread/shared_mutex.hpp>
44 47
45 #define ORTHANC_PLUGIN_NAME "stl" 48 #define ORTHANC_PLUGIN_NAME "stl"
46 49
47 50
145 } 148 }
146 149
147 150
148 151
149 152
150 #include <dcmtk/dcmdata/dcdeftag.h>
151 #include <dcmtk/dcmdata/dcfilefo.h> 153 #include <dcmtk/dcmdata/dcfilefo.h>
152 #include <dcmtk/dcmdata/dcitem.h>
153 #include <dcmtk/dcmdata/dcsequen.h> 154 #include <dcmtk/dcmdata/dcsequen.h>
154 #include <dcmtk/dcmdata/dcuid.h> 155 #include <dcmtk/dcmdata/dcuid.h>
155
156
157 static std::string GetStringValue(DcmItem& item,
158 const DcmTagKey& key)
159 {
160 const char* s = NULL;
161 if (!item.findAndGetString(key, s).good() ||
162 s == NULL)
163 {
164 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
165 }
166 else
167 {
168 return Orthanc::Toolbox::StripSpaces(s);
169 }
170 }
171
172
173 static void ListStructuresNames(std::set<std::string>& target,
174 Orthanc::ParsedDicomFile& source)
175 {
176 target.clear();
177
178 DcmSequenceOfItems* sequence = NULL;
179 if (!source.GetDcmtkObject().getDataset()->findAndGetSequence(DCM_StructureSetROISequence, sequence).good() ||
180 sequence == NULL)
181 {
182 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
183 }
184
185 for (unsigned long i = 0; i < sequence->card(); i++)
186 {
187 DcmItem* item = sequence->getItem(i);
188 if (item == NULL)
189 {
190 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
191 }
192 else
193 {
194 target.insert(GetStringValue(*item, DCM_ROIName));
195 }
196 }
197 }
198
199
200 class StructurePolygon : public boost::noncopyable
201 {
202 private:
203 std::string roiName_;
204 std::string referencedSopInstanceUid_;
205 uint8_t red_;
206 uint8_t green_;
207 uint8_t blue_;
208 std::vector<Vector3D> points_;
209
210 public:
211 StructurePolygon(Orthanc::ParsedDicomFile& dicom,
212 unsigned long roiIndex,
213 unsigned long contourIndex)
214 {
215 DcmDataset& dataset = *dicom.GetDcmtkObject().getDataset();
216
217 DcmItem* structure = NULL;
218 DcmItem* roi = NULL;
219 DcmItem* contour = NULL;
220 DcmSequenceOfItems* referenced = NULL;
221
222 if (!dataset.findAndGetSequenceItem(DCM_StructureSetROISequence, structure, roiIndex).good() ||
223 structure == NULL ||
224 !dataset.findAndGetSequenceItem(DCM_ROIContourSequence, roi, roiIndex).good() ||
225 roi == NULL ||
226 !roi->findAndGetSequenceItem(DCM_ContourSequence, contour, contourIndex).good() ||
227 contour == NULL ||
228 !contour->findAndGetSequence(DCM_ContourImageSequence, referenced).good() ||
229 referenced == NULL ||
230 referenced->card() != 1)
231 {
232 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
233 }
234
235 roiName_ = GetStringValue(*structure, DCM_ROIName);
236 referencedSopInstanceUid_ = GetStringValue(*referenced->getItem(0), DCM_ReferencedSOPInstanceUID);
237
238 if (GetStringValue(*contour, DCM_ContourGeometricType) != "CLOSED_PLANAR")
239 {
240 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
241 }
242
243 {
244 std::vector<std::string> tokens;
245 Orthanc::Toolbox::TokenizeString(tokens, GetStringValue(*roi, DCM_ROIDisplayColor), '\\');
246
247 uint32_t r, g, b;
248 if (tokens.size() != 3 ||
249 !Orthanc::SerializationToolbox::ParseFirstUnsignedInteger32(r, tokens[0]) ||
250 !Orthanc::SerializationToolbox::ParseFirstUnsignedInteger32(g, tokens[1]) ||
251 !Orthanc::SerializationToolbox::ParseFirstUnsignedInteger32(b, tokens[2]) ||
252 r > 255 ||
253 g > 255 ||
254 b > 255)
255 {
256 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
257 }
258
259 red_ = r;
260 green_ = g;
261 blue_ = b;
262 }
263
264 {
265 std::vector<std::string> tokens;
266 Orthanc::Toolbox::TokenizeString(tokens, GetStringValue(*contour, DCM_ContourData), '\\');
267
268 const std::string s = GetStringValue(*contour, DCM_NumberOfContourPoints);
269
270 uint32_t countPoints;
271 if (!Orthanc::SerializationToolbox::ParseUnsignedInteger32(countPoints, s) ||
272 tokens.size() != 3 * countPoints)
273 {
274 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
275 }
276
277 points_.reserve(countPoints);
278
279 for (size_t i = 0; i < tokens.size(); i += 3)
280 {
281 double x, y, z;
282 if (!Toolbox::MyParseDouble(x, tokens[i]) ||
283 !Toolbox::MyParseDouble(y, tokens[i + 1]) ||
284 !Toolbox::MyParseDouble(z, tokens[i + 2]))
285 {
286 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
287 }
288
289 points_.push_back(Vector3D(x, y, z));
290 }
291
292 assert(points_.size() == countPoints);
293 }
294 }
295
296 const std::string& GetRoiName() const
297 {
298 return roiName_;
299 }
300
301 const std::string& GetReferencedSopInstanceUid() const
302 {
303 return referencedSopInstanceUid_;
304 }
305
306 size_t GetPointsCount() const
307 {
308 return points_.size();
309 }
310
311 const Vector3D& GetPoint(size_t i) const
312 {
313 if (i >= points_.size())
314 {
315 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
316 }
317 else
318 {
319 return points_[i];
320 }
321 }
322
323 bool IsCoplanar(Vector3D& normal) const
324 {
325 if (points_.size() < 3)
326 {
327 return false;
328 }
329
330 bool hasNormal = false;
331
332 for (size_t i = 0; i < points_.size(); i++)
333 {
334 normal = Vector3D::CrossProduct(Vector3D(points_[1], points_[0]),
335 Vector3D(points_[2], points_[0]));
336 if (!Toolbox::IsNear(normal.ComputeNorm(), 0))
337 {
338 normal.Normalize();
339 hasNormal = true;
340 }
341 }
342
343 if (!hasNormal)
344 {
345 return false;
346 }
347
348 double a = Vector3D::DotProduct(points_[0], normal);
349
350 for (size_t i = 1; i < points_.size(); i++)
351 {
352 double b = Vector3D::DotProduct(points_[i], normal);
353 if (!Toolbox::IsNear(a, b))
354 {
355 return false;
356 }
357 }
358
359 return true;
360 }
361
362 void Add(Extent2D& extent,
363 const Vector3D& axisX,
364 const Vector3D& axisY) const
365 {
366 assert(Toolbox::IsNear(1, axisX.ComputeNorm()));
367 assert(Toolbox::IsNear(1, axisY.ComputeNorm()));
368
369 for (size_t i = 0; i < points_.size(); i++)
370 {
371 extent.Add(Vector3D::DotProduct(axisX, points_[i]),
372 Vector3D::DotProduct(axisY, points_[i]));
373 }
374 }
375 };
376
377 156
378 157
379 class StructureSet : public boost::noncopyable 158 class StructureSet : public boost::noncopyable
380 { 159 {
381 private: 160 private:
390 public: 169 public:
391 explicit StructureSet(Orthanc::ParsedDicomFile& dicom) : 170 explicit StructureSet(Orthanc::ParsedDicomFile& dicom) :
392 hasFrameOfReferenceUid_(false) 171 hasFrameOfReferenceUid_(false)
393 { 172 {
394 DcmDataset& dataset = *dicom.GetDcmtkObject().getDataset(); 173 DcmDataset& dataset = *dicom.GetDcmtkObject().getDataset();
395 patientId_ = GetStringValue(dataset, DCM_PatientID); 174 patientId_ = STLToolbox::GetStringValue(dataset, DCM_PatientID);
396 studyInstanceUid_ = GetStringValue(dataset, DCM_StudyInstanceUID); 175 studyInstanceUid_ = STLToolbox::GetStringValue(dataset, DCM_StudyInstanceUID);
397 seriesInstanceUid_ = GetStringValue(dataset, DCM_SeriesInstanceUID); 176 seriesInstanceUid_ = STLToolbox::GetStringValue(dataset, DCM_SeriesInstanceUID);
398 sopInstanceUid_ = GetStringValue(dataset, DCM_SOPInstanceUID); 177 sopInstanceUid_ = STLToolbox::GetStringValue(dataset, DCM_SOPInstanceUID);
399 178
400 DcmSequenceOfItems* frame = NULL; 179 DcmSequenceOfItems* frame = NULL;
401 if (!dataset.findAndGetSequence(DCM_ReferencedFrameOfReferenceSequence, frame).good() || 180 if (!dataset.findAndGetSequence(DCM_ReferencedFrameOfReferenceSequence, frame).good() ||
402 frame == NULL) 181 frame == NULL)
403 { 182 {
519 return frameOfReferenceUid_; 298 return frameOfReferenceUid_;
520 } 299 }
521 else 300 else
522 { 301 {
523 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); 302 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
303 }
304 }
305
306 // This static method is faster than constructing the full "StructureSet" object
307 static void ListStructuresNames(std::set<std::string>& target,
308 Orthanc::ParsedDicomFile& source)
309 {
310 target.clear();
311
312 DcmSequenceOfItems* sequence = NULL;
313 if (!source.GetDcmtkObject().getDataset()->findAndGetSequence(DCM_StructureSetROISequence, sequence).good() ||
314 sequence == NULL)
315 {
316 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
317 }
318
319 for (unsigned long i = 0; i < sequence->card(); i++)
320 {
321 DcmItem* item = sequence->getItem(i);
322 if (item == NULL)
323 {
324 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
325 }
326 else
327 {
328 target.insert(STLToolbox::GetStringValue(*item, DCM_ROIName));
329 }
524 } 330 }
525 } 331 }
526 }; 332 };
527 333
528 334
553 assert(slicesSpacing_ > 0 && 359 assert(slicesSpacing_ > 0 &&
554 minProjectionAlongNormal_ < maxProjectionAlongNormal_); 360 minProjectionAlongNormal_ < maxProjectionAlongNormal_);
555 361
556 double d = (z - minProjectionAlongNormal_) / slicesSpacing_; 362 double d = (z - minProjectionAlongNormal_) / slicesSpacing_;
557 363
558 if (Toolbox::IsNear(d, round(d))) 364 if (STLToolbox::IsNear(d, round(d)))
559 { 365 {
560 if (d < 0.0 || 366 if (d < 0.0 ||
561 d > static_cast<double>(slicesCount_) - 1.0) 367 d > static_cast<double>(slicesCount_) - 1.0)
562 { 368 {
563 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); 369 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
621 } 427 }
622 428
623 // Only keep unique projections 429 // Only keep unique projections
624 430
625 std::sort(projections.begin(), projections.end()); 431 std::sort(projections.begin(), projections.end());
626 Toolbox::RemoveDuplicateValues(projections); 432 STLToolbox::RemoveDuplicateValues(projections);
627 assert(!projections.empty()); 433 assert(!projections.empty());
628 434
629 if (projections.size() == 1) 435 if (projections.size() == 1)
630 { 436 {
631 // Volume with one single slice 437 // Volume with one single slice
648 spacings[i] = projections[i + 1] - projections[i]; 454 spacings[i] = projections[i + 1] - projections[i];
649 assert(spacings[i] > 0); 455 assert(spacings[i] > 0);
650 } 456 }
651 457
652 std::sort(spacings.begin(), spacings.end()); 458 std::sort(spacings.begin(), spacings.end());
653 Toolbox::RemoveDuplicateValues(spacings); 459 STLToolbox::RemoveDuplicateValues(spacings);
654 460
655 if (spacings.empty()) 461 if (spacings.empty())
656 { 462 {
657 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); 463 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
658 } 464 }
685 it++; 491 it++;
686 492
687 while (it != candidates.end()) 493 while (it != candidates.end())
688 { 494 {
689 double d = (projections[*it] - projections[reference]) / slicesSpacing_; 495 double d = (projections[*it] - projections[reference]) / slicesSpacing_;
690 if (Toolbox::IsNear(d, round(d))) 496 if (STLToolbox::IsNear(d, round(d)))
691 { 497 {
692 countSupport ++; 498 countSupport ++;
693 } 499 }
694 else 500 else
695 { 501 {
724 maxProjectionAlongNormal_ = bestProjection; 530 maxProjectionAlongNormal_ = bestProjection;
725 531
726 for (size_t i = 0; i < projections.size(); i++) 532 for (size_t i = 0; i < projections.size(); i++)
727 { 533 {
728 double d = (projections[i] - bestProjection) / slicesSpacing_; 534 double d = (projections[i] - bestProjection) / slicesSpacing_;
729 if (Toolbox::IsNear(d, round(d))) 535 if (STLToolbox::IsNear(d, round(d)))
730 { 536 {
731 minProjectionAlongNormal_ = std::min(minProjectionAlongNormal_, projections[i]); 537 minProjectionAlongNormal_ = std::min(minProjectionAlongNormal_, projections[i]);
732 maxProjectionAlongNormal_ = std::max(maxProjectionAlongNormal_, projections[i]); 538 maxProjectionAlongNormal_ = std::max(maxProjectionAlongNormal_, projections[i]);
733 } 539 }
734 } 540 }
735 541
736 double d = (maxProjectionAlongNormal_ - minProjectionAlongNormal_) / slicesSpacing_; 542 double d = (maxProjectionAlongNormal_ - minProjectionAlongNormal_) / slicesSpacing_;
737 if (Toolbox::IsNear(d, round(d))) 543 if (STLToolbox::IsNear(d, round(d)))
738 { 544 {
739 slicesCount_ = static_cast<size_t>(round(d)) + 1; 545 slicesCount_ = static_cast<size_t>(round(d)) + 1;
740 } 546 }
741 else 547 else
742 { 548 {
900 Orthanc::Toolbox::TokenizeString(items, imageOrientation, '\\'); 706 Orthanc::Toolbox::TokenizeString(items, imageOrientation, '\\');
901 707
902 double x1, x2, x3, y1, y2, y3; 708 double x1, x2, x3, y1, y2, y3;
903 709
904 if (items.size() == 6 && 710 if (items.size() == 6 &&
905 Toolbox::MyParseDouble(x1, items[0]) && 711 STLToolbox::MyParseDouble(x1, items[0]) &&
906 Toolbox::MyParseDouble(x2, items[1]) && 712 STLToolbox::MyParseDouble(x2, items[1]) &&
907 Toolbox::MyParseDouble(x3, items[2]) && 713 STLToolbox::MyParseDouble(x3, items[2]) &&
908 Toolbox::MyParseDouble(y1, items[3]) && 714 STLToolbox::MyParseDouble(y1, items[3]) &&
909 Toolbox::MyParseDouble(y2, items[4]) && 715 STLToolbox::MyParseDouble(y2, items[4]) &&
910 Toolbox::MyParseDouble(y3, items[5])) 716 STLToolbox::MyParseDouble(y3, items[5]))
911 { 717 {
912 axisX = Vector3D(x1, x2, x3); 718 axisX = Vector3D(x1, x2, x3);
913 axisY = Vector3D(y1, y2, y3); 719 axisY = Vector3D(y1, y2, y3);
914 } 720 }
915 } 721 }
936 structureSet.GetPolygonsCount() < 1) 742 structureSet.GetPolygonsCount() < 1)
937 { 743 {
938 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); 744 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
939 } 745 }
940 746
941 if (!Toolbox::IsNear(1, geometry.GetSlicesNormal().ComputeNorm())) 747 if (!STLToolbox::IsNear(1, geometry.GetSlicesNormal().ComputeNorm()))
942 { 748 {
943 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); 749 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
944 } 750 }
945 751
946 Vector3D axisX, axisY; 752 Vector3D axisX, axisY;
947 GetReferencedVolumeAxes(axisX, axisY, structureSet); 753 GetReferencedVolumeAxes(axisX, axisY, structureSet);
948 754
949 Vector3D axisZ = Vector3D::CrossProduct(axisX, axisY); 755 Vector3D axisZ = Vector3D::CrossProduct(axisX, axisY);
950 756
951 if (!Toolbox::IsNear(1, axisX.ComputeNorm()) || 757 if (!STLToolbox::IsNear(1, axisX.ComputeNorm()) ||
952 !Toolbox::IsNear(1, axisY.ComputeNorm()) || 758 !STLToolbox::IsNear(1, axisY.ComputeNorm()) ||
953 !Vector3D::AreParallel(axisZ, geometry.GetSlicesNormal())) 759 !Vector3D::AreParallel(axisZ, geometry.GetSlicesNormal()))
954 { 760 {
955 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); 761 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
956 } 762 }
957 763
1035 const std::string instanceId(request->groups[0]); 841 const std::string instanceId(request->groups[0]);
1036 842
1037 std::unique_ptr<Orthanc::ParsedDicomFile> dicom(LoadInstance(instanceId)); 843 std::unique_ptr<Orthanc::ParsedDicomFile> dicom(LoadInstance(instanceId));
1038 844
1039 std::set<std::string> names; 845 std::set<std::string> names;
1040 ListStructuresNames(names, *dicom); 846 StructureSet::ListStructuresNames(names, *dicom);
1041 847
1042 Json::Value answer = Json::arrayValue; 848 Json::Value answer = Json::arrayValue;
1043 849
1044 for (std::set<std::string>::const_iterator it = names.begin(); it != names.end(); ++it) 850 for (std::set<std::string>::const_iterator it = names.begin(); it != names.end(); ++it)
1045 { 851 {
1258 1064
1259 std::unique_ptr<Orthanc::ParsedDicomFile> dicom(LoadInstance(instanceId)); 1065 std::unique_ptr<Orthanc::ParsedDicomFile> dicom(LoadInstance(instanceId));
1260 DcmDataset& dataset = *dicom->GetDcmtkObject().getDataset(); 1066 DcmDataset& dataset = *dicom->GetDcmtkObject().getDataset();
1261 1067
1262 std::string stl; 1068 std::string stl;
1263 if (GetStringValue(dataset, DCM_MIMETypeOfEncapsulatedDocument) != Orthanc::MIME_STL || 1069 if (STLToolbox::GetStringValue(dataset, DCM_MIMETypeOfEncapsulatedDocument) != Orthanc::MIME_STL ||
1264 GetStringValue(dataset, DCM_SOPClassUID) != UID_EncapsulatedSTLStorage || 1070 STLToolbox::GetStringValue(dataset, DCM_SOPClassUID) != UID_EncapsulatedSTLStorage ||
1265 !dicom->GetTagValue(stl, Orthanc::DICOM_TAG_ENCAPSULATED_DOCUMENT)) 1071 !dicom->GetTagValue(stl, Orthanc::DICOM_TAG_ENCAPSULATED_DOCUMENT))
1266 { 1072 {
1267 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadRequest, "DICOM instance not encapsulating a STL model: " + instanceId); 1073 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadRequest, "DICOM instance not encapsulating a STL model: " + instanceId);
1268 } 1074 }
1269 else 1075 else