RTK  2.0.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  * http://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 {
125 template< typename TOutputImage, typename TPhotonCounts, typename TSpectrum >
127 {
128 public:
129  ITK_DISALLOW_COPY_AND_ASSIGN(MechlemOneStepSpectralReconstructionFilter);
130 
135 
137  itkNewMacro(Self)
138 
139 
141 
143  static constexpr unsigned int nBins = TPhotonCounts::PixelType::Dimension;
144  static constexpr unsigned int nMaterials = TOutputImage::PixelType::Dimension;
145  using dataType = typename TOutputImage::PixelType::ValueType;
146 
148  using CPUOutputImageType = typename itk::Image< typename TOutputImage::PixelType,
149  TOutputImage::ImageDimension>;
150 #ifdef RTK_USE_CUDA
151  typedef typename std::conditional< std::is_same< TOutputImage, CPUOutputImageType >::value,
152  itk::Image< itk::Vector<dataType, nMaterials * nMaterials>, TOutputImage::ImageDimension >,
153  itk::CudaImage< itk::Vector<dataType, nMaterials * nMaterials>, TOutputImage::ImageDimension > >::type
155  typedef typename std::conditional< std::is_same< TOutputImage, CPUOutputImageType >::value,
157  itk::CudaImage<dataType, TOutputImage::ImageDimension> >::type
159 #else
160  using THessiansImage = typename itk::Image< itk::Vector<dataType, nMaterials * nMaterials>,
161  TOutputImage::ImageDimension >;
162  using SingleComponentImageType = typename itk::Image<dataType, TOutputImage::ImageDimension>;
163 #endif
164 
165 #if !defined( ITK_WRAPPING_PARSER )
166 #ifdef RTK_USE_CUDA
167  typedef typename std::conditional< std::is_same< TOutputImage, CPUOutputImageType >::value,
169  CudaWeidingerForwardModelImageFilter<TOutputImage, TPhotonCounts, TSpectrum> >::type
171  typedef typename std::conditional< std::is_same< TOutputImage, CPUOutputImageType >::value,
173  CudaForwardProjectionImageFilter<SingleComponentImageType, SingleComponentImageType> >::type
175  typedef typename std::conditional< std::is_same< TOutputImage, CPUOutputImageType >::value,
177  CudaBackProjectionImageFilter<THessiansImage> >::type
179 #else
180  using WeidingerForwardModelType = WeidingerForwardModelImageFilter<TOutputImage, TPhotonCounts, TSpectrum>;
181  using CudaSingleComponentForwardProjectionImageFilterType = JosephForwardProjectionImageFilter<SingleComponentImageType,
182  SingleComponentImageType>;
183  using CudaHessiansBackProjectionImageFilterType = BackProjectionImageFilter<THessiansImage, THessiansImage>;
184 #endif
185  using TGradientsImage = TOutputImage;
186 #endif
187 
190 
191 #if !defined( ITK_WRAPPING_PARSER )
192 
196  SingleComponentImageType >;
209 #endif
210 
213 
215  void SetBackProjectionFilter (BackProjectionType _arg) override;
216 
218  itkSetConstObjectMacro(Geometry, ThreeDCircularProjectionGeometry)
219 
220  itkSetMacro(NumberOfIterations, int)
221  itkGetMacro(NumberOfIterations, int)
222 
224  itkSetMacro(NumberOfSubsets, int)
225  itkGetMacro(NumberOfSubsets, int)
226 
230  itkSetMacro(ResetNesterovEvery, int)
231  itkGetMacro(ResetNesterovEvery, int)
232 
234  void SetInputMaterialVolumes(const TOutputImage* materialVolumes);
235  void SetInputPhotonCounts(const TPhotonCounts* photonCounts);
236  void SetInputSpectrum(const TSpectrum* spectrum);
237  void SetSupportMask(const SingleComponentImageType* support);
239 
241  itkSetMacro(RegularizationWeights, typename TOutputImage::PixelType)
242  itkGetMacro(RegularizationWeights, typename TOutputImage::PixelType)
243 
245  itkSetMacro(RegularizationRadius, typename TOutputImage::RegionType::SizeType)
246  itkGetMacro(RegularizationRadius, typename TOutputImage::RegionType::SizeType)
247 
248 // itkSetMacro(IterationCosts, bool)
249 // itkGetMacro(IterationCosts, bool)
250 
254  using MaterialAttenuationsType = vnl_matrix<dataType>;
255  virtual void SetBinnedDetectorResponse(const BinnedDetectorResponseType & detResp);
256  virtual void SetMaterialAttenuations(const MaterialAttenuationsType & matAtt);
258 
259 protected:
260  MechlemOneStepSpectralReconstructionFilter();
261  ~MechlemOneStepSpectralReconstructionFilter() override = default;
262 
264  void GenerateData() override;
265 
266 #if !defined( ITK_WRAPPING_PARSER )
267 
268  typename ExtractPhotonCountsFilterType::Pointer m_ExtractPhotonCountsFilter;
285 #endif
286 
290 #if ITK_VERSION_MAJOR<5
291  void VerifyInputInformation() override {}
292 #else
293  void VerifyInputInformation() const override {}
294 #endif
295 
296 
299  void GenerateInputRequestedRegion() override;
300  void GenerateOutputInformation() override;
302 
304  typename TOutputImage::ConstPointer GetInputMaterialVolumes();
305  typename TPhotonCounts::ConstPointer GetInputPhotonCounts();
306  typename TSpectrum::ConstPointer GetInputSpectrum();
307  typename SingleComponentImageType::ConstPointer GetSupportMask();
309 
310 #if !defined( ITK_WRAPPING_PARSER )
311 
315 #endif
316 
317 
319 
326 
327  typename TOutputImage::PixelType m_RegularizationWeights;
328  typename TOutputImage::RegionType::SizeType m_RegularizationRadius;
329 
330 // private:
331 // bool m_IterationCosts;
332 };
333 } //namespace ITK
334 
335 
336 #ifndef ITK_MANUAL_INSTANTIATION
337 #include "rtkMechlemOneStepSpectralReconstructionFilter.hxx"
338 #endif
339 
340 #endif
typename itk::Image< dataType, TOutputImage::ImageDimension > SingleComponentImageType
void SetInputPhotonCounts(const TPhotonCounts *photonCounts)
Base class for forward projection, i.e. accumulation along x-ray lines.
SingleComponentForwardProjectionFilterType::Pointer InstantiateSingleComponentForwardProjectionFilter(int fwtype)
HessiansBackProjectionFilterType::Pointer InstantiateHessiansBackProjectionFilter(int bptype)
Applies Nesterov&#39;s momentum technique.
JosephForwardProjectionImageFilter< SingleComponentImageType, SingleComponentImageType > CudaSingleComponentForwardProjectionImageFilterType
Generate an n-dimensional image with constant pixel values.
void SetSupportMask(const SingleComponentImageType *support)
typename itk::Image< typename TOutputImage::PixelType, TOutputImage::ImageDimension > CPUOutputImageType
void SetBackProjectionFilter(BackProjectionType _arg) override
BackProjectionImageFilter< THessiansImage, THessiansImage > CudaHessiansBackProjectionImageFilterType
virtual void SetMaterialAttenuations(const MaterialAttenuationsType &matAtt)
Projection geometry for a source and a 2-D flat panel.
ImageBaseType::SizeType SizeType
Computes update from gradient and Hessian in Newton&#39;s method.
void SetInputMaterialVolumes(const TOutputImage *materialVolumes)
typename itk::Image< itk::Vector< dataType, nMaterials *nMaterials >, TOutputImage::ImageDimension > THessiansImage
void SetInputSpectrum(const TSpectrum *spectrum)
Performs intermediate computations in Weidinger2016.
WeidingerForwardModelImageFilter< TOutputImage, TPhotonCounts, TSpectrum > WeidingerForwardModelType
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...
virtual void SetBinnedDetectorResponse(const BinnedDetectorResponseType &detResp)
SingleComponentImageType::ConstPointer GetSupportMask()
Implements the one-step spectral CT inversion method described by Mechlem et al.
ImageBaseType::RegionType RegionType
SingleComponentForwardProjectionFilterType::Pointer m_SingleComponentForwardProjectionFilter
TPhotonCounts::ConstPointer GetInputPhotonCounts()
TOutputImage::ConstPointer GetInputMaterialVolumes()
void SetForwardProjectionFilter(ForwardProjectionType _arg) override
#define itkSetMacro(name, type)