Line data    Source code

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