rtkSpectralForwardModelImageFilter.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright RTK Consortium
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         https://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/

#ifndef rtkSpectralForwardModelImageFilter_h
#define rtkSpectralForwardModelImageFilter_h

#include "rtkSchlomka2008NegativeLogLikelihood.h"
#include "rtkDualEnergyNegativeLogLikelihood.h"
#include "rtkVectorImageToImageFilter.h"

#include <itkPermuteAxesImageFilter.h>
#include <itkInPlaceImageFilter.h>
#include <itkCastImageFilter.h>

namespace rtk
{
template <typename DecomposedProjectionsType,
          typename MeasuredProjectionsType,
          typename IncidentSpectrumImageType = itk::Image<float, 3>,
          typename DetectorResponseImageType = itk::Image<float, 2>,
          typename MaterialAttenuationsImageType = itk::Image<float, 2>>
class ITK_TEMPLATE_EXPORT SpectralForwardModelImageFilter
  : public itk::InPlaceImageFilter<MeasuredProjectionsType, MeasuredProjectionsType>
{
public:
  ITK_DISALLOW_COPY_AND_MOVE(SpectralForwardModelImageFilter);

  using Self = SpectralForwardModelImageFilter;
  using Superclass = itk::ImageToImageFilter<MeasuredProjectionsType, MeasuredProjectionsType>;
  using Pointer = itk::SmartPointer<Self>;
  using ConstPointer = itk::SmartPointer<const Self>;

  using InputImageType = MeasuredProjectionsType;
  using OutputImageType = MeasuredProjectionsType;

  using ThresholdsType = itk::VariableLengthVector<double>;
  using DetectorResponseType = vnl_matrix<double>;
  using MaterialAttenuationsType = vnl_matrix<double>;
  using DecomposedProjectionsDataType = typename DecomposedProjectionsType::PixelType::ValueType;
  using MeasuredProjectionsDataType = typename MeasuredProjectionsType::PixelType::ValueType;

#ifndef ITK_FUTURE_LEGACY_REMOVE

  using VectorSpectrumImageType = itk::VectorImage<float, 2>;
  using FlattenVectorFilterType = rtk::VectorImageToImageFilter<VectorSpectrumImageType, IncidentSpectrumImageType>;
  using PermuteFilterType = itk::PermuteAxesImageFilter<IncidentSpectrumImageType>;
#endif

  itkNewMacro(Self);

  itkOverrideGetNameOfClassMacro(SpectralForwardModelImageFilter);

  void
  SetInputIncidentSpectrum(const IncidentSpectrumImageType * IncidentSpectrum);
  void
  SetInputSecondIncidentSpectrum(const IncidentSpectrumImageType * SecondIncidentSpectrum);
#ifndef ITK_FUTURE_LEGACY_REMOVE
  void
  SetInputIncidentSpectrum(const VectorSpectrumImageType * IncidentSpectrum);
  void
  SetInputSecondIncidentSpectrum(const VectorSpectrumImageType * SecondIncidentSpectrum);
#endif
  typename IncidentSpectrumImageType::ConstPointer
  GetInputIncidentSpectrum();
  typename IncidentSpectrumImageType::ConstPointer
  GetInputSecondIncidentSpectrum();

  void
  SetInputDecomposedProjections(
    const typename itk::ImageBase<DecomposedProjectionsType::ImageDimension> * DecomposedProjections);
  template <unsigned int VNumberOfMaterials>
  void
  SetInputFixedVectorLengthDecomposedProjections(
    const itk::Image<itk::Vector<DecomposedProjectionsDataType, VNumberOfMaterials>,
                     DecomposedProjectionsType::ImageDimension> * DecomposedProjections);
  typename DecomposedProjectionsType::ConstPointer
  GetInputDecomposedProjections();

  void
  SetInputMeasuredProjections(
    const typename itk::ImageBase<MeasuredProjectionsType::ImageDimension> * MeasuredProjections);
  template <unsigned int VNumberOfSpectralBins>
  void
  SetInputFixedVectorLengthMeasuredProjections(
    const itk::Image<itk::Vector<MeasuredProjectionsDataType, VNumberOfSpectralBins>,
                     MeasuredProjectionsType::ImageDimension> * MeasuredProjections);
  typename MeasuredProjectionsType::ConstPointer
  GetInputMeasuredProjections();

  void
  SetDetectorResponse(const DetectorResponseImageType * DetectorResponse);
  typename DetectorResponseImageType::ConstPointer
  GetDetectorResponse();

  void
  SetMaterialAttenuations(const MaterialAttenuationsImageType * MaterialAttenuations);
  typename MaterialAttenuationsImageType::ConstPointer
  GetMaterialAttenuations();

  typename DecomposedProjectionsType::ConstPointer
  GetOutputCramerRaoLowerBound();

  typename MeasuredProjectionsType::ConstPointer
  GetOutputVariances();

  itkSetMacro(Thresholds, ThresholdsType);
  itkGetMacro(Thresholds, ThresholdsType);

  itkSetMacro(NumberOfSpectralBins, unsigned int);
  itkGetMacro(NumberOfSpectralBins, unsigned int);

  itkSetMacro(NumberOfMaterials, unsigned int);
  itkGetMacro(NumberOfMaterials, unsigned int);

  itkSetMacro(NumberOfEnergies, unsigned int);
  itkGetMacro(NumberOfEnergies, unsigned int);

  itkSetMacro(IsSpectralCT, bool);
  itkGetMacro(IsSpectralCT, bool);

  itkSetMacro(ComputeVariances, bool);
  itkGetMacro(ComputeVariances, bool);

  itkSetMacro(ComputeCramerRaoLowerBound, bool);
  itkGetMacro(ComputeCramerRaoLowerBound, bool);

protected:
  SpectralForwardModelImageFilter();
  ~SpectralForwardModelImageFilter() override = default;

  using DataObjectPointerArraySizeType = itk::ProcessObject::DataObjectPointerArraySizeType;
  using Superclass::MakeOutput;
  itk::DataObject::Pointer
  MakeOutput(DataObjectPointerArraySizeType idx) override;

  void
  GenerateOutputInformation() override;

  void
  GenerateInputRequestedRegion() override;

  void
  BeforeThreadedGenerateData() override;
  void
  DynamicThreadedGenerateData(const typename OutputImageType::RegionType & outputRegionForThread) override;

  void
  VerifyInputInformation() const override
  {}

  ThresholdsType m_Thresholds;
  unsigned int   m_NumberOfSpectralBins;

  MaterialAttenuationsType m_MaterialAttenuations;
  DetectorResponseType     m_DetectorResponse;
  unsigned int             m_NumberOfEnergies;

  unsigned int m_NumberOfIterations;
  unsigned int m_NumberOfMaterials;
  bool         m_OptimizeWithRestarts;
  bool         m_IsSpectralCT;               // If not, it is dual energy CT
  bool         m_ComputeVariances;           // Only implemented for dual energy CT
  bool         m_ComputeCramerRaoLowerBound; // Only implemented for spectral CT

#ifndef ITK_FUTURE_LEGACY_REMOVE

  typename FlattenVectorFilterType::Pointer m_FlattenFilter;
  typename FlattenVectorFilterType::Pointer m_FlattenSecondFilter;
  typename PermuteFilterType::Pointer       m_PermuteFilter;
  typename PermuteFilterType::Pointer       m_PermuteSecondFilter;
#endif

}; // end of class

// Function to bin a detector response matrix according to given energy thresholds
template <typename OutputElementType, typename DetectorResponseImageType, typename ThresholdsType>
vnl_matrix<OutputElementType>
SpectralBinDetectorResponse(const DetectorResponseImageType * drm,
                            const ThresholdsType &            thresholds,
                            const unsigned int                numberOfEnergies);

} // end namespace rtk


#ifndef ITK_MANUAL_INSTANTIATION
#  include "rtkSpectralForwardModelImageFilter.hxx"
#endif

#endif