Mercurial > hg > orthanc-stone

comparison Framework/Toolbox/DicomInstanceParameters.cpp @ 860:238693c3bc51 am-dev

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merge default -> am-dev
author Alain Mazy <alain@mazy.be>
date Mon, 24 Jun 2019 14:35:00 +0200
parents 55411e7da2f7
children 32eaf4929b08
comparison
equal deleted inserted replaced
856:a6e17a5a39e7 860:238693c3bc51
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-2019 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
28 #include <Core/Images/Image.h>
29 #include <Core/Images/ImageProcessing.h>
30 #include <Core/Logging.h>
31 #include <Core/OrthancException.h>
32 #include <Core/Toolbox.h>
33
34
35 namespace OrthancStone
36 {
37 void DicomInstanceParameters::Data::ComputeDoseOffsets(const Orthanc::DicomMap& dicom)
38 {
39 // http://dicom.nema.org/medical/Dicom/2016a/output/chtml/part03/sect_C.8.8.3.2.html
40
41 {
42 std::string increment;
43
44 if (dicom.CopyToString(increment, Orthanc::DICOM_TAG_FRAME_INCREMENT_POINTER, false))
45 {
46 Orthanc::Toolbox::ToUpperCase(increment);
47 if (increment != "3004,000C") // This is the "Grid Frame Offset Vector" tag
48 {
49 LOG(ERROR) << "RT-DOSE: Bad value for the \"FrameIncrementPointer\" tag";
50 return;
51 }
52 }
53 }
54
55 if (!LinearAlgebra::ParseVector(frameOffsets_, dicom, Orthanc::DICOM_TAG_GRID_FRAME_OFFSET_VECTOR) ||
56 frameOffsets_.size() < imageInformation_.GetNumberOfFrames())
57 {
58 LOG(ERROR) << "RT-DOSE: No information about the 3D location of some slice(s)";
59 frameOffsets_.clear();
60 }
61 else
62 {
63 if (frameOffsets_.size() >= 2)
64 {
65 thickness_ = std::abs(frameOffsets_[1] - frameOffsets_[0]);
66 }
67 }
68 }
69
70
71 DicomInstanceParameters::Data::Data(const Orthanc::DicomMap& dicom) :
72 imageInformation_(dicom)
73 {
74 if (imageInformation_.GetNumberOfFrames() <= 0)
75 {
76 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
77 }
78
79 if (!dicom.CopyToString(studyInstanceUid_, Orthanc::DICOM_TAG_STUDY_INSTANCE_UID, false) ||
80 !dicom.CopyToString(seriesInstanceUid_, Orthanc::DICOM_TAG_SERIES_INSTANCE_UID, false) ||
81 !dicom.CopyToString(sopInstanceUid_, Orthanc::DICOM_TAG_SOP_INSTANCE_UID, false))
82 {
83 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
84 }
85
86 std::string s;
87 if (!dicom.CopyToString(s, Orthanc::DICOM_TAG_SOP_CLASS_UID, false))
88 {
89 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
90 }
91 else
92 {
93 sopClassUid_ = StringToSopClassUid(s);
94 }
95
96 if (!dicom.ParseDouble(thickness_, Orthanc::DICOM_TAG_SLICE_THICKNESS))
97 {
98 thickness_ = 100.0 * std::numeric_limits<double>::epsilon();
99 }
100
101 GeometryToolbox::GetPixelSpacing(pixelSpacingX_, pixelSpacingY_, dicom);
102
103 std::string position, orientation;
104 if (dicom.CopyToString(position, Orthanc::DICOM_TAG_IMAGE_POSITION_PATIENT, false) &&
105 dicom.CopyToString(orientation, Orthanc::DICOM_TAG_IMAGE_ORIENTATION_PATIENT, false))
106 {
107 geometry_ = CoordinateSystem3D(position, orientation);
108 }
109
110 if (sopClassUid_ == SopClassUid_RTDose)
111 {
112 ComputeDoseOffsets(dicom);
113 }
114
115 isColor_ = (imageInformation_.GetPhotometricInterpretation() != Orthanc::PhotometricInterpretation_Monochrome1 &&
116 imageInformation_.GetPhotometricInterpretation() != Orthanc::PhotometricInterpretation_Monochrome2);
117
118 double doseGridScaling;
119
120 if (dicom.ParseDouble(rescaleIntercept_, Orthanc::DICOM_TAG_RESCALE_INTERCEPT) &&
121 dicom.ParseDouble(rescaleSlope_, Orthanc::DICOM_TAG_RESCALE_SLOPE))
122 {
123 hasRescale_ = true;
124 }
125 else if (dicom.ParseDouble(doseGridScaling, Orthanc::DICOM_TAG_DOSE_GRID_SCALING))
126 {
127 hasRescale_ = true;
128 rescaleIntercept_ = 0;
129 rescaleSlope_ = doseGridScaling;
130 }
131 else
132 {
133 hasRescale_ = false;
134 }
135
136 Vector c, w;
137 if (LinearAlgebra::ParseVector(c, dicom, Orthanc::DICOM_TAG_WINDOW_CENTER) &&
138 LinearAlgebra::ParseVector(w, dicom, Orthanc::DICOM_TAG_WINDOW_WIDTH) &&
139 c.size() > 0 &&
140 w.size() > 0)
141 {
142 hasDefaultWindowing_ = true;
143 defaultWindowingCenter_ = static_cast<float>(c[0]);
144 defaultWindowingWidth_ = static_cast<float>(w[0]);
145 }
146 else
147 {
148 hasDefaultWindowing_ = false;
149 }
150
151 if (sopClassUid_ == SopClassUid_RTDose)
152 {
153 switch (imageInformation_.GetBitsStored())
154 {
155 case 16:
156 expectedPixelFormat_ = Orthanc::PixelFormat_Grayscale16;
157 break;
158
159 case 32:
160 expectedPixelFormat_ = Orthanc::PixelFormat_Grayscale32;
161 break;
162
163 default:
164 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
165 }
166 }
167 else if (isColor_)
168 {
169 expectedPixelFormat_ = Orthanc::PixelFormat_RGB24;
170 }
171 else if (imageInformation_.IsSigned())
172 {
173 expectedPixelFormat_ = Orthanc::PixelFormat_SignedGrayscale16;
174 }
175 else
176 {
177 expectedPixelFormat_ = Orthanc::PixelFormat_Grayscale16;
178 }
179 }
180
181
182 CoordinateSystem3D DicomInstanceParameters::Data::GetFrameGeometry(unsigned int frame) const
183 {
184 if (frame == 0)
185 {
186 return geometry_;
187 }
188 else if (frame >= imageInformation_.GetNumberOfFrames())
189 {
190 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
191 }
192 else if (sopClassUid_ == SopClassUid_RTDose)
193 {
194 if (frame >= frameOffsets_.size())
195 {
196 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
197 }
198
199 return CoordinateSystem3D(
200 geometry_.GetOrigin() + frameOffsets_[frame] * geometry_.GetNormal(),
201 geometry_.GetAxisX(),
202 geometry_.GetAxisY());
203 }
204 else
205 {
206 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
207 }
208 }
209
210
211 bool DicomInstanceParameters::Data::IsPlaneWithinSlice(unsigned int frame,
212 const CoordinateSystem3D& plane) const
213 {
214 if (frame >= imageInformation_.GetNumberOfFrames())
215 {
216 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
217 }
218
219 CoordinateSystem3D tmp = geometry_;
220
221 if (frame != 0)
222 {
223 tmp = GetFrameGeometry(frame);
224 }
225
226 double distance;
227
228 return (CoordinateSystem3D::ComputeDistance(distance, tmp, plane) &&
229 distance <= thickness_ / 2.0);
230 }
231
232
233 void DicomInstanceParameters::Data::ApplyRescale(Orthanc::ImageAccessor& image,
234 bool useDouble) const
235 {
236 if (image.GetFormat() != Orthanc::PixelFormat_Float32)
237 {
238 throw Orthanc::OrthancException(Orthanc::ErrorCode_IncompatibleImageFormat);
239 }
240
241 if (hasRescale_)
242 {
243 const unsigned int width = image.GetWidth();
244 const unsigned int height = image.GetHeight();
245
246 for (unsigned int y = 0; y < height; y++)
247 {
248 float* p = reinterpret_cast<float*>(image.GetRow(y));
249
250 if (useDouble)
251 {
252 // Slower, accurate implementation using double
253 for (unsigned int x = 0; x < width; x++, p++)
254 {
255 double value = static_cast<double>(*p);
256 *p = static_cast<float>(value * rescaleSlope_ + rescaleIntercept_);
257 }
258 }
259 else
260 {
261 // Fast, approximate implementation using float
262 for (unsigned int x = 0; x < width; x++, p++)
263 {
264 *p = (*p) * static_cast<float>(rescaleSlope_) + static_cast<float>(rescaleIntercept_);
265 }
266 }
267 }
268 }
269 }
270
271 double DicomInstanceParameters::GetRescaleIntercept() const
272 {
273 if (data_.hasRescale_)
274 {
275 return data_.rescaleIntercept_;
276 }
277 else
278 {
279 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
280 }
281 }
282
283
284 double DicomInstanceParameters::GetRescaleSlope() const
285 {
286 if (data_.hasRescale_)
287 {
288 return data_.rescaleSlope_;
289 }
290 else
291 {
292 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
293 }
294 }
295
296
297 float DicomInstanceParameters::GetDefaultWindowingCenter() const
298 {
299 if (data_.hasDefaultWindowing_)
300 {
301 return data_.defaultWindowingCenter_;
302 }
303 else
304 {
305 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
306 }
307 }
308
309
310 float DicomInstanceParameters::GetDefaultWindowingWidth() const
311 {
312 if (data_.hasDefaultWindowing_)
313 {
314 return data_.defaultWindowingWidth_;
315 }
316 else
317 {
318 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
319 }
320 }
321
322
323 Orthanc::ImageAccessor* DicomInstanceParameters::ConvertToFloat(const Orthanc::ImageAccessor& pixelData) const
324 {
325 std::auto_ptr<Orthanc::Image> converted(new Orthanc::Image(Orthanc::PixelFormat_Float32,
326 pixelData.GetWidth(),
327 pixelData.GetHeight(),
328 false));
329 Orthanc::ImageProcessing::Convert(*converted, pixelData);
330
331 // Correct rescale slope/intercept if need be
332 data_.ApplyRescale(*converted, (pixelData.GetFormat() == Orthanc::PixelFormat_Grayscale32));
333
334 return converted.release();
335 }
336
337
338
339 TextureBaseSceneLayer* DicomInstanceParameters::CreateTexture
340 (const Orthanc::ImageAccessor& pixelData) const
341 {
342 assert(sizeof(float) == 4);
343
344 Orthanc::PixelFormat sourceFormat = pixelData.GetFormat();
345
346 if (sourceFormat != GetExpectedPixelFormat())
347 {
348 throw Orthanc::OrthancException(Orthanc::ErrorCode_IncompatibleImageFormat);
349 }
350
351 if (sourceFormat == Orthanc::PixelFormat_RGB24)
352 {
353 // This is the case of a color image. No conversion has to be done.
354 return new ColorTextureSceneLayer(pixelData);
355 }
356 else
357 {
358 // This is the case of a grayscale frame. Convert it to Float32.
359 std::auto_ptr<FloatTextureSceneLayer> texture;
360
361 if (pixelData.GetFormat() == Orthanc::PixelFormat_Float32)
362 {
363 texture.reset(new FloatTextureSceneLayer(pixelData));
364 }
365 else
366 {
367 std::auto_ptr<Orthanc::ImageAccessor> converted(ConvertToFloat(pixelData));
368 texture.reset(new FloatTextureSceneLayer(*converted));
369 }
370
371 if (data_.hasDefaultWindowing_)
372 {
373 texture->SetCustomWindowing(data_.defaultWindowingCenter_,
374 data_.defaultWindowingWidth_);
375 }
376
377 return texture.release();
378 }
379 }
380
381
382 LookupTableTextureSceneLayer* DicomInstanceParameters::CreateLookupTableTexture
383 (const Orthanc::ImageAccessor& pixelData) const
384 {
385 std::auto_ptr<FloatTextureSceneLayer> texture;
386
387 if (pixelData.GetFormat() == Orthanc::PixelFormat_Float32)
388 {
389 return new LookupTableTextureSceneLayer(pixelData);
390 }
391 else
392 {
393 std::auto_ptr<Orthanc::ImageAccessor> converted(ConvertToFloat(pixelData));
394 return new LookupTableTextureSceneLayer(*converted);
395 }
396 }
397 }