RTK  2.5.0
Reconstruction Toolkit
rtkIterativeConeBeamReconstructionFilter.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 rtkIterativeConeBeamReconstructionFilter_h
20 #define rtkIterativeConeBeamReconstructionFilter_h
21 
22 #include <itkPixelTraits.h>
23 
24 // Forward projection filters
25 #include "rtkConfiguration.h"
28 // Back projection filters
31 
32 #ifdef RTK_USE_CUDA
36 #endif
37 
38 #include <random>
39 #include <algorithm>
40 
41 namespace rtk
42 {
43 
56 template <class TOutputImage, class ProjectionStackType = TOutputImage>
57 class ITK_TEMPLATE_EXPORT IterativeConeBeamReconstructionFilter
58  : public itk::ImageToImageFilter<TOutputImage, TOutputImage>
59 {
60 public:
61  ITK_DISALLOW_COPY_AND_MOVE(IterativeConeBeamReconstructionFilter);
62 
68 
70  using VolumeType = ProjectionStackType;
72  typedef enum
73  {
74  FP_JOSEPH = 0,
75  FP_CUDARAYCAST = 2,
76  FP_JOSEPHATTENUATED = 3,
77  FP_ZENG = 4
78  } ForwardProjectionType;
79  typedef enum
80  {
81  BP_VOXELBASED = 0,
82  BP_JOSEPH = 1,
83  BP_CUDAVOXELBASED = 2,
84  BP_CUDARAYCAST = 4,
85  BP_JOSEPHATTENUATED = 5,
86  BP_ZENG = 6
87  } BackProjectionType;
88 
94 
96  itkNewMacro(Self);
97 
99 #ifdef itkOverrideGetNameOfClassMacro
100  itkOverrideGetNameOfClassMacro(IterativeConeBeamReconstructionFilter);
101 #else
103 #endif
104 
105 
107  virtual void
108  SetForwardProjectionFilter(ForwardProjectionType fwtype);
111  {
112  return m_CurrentForwardProjectionConfiguration;
113  }
114  virtual void
115  SetBackProjectionFilter(BackProjectionType bptype);
116  BackProjectionType
118  {
119  return m_CurrentBackProjectionConfiguration;
120  }
122 
125  void
126  SetAttenuationMap(const VolumeType * attenuationMap)
127  {
128  // Process object is not const-correct so the const casting is required.
129  this->SetNthInput(2, const_cast<VolumeType *>(attenuationMap));
130  }
131  typename VolumeType::ConstPointer
133  {
134  return static_cast<const VolumeType *>(this->itk::ProcessObject::GetInput(2));
135  }
137 
141  void
142  SetInferiorClipImage(const TClipImageType * inferiorClipImage)
143  {
144  // Process object is not const-correct so the const casting is required.
145  this->SetInput("InferiorClipImage", const_cast<TClipImageType *>(inferiorClipImage));
146  }
147  typename TClipImageType::ConstPointer
149  {
150  return static_cast<const TClipImageType *>(this->itk::ProcessObject::GetInput("InferiorClipImage"));
151  }
153 
157  void
158  SetSuperiorClipImage(const TClipImageType * superiorClipImage)
159  {
160  // Process object is not const-correct so the const casting is required.
161  this->SetInput("SuperiorClipImage", const_cast<TClipImageType *>(superiorClipImage));
162  }
163  typename TClipImageType::ConstPointer
165  {
166  return static_cast<const TClipImageType *>(this->itk::ProcessObject::GetInput("SuperiorClipImage"));
167  }
169 
171  itkGetMacro(SigmaZero, double);
172  itkSetMacro(SigmaZero, double);
174 
176  itkGetMacro(AlphaPSF, double);
177  itkSetMacro(AlphaPSF, double);
179 
181  itkGetConstMacro(StepSize, double);
182  itkSetMacro(StepSize, double);
184 
185 protected:
187  ~IterativeConeBeamReconstructionFilter() override = default;
188 
191  virtual BackProjectionPointerType
192  InstantiateBackProjectionFilter(int bptype);
193 
196  virtual ForwardProjectionPointerType
197  InstantiateForwardProjectionFilter(int fwtype);
198 
203 
206  std::default_random_engine m_DefaultRandomEngine = std::default_random_engine{};
207 
209  double m_SigmaZero{ 1.5417233052142099 };
210  double m_AlphaPSF{ 0.016241189545787734 };
211 
213  double m_StepSize{ 1.0 };
214 
216  using CPUImageType =
218  template <typename ImageType>
219  using EnableCudaScalarAndVectorType = typename std::enable_if<
220  !std::is_same<CPUImageType, ImageType>::value &&
221  std::is_same<typename itk::PixelTraits<typename ImageType::PixelType>::ValueType, float>::value &&
225  template <typename ImageType>
226  using DisableCudaScalarAndVectorType = typename std::enable_if<
227  std::is_same<CPUImageType, ImageType>::value ||
228  !std::is_same<typename itk::PixelTraits<typename ImageType::PixelType>::ValueType, float>::value ||
232  template <typename ImageType>
233  using EnableCudaScalarType = typename std::enable_if<
234  !std::is_same<CPUImageType, ImageType>::value &&
235  std::is_same<typename itk::PixelTraits<typename ImageType::PixelType>::ValueType, float>::value &&
237  template <typename ImageType>
238  using DisableCudaScalarType = typename std::enable_if<
239  std::is_same<CPUImageType, ImageType>::value ||
240  !std::is_same<typename itk::PixelTraits<typename ImageType::PixelType>::ValueType, float>::value ||
242  template <typename ImageType>
243  using EnableVectorType =
244  typename std::enable_if<itk::PixelTraits<typename ImageType::PixelType>::Dimension != 1>::type;
245  template <typename ImageType>
246  using DisableVectorType =
247  typename std::enable_if<itk::PixelTraits<typename ImageType::PixelType>::Dimension == 1>::type;
249 
250  template <typename ImageType, EnableCudaScalarAndVectorType<ImageType> * = nullptr>
253  {
255 #ifdef RTK_USE_CUDA
256  fw = CudaForwardProjectionImageFilter<ImageType, ImageType>::New();
257  dynamic_cast<rtk::CudaForwardProjectionImageFilter<ImageType, ImageType> *>(fw.GetPointer())
258  ->SetStepSize(m_StepSize);
259 #endif
260  return fw;
261  }
262 
263 
264  template <typename ImageType, DisableCudaScalarAndVectorType<ImageType> * = nullptr>
265  ForwardProjectionPointerType
267  {
268  itkGenericExceptionMacro(
269  << "CudaRayCastBackProjectionImageFilter only available with 3D CudaImage of float or itk::Vector<float,3>.");
270  return nullptr;
271  }
272 
273 
274  template <typename ImageType, EnableVectorType<ImageType> * = nullptr>
275  ForwardProjectionPointerType
277  {
278  itkGenericExceptionMacro(<< "JosephForwardAttenuatedProjectionImageFilter only available with scalar pixel types.");
279  return nullptr;
280  }
281 
282 
283  template <typename ImageType, DisableVectorType<ImageType> * = nullptr>
284  ForwardProjectionPointerType
286  {
289  if (this->GetAttenuationMap().IsNotNull())
290  {
291  fw->SetInput(2, this->GetAttenuationMap());
292  }
293  else
294  {
295  itkExceptionMacro(<< "Set Joseph attenuated forward projection filter but no attenuation map is given");
296  return nullptr;
297  }
298  if (this->GetSuperiorClipImage().IsNotNull())
299  {
301  fw.GetPointer())
302  ->SetSuperiorClipImage(this->GetSuperiorClipImage());
303  }
304  if (this->GetInferiorClipImage().IsNotNull())
305  {
307  fw.GetPointer())
308  ->SetInferiorClipImage(this->GetInferiorClipImage());
309  }
310  return fw;
311  }
312 
313  template <typename ImageType, EnableVectorType<ImageType> * = nullptr>
314  ForwardProjectionPointerType
316  {
317  itkGenericExceptionMacro(<< "JosephForwardAttenuatedProjectionImageFilter only available with scalar pixel types.");
318  return nullptr;
319  }
320 
321 
322  template <typename ImageType, DisableVectorType<ImageType> * = nullptr>
323  ForwardProjectionPointerType
325  {
328  if (this->GetAttenuationMap().IsNotNull())
329  {
330  fw->SetInput(2, this->GetAttenuationMap());
331  }
333  ->SetSigmaZero(m_SigmaZero);
335  ->SetAlpha(m_AlphaPSF);
336  return fw;
337  }
338 
339  template <typename ImageType, EnableCudaScalarAndVectorType<ImageType> * = nullptr>
340  BackProjectionPointerType
342  {
344 #ifdef RTK_USE_CUDA
345  bp = CudaBackProjectionImageFilter<ImageType>::New();
346 #endif
347  return bp;
348  }
349 
350 
351  template <typename ImageType, DisableCudaScalarAndVectorType<ImageType> * = nullptr>
352  BackProjectionPointerType
354  {
355  itkGenericExceptionMacro(
356  << "CudaBackProjectionImageFilter only available with 3D CudaImage of float or itk::Vector<float,3>.");
357  return nullptr;
358  }
359 
360 
361  template <typename ImageType, EnableCudaScalarType<ImageType> * = nullptr>
362  BackProjectionPointerType
364  {
366 #ifdef RTK_USE_CUDA
367  bp = CudaRayCastBackProjectionImageFilter::New();
368  dynamic_cast<rtk::CudaRayCastBackProjectionImageFilter *>(bp.GetPointer())->SetStepSize(m_StepSize);
369 #endif
370  return bp;
371  }
372 
373 
374  template <typename ImageType, DisableCudaScalarType<ImageType> * = nullptr>
375  BackProjectionPointerType
377  {
378  itkGenericExceptionMacro(<< "CudaRayCastBackProjectionImageFilter only available with 3D CudaImage of float.");
379  return nullptr;
380  }
381 
382 
383  template <typename ImageType, EnableVectorType<ImageType> * = nullptr>
384  BackProjectionPointerType
386  {
387  itkGenericExceptionMacro(<< "JosephBackAttenuatedProjectionImageFilter only available with scalar pixel types.");
388  return nullptr;
389  }
390 
391 
392  template <typename ImageType, DisableVectorType<ImageType> * = nullptr>
393  BackProjectionPointerType
395  {
398  if (this->GetAttenuationMap().IsNotNull())
399  {
400  bp->SetInput(2, this->GetAttenuationMap());
401  return bp;
402  }
403  else
404  {
405  itkExceptionMacro(<< "Set Joseph attenuated backprojection filter but no attenuation map is given");
406  return nullptr;
407  }
408  }
409 
410  template <typename ImageType, EnableVectorType<ImageType> * = nullptr>
411  BackProjectionPointerType
413  {
414  itkGenericExceptionMacro(<< "JosephBackAttenuatedProjectionImageFilter only available with scalar pixel types.");
415  return nullptr;
416  }
417 
418 
419  template <typename ImageType, DisableVectorType<ImageType> * = nullptr>
420  BackProjectionPointerType
422  {
425  if (this->GetAttenuationMap().IsNotNull())
426  {
427  bp->SetInput(2, this->GetAttenuationMap());
428  }
430  ->SetSigmaZero(m_SigmaZero);
432  ->SetAlpha(m_AlphaPSF);
433  return bp;
434  }
435 
436 }; // end of class
437 
438 } // end namespace rtk
439 
440 #ifndef ITK_MANUAL_INSTANTIATION
441 # include "rtkIterativeConeBeamReconstructionFilter.hxx"
442 #endif
443 
444 #endif
Base class for forward projection, i.e. accumulation along x-ray lines.
void SetSuperiorClipImage(const TClipImageType *superiorClipImage)
typename std::enable_if< itk::PixelTraits< typename ImageType::PixelType >::Dimension==1 >::type DisableVectorType
typename std::enable_if< !std::is_same< CPUImageType, ImageType >::value &&std::is_same< typename itk::PixelTraits< typename ImageType::PixelType >::ValueType, float >::value &&(itk::PixelTraits< typename ImageType::PixelType >::Dimension==1||itk::PixelTraits< typename ImageType::PixelType >::Dimension==2||itk::PixelTraits< typename ImageType::PixelType >::Dimension==3)>::type EnableCudaScalarAndVectorType
typename std::enable_if< !std::is_same< CPUImageType, ImageType >::value &&std::is_same< typename itk::PixelTraits< typename ImageType::PixelType >::ValueType, float >::value &&itk::PixelTraits< typename ImageType::PixelType >::Dimension==1 >::type EnableCudaScalarType
#define itkSetMacro(name, type)
DataObject * GetInput(const DataObjectIdentifierType &key)
void SetInferiorClipImage(const TClipImageType *inferiorClipImage)
Mother class for cone beam reconstruction filters which need runtime selection of their forward and b...
typename std::enable_if< std::is_same< CPUImageType, ImageType >::value||!std::is_same< typename itk::PixelTraits< typename ImageType::PixelType >::ValueType, float >::value||(itk::PixelTraits< typename ImageType::PixelType >::Dimension !=1 &&itk::PixelTraits< typename ImageType::PixelType >::Dimension !=2 &&itk::PixelTraits< typename ImageType::PixelType >::Dimension !=3)>::type DisableCudaScalarAndVectorType
typename std::enable_if< itk::PixelTraits< typename ImageType::PixelType >::Dimension !=1 >::type EnableVectorType
typename std::enable_if< std::is_same< CPUImageType, ImageType >::value||!std::is_same< typename itk::PixelTraits< typename ImageType::PixelType >::ValueType, float >::value||itk::PixelTraits< typename ImageType::PixelType >::Dimension !=1 >::type DisableCudaScalarType