RTK  2.4.1
Reconstruction Toolkit
rtkMechlemOneStepSpectralReconstructionFilter.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright RTK Consortium
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * https://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
18 
19 #ifndef rtkMechlemOneStepSpectralReconstructionFilter_h
20 #define rtkMechlemOneStepSpectralReconstructionFilter_h
21 
26 #include "rtkConstantImageSource.h"
30 
31 #include <itkExtractImageFilter.h>
32 #include <itkAddImageFilter.h>
33 #include <itkMultiplyImageFilter.h>
34 
35 #ifdef RTK_USE_CUDA
37 #endif
38 
39 namespace rtk
40 {
145 template <typename TOutputImage, typename TPhotonCounts, typename TSpectrum>
147  : public rtk::IterativeConeBeamReconstructionFilter<TOutputImage, TOutputImage>
148 {
149 public:
150  ITK_DISALLOW_COPY_AND_MOVE(MechlemOneStepSpectralReconstructionFilter);
151 
156 
158  itkNewMacro(Self);
159 
162 
164  static constexpr unsigned int nBins = TPhotonCounts::PixelType::Dimension;
165  static constexpr unsigned int nMaterials = TOutputImage::PixelType::Dimension;
166  using dataType = typename TOutputImage::PixelType::ValueType;
167 
170 #ifdef RTK_USE_CUDA
171  typedef
172  typename std::conditional<std::is_same<TOutputImage, CPUOutputImageType>::value,
173  itk::Image<itk::Vector<dataType, nMaterials * nMaterials>, TOutputImage::ImageDimension>,
174  itk::CudaImage<itk::Vector<dataType, nMaterials * nMaterials>,
175  TOutputImage::ImageDimension>>::type HessiansImageType;
176  typedef
177  typename std::conditional<std::is_same<TOutputImage, CPUOutputImageType>::value,
179  itk::CudaImage<dataType, TOutputImage::ImageDimension>>::type SingleComponentImageType;
180 #else
181  using HessiansImageType =
182  typename itk::Image<itk::Vector<dataType, nMaterials * nMaterials>, TOutputImage::ImageDimension>;
184 #endif
185 
186 #if !defined(ITK_WRAPPING_PARSER)
187 # ifdef RTK_USE_CUDA
188  typedef typename std::conditional<std::is_same<TOutputImage, CPUOutputImageType>::value,
190  CudaWeidingerForwardModelImageFilter<TOutputImage, TPhotonCounts, TSpectrum>>::type
192  typedef
193  typename std::conditional<std::is_same<TOutputImage, CPUOutputImageType>::value,
195  CudaForwardProjectionImageFilter<SingleComponentImageType, SingleComponentImageType>>::
197  typedef typename std::conditional<std::is_same<TOutputImage, CPUOutputImageType>::value,
199  CudaBackProjectionImageFilter<HessiansImageType>>::type
201 # else
206 # endif
207  using GradientsImageType = TOutputImage;
208 #endif
209 
210  using ForwardProjectionType = typename Superclass::ForwardProjectionType;
211  using BackProjectionType = typename Superclass::BackProjectionType;
212 
213 #if !defined(ITK_WRAPPING_PARSER)
214 
216  using AddFilterType = itk::AddImageFilter<GradientsImageType>;
232 #endif
233 
235  itkSetConstObjectMacro(Geometry, ThreeDCircularProjectionGeometry);
236 
237  itkSetMacro(NumberOfIterations, int);
238  itkGetMacro(NumberOfIterations, int);
239 
241  itkSetMacro(NumberOfSubsets, int);
242  itkGetMacro(NumberOfSubsets, int);
244 
248  itkSetMacro(ResetNesterovEvery, int);
249  itkGetMacro(ResetNesterovEvery, int);
251 
253  void
254  SetInputMaterialVolumes(const TOutputImage * materialVolumes);
255  void
256  SetInputPhotonCounts(const TPhotonCounts * photonCounts);
257  void
258  SetInputSpectrum(const TSpectrum * spectrum);
259  void
260  SetSupportMask(const SingleComponentImageType * support);
261  void
262  SetSpatialRegularizationWeights(const SingleComponentImageType * regweights);
263  void
264  SetProjectionWeights(const SingleComponentImageType * weiprojections);
266 
268  itkSetMacro(RegularizationWeights, typename TOutputImage::PixelType);
269  itkGetMacro(RegularizationWeights, typename TOutputImage::PixelType);
271 
273  itkSetMacro(RegularizationRadius, typename TOutputImage::RegionType::SizeType);
274  itkGetMacro(RegularizationRadius, typename TOutputImage::RegionType::SizeType);
276 
279  using BinnedDetectorResponseType = vnl_matrix<dataType>;
280  using MaterialAttenuationsType = vnl_matrix<dataType>;
281  virtual void
282  SetBinnedDetectorResponse(const BinnedDetectorResponseType & detResp);
283  virtual void
284  SetMaterialAttenuations(const MaterialAttenuationsType & matAtt);
286 
287 protected:
289  ~MechlemOneStepSpectralReconstructionFilter() override = default;
290 
292  void
293  VerifyPreconditions() ITKv5_CONST override;
294 
296  void
297  GenerateData() override;
298 
299 #if !defined(ITK_WRAPPING_PARSER)
300 
302  typename AddFilterType::Pointer m_AddGradients;
321 #endif
322 
326  void
327  VerifyInputInformation() const override
328  {}
329 
332  void
333  GenerateInputRequestedRegion() override;
334  void
335  GenerateOutputInformation() override;
337 
339  typename TOutputImage::ConstPointer
340  GetInputMaterialVolumes();
341  typename TPhotonCounts::ConstPointer
342  GetInputPhotonCounts();
343  typename TSpectrum::ConstPointer
344  GetInputSpectrum();
345  typename SingleComponentImageType::ConstPointer
346  GetSupportMask();
347  typename SingleComponentImageType::ConstPointer
348  GetSpatialRegularizationWeights();
349  typename SingleComponentImageType::ConstPointer
350  GetProjectionWeights();
352 
353 #if !defined(ITK_WRAPPING_PARSER)
354 
356  typename SingleComponentForwardProjectionFilterType::Pointer
357  InstantiateSingleComponentForwardProjectionFilter(int fwtype);
358  typename HessiansBackProjectionFilterType::Pointer
359  InstantiateHessiansBackProjectionFilter(int bptype);
360 #endif
361 
362 
364 
371 
372  typename TOutputImage::PixelType m_RegularizationWeights;
373  typename TOutputImage::RegionType::SizeType m_RegularizationRadius;
374 };
375 } // namespace rtk
376 
377 
378 #ifndef ITK_MANUAL_INSTANTIATION
379 # include "rtkMechlemOneStepSpectralReconstructionFilter.hxx"
380 #endif
381 
382 #endif
typename itk::Image< dataType, TOutputImage::ImageDimension > SingleComponentImageType
Base class for forward projection, i.e. accumulation along x-ray lines.
Applies Nesterov&#39;s momentum technique.
Generate an n-dimensional image with constant pixel values.
typename itk::Image< itk::Vector< dataType, nMaterials *nMaterials >, TOutputImage::ImageDimension > HessiansImageType
Projection geometry for a source and a 2-D flat panel.
#define itkSetMacro(name, type)
Computes update from gradient and Hessian in Newton&#39;s method.
Performs intermediate computations in Weidinger2016.
typename itk::Image< typename TOutputImage::PixelType, TOutputImage::ImageDimension > CPUOutputImageType
For each vector-valued pixel, adds a vector to the diagonal of a matrix.
For one-step inversion of spectral CT data by the method Mechlem2017, computes regularization term&#39;s ...
Mother class for cone beam reconstruction filters which need runtime selection of their forward and b...
Implements the one-step spectral CT inversion method described by Mechlem et al.
SingleComponentForwardProjectionFilterType::Pointer m_SingleComponentForwardProjectionFilter