rtkSimplexSpectralProjectionsDecompositionImageFilter.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 rtkSimplexSpectralProjectionsDecompositionImageFilter_h
#define rtkSimplexSpectralProjectionsDecompositionImageFilter_h

#include <itkImageToImageFilter.h>
#include <itkAmoebaOptimizer.h>
#include "rtkSchlomka2008NegativeLogLikelihood.h"
#include "rtkDualEnergyNegativeLogLikelihood.h"
#include "rtkVectorImageToImageFilter.h"

#include <itkPermuteAxesImageFilter.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 SimplexSpectralProjectionsDecompositionImageFilter
  : public itk::ImageToImageFilter<DecomposedProjectionsType, DecomposedProjectionsType>
{
public:
  ITK_DISALLOW_COPY_AND_MOVE(SimplexSpectralProjectionsDecompositionImageFilter);

  using Self = SimplexSpectralProjectionsDecompositionImageFilter;
  using Superclass = itk::ImageToImageFilter<DecomposedProjectionsType, DecomposedProjectionsType>;
  using Pointer = itk::SmartPointer<Self>;
  using ConstPointer = itk::SmartPointer<const Self>;

  using InputImageType = DecomposedProjectionsType;
  using OutputImageType = DecomposedProjectionsType;

  using ThresholdsType = itk::VariableLengthVector<double>;
  using MeanAttenuationInBinType = itk::VariableSizeMatrix<double>;
  using DetectorResponseType = vnl_matrix<double>;
  using MaterialAttenuationsType = vnl_matrix<double>;
  using CostFunctionType = ProjectionsDecompositionNegativeLogLikelihood;
  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(SimplexSpectralProjectionsDecompositionImageFilter);

  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();

  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();

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

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

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

  itkSetMacro(OptimizeWithRestarts, bool);
  itkGetMacro(OptimizeWithRestarts, bool);

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

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

  itkSetMacro(OutputInverseCramerRaoLowerBound, bool);
  itkGetMacro(OutputInverseCramerRaoLowerBound, bool);

  itkSetMacro(OutputFischerMatrix, bool);
  itkGetMacro(OutputFischerMatrix, bool);

  itkSetMacro(LogTransformEachBin, bool);
  itkGetMacro(LogTransformEachBin, bool);

  itkSetMacro(GuessInitialization, bool);
  itkGetMacro(GuessInitialization, bool);

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

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

  void
  GenerateOutputInformation() override;

  void
  GenerateInputRequestedRegion() override;

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

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

  void
  VerifyInputInformation() const override
  {}

  MaterialAttenuationsType m_MaterialAttenuations;
  DetectorResponseType     m_DetectorResponse;
  ThresholdsType           m_Thresholds;
  MeanAttenuationInBinType m_MeanAttenuationInBin;
  bool                     m_OutputInverseCramerRaoLowerBound;
  bool                     m_OutputFischerMatrix;
  bool                     m_LogTransformEachBin;
  bool                     m_GuessInitialization;
  bool                     m_IsSpectralCT; // If not, it is dual energy CT
  bool                     m_OptimizeWithRestarts;
  unsigned int             m_NumberOfIterations;
  unsigned int             m_NumberOfMaterials;
  unsigned int             m_NumberOfEnergies;
  unsigned int             m_NumberOfSpectralBins;

#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

} // end namespace rtk


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

#endif