Commit 13334211 authored by Andi Braimllari's avatar Andi Braimllari Committed by David Frank
Browse files

Add 2D circular limited angle trajectory generator

The new trajectory generator creates a circular trajectory,
where a certain portion of the trajectory is missing.

This commit also refactors out common code from both generators
and fixes a minor bug, where the trajectory would overshoot slightly
(add a pose at 360.2 or similar).
parent 5c098e26
Pipeline #703454 passed with stages
in 23 minutes and 3 seconds
......@@ -28,6 +28,7 @@
// Generators headers
#include "PhantomGenerator.h"
#include "CircleTrajectoryGenerator.h"
#include "LimitedAngleTrajectoryGenerator.h"
#include "SphereTrajectoryGenerator.h"
// IO headers
......
......@@ -3,16 +3,20 @@ set(MODULE_HEADERS
EllipseGenerator.h
PhantomGenerator.h
CircleTrajectoryGenerator.h
LimitedAngleTrajectoryGenerator.h
NoiseGenerators.h
SphereTrajectoryGenerator.h)
SphereTrajectoryGenerator.h
TrajectoryGenerator.h)
# list all the code files of the module
set(MODULE_SOURCES
EllipseGenerator.cpp
PhantomGenerator.cpp
CircleTrajectoryGenerator.cpp
LimitedAngleTrajectoryGenerator.cpp
NoiseGenerators.cpp
SphereTrajectoryGenerator.cpp)
SphereTrajectoryGenerator.cpp
TrajectoryGenerator.cpp)
list(APPEND MODULE_PUBLIC_DEPS elsa_core elsa_logging)
list(APPEND MODULE_PRIVATE_DEPS)
......
......@@ -25,30 +25,8 @@ namespace elsa
arcDegrees);
// Calculate size and spacing for each geometry pose using a IIFE
const auto [coeffs, spacing] = [&] {
IndexVector_t coeffs(dim);
RealVector_t spacing(dim);
// Scale coeffsPerDim by sqrt(2), this reduces undersampling of the corners, as the
// detector is larger than the volume. Cast back and forthe to reduce warnings...
// This has to be a RealVector_t, most likely that the cast happens, anyway we get
// errors down the line see #86 in Gitlab
const RealVector_t coeffsPerDim =
volumeDescriptor.getNumberOfCoefficientsPerDimension().template cast<real_t>();
const real_t sqrt2 = std::sqrt(2.f);
const auto coeffsPerDimScaled = (coeffsPerDim * sqrt2).template cast<index_t>();
const auto spacingPerDim = volumeDescriptor.getSpacingPerDimension();
coeffs.head(dim - 1) = coeffsPerDimScaled.head(dim - 1);
coeffs[dim - 1] = numberOfPoses; // TODO: with eigen 3.4: `coeffs(Eigen::last) = 1`
spacing.head(dim - 1) = spacingPerDim.head(dim - 1);
spacing[dim - 1] = 1; // TODO: same as coeffs
// Return a pair, then split it using structured bindings
return std::pair{coeffs, spacing};
}();
const auto [coeffs, spacing] =
calculateSizeAndSpacingPerGeometry(volumeDescriptor, numberOfPoses);
// Create vector and reserve the necessary size, minor optimization such that no new
// allocations are necessary in the loop
......
......@@ -3,6 +3,7 @@
#include "Geometry.h"
#include "DetectorDescriptor.h"
#include "TrajectoryGenerator.h"
#include <vector>
#include <utility>
......@@ -16,7 +17,7 @@ namespace elsa
* @author Maximilan Hornung - initial code
* @author Tobias Lasser - modernization, fixes
*/
class CircleTrajectoryGenerator
class CircleTrajectoryGenerator : public TrajectoryGenerator
{
public:
/**
......@@ -26,8 +27,10 @@ namespace elsa
* @param numberOfPoses the number of (equally spaced) acquisition poses to be generated
* @param volumeDescriptor the volume around which the trajectory should go
* @param arcDegrees the size of the arc of the circle covered by the trajectory (in
* degrees, 360 for full circle) @param sourceToCenter the distance of the X-ray source to
* the center of the volume @param centerToDetector the distance of the center of the volume
* degrees, 360 for full circle)
* @param sourceToCenter the distance of the X-ray source to
* the center of the volume
* @param centerToDetector the distance of the center of the volume
* to the X-ray detector
*
* @returns a pair containing the list of geometries with a circular trajectory, and the
......
#include "LimitedAngleTrajectoryGenerator.h"
#include "VolumeDescriptor.h"
#include "PlanarDetectorDescriptor.h"
#include "Logger.h"
#include <stdexcept>
namespace elsa
{
std::unique_ptr<DetectorDescriptor> LimitedAngleTrajectoryGenerator::createTrajectory(
index_t numberOfPoses,
std::pair<elsa::geometry::Degree, elsa::geometry::Degree> missingWedgeAngles,
const DataDescriptor& volumeDescriptor, index_t arcDegrees, real_t sourceToCenter,
real_t centerToDetector, bool mirrored)
{
// pull in geometry namespace, to reduce cluttering
using namespace geometry;
// sanity check
const auto dim = volumeDescriptor.getNumberOfDimensions();
if (dim != 2) {
throw InvalidArgumentError("LimitedAngleTrajectoryGenerator: can only handle 2D");
}
Logger::get("LimitedAngleTrajectoryGenerator")
->info("creating 2D trajectory with {} poses in an {} degree arc", numberOfPoses,
arcDegrees);
const auto [coeffs, spacing] =
calculateSizeAndSpacingPerGeometry(volumeDescriptor, numberOfPoses);
// create vector and reserve the necessary size, minor optimization such that no new
// allocations are necessary in the loop
std::vector<Geometry> geometryList;
geometryList.reserve(static_cast<std::size_t>(numberOfPoses));
real_t wedgeArc = mirrored ? 2 * (missingWedgeAngles.second - missingWedgeAngles.first)
: missingWedgeAngles.second - missingWedgeAngles.first;
const real_t angleIncrement =
(static_cast<real_t>(arcDegrees) - wedgeArc) / (static_cast<real_t>(numberOfPoses));
for (index_t i = 0;; ++i) {
Radian angle = Degree{static_cast<real_t>(i) * angleIncrement};
if (angle.to_degree() >= static_cast<real_t>(arcDegrees)) {
break;
}
if (notInMissingWedge(angle, missingWedgeAngles, mirrored)) {
// use emplace_back, then no copy is created
geometryList.emplace_back(SourceToCenterOfRotation{sourceToCenter},
CenterOfRotationToDetector{centerToDetector},
Radian{angle},
VolumeData2D{volumeDescriptor.getSpacingPerDimension(),
volumeDescriptor.getLocationOfOrigin()},
SinogramData2D{Size2D{coeffs}, Spacing2D{spacing}});
}
}
return std::make_unique<PlanarDetectorDescriptor>(coeffs, spacing, geometryList);
}
bool LimitedAngleTrajectoryGenerator::notInMissingWedge(
elsa::geometry::Radian angle,
std::pair<elsa::geometry::Degree, elsa::geometry::Degree> missingWedgeAngles, bool mirrored)
{
if (!mirrored) {
return !(angle.to_degree() >= missingWedgeAngles.first
&& angle.to_degree() <= missingWedgeAngles.second);
} else {
return !((angle.to_degree() >= missingWedgeAngles.first
&& angle.to_degree() <= missingWedgeAngles.second)
|| (angle.to_degree() >= (missingWedgeAngles.first + 180)
&& angle.to_degree() <= (missingWedgeAngles.second + 180)));
}
}
} // namespace elsa
#pragma once
#include "Geometry.h"
#include "DetectorDescriptor.h"
#include "TrajectoryGenerator.h"
#include <vector>
#include <utility>
namespace elsa
{
/**
* @brief Generator for limited angle trajectories as used in X-ray Computed Tomography
* (for 2D/3D).
*
* @author Andi Braimllari - initial code
*/
class LimitedAngleTrajectoryGenerator : public TrajectoryGenerator
{
public:
/**
* @brief Generate a list of geometries corresponding to a limited angle trajectory around a
* volume.
*
* @param numberOfPoses the number of (equally spaced) acquisition poses to be generated
* @param missingWedgeAngles the two angles between which the missing wedge is located
* @param volumeDescriptor the volume around which the trajectory should go
* @param arcDegrees the size of the arc of the circle covered by the trajectory (in
* degrees, 360 for full circle)
* @param sourceToCenter the distance of the X-ray source to
* the center of the volume
* @param centerToDetector the distance of the center of the volume
* to the X-ray detector
* @param mirrored flag indicating whether the missing wedge will be mirrored or not,
* defaults to true
*
* @returns a pair containing the list of geometries with a circular trajectory, and the
* sinogram data descriptor
*
* Please note: the first pose will be at 0 degrees, the last pose will be at arcDegrees
* For example: 3 poses over a 180 arc will yield: 0, 90, 180 degrees.
*
* Please note: the sinogram size/spacing will match the volume size/spacing.
*
* TODO: Make it possible to return either PlanarDetectorDescriptor, or
* CurvedDetectorDescriptor
*/
static std::unique_ptr<DetectorDescriptor> createTrajectory(
index_t numberOfPoses,
std::pair<elsa::geometry::Degree, elsa::geometry::Degree> missingWedgeAngles,
const DataDescriptor& volumeDescriptor, index_t arcDegrees, real_t sourceToCenter,
real_t centerToDetector, bool mirrored = true);
private:
static bool notInMissingWedge(
elsa::geometry::Radian angle,
std::pair<elsa::geometry::Degree, elsa::geometry::Degree> missingWedgeAngles,
bool mirrored = true);
};
} // namespace elsa
......@@ -23,31 +23,8 @@ namespace elsa
"trajectories",
dim, numberOfPoses, numberOfCircles);
// Calculate size and spacing for each geometry pose using a IIFE
const auto [coeffs, spacing] = [&] {
IndexVector_t coeffs(dim);
RealVector_t spacing(dim);
// Scale coeffsPerDim by sqrt(2), this reduces undersampling of the corners, as the
// detector is larger than the volume. Cast back and forth to reduce warnings...
// This has to be a RealVector_t, most likely that the cast happens, anyway we get
// errors down the line see #86 in Gitlab
const RealVector_t coeffsPerDim =
volumeDescriptor.getNumberOfCoefficientsPerDimension().template cast<real_t>();
const real_t sqrt2 = std::sqrt(2.f);
const auto coeffsPerDimScaled = (coeffsPerDim * sqrt2).template cast<index_t>();
const auto spacingPerDim = volumeDescriptor.getSpacingPerDimension();
coeffs.head(dim - 1) = coeffsPerDimScaled.head(dim - 1);
coeffs[dim - 1] = numberOfPoses; // TODO: with eigen 3.4: `coeffs(Eigen::last) = 1`
spacing.head(dim - 1) = spacingPerDim.head(dim - 1);
spacing[dim - 1] = 1; // TODO: same as coeffs
// Return a pair, then split it using structured bindings
return std::pair{coeffs, spacing};
}();
const auto [coeffs, spacing] =
calculateSizeAndSpacingPerGeometry(volumeDescriptor, numberOfPoses);
// Create vector and reserve the necessary size, minor optimization such that no new
// allocations are necessary in the loop
......
......@@ -2,6 +2,7 @@
#include "Geometry.h"
#include "DetectorDescriptor.h"
#include "TrajectoryGenerator.h"
#include <vector>
#include <utility>
......@@ -15,7 +16,7 @@ namespace elsa
* @author Michael Loipführer - initial code
*/
class SphereTrajectoryGenerator
class SphereTrajectoryGenerator : public TrajectoryGenerator
{
public:
/**
......
#include "TrajectoryGenerator.h"
#include "VolumeDescriptor.h"
namespace elsa
{
std::pair<IndexVector_t, RealVector_t>
TrajectoryGenerator::calculateSizeAndSpacingPerGeometry(const DataDescriptor& volDescr,
index_t numberOfPoses)
{
const auto dim = volDescr.getNumberOfDimensions();
IndexVector_t coeffs(dim);
RealVector_t spacing(dim);
// Scale coeffsPerDim by sqrt(2), this reduces undersampling of the corners, as the
// detector is larger than the volume. Cast back and forthe to reduce warnings...
// This has to be a RealVector_t, most likely that the cast happens, anyway we get
// errors down the line see #86 in Gitlab
const RealVector_t coeffsPerDim =
volDescr.getNumberOfCoefficientsPerDimension().template cast<real_t>();
const real_t sqrt2 = std::sqrt(2.f);
const auto coeffsPerDimScaled = (coeffsPerDim * sqrt2).template cast<index_t>();
const auto spacingPerDim = volDescr.getSpacingPerDimension();
coeffs.head(dim - 1) = coeffsPerDimScaled.head(dim - 1);
coeffs[dim - 1] = numberOfPoses; // TODO: with eigen 3.4: `coeffs(Eigen::last) = 1`
spacing.head(dim - 1) = spacingPerDim.head(dim - 1);
spacing[dim - 1] = 1; // TODO: same as coeffs
// return a pair, then split it using structured bindings
return std::pair{coeffs, spacing};
}
} // namespace elsa
#pragma once
#include "VolumeDescriptor.h"
#include "PlanarDetectorDescriptor.h"
namespace elsa
{
/**
* @brief Parent class for other trajectory generator classes. Currently contains extracted
* duplicate code.
*
* @author Andi Braimllari - initial code
*/
class TrajectoryGenerator
{
protected:
static std::pair<IndexVector_t, RealVector_t>
calculateSizeAndSpacingPerGeometry(const DataDescriptor& volDescr,
index_t numberOfPoses);
};
} // namespace elsa
......@@ -17,4 +17,5 @@ add_custom_target(build-tests-generators)
ELSA_DOCTEST(EllipseGenerator)
ELSA_DOCTEST(PhantomGenerator)
ELSA_DOCTEST(CircleTrajectoryGenerator)
ELSA_DOCTEST(LimitedAngleTrajectoryGenerator)
ELSA_DOCTEST(SphereTrajectoryGenerator)
......@@ -81,12 +81,12 @@ TEST_CASE("CircleTrajectoryGenerator: Create a Circular Trajectory")
WHEN("We create a full circular trajectory for this scenario")
{
index_t halfCircular = 359;
index_t fullyCircular = 359;
real_t diffCenterSource{s * 100};
real_t diffCenterDetector{s};
auto sdesc = CircleTrajectoryGenerator::createTrajectory(
numberOfAngles, desc, halfCircular, diffCenterSource, diffCenterDetector);
numberOfAngles, desc, fullyCircular, diffCenterSource, diffCenterDetector);
// Check that the detector size is correct
REQUIRE_EQ(sdesc->getNumberOfCoefficientsPerDimension()[0], expectedDetectorSize);
......@@ -97,7 +97,7 @@ TEST_CASE("CircleTrajectoryGenerator: Create a Circular Trajectory")
const real_t sourceToCenter = diffCenterSource;
const real_t centerToDetector = diffCenterDetector;
real_t angle = static_cast<real_t>(1.0) * static_cast<real_t>(halfCircular)
real_t angle = static_cast<real_t>(1.0) * static_cast<real_t>(fullyCircular)
/ static_cast<real_t>(numberOfAngles - 1);
for (index_t i = 0; i < numberOfAngles; ++i) {
real_t currAngle = static_cast<real_t>(i) * angle * pi_t / 180.0f;
......@@ -182,12 +182,12 @@ TEST_CASE("CircleTrajectoryGenerator: Create a Circular Trajectory")
}
WHEN("We create a full circular trajectory for this scenario")
{
const index_t halfCircular = 359;
const index_t fullyCircular = 359;
real_t diffCenterSource{s * 100};
real_t diffCenterDetector{s};
auto sdesc = CircleTrajectoryGenerator::createTrajectory(
numberOfAngles, desc, halfCircular, diffCenterSource, diffCenterDetector);
numberOfAngles, desc, fullyCircular, diffCenterSource, diffCenterDetector);
// Check that the detector size is correct
REQUIRE_EQ(sdesc->getNumberOfCoefficientsPerDimension()[0], expectedDetectorSize);
......@@ -199,7 +199,7 @@ TEST_CASE("CircleTrajectoryGenerator: Create a Circular Trajectory")
const auto sourceToCenter = diffCenterSource;
const auto centerToDetector = diffCenterDetector;
real_t angleInc = 1.0f * static_cast<real_t>(halfCircular)
real_t angleInc = 1.0f * static_cast<real_t>(fullyCircular)
/ static_cast<real_t>(numberOfAngles - 1);
for (index_t i = 0; i < numberOfAngles; ++i) {
real_t angle = static_cast<real_t>(i) * angleInc * pi_t / 180.0f;
......
/**
* @file test_LimitedAngleTrajectoryGenerator.cpp
*
* @brief Tests for the LimitedAngleTrajectoryGenerator class
*
* @author Andi Braimllari
*/
#include "LimitedAngleTrajectoryGenerator.h"
#include "Logger.h"
#include "VolumeDescriptor.h"
#include "testHelpers.h"
#include "doctest/doctest.h"
using namespace elsa;
using namespace doctest;
TEST_CASE("LimitedAngleTrajectoryGenerator: Create a mirrored Limited Angle Trajectory")
{
using namespace geometry;
const index_t s = 64;
// Detector size is the volume size scaled by the square root of 2
const auto expectedDetectorSize = static_cast<index_t>(s * std::sqrt(2));
GIVEN("A 2D descriptor and 256 angles")
{
index_t numberOfAngles = 256;
IndexVector_t volSize(2);
volSize << s, s;
VolumeDescriptor desc{volSize};
WHEN("We create a half circular limited angle trajectory for this scenario")
{
index_t halfCircular = 180;
real_t diffCenterSource{s * 100};
real_t diffCenterDetector{s};
std::pair missingWedgeAngles(geometry::Degree(110), geometry::Degree(170));
auto sinoDescriptor = LimitedAngleTrajectoryGenerator::createTrajectory(
numberOfAngles, missingWedgeAngles, desc, halfCircular, diffCenterSource,
diffCenterDetector);
// check that the detector size is correct
REQUIRE_EQ(sinoDescriptor->getNumberOfCoefficientsPerDimension()[0],
expectedDetectorSize);
THEN("Every geomList in our list has the same camera center and the same projection "
"matrix")
{
const real_t sourceToCenter = diffCenterSource;
const real_t centerToDetector = diffCenterDetector;
real_t wedgeArc = 2 * (missingWedgeAngles.second - missingWedgeAngles.first);
const real_t angleIncrement = (static_cast<real_t>(halfCircular) - wedgeArc)
/ (static_cast<real_t>(numberOfAngles));
index_t j = 0;
for (index_t i = 0;; ++i) {
Radian angle = Degree{static_cast<real_t>(i) * angleIncrement};
if (angle.to_degree() >= static_cast<real_t>(halfCircular)) {
break;
}
if (!((angle.to_degree() >= missingWedgeAngles.first
&& angle.to_degree() <= missingWedgeAngles.second)
|| (angle.to_degree() >= (missingWedgeAngles.first + 180)
&& angle.to_degree() <= (missingWedgeAngles.second + 180)))) {
Geometry tmpGeom(SourceToCenterOfRotation{sourceToCenter},
CenterOfRotationToDetector{centerToDetector},
Radian{angle}, VolumeData2D{volSize},
SinogramData2D{sinoDescriptor->getSpacingPerDimension(),
sinoDescriptor->getLocationOfOrigin()});
auto geom = sinoDescriptor->getGeometryAt(j++);
CHECK(geom);
const auto centerNorm =
(tmpGeom.getCameraCenter() - geom->getCameraCenter()).norm();
const auto projMatNorm =
(tmpGeom.getProjectionMatrix() - geom->getProjectionMatrix()).norm();
const auto invProjMatNorm = (tmpGeom.getInverseProjectionMatrix()
- geom->getInverseProjectionMatrix())
.norm();
REQUIRE_UNARY(checkApproxEq(centerNorm, 0));
REQUIRE_UNARY(checkApproxEq(projMatNorm, 0, 0.0000001));
REQUIRE_UNARY(checkApproxEq(invProjMatNorm, 0, 0.0000001));
}
}
}
}
WHEN("We create a full circular limited angle trajectory for this scenario")
{
index_t fullyCircular = 359;
real_t diffCenterSource{s * 100};
real_t diffCenterDetector{s};
std::pair missingWedgeAngles(geometry::Degree(30), geometry::Degree(70));
auto sinoDescriptor = LimitedAngleTrajectoryGenerator::createTrajectory(
numberOfAngles, missingWedgeAngles, desc, fullyCircular, diffCenterSource,
diffCenterDetector);
// check that the detector size is correct
REQUIRE_EQ(sinoDescriptor->getNumberOfCoefficientsPerDimension()[0],
expectedDetectorSize);
THEN("Every geomList in our list has the same camera center and the same projection "
"matrix")
{
const real_t sourceToCenter = diffCenterSource;
const real_t centerToDetector = diffCenterDetector;
real_t wedgeArc = 2 * (missingWedgeAngles.second - missingWedgeAngles.first);
const real_t angleIncrement = (static_cast<real_t>(fullyCircular) - wedgeArc)
/ (static_cast<real_t>(numberOfAngles));
index_t j = 0;
for (index_t i = 0;; ++i) {
Radian angle = Degree{static_cast<real_t>(i) * angleIncrement};
if (angle.to_degree() > static_cast<real_t>(fullyCircular)) {
break;
}
if (!((angle.to_degree() >= missingWedgeAngles.first
&& angle.to_degree() <= missingWedgeAngles.second)
|| (angle.to_degree() >= (missingWedgeAngles.first + 180)
&& angle.to_degree() <= (missingWedgeAngles.second + 180)))) {
Geometry tmpGeom(SourceToCenterOfRotation{sourceToCenter},
CenterOfRotationToDetector{centerToDetector},
Radian{angle}, VolumeData2D{volSize},
SinogramData2D{sinoDescriptor->getSpacingPerDimension(),
sinoDescriptor->getLocationOfOrigin()});
auto geom = sinoDescriptor->getGeometryAt(j++);
CHECK(geom);
const auto centerNorm =
(tmpGeom.getCameraCenter() - geom->getCameraCenter()).norm();
const auto projMatNorm =
(tmpGeom.getProjectionMatrix() - geom->getProjectionMatrix()).norm();
const auto invProjMatNorm = (tmpGeom.getInverseProjectionMatrix()
- geom->getInverseProjectionMatrix())
.norm();
REQUIRE_UNARY(checkApproxEq(centerNorm, 0));
REQUIRE_UNARY(checkApproxEq(projMatNorm, 0, 0.0000001));
REQUIRE_UNARY(checkApproxEq(invProjMatNorm, 0, 0.0000001));
}
}
}
}
}
}
TEST_CASE("LimitedAngleTrajectoryGenerator: Create a non-mirrored Limited Angle Trajectory")
{
using namespace geometry;
const index_t s = 64;
// Detector size is the volume size scaled by the square root of 2
const auto expectedDetectorSize = static_cast<index_t>(s * std::sqrt(2));
GIVEN("A 2D descriptor and 256 angles")
{
index_t numberOfAngles = 256;
IndexVector_t volSize(2);
volSize << s, s;
VolumeDescriptor desc{volSize};
WHEN("We create a half circular limited angle trajectory for this scenario")
{
index_t halfCircular = 180;
real_t diffCenterSource{s * 100};
real_t diffCenterDetector{s};
std::pair missingWedgeAngles(geometry::Degree(40), geometry::Degree(45));
auto sinoDescriptor = LimitedAngleTrajectoryGenerator::createTrajectory(