Mercurial > hg > orthanc-stone
comparison Framework/Volumes/VolumeReslicer.cpp @ 177:83200c4d07ca wasm
fix interpolation
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Fri, 09 Mar 2018 12:32:01 +0100 |
parents | 316324f42848 |
children | db21c1810c89 |
comparison
equal
deleted
inserted
replaced
176:ab9c799f5de1 | 177:83200c4d07ca |
---|---|
232 } | 232 } |
233 }; | 233 }; |
234 | 234 |
235 | 235 |
236 | 236 |
237 class VoxelReaderBase | 237 class VoxelReaderBase : public boost::noncopyable |
238 { | 238 { |
239 protected: | 239 private: |
240 const ImageBuffer3D& source_; | 240 const ImageBuffer3D& image_; |
241 unsigned int sourceWidth_; | 241 unsigned int width_; |
242 unsigned int sourceHeight_; | 242 unsigned int height_; |
243 unsigned int sourceDepth_; | 243 unsigned int depth_; |
244 | 244 float widthFloat_; |
245 public: | 245 float heightFloat_; |
246 VoxelReaderBase(const ImageBuffer3D& source) : | 246 float depthFloat_; |
247 source_(source), | 247 |
248 sourceWidth_(source.GetWidth()), | 248 public: |
249 sourceHeight_(source.GetHeight()), | 249 VoxelReaderBase(const ImageBuffer3D& image) : |
250 sourceDepth_(source.GetDepth()) | 250 image_(image), |
251 { | 251 width_(image.GetWidth()), |
252 } | 252 height_(image.GetHeight()), |
253 | 253 depth_(image.GetDepth()), |
254 unsigned int GetSourceWidth() const | 254 widthFloat_(static_cast<float>(image.GetWidth())), |
255 { | 255 heightFloat_(static_cast<float>(image.GetHeight())), |
256 return sourceWidth_; | 256 depthFloat_(static_cast<float>(image.GetDepth())) |
257 } | 257 { |
258 | 258 } |
259 unsigned int GetSourceHeight() const | 259 |
260 { | 260 const ImageBuffer3D& GetImage() const |
261 return sourceHeight_; | 261 { |
262 } | 262 return image_; |
263 | 263 } |
264 unsigned int GetSourceDepth() const | 264 |
265 { | 265 const unsigned int GetImageWidth() const |
266 return sourceDepth_; | 266 { |
267 } | 267 return width_; |
268 | 268 } |
269 bool GetNearestCoordinates(unsigned int& sourceX, | 269 |
270 unsigned int& sourceY, | 270 const unsigned int GetImageHeight() const |
271 unsigned int& sourceZ, | 271 { |
272 float worldX, | 272 return height_; |
273 float worldY, | 273 } |
274 float worldZ) const | 274 |
275 { | 275 const unsigned int GetImageDepth() const |
276 if (worldX >= 0 && | 276 { |
277 worldY >= 0 && | 277 return depth_; |
278 worldZ >= 0) | 278 } |
279 { | 279 |
280 sourceX = static_cast<unsigned int>(worldX * static_cast<float>(sourceWidth_)); | 280 bool GetNearestCoordinates(unsigned int& imageX, |
281 sourceY = static_cast<unsigned int>(worldY * static_cast<float>(sourceHeight_)); | 281 unsigned int& imageY, |
282 sourceZ = static_cast<unsigned int>(worldZ * static_cast<float>(sourceDepth_)); | 282 unsigned int& imageZ, |
283 | 283 float& fractionalX, |
284 return (sourceX < sourceWidth_ && | 284 float& fractionalY, |
285 sourceY < sourceHeight_ && | 285 float& fractionalZ, |
286 sourceZ < sourceDepth_); | 286 float volumeX, |
287 float volumeY, | |
288 float volumeZ) const | |
289 { | |
290 if (volumeX >= 0 && | |
291 volumeY >= 0 && | |
292 volumeZ >= 0) | |
293 { | |
294 const float x = volumeX * widthFloat_; | |
295 const float y = volumeY * heightFloat_; | |
296 const float z = volumeZ * depthFloat_; | |
297 | |
298 imageX = static_cast<unsigned int>(std::floor(x)); | |
299 imageY = static_cast<unsigned int>(std::floor(y)); | |
300 imageZ = static_cast<unsigned int>(std::floor(z)); | |
301 | |
302 if (imageX < width_ && | |
303 imageY < height_ && | |
304 imageZ < depth_) | |
305 { | |
306 fractionalX = x - static_cast<float>(imageX); | |
307 fractionalY = y - static_cast<float>(imageY); | |
308 fractionalZ = z - static_cast<float>(imageZ); | |
309 return true; | |
310 } | |
311 else | |
312 { | |
313 return false; | |
314 } | |
287 } | 315 } |
288 else | 316 else |
289 { | 317 { |
290 return false; | 318 return false; |
291 } | 319 } |
303 public VoxelReaderBase | 331 public VoxelReaderBase |
304 { | 332 { |
305 public: | 333 public: |
306 typedef typename PixelTraits<InputFormat>::PixelType InputPixelType; | 334 typedef typename PixelTraits<InputFormat>::PixelType InputPixelType; |
307 | 335 |
308 VoxelReader(const ImageBuffer3D& source) : | 336 VoxelReader(const ImageBuffer3D& image) : |
309 VoxelReaderBase(source) | 337 VoxelReaderBase(image) |
310 { | 338 { |
311 } | 339 } |
312 | 340 |
313 ORTHANC_STONE_FORCE_INLINE | 341 ORTHANC_STONE_FORCE_INLINE |
314 float GetFloatValue(float worldX, | 342 float GetFloatValue(float volumeX, |
315 float worldY, | 343 float volumeY, |
316 float worldZ) const | 344 float volumeZ) const |
317 { | 345 { |
318 InputPixelType value; | 346 InputPixelType value; |
319 GetValue(value, worldX, worldY, worldZ); | 347 GetValue(value, volumeX, volumeY, volumeZ); |
320 return static_cast<float>(value); | 348 return static_cast<float>(value); |
321 } | 349 } |
322 | 350 |
323 ORTHANC_STONE_FORCE_INLINE | 351 ORTHANC_STONE_FORCE_INLINE |
324 void GetValue(InputPixelType& target, | 352 void GetValue(InputPixelType& target, |
325 float worldX, | 353 float volumeX, |
326 float worldY, | 354 float volumeY, |
327 float worldZ) const | 355 float volumeZ) const |
328 { | 356 { |
329 unsigned int sourceX, sourceY, sourceZ; | 357 unsigned int imageX, imageY, imageZ; |
330 | 358 float fractionalX, fractionalY, fractionalZ; // unused |
331 if (GetNearestCoordinates(sourceX, sourceY, sourceZ, worldX, worldY, worldZ)) | 359 |
332 { | 360 if (GetNearestCoordinates(imageX, imageY, imageZ, |
333 PixelTraits<InputFormat>::GetVoxel(target, source_, sourceX, sourceY, sourceZ); | 361 fractionalX, fractionalY, fractionalZ, |
362 volumeX, volumeY, volumeZ)) | |
363 { | |
364 PixelTraits<InputFormat>::GetVoxel(target, GetImage(), imageX, imageY, imageZ); | |
334 } | 365 } |
335 else | 366 else |
336 { | 367 { |
337 PixelTraits<InputFormat>::SetOutOfVolume(target); | 368 PixelTraits<InputFormat>::SetOutOfVolume(target); |
338 } | 369 } |
346 { | 377 { |
347 private: | 378 private: |
348 float outOfVolume_; | 379 float outOfVolume_; |
349 | 380 |
350 public: | 381 public: |
351 VoxelReader(const ImageBuffer3D& source) : | 382 VoxelReader(const ImageBuffer3D& image) : |
352 VoxelReaderBase(source) | 383 VoxelReaderBase(image) |
353 { | 384 { |
354 typename PixelTraits<InputFormat>::PixelType value; | 385 typename PixelTraits<InputFormat>::PixelType value; |
355 PixelTraits<InputFormat>::SetOutOfVolume(value); | 386 PixelTraits<InputFormat>::SetOutOfVolume(value); |
356 outOfVolume_ = static_cast<float>(value); | 387 outOfVolume_ = static_cast<float>(value); |
357 } | 388 } |
358 | 389 |
359 void SampleVoxels(float& f00, | 390 void SampleVoxels(float& f00, |
391 float& f01, | |
360 float& f10, | 392 float& f10, |
361 float& f01, | |
362 float& f11, | 393 float& f11, |
363 unsigned int sourceX, | 394 unsigned int imageX, |
364 unsigned int sourceY, | 395 unsigned int imageY, |
365 unsigned int sourceZ) const | 396 unsigned int imageZ) const |
366 { | 397 { |
367 f00 = PixelTraits<InputFormat>::GetFloatVoxel(source_, sourceX, sourceY, sourceZ); | 398 f00 = PixelTraits<InputFormat>::GetFloatVoxel(GetImage(), imageX, imageY, imageZ); |
368 | 399 |
369 if (sourceX + 1 < sourceWidth_) | 400 if (imageX + 1 < GetImageWidth()) |
370 { | 401 { |
371 f01 = PixelTraits<InputFormat>::GetFloatVoxel(source_, sourceX + 1, sourceY, sourceZ); | 402 f01 = PixelTraits<InputFormat>::GetFloatVoxel(GetImage(), imageX + 1, imageY, imageZ); |
372 } | 403 } |
373 else | 404 else |
374 { | 405 { |
375 f01 = f00; | 406 f01 = f00; |
376 } | 407 } |
377 | 408 |
378 if (sourceY + 1 < sourceHeight_) | 409 if (imageY + 1 < GetImageWidth()) |
379 { | 410 { |
380 f10 = PixelTraits<InputFormat>::GetFloatVoxel(source_, sourceX, sourceY + 1, sourceZ); | 411 f10 = PixelTraits<InputFormat>::GetFloatVoxel(GetImage(), imageX, imageY + 1, imageZ); |
381 } | 412 } |
382 else | 413 else |
383 { | 414 { |
384 f10 = f00; | 415 f10 = f00; |
385 } | 416 } |
386 | 417 |
387 if (sourceX + 1 < sourceWidth_ && | 418 if (imageX + 1 < GetImageWidth() && |
388 sourceY + 1 < sourceHeight_) | 419 imageY + 1 < GetImageHeight()) |
389 { | 420 { |
390 f11 = PixelTraits<InputFormat>::GetFloatVoxel(source_, sourceX + 1, sourceY + 1, sourceZ); | 421 f11 = PixelTraits<InputFormat>::GetFloatVoxel(GetImage(), imageX + 1, imageY + 1, imageZ); |
391 } | 422 } |
392 else | 423 else |
393 { | 424 { |
394 f11 = f00; | 425 f11 = f00; |
395 } | 426 } |
398 float GetOutOfVolume() const | 429 float GetOutOfVolume() const |
399 { | 430 { |
400 return outOfVolume_; | 431 return outOfVolume_; |
401 } | 432 } |
402 | 433 |
403 float GetFloatValue(float worldX, | 434 float GetFloatValue(float volumeX, |
404 float worldY, | 435 float volumeY, |
405 float worldZ) const | 436 float volumeZ) const |
406 { | 437 { |
407 unsigned int sourceX, sourceY, sourceZ; | 438 unsigned int imageX, imageY, imageZ; |
408 | 439 float fractionalX, fractionalY, fractionalZ; |
409 if (GetNearestCoordinates(sourceX, sourceY, sourceZ, worldX, worldY, worldZ)) | 440 |
410 { | 441 if (GetNearestCoordinates(imageX, imageY, imageZ, |
411 float f00, f10, f01, f11; | 442 fractionalX, fractionalY, fractionalZ, |
412 SampleVoxels(f00, f10, f01, f11, sourceX, sourceY, sourceZ); | 443 volumeX, volumeY, volumeZ)) |
413 return GeometryToolbox::ComputeBilinearInterpolation(worldX, worldY, f00, f01, f10, f11); | 444 { |
445 float f00, f01, f10, f11; | |
446 SampleVoxels(f00, f01, f10, f11, imageX, imageY, imageZ); | |
447 return GeometryToolbox::ComputeBilinearInterpolationUnitSquare | |
448 (fractionalX, fractionalY, f00, f01, f10, f11); | |
414 } | 449 } |
415 else | 450 else |
416 { | 451 { |
417 return outOfVolume_; | 452 return outOfVolume_; |
418 } | 453 } |
419 } | 454 } |
420 }; | 455 }; |
421 | 456 |
422 | 457 |
423 template <Orthanc::PixelFormat InputFormat> | 458 template <Orthanc::PixelFormat InputFormat> |
424 class VoxelReader<InputFormat, ImageInterpolation_Trilinear> | 459 class VoxelReader<InputFormat, ImageInterpolation_Trilinear> : |
460 public VoxelReaderBase | |
425 { | 461 { |
426 private: | 462 private: |
427 typedef VoxelReader<InputFormat, ImageInterpolation_Bilinear> Bilinear; | 463 typedef VoxelReader<InputFormat, ImageInterpolation_Bilinear> Bilinear; |
428 | 464 |
429 Bilinear bilinear_; | 465 Bilinear bilinear_; |
430 unsigned int sourceDepth_; | 466 unsigned int imageDepth_; |
431 | 467 |
432 public: | 468 public: |
433 VoxelReader(const ImageBuffer3D& source) : | 469 VoxelReader(const ImageBuffer3D& image) : |
434 bilinear_(source), | 470 VoxelReaderBase(image), |
435 sourceDepth_(source.GetDepth()) | 471 bilinear_(image), |
436 { | 472 imageDepth_(image.GetDepth()) |
437 } | 473 { |
438 | 474 } |
439 float GetFloatValue(float worldX, | 475 |
440 float worldY, | 476 float GetFloatValue(float volumeX, |
441 float worldZ) const | 477 float volumeY, |
442 { | 478 float volumeZ) const |
443 unsigned int sourceX, sourceY, sourceZ; | 479 { |
444 | 480 unsigned int imageX, imageY, imageZ; |
445 if (bilinear_.GetNearestCoordinates(sourceX, sourceY, sourceZ, worldX, worldY, worldZ)) | 481 float fractionalX, fractionalY, fractionalZ; |
446 { | 482 |
447 float f000, f010, f001, f011, f100, f110, f101, f111; | 483 if (GetNearestCoordinates(imageX, imageY, imageZ, |
448 | 484 fractionalX, fractionalY, fractionalZ, |
449 bilinear_.SampleVoxels(f000, f010, f001, f011, sourceX, sourceY, sourceZ); | 485 volumeX, volumeY, volumeZ)) |
450 | 486 { |
451 if (sourceZ + 1 < sourceDepth_) | 487 float f000, f001, f010, f011; |
488 bilinear_.SampleVoxels(f000, f001, f010, f011, imageX, imageY, imageZ); | |
489 | |
490 if (imageZ + 1 < imageDepth_) | |
452 { | 491 { |
453 bilinear_.SampleVoxels(f100, f110, f101, f111, sourceX, sourceY, sourceZ + 1); | 492 float f100, f101, f110, f111; |
454 return GeometryToolbox::ComputeTrilinearInterpolation | 493 bilinear_.SampleVoxels(f100, f101, f110, f111, imageX, imageY, imageZ + 1); |
455 (worldX, worldY, worldZ, f000, f001, f010, f011, f100, f110, f101, f111); | 494 return GeometryToolbox::ComputeTrilinearInterpolationUnitSquare |
495 (fractionalX, fractionalY, fractionalZ, | |
496 f000, f001, f010, f011, f100, f101, f110, f111); | |
456 } | 497 } |
457 else | 498 else |
458 { | 499 { |
459 return GeometryToolbox::ComputeBilinearInterpolation(worldX, worldY, f000, f001, f010, f011); | 500 return GeometryToolbox::ComputeBilinearInterpolationUnitSquare |
501 (fractionalX, fractionalY, f000, f001, f010, f011); | |
460 } | 502 } |
461 } | 503 } |
462 else | 504 else |
463 { | 505 { |
464 return bilinear_.GetOutOfVolume(); | 506 return bilinear_.GetOutOfVolume(); |
489 { | 531 { |
490 } | 532 } |
491 | 533 |
492 ORTHANC_STONE_FORCE_INLINE | 534 ORTHANC_STONE_FORCE_INLINE |
493 void Apply(typename PixelWriter::OutputPixelType* pixel, | 535 void Apply(typename PixelWriter::OutputPixelType* pixel, |
494 float worldX, | 536 float volumeX, |
495 float worldY, | 537 float volumeY, |
496 float worldZ) | 538 float volumeZ) |
497 { | 539 { |
498 typename VoxelReader::InputPixelType source; | 540 typename VoxelReader::InputPixelType image; |
499 reader_.GetValue(source, worldX, worldY, worldZ); | 541 reader_.GetValue(image, volumeX, volumeY, volumeZ); |
500 writer_.SetValue(pixel, source); | 542 writer_.SetValue(pixel, image); |
501 } | 543 } |
502 }; | 544 }; |
503 | 545 |
504 | 546 |
505 template <typename VoxelReader, | 547 template <typename VoxelReader, |
518 { | 560 { |
519 } | 561 } |
520 | 562 |
521 ORTHANC_STONE_FORCE_INLINE | 563 ORTHANC_STONE_FORCE_INLINE |
522 void Apply(typename PixelWriter::OutputPixelType* pixel, | 564 void Apply(typename PixelWriter::OutputPixelType* pixel, |
523 float worldX, | 565 float volumeX, |
524 float worldY, | 566 float volumeY, |
525 float worldZ) | 567 float volumeZ) |
526 { | 568 { |
527 writer_.SetFloatValue(pixel, reader_.GetFloatValue(worldX, worldY, worldZ)); | 569 writer_.SetFloatValue(pixel, reader_.GetFloatValue(volumeX, volumeY, volumeZ)); |
528 } | 570 } |
529 }; | 571 }; |
530 | 572 |
531 | 573 |
532 template <typename VoxelReader, | 574 template <typename VoxelReader, |
549 { | 591 { |
550 } | 592 } |
551 | 593 |
552 ORTHANC_STONE_FORCE_INLINE | 594 ORTHANC_STONE_FORCE_INLINE |
553 void Apply(typename PixelWriter::OutputPixelType* pixel, | 595 void Apply(typename PixelWriter::OutputPixelType* pixel, |
554 float worldX, | 596 float volumeX, |
555 float worldY, | 597 float volumeY, |
556 float worldZ) | 598 float volumeZ) |
557 { | 599 { |
558 writer_.SetFloatValue(pixel, scaling_ * reader_.GetFloatValue(worldX, worldY, worldZ) + offset_); | 600 writer_.SetFloatValue(pixel, scaling_ * reader_.GetFloatValue(volumeX, volumeY, volumeZ) + offset_); |
559 } | 601 } |
560 }; | 602 }; |
561 | 603 |
562 | 604 |
563 | 605 |
607 position_[1] += offset_[1]; | 649 position_[1] += offset_[1]; |
608 position_[2] += offset_[2]; | 650 position_[2] += offset_[2]; |
609 } | 651 } |
610 | 652 |
611 ORTHANC_STONE_FORCE_INLINE | 653 ORTHANC_STONE_FORCE_INLINE |
612 void GetNearestCoordinates(float& x, | 654 void GetVolumeCoordinates(float& x, |
613 float& y, | 655 float& y, |
614 float& z) const | 656 float& z) const |
615 { | 657 { |
616 x = position_[0]; | 658 x = position_[0]; |
617 y = position_[1]; | 659 y = position_[1]; |
618 z = position_[2]; | 660 z = position_[2]; |
619 } | 661 } |
649 void Next() | 691 void Next() |
650 { | 692 { |
651 x_++; | 693 x_++; |
652 } | 694 } |
653 | 695 |
654 void GetNearestCoordinates(float& x, | 696 void GetVolumeCoordinates(float& x, |
655 float& y, | 697 float& y, |
656 float& z) const | 698 float& z) const |
657 { | 699 { |
658 assert(x_ < slice_.GetWidth()); | 700 assert(x_ < slice_.GetWidth()); |
659 | 701 |
660 const double width = static_cast<double>(slice_.GetWidth()); | 702 const double width = static_cast<double>(slice_.GetWidth()); |
661 const double height = static_cast<double>(slice_.GetHeight()); | 703 const double height = static_cast<double>(slice_.GetHeight()); |
702 | 744 |
703 RowIterator it(slice, extent, plane, box, y); | 745 RowIterator it(slice, extent, plane, box, y); |
704 | 746 |
705 for (unsigned int x = 0; x < outputWidth; x++, p++) | 747 for (unsigned int x = 0; x < outputWidth; x++, p++) |
706 { | 748 { |
707 float worldX, worldY, worldZ; | 749 float volumeX, volumeY, volumeZ; |
708 it.GetNearestCoordinates(worldX, worldY, worldZ); | 750 it.GetVolumeCoordinates(volumeX, volumeY, volumeZ); |
709 shader.Apply(p, worldX, worldY, worldZ); | 751 shader.Apply(p, volumeX, volumeY, volumeZ); |
710 it.Next(); | 752 it.Next(); |
711 } | 753 } |
712 } | 754 } |
713 } | 755 } |
714 | 756 |
843 SlowRowIterator slow(*slice_, extent_, plane, box, y); | 885 SlowRowIterator slow(*slice_, extent_, plane, box, y); |
844 | 886 |
845 for (unsigned int x = 0; x < slice_->GetWidth(); x++) | 887 for (unsigned int x = 0; x < slice_->GetWidth(); x++) |
846 { | 888 { |
847 float px, py, pz; | 889 float px, py, pz; |
848 fast.GetNearestCoordinates(px, py, pz); | 890 fast.GetVolumeCoordinates(px, py, pz); |
849 | 891 |
850 float qx, qy, qz; | 892 float qx, qy, qz; |
851 slow.GetNearestCoordinates(qx, qy, qz); | 893 slow.GetVolumeCoordinates(qx, qy, qz); |
852 | 894 |
853 Vector d; | 895 Vector d; |
854 LinearAlgebra::AssignVector(d, px - qx, py - qy, pz - qz); | 896 LinearAlgebra::AssignVector(d, px - qx, py - qy, pz - qz); |
855 double norm = boost::numeric::ublas::norm_2(d); | 897 double norm = boost::numeric::ublas::norm_2(d); |
856 if (norm > 0.0001) | 898 if (norm > 0.0001) |