Mercurial > hg > orthanc-stone
comparison Framework/dev.h @ 102:fcec0ab44054 wasm
display volumes
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Wed, 31 May 2017 17:01:18 +0200 |
parents | |
children | eccd64f8e297 |
comparison
equal
deleted
inserted
replaced
101:af312ce4fe59 | 102:fcec0ab44054 |
---|---|
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 Osimis, 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 #pragma once | |
23 | |
24 #include "Layers/FrameRenderer.h" | |
25 #include "Layers/LayerSourceBase.h" | |
26 #include "Layers/SliceOutlineRenderer.h" | |
27 #include "Toolbox/DownloadStack.h" | |
28 #include "Toolbox/OrthancSlicesLoader.h" | |
29 #include "Volumes/ImageBuffer3D.h" | |
30 #include "Volumes/SlicedVolumeBase.h" | |
31 | |
32 #include "../Resources/Orthanc/Core/Logging.h" | |
33 #include "../Resources/Orthanc/Core/Images/ImageProcessing.h" | |
34 | |
35 #include <boost/math/special_functions/round.hpp> | |
36 | |
37 | |
38 namespace OrthancStone | |
39 { | |
40 class OrthancVolumeImage : | |
41 public SlicedVolumeBase, | |
42 private OrthancSlicesLoader::ICallback | |
43 { | |
44 private: | |
45 OrthancSlicesLoader loader_; | |
46 std::auto_ptr<ImageBuffer3D> image_; | |
47 std::auto_ptr<DownloadStack> downloadStack_; | |
48 | |
49 | |
50 void ScheduleSliceDownload() | |
51 { | |
52 assert(downloadStack_.get() != NULL); | |
53 | |
54 unsigned int slice; | |
55 if (downloadStack_->Pop(slice)) | |
56 { | |
57 loader_.ScheduleLoadSliceImage(slice, SliceImageQuality_Jpeg90); | |
58 } | |
59 } | |
60 | |
61 | |
62 static bool IsCompatible(const Slice& a, | |
63 const Slice& b) | |
64 { | |
65 if (!GeometryToolbox::IsParallel(a.GetGeometry().GetNormal(), | |
66 b.GetGeometry().GetNormal())) | |
67 { | |
68 LOG(ERROR) << "Some slice in the volume image is not parallel to the others"; | |
69 return false; | |
70 } | |
71 | |
72 if (a.GetConverter().GetExpectedPixelFormat() != b.GetConverter().GetExpectedPixelFormat()) | |
73 { | |
74 LOG(ERROR) << "The pixel format changes across the slices of the volume image"; | |
75 return false; | |
76 } | |
77 | |
78 if (a.GetWidth() != b.GetWidth() || | |
79 a.GetHeight() != b.GetHeight()) | |
80 { | |
81 LOG(ERROR) << "The width/height of the slices change across the volume image"; | |
82 return false; | |
83 } | |
84 | |
85 if (!GeometryToolbox::IsNear(a.GetPixelSpacingX(), b.GetPixelSpacingX()) || | |
86 !GeometryToolbox::IsNear(a.GetPixelSpacingY(), b.GetPixelSpacingY())) | |
87 { | |
88 LOG(ERROR) << "The pixel spacing of the slices change across the volume image"; | |
89 return false; | |
90 } | |
91 | |
92 return true; | |
93 } | |
94 | |
95 | |
96 static double GetDistance(const Slice& a, | |
97 const Slice& b) | |
98 { | |
99 return fabs(a.GetGeometry().ProjectAlongNormal(a.GetGeometry().GetOrigin()) - | |
100 a.GetGeometry().ProjectAlongNormal(b.GetGeometry().GetOrigin())); | |
101 } | |
102 | |
103 | |
104 virtual void NotifyGeometryReady(const OrthancSlicesLoader& loader) | |
105 { | |
106 if (loader.GetSliceCount() == 0) | |
107 { | |
108 LOG(ERROR) << "Empty volume image"; | |
109 SlicedVolumeBase::NotifyGeometryError(); | |
110 return; | |
111 } | |
112 | |
113 for (size_t i = 1; i < loader.GetSliceCount(); i++) | |
114 { | |
115 if (!IsCompatible(loader.GetSlice(0), loader.GetSlice(i))) | |
116 { | |
117 SlicedVolumeBase::NotifyGeometryError(); | |
118 return; | |
119 } | |
120 } | |
121 | |
122 double spacingZ; | |
123 | |
124 if (loader.GetSliceCount() > 1) | |
125 { | |
126 spacingZ = GetDistance(loader.GetSlice(0), loader.GetSlice(1)); | |
127 } | |
128 else | |
129 { | |
130 // This is a volume with one single slice: Choose a dummy | |
131 // z-dimension for voxels | |
132 spacingZ = 1; | |
133 } | |
134 | |
135 for (size_t i = 1; i < loader.GetSliceCount(); i++) | |
136 { | |
137 if (!GeometryToolbox::IsNear(spacingZ, GetDistance(loader.GetSlice(i - 1), loader.GetSlice(i)), | |
138 0.001 /* this is expressed in mm */)) | |
139 { | |
140 LOG(ERROR) << "The distance between successive slices is not constant in a volume image"; | |
141 SlicedVolumeBase::NotifyGeometryError(); | |
142 return; | |
143 } | |
144 } | |
145 | |
146 unsigned int width = loader.GetSlice(0).GetWidth(); | |
147 unsigned int height = loader.GetSlice(0).GetHeight(); | |
148 Orthanc::PixelFormat format = loader.GetSlice(0).GetConverter().GetExpectedPixelFormat(); | |
149 LOG(INFO) << "Creating a volume image of size " << width << "x" << height | |
150 << "x" << loader.GetSliceCount() << " in " << Orthanc::EnumerationToString(format); | |
151 | |
152 image_.reset(new ImageBuffer3D(format, width, height, loader.GetSliceCount())); | |
153 image_->SetAxialGeometry(loader.GetSlice(0).GetGeometry()); | |
154 image_->SetVoxelDimensions(loader.GetSlice(0).GetPixelSpacingX(), | |
155 loader.GetSlice(0).GetPixelSpacingY(), spacingZ); | |
156 image_->Clear(); | |
157 | |
158 downloadStack_.reset(new DownloadStack(loader.GetSliceCount())); | |
159 | |
160 for (unsigned int i = 0; i < 4; i++) // Limit to 4 simultaneous downloads | |
161 { | |
162 ScheduleSliceDownload(); | |
163 } | |
164 | |
165 // TODO Check the DicomFrameConverter are constant | |
166 | |
167 SlicedVolumeBase::NotifyGeometryReady(); | |
168 } | |
169 | |
170 virtual void NotifyGeometryError(const OrthancSlicesLoader& loader) | |
171 { | |
172 LOG(ERROR) << "Unable to download a volume image"; | |
173 SlicedVolumeBase::NotifyGeometryError(); | |
174 } | |
175 | |
176 virtual void NotifySliceImageReady(const OrthancSlicesLoader& loader, | |
177 unsigned int sliceIndex, | |
178 const Slice& slice, | |
179 std::auto_ptr<Orthanc::ImageAccessor>& image, | |
180 SliceImageQuality quality) | |
181 { | |
182 { | |
183 ImageBuffer3D::SliceWriter writer(*image_, VolumeProjection_Axial, sliceIndex); | |
184 Orthanc::ImageProcessing::Copy(writer.GetAccessor(), *image); | |
185 } | |
186 | |
187 SlicedVolumeBase::NotifySliceChange(sliceIndex, slice); | |
188 | |
189 ScheduleSliceDownload(); | |
190 } | |
191 | |
192 virtual void NotifySliceImageError(const OrthancSlicesLoader& loader, | |
193 unsigned int sliceIndex, | |
194 const Slice& slice, | |
195 SliceImageQuality quality) | |
196 { | |
197 LOG(ERROR) << "Cannot download slice " << sliceIndex << " in a volume image"; | |
198 ScheduleSliceDownload(); | |
199 } | |
200 | |
201 public: | |
202 OrthancVolumeImage(IWebService& orthanc) : | |
203 loader_(*this, orthanc) | |
204 { | |
205 } | |
206 | |
207 void ScheduleLoadSeries(const std::string& seriesId) | |
208 { | |
209 loader_.ScheduleLoadSeries(seriesId); | |
210 } | |
211 | |
212 void ScheduleLoadInstance(const std::string& instanceId, | |
213 unsigned int frame) | |
214 { | |
215 loader_.ScheduleLoadInstance(instanceId, frame); | |
216 } | |
217 | |
218 virtual size_t GetSliceCount() const | |
219 { | |
220 return loader_.GetSliceCount(); | |
221 } | |
222 | |
223 virtual const Slice& GetSlice(size_t index) const | |
224 { | |
225 return loader_.GetSlice(index); | |
226 } | |
227 | |
228 ImageBuffer3D& GetImage() const | |
229 { | |
230 if (image_.get() == NULL) | |
231 { | |
232 // The geometry is not ready yet | |
233 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); | |
234 } | |
235 else | |
236 { | |
237 return *image_; | |
238 } | |
239 } | |
240 }; | |
241 | |
242 | |
243 class VolumeImageGeometry | |
244 { | |
245 private: | |
246 unsigned int width_; | |
247 unsigned int height_; | |
248 size_t depth_; | |
249 double pixelSpacingX_; | |
250 double pixelSpacingY_; | |
251 double sliceThickness_; | |
252 SliceGeometry reference_; | |
253 DicomFrameConverter converter_; | |
254 | |
255 double ComputeAxialThickness(const OrthancVolumeImage& volume) const | |
256 { | |
257 double thickness; | |
258 | |
259 size_t n = volume.GetSliceCount(); | |
260 if (n > 1) | |
261 { | |
262 const Slice& a = volume.GetSlice(0); | |
263 const Slice& b = volume.GetSlice(n - 1); | |
264 thickness = ((reference_.ProjectAlongNormal(b.GetGeometry().GetOrigin()) - | |
265 reference_.ProjectAlongNormal(a.GetGeometry().GetOrigin())) / | |
266 (static_cast<double>(n) - 1.0)); | |
267 } | |
268 else | |
269 { | |
270 thickness = volume.GetSlice(0).GetThickness(); | |
271 } | |
272 | |
273 if (thickness <= 0) | |
274 { | |
275 // The slices should have been sorted with increasing Z | |
276 // (along the normal) by the OrthancSlicesLoader | |
277 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); | |
278 } | |
279 else | |
280 { | |
281 return thickness; | |
282 } | |
283 } | |
284 | |
285 void SetupAxial(const OrthancVolumeImage& volume) | |
286 { | |
287 const Slice& axial = volume.GetSlice(0); | |
288 | |
289 width_ = axial.GetWidth(); | |
290 height_ = axial.GetHeight(); | |
291 depth_ = volume.GetSliceCount(); | |
292 | |
293 pixelSpacingX_ = axial.GetPixelSpacingX(); | |
294 pixelSpacingY_ = axial.GetPixelSpacingY(); | |
295 sliceThickness_ = ComputeAxialThickness(volume); | |
296 | |
297 reference_ = axial.GetGeometry(); | |
298 } | |
299 | |
300 void SetupCoronal(const OrthancVolumeImage& volume) | |
301 { | |
302 const Slice& axial = volume.GetSlice(0); | |
303 double axialThickness = ComputeAxialThickness(volume); | |
304 | |
305 width_ = axial.GetWidth(); | |
306 height_ = volume.GetSliceCount(); | |
307 depth_ = axial.GetHeight(); | |
308 | |
309 pixelSpacingX_ = axial.GetPixelSpacingX(); | |
310 pixelSpacingY_ = axialThickness; | |
311 sliceThickness_ = axial.GetPixelSpacingY(); | |
312 | |
313 Vector origin = axial.GetGeometry().GetOrigin(); | |
314 origin += (static_cast<double>(volume.GetSliceCount() - 1) * | |
315 axialThickness * axial.GetGeometry().GetNormal()); | |
316 | |
317 reference_ = SliceGeometry(origin, | |
318 axial.GetGeometry().GetAxisX(), | |
319 -axial.GetGeometry().GetNormal()); | |
320 } | |
321 | |
322 void SetupSagittal(const OrthancVolumeImage& volume) | |
323 { | |
324 const Slice& axial = volume.GetSlice(0); | |
325 double axialThickness = ComputeAxialThickness(volume); | |
326 | |
327 width_ = axial.GetHeight(); | |
328 height_ = volume.GetSliceCount(); | |
329 depth_ = axial.GetWidth(); | |
330 | |
331 pixelSpacingX_ = axial.GetPixelSpacingY(); | |
332 pixelSpacingY_ = axialThickness; | |
333 sliceThickness_ = axial.GetPixelSpacingX(); | |
334 | |
335 Vector origin = axial.GetGeometry().GetOrigin(); | |
336 origin += (static_cast<double>(volume.GetSliceCount() - 1) * | |
337 axialThickness * axial.GetGeometry().GetNormal()); | |
338 | |
339 reference_ = SliceGeometry(origin, | |
340 axial.GetGeometry().GetAxisY(), | |
341 axial.GetGeometry().GetNormal()); | |
342 } | |
343 | |
344 public: | |
345 VolumeImageGeometry(const OrthancVolumeImage& volume, | |
346 VolumeProjection projection) | |
347 { | |
348 if (volume.GetSliceCount() == 0) | |
349 { | |
350 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); | |
351 } | |
352 | |
353 converter_ = volume.GetSlice(0).GetConverter(); | |
354 | |
355 switch (projection) | |
356 { | |
357 case VolumeProjection_Axial: | |
358 SetupAxial(volume); | |
359 break; | |
360 | |
361 case VolumeProjection_Coronal: | |
362 SetupCoronal(volume); | |
363 break; | |
364 | |
365 case VolumeProjection_Sagittal: | |
366 SetupSagittal(volume); | |
367 break; | |
368 | |
369 default: | |
370 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); | |
371 } | |
372 } | |
373 | |
374 size_t GetSliceCount() const | |
375 { | |
376 return depth_; | |
377 } | |
378 | |
379 const Vector& GetNormal() const | |
380 { | |
381 return reference_.GetNormal(); | |
382 } | |
383 | |
384 bool LookupSlice(size_t& index, | |
385 const SliceGeometry& slice) const | |
386 { | |
387 bool opposite; | |
388 if (!GeometryToolbox::IsParallelOrOpposite(opposite, | |
389 reference_.GetNormal(), | |
390 slice.GetNormal())) | |
391 { | |
392 return false; | |
393 } | |
394 | |
395 double z = (reference_.ProjectAlongNormal(slice.GetOrigin()) - | |
396 reference_.ProjectAlongNormal(reference_.GetOrigin())) / sliceThickness_; | |
397 | |
398 int s = static_cast<int>(boost::math::iround(z)); | |
399 | |
400 if (s < 0 || | |
401 s >= static_cast<int>(depth_)) | |
402 { | |
403 return false; | |
404 } | |
405 else | |
406 { | |
407 index = static_cast<size_t>(s); | |
408 return true; | |
409 } | |
410 } | |
411 | |
412 Slice GetSlice(size_t slice) const | |
413 { | |
414 if (slice < 0 || | |
415 slice >= depth_) | |
416 { | |
417 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); | |
418 } | |
419 else | |
420 { | |
421 SliceGeometry origin(reference_.GetOrigin() + | |
422 static_cast<double>(slice) * sliceThickness_ * reference_.GetNormal(), | |
423 reference_.GetAxisX(), | |
424 reference_.GetAxisY()); | |
425 | |
426 return Slice(origin, pixelSpacingX_, pixelSpacingY_, sliceThickness_, | |
427 width_, height_, converter_); | |
428 } | |
429 } | |
430 }; | |
431 | |
432 | |
433 | |
434 class VolumeImageSource : | |
435 public LayerSourceBase, | |
436 private ISlicedVolume::IObserver | |
437 { | |
438 private: | |
439 OrthancVolumeImage& volume_; | |
440 std::auto_ptr<VolumeImageGeometry> axialGeometry_; | |
441 std::auto_ptr<VolumeImageGeometry> coronalGeometry_; | |
442 std::auto_ptr<VolumeImageGeometry> sagittalGeometry_; | |
443 | |
444 | |
445 bool IsGeometryReady() const | |
446 { | |
447 return axialGeometry_.get() != NULL; | |
448 } | |
449 | |
450 | |
451 virtual void NotifyGeometryReady(const ISlicedVolume& volume) | |
452 { | |
453 // These 3 values are only used to speed up the ILayerSource | |
454 axialGeometry_.reset(new VolumeImageGeometry(volume_, VolumeProjection_Axial)); | |
455 coronalGeometry_.reset(new VolumeImageGeometry(volume_, VolumeProjection_Coronal)); | |
456 sagittalGeometry_.reset(new VolumeImageGeometry(volume_, VolumeProjection_Sagittal)); | |
457 | |
458 LayerSourceBase::NotifyGeometryReady(); | |
459 } | |
460 | |
461 virtual void NotifyGeometryError(const ISlicedVolume& volume) | |
462 { | |
463 LayerSourceBase::NotifyGeometryError(); | |
464 } | |
465 | |
466 virtual void NotifyContentChange(const ISlicedVolume& volume) | |
467 { | |
468 LayerSourceBase::NotifyContentChange(); | |
469 } | |
470 | |
471 virtual void NotifySliceChange(const ISlicedVolume& volume, | |
472 const size_t& sliceIndex, | |
473 const Slice& slice) | |
474 { | |
475 //LayerSourceBase::NotifySliceChange(slice); | |
476 | |
477 // TODO Improve this? | |
478 LayerSourceBase::NotifyContentChange(); | |
479 } | |
480 | |
481 | |
482 const VolumeImageGeometry& GetProjectionGeometry(VolumeProjection projection) | |
483 { | |
484 if (!IsGeometryReady()) | |
485 { | |
486 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); | |
487 } | |
488 | |
489 switch (projection) | |
490 { | |
491 case VolumeProjection_Axial: | |
492 return *axialGeometry_; | |
493 | |
494 case VolumeProjection_Sagittal: | |
495 return *sagittalGeometry_; | |
496 | |
497 case VolumeProjection_Coronal: | |
498 return *coronalGeometry_; | |
499 | |
500 default: | |
501 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); | |
502 } | |
503 } | |
504 | |
505 | |
506 bool DetectProjection(VolumeProjection& projection, | |
507 const SliceGeometry& viewportSlice) | |
508 { | |
509 bool isOpposite; // Ignored | |
510 | |
511 if (GeometryToolbox::IsParallelOrOpposite(isOpposite, | |
512 viewportSlice.GetNormal(), | |
513 axialGeometry_->GetNormal())) | |
514 { | |
515 projection = VolumeProjection_Axial; | |
516 return true; | |
517 } | |
518 else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, | |
519 viewportSlice.GetNormal(), | |
520 sagittalGeometry_->GetNormal())) | |
521 { | |
522 projection = VolumeProjection_Sagittal; | |
523 return true; | |
524 } | |
525 else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, | |
526 viewportSlice.GetNormal(), | |
527 coronalGeometry_->GetNormal())) | |
528 { | |
529 projection = VolumeProjection_Coronal; | |
530 return true; | |
531 } | |
532 else | |
533 { | |
534 return false; | |
535 } | |
536 } | |
537 | |
538 | |
539 public: | |
540 VolumeImageSource(OrthancVolumeImage& volume) : | |
541 volume_(volume) | |
542 { | |
543 volume_.Register(*this); | |
544 } | |
545 | |
546 virtual bool GetExtent(std::vector<Vector>& points, | |
547 const SliceGeometry& viewportSlice) | |
548 { | |
549 VolumeProjection projection; | |
550 | |
551 if (!IsGeometryReady() || | |
552 !DetectProjection(projection, viewportSlice)) | |
553 { | |
554 return false; | |
555 } | |
556 else | |
557 { | |
558 // As the slices of the volumic image are arranged in a box, | |
559 // we only consider one single reference slice (the one with index 0). | |
560 GetProjectionGeometry(projection).GetSlice(0).GetExtent(points); | |
561 | |
562 return true; | |
563 } | |
564 } | |
565 | |
566 | |
567 virtual void ScheduleLayerCreation(const SliceGeometry& viewportSlice) | |
568 { | |
569 VolumeProjection projection; | |
570 | |
571 if (IsGeometryReady() && | |
572 DetectProjection(projection, viewportSlice)) | |
573 { | |
574 const VolumeImageGeometry& geometry = GetProjectionGeometry(projection); | |
575 | |
576 size_t closest; | |
577 | |
578 if (geometry.LookupSlice(closest, viewportSlice)) | |
579 { | |
580 bool isFullQuality = true; // TODO | |
581 | |
582 std::auto_ptr<Orthanc::Image> frame; | |
583 | |
584 { | |
585 ImageBuffer3D::SliceReader reader(volume_.GetImage(), projection, closest); | |
586 | |
587 // TODO Transfer ownership if non-axial, to avoid memcpy | |
588 frame.reset(Orthanc::Image::Clone(reader.GetAccessor())); | |
589 } | |
590 | |
591 Slice slice = geometry.GetSlice(closest); | |
592 LayerSourceBase::NotifyLayerReady( | |
593 FrameRenderer::CreateRenderer(frame.release(), slice, isFullQuality), | |
594 //new SliceOutlineRenderer(slice), | |
595 slice, false); | |
596 return; | |
597 } | |
598 } | |
599 | |
600 // Error | |
601 Slice slice; | |
602 LayerSourceBase::NotifyLayerReady(NULL, slice, true); | |
603 } | |
604 }; | |
605 } |