11.08., 9:00 - 11:00: Due to updates GitLab will be unavailable for some minutes between 09:00 and 11:00.

Commit 7397e31a authored by David Frank's avatar David Frank Committed by Tobias Lasser

#63 Add DetectorDescriptor and PlanarDetectorDescriptor

- Add Abstract class DetectorDescriptor, which is derived from DataDescriptor (Support ray generation from a given pose and detector pixel and given pose and voxel)
- Add first derived class of DetectorDescriptor: PlanarDetectorDescriptor
- Add usage of DetectorDescriptor in Projectors (remove dependency to geometry, ray generation is only handled by DetectorDescriptor, adapt tests)
- Restructure tests of CUDA projectors to make failing tests more readable
parent f6e6f3d7
Pipeline #287374 failed with stages
in 3 minutes and 46 seconds
......@@ -13,6 +13,8 @@ set(MODULE_HEADERS
Descriptors/DataDescriptor.h
Descriptors/DescriptorUtils.h
Descriptors/VolumeDescriptor.h
Descriptors/DetectorDescriptor.h
Descriptors/PlanarDetectorDescriptor.h
Descriptors/BlockDescriptor.h
Descriptors/IdenticalBlocksDescriptor.h
Descriptors/PartitionDescriptor.h
......@@ -32,9 +34,11 @@ set(MODULE_HEADERS
set(MODULE_SOURCES
Descriptors/DataDescriptor.cpp
Descriptors/VolumeDescriptor.cpp
Descriptors/PlanarDetectorDescriptor.cpp
Descriptors/RandomBlocksDescriptor.cpp
Descriptors/DescriptorUtils.cpp
Descriptors/IdenticalBlocksDescriptor.cpp
Descriptors/DetectorDescriptor.cpp
Descriptors/PartitionDescriptor.cpp
DataContainer.cpp
Handlers/DataHandlerCPU.cpp
......
#include "DetectorDescriptor.h"
#include <iostream>
namespace elsa
{
DetectorDescriptor::DetectorDescriptor(const IndexVector_t& numOfCoeffsPerDim,
const std::vector<Geometry>& geometryList)
: DataDescriptor(numOfCoeffsPerDim), _geometry(geometryList)
{
// TODO Clarify: What about empty geometryList? Do we want to support it, or throw an
// execption?
}
DetectorDescriptor::DetectorDescriptor(const IndexVector_t& numOfCoeffsPerDim,
const RealVector_t& spacingPerDim,
const std::vector<Geometry>& geometryList)
: DataDescriptor(numOfCoeffsPerDim, spacingPerDim), _geometry(geometryList)
{
}
DetectorDescriptor::Ray
DetectorDescriptor::computeRayFromDetectorCoord(const index_t detectorIndex) const
{
// Return empty, if access out of bounds
assert(detectorIndex < getNumberOfCoefficients()
&& "PlanarDetectorDescriptor::computeRayToDetector(index_t): Assumption "
"detectorIndex smaller than number of coeffs, broken");
auto coord = getCoordinateFromIndex(detectorIndex);
return computeRayFromDetectorCoord(coord);
}
DetectorDescriptor::Ray
DetectorDescriptor::computeRayFromDetectorCoord(const IndexVector_t coord) const
{
// Assume all of the coordinates are inside of the volume
// auto tmp = (coord.array() < getNumberOfCoefficientsPerDimension().array());
// assert(tmp.all()
// && "DetectorDescriptor::computeRayToDetector(IndexVector_t): Assumption coord "
// "in bound wrong");
auto dim = getNumberOfDimensions();
// Cast to real_t and shift to center of pixel
auto detectorCoord = coord.head(dim - 1).template cast<real_t>().array() + 0.5;
// Last dimension is always the pose index
auto poseIndex = coord[dim - 1];
return computeRayFromDetectorCoord(detectorCoord, poseIndex);
}
index_t DetectorDescriptor::getNumberOfGeometryPoses() const { return _geometry.size(); }
std::optional<Geometry> DetectorDescriptor::getGeometryAt(const index_t index) const
{
if (_geometry.size() <= std::make_unsigned_t<std::size_t>(index))
return {};
return _geometry[index];
}
bool DetectorDescriptor::isEqual(const DataDescriptor& other) const
{
if (!DataDescriptor::isEqual(other))
return false;
// static cast as type checked in base comparison
auto otherBlock = static_cast<const DetectorDescriptor*>(&other);
if (getNumberOfGeometryPoses() != otherBlock->getNumberOfGeometryPoses())
return false;
return std::equal(std::cbegin(_geometry), std::cend(_geometry),
std::cbegin(otherBlock->_geometry));
}
} // namespace elsa
#pragma once
#include "elsaDefines.h"
#include "DataDescriptor.h"
#include "Geometry.h"
#include <optional>
#include "Eigen/Geometry"
namespace elsa
{
/**
* @brief Class representing metadata for lineraized n-dimensional signal stored in memory. It
* is a base class for different type signals caputred by some kind of detectors (i.e. a planar
* detector, curved or some other shaped detector). Is additionally stored information about the
* different poses of a trajectory.
*/
class DetectorDescriptor : public DataDescriptor
{
public:
using Ray = Eigen::ParametrizedLine<real_t, Eigen::Dynamic>;
public:
/// There is not default signal
DetectorDescriptor() = delete;
/// Default destructor
~DetectorDescriptor() override = default;
/**
* @brief Construct a DetectorDescriptor with a number of coefficients for each dimension
* and a list of geometry poses
*/
DetectorDescriptor(const IndexVector_t& numOfCoeffsPerDim,
const std::vector<Geometry>& geometryList);
/**
* @brief Construct a DetectorDescriptor with a number of coefficients and spacing for each
* dimension and a list of geometry poses
*/
DetectorDescriptor(const IndexVector_t& numOfCoeffsPerDim,
const RealVector_t& spacingPerDim,
const std::vector<Geometry>& geometryList);
/**
* @brief Overload of computeRayToDetector with a single detectorIndex. Compute the pose and
* coord index using getCoordinateFromIndex and call overload
*/
Ray computeRayFromDetectorCoord(const index_t detectorIndex) const;
/**
* @brief Overload of computeRayToDetector with a single coord vector. This vector encodes
* the pose index and the detector coordinate. So for a 1D detector, this will be 2D and the
* second dimension, will reference the pose index
*/
Ray computeRayFromDetectorCoord(const IndexVector_t coord) const;
/**
* @brief Compute a ray from the source from a pose to the given detector coordinate
*
* @param[in] detectorCoord Vector of size dim - 1, specifying the coordinate the ray should
* hit
* @param[in] poseIndex index into geometryList array, which pose to use for ray computation
*
*/
virtual Ray computeRayFromDetectorCoord(const RealVector_t& detectorCoord,
const index_t poseIndex) const = 0;
/**
* @brief Compute a ray from the source to trougth a pixel/voxel
*/
virtual RealVector_t computeDetectorCoordFromRay(const Ray& ray,
const index_t poseIndex) const = 0;
/// Get the number of poses used in the geometry
index_t getNumberOfGeometryPoses() const;
/// Get the i-th geometry in the trajectory.
std::optional<Geometry> getGeometryAt(const index_t index) const;
protected:
/// implement the polymorphic comparison operation
bool isEqual(const DataDescriptor& other) const override;
/// List of geometry poses
std::vector<Geometry> _geometry;
};
} // namespace elsa
#include "PlanarDetectorDescriptor.h"
#include <iostream>
namespace elsa
{
PlanarDetectorDescriptor::PlanarDetectorDescriptor(const IndexVector_t& numOfCoeffsPerDim,
const std::vector<Geometry>& geometryList)
: DetectorDescriptor(numOfCoeffsPerDim, geometryList)
{
}
PlanarDetectorDescriptor::PlanarDetectorDescriptor(const IndexVector_t& numOfCoeffsPerDim,
const RealVector_t& spacingPerDim,
const std::vector<Geometry>& geometryList)
: DetectorDescriptor(numOfCoeffsPerDim, spacingPerDim, geometryList)
{
}
DetectorDescriptor::Ray
PlanarDetectorDescriptor::computeRayFromDetectorCoord(const RealVector_t& detectorCoord,
const index_t poseIndex) const
{
// Assert that for all dimension of detectorCoord is in bounds and poseIndex can
// be index in the _geometry. If not the calculation will not be correct, but
// as this is the hot path, I don't want execptions and unpacking everything
// We'll just have to ensure, that we don't mess up in our hot path! :-)
assert((detectorCoord.block(0, 0, getNumberOfDimensions() - 1, 0).array()
< getNumberOfCoefficientsPerDimension()
.block(0, 0, getNumberOfDimensions() - 1, 0)
.template cast<real_t>()
.array())
.all()
&& "PlanarDetectorDescriptor::computeRayToDetector: Assumption detectorCoord in "
"bounds, wrong");
assert(std::make_unsigned_t<std::size_t>(poseIndex) < _geometry.size()
&& "PlanarDetectorDescriptor::computeRayToDetector: Assumption poseIndex smaller "
"than number of poses, wrong");
auto dim = getNumberOfDimensions();
// get the pose of trajectory
auto geometry = _geometry[std::make_unsigned_t<index_t>(poseIndex)];
auto projInvMatrix = geometry.getInverseProjectionMatrix();
// homogeneous coordinates [p;1], with p in detector space
RealVector_t homogeneousPixelCoord(dim);
homogeneousPixelCoord << detectorCoord, 1;
// Camera center is always the ray origin
auto ro = geometry.getCameraCenter();
auto rd = (projInvMatrix * homogeneousPixelCoord) // Matrix-Vector multiplication
.head(dim) // Transform to non-homogeneous
.normalized(); // normalize vector
return Ray(ro, rd);
}
RealVector_t
PlanarDetectorDescriptor::computeDetectorCoordFromRay(const Ray& ray,
const index_t poseIndex) const
{
auto dim = getNumberOfDimensions();
auto geometry = _geometry[(poseIndex)];
auto projMatrix = geometry.getProjectionMatrix();
// Only take the square matrix part
auto pixel = (projMatrix.block(0, 0, dim, dim) * ray.direction()).head(dim - 1);
return pixel;
}
bool PlanarDetectorDescriptor::isEqual(const DataDescriptor& other) const
{
// PlanarDetectorDescriptor has no data, so just deligate it to base class
return DetectorDescriptor::isEqual(other);
}
PlanarDetectorDescriptor* PlanarDetectorDescriptor::cloneImpl() const
{
return new PlanarDetectorDescriptor(getNumberOfCoefficientsPerDimension(),
getSpacingPerDimension(), _geometry);
}
} // namespace elsa
#pragma once
#include "DetectorDescriptor.h"
namespace elsa
{
/**
* @brief Class representing metadata for lineraized n-dimensional signal stored in memory. It
* specifically describes signals, which were captured by a planar detector and stores
* additional information such as different poses
*
* @author David Frank - initial code
*/
class PlanarDetectorDescriptor : public DetectorDescriptor
{
using DetectorDescriptor::Ray;
public:
PlanarDetectorDescriptor() = delete;
~PlanarDetectorDescriptor() = default;
/**
* @brief Construct a PlanatDetectorDescriptor with given number of coefficients and spacing
* per dimension and a list of geometry poses in the trajectory
*/
PlanarDetectorDescriptor(const IndexVector_t& numOfCoeffsPerDim,
const RealVector_t& spacingPerDim,
const std::vector<Geometry>& geometryList);
/**
* @brief Construct a PlanatDetectorDescriptor with given number of coefficients
* per dimension and a list of geometry poses in the trajectory
*/
PlanarDetectorDescriptor(const IndexVector_t& numOfCoeffsPerDim,
const std::vector<Geometry>& geometryList);
using DetectorDescriptor::computeRayFromDetectorCoord;
/// Override function to compute rays for a planar detector
Ray computeRayFromDetectorCoord(const RealVector_t& detectorCoord,
const index_t poseIndex) const override;
RealVector_t computeDetectorCoordFromRay(const Ray& ray,
const index_t poseIndex) const override;
private:
PlanarDetectorDescriptor* cloneImpl() const override;
bool isEqual(const DataDescriptor& other) const override;
};
} // namespace elsa
......@@ -147,20 +147,6 @@ namespace elsa
buildMatrices();
}
std::pair<RealVector_t, RealVector_t> Geometry::computeRayTo(const RealVector_t& p) const
{
// the ray origin is always the camera center
RealVector_t ro = _C;
// homogeneous coordinates [p;1] - p is in detector space
RealVector_t homP(_objectDimension);
homP << p, 1;
// multiplication of inverse projection matrix and homogeneous detector coordinate
auto rd = (_Pinv * homP).normalized();
return std::make_pair(ro, rd.head(_objectDimension));
}
const RealMatrix_t& Geometry::getProjectionMatrix() const { return _P; }
const RealMatrix_t& Geometry::getInverseProjectionMatrix() const { return _Pinv; }
......
......@@ -109,18 +109,6 @@ namespace elsa
real_t centerOfRotationOffsetY = static_cast<real_t>(0.0),
real_t centerOfRotationOffsetZ = static_cast<real_t>(0.0));
/**
* @brief Compute a ray (ray origin and ray direction) that hits the specified point p on
* the detector
*
* @param[in] p point p on detector
*
* @returns pair of <ro, rd>, where ro = ray origin and rd = ray direction
*
* Computation are done using the projection matrix.
*/
std::pair<RealVector_t, RealVector_t> computeRayTo(const RealVector_t& p) const;
/**
* @brief Return the projection matrix
*
......
......@@ -248,9 +248,9 @@ namespace elsa
}
/// Access to gamma
constexpr Radian gamma() const { return operator[](0); }
constexpr Radian gamma() const { return operator[](0u); }
/// Access to beta
constexpr Radian beta() const { return operator[](1); }
constexpr Radian beta() const { return operator[](1u); }
/// Access to alpha
constexpr Radian alpha() const { return operator[](2); }
};
......
......@@ -17,6 +17,7 @@ ELSA_TEST(LinearOperator)
ELSA_TEST(ExpressionTemplates)
ELSA_TEST(Geometry)
ELSA_TEST(StrongTypes)
ELSA_TEST(PlanarDetectorDescriptor)
if(ELSA_BENCHMARKS)
ELSA_TEST(BenchmarkExpressionTemplates)
......
This diff is collapsed.
/**
* @file test_PlanarDetectorDescriptor.cpp
*
* @brief Test for PlanarDetectorDescriptor
*
* @author David Frank - initial code
*/
#include <catch2/catch.hpp>
#include "PlanarDetectorDescriptor.h"
#include "VolumeDescriptor.h"
#include <iostream>
using namespace elsa;
using namespace elsa::geometry;
using Ray = DetectorDescriptor::Ray;
SCENARIO("Testing 2D PlanarDetectorDescriptor")
{
GIVEN("Given a 5x5 Volume and a single 5 wide detector pose")
{
IndexVector_t volSize(2);
volSize << 5, 5;
VolumeDescriptor ddVol(volSize);
IndexVector_t sinoSize(2);
sinoSize << 5, 1;
real_t s2c = 10;
real_t c2d = 4;
Geometry g(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Radian{0},
VolumeData2D{Size2D{volSize}}, SinogramData2D{Size2D{sinoSize}});
PlanarDetectorDescriptor desc(sinoSize, {g});
WHEN("Retreiving the single geometry pose")
{
auto geom = desc.getGeometryAt(0);
CHECK(desc.getNumberOfGeometryPoses() == 1);
THEN("Geometry is equal") { CHECK((geom) == g); }
}
WHEN("Generating rays for detecor pixels 0, 2 and 4")
{
for (real_t detPixel : std::initializer_list<real_t>{0, 2, 4}) {
RealVector_t pixel(1);
pixel << detPixel + 0.5f;
// Check that ray for IndexVector_t is equal to previous one
auto ray = desc.computeRayFromDetectorCoord(pixel, 0);
// Create variables, which make typing quicker
auto ro = ray.origin();
auto rd = ray.direction();
// Check that ray origin is camera center
auto c = g.getCameraCenter();
CHECK((ro - c).sum() == Approx(0));
// compute intersection manually
real_t t = Approx(rd[0]) == 0 ? (s2c + c2d) : ((pixel[0] - ro[0]) / rd[0]);
auto detCoordY = ro[1] + t * rd[1];
CHECK(detCoordY == Approx(ddVol.getLocationOfOrigin()[1] + c2d));
}
}
WHEN("Computing detector coord from ray")
{
auto ro = g.getCameraCenter();
auto rd = RealVector_t(2);
THEN("The detector coord is correct")
{
rd << -0.141421f, 0.989949f;
rd.normalize();
auto detCoord = desc.computeDetectorCoordFromRay(Ray(ro, rd), 0);
CHECK(detCoord[0] == Approx(0.5).margin(0.05));
}
THEN("The detector coord is correct asd")
{
rd << 0.f, 1.f;
rd.normalize();
auto detCoord = desc.computeDetectorCoordFromRay(Ray(ro, rd), 0);
CHECK(detCoord[0] == Approx(2.5));
}
THEN("The detector coord is correct asd")
{
rd << 0.141421f, 0.989949f;
rd.normalize();
auto detCoord = desc.computeDetectorCoordFromRay(Ray(ro, rd), 0);
CHECK(detCoord[0] == Approx(4.5).margin(0.05));
}
}
GIVEN("Given a 5x5 Volume and a multiple 5 wide detector pose")
{
IndexVector_t volSize(2);
volSize << 5, 5;
VolumeDescriptor ddVol(volSize);
IndexVector_t sinoSize(2);
sinoSize << 5, 4;
real_t s2c = 10;
real_t c2d = 4;
Geometry g1(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Degree{0},
VolumeData2D{Size2D{volSize}}, SinogramData2D{Size2D{sinoSize}});
Geometry g2(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Degree{90},
VolumeData2D{Size2D{volSize}}, SinogramData2D{Size2D{sinoSize}});
Geometry g3(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Degree{180},
VolumeData2D{Size2D{volSize}}, SinogramData2D{Size2D{sinoSize}});
Geometry g4(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Degree{270},
VolumeData2D{Size2D{volSize}}, SinogramData2D{Size2D{sinoSize}});
PlanarDetectorDescriptor desc(sinoSize, {g1, g2, g3, g4});
WHEN("Retreiving the single geometry pose")
{
auto geom = desc.getGeometryAt(0);
CHECK(desc.getNumberOfGeometryPoses() == 4);
THEN("Geometry is equal") { CHECK((desc.getGeometryAt(0)) == g1); }
THEN("Geometry is equal") { CHECK((desc.getGeometryAt(1)) == g2); }
THEN("Geometry is equal") { CHECK((desc.getGeometryAt(2)) == g3); }
THEN("Geometry is equal") { CHECK((desc.getGeometryAt(3)) == g4); }
}
WHEN("Check for multiple poses, that all the overloads compute the same rays")
{
for (index_t pose : {0, 1, 2, 3}) {
for (index_t detPixel : {0, 2, 4}) {
IndexVector_t pixel(2);
pixel << detPixel, pose;
RealVector_t pixelReal(1);
pixelReal << static_cast<real_t>(detPixel) + 0.5f;
auto ray1 =
desc.computeRayFromDetectorCoord(desc.getIndexFromCoordinate(pixel));
auto ray2 = desc.computeRayFromDetectorCoord(pixel);
auto ray3 = desc.computeRayFromDetectorCoord(pixelReal, pose);
auto ro1 = ray1.origin();
auto rd1 = ray1.direction();
auto ro2 = ray2.origin();
auto rd2 = ray2.direction();
auto ro3 = ray3.origin();
auto rd3 = ray3.direction();
CHECK(ro1 == ro2);
CHECK(ro1 == ro3);
CHECK(rd1 == rd2);
CHECK(rd1 == rd3);
// Shouldn't be necessary, but whatever
CHECK(ro2 == ro3);
CHECK(rd2 == rd3);
}
}
}
}
}
}
SCENARIO("Testing 3D PlanarDetectorDescriptor")
{
GIVEN("Given a 5x5x5 Volume and a single 5x5 wide detector pose")
{
IndexVector_t volSize(3);
volSize << 5, 5, 5;
VolumeDescriptor ddVol(volSize);
IndexVector_t sinoSize(3);
sinoSize << 5, 5, 1;
real_t s2c = 10;
real_t c2d = 4;
Geometry g(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d},
VolumeData3D{Size3D{volSize}}, SinogramData3D{Size3D{sinoSize}},
RotationAngles3D{Gamma{0}});
PlanarDetectorDescriptor desc(sinoSize, {g});
WHEN("Retreiving the single geometry pose")
{
auto geom = desc.getGeometryAt(0);
CHECK(desc.getNumberOfGeometryPoses() == 1);
THEN("Geometry is equal") { CHECK((geom) == g); }
}
WHEN("Generating rays for detecor pixels 0, 2 and 4 for each dim")
{
for (index_t detPixel1 : {0, 2, 4}) {
for (index_t detPixel2 : {0, 2, 4}) {
RealVector_t pixel(2);
pixel << static_cast<real_t>(detPixel1) + 0.5f,
static_cast<real_t>(detPixel2) + 0.5f;
// Check that ray for IndexVector_t is equal to previous one
auto ray = desc.computeRayFromDetectorCoord(pixel, 0);
// Create variables, which make typing quicker
auto ro = ray.origin();
auto rd = ray.direction();
// Check that ray origin is camera center
auto c = g.getCameraCenter();
CHECK((ro - c).sum() == Approx(0));
auto o = ddVol.getLocationOfOrigin();
RealVector_t detCoordWorld(3);
detCoordWorld << pixel[0] - o[0], pixel[1] - o[1], c2d;
RealVector_t rotD = g.getRotationMatrix().transpose() * detCoordWorld + o;
real_t factor = 0;
if (std::abs(rd[0]) > 0)
factor = (rotD[0] - ro[0]) / rd[0];
else if (std::abs(rd[1]) > 0)
factor = (rotD[1] - ro[1]) / rd[1];
else if (std::abs(rd[2]) > 0)
factor = (rotD[2] - ro[2] / rd[2]);
CHECK((ro[0] + factor * rd[0]) == Approx(rotD[0]));
CHECK((ro[1] + factor * rd[1]) == Approx(rotD[1]));
CHECK((ro[2] + factor * rd[2]) == Approx(rotD[2]));
}
}
}
}
}
#include "CircleTrajectoryGenerator.h"
#include "Logger.h"
#include "VolumeDescriptor.h"
#include "PlanarDetectorDescriptor.h"
#include <stdexcept>
namespace elsa
{
std::pair<std::vector<Geometry>, std::unique_ptr<DataDescriptor>>
CircleTrajectoryGenerator::createTrajectory(index_t numberOfPoses,
const DataDescriptor& volumeDescriptor,
index_t arcDegrees, real_t sourceToCenter,
real_t centerToDetector)
std::unique_ptr<DetectorDescriptor> CircleTrajectoryGenerator::createTrajectory(
index_t numberOfPoses, const DataDescriptor& volumeDescriptor, index_t arcDegrees,
real_t sourceToCenter, real_t centerToDetector)
{
// pull in geometry namespace, to reduce cluttering
using namespace geometry;
......@@ -37,8 +36,6 @@ namespace elsa
volumeDescriptor.getSpacingPerDimension()[1], 1;
}
VolumeDescriptor sinoDescriptor(coeffs, spacing);
std::vector<Geometry> geometryList;
real_t angleIncrement = static_cast<real_t>(1.0f) * static_cast<real_t>(arcDegrees)
......@@ -64,7 +61,7 @@ namespace elsa
}
}
return std::make_pair(geometryList, sinoDescriptor.clone());
return std::make_unique<PlanarDetectorDescriptor>(coeffs, spacing, geometryList);
}
} // namespace elsa