Mercurial > hg > orthanc-stone
comparison Framework/Toolbox/ShearWarpProjectiveTransform.cpp @ 193:4abddd083374 wasm
ShearWarpProjectiveTransform::ApplyAxial()
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Tue, 20 Mar 2018 14:05:39 +0100 |
parents | 46cb2eedc2e0 |
children | dabe9982fca3 |
comparison
equal
deleted
inserted
replaced
192:371da7fe2c0e | 193:4abddd083374 |
---|---|
24 #include "ImageGeometry.h" | 24 #include "ImageGeometry.h" |
25 #include "Extent2D.h" | 25 #include "Extent2D.h" |
26 #include "FiniteProjectiveCamera.h" | 26 #include "FiniteProjectiveCamera.h" |
27 #include "GeometryToolbox.h" | 27 #include "GeometryToolbox.h" |
28 | 28 |
29 #include <Core/Images/PixelTraits.h> | |
30 #include <Core/Images/ImageProcessing.h> | |
29 #include <Core/OrthancException.h> | 31 #include <Core/OrthancException.h> |
30 #include <Core/Logging.h> | 32 #include <Core/Logging.h> |
31 | 33 |
32 #include <boost/numeric/ublas/matrix_proxy.hpp> | 34 #include <boost/numeric/ublas/matrix_proxy.hpp> |
35 #include <boost/math/special_functions/round.hpp> | |
33 #include <cassert> | 36 #include <cassert> |
34 | 37 |
35 | 38 |
36 namespace OrthancStone | 39 namespace OrthancStone |
37 { | 40 { |
324 } | 327 } |
325 } | 328 } |
326 | 329 |
327 return M_view; | 330 return M_view; |
328 } | 331 } |
332 | |
333 | |
334 template <Orthanc::PixelFormat SourceFormat, | |
335 Orthanc::PixelFormat TargetFormat, | |
336 bool MIP> | |
337 static void ApplyAxialInternal(Orthanc::ImageAccessor& target, | |
338 float& maxValue, | |
339 const Matrix& M_view, | |
340 const ImageBuffer3D& source, | |
341 double pixelSpacing, | |
342 unsigned int countSlices, | |
343 ImageInterpolation shearInterpolation, | |
344 ImageInterpolation warpInterpolation) | |
345 { | |
346 typedef Orthanc::PixelTraits<SourceFormat> SourceTraits; | |
347 typedef Orthanc::PixelTraits<TargetFormat> TargetTraits; | |
348 | |
349 /** | |
350 * Step 1: Precompute some information. | |
351 **/ | |
352 | |
353 if (target.GetFormat() != TargetFormat || | |
354 source.GetFormat() != SourceFormat || | |
355 !std::numeric_limits<float>::is_iec559 || | |
356 sizeof(float) != 4) | |
357 { | |
358 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); | |
359 } | |
360 | |
361 if (countSlices > source.GetDepth()) | |
362 { | |
363 countSlices = source.GetDepth(); | |
364 } | |
365 | |
366 if (countSlices == 0) | |
367 { | |
368 maxValue = 0; | |
369 Orthanc::ImageProcessing::Set(target, 0); | |
370 return; | |
371 } | |
372 | |
373 LOG(INFO) << "Number of rendered slices: " << countSlices; | |
374 | |
375 | |
376 /** | |
377 * Step 2: Extract the shear-warp transform corresponding to | |
378 * M_view. | |
379 **/ | |
380 | |
381 // Compute the "world" matrix that maps the source volume to the | |
382 // (0,0,0)->(1,1,1) unit cube | |
383 Vector origin = source.GetCoordinates(0, 0, 0); | |
384 Vector ps = source.GetVoxelDimensions(VolumeProjection_Axial); | |
385 Matrix world = LinearAlgebra::Product( | |
386 GeometryToolbox::CreateScalingMatrix(1.0 / ps[0], 1.0 / ps[1], 1.0 / ps[2]), | |
387 GeometryToolbox::CreateTranslationMatrix(-origin[0], -origin[1], -origin[2])); | |
388 | |
389 Matrix worldInv; | |
390 LinearAlgebra::InvertMatrix(worldInv, world); | |
391 | |
392 ShearWarpProjectiveTransform shearWarp(LinearAlgebra::Product(M_view, worldInv), | |
393 /*LinearAlgebra::IdentityMatrix(4),*/ | |
394 source.GetWidth(), | |
395 source.GetHeight(), | |
396 source.GetDepth(), | |
397 pixelSpacing, pixelSpacing, | |
398 target.GetWidth(), target.GetHeight()); | |
399 | |
400 const unsigned int intermediateWidth = shearWarp.GetIntermediateWidth(); | |
401 const unsigned int intermediateHeight = shearWarp.GetIntermediateHeight(); | |
402 | |
403 | |
404 /** | |
405 * Step 3: Apply the "shear" part of the transform to form the | |
406 * intermediate image. The sheared images are accumulated into the | |
407 * Float32 image "accumulator". The number of samples available | |
408 * for each pixel is stored in the "counter" image. | |
409 **/ | |
410 | |
411 std::auto_ptr<Orthanc::ImageAccessor> accumulator, counter, intermediate; | |
412 | |
413 accumulator.reset(new Orthanc::Image(Orthanc::PixelFormat_Float32, | |
414 intermediateWidth, intermediateHeight, false)); | |
415 counter.reset(new Orthanc::Image(Orthanc::PixelFormat_Grayscale16, | |
416 intermediateWidth, intermediateHeight, false)); | |
417 intermediate.reset(new Orthanc::Image(SourceFormat, intermediateWidth, intermediateHeight, false)); | |
418 | |
419 Orthanc::ImageProcessing::Set(*accumulator, 0); | |
420 Orthanc::ImageProcessing::Set(*counter, 0); | |
421 | |
422 // Loop around the slices of the volume | |
423 for (unsigned int i = 0; i <= countSlices; i++) | |
424 { | |
425 // (3.a) Compute the shear for this specific slice | |
426 unsigned int z = static_cast<unsigned int>( | |
427 boost::math::iround(static_cast<double>(i) / | |
428 static_cast<double>(countSlices) * | |
429 static_cast<double>(source.GetDepth() - 1))); | |
430 | |
431 double a11, b1, a22, b2, vz; | |
432 shearWarp.ComputeShearOnSlice(a11, b1, a22, b2, vz, static_cast<double>(z) + 0.5); | |
433 | |
434 | |
435 { | |
436 // (3.b) Detect the "useful" portion of the intermediate image | |
437 // for this slice (i.e. the bounding box where the source | |
438 // slice is mapped to by the shear), so as to update "counter" | |
439 Matrix a = LinearAlgebra::ZeroMatrix(3, 3); | |
440 a(0,0) = a11; | |
441 a(0,2) = b1; | |
442 a(1,1) = a22; | |
443 a(1,2) = b2; | |
444 a(2,2) = 1; | |
445 | |
446 unsigned int x1, y1, x2, y2; | |
447 if (GetProjectiveTransformExtent(x1, y1, x2, y2, a, | |
448 source.GetWidth(), source.GetHeight(), | |
449 intermediateWidth, intermediateHeight)) | |
450 { | |
451 for (unsigned int y = y1; y <= y2; y++) | |
452 { | |
453 uint16_t* p = reinterpret_cast<uint16_t*>(counter->GetRow(y)) + x1; | |
454 for (unsigned int x = x1; x <= x2; x++, p++) | |
455 { | |
456 if (MIP) | |
457 { | |
458 // TODO - In the case of MIP, "counter" could be | |
459 // reduced to "PixelFormat_Grayscale8" to reduce | |
460 // memory usage | |
461 *p = 1; | |
462 } | |
463 else | |
464 { | |
465 *p += 1; | |
466 } | |
467 } | |
468 } | |
469 } | |
470 } | |
471 | |
472 | |
473 { | |
474 // (3.c) Shear the source slice into a temporary image | |
475 ImageBuffer3D::SliceReader reader(source, VolumeProjection_Axial, z); | |
476 ApplyAffineTransform(*intermediate, reader.GetAccessor(), | |
477 a11, 0, b1, | |
478 0, a22, b2, | |
479 shearInterpolation); | |
480 } | |
481 | |
482 | |
483 for (unsigned int y = 0; y < intermediateHeight; y++) | |
484 { | |
485 // (3.d) Accumulate the pixels of the sheared image into "accumulator" | |
486 const typename SourceTraits::PixelType* p = | |
487 reinterpret_cast<const typename SourceTraits::PixelType*>(intermediate->GetConstRow(y)); | |
488 | |
489 float* q = reinterpret_cast<float*>(accumulator->GetRow(y)); | |
490 | |
491 for (unsigned int x = 0; x < intermediateWidth; x++) | |
492 { | |
493 float pixel = SourceTraits::PixelToFloat(*p); | |
494 | |
495 if (MIP) | |
496 { | |
497 // Get maximum for MIP | |
498 if (*q < pixel) | |
499 { | |
500 *q = pixel; | |
501 } | |
502 } | |
503 else | |
504 { | |
505 *q += pixel; | |
506 } | |
507 | |
508 p++; | |
509 q++; | |
510 } | |
511 } | |
512 } | |
513 | |
514 | |
515 /** | |
516 * Step 4: The intermediate image (that will be transformed by the | |
517 * "warp") is now available as an accumulator image together with | |
518 * a counter image. "Flatten" these two images into one. | |
519 **/ | |
520 | |
521 intermediate.reset(new Orthanc::Image | |
522 (TargetFormat, intermediateWidth, intermediateHeight, false)); | |
523 | |
524 maxValue = 0; | |
525 | |
526 for (unsigned int y = 0; y < intermediateHeight; y++) | |
527 { | |
528 const float *qacc = reinterpret_cast<const float*>(accumulator->GetConstRow(y)); | |
529 const uint16_t *qcount = reinterpret_cast<const uint16_t*>(counter->GetConstRow(y)); | |
530 typename TargetTraits::PixelType *p = | |
531 reinterpret_cast<typename TargetTraits::PixelType*>(intermediate->GetRow(y)); | |
532 | |
533 for (unsigned int x = 0; x < intermediateWidth; x++) | |
534 { | |
535 if (*qcount == 0) | |
536 { | |
537 TargetTraits::SetZero(*p); | |
538 } | |
539 else | |
540 { | |
541 *p = *qacc / static_cast<float>(*qcount); | |
542 | |
543 if (*p > maxValue) | |
544 { | |
545 maxValue = *p; | |
546 } | |
547 } | |
548 | |
549 p++; | |
550 qacc++; | |
551 qcount++; | |
552 } | |
553 } | |
554 | |
555 // We don't need the accumulator images anymore | |
556 accumulator.reset(NULL); | |
557 counter.reset(NULL); | |
558 | |
559 | |
560 /** | |
561 * Step 6: Apply the "warp" part of the transform to map the | |
562 * intermediate image to the final image. | |
563 **/ | |
564 | |
565 Matrix warp; | |
566 | |
567 { | |
568 // (5.a) Compute the "warp" matrix by removing the 3rd row and | |
569 // 3rd column from the GetWarp() matrix | |
570 // Check out: ../../Resources/Computations/ComputeWarp.py | |
571 | |
572 Matrix fullWarp = LinearAlgebra::Product | |
573 (shearWarp.GetIntrinsicParameters(), shearWarp.GetWarp()); | |
574 | |
575 const double v[] = { | |
576 fullWarp(0,0), fullWarp(0,1), fullWarp(0,3), | |
577 fullWarp(1,0), fullWarp(1,1), fullWarp(1,3), | |
578 fullWarp(2,0), fullWarp(2,1), fullWarp(2,3) | |
579 }; | |
580 | |
581 LinearAlgebra::FillMatrix(warp, 3, 3, v); | |
582 } | |
583 | |
584 // (5.b) Apply the projective transform to the image | |
585 ApplyProjectiveTransform(target, *intermediate, warp, warpInterpolation); | |
586 } | |
587 | |
588 | |
589 template <Orthanc::PixelFormat SourceFormat, | |
590 Orthanc::PixelFormat TargetFormat> | |
591 static void ApplyAxialInternal2(Orthanc::ImageAccessor& target, | |
592 float& maxValue, | |
593 const Matrix& M_view, | |
594 const ImageBuffer3D& source, | |
595 bool mip, | |
596 double pixelSpacing, | |
597 unsigned int countSlices, | |
598 ImageInterpolation shearInterpolation, | |
599 ImageInterpolation warpInterpolation) | |
600 { | |
601 if (mip) | |
602 { | |
603 ApplyAxialInternal<SourceFormat, TargetFormat, true> | |
604 (target, maxValue, M_view, source, pixelSpacing, | |
605 countSlices, shearInterpolation, warpInterpolation); | |
606 } | |
607 else | |
608 { | |
609 ApplyAxialInternal<SourceFormat, TargetFormat, false> | |
610 (target, maxValue, M_view, source, pixelSpacing, | |
611 countSlices, shearInterpolation, warpInterpolation); | |
612 } | |
613 } | |
614 | |
615 | |
616 Orthanc::ImageAccessor* | |
617 ShearWarpProjectiveTransform::ApplyAxial(float& maxValue, | |
618 const Matrix& M_view, | |
619 const ImageBuffer3D& source, | |
620 Orthanc::PixelFormat targetFormat, | |
621 unsigned int targetWidth, | |
622 unsigned int targetHeight, | |
623 bool mip, | |
624 double pixelSpacing, | |
625 unsigned int countSlices, | |
626 ImageInterpolation shearInterpolation, | |
627 ImageInterpolation warpInterpolation) | |
628 { | |
629 std::auto_ptr<Orthanc::ImageAccessor> target | |
630 (new Orthanc::Image(targetFormat, targetWidth, targetHeight, false)); | |
631 | |
632 if (source.GetFormat() == Orthanc::PixelFormat_Grayscale16 && | |
633 targetFormat == Orthanc::PixelFormat_Grayscale16) | |
634 { | |
635 ApplyAxialInternal2<Orthanc::PixelFormat_Grayscale16, | |
636 Orthanc::PixelFormat_Grayscale16> | |
637 (*target, maxValue, M_view, source, mip, pixelSpacing, | |
638 countSlices, shearInterpolation, warpInterpolation); | |
639 } | |
640 else | |
641 { | |
642 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); | |
643 } | |
644 | |
645 return target.release(); | |
646 } | |
329 } | 647 } |