Mercurial > hg > orthanc-stone

comparison OrthancStone/Sources/Loaders/OrthancMultiframeVolumeLoader.cpp @ 1512:244ad1e4e76a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
reorganization of folders
author Sebastien Jodogne <s.jodogne@gmail.com>
date Tue, 07 Jul 2020 16:21:02 +0200
parents Framework/Loaders/OrthancMultiframeVolumeLoader.cpp@30deba7bc8e2
children 85e117739eca
comparison
equal deleted inserted replaced
1511:9dfeee74c1e6 1512:244ad1e4e76a
1 /**
2 * Stone of Orthanc
3 * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics
4 * Department, University Hospital of Liege, Belgium
5 * Copyright (C) 2017-2020 Osimis S.A., Belgium
6 *
7 * This program is free software: you can redistribute it and/or
8 * modify it under the terms of the GNU Affero General Public License
9 * as published by the Free Software Foundation, either version 3 of
10 * the License, or (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful, but
13 * WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Affero General Public License for more details.
16 *
17 * You should have received a copy of the GNU Affero General Public License
18 * along with this program. If not, see <http://www.gnu.org/licenses/>.
19 **/
20
21
22 #include "OrthancMultiframeVolumeLoader.h"
23
24 #include <Endianness.h>
25 #include <Toolbox.h>
26
27 namespace OrthancStone
28 {
29 class OrthancMultiframeVolumeLoader::LoadRTDoseGeometry : public LoaderStateMachine::State
30 {
31 private:
32 std::unique_ptr<Orthanc::DicomMap> dicom_;
33
34 public:
35 LoadRTDoseGeometry(OrthancMultiframeVolumeLoader& that,
36 Orthanc::DicomMap* dicom) :
37 State(that),
38 dicom_(dicom)
39 {
40 if (dicom == NULL)
41 {
42 throw Orthanc::OrthancException(Orthanc::ErrorCode_NullPointer);
43 }
44
45 }
46
47 virtual void Handle(const OrthancStone::OrthancRestApiCommand::SuccessMessage& message)
48 {
49 // Complete the DICOM tags with just-received "Grid Frame Offset Vector"
50 std::string s = Orthanc::Toolbox::StripSpaces(message.GetAnswer());
51 dicom_->SetValue(Orthanc::DICOM_TAG_GRID_FRAME_OFFSET_VECTOR, s, false);
52
53 GetLoader<OrthancMultiframeVolumeLoader>().SetGeometry(*dicom_);
54 }
55 };
56
57
58 static std::string GetSopClassUid(const Orthanc::DicomMap& dicom)
59 {
60 std::string s;
61 if (!dicom.LookupStringValue(s, Orthanc::DICOM_TAG_SOP_CLASS_UID, false))
62 {
63 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat,
64 "DICOM file without SOP class UID");
65 }
66 else
67 {
68 return s;
69 }
70 }
71
72
73 class OrthancMultiframeVolumeLoader::LoadGeometry : public State
74 {
75 public:
76 LoadGeometry(OrthancMultiframeVolumeLoader& that) :
77 State(that)
78 {
79 }
80
81 virtual void Handle(const OrthancStone::OrthancRestApiCommand::SuccessMessage& message)
82 {
83 OrthancMultiframeVolumeLoader& loader = GetLoader<OrthancMultiframeVolumeLoader>();
84
85 Json::Value body;
86 message.ParseJsonBody(body);
87
88 if (body.type() != Json::objectValue)
89 {
90 throw Orthanc::OrthancException(Orthanc::ErrorCode_NetworkProtocol);
91 }
92
93 std::unique_ptr<Orthanc::DicomMap> dicom(new Orthanc::DicomMap);
94 dicom->FromDicomAsJson(body);
95
96 if (OrthancStone::StringToSopClassUid(GetSopClassUid(*dicom)) == OrthancStone::SopClassUid_RTDose)
97 {
98 // Download the "Grid Frame Offset Vector" DICOM tag, that is
99 // mandatory for RT-DOSE, but is too long to be returned by default
100
101 std::unique_ptr<OrthancStone::OrthancRestApiCommand> command(new OrthancStone::OrthancRestApiCommand);
102 command->SetUri("/instances/" + loader.GetInstanceId() + "/content/" +
103 Orthanc::DICOM_TAG_GRID_FRAME_OFFSET_VECTOR.Format());
104 command->AcquirePayload(new LoadRTDoseGeometry(loader, dicom.release()));
105
106 Schedule(command.release());
107 }
108 else
109 {
110 loader.SetGeometry(*dicom);
111 }
112 }
113 };
114
115 class OrthancMultiframeVolumeLoader::LoadTransferSyntax : public State
116 {
117 public:
118 LoadTransferSyntax(OrthancMultiframeVolumeLoader& that) :
119 State(that)
120 {
121 }
122
123 virtual void Handle(const OrthancStone::OrthancRestApiCommand::SuccessMessage& message)
124 {
125 GetLoader<OrthancMultiframeVolumeLoader>().SetTransferSyntax(message.GetAnswer());
126 }
127 };
128
129 class OrthancMultiframeVolumeLoader::LoadUncompressedPixelData : public State
130 {
131 public:
132 LoadUncompressedPixelData(OrthancMultiframeVolumeLoader& that) :
133 State(that)
134 {
135 }
136
137 virtual void Handle(const OrthancStone::OrthancRestApiCommand::SuccessMessage& message)
138 {
139 GetLoader<OrthancMultiframeVolumeLoader>().SetUncompressedPixelData(message.GetAnswer());
140 }
141 };
142
143 const std::string& OrthancMultiframeVolumeLoader::GetInstanceId() const
144 {
145 if (IsActive())
146 {
147 return instanceId_;
148 }
149 else
150 {
151 LOG(ERROR) << "OrthancMultiframeVolumeLoader::GetInstanceId(): (!IsActive())";
152 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
153 }
154 }
155
156 void OrthancMultiframeVolumeLoader::ScheduleFrameDownloads()
157 {
158 if (transferSyntaxUid_.empty() ||
159 !volume_->HasGeometry())
160 {
161 return;
162 }
163 /*
164 1.2.840.10008.1.2 Implicit VR Endian: Default Transfer Syntax for DICOM
165 1.2.840.10008.1.2.1 Explicit VR Little Endian
166 1.2.840.10008.1.2.2 Explicit VR Big Endian
167
168 See https://www.dicomlibrary.com/dicom/transfer-syntax/
169 */
170 if (transferSyntaxUid_ == "1.2.840.10008.1.2" ||
171 transferSyntaxUid_ == "1.2.840.10008.1.2.1" ||
172 transferSyntaxUid_ == "1.2.840.10008.1.2.2")
173 {
174 std::unique_ptr<OrthancStone::OrthancRestApiCommand> command(new OrthancStone::OrthancRestApiCommand);
175 command->SetHttpHeader("Accept-Encoding", "gzip");
176 command->SetUri("/instances/" + instanceId_ + "/content/" +
177 Orthanc::DICOM_TAG_PIXEL_DATA.Format() + "/0");
178 command->AcquirePayload(new LoadUncompressedPixelData(*this));
179 Schedule(command.release());
180 }
181 else
182 {
183 throw Orthanc::OrthancException(
184 Orthanc::ErrorCode_NotImplemented,
185 "No support for multiframe instances with transfer syntax: " + transferSyntaxUid_);
186 }
187 }
188
189 void OrthancMultiframeVolumeLoader::SetTransferSyntax(const std::string& transferSyntax)
190 {
191 transferSyntaxUid_ = Orthanc::Toolbox::StripSpaces(transferSyntax);
192 ScheduleFrameDownloads();
193 }
194
195 void OrthancMultiframeVolumeLoader::SetGeometry(const Orthanc::DicomMap& dicom)
196 {
197 OrthancStone::DicomInstanceParameters parameters(dicom);
198 volume_->SetDicomParameters(parameters);
199
200 Orthanc::PixelFormat format;
201 if (!parameters.GetImageInformation().ExtractPixelFormat(format, true))
202 {
203 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
204 }
205
206 double spacingZ;
207 switch (parameters.GetSopClassUid())
208 {
209 case OrthancStone::SopClassUid_RTDose:
210 spacingZ = parameters.GetThickness();
211 break;
212
213 default:
214 throw Orthanc::OrthancException(
215 Orthanc::ErrorCode_NotImplemented,
216 "No support for multiframe instances with SOP class UID: " + GetSopClassUid(dicom));
217 }
218
219 const unsigned int width = parameters.GetImageInformation().GetWidth();
220 const unsigned int height = parameters.GetImageInformation().GetHeight();
221 const unsigned int depth = parameters.GetImageInformation().GetNumberOfFrames();
222
223 {
224 OrthancStone::VolumeImageGeometry geometry;
225 geometry.SetSizeInVoxels(width, height, depth);
226 geometry.SetAxialGeometry(parameters.GetGeometry());
227 geometry.SetVoxelDimensions(parameters.GetPixelSpacingX(),
228 parameters.GetPixelSpacingY(), spacingZ);
229 volume_->Initialize(geometry, format, true /* Do compute range */);
230 }
231
232 volume_->GetPixelData().Clear();
233
234 ScheduleFrameDownloads();
235
236
237
238 BroadcastMessage(OrthancStone::DicomVolumeImage::GeometryReadyMessage(*volume_));
239 }
240
241
242 ORTHANC_FORCE_INLINE
243 static void CopyPixel(uint32_t& target, const void* source)
244 {
245 // TODO - check alignement?
246 target = le32toh(*reinterpret_cast<const uint32_t*>(source));
247 }
248
249 ORTHANC_FORCE_INLINE
250 static void CopyPixel(uint16_t& target, const void* source)
251 {
252 // TODO - check alignement?
253 target = le16toh(*reinterpret_cast<const uint16_t*>(source));
254 }
255
256 ORTHANC_FORCE_INLINE
257 static void CopyPixel(int16_t& target, const void* source)
258 {
259 // byte swapping is the same for unsigned and signed integers
260 // (the sign bit is always stored with the MSByte)
261 uint16_t* targetUp = reinterpret_cast<uint16_t*>(&target);
262 CopyPixel(*targetUp, source);
263 }
264
265 template <typename T>
266 void OrthancMultiframeVolumeLoader::CopyPixelDataAndComputeDistribution(
267 const std::string& pixelData, std::map<T,uint64_t>& distribution)
268 {
269 OrthancStone::ImageBuffer3D& target = volume_->GetPixelData();
270
271 const unsigned int bpp = target.GetBytesPerPixel();
272 const unsigned int width = target.GetWidth();
273 const unsigned int height = target.GetHeight();
274 const unsigned int depth = target.GetDepth();
275
276 if (pixelData.size() != bpp * width * height * depth)
277 {
278 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat,
279 "The pixel data has not the proper size");
280 }
281
282 if (pixelData.empty())
283 {
284 return;
285 }
286
287 // first pass to initialize map
288 {
289 const uint8_t* source = reinterpret_cast<const uint8_t*>(pixelData.c_str());
290
291 for (unsigned int z = 0; z < depth; z++)
292 {
293 for (unsigned int y = 0; y < height; y++)
294 {
295 for (unsigned int x = 0; x < width; x++)
296 {
297 T value;
298 CopyPixel(value, source);
299 distribution[value] = 0;
300 source += bpp;
301 }
302 }
303 }
304 }
305
306 {
307 const uint8_t* source = reinterpret_cast<const uint8_t*>(pixelData.c_str());
308
309 for (unsigned int z = 0; z < depth; z++)
310 {
311 OrthancStone::ImageBuffer3D::SliceWriter writer(target, OrthancStone::VolumeProjection_Axial, z);
312
313 assert(writer.GetAccessor().GetWidth() == width &&
314 writer.GetAccessor().GetHeight() == height);
315 #if 0
316 for (unsigned int y = 0; y < height; y++)
317 {
318 assert(sizeof(T) == Orthanc::GetBytesPerPixel(target.GetFormat()));
319
320 T* target = reinterpret_cast<T*>(writer.GetAccessor().GetRow(y));
321
322 for (unsigned int x = 0; x < width; x++)
323 {
324 CopyPixel(*target, source);
325
326 distribution[*target] += 1;
327
328 target++;
329 source += bpp;
330 }
331 }
332 #else
333 // optimized version (fixed) as of 2020-04-15
334 unsigned int pitch = writer.GetAccessor().GetPitch();
335 T* targetAddrLine = reinterpret_cast<T*>(writer.GetAccessor().GetRow(0));
336 assert(sizeof(T) == Orthanc::GetBytesPerPixel(target.GetFormat()));
337
338 for (unsigned int y = 0; y < height; y++)
339 {
340 T* targetAddrPix = targetAddrLine;
341 for (unsigned int x = 0; x < width; x++)
342 {
343 CopyPixel(*targetAddrPix, source);
344
345 distribution[*targetAddrPix] += 1;
346
347 targetAddrPix++;
348 source += bpp;
349 }
350 uint8_t* targetAddrLineBytes = reinterpret_cast<uint8_t*>(targetAddrLine) + pitch;
351 targetAddrLine = reinterpret_cast<T*>(targetAddrLineBytes);
352 }
353 #endif
354 }
355 }
356 }
357
358 template <typename T>
359 void OrthancMultiframeVolumeLoader::ComputeMinMaxWithOutlierRejection(
360 const std::map<T, uint64_t>& distribution)
361 {
362 if (distribution.size() == 0)
363 {
364 LOG(ERROR) << "ComputeMinMaxWithOutlierRejection -- Volume image empty.";
365 }
366 else
367 {
368 OrthancStone::ImageBuffer3D& target = volume_->GetPixelData();
369
370 const uint64_t width = target.GetWidth();
371 const uint64_t height = target.GetHeight();
372 const uint64_t depth = target.GetDepth();
373 const uint64_t voxelCount = width * height * depth;
374
375 // now that we have distribution[pixelValue] == numberOfPixelsWithValue
376 // compute number of values and check (assertion) that it is equal to
377 // width * height * depth
378 {
379 typename std::map<T, uint64_t>::const_iterator it = distribution.begin();
380 uint64_t totalCount = 0;
381 distributionRawMin_ = static_cast<float>(it->first);
382
383 while (it != distribution.end())
384 {
385 T pixelValue = it->first;
386 uint64_t count = it->second;
387 totalCount += count;
388 it++;
389 if (it == distribution.end())
390 distributionRawMax_ = static_cast<float>(pixelValue);
391 }
392 LOG(INFO) << "Volume image. First distribution value = "
393 << static_cast<float>(distributionRawMin_)
394 << " | Last distribution value = "
395 << static_cast<float>(distributionRawMax_);
396
397 if (totalCount != voxelCount)
398 {
399 LOG(ERROR) << "Internal error in dose distribution computation. TC ("
400 << totalCount << ") != VoxC (" << voxelCount;
401 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
402 }
403 }
404
405 // compute the number of voxels to reject at each end of the distribution
406 uint64_t endRejectionCount = static_cast<uint64_t>(
407 outliersHalfRejectionRate_ * voxelCount);
408
409 if (endRejectionCount > voxelCount)
410 {
411 LOG(ERROR) << "Internal error in dose distribution computation."
412 << " endRejectionCount = " << endRejectionCount
413 << " | voxelCount = " << voxelCount;
414 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
415 }
416
417 // this will contain the actual distribution minimum after outlier
418 // rejection
419 T resultMin = 0;
420
421 // then start from start and remove pixel values up to
422 // endRejectionCount voxels rejected
423 {
424 typename std::map<T, uint64_t>::const_iterator it = distribution.begin();
425
426 uint64_t currentCount = 0;
427
428 while (it != distribution.end())
429 {
430 T pixelValue = it->first;
431 uint64_t count = it->second;
432
433 // if this pixelValue crosses the rejection threshold, let's set it
434 // and exit the loop
435 if ((currentCount <= endRejectionCount) &&
436 (currentCount + count > endRejectionCount))
437 {
438 resultMin = pixelValue;
439 break;
440 }
441 else
442 {
443 currentCount += count;
444 }
445 // and continue walking along the distribution
446 it++;
447 }
448 }
449
450 // this will contain the actual distribution maximum after outlier
451 // rejection
452 T resultMax = 0;
453 // now start from END and remove pixel values up to
454 // endRejectionCount voxels rejected
455 {
456 typename std::map<T, uint64_t>::const_reverse_iterator it = distribution.rbegin();
457
458 uint64_t currentCount = 0;
459
460 while (it != distribution.rend())
461 {
462 T pixelValue = it->first;
463 uint64_t count = it->second;
464
465 if ((currentCount <= endRejectionCount) &&
466 (currentCount + count > endRejectionCount))
467 {
468 resultMax = pixelValue;
469 break;
470 }
471 else
472 {
473 currentCount += count;
474 }
475 // and continue walking along the distribution
476 it++;
477 }
478 }
479 if (resultMin > resultMax)
480 {
481 LOG(ERROR) << "Internal error in dose distribution computation! " <<
482 "resultMin (" << resultMin << ") > resultMax (" << resultMax << ")";
483 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
484 }
485 computedDistributionMin_ = static_cast<float>(resultMin);
486 computedDistributionMax_ = static_cast<float>(resultMax);
487 }
488 }
489
490 template <typename T>
491 void OrthancMultiframeVolumeLoader::CopyPixelDataAndComputeMinMax(
492 const std::string& pixelData)
493 {
494 std::map<T, uint64_t> distribution;
495 CopyPixelDataAndComputeDistribution(pixelData, distribution);
496 ComputeMinMaxWithOutlierRejection(distribution);
497 }
498
499 void OrthancMultiframeVolumeLoader::SetUncompressedPixelData(const std::string& pixelData)
500 {
501 switch (volume_->GetPixelData().GetFormat())
502 {
503 case Orthanc::PixelFormat_Grayscale32:
504 CopyPixelDataAndComputeMinMax<uint32_t>(pixelData);
505 break;
506 case Orthanc::PixelFormat_Grayscale16:
507 CopyPixelDataAndComputeMinMax<uint16_t>(pixelData);
508 break;
509 case Orthanc::PixelFormat_SignedGrayscale16:
510 CopyPixelDataAndComputeMinMax<int16_t>(pixelData);
511 break;
512 default:
513 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
514 }
515
516 volume_->IncrementRevision();
517
518 pixelDataLoaded_ = true;
519 BroadcastMessage(OrthancStone::DicomVolumeImage::ContentUpdatedMessage(*volume_));
520 }
521
522 bool OrthancMultiframeVolumeLoader::HasGeometry() const
523 {
524 return volume_->HasGeometry();
525 }
526
527 const OrthancStone::VolumeImageGeometry& OrthancMultiframeVolumeLoader::GetImageGeometry() const
528 {
529 return volume_->GetGeometry();
530 }
531
532 OrthancMultiframeVolumeLoader::OrthancMultiframeVolumeLoader(
533 OrthancStone::ILoadersContext& loadersContext,
534 boost::shared_ptr<OrthancStone::DicomVolumeImage> volume,
535 float outliersHalfRejectionRate)
536 : LoaderStateMachine(loadersContext)
537 , volume_(volume)
538 , pixelDataLoaded_(false)
539 , outliersHalfRejectionRate_(outliersHalfRejectionRate)
540 , distributionRawMin_(0)
541 , distributionRawMax_(0)
542 , computedDistributionMin_(0)
543 , computedDistributionMax_(0)
544 {
545 if (volume.get() == NULL)
546 {
547 throw Orthanc::OrthancException(Orthanc::ErrorCode_NullPointer);
548 }
549 }
550
551
552 boost::shared_ptr<OrthancMultiframeVolumeLoader>
553 OrthancMultiframeVolumeLoader::Create(
554 OrthancStone::ILoadersContext& loadersContext,
555 boost::shared_ptr<OrthancStone::DicomVolumeImage> volume,
556 float outliersHalfRejectionRate /*= 0.0005*/)
557 {
558 boost::shared_ptr<OrthancMultiframeVolumeLoader> obj(
559 new OrthancMultiframeVolumeLoader(
560 loadersContext,
561 volume,
562 outliersHalfRejectionRate));
563 obj->LoaderStateMachine::PostConstructor();
564 return obj;
565 }
566
567 OrthancMultiframeVolumeLoader::~OrthancMultiframeVolumeLoader()
568 {
569 LOG(TRACE) << "OrthancMultiframeVolumeLoader::~OrthancMultiframeVolumeLoader()";
570 }
571
572 void OrthancMultiframeVolumeLoader::GetDistributionMinMax
573 (float& minValue, float& maxValue) const
574 {
575 if (distributionRawMin_ == 0 && distributionRawMax_ == 0)
576 {
577 LOG(WARNING) << "GetDistributionMinMaxWithOutliersRejection called before computation!";
578 }
579 minValue = distributionRawMin_;
580 maxValue = distributionRawMax_;
581 }
582
583 void OrthancMultiframeVolumeLoader::GetDistributionMinMaxWithOutliersRejection
584 (float& minValue, float& maxValue) const
585 {
586 if (computedDistributionMin_ == 0 && computedDistributionMax_ == 0)
587 {
588 LOG(WARNING) << "GetDistributionMinMaxWithOutliersRejection called before computation!";
589 }
590 minValue = computedDistributionMin_;
591 maxValue = computedDistributionMax_;
592 }
593
594 void OrthancMultiframeVolumeLoader::LoadInstance(const std::string& instanceId)
595 {
596 Start();
597
598 instanceId_ = instanceId;
599
600 {
601 std::unique_ptr<OrthancStone::OrthancRestApiCommand> command(new OrthancStone::OrthancRestApiCommand);
602 command->SetHttpHeader("Accept-Encoding", "gzip");
603 command->SetUri("/instances/" + instanceId + "/tags");
604 command->AcquirePayload(new LoadGeometry(*this));
605 Schedule(command.release());
606 }
607
608 {
609 std::unique_ptr<OrthancStone::OrthancRestApiCommand> command(new OrthancStone::OrthancRestApiCommand);
610 command->SetUri("/instances/" + instanceId + "/metadata/TransferSyntax");
611 command->AcquirePayload(new LoadTransferSyntax(*this));
612 Schedule(command.release());
613 }
614 }
615 }