Commit 7004ac58 authored by Nikola Dinev's avatar Nikola Dinev Committed by Tobias Lasser

Protect copy constructors in non-leaf and leaf classes, to avoid the...

Protect copy constructors in non-leaf and leaf classes, to avoid the implicitly generated ones by the compiler that could cause slicing. In Cloneable, the assignment operator was deleted to disable issues with implicitly generated assingment operators in derived classes to prevent slicing. (resolves #31)
parent aea3cb41
Pipeline #190113 passed with stages
in 7 minutes and 20 seconds
......@@ -31,6 +31,9 @@ namespace elsa
bool operator!=(const Derived& other) const { return !(*this == other); }
/// delete implicitly declared copy assignment to prevent copy assignment of derived classes
Cloneable& operator=(const Cloneable&) = delete;
protected:
/// actual clone implementation method, abstract to force override in derived classes
virtual Derived* cloneImpl() const = 0;
......@@ -39,10 +42,8 @@ namespace elsa
virtual bool isEqual(const Derived& other) const = 0;
/// default copy constructor, protected to not be publicly available (but available for
/// cloneImpl)
/// cloneImpl()
Cloneable(const Cloneable&) = default;
/// default copy assignment, protected to not be publicly available
Cloneable& operator=(const Cloneable&) = default;
};
} // namespace elsa
......@@ -105,6 +105,9 @@ namespace elsa
/// vector containing the origin of the described volume (typically the center)
RealVector_t _locationOfOrigin;
/// default copy constructor, hidden from non-derived classes to prevent potential slicing
DataDescriptor(const DataDescriptor&) = default;
/// implement the polymorphic clone operation
DataDescriptor* cloneImpl() const override;
......
......@@ -68,6 +68,9 @@ namespace elsa
~EmissionLogLikelihood() override = default;
protected:
/// default copy constructor, hidden from non-derived classes to prevent potential slicing
EmissionLogLikelihood(const EmissionLogLikelihood<data_t>&) = default;
/// the evaluation of the emission log-likelihood
data_t evaluateImpl(const DataContainer<data_t>& Rx) override;
......
......@@ -46,6 +46,9 @@ namespace elsa
~Huber() override = default;
protected:
/// default copy constructor, hidden from non-derived classes to prevent potential slicing
Huber(const Huber<data_t>&) = default;
/// the evaluation of the Huber norm
data_t evaluateImpl(const DataContainer<data_t>& Rx) override;
......
......@@ -42,6 +42,9 @@ namespace elsa
~L1Norm() override = default;
protected:
/// default copy constructor, hidden from non-derived classes to prevent potential slicing
L1Norm(const L1Norm<data_t>&) = default;
/// the evaluation of the l1 norm
data_t evaluateImpl(const DataContainer<data_t>& Rx) override;
......
......@@ -41,6 +41,9 @@ namespace elsa
~L2NormPow2() override = default;
protected:
/// default copy constructor, hidden from non-derived classes to prevent potential slicing
L2NormPow2(const L2NormPow2<data_t>&) = default;
/// the evaluation of the l2 norm (squared)
data_t evaluateImpl(const DataContainer<data_t>& Rx) override;
......
......@@ -42,6 +42,9 @@ namespace elsa
~LInfNorm() override = default;
protected:
/// default copy constructor, hidden from non-derived classes to prevent potential slicing
LInfNorm(const LInfNorm<data_t>&) = default;
/// the evaluation of the linf norm
data_t evaluateImpl(const DataContainer<data_t>& Rx) override;
......
......@@ -65,6 +65,9 @@ namespace elsa
const DataContainer<data_t>& getDataVector() const;
protected:
/// default copy constructor, hidden from non-derived classes to prevent potential slicing
LinearResidual(const LinearResidual<data_t>&) = default;
/// implement the polymorphic clone operation
LinearResidual<data_t>* cloneImpl() const override;
......
......@@ -47,6 +47,9 @@ namespace elsa
~PseudoHuber() override = default;
protected:
/// default copy constructor, hidden from non-derived classes to prevent potential slicing
PseudoHuber(const PseudoHuber<data_t>&) = default;
/// the evaluation of the Huber norm
data_t evaluateImpl(const DataContainer<data_t>& Rx) override;
......
......@@ -67,6 +67,9 @@ namespace elsa
const LinearResidual<data_t>& getGradientExpression() const;
protected:
/// default copy constructor, hidden from non-derived classes to prevent potential slicing
Quadric(const Quadric<data_t>&) = default;
/// the evaluation of the Quadric functional
data_t evaluateImpl(const DataContainer<data_t>& Rx) override;
......
......@@ -77,6 +77,9 @@ namespace elsa
~TransmissionLogLikelihood() override = default;
protected:
/// default copy constructor, hidden from non-derived classes to prevent potential slicing
TransmissionLogLikelihood(const TransmissionLogLikelihood<data_t>&) = default;
/// the evaluation of the transmission log-likelihood
data_t evaluateImpl(const DataContainer<data_t>& Rx) override;
......
......@@ -46,6 +46,9 @@ namespace elsa
const Scaling<data_t>& getWeightingOperator() const;
protected:
/// default copy constructor, hidden from non-derived classes to prevent potential slicing
WeightedL2NormPow2(const WeightedL2NormPow2<data_t>&) = default;
/// the evaluation of the weighted, squared l2 norm
data_t evaluateImpl(const DataContainer<data_t>& Rx) override;
......
......@@ -5,9 +5,11 @@
namespace elsa
{
std::pair<std::vector<Geometry>, DataDescriptor> CircleTrajectoryGenerator::createTrajectory(
index_t numberOfPoses, const DataDescriptor& volumeDescriptor, index_t arcDegrees,
real_t sourceToCenter, real_t centerToDetector)
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)
{
// sanity check
auto dim = volumeDescriptor.getNumberOfDimensions();
......@@ -48,7 +50,7 @@ namespace elsa
}
}
return std::make_pair(geometryList, sinoDescriptor);
return std::make_pair(geometryList, sinoDescriptor.clone());
}
} // namespace elsa
......@@ -36,7 +36,7 @@ namespace elsa
*
* Please note: the sinogram size/spacing will match the volume size/spacing.
*/
static std::pair<std::vector<Geometry>, DataDescriptor>
static std::pair<std::vector<Geometry>, std::unique_ptr<DataDescriptor>>
createTrajectory(index_t numberOfPoses, const DataDescriptor& volumeDescriptor,
index_t arcDegrees, real_t sourceToCenter, real_t centerToDetector);
};
......
......@@ -44,7 +44,7 @@ SCENARIO("Create a Circular Trajectory")
real_t angle = (1.0 / (numberOfAngles - 1)) * halfCircular;
for (int i = 0; i < numberOfAngles; ++i) {
real_t currAngle = i * angle * pi / 180.0;
Geometry tmpGeom(sourceToCenter, centerToDetector, currAngle, desc, sdesc);
Geometry tmpGeom(sourceToCenter, centerToDetector, currAngle, desc, *sdesc);
REQUIRE((tmpGeom.getCameraCenter() - geomList[i].getCameraCenter()).norm()
== Approx(0));
......@@ -77,7 +77,7 @@ SCENARIO("Create a Circular Trajectory")
real_t angle = (1.0 / (numberOfAngles - 1)) * halfCircular;
for (int i = 0; i < numberOfAngles; ++i) {
real_t currAngle = i * angle * pi / 180.0;
Geometry tmpGeom(sourceToCenter, centerToDetector, currAngle, desc, sdesc);
Geometry tmpGeom(sourceToCenter, centerToDetector, currAngle, desc, *sdesc);
REQUIRE((tmpGeom.getCameraCenter() - geomList[i].getCameraCenter()).norm()
== Approx(0));
......@@ -118,7 +118,7 @@ SCENARIO("Create a Circular Trajectory")
real_t angleInc = 1.0 * halfCircular / (numberOfAngles - 1);
for (int i = 0; i < numberOfAngles; ++i) {
real_t angle = i * angleInc * pi / 180.0;
Geometry tmpGeom(sourceToCenter, centerToDetector, desc, sdesc, angle);
Geometry tmpGeom(sourceToCenter, centerToDetector, desc, *sdesc, angle);
REQUIRE((tmpGeom.getCameraCenter() - geomList[i].getCameraCenter()).norm()
== Approx(0));
......@@ -151,7 +151,7 @@ SCENARIO("Create a Circular Trajectory")
real_t angleInc = 1.0 * halfCircular / (numberOfAngles - 1);
for (int i = 0; i < numberOfAngles; ++i) {
real_t angle = i * angleInc * pi / 180.0;
Geometry tmpGeom(sourceToCenter, centerToDetector, desc, sdesc, angle);
Geometry tmpGeom(sourceToCenter, centerToDetector, desc, *sdesc, angle);
REQUIRE((tmpGeom.getCameraCenter() - geomList[i].getCameraCenter()).norm()
== Approx(0));
......
......@@ -20,7 +20,7 @@ namespace elsa
auto [descriptor, dataType] = parseHeader(properties);
// read in the data
DataContainer<data_t> dataContainer(descriptor);
DataContainer<data_t> dataContainer(*descriptor);
if (dataType == DataUtils::DataType::UINT16)
DataUtils::parseRawData<uint16_t, data_t>(file, dataContainer);
......@@ -117,7 +117,7 @@ namespace elsa
return properties;
}
std::pair<DataDescriptor, DataUtils::DataType>
std::pair<std::unique_ptr<DataDescriptor>, DataUtils::DataType>
EDF::parseHeader(const std::map<std::string, std::string>& properties)
{
// read the dimensions
......@@ -208,10 +208,8 @@ namespace elsa
dimSpacingVec[i] = spacing[i];
}
// the data descriptor condensed from the info
DataDescriptor dataDescriptor(dimSizeVec, dimSpacingVec);
return std::make_pair(dataDescriptor, dataType);
return std::make_pair(std::make_unique<DataDescriptor>(dimSizeVec, dimSpacingVec),
dataType);
}
template <typename data_t>
......@@ -225,7 +223,7 @@ namespace elsa
file << "ByteOrder = LowByteFirst;\n";
file << "DataType = " << getDataTypeName(data) << ";\n";
auto descriptor = data.getDataDescriptor();
auto& descriptor = data.getDataDescriptor();
// write dimension and size
for (std::size_t i = 0; i < descriptor.getNumberOfDimensions(); ++i)
......
......@@ -43,7 +43,7 @@ namespace elsa
static std::map<std::string, std::string> readHeader(std::ifstream& file);
/// parse the EDF header property map into a DataDescriptor and DataType
static std::pair<DataDescriptor, DataUtils::DataType>
static std::pair<std::unique_ptr<DataDescriptor>, DataUtils::DataType>
parseHeader(const std::map<std::string, std::string>& properties);
/// write the EDF header to file
......
......@@ -25,7 +25,7 @@ namespace elsa
throw std::runtime_error("MHD::read: can not read from '" + dataPath + "'");
// read in the data
DataContainer<data_t> dataContainer(descriptor);
DataContainer<data_t> dataContainer(*descriptor);
if (dataType == DataUtils::DataType::UINT16)
DataUtils::parseRawData<uint16_t, data_t>(dataFile, dataContainer);
......@@ -97,7 +97,7 @@ namespace elsa
return properties;
}
std::tuple<DataDescriptor, std::string, DataUtils::DataType>
std::tuple<std::unique_ptr<DataDescriptor>, std::string, DataUtils::DataType>
MHD::parseHeader(const std::map<std::string, std::string>& properties)
{
// check the dimensions
......@@ -183,17 +183,15 @@ namespace elsa
dimSpacing[i] = dimSpacingVec[i];
}
// the data descriptor condensed form the info
DataDescriptor dataDescriptor(dimSizes, dimSpacing);
return std::make_tuple(dataDescriptor, rawDataPath, dataType);
return std::make_tuple(std::make_unique<DataDescriptor>(dimSizes, dimSpacing), rawDataPath,
dataType);
}
template <typename data_t>
void MHD::writeHeader(std::ofstream& metaFile, const DataContainer<data_t>& data,
std::string rawFilename)
{
auto descriptor = data.getDataDescriptor();
auto& descriptor = data.getDataDescriptor();
// write dimension, size and spacing
metaFile << "NDims = " << descriptor.getNumberOfDimensions() << "\n";
......
......@@ -46,7 +46,7 @@ namespace elsa
static std::map<std::string, std::string> readHeader(std::ifstream& metaFile);
/// parse the MHD header property map into a DataDescriptor and DataType
static std::tuple<DataDescriptor, std::string, DataUtils::DataType>
static std::tuple<std::unique_ptr<DataDescriptor>, std::string, DataUtils::DataType>
parseHeader(const std::map<std::string, std::string>& properties);
/// write the MHD header to file
......
......@@ -64,6 +64,9 @@ namespace elsa
~FiniteDifferences() override = default;
protected:
/// default copy constructor, hidden from non-derived classes to prevent potential slicing
FiniteDifferences(const FiniteDifferences<data_t>&) = default;
/// apply the finite differences operator
void applyImpl(const DataContainer<data_t>& x, DataContainer<data_t>& Ax) const override;
......
......@@ -29,6 +29,9 @@ namespace elsa
~Identity() override = default;
protected:
/// default copy constructor, hidden from non-derived classes to prevent potential slicing
Identity(const Identity<data_t>&) = default;
/**
* \brief apply the identity operator A to x, i.e. Ax = x
*
......
......@@ -50,6 +50,9 @@ namespace elsa
const DataContainer<data_t>& getScaleFactors() const;
protected:
/// default copy constructor, hidden from non-derived classes to prevent potential slicing
Scaling(const Scaling<data_t>&) = default;
/// apply the scaling operation
void applyImpl(const DataContainer<data_t>& x, DataContainer<data_t>& Ax) const override;
......
......@@ -55,6 +55,9 @@ namespace elsa
~BinaryMethod() override = default;
protected:
/// default copy constructor, hidden from non-derived classes to prevent potential slicing
BinaryMethod(const BinaryMethod<data_t>&) = default;
/// apply the binary method (i.e. forward projection)
void applyImpl(const DataContainer<data_t>& x, DataContainer<data_t>& Ax) const override;
......
......@@ -62,6 +62,9 @@ namespace elsa
~JosephsMethod() = default;
protected:
/// default copy constructor, hidden from non-derived classes to prevent potential slicing
JosephsMethod(const JosephsMethod<data_t>&) = default;
/// apply Joseph's method (i.e. forward projection)
void applyImpl(const DataContainer<data_t>& x, DataContainer<data_t>& Ax) const override;
......
......@@ -51,6 +51,9 @@ namespace elsa
~SiddonsMethod() override = default;
protected:
/// default copy constructor, hidden from non-derived classes to prevent potential slicing
SiddonsMethod(const SiddonsMethod<data_t>&) = default;
/// apply Siddon's method (i.e. forward projection)
void applyImpl(const DataContainer<data_t>& x, DataContainer<data_t>& Ax) const override;
......
......@@ -277,7 +277,7 @@ SCENARIO("Testing BinaryMethod")
THEN("Domain descriptor is still the same")
{
auto retDescriptor = op.getDomainDescriptor();
auto& retDescriptor = op.getDomainDescriptor();
CHECK(retDescriptor.getNumberOfCoefficientsPerDimension()(0) == 5);
CHECK(retDescriptor.getNumberOfCoefficientsPerDimension()(1) == 5);
......@@ -285,7 +285,7 @@ SCENARIO("Testing BinaryMethod")
THEN("Domain descriptor is still the same")
{
auto retDescriptor = op.getRangeDescriptor();
auto& retDescriptor = op.getRangeDescriptor();
CHECK(retDescriptor.getNumberOfCoefficientsPerDimension()(0) == 5);
CHECK(retDescriptor.getNumberOfCoefficientsPerDimension()(1) == 1);
......
......@@ -64,6 +64,9 @@ namespace elsa
virtual ~JosephsMethodCUDA();
protected:
/// default copy constructor, hidden from non-derived classes to prevent potential slicing
JosephsMethodCUDA(const JosephsMethodCUDA<data_t>&) = default;
/// apply Joseph's method (i.e. forward projection)
void applyImpl(const DataContainer<data_t>& x, DataContainer<data_t>& Ax) const override;
......
......@@ -55,6 +55,9 @@ namespace elsa
~SiddonsMethodCUDA();
protected:
/// default copy constructor, hidden from non-derived classes to prevent potential slicing
SiddonsMethodCUDA(const SiddonsMethodCUDA<data_t>&) = default;
/// apply Siddon's method (i.e. forward projection)
void applyImpl(const DataContainer<data_t>& x, DataContainer<data_t>& Ax) const override;
......
......@@ -64,6 +64,10 @@ namespace elsa
/// default destructor
~CG() override = default;
protected:
/// default copy constructor, hidden from non-derived classes to prevent potential slicing
CG(const CG<data_t>&) = default;
private:
/// the default number of iterations
const index_t _defaultIterations{100};
......
......@@ -30,6 +30,10 @@ namespace elsa
/// default destructor
~GradientDescent() override = default;
protected:
/// default copy constructor, hidden from non-derived classes to prevent potential slicing
GradientDescent(const GradientDescent<data_t>&) = default;
private:
/// the step size
real_t _stepSize;
......
......@@ -23,7 +23,7 @@ if(ELSA_BUILD_CUDA_PROJECTORS)
# build the GPU projector speed test program
add_executable(speed_test speed_test.cpp)
target_link_libraries(speed_test elsa_generators elsa_projectors_cuda elsa_logging)
target_link_libraries(speed_test elsa::all)
target_compile_features(speed_test PUBLIC cxx_std_17)
add_dependencies(examples speed_test)
endif()
......
......@@ -9,7 +9,8 @@ using namespace elsa;
void example2d()
{
// generate 2d phantom
IndexVector_t size(2); size << 128, 128;
IndexVector_t size(2);
size << 128, 128;
auto phantom = PhantomGenerator<real_t>::createModifiedSheppLogan(size);
// write the phantom out
......@@ -17,12 +18,12 @@ void example2d()
// generate circular trajectory
index_t noAngles{100}, arc{360};
auto [geometry, sinoDescriptor] = CircleTrajectoryGenerator::createTrajectory(noAngles, phantom.getDataDescriptor(),
arc, size(0)*100, size(0));
auto [geometry, sinoDescriptor] = CircleTrajectoryGenerator::createTrajectory(
noAngles, phantom.getDataDescriptor(), arc, size(0) * 100, size(0));
// setup operator for 2d X-ray transform
Logger::get("Info")->info("Simulating sinogram using Siddon's method");
SiddonsMethod projector(phantom.getDataDescriptor(), sinoDescriptor, geometry);
SiddonsMethod projector(phantom.getDataDescriptor(), *sinoDescriptor, geometry);
// simulate the sinogram
auto sinogram = projector.apply(phantom);
......@@ -30,28 +31,26 @@ void example2d()
// write the sinogram out
EDF::write(sinogram, "2dsinogram.edf");
// setup reconstruction problem
WLSProblem problem(projector, sinogram);
// solve the reconstruction problem
GradientDescent solver(problem, 1.0/size.prod());
GradientDescent solver(problem, 1.0 / size.prod());
index_t noIterations{50};
Logger::get("Info")->info("Solving reconstruction using {} iterations of gradient descent", noIterations);
Logger::get("Info")->info("Solving reconstruction using {} iterations of gradient descent",
noIterations);
auto reconstruction = solver.solve(noIterations);
// write the reconstruction out
EDF::write(reconstruction, "2dreconstruction.edf");
}
int main()
{
try {
example2d();
}
catch (std::exception& e) {
} catch (std::exception& e) {
std::cerr << "An exception occurred: " << e.what() << "\n";
}
}
......@@ -6,10 +6,11 @@
using namespace elsa;
void example2d()
void example3d()
{
// generate 3d phantom
IndexVector_t size(3); size << 128, 128, 128;
IndexVector_t size(3);
size << 128, 128, 128;
auto phantom = PhantomGenerator<real_t>::createModifiedSheppLogan(size);
// write the phantom out
......@@ -17,12 +18,12 @@ void example2d()
// generate circular trajectory
index_t noAngles{100}, arc{360};
auto [geometry, sinoDescriptor] = CircleTrajectoryGenerator::createTrajectory(noAngles, phantom.getDataDescriptor(),
arc, size(0)*100, size(0));
auto [geometry, sinoDescriptor] = CircleTrajectoryGenerator::createTrajectory(
noAngles, phantom.getDataDescriptor(), arc, size(0) * 100, size(0));
// setup operator for 2d X-ray transform
Logger::get("Info")->info("Simulating sinogram using Siddon's method");
JosephsMethodCUDA projector(phantom.getDataDescriptor(), sinoDescriptor, geometry);
JosephsMethodCUDA projector(phantom.getDataDescriptor(), *sinoDescriptor, geometry);
// simulate the sinogram
auto sinogram = projector.apply(phantom);
......@@ -30,28 +31,26 @@ void example2d()
// write the sinogram out
EDF::write(sinogram, "3dsinogram.edf");
// setup reconstruction problem
WLSProblem problem(projector, sinogram);
// solve the reconstruction problem
GradientDescent solver(problem, 1.0/size.prod());
GradientDescent solver(problem, 1.0 / size.prod());
index_t noIterations{50};
Logger::get("Info")->info("Solving reconstruction using {} iterations of gradient descent", noIterations);
Logger::get("Info")->info("Solving reconstruction using {} iterations of gradient descent",
noIterations);
auto reconstruction = solver.solve(noIterations);
// write the reconstruction out
EDF::write(reconstruction, "3dreconstruction.edf");
}
int main()
{
try {
example2d();
}
catch (std::exception& e) {
example3d();
} catch (std::exception& e) {
std::cerr << "An exception occurred: " << e.what() << "\n";
}
}
/// Elsa example program: test execution speed of GPU projectors
#include "CircleTrajectoryGenerator.h"
#include "PhantomGenerator.h"
#include "Logger.h"
#include "JosephsMethodCUDA.h"
#include "SiddonsMethodCUDA.h"
#include "elsa.h"
#include <chrono>
#include <iostream>
......@@ -13,7 +8,8 @@
using namespace elsa;
void testExecutionSpeed(LinearOperator<real_t>& projector, DataContainer<real_t>& volume, int numIters)
void testExecutionSpeed(LinearOperator<real_t>& projector, DataContainer<real_t>& volume,
int numIters)
{
real_t timer = 0;
......@@ -22,30 +18,30 @@ void testExecutionSpeed(LinearOperator<real_t>& projector, DataContainer<real_t>
DataContainer backproj(projector.getDomainDescriptor());
// run test for forward projection
for (int i = 0 ; i < numIters; ++i) {
for (int i = 0; i < numIters; ++i