Commit e3b20eda authored by David Frank's avatar David Frank
Browse files

Add global ray type alias

parent 476b1c97
......@@ -8,7 +8,7 @@ namespace elsa
: DataDescriptor(numOfCoeffsPerDim), _geometry(geometryList)
{
// TODO Clarify: What about empty geometryList? Do we want to support it, or throw an
// execption?
// exception?
}
DetectorDescriptor::DetectorDescriptor(const IndexVector_t& numOfCoeffsPerDim,
......@@ -18,8 +18,7 @@ namespace elsa
{
}
DetectorDescriptor::Ray
DetectorDescriptor::computeRayFromDetectorCoord(const index_t detectorIndex) const
RealRay_t DetectorDescriptor::computeRayFromDetectorCoord(const index_t detectorIndex) const
{
// Return empty, if access out of bounds
......@@ -31,8 +30,7 @@ namespace elsa
return computeRayFromDetectorCoord(coord);
}
DetectorDescriptor::Ray
DetectorDescriptor::computeRayFromDetectorCoord(const IndexVector_t coord) const
RealRay_t DetectorDescriptor::computeRayFromDetectorCoord(const IndexVector_t coord) const
{
// Assume all of the coordinates are inside of the volume
// auto tmp = (coord.array() < getNumberOfCoefficientsPerDimension().array());
......
......@@ -17,9 +17,6 @@ namespace elsa
*/
class DetectorDescriptor : public DataDescriptor
{
public:
using Ray = Eigen::ParametrizedLine<real_t, Eigen::Dynamic>;
public:
/// There is not default signal
DetectorDescriptor() = delete;
......@@ -45,14 +42,14 @@ namespace elsa
* @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;
RealRay_t 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;
RealRay_t computeRayFromDetectorCoord(const IndexVector_t coord) const;
/**
* @brief Compute a ray from the source from a pose to the given detector coordinate
......@@ -62,13 +59,13 @@ namespace elsa
* @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;
virtual RealRay_t 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,
virtual RealVector_t computeDetectorCoordFromRay(const RealRay_t& ray,
const index_t poseIndex) const = 0;
/// Get the number of poses used in the geometry
......
......@@ -15,13 +15,13 @@ namespace elsa
{
}
DetectorDescriptor::Ray
RealRay_t
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
// as this is the hot path, I don't want exceptions 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()
......@@ -53,11 +53,11 @@ namespace elsa
.head(dim) // Transform to non-homogeneous
.normalized(); // normalize vector
return Ray(ro, rd);
return RealRay_t(ro, rd);
}
RealVector_t
PlanarDetectorDescriptor::computeDetectorCoordFromRay(const Ray& ray,
PlanarDetectorDescriptor::computeDetectorCoordFromRay(const RealRay_t& ray,
const index_t poseIndex) const
{
auto dim = getNumberOfDimensions();
......
......@@ -2,6 +2,7 @@
#include "DetectorDescriptor.h"
#include "TypeCasts.hpp"
#include "elsaDefines.h"
namespace elsa
{
......@@ -14,8 +15,6 @@ namespace elsa
*/
class PlanarDetectorDescriptor : public DetectorDescriptor
{
using DetectorDescriptor::Ray;
public:
PlanarDetectorDescriptor() = delete;
......@@ -39,10 +38,10 @@ namespace elsa
using DetectorDescriptor::computeRayFromDetectorCoord;
/// Override function to compute rays for a planar detector
Ray computeRayFromDetectorCoord(const RealVector_t& detectorCoord,
const index_t poseIndex) const override;
RealRay_t computeRayFromDetectorCoord(const RealVector_t& detectorCoord,
const index_t poseIndex) const override;
RealVector_t computeDetectorCoordFromRay(const Ray& ray,
RealVector_t computeDetectorCoordFromRay(const RealRay_t& ray,
const index_t poseIndex) const override;
private:
......
......@@ -3,6 +3,7 @@
#include <complex>
#include <cstddef>
#include <Eigen/Core>
#include <Eigen/Geometry>
#include <type_traits>
#ifdef ELSA_CUDA_VECTOR
......@@ -36,6 +37,13 @@ namespace elsa
/// global type for matrices of real numbers
using RealMatrix_t = Eigen::Matrix<real_t, Eigen::Dynamic, Eigen::Dynamic>;
/// global type alias for rays
using RealRay_t = Eigen::ParametrizedLine<real_t, Eigen::Dynamic>;
/// global type alias for rays
template <typename data_t>
using Ray_t = Eigen::ParametrizedLine<data_t, Eigen::Dynamic>;
/// template global constexpr for the number pi
template <typename T>
constexpr auto pi = static_cast<T>(3.14159265358979323846);
......
......@@ -15,8 +15,6 @@ using namespace elsa;
using namespace elsa::geometry;
using namespace doctest;
using Ray = DetectorDescriptor::Ray;
TEST_SUITE_BEGIN("core");
TEST_CASE("PlanarDetectorDescriptor: Testing 2D PlanarDetectorDescriptor")
......@@ -88,7 +86,7 @@ TEST_CASE("PlanarDetectorDescriptor: Testing 2D PlanarDetectorDescriptor")
rd << -0.141421f, 0.989949f;
rd.normalize();
auto detCoord = desc.computeDetectorCoordFromRay(Ray(ro, rd), 0);
auto detCoord = desc.computeDetectorCoordFromRay(RealRay_t(ro, rd), 0);
CHECK_EQ(detCoord[0], Approx(0.5).epsilon(0.05));
}
THEN("The detector coord is correct asd")
......@@ -96,7 +94,7 @@ TEST_CASE("PlanarDetectorDescriptor: Testing 2D PlanarDetectorDescriptor")
rd << 0.f, 1.f;
rd.normalize();
auto detCoord = desc.computeDetectorCoordFromRay(Ray(ro, rd), 0);
auto detCoord = desc.computeDetectorCoordFromRay(RealRay_t(ro, rd), 0);
CHECK_EQ(detCoord[0], Approx(2.5));
}
......@@ -105,7 +103,7 @@ TEST_CASE("PlanarDetectorDescriptor: Testing 2D PlanarDetectorDescriptor")
rd << 0.141421f, 0.989949f;
rd.normalize();
auto detCoord = desc.computeDetectorCoordFromRay(Ray(ro, rd), 0);
auto detCoord = desc.computeDetectorCoordFromRay(RealRay_t(ro, rd), 0);
CHECK_EQ(detCoord[0], Approx(4.5).epsilon(0.05));
}
}
......
......@@ -91,9 +91,6 @@ namespace elsa
void traverseVolume(const DataContainer<data_t>& vector,
DataContainer<data_t>& result) const;
/// convenience typedef for ray
using Ray = Eigen::ParametrizedLine<real_t, Eigen::Dynamic>;
/// lift from base class
using LinearOperator<data_t>::_domainDescriptor;
using LinearOperator<data_t>::_rangeDescriptor;
......
......@@ -2,7 +2,8 @@
namespace elsa
{
std::optional<IntersectionResult> Intersection::withRay(const BoundingBox& aabb, const Ray& r)
std::optional<IntersectionResult> Intersection::withRay(const BoundingBox& aabb,
const RealRay_t& r)
{
real_t invDir = 1 / r.direction()(0);
......@@ -28,5 +29,4 @@ namespace elsa
return std::make_optional<IntersectionResult>(tmin, tmax);
return std::nullopt;
}
} // namespace elsa
......@@ -39,10 +39,6 @@ namespace elsa
*/
class Intersection
{
private:
/// the type for a ray using Eigen
using Ray = Eigen::ParametrizedLine<real_t, Eigen::Dynamic>;
public:
/**
* @brief Compute entry and exit point of ray in a volume (given as an AABB)
......@@ -59,7 +55,8 @@ namespace elsa
* Method adapted from
https://tavianator.com/fast-branchless-raybounding-box-intersections-part-2-nans/
*/
static std::optional<IntersectionResult> withRay(const BoundingBox& aabb, const Ray& r);
static std::optional<IntersectionResult> withRay(const BoundingBox& aabb,
const RealRay_t& r);
};
/**
......
......@@ -97,9 +97,6 @@ namespace elsa
void traverseVolume(const DataContainer<data_t>& vector,
DataContainer<data_t>& result) const;
/// convenience typedef for ray
using Ray = Eigen::ParametrizedLine<real_t, Eigen::Dynamic>;
/**
* @brief Linear interpolation, works in any dimension
*
......
......@@ -83,9 +83,6 @@ namespace elsa
void traverseVolume(const DataContainer<data_t>& vector,
DataContainer<data_t>& result) const;
/// convenience typedef for ray
using Ray = Eigen::ParametrizedLine<real_t, Eigen::Dynamic>;
/// lift from base class
using LinearOperator<data_t>::_domainDescriptor;
using LinearOperator<data_t>::_rangeDescriptor;
......
......@@ -3,7 +3,7 @@
namespace elsa
{
TraverseAABB::TraverseAABB(const BoundingBox& aabb, const Ray& r) : _aabb{aabb}
TraverseAABB::TraverseAABB(const BoundingBox& aabb, const RealRay_t& r) : _aabb{aabb}
{
// compute the first intersection
calculateAABBIntersections(r);
......@@ -62,7 +62,7 @@ namespace elsa
return _currentPos.template cast<IndexVector_t::Scalar>();
}
void TraverseAABB::calculateAABBIntersections(const TraverseAABB::Ray& r)
void TraverseAABB::calculateAABBIntersections(const RealRay_t& r)
{
// entry and exit point parameters
real_t tmin;
......
......@@ -22,10 +22,6 @@ namespace elsa
*/
class TraverseAABB
{
private:
/// convenience type def for ray
using Ray = Eigen::ParametrizedLine<real_t, Eigen::Dynamic>;
public:
/**
* @brief Constructor for traversal, accepting bounding box and ray
......@@ -33,7 +29,7 @@ namespace elsa
* @param[in] aabb axis-aligned boundary box describing the volume
* @param[in] r the ray to be traversed
*/
TraverseAABB(const BoundingBox& aabb, const Ray& r);
TraverseAABB(const BoundingBox& aabb, const RealRay_t& r);
/**
* @brief Update the traverser status by taking the next traversal step
......@@ -88,7 +84,7 @@ namespace elsa
RealVector_t(_aabb._dim).setConstant(std::numeric_limits<real_t>::max())};
/// compute the entry and exit points of ray r with the volume (aabb)
void calculateAABBIntersections(const Ray& r);
void calculateAABBIntersections(const RealRay_t& r);
/// setup the step directions (which is basically the sign of the ray direction rd)
void initStepDirection(const RealVector_t& rd);
/// select the closest voxel to the entry point
......
......@@ -3,12 +3,13 @@
namespace elsa
{
TraverseAABBJosephsMethod::TraverseAABBJosephsMethod(const BoundingBox& aabb, const Ray& r)
TraverseAABBJosephsMethod::TraverseAABBJosephsMethod(const BoundingBox& aabb,
const RealRay_t& r)
: _aabb{aabb}
{
initStepDirection(r.direction());
Ray rt(r.origin(), r.direction());
RealRay_t rt(r.origin(), r.direction());
// determinge length of entire intersection with AABB
real_t intersectionLength = calculateAABBIntersections(rt);
......@@ -108,7 +109,7 @@ namespace elsa
// determine main direction
rd.cwiseAbs().maxCoeff(&_ignoreDirection);
// intialize _nextStep as the step for interior pixels
// initialize _nextStep as the step for interior pixels
_nextStep(_ignoreDirection) = _stepDirection(_ignoreDirection);
for (int i = 0; i < _aabb._dim; ++i) {
if (i != _ignoreDirection) {
......@@ -185,7 +186,7 @@ namespace elsa
_stepDirection = rd.array().sign();
}
real_t TraverseAABBJosephsMethod::calculateAABBIntersections(const Ray& ray)
real_t TraverseAABBJosephsMethod::calculateAABBIntersections(const RealRay_t& ray)
{
real_t tmin;
......
......@@ -15,9 +15,6 @@ namespace elsa
class TraverseAABBJosephsMethod
{
private:
/// convenience type def for ray
using Ray = Eigen::ParametrizedLine<real_t, Eigen::Dynamic>;
/**
* @brief Possible stages of the traversal
*
......@@ -38,7 +35,7 @@ namespace elsa
* @param[in] aabb axis-aligned boundary box describing the volume
* @param[in] r the ray to be traversed
*/
explicit TraverseAABBJosephsMethod(const BoundingBox& aabb, const Ray& r);
explicit TraverseAABBJosephsMethod(const BoundingBox& aabb, const RealRay_t& r);
/**
* @brief Update the traverser status by taking the next traversal step
......@@ -124,7 +121,7 @@ namespace elsa
void initStepDirection(const RealVector_t& rd);
/// compute the entry and exit points of ray r with the volume (aabb), returns the length of
/// the intersection
real_t calculateAABBIntersections(const Ray& ray);
real_t calculateAABBIntersections(const RealRay_t& ray);
/// select the closest voxel to the entry point (considering the current position)
void selectClosestVoxel(const RealVector_t& currentPosition);
/// check if the current index is still in the bounding box
......
......@@ -15,8 +15,6 @@
using namespace elsa;
using namespace doctest;
using Ray = Eigen::ParametrizedLine<real_t, Eigen::Dynamic>;
TEST_CASE("Intersection: Intersect corners of pixels")
{
size_t dim = 2;
......@@ -31,7 +29,7 @@ TEST_CASE("Intersection: Intersect corners of pixels")
RealVector_t rd(dim);
rd << 1.0, 1.0;
rd.normalize();
Ray r(ro, rd);
RealRay_t r(ro, rd);
REQUIRE_UNARY(Intersection::withRay(aabb, r));
......@@ -39,7 +37,7 @@ TEST_CASE("Intersection: Intersect corners of pixels")
ro << 1, 2;
rd << 1.0, -1.0;
rd.normalize();
r = Ray(ro, rd);
r = RealRay_t(ro, rd);
REQUIRE_UNARY_FALSE(Intersection::withRay(aabb, r));
......@@ -47,7 +45,7 @@ TEST_CASE("Intersection: Intersect corners of pixels")
ro << 3, -2;
rd << -1.0, 1.0;
rd.normalize();
r = Ray(ro, rd);
r = RealRay_t(ro, rd);
REQUIRE_UNARY(Intersection::withRay(aabb, r));
......@@ -55,7 +53,7 @@ TEST_CASE("Intersection: Intersect corners of pixels")
ro << 3, 1;
rd << -1.0, -1.0;
rd.normalize();
r = Ray(ro, rd);
r = RealRay_t(ro, rd);
REQUIRE_UNARY_FALSE(Intersection::withRay(aabb, r));
}
......@@ -69,7 +67,7 @@ TEST_CASE("Intersection: Intersect edges of voxels")
ro << 132, 30;
RealVector_t rd(dim);
rd << -1.0, 0;
Ray r(ro, rd);
RealRay_t r(ro, rd);
IndexVector_t voxel(dim);
......@@ -98,7 +96,7 @@ TEST_CASE("Intersection: Intersect edges of voxels")
// vertical check
ro << 30, -35;
rd << 0, 1.0;
r = Ray(ro, rd);
r = RealRay_t(ro, rd);
voxel << 31, 30;
THEN("the ray intersects")
......@@ -123,7 +121,7 @@ TEST_CASE("Intersection: Intersect edges of voxels")
rd << 0.0, 1.0;
ro << 1.5, -10;
r = Ray(ro, rd);
r = RealRay_t(ro, rd);
voxel << 0, 1;
THEN("the ray does not intersect")
......
......@@ -16,9 +16,7 @@
using namespace elsa;
using namespace doctest;
using Ray = Eigen::ParametrizedLine<real_t, Eigen::Dynamic>;
bool intersect(const RealVector_t& voxel, const Ray& r)
bool intersect(const RealVector_t& voxel, const RealRay_t& r)
{
// pre-check parallel rays
for (index_t i = 0; i < r.dim(); ++i) {
......@@ -64,7 +62,7 @@ TEST_CASE("TraverseAABB: Construction of a 2D traversal object")
{
ro << 0.5, -0.5;
rd << 0.0, 1.0;
Ray r(ro, rd);
RealRay_t r(ro, rd);
TraverseAABB traverse(aabb, r);
CHECK_UNARY(traverse.isInBoundingBox());
......@@ -80,7 +78,7 @@ TEST_CASE("TraverseAABB: Construction of a 2D traversal object")
{
ro << 1.0, -0.5;
rd << 0.0, 1.0;
Ray r(ro, rd);
RealRay_t r(ro, rd);
TraverseAABB traverse(aabb, r);
CHECK_UNARY(traverse.isInBoundingBox());
......@@ -96,7 +94,7 @@ TEST_CASE("TraverseAABB: Construction of a 2D traversal object")
{
ro << 1.5, -0.5;
rd << 0.0, 1.0;
Ray r(ro, rd);
RealRay_t r(ro, rd);
TraverseAABB traverse(aabb, r);
CHECK_UNARY(traverse.isInBoundingBox());
......@@ -112,7 +110,7 @@ TEST_CASE("TraverseAABB: Construction of a 2D traversal object")
{
ro << 2.0, -0.5;
rd << 0.0, 1.0;
Ray r(ro, rd);
RealRay_t r(ro, rd);
TraverseAABB traverse(aabb, r);
CHECK_UNARY(traverse.isInBoundingBox());
......@@ -128,7 +126,7 @@ TEST_CASE("TraverseAABB: Construction of a 2D traversal object")
{
ro << 2.5, -0.5;
rd << 0.0, 1.0;
Ray r(ro, rd);
RealRay_t r(ro, rd);
TraverseAABB traverse(aabb, r);
CHECK_UNARY(traverse.isInBoundingBox());
......@@ -147,7 +145,7 @@ TEST_CASE("TraverseAABB: Construction of a 2D traversal object")
{
ro << -0.5, 0.5;
rd << 1.0, 0.0;
Ray r(ro, rd);
RealRay_t r(ro, rd);
TraverseAABB traverse(aabb, r);
CHECK_UNARY(traverse.isInBoundingBox());
......@@ -163,7 +161,7 @@ TEST_CASE("TraverseAABB: Construction of a 2D traversal object")
{
ro << -0.5, 1.0;
rd << 1.0, 0.0;
Ray r(ro, rd);
RealRay_t r(ro, rd);
TraverseAABB traverse(aabb, r);
CHECK_UNARY(traverse.isInBoundingBox());
......@@ -179,7 +177,7 @@ TEST_CASE("TraverseAABB: Construction of a 2D traversal object")
{
ro << -0.5, 1.5;
rd << 1.0, 0.0;
Ray r(ro, rd);
RealRay_t r(ro, rd);
TraverseAABB traverse(aabb, r);
CHECK_UNARY(traverse.isInBoundingBox());
......@@ -195,7 +193,7 @@ TEST_CASE("TraverseAABB: Construction of a 2D traversal object")
{
ro << -0.5, 2.0;
rd << 1.0, 0.0;
Ray r(ro, rd);
RealRay_t r(ro, rd);
TraverseAABB traverse(aabb, r);
CHECK_UNARY(traverse.isInBoundingBox());
......@@ -211,7 +209,7 @@ TEST_CASE("TraverseAABB: Construction of a 2D traversal object")
{
ro << -0.5, 2.5;
rd << 1.0, 0.0;
Ray r(ro, rd);
RealRay_t r(ro, rd);
TraverseAABB traverse(aabb, r);
CHECK_UNARY(traverse.isInBoundingBox());
......@@ -230,7 +228,7 @@ TEST_CASE("TraverseAABB: Construction of a 2D traversal object")
{
ro << 3.5, 0.5;
rd << -1.0, 0.0;
Ray r(ro, rd);
RealRay_t r(ro, rd);
TraverseAABB traverse(aabb, r);
CHECK_UNARY(traverse.isInBoundingBox());
......@@ -246,7 +244,7 @@ TEST_CASE("TraverseAABB: Construction of a 2D traversal object")
{
ro << 3.5, 1.0;
rd << -1.0, 0.0;
Ray r(ro, rd);
RealRay_t r(ro, rd);
TraverseAABB traverse(aabb, r);
CHECK_UNARY(traverse.isInBoundingBox());
......@@ -262,7 +260,7 @@ TEST_CASE("TraverseAABB: Construction of a 2D traversal object")
{
ro << 3.5, 1.5;
rd << -1.0, 0.0;
Ray r(ro, rd);
RealRay_t r(ro, rd);
TraverseAABB traverse(aabb, r);
CHECK_UNARY(traverse.isInBoundingBox());
......@@ -278,7 +276,7 @@ TEST_CASE("TraverseAABB: Construction of a 2D traversal object")
{
ro << 3.5, 2.0;
rd << -1.0, 0.0;
Ray r(ro, rd);
RealRay_t r(ro, rd);
TraverseAABB traverse(aabb, r);
CHECK_UNARY(traverse.isInBoundingBox());
......@@ -294,7 +292,7 @@ TEST_CASE("TraverseAABB: Construction of a 2D traversal object")
{
ro << 3.5, 2.5;
rd << -1.0, 0.0;
Ray r(ro, rd);
RealRay_t r(ro, rd);
TraverseAABB traverse(aabb, r);
CHECK_UNARY(traverse.isInBoundingBox());
......@@ -313,7 +311,7 @@ TEST_CASE("TraverseAABB: Construction of a 2D traversal object")
{
ro << 0.5, 3.5;
rd << 0.0, -1.0;
Ray r(ro, rd);
RealRay_t r(ro, rd);
TraverseAABB traverse(aabb, r);
CHECK_UNARY(traverse.isInBoundingBox());
......@@ -329,7 +327,7 @@ TEST_CASE("TraverseAABB: Construction of a 2D traversal object")
{
ro << 1.0, 3.5;
rd << 0.0, -1.0;
Ray r(ro, rd);
RealRay_t r(ro, rd);
TraverseAABB traverse(aabb, r);
CHECK_UNARY(traverse.isInBoundingBox());
......@@ -345,7 +343,7 @@ TEST_CASE("TraverseAABB: Construction of a 2D traversal object")
{
ro << 1.5, 3.5;
rd << 0.0, -1.0;
Ray r(ro, rd);
RealRay_t r(ro, rd);
TraverseAABB traverse(aabb, r);
CHECK_UNARY(traverse.isInBoundingBox());
......@@ -361,7 +359,7 @@ TEST_CASE("TraverseAABB: Construction of a 2D traversal object")
{
ro << 2.0, 3.5;
rd << 0.0, -1.0;
Ray r(ro, rd);
RealRay_t r(ro, rd);
TraverseAABB traverse(aabb, r);
CHECK_UNARY(traverse.isInBoundingBox());
......@@ -377,7 +375,7 @@ TEST_CASE("TraverseAABB: Construction of a 2D traversal object")
{