Commit 3554361d authored by David Frank's avatar David Frank
Browse files

Add BlobProjector prototype

parent 3268337b
#include "Blobs.h"
#include "Bessel.h"
namespace elsa
{
template <typename data_t>
constexpr Blob<data_t>::Blob(data_t radius, SelfType_t<data_t> alpha, SelfType_t<data_t> order)
constexpr Blob<data_t>::Blob(data_t radius, SelfType_t<data_t> alpha, int order)
: radius_(radius), alpha_(alpha), order_(order)
{
}
......@@ -27,7 +28,7 @@ namespace elsa
}
template <typename data_t>
constexpr data_t Blob<data_t>::order() const
constexpr int Blob<data_t>::order() const
{
return order_;
}
......@@ -35,19 +36,16 @@ namespace elsa
// explicit template instantiation
namespace blobs
{
template float blob_evaluate<float>(float, SelfType_t<float>, SelfType_t<float>,
SelfType_t<float>);
template double blob_evaluate<double>(double, SelfType_t<double>, SelfType_t<double>,
SelfType_t<double>);
template float blob_evaluate<float>(float, SelfType_t<float>, SelfType_t<float>, int);
template double blob_evaluate<double>(double, SelfType_t<double>, SelfType_t<double>, int);
template float blob_projected<float>(float, SelfType_t<float>, SelfType_t<float>,
SelfType_t<float>);
template double blob_projected<double>(double, SelfType_t<double>, SelfType_t<double>,
SelfType_t<double>);
template float blob_projected<float>(float, SelfType_t<float>, SelfType_t<float>, int);
template double blob_projected<double>(double, SelfType_t<double>, SelfType_t<double>, int);
template float blob_projected<float>(float);
template double blob_projected<double>(double);
} // namespace blobs
template class Blob<float>;
template class Blob<double>;
} // namespace elsa
#pragma once
#include "Error.h"
#include "elsaDefines.h"
#include "Bessel.h"
namespace elsa
{
......@@ -14,15 +16,19 @@ namespace elsa
namespace blobs
{
// The devinition given by Lewitt (1992):
// \f$
// b_{m, alpha, a}(r) = \frac{I_m(\alpha * \sqrt{1 - (r/a)^2)}{I_m(\alpha)} * (\sqrt{1 -
// (r/a)^2})^m \f$ for any \f$ 0 <= r <= a \f$.
template <typename data_t>
constexpr data_t blob_evaluate(data_t r, SelfType_t<data_t> a, SelfType_t<data_t> alpha,
SelfType_t<data_t> m)
int m) noexcept
{
const auto w = static_cast<data_t>(1) - std::pow(r / a, static_cast<data_t>(2));
if (w >= 0) {
const data_t Im1 = std::cyl_bessel_i(m, alpha);
const data_t Im1 = math::bessi(m, alpha);
const data_t arg = std::sqrt(w);
const data_t Im2 = std::cyl_bessel_i(m, alpha * arg);
const data_t Im2 = math::bessi(m, alpha * arg);
return Im2 / Im1 * std::pow(arg, m);
}
......@@ -30,13 +36,35 @@ namespace elsa
}
/// @brief Compute line integral of blob through a straight line
/// The exact formulations are quite easily derived using the fact, that the line
/// integral for a single blob is given by:
///
/// \f$
/// \frac{a}{I_m(\alpha)} \sqrt{\frac{2 \pi}{\alpha}} \left(\sqrt{1 -
/// \left(\frac{s}{a}\right)^2}\right)^{m + 0.5} I_{m + 0.5}\left(\alpha \sqrt{1 -
/// \left(\frac{s}{a}\right)^2}\right)
/// \f$
///
/// Then for a given order substitute I_{m + 0.5}(x) with the elementary function
/// representations, as they can be found in Abramowitz and Stegun's Handbook of
/// mathematical functions (1972):
/// \f$ I_{0.5}(x) = \sqrt{\frac{2}{\pi x}} \sinh(x)\$f
/// \f$ I_{1.5}(x) = \sqrt{\frac{2}{\pi x}} \left( \cosh(x) \frac{\sinh(x)}{x} \right) \$f
/// \f$ I_{2.5}(x) = \sqrt{\frac{2}{\pi x}} \left(\left(\frac{3}{x^2} +
/// \frac{1}{x}\right)\sinh(x) - \frac{3}{x^2} \cosh(x)\right)\$f
///
/// Which will result in the below formulations. In theory using the recurrent relations,
/// this could be extended to further orders, but is deemed unnecessary for now.
///
/// TODO: What about alpha = 0? Check if you find anything in the literature for that.
///
/// @param distance distance of blob center to straight line, in literature often referred
/// to as `r`
/// @param radius radius of blob, often referred to as `a`
/// @param alpha smoothness factor of blob
/// @param alpha smoothness factor of blob, expected to be larger than 0
/// @param order order of Bessel function, often referred to as `m`
///
/// Ref:
/// https://github.com/I2PC/xmipp/blob/3d4cc3f430cbc99a337635edbd95ebbcef33fc44/src/xmipp/libraries/data/blobs.cpp#L91A
/// Distance-Driven Projection and Backprojection for Spherically Symmetric Basis Functions
/// in CT - Levakhina
/// Spherically symmetric volume elements as basis functions for image reconstructions in
......@@ -44,25 +72,32 @@ namespace elsa
/// Semi-Discrete Iteration Methods in X-Ray Tomography - Jonas Vogelgesang
template <typename data_t>
constexpr data_t blob_projected(data_t s, SelfType_t<data_t> a, SelfType_t<data_t> alpha,
SelfType_t<data_t> m)
int m)
{
// Equation derived in Lewitt 1990
const data_t w = static_cast<data_t>(1) - ((s * s) / (a * a));
// If `w` is close to zero or negative, `s` > `a`, and therefore just return 0
if (w > 1e-10) {
const data_t root = std::sqrt(w);
// First three terms of equation
const data_t q1 = a / std::cyl_bessel_i(m, alpha);
const data_t q2 = std::sqrt(2 * pi<data_t> / alpha);
const data_t q3 = std::pow(root, m + static_cast<data_t>(0.5));
const data_t q4 = std::cyl_bessel_i(m + static_cast<data_t>(0.5), alpha * root);
return q1 * q2 * q3 * q4;
// expect alpha > 0
using namespace elsa::math;
const data_t sda = s / a;
const data_t sdas = std::pow(sda, 2);
const data_t w = 1.0 - sdas;
if (w > 1.0e-10) {
const auto arg = alpha * std::sqrt(w);
if (m == 0) {
return (2.0 * a / alpha) * std::sinh(arg) / math::bessi0(alpha);
} else if (m == 1) {
return (2.0 * a / alpha) * std::sqrt(w)
* (std::cosh(arg) - std::sinh(arg) / arg) / math::bessi1(alpha);
} else if (m == 2) {
return (2.0 * a / alpha) * w
* ((3.0 / (arg * arg) + 1.0) * std::sinh(arg) - (3.0 / arg) * cosh(arg))
/ bessi2(alpha);
} else {
throw Error("m out of range in blob_projected()");
}
}
return 0;
return 0.0;
}
template <typename data_t>
......@@ -76,7 +111,7 @@ namespace elsa
class Blob
{
public:
constexpr Blob(data_t radius, SelfType_t<data_t> alpha, SelfType_t<data_t> order);
constexpr Blob(data_t radius, SelfType_t<data_t> alpha, int order);
constexpr data_t operator()(data_t s);
......@@ -84,19 +119,19 @@ namespace elsa
constexpr data_t alpha() const;
constexpr data_t order() const;
constexpr int order() const;
private:
data_t radius_;
data_t alpha_;
data_t order_;
int order_;
};
template <typename data_t>
class ProjectedBlob
{
public:
constexpr ProjectedBlob(data_t radius, SelfType_t<data_t> alpha, SelfType_t<data_t> order)
constexpr ProjectedBlob(data_t radius, SelfType_t<data_t> alpha, int order)
: radius_(radius), alpha_(alpha), order_(order)
{
}
......@@ -105,17 +140,17 @@ namespace elsa
{
return blobs::blob_projected(s, radius_, alpha_, order_);
}
#
data_t radius() const { return radius_; }
data_t alpha() const { return alpha_; }
data_t order() const { return order_; }
int order() const { return order_; }
private:
data_t radius_;
data_t alpha_;
data_t order_;
int order_;
};
} // namespace elsa
......@@ -11,6 +11,11 @@
#include "BoundingBox.h"
#include "Logger.h"
#include "Blobs.h"
#include "CartesianIndices.h"
#include "spdlog/fmt/bundled/core.h"
#include "spdlog/fmt/fmt.h"
#include "spdlog/fmt/ostr.h"
namespace elsa
{
......@@ -60,8 +65,12 @@ namespace elsa
const auto sizeRange = Ax.getSize();
const auto volume_shape = x.getDataDescriptor().getNumberOfCoefficientsPerDimension();
// Loop over all the poses, and for each pose loop over all detector pixels
const IndexVector_t lower = _boundingBox._min.template cast<index_t>();
const IndexVector_t upper = _boundingBox._max.template cast<index_t>();
const auto support = self().support();
#pragma omp parallel for
// Loop over all the poses, and for each pose loop over all detector pixels
for (index_t rangeIndex = 0; rangeIndex < sizeRange; ++rangeIndex) {
// --> get the current ray to the detector center
auto ray = _detectorDescriptor.computeRayFromDetectorCoord(rangeIndex);
......@@ -69,7 +78,7 @@ namespace elsa
index_t leadingdir = 0;
ray.direction().array().cwiseAbs().maxCoeff(&leadingdir);
auto& rangeVal = Ax[rangeIndex];
auto rangeVal = Ax[rangeIndex];
// Expand bounding box as rays have larger support now
auto aabb = _boundingBox;
......@@ -88,10 +97,14 @@ namespace elsa
// Correct position, such that the distance is still correct
auto correctedPos = neighbour.template cast<real_t>().array() + 0.5;
const auto distance = ray.distance(curPos);
const auto weight = self().weight(distance);
rangeVal += weight * x.at(curVoxel);
const auto distance = ray.distance(correctedPos);
const auto weight = self().weight(distance);
rangeVal += weight * x.at(neighbour);
}
}
Ax[rangeIndex] = rangeVal;
}
}
......@@ -106,6 +119,10 @@ namespace elsa
const auto shape = _domainDescriptor->getNumberOfCoefficientsPerDimension();
const IndexVector_t lower = _boundingBox._min.template cast<index_t>();
const IndexVector_t upper = _boundingBox._max.template cast<index_t>();
const auto support = self().support();
#pragma omp parallel for
// Loop over all the poses, and for each pose loop over all detector pixels
for (index_t rangeIndex = 0; rangeIndex < sizeRange; ++rangeIndex) {
......@@ -191,6 +208,8 @@ namespace elsa
data_t weight(data_t distance) const { return lut_(distance); }
index_t support() const { return static_cast<index_t>(std::ceil(lut_.radius())); }
/// implement the polymorphic clone operation
BlobProjector<data_t>* cloneImpl() const override
{
......
......@@ -74,8 +74,10 @@ namespace elsa
const auto b = static_cast<std::size_t>(std::ceil(index));
// Get the function values
const auto fa = data_[a];
const auto fb = data_[b];
const auto fa = a == N ? 0 : data_[a];
const auto fb = b == N ? 0 : data_[b];
auto ret = detail::lerp(fa, fb, index - static_cast<data_t>(a));
// Bilinear interpolation
return detail::lerp(fa, fb, index - static_cast<data_t>(a));
......@@ -99,6 +101,12 @@ namespace elsa
{
}
data_t radius() const { return blob_.radius(); }
data_t alpha() const { return blob_.alpha(); }
data_t order() const { return blob_.order(); }
data_t operator()(data_t distance) const { return lut_((distance / blob_.radius()) * N); }
private:
......
#include "SliceTraversal.h"
#include "Error.h"
#include "elsaDefines.h"
#include "spdlog/fmt/fmt.h"
#include "spdlog/fmt/ostr.h"
......@@ -32,6 +33,49 @@ namespace elsa
return createRotationMatrix(0);
}
RealMatrix_t TransformToTraversal::create3DRotation(const real_t leadingCoeff,
const index_t leadingAxisIndex)
{
auto identity = []() { return RealMatrix_t({{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}); };
auto rotationX = [](real_t radian) {
real_t c = std::cos(radian);
real_t s = std::sin(radian);
return RealMatrix_t({{1, 0, 0}, {0, c, -s}, {0, s, c}});
};
auto rotationY = [](real_t radian) {
real_t c = std::cos(radian);
real_t s = std::sin(radian);
return RealMatrix_t({{c, 0, s}, {0, 1, 0}, {-s, 0, c}});
};
auto rotationZ = [](real_t radian) {
real_t c = std::cos(radian);
real_t s = std::sin(radian);
return RealMatrix_t({{c, -s, 0}, {s, c, 0}, {0, 0, 1}});
};
// TODO: Does a closed form solution exist?
// Use conversion to polar coordinates
if (leadingAxisIndex == 0 && leadingCoeff >= 0) {
// Already identity, so do nothing
return identity();
} else if (leadingAxisIndex == 0 && leadingCoeff < 0) {
return rotationY(pi_t);
} else if (leadingAxisIndex == 1 && leadingCoeff >= 0) {
return rotationZ(-0.5 * pi_t);
} else if (leadingAxisIndex == 1 && leadingCoeff <= 0) {
return rotationZ(0.5 * pi_t);
} else if (leadingAxisIndex == 2 && leadingCoeff >= 0) {
return rotationY(0.5 * pi_t);
} else if (leadingAxisIndex == 2 && leadingCoeff <= 0) {
return rotationY(-0.5 * pi_t);
}
return identity();
}
RealMatrix_t TransformToTraversal::createRotation(RealRay_t ray)
{
// Get the leading axis, by absolute value
......@@ -41,7 +85,13 @@ namespace elsa
// Get the signed leading coefficient
const auto leadingCoeff = ray.direction()(leadingAxisIndex);
return create2DRotation(leadingCoeff, leadingAxisIndex);
if (ray.dim() == 2) {
return create2DRotation(leadingCoeff, leadingAxisIndex);
} else if (ray.dim() == 3) {
return create3DRotation(leadingCoeff, leadingAxisIndex);
}
throw Error("Can not create a {}-dimensional transformation for the slice traversal",
ray.dim());
}
RealMatrix_t TransformToTraversal::createTransformation(const RealRay_t& ray,
......@@ -118,25 +168,9 @@ namespace elsa
SliceTraversal::SliceTraversal(BoundingBox aabb, RealRay_t ray)
: transformation_(ray, aabb._max.array() / 2), ray_(ray)
{
Eigen::IOFormat vecfmt(10, 0, ", ", ", ", "", "", "[", "]");
// {
// auto tmphit = Intersection::xPlanesWithRay(aabb, ray);
//
// fmt::print("aabb: {}\n", aabb);
// fmt::print("ro: {}\n", ray.origin().format(vecfmt));
// fmt::print("rd: {}\n", ray.direction().format(vecfmt));
// if (tmphit) {
// fmt::print("hit: {}\n", *tmphit);
// fmt::print("diff: {}\n", tmphit->_tmax - tmphit->_tmin);
// fmt::print("firstVoxel: {}\n", ray.pointAt(tmphit->_tmin).format(vecfmt));
// fmt::print("lastVoxel: {}\n", ray.pointAt(tmphit->_tmax).format(vecfmt));
// }
// }
// Transform ray and aabb to traversal coordinate space
ray = transformation_.toTraversalCoordinates(ray);
aabb = transformation_.toTraversalCoordinates(aabb);
// fmt::print("aabb: {}\n", aabb);
// TODO: We only want to shift in x direction by 0.5, not in y
auto hit = Intersection::xPlanesWithRay(aabb, ray);
......@@ -146,13 +180,6 @@ namespace elsa
const auto firstVoxel = ray.pointAt(hit->_tmin);
const auto lastVoxel = ray.pointAt(hit->_tmax);
// fmt::print("ro: {}\n", ray.origin().format(vecfmt));
// fmt::print("rd: {}\n", ray.direction().format(vecfmt));
// fmt::print("hit: {}\n", *hit);
// fmt::print("diff: {}\n", hit->_tmax - hit->_tmin);
// fmt::print("firstVoxel: {}\n", firstVoxel.format(vecfmt));
// fmt::print("lastVoxel: {}\n", lastVoxel.format(vecfmt));
endIndex_ = std::max<index_t>(
1, static_cast<index_t>(std::ceil(lastVoxel[0] - firstVoxel[0] + 0.5f)));
......
......@@ -76,6 +76,9 @@ namespace elsa
static RealMatrix_t create2DRotation(const real_t leadingCoeff,
const index_t leadingAxisIndex);
static RealMatrix_t create3DRotation(const real_t leadingCoeff,
const index_t leadingAxisIndex);
static RealMatrix_t createRotation(RealRay_t ray);
static RealMatrix_t createTransformation(const RealRay_t& ray,
......
......@@ -84,7 +84,7 @@ TEST_CASE_TEMPLATE("Blobs: Test projected values", data_t, float, double)
constexpr std::array<data_t, 3> scales{100, 50, 30};
constexpr std::array<data_t, 3> as{1, 2, 3};
constexpr std::array<data_t, 4> alphas{2.5, 5., 7.5, 10.83};
constexpr std::array<index_t, 4> ms{0, 1, 2, 3};
constexpr std::array<index_t, 3> ms{0, 1, 2};
for (const auto a : as) {
const auto scale = scales[static_cast<index_t>(a - 1)];
......
......@@ -307,4 +307,255 @@ TEST_CASE("Intersection: Quick bug tests")
}
}
TEST_CASE("Intersection: Quick bug tests")
{
const IndexVector_t size({{3, 3}});
BoundingBox aabb(size);
aabb._min = RealVector_t({{-1.5f, -1.5f}});
aabb._max = RealVector_t({{1.5f, 1.5f}});
const RealVector_t ro({{-2, -3}});
RealVector_t rd({{1, 1}});
rd.normalize();
const RealRay_t ray(ro, rd);
auto hit = Intersection::xPlanesWithRay(aabb, ray);
CAPTURE(hit);
REQUIRE(hit);
auto dist_to_integer = [&](auto f) {
real_t aabbmin = static_cast<int>(std::round(f));
auto frac = std::abs(f - aabbmin);
return frac;
};
THEN("entries x-coord lies in the middle of a voxel")
{
RealVector_t entry = ray.pointAt(hit->_tmin);
CAPTURE(entry.format(vecfmt));
// real_t trash = 0;
CHECK_EQ(0, doctest::Approx(dist_to_integer(entry[0])));
CHECK((entry.array() < 1.5).all());
CHECK((entry.array() > -1.5).all());
}
THEN("x-coord lies in the middle of a voxel")
{
RealVector_t leave = ray.pointAt(hit->_tmax);
CAPTURE(leave.format(vecfmt));
// real_t trash = 0;
// CHECK_EQ(0, doctest::Approx(std::modf(leave[0], &trash)));
CHECK_EQ(0, doctest::Approx(dist_to_integer(leave[0])));
CHECK((leave.array() < 1.5).all());
CHECK((leave.array() > -1.5).all());
}
}
// Redefine GIVEN such that it's nicely usable inside an loop
#undef GIVEN
#define GIVEN(...) DOCTEST_SUBCASE((std::string(" Given: ") + std::string(__VA_ARGS__)).c_str())
TEST_CASE("Intersection: 3D rays with bounding box")
{
const IndexVector_t upperbound({{5, 5, 5}});
const BoundingBox aabb(upperbound);
for (int z = 0; z < 25; ++z) {
for (int y = 0; y < 25; ++y) {
const auto zinc = static_cast<real_t>(z) / 5.f;
const auto yinc = static_cast<real_t>(y) / 5.f;
GIVEN("ray starting at (6, " + std::to_string(yinc) + ", " + std::to_string(zinc) + ")")
{
const RealVector_t origin({{-1, yinc, zinc}});
RealVector_t dir({{1, 0, 0}});
CAPTURE(yinc);
CAPTURE(zinc);
CAPTURE(origin.format(vecfmt));
CAPTURE(dir.format(vecfmt));
RealRay_t ray(origin, dir);
auto hit = Intersection::xPlanesWithRay(aabb, ray);
CAPTURE(hit);
REQUIRE(hit);
THEN("entries x-coord lies in the middle of a voxel")
{
RealVector_t enter = ray.pointAt(hit->_tmin);
CAPTURE(enter.format(vecfmt));
real_t trash = 0;
CHECK_EQ(0.5, doctest::Approx(std::modf(enter[0], &trash)));
CHECK((enter.array() <= upperbound.template cast<real_t>().array()).all());
CHECK((enter.array() >= 0).all());
}
THEN("x-coord lies in the middle of a voxel")
{
RealVector_t leave = ray.pointAt(hit->_tmax);
CAPTURE(leave.format(vecfmt));
real_t trash = 0;
CHECK_EQ(0.5, doctest::Approx(std::modf(leave[0], &trash)));
CHECK((leave.array() <= upperbound.template cast<real_t>().array()).all());
CHECK((leave.array() >= 0).all());
}
}
}
}
for (int z = 0; z < 25; ++z) {
for (int y = 0; y < 25; ++y) {
const auto zinc = static_cast<real_t>(z) / 5.f;
const auto yinc = static_cast<real_t>(y) / 5.f;
GIVEN("ray starting at (6, " + std::to_string(yinc) + ", " + std::to_string(zinc) + ")")
{
const RealVector_t origin({{6, yinc, zinc}});
RealVector_t dir({{-1, 0, 0}});
CAPTURE(yinc);
CAPTURE(zinc);
CAPTURE(origin.format(vecfmt));
CAPTURE(dir.format(vecfmt));
RealRay_t ray(origin, dir);
auto hit = Intersection::xPlanesWithRay(aabb, ray);
CAPTURE(hit);
REQUIRE(hit);
THEN("entries x-coord lies in the middle of a voxel")
{
RealVector_t enter = ray.pointAt(hit->_tmin);
CAPTURE(enter.format(vecfmt));
real_t trash = 0;
CHECK_EQ(0.5, doctest::Approx(std::modf(enter[0], &trash)));
CHECK((enter.array() <= upperbound.template cast<real_t>().array()).all());
CHECK((enter.array() >= 0).all());
}
THEN("x-coord lies in the middle of a voxel")
{
RealVector_t leave = ray.pointAt(hit->_tmax);
CAPTURE(leave.format(vecfmt));
real_t trash = 0;
CHECK_EQ(0.5, doctest::Approx(std::modf(leave[0], &trash)));
CHECK((leave.array() <= upperbound.template cast<real_t>().array()).all());
CHECK((leave.array() >= 0).all());
}
}
}
}
}