artreconstructor.h

#ifndef CTL_ARTRECONSTRUCTOR_H
#define CTL_ARTRECONSTRUCTOR_H

#include "abstractreconstructor.h"
#include "abstractsubsetgenerator.h"
#include "acquisition/acquisitionsetup.h"
#include "processing/abstractvolumefilter.h"
#include "projectors/abstractprojector.h"
#include "defaultsubsetgenerator.h"

#include <memory>

namespace CTL {

class ARTReconstructor : public AbstractReconstructor
{
    CTL_TYPE_ID(200)
    Q_GADGET

public:
    enum StoppingCriterion
    {
        NoStoppingCriterion     = 0,

        MaximumNbIterations     = 1 << 0,
        MaximumTime             = 1 << 1,
        ProjectionErrorChange   = 1 << 2,
        VolumeDomainChange      = 1 << 3,
        RelativeProjectionError = 1 << 4,
        NormalEquationSatisfied = 1 << 5,

        AllStoppingCriteria = MaximumNbIterations |
                              MaximumTime |
                              ProjectionErrorChange |
                              VolumeDomainChange |
                              RelativeProjectionError |
                              NormalEquationSatisfied
    };
    Q_DECLARE_FLAGS(StoppingCriteria, StoppingCriterion)
    Q_FLAG(StoppingCriteria)

    ARTReconstructor();
    ARTReconstructor(std::unique_ptr<AbstractProjector> forwardProjector,
                     std::unique_ptr<AbstractReconstructor> backprojector);
    ARTReconstructor(AbstractProjector* forwardProjector,
                     AbstractReconstructor* backprojector);

    public: void configure(const AcquisitionSetup& setup) override;
    public: bool reconstructToPlain(const ProjectionDataView& projections,
                                    VoxelVolume<float>& targetVolume) override;
    public: bool reconstructToSpectral(const ProjectionDataView& projections,
                                       SpectralVolumeData& targetVolume) override;
public:
    // SerializationInterface interface
    QVariant parameter() const override;
    void setParameter(const QVariant& parameter) override;
    QVariant toVariant() const override;

    static float computeRelaxationEstimate(const AcquisitionSetup& setup,
                                           const VoxelVolume<float>& targetVolume,
                                           std::shared_ptr<AbstractProjector> forwardProjector,
                                           std::shared_ptr<AbstractReconstructor> backprojector,
                                           uint nbPowerIter = 1u);
    static float computeRelaxationEstimate(const AcquisitionSetup& setup,
                                           const VoxelVolume<float>& targetVolume,
                                           AbstractProjector* forwardProjector = nullptr,
                                           AbstractReconstructor* backprojector = nullptr,
                                           uint nbPowerIter = 1u);

    // setter
    void setPositivityConstraintEnabled(bool enabled);
    void setRegularizationEnabled(bool enabled);
    void setRelaxationEstimationEnabled(bool enabled);

    // relaxation setters -> deactivate automatic estimation when used
    void setRelaxation(float relax);
    void setRelaxationByEstimation(const VoxelVolume<float>& targetVolume,
                                   uint nbPowerIter = 1u);
    void setRelaxationByEstimation(const AcquisitionSetup& setup,
                                   const VoxelVolume<float>& targetVolume,
                                   uint nbPowerIter = 1u);

    void setRegularizer(AbstractVolumeFilter* regularizer, bool enableRegularization = true);
    void setRegularizer(std::unique_ptr<AbstractVolumeFilter> regularizer,
                        bool enableRegularization = true);
    void setSubsetGenerator(AbstractSubsetGenerator* generator);
    void setSubsetGenerator(std::unique_ptr<AbstractSubsetGenerator> generator);
    void setForwardProjector(std::unique_ptr<AbstractProjector> projector);
    void setBackprojector(std::unique_ptr<AbstractReconstructor> projector);
    void setForwardProjector(AbstractProjector* projector);
    void setBackprojector(AbstractReconstructor* projector);
    void useCustomSubsetGenerator();
    void useDefaultSubsetGenerator();

    // stopping criteria
    void setMaxNbIterations(uint nbIt, bool enableCriterion = true);
    void setMaxTime(float seconds, bool enableCriterion = true);
    void setMinChangeInVolumeDomain(float minRelativeChange, bool enableCriterion = true);
    void setMinChangeInProjectionError(float minRelativeChange, bool enableCriterion = true);
    void setMinRelativeProjectionError(float minRelativeError, bool enableCriterion = true);
    void setNormalEqTolerance(float relativeTol, bool enableCriterion = true);
    void setStoppingCriteria(int criteria);
    void setStoppingCriterionEnabled(StoppingCriterion criterion, bool enabled = true);


    // getter
    bool isCustomSubsetGeneratorInUse() const;
    bool isDefaultSubsetGeneratorInUse() const;
    bool isPositivityConstraintEnabled() const;
    bool isRelaxationEstimationEnabled() const;
    bool isRegularizationEnabled() const;
    bool isStoppingCriterionEnabled(int criteria) const;
    float relaxation() const;
    int stoppingCriteria() const;
    AbstractVolumeFilter& regularizer() const;
    AbstractSubsetGenerator& customSubsetGenerator() const;
    DefaultSubsetGenerator& defaultSubsetGenerator() const;
    AbstractSubsetGenerator& subsetGenerator() const;

private:
    AcquisitionSetup _setup;
    std::unique_ptr<AbstractProjector> _fp;
    std::unique_ptr<AbstractReconstructor> _bp;
    std::unique_ptr<AbstractVolumeFilter> _regularizer;
    std::unique_ptr<AbstractSubsetGenerator> _customSubsetGen;
    std::unique_ptr<DefaultSubsetGenerator> _defaultSubsetGen;
    AbstractSubsetGenerator* _activeSubsetGen;

    // parameters
    uint _maxNbIterations         = 5;
    bool _usePositivityConstraint = true;
    bool _useRegularization       = false;
    bool _useRelaxationEstimation = true;
    float _relax                  = 1.0e-6f; 
    float _terminateNormalEqTol   = 0.001f;  
    float _terminateVolChange     = 0.001f;  
    float _terminateProjErrChange = 0.01f;   
    float _terminateProjErrRel    = 0.01f;   
    uint  _terminateMaxTime       = 60000u;  
    int _stoppingCriteria = MaximumNbIterations;

    bool consistencyChecks(const ProjectionDataView& projections) const;

    // stopping criteria methods
    bool hasAnyStoppingCriterion() const;
    bool stoppingCriterionReached(const VoxelVolume<float>& oldVol,
                                  const VoxelVolume<float>& newVol,
                                  float oldProjError,
                                  float newProjError,
                                  float inputProjNorm,
                                  float normalEqNormalization,
                                  int msSpent,
                                  uint iteration) const;
    bool terminateByComputationTime(int msSpent) const;
    bool terminateByIterationCount(uint iteration) const;
    bool terminateByProjError(float projError, float inputProjNorm) const;
    bool terminateByProjErrorChange(float oldProjError, float newProjError) const;
    bool terminateByVolumeChange(const VoxelVolume<float>& oldVol,
                                 const VoxelVolume<float>& newVol) const;
    bool terminateByNormalEqTol(const VoxelVolume<float>& oldVol,
                                const VoxelVolume<float>& newVol,
                                float normalization) const;

    bool reconstructData(const ProjectionDataView& projections,
                         SpectralVolumeData& targetVolume,
                         AbstractProjector& fp,
                         AbstractReconstructor& bp);

    static float computeRelaxationEstimate(const AcquisitionSetup& setup,
                                           const VoxelVolume<float>& targetVolume,
                                           AbstractProjector& fp,
                                           AbstractReconstructor& bp,
                                           uint nbPowerIter);
    static QMetaEnum metaEnum();
};

Q_DECLARE_OPERATORS_FOR_FLAGS(ARTReconstructor::StoppingCriteria)


class ARTReconstructorSFP : public ARTReconstructor
{
public:
    enum FootprintType { TR, TT, TT_Generic };
    explicit ARTReconstructorSFP(FootprintType footprintType = TR);
    ARTReconstructorSFP(FootprintType footprintTypeFP, FootprintType footprintTypeBP);
};

namespace assist {
void enforcePositivity(VoxelVolume<float>& volume);
} // namespace assist
} // namespace CTL

#endif // CTL_ARTRECONSTRUCTOR_H