Mercurial > hg > orthanc-stone

comparison OrthancStone/Sources/Toolbox/DicomInstanceParameters.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/Toolbox/DicomInstanceParameters.cpp@970ee51fe01f
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 "DicomInstanceParameters.h"
23
24 #include "../Scene2D/ColorTextureSceneLayer.h"
25 #include "../Scene2D/FloatTextureSceneLayer.h"
26 #include "../Toolbox/GeometryToolbox.h"
27 #include "../Toolbox/ImageToolbox.h"
28
29 #include <Images/Image.h>
30 #include <Images/ImageProcessing.h>
31 #include <Logging.h>
32 #include <OrthancException.h>
33 #include <Toolbox.h>
34
35
36 namespace OrthancStone
37 {
38 void DicomInstanceParameters::Data::ComputeDoseOffsets(const Orthanc::DicomMap& dicom)
39 {
40 // http://dicom.nema.org/medical/Dicom/2016a/output/chtml/part03/sect_C.8.8.3.2.html
41
42 {
43 std::string increment;
44
45 if (dicom.LookupStringValue(increment, Orthanc::DICOM_TAG_FRAME_INCREMENT_POINTER, false))
46 {
47 Orthanc::Toolbox::ToUpperCase(increment);
48 if (increment != "3004,000C") // This is the "Grid Frame Offset Vector" tag (DICOM_TAG_GRID_FRAME_OFFSET_VECTOR)
49 {
50 LOG(ERROR) << "RT-DOSE: Bad value for the \"FrameIncrementPointer\" tag";
51 return;
52 }
53 }
54 }
55
56 if (!LinearAlgebra::ParseVector(frameOffsets_, dicom, Orthanc::DICOM_TAG_GRID_FRAME_OFFSET_VECTOR) ||
57 frameOffsets_.size() < imageInformation_.GetNumberOfFrames())
58 {
59 LOG(ERROR) << "RT-DOSE: No information about the 3D location of some slice(s)";
60 frameOffsets_.clear();
61 }
62 else
63 {
64 if (frameOffsets_.size() >= 2)
65 {
66 thickness_ = std::abs(frameOffsets_[1] - frameOffsets_[0]);
67 }
68 }
69 }
70
71
72 DicomInstanceParameters::Data::Data(const Orthanc::DicomMap& dicom) :
73 imageInformation_(dicom)
74 {
75 if (imageInformation_.GetNumberOfFrames() <= 0)
76 {
77 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
78 }
79
80 if (!dicom.LookupStringValue(studyInstanceUid_, Orthanc::DICOM_TAG_STUDY_INSTANCE_UID, false) ||
81 !dicom.LookupStringValue(seriesInstanceUid_, Orthanc::DICOM_TAG_SERIES_INSTANCE_UID, false) ||
82 !dicom.LookupStringValue(sopInstanceUid_, Orthanc::DICOM_TAG_SOP_INSTANCE_UID, false))
83 {
84 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
85 }
86
87 std::string s;
88 if (!dicom.LookupStringValue(s, Orthanc::DICOM_TAG_SOP_CLASS_UID, false))
89 {
90 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
91 }
92 else
93 {
94 sopClassUid_ = StringToSopClassUid(s);
95 }
96
97 if (!dicom.ParseDouble(thickness_, Orthanc::DICOM_TAG_SLICE_THICKNESS))
98 {
99 thickness_ = 100.0 * std::numeric_limits<double>::epsilon();
100 }
101
102 GeometryToolbox::GetPixelSpacing(pixelSpacingX_, pixelSpacingY_, dicom);
103
104 std::string position, orientation;
105 if (dicom.LookupStringValue(position, Orthanc::DICOM_TAG_IMAGE_POSITION_PATIENT, false) &&
106 dicom.LookupStringValue(orientation, Orthanc::DICOM_TAG_IMAGE_ORIENTATION_PATIENT, false))
107 {
108 geometry_ = CoordinateSystem3D(position, orientation);
109 }
110
111 if (sopClassUid_ == SopClassUid_RTDose)
112 {
113 ComputeDoseOffsets(dicom);
114
115 static const Orthanc::DicomTag DICOM_TAG_DOSE_UNITS(0x3004, 0x0002);
116
117 if (!dicom.LookupStringValue(doseUnits_, DICOM_TAG_DOSE_UNITS, false))
118 {
119 LOG(ERROR) << "Tag DoseUnits (0x3004, 0x0002) is missing in " << sopInstanceUid_;
120 doseUnits_ = "";
121 }
122 }
123
124 isColor_ = (imageInformation_.GetPhotometricInterpretation() != Orthanc::PhotometricInterpretation_Monochrome1 &&
125 imageInformation_.GetPhotometricInterpretation() != Orthanc::PhotometricInterpretation_Monochrome2);
126
127 if (dicom.ParseDouble(rescaleIntercept_, Orthanc::DICOM_TAG_RESCALE_INTERCEPT) &&
128 dicom.ParseDouble(rescaleSlope_, Orthanc::DICOM_TAG_RESCALE_SLOPE))
129 {
130 if (sopClassUid_ == SopClassUid_RTDose)
131 {
132 LOG(INFO) << "DOSE HAS Rescale*: rescaleIntercept_ = " << rescaleIntercept_ << " rescaleSlope_ = " << rescaleSlope_;
133 // WE SHOULD NOT TAKE THE RESCALE VALUE INTO ACCOUNT IN THE CASE OF DOSES
134 hasRescale_ = false;
135 }
136 else
137 {
138 hasRescale_ = true;
139 }
140
141 }
142 else
143 {
144 hasRescale_ = false;
145 }
146
147 if (dicom.ParseDouble(doseGridScaling_, Orthanc::DICOM_TAG_DOSE_GRID_SCALING))
148 {
149 if (sopClassUid_ == SopClassUid_RTDose)
150 {
151 LOG(INFO) << "DOSE HAS DoseGridScaling: doseGridScaling_ = " << doseGridScaling_;
152 }
153 }
154 else
155 {
156 doseGridScaling_ = 1.0;
157 if (sopClassUid_ == SopClassUid_RTDose)
158 {
159 LOG(ERROR) << "Tag DoseGridScaling (0x3004, 0x000e) is missing in " << sopInstanceUid_ << " doseGridScaling_ will be set to 1.0";
160 }
161 }
162
163 Vector c, w;
164 if (LinearAlgebra::ParseVector(c, dicom, Orthanc::DICOM_TAG_WINDOW_CENTER) &&
165 LinearAlgebra::ParseVector(w, dicom, Orthanc::DICOM_TAG_WINDOW_WIDTH) &&
166 c.size() > 0 &&
167 w.size() > 0)
168 {
169 hasDefaultWindowing_ = true;
170 defaultWindowingCenter_ = static_cast<float>(c[0]);
171 defaultWindowingWidth_ = static_cast<float>(w[0]);
172 }
173 else
174 {
175 hasDefaultWindowing_ = false;
176 defaultWindowingCenter_ = 0;
177 defaultWindowingWidth_ = 0;
178 }
179
180 if (sopClassUid_ == SopClassUid_RTDose)
181 {
182 switch (imageInformation_.GetBitsStored())
183 {
184 case 16:
185 expectedPixelFormat_ = Orthanc::PixelFormat_Grayscale16;
186 break;
187
188 case 32:
189 expectedPixelFormat_ = Orthanc::PixelFormat_Grayscale32;
190 break;
191
192 default:
193 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
194 }
195 }
196 else if (isColor_)
197 {
198 expectedPixelFormat_ = Orthanc::PixelFormat_RGB24;
199 }
200 else if (imageInformation_.IsSigned())
201 {
202 expectedPixelFormat_ = Orthanc::PixelFormat_SignedGrayscale16;
203 }
204 else
205 {
206 expectedPixelFormat_ = Orthanc::PixelFormat_Grayscale16;
207 }
208
209 // This computes the "IndexInSeries" metadata from Orthanc (check
210 // out "Orthanc::ServerIndex::Store()")
211 hasIndexInSeries_ = (
212 dicom.ParseUnsignedInteger32(indexInSeries_, Orthanc::DICOM_TAG_INSTANCE_NUMBER) ||
213 dicom.ParseUnsignedInteger32(indexInSeries_, Orthanc::DICOM_TAG_IMAGE_INDEX));
214 }
215
216
217 CoordinateSystem3D DicomInstanceParameters::Data::GetFrameGeometry(unsigned int frame) const
218 {
219 if (frame == 0)
220 {
221 return geometry_;
222 }
223 else if (frame >= imageInformation_.GetNumberOfFrames())
224 {
225 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
226 }
227 else if (sopClassUid_ == SopClassUid_RTDose)
228 {
229 if (frame >= frameOffsets_.size())
230 {
231 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
232 }
233
234 return CoordinateSystem3D(
235 geometry_.GetOrigin() + frameOffsets_[frame] * geometry_.GetNormal(),
236 geometry_.GetAxisX(),
237 geometry_.GetAxisY());
238 }
239 else
240 {
241 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
242 }
243 }
244
245
246 bool DicomInstanceParameters::Data::IsPlaneWithinSlice(unsigned int frame,
247 const CoordinateSystem3D& plane) const
248 {
249 if (frame >= imageInformation_.GetNumberOfFrames())
250 {
251 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
252 }
253
254 CoordinateSystem3D tmp = geometry_;
255
256 if (frame != 0)
257 {
258 tmp = GetFrameGeometry(frame);
259 }
260
261 double distance;
262
263 return (CoordinateSystem3D::ComputeDistance(distance, tmp, plane) &&
264 distance <= thickness_ / 2.0);
265 }
266
267 void DicomInstanceParameters::Data::ApplyRescaleAndDoseScaling(Orthanc::ImageAccessor& image,
268 bool useDouble) const
269 {
270 if (image.GetFormat() != Orthanc::PixelFormat_Float32)
271 {
272 throw Orthanc::OrthancException(Orthanc::ErrorCode_IncompatibleImageFormat);
273 }
274
275 double factor = doseGridScaling_;
276 double offset = 0.0;
277
278 if (hasRescale_)
279 {
280 factor *= rescaleSlope_;
281 offset = rescaleIntercept_;
282 }
283
284 if ( (factor != 1.0) || (offset != 0.0) )
285 {
286 const unsigned int width = image.GetWidth();
287 const unsigned int height = image.GetHeight();
288
289 for (unsigned int y = 0; y < height; y++)
290 {
291 float* p = reinterpret_cast<float*>(image.GetRow(y));
292
293 if (useDouble)
294 {
295 // Slower, accurate implementation using double
296 for (unsigned int x = 0; x < width; x++, p++)
297 {
298 double value = static_cast<double>(*p);
299 *p = static_cast<float>(value * factor + offset);
300 }
301 }
302 else
303 {
304 // Fast, approximate implementation using float
305 for (unsigned int x = 0; x < width; x++, p++)
306 {
307 *p = (*p) * static_cast<float>(factor) + static_cast<float>(offset);
308 }
309 }
310 }
311 }
312 }
313
314 double DicomInstanceParameters::GetRescaleIntercept() const
315 {
316 if (data_.hasRescale_)
317 {
318 return data_.rescaleIntercept_;
319 }
320 else
321 {
322 LOG(ERROR) << "DicomInstanceParameters::GetRescaleIntercept(): !data_.hasRescale_";
323 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
324 }
325 }
326
327
328 double DicomInstanceParameters::GetRescaleSlope() const
329 {
330 if (data_.hasRescale_)
331 {
332 return data_.rescaleSlope_;
333 }
334 else
335 {
336 LOG(ERROR) << "DicomInstanceParameters::GetRescaleSlope(): !data_.hasRescale_";
337 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
338 }
339 }
340
341
342 float DicomInstanceParameters::GetDefaultWindowingCenter() const
343 {
344 if (data_.hasDefaultWindowing_)
345 {
346 return data_.defaultWindowingCenter_;
347 }
348 else
349 {
350 LOG(ERROR) << "DicomInstanceParameters::GetDefaultWindowingCenter(): no default windowing";
351 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
352 }
353 }
354
355
356 float DicomInstanceParameters::GetDefaultWindowingWidth() const
357 {
358 if (data_.hasDefaultWindowing_)
359 {
360 return data_.defaultWindowingWidth_;
361 }
362 else
363 {
364 LOG(ERROR) << "DicomInstanceParameters::GetDefaultWindowingWidth(): no default windowing";
365 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
366 }
367 }
368
369
370 Orthanc::ImageAccessor* DicomInstanceParameters::ConvertToFloat(const Orthanc::ImageAccessor& pixelData) const
371 {
372 std::unique_ptr<Orthanc::Image> converted(new Orthanc::Image(Orthanc::PixelFormat_Float32,
373 pixelData.GetWidth(),
374 pixelData.GetHeight(),
375 false));
376 Orthanc::ImageProcessing::Convert(*converted, pixelData);
377
378
379 // Correct rescale slope/intercept if need be
380 //data_.ApplyRescaleAndDoseScaling(*converted, (pixelData.GetFormat() == Orthanc::PixelFormat_Grayscale32));
381 data_.ApplyRescaleAndDoseScaling(*converted, false);
382
383 return converted.release();
384 }
385
386
387 TextureBaseSceneLayer* DicomInstanceParameters::CreateTexture
388 (const Orthanc::ImageAccessor& pixelData) const
389 {
390 // {
391 // const Orthanc::ImageAccessor& source = pixelData;
392 // const void* sourceBuffer = source.GetConstBuffer();
393 // intptr_t sourceBufferInt = reinterpret_cast<intptr_t>(sourceBuffer);
394 // int sourceWidth = source.GetWidth();
395 // int sourceHeight = source.GetHeight();
396 // int sourcePitch = source.GetPitch();
397
398 // // TODO: turn error into trace below
399 // LOG(ERROR) << "ConvertGrayscaleToFloat | source:"
400 // << " W = " << sourceWidth << " H = " << sourceHeight
401 // << " P = " << sourcePitch << " B = " << sourceBufferInt
402 // << " B % 4 == " << sourceBufferInt % 4;
403 // }
404
405 assert(sizeof(float) == 4);
406
407 Orthanc::PixelFormat sourceFormat = pixelData.GetFormat();
408
409 if (sourceFormat != GetExpectedPixelFormat())
410 {
411 throw Orthanc::OrthancException(Orthanc::ErrorCode_IncompatibleImageFormat);
412 }
413
414 if (sourceFormat == Orthanc::PixelFormat_RGB24)
415 {
416 // This is the case of a color image. No conversion has to be done.
417 return new ColorTextureSceneLayer(pixelData);
418 }
419 else
420 {
421 // This is the case of a grayscale frame. Convert it to Float32.
422 std::unique_ptr<FloatTextureSceneLayer> texture;
423
424 if (pixelData.GetFormat() == Orthanc::PixelFormat_Float32)
425 {
426 texture.reset(new FloatTextureSceneLayer(pixelData));
427 }
428 else
429 {
430 std::unique_ptr<Orthanc::ImageAccessor> converted(ConvertToFloat(pixelData));
431 texture.reset(new FloatTextureSceneLayer(*converted));
432 }
433
434 if (data_.hasDefaultWindowing_)
435 {
436 texture->SetCustomWindowing(data_.defaultWindowingCenter_,
437 data_.defaultWindowingWidth_);
438 }
439
440
441 if (data_.imageInformation_.GetPhotometricInterpretation()
442 == Orthanc::PhotometricInterpretation_Monochrome1)
443 {
444 texture->SetInverted(true);
445 }
446 else if (data_.imageInformation_.GetPhotometricInterpretation()
447 == Orthanc::PhotometricInterpretation_Monochrome2)
448 {
449 texture->SetInverted(false);
450 }
451
452 return texture.release();
453 }
454 }
455
456
457 LookupTableTextureSceneLayer* DicomInstanceParameters::CreateLookupTableTexture
458 (const Orthanc::ImageAccessor& pixelData) const
459 {
460 std::unique_ptr<FloatTextureSceneLayer> texture;
461
462 if (pixelData.GetFormat() == Orthanc::PixelFormat_Float32)
463 {
464 return new LookupTableTextureSceneLayer(pixelData);
465 }
466 else
467 {
468 std::unique_ptr<Orthanc::ImageAccessor> converted(ConvertToFloat(pixelData));
469 return new LookupTableTextureSceneLayer(*converted);
470 }
471 }
472
473
474 unsigned int DicomInstanceParameters::GetIndexInSeries() const
475 {
476 if (data_.hasIndexInSeries_)
477 {
478 return data_.indexInSeries_;
479 }
480 else
481 {
482 LOG(ERROR) << "DicomInstanceParameters::GetIndexInSeries(): !data_.hasIndexInSeries_";
483 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
484 }
485 }
486
487
488 double DicomInstanceParameters::Data::ApplyRescale(double value) const
489 {
490 double factor = doseGridScaling_;
491 double offset = 0.0;
492
493 if (hasRescale_)
494 {
495 factor *= rescaleSlope_;
496 offset = rescaleIntercept_;
497 }
498
499 return (value * factor + offset);
500 }
501
502
503 bool DicomInstanceParameters::Data::ComputeRegularSpacing(double& spacing) const
504 {
505 if (frameOffsets_.size() == 0) // Not a RT-DOSE
506 {
507 return false;
508 }
509 else if (frameOffsets_.size() == 1)
510 {
511 spacing = 1; // Edge case: RT-DOSE with one single frame
512 return true;
513 }
514 else
515 {
516 spacing = std::abs(frameOffsets_[1] - frameOffsets_[0]);
517
518 for (size_t i = 1; i + 1 < frameOffsets_.size(); i++)
519 {
520 double s = frameOffsets_[i + 1] - frameOffsets_[i];
521 if (!LinearAlgebra::IsNear(spacing, s, 0.001))
522 {
523 return false;
524 }
525 }
526
527 return true;
528 }
529 }
530 }