Line data    Source code

       1             : #include "SphereTrajectoryGenerator.h"
       2             : #include "PlanarDetectorDescriptor.h"
       3             : #include "Logger.h"
       4             : 
       5             : namespace elsa
       6             : {
       7             :     std::unique_ptr<DetectorDescriptor> SphereTrajectoryGenerator::createTrajectory(
       8             :         index_t numberOfPoses, const DataDescriptor& volumeDescriptor, index_t numberOfCircles,
       9             :         geometry::SourceToCenterOfRotation sourceToCenter,
      10             :         geometry::CenterOfRotationToDetector centerToDetector)
      11          11 :     {
      12             :         // pull in geometry namespace, to reduce cluttering
      13          11 :         using namespace geometry;
      14             : 
      15             :         // sanity check
      16          11 :         const auto dim = volumeDescriptor.getNumberOfDimensions();
      17             : 
      18          11 :         if (dim != 3)
      19           1 :             throw InvalidArgumentError("SphereTrajectoryGenerator: can only handle 3d");
      20             : 
      21          10 :         Logger::get("SphereTrajectoryGenerator")
      22          10 :             ->info("creating {}D trajectory with {} poses consisting of {} circular "
      23          10 :                    "trajectories",
      24          10 :                    dim, numberOfPoses, numberOfCircles);
      25             : 
      26          10 :         const auto [coeffs, spacing] =
      27          10 :             calculateSizeAndSpacingPerGeometry(volumeDescriptor, numberOfPoses);
      28             : 
      29             :         // Create vector and reserve the necessary size, minor optimization such that no new
      30             :         // allocations are necessary in the loop
      31          10 :         std::vector<Geometry> geometryList;
      32          10 :         geometryList.reserve(static_cast<std::size_t>(numberOfPoses));
      33             : 
      34             :         // compute the total circumference of all circular trajectories we are generating
      35          10 :         double circumference = 0;
      36          57 :         for (index_t i = 0; i < numberOfCircles; ++i) {
      37          47 :             circumference +=
      38          47 :                 M_PI * 2.0 * centerToDetector
      39          47 :                 * cos(M_PI * static_cast<double>(i + 1) / static_cast<double>(numberOfCircles + 2)
      40          47 :                       - M_PI_2);
      41          47 :         }
      42             : 
      43          10 :         auto circumferenceChangePerPose = circumference / static_cast<double>(numberOfPoses);
      44          10 :         index_t createdPoses = 0;
      45             : 
      46          57 :         for (index_t i = 0; i < numberOfCircles; ++i) {
      47          47 :             auto beta = M_PI * static_cast<double>(i + 1) / static_cast<double>(numberOfCircles + 2)
      48          47 :                         - M_PI_2;
      49             : 
      50          47 :             index_t posesInCircle;
      51          47 :             if (i == numberOfCircles - 1) {
      52          10 :                 posesInCircle = numberOfPoses - createdPoses;
      53          37 :             } else {
      54          37 :                 posesInCircle = static_cast<index_t>(
      55          37 :                     (2.0 * M_PI * static_cast<double>(centerToDetector) * cos(beta))
      56          37 :                     / circumferenceChangePerPose);
      57          37 :             }
      58             : 
      59         495 :             for (index_t poseIndex = 0; poseIndex < posesInCircle; ++poseIndex) {
      60         448 :                 auto gamma = M_PI * 2.0 * static_cast<double>(poseIndex)
      61         448 :                              / static_cast<double>(posesInCircle);
      62             : 
      63         448 :                 geometryList.emplace_back(sourceToCenter, centerToDetector,
      64         448 :                                           VolumeData3D{volumeDescriptor.getSpacingPerDimension(),
      65         448 :                                                        volumeDescriptor.getLocationOfOrigin()},
      66         448 :                                           SinogramData3D{Size3D{coeffs}, Spacing3D{spacing}},
      67         448 :                                           RotationAngles3D{Radian{static_cast<real_t>(gamma)},
      68         448 :                                                            Radian{static_cast<real_t>(beta)},
      69         448 :                                                            Radian{0}});
      70         448 :             }
      71             : 
      72          47 :             createdPoses += posesInCircle;
      73          47 :         }
      74             : 
      75          10 :         return std::make_unique<PlanarDetectorDescriptor>(coeffs, spacing, geometryList);
      76          10 :     }
      77             : 
      78             : } // namespace elsa