The name of the initial branch for new projects is now "main" instead of "master". Existing projects remain unchanged. More information: https://doku.lrz.de/display/PUBLIC/GitLab

Commit cb408faa authored by Michael Loipführer's avatar Michael Loipführer
Browse files

add dummy tests for subset problems and sampler, fix some issues

parent ebd501c8
Pipeline #458020 failed with stages
in 25 minutes and 15 seconds
......@@ -15,6 +15,7 @@ set(MODULE_HEADERS
QuadricProblem.h
LASSOProblem.h
SplittingProblem.h
SubsetProblem.h
WLSSubsetProblem.h)
# list all the code files of the module
......@@ -26,6 +27,7 @@ set(MODULE_SOURCES
QuadricProblem.cpp
LASSOProblem.cpp
SplittingProblem.cpp
SubsetProblem.cpp
WLSSubsetProblem.cpp)
......
......@@ -4,10 +4,11 @@ namespace elsa
{
template <typename data_t>
WLSSubsetProblem<data_t>::WLSSubsetProblem(LinearOperator<data_t> A,
WLSSubsetProblem<data_t>::WLSSubsetProblem(LinearOperator<data_t> A, DataContainer<data_t> b,
std::vector<LinearOperator<data_t>> subsetAs,
DataContainer<data_t> b)
: SubsetProblem<data_t>(WLSProblem<data_t>(A, b), *wslProblemsFromOperators(subsetAs, b))
DataContainer<data_t> subsetBs)
: SubsetProblem<data_t>(WLSProblem<data_t>(A, b),
*wslProblemsFromOperators(subsetAs, subsetBs))
{
}
......
......@@ -25,11 +25,14 @@ namespace elsa
* \brief Constructor for a subset problem
*
* \param[in] A the full system matrix of the whole WSL problem
* \param[in] b data vector
* \param[in] subsetAs the system matrices corresponding to each subset
* \param[in] b data vector with a BlockDescriptor for each subset
* \param[in] subsetBs a data vector with a Block Descriptor containing the data vector vor
* each subset
*/
WLSSubsetProblem(LinearOperator<data_t> A, std::vector<LinearOperator<data_t>> subsetAs,
DataContainer<data_t> b);
WLSSubsetProblem(LinearOperator<data_t> A, DataContainer<data_t> b,
std::vector<LinearOperator<data_t>> subsetAs,
DataContainer<data_t> subsetBs);
/// default destructor
~WLSSubsetProblem() override = default;
......
......@@ -9,4 +9,6 @@ ELSA_TEST(WLSProblem)
ELSA_TEST(QuadricProblem)
ELSA_TEST(TikhonovProblem)
ELSA_TEST(LASSOProblem)
ELSA_TEST(SplittingProblem)
\ No newline at end of file
ELSA_TEST(SplittingProblem)
ELSA_TEST(SubsetProblem)
ELSA_TEST(WLSSubsetProblem)
/**
* \file test_SubsetProblem.cpp
*
* \brief Tests for the SubsetProblem class
*
* \author Michael Loipführer - initial code
*/
#include <catch2/catch.hpp>
#include "SubsetProblem.h"
using namespace elsa;
TEMPLATE_TEST_CASE("Scenario: Testing SubsetProblem", "", float, double)
{
}
/**
* \file test_WLSSubsetProblem.cpp
*
* \brief Tests for the WLSSubsetProblem class
*
* \author Michael Loipführer - initial code
*/
#include <catch2/catch.hpp>
#include "WLSSubsetProblem.h"
using namespace elsa;
TEMPLATE_TEST_CASE("Scenario: Testing WLSSubsetProblem", "", float, double)
{
}
......@@ -34,7 +34,7 @@ set(MODULE_SOURCES
add_library(${ELSA_MODULE_TARGET_NAME} ${MODULE_HEADERS} ${MODULE_SOURCES})
add_library(elsa::${ELSA_MODULE_NAME} ALIAS ${ELSA_MODULE_TARGET_NAME})
target_link_libraries(${ELSA_MODULE_TARGET_NAME} PUBLIC elsa_core elsa_logging elsa_problems elsa_proximity_operators)
target_link_libraries(${ELSA_MODULE_TARGET_NAME} PUBLIC elsa_core elsa_logging elsa_problems elsa_proximity_operators elsa_projectors)
target_include_directories(${ELSA_MODULE_TARGET_NAME}
PUBLIC
......
#include "SubsetSampler.h"
#include "RandomBlocksDescriptor.h"
#include "PlanarDetectorDescriptor.h"
#include "SiddonsMethod.h"
namespace elsa
{
......@@ -7,18 +9,48 @@ namespace elsa
SubsetSampler<detectorDescriptor_t, data_t>::SubsetSampler(
const VolumeDescriptor& volumeDescriptor, const detectorDescriptor_t& detectorDescriptor,
const DataContainer<data_t>& sinogram, index_t nSubsets)
: _nSubsets{nSubsets}, _fullDetectorDescriptor(detectorDescriptor)
// : _fullDataDescriptor(sinogram.getDataDescriptor),
: _volumeDescriptor(volumeDescriptor),
_fullDetectorDescriptor(detectorDescriptor),
_nSubsets{nSubsets}
{
// the individual data descriptors for each block
std::vector<DataDescriptor> dataDescriptors;
std::vector<std::unique_ptr<DataDescriptor>> dataDescriptors;
std::vector<IndexVector_t> subsetIndices;
// determine the mapping of indices to subsets
index_t subsetSize = static_cast<index_t>(
sinogram.getDataDescriptor().getNumberOfCoefficientsPerDimension()[0] / _nSubsets);
for (index_t i = 0; i < _nSubsets - 1; i++) {
IndexVector_t indices(subsetSize);
for (index_t j = 0; j < subsetSize; j++) {
indices[j] = i + j * _nSubsets;
}
subsetIndices.emplace_back(indices);
}
index_t lastSubsetSize =
sinogram.getDataDescriptor().getNumberOfCoefficientsPerDimension()[0]
- subsetSize * (_nSubsets - 1);
IndexVector_t lastSubsetIndices(lastSubsetSize);
for (index_t j = 0; j < subsetSize; j++) {
lastSubsetIndices[j] = _nSubsets - 1 + j * _nSubsets;
}
// TODO: this is not quite right, better would be for the first subset to be larger
lastSubsetIndices[lastSubsetIndices.size() - 1] =
sinogram.getDataDescriptor().getNumberOfCoefficientsPerDimension()[0] - 1;
subsetIndices.emplace_back(lastSubsetIndices);
for (index_t i = 0; i < _nSubsets; i++) {
IndexVector_t indices =
subsetIndices[i]; // the indices of the data rows belonging to this subset
std::vector<Geometry> geometry;
for (auto index : indices) {
geometry.emplace_back(detectorDescriptor.getGeometryAt(index));
auto geo = detectorDescriptor.getGeometryAt(index);
if (geo.has_value()) {
geometry.emplace_back(*geo);
}
}
IndexVector_t numOfCoeffsPerDim =
detectorDescriptor.getNumberOfCoefficientsPerDimension();
......@@ -34,35 +66,45 @@ namespace elsa
// and then replacing the data descriptor or create a full new DataContainer and move
// data to it row by row (would also need support in the DataContainer itself to some
// degree)
dataDescriptors.emplace_back(DataDescriptor(numOfCoeffsPerDim));
VolumeDescriptor dataDescriptor(numOfCoeffsPerDim);
dataDescriptors.emplace_back(std::make_unique<VolumeDescriptor>(dataDescriptor));
}
RandomBlocksDescriptor dataDescriptor{dataDescriptors};
_data(dataDescriptor, sinogram.getDataHandlerType());
RandomBlocksDescriptor dataDescriptor(dataDescriptors);
// TODO: make this initialization better to not initialize _data twice (if neccessary)
_data =
std::make_unique<DataContainer<data_t>>(dataDescriptor, sinogram.getDataHandlerType());
// resort the actual measurement data
IndexVector_t numOfCoeffsPerDim =
_data.getDataDescriptor().getNumberOfCoefficientsPerDimension();
sinogram.getDataDescriptor().getNumberOfCoefficientsPerDimension();
index_t coeffsPerRow = numOfCoeffsPerDim.tail(numOfCoeffsPerDim.size() - 1).prod();
for (index_t i = 0; i < _nSubsets; i++) {
IndexVector_t indices =
subsetIndices[i]; // the indices of the data rows belonging to this subset
// the indices of the data rows belonging to this subset
IndexVector_t indices = subsetIndices[i];
// TODO: directly perform this operation on the linearized data
auto block = (*_data).getBlock(i);
auto coeffsPerColumnBlock =
block.getDataDescriptor().getNumberOfCoefficientsPerDimension()[0];
// FIXME: index computation foo
// TODO: determine if this approach works for >= 3 dimensions
for (int j = 0; j < indices.size(); j++) {
for (int k = 0; i < coeffsPerRow; k++) {
for (int k = 0; k < coeffsPerRow; k++) {
index_t rowIndex = numOfCoeffsPerDim[0] * k;
_data[i * indices.size() * coeffsPerRow + j + rowIndex] =
sinogram[indices[j] + rowIndex];
block[j + coeffsPerColumnBlock * k] = sinogram[indices[j] + rowIndex];
}
}
}
}
template <typename detectorDescriptor_t, typename data_t>
DataContainer<data_t> SubsetSampler<detectorDescriptor_t, data_t>::getData()
DataContainer<data_t> SubsetSampler<detectorDescriptor_t, data_t>::getData(bool viewAsBlocks)
{
return _data;
// if (viewAsBlocks)
return (*_data);
// return (*_data).viewAs(_fullDataDescriptor);
}
template <typename detectorDescriptor_t, typename data_t>
......@@ -78,7 +120,7 @@ namespace elsa
if (i < 0 || i >= _nSubsets) {
throw std::invalid_argument("SubsetSampler: index out of bounds for number of subsets");
}
return _data.getBlock(i);
return (*_data).getBlock(i);
}
template <typename detectorDescriptor_t, typename data_t>
......@@ -92,6 +134,19 @@ namespace elsa
return projector_t(_volumeDescriptor, _detectorDescriptors[i]);
}
template <typename detectorDescriptor_t, typename data_t>
template <typename projector_t>
std::vector<projector_t> SubsetSampler<detectorDescriptor_t, data_t>::getSubsetProjectors()
{
std::vector<projector_t> projectors;
for (const auto& detectorDescriptor : _detectorDescriptors) {
projectors.emplace_back(projector_t(_volumeDescriptor, detectorDescriptor));
}
return projectors;
}
template <typename detectorDescriptor_t, typename data_t>
template <typename projector_t>
std::pair<projector_t, DataContainer<data_t>>
......@@ -102,6 +157,13 @@ namespace elsa
"SubsetSampler: index out of bounds for number of projectors");
}
return std::pair(projector_t(_volumeDescriptor, _detectorDescriptors[i]),
_data.getBlock(i));
(*_data).getBlock(i));
}
// ------------------------------------------
// explicit template instantiation
template class SubsetSampler<PlanarDetectorDescriptor, float>;
template class SubsetSampler<PlanarDetectorDescriptor, double>;
template class SubsetSampler<PlanarDetectorDescriptor, std::complex<float>>;
template class SubsetSampler<PlanarDetectorDescriptor, std::complex<double>>;
} // namespace elsa
\ No newline at end of file
......@@ -12,6 +12,7 @@ namespace elsa
*
* \author Michael Loipführer - initial code
*
* \tparam detectorDescriptor_t
* \tparam data_t data type for the domain and range of the problem, defaulting to real_t
*
*
......@@ -41,7 +42,7 @@ namespace elsa
/**
* \brief return the full sinogram
*/
DataContainer<data_t> getData();
DataContainer<data_t> getData(bool viewAsBlocks);
/**
* \brief return the full projector
......@@ -51,36 +52,33 @@ namespace elsa
/**
* \brief return the i-th subset generated by the sampler
*
* \param i is the index of the subset to be returned
*
* \returns the i-th subset of the measurement data
*/
DataContainer<data_t> getSubsetData(index_t i);
/**
* \brief return the i-th projector generated by the sampler
*
* \param i is the index of the projector to be returned
*
* \returns the i-th projector
*/
template <typename projector_t>
projector_t getSubsetProjector(index_t i);
/**
* \brief return a list of projectors that correspond to each subset
*/
template <typename projector_t>
std::vector<projector_t> getSubsetProjectors();
/**
* \brief return both the i-th projector and subset generated by the sampler
*
* \param i is the index of the projector and subset to be returned
*
* \returns pair containing the i-th projector and corresponding subset
*/
template <typename projector_t>
std::pair<projector_t, DataContainer<data_t>> getSubsetProjectorAndData(index_t i);
private:
/// measurement data split into subsets using a RandomBlockDescriptor
DataContainer<data_t> _data;
std::unique_ptr<DataContainer<data_t>> _data;
// /// Data descriptor of the full data
// DataDescriptor _fullDataDescriptor;
/// volume descriptor of the problem
VolumeDescriptor _volumeDescriptor;
......
......@@ -10,6 +10,9 @@
#include <iostream>
#include "SQS.h"
#include "WLSProblem.h"
#include "WLSSubsetProblem.h"
#include "SubsetSampler.h"
#include "PlanarDetectorDescriptor.h"
#include "Problem.h"
#include "Identity.h"
#include "LinearResidual.h"
......@@ -179,3 +182,63 @@ TEST_CASE("Scenario: Solving a simple phantom reconstruction", "")
}
}
}
TEST_CASE("Scenario: Solving a simple phantom problem using ordered subsets", "")
{
// eliminate the timing info from console for the tests
Logger::setLevel(Logger::LogLevel::INFO);
GIVEN("a Phantom reconstruction problem")
{
IndexVector_t size(2);
size << 16, 16; // TODO: determine optimal phantom size for efficient testing
auto phantom = PhantomGenerator<real_t>::createModifiedSheppLogan(size);
auto& volumeDescriptor = phantom.getDataDescriptor();
index_t numAngles{90}, arc{360};
auto sinoDescriptor = CircleTrajectoryGenerator::createTrajectory(
numAngles, phantom.getDataDescriptor(), arc, size(0) * 100, size(0));
SiddonsMethod projector(static_cast<const VolumeDescriptor&>(volumeDescriptor),
*sinoDescriptor);
auto sinogram = projector.apply(phantom);
index_t nSubsets{4};
SubsetSampler<PlanarDetectorDescriptor, real_t> subsetSampler(
static_cast<const VolumeDescriptor&>(volumeDescriptor),
static_cast<const PlanarDetectorDescriptor&>(*sinoDescriptor), sinogram, nSubsets);
subsetSampler.getSubsetProjectors<SiddonsMethod>();
WLSSubsetProblem problem(
subsetSampler.getProjector<SiddonsMethod>(), subsetSampler.getData(),
subsetSampler.getSubsetProjectors<SiddonsMethod>(), subsetSampler.getData(true));
real_t epsilon = std::numeric_limits<real_t>::epsilon();
WHEN("setting up a SQS solver")
{
SQS solver{problem, true, epsilon};
THEN("the clone works correctly")
{
auto sqsClone = solver.clone();
REQUIRE(sqsClone.get() != &solver);
REQUIRE(*sqsClone == solver);
AND_THEN("it works as expected")
{
auto reconstruction = solver.solve(40);
DataContainer resultsDifference = reconstruction - phantom;
// should have converged for the given number of iterations
// does not converge to the optimal solution because of the regularization term
REQUIRE(Approx(resultsDifference.squaredL2Norm()).margin(1)
<= epsilon * epsilon * phantom.squaredL2Norm());
}
}
}
}
}
......@@ -13,6 +13,8 @@
#include "SiddonsMethod.h"
#include "PlanarDetectorDescriptor.h"
#include "PhantomGenerator.h"
#include "SiddonsMethod.h"
#include "JosephsMethod.h"
using namespace elsa;
......@@ -21,27 +23,47 @@ SCENARIO("Testing SubsetSampler with PlanarDetectorDescriptor")
Logger::setLevel(Logger::LogLevel::WARN);
IndexVector_t size(2);
size << 13, 24;
VolumeDescriptor dd{size};
Eigen::Matrix<real_t, -1, 1> bVec(dd.getNumberOfCoefficients());
bVec.setRandom();
DataContainer<real_t> dcB{dd, bVec};
size << 128, 128;
auto phantom = PhantomGenerator<real_t>::createModifiedSheppLogan(size);
auto& volumeDescriptor = phantom.getDataDescriptor();
index_t numAngles{180}, arc{360};
auto sinoDescriptor = CircleTrajectoryGenerator::createTrajectory(
numAngles, dcB.getDataDescriptor(), arc, size(0) * 100, size(0));
numAngles, phantom.getDataDescriptor(), arc, size(0) * 100, size(0));
SiddonsMethod projector(static_cast<const VolumeDescriptor&>(volumeDescriptor),
*sinoDescriptor);
auto sinogram = projector.apply(phantom);
const auto coeffsPerDimSinogram =
sinogram.getDataDescriptor().getNumberOfCoefficientsPerDimension();
GIVEN("A SubsetSampler with 4 subsets")
{
index_t nSubsets{4};
SubsetSampler<PlanarDetectorDescriptor, real_t> subsetSampler(
dd, static_cast<const PlanarDetectorDescriptor&>(*sinoDescriptor), dcB, nSubsets);
static_cast<const VolumeDescriptor&>(volumeDescriptor),
static_cast<const PlanarDetectorDescriptor&>(*sinoDescriptor), sinogram, nSubsets);
WHEN("Returning the individual subsets and corresponding projectors")
{
THEN("The correct ones are returned") {}
THEN("The correct ones are returned")
{
for (index_t i = 0; i < nSubsets; i++) {
const auto subset = subsetSampler.getSubsetData(i);
const auto coeffsPerDimension =
subset.getDataDescriptor().getNumberOfCoefficientsPerDimension();
REQUIRE(coeffsPerDimension[0] == coeffsPerDimSinogram[0] / nSubsets);
REQUIRE(coeffsPerDimension[1] == coeffsPerDimSinogram[1]);
}
subsetSampler.getProjector<SiddonsMethod>();
subsetSampler.getProjector<JosephsMethod>();
subsetSampler.getSubsetProjectors<SiddonsMethod>();
subsetSampler.getSubsetProjectors<JosephsMethod>();
}
}
}
}
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment