RTK  2.0.1
Reconstruction Toolkit
rtkReconstructionConjugateGradientOperator.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 rtkReconstructionConjugateGradientOperator_h
20 #define rtkReconstructionConjugateGradientOperator_h
21 
22 #include <itkMultiplyImageFilter.h>
23 #include <itkAddImageFilter.h>
24 
25 #include "rtkConstantImageSource.h"
26 
30 
34 
35 #ifdef RTK_USE_CUDA
38 #endif
39 
40 namespace rtk
41 {
42 
117 template< typename TOutputImage,
118  typename TSingleComponentImage = TOutputImage,
119  typename TWeightsImage = TOutputImage>
121 {
122 public:
123  ITK_DISALLOW_COPY_AND_ASSIGN(ReconstructionConjugateGradientOperator);
124 
129 #ifdef RTK_USE_CUDA
130  using GradientImageType = itk::CudaImage<itk::CovariantVector<typename TOutputImage::PixelType, TOutputImage::ImageDimension>, TOutputImage::ImageDimension >;
131 #else
133 #endif
134 
136  itkNewMacro(Self)
137 
138 
139  void SetInputVolume(const TOutputImage* vol);
140  void SetInputProjectionStack(const TOutputImage* projs);
141  void SetInputWeights(const TWeightsImage* weights);
143 
145  itkTypeMacro(rtkReconstructionConjugateGradientOperator, ConjugateGradientOperator)
146 
147  using BackProjectionFilterType = rtk::BackProjectionImageFilter< TOutputImage, TOutputImage >;
149 
150  using ForwardProjectionFilterType = rtk::ForwardProjectionImageFilter< TOutputImage, TOutputImage >;
152 
154  using MultiplyFilterType = itk::MultiplyImageFilter<TOutputImage, TSingleComponentImage>;
155  using AddFilterType = itk::AddImageFilter<TOutputImage>;
156 
157  // If TOutputImage is an itk::Image of floats or double, so are the weights, and a simple Multiply filter is required
158  // If TOutputImage is an itk::Image of itk::Vector<float (or double)>, a BlockDiagonalMatrixVectorMultiply filter
159  // is needed. Thus the meta-programming construct
161  using PlainMultiplyFilterType = itk::MultiplyImageFilter<TOutputImage, TOutputImage, TOutputImage>;
162  typedef typename std::conditional<std::is_same< TSingleComponentImage, TOutputImage>::value,
165 
166  using OutputImagePointer = typename TOutputImage::Pointer;
167 
170 
173 
175  void SetSupportMask(const TSingleComponentImage *SupportMask);
176  typename TSingleComponentImage::ConstPointer GetSupportMask();
178 
180  itkSetConstObjectMacro(Geometry, ThreeDCircularProjectionGeometry)
181 
184  itkSetMacro(Gamma, float)
185  itkGetMacro(Gamma, float)
186  itkSetMacro(Tikhonov, float)
187  itkGetMacro(Tikhonov, float)
188 
189 protected:
191  ~ReconstructionConjugateGradientOperator() override = default;
192 
194  void GenerateData() override;
195 
196  template < typename ImageType >
197  typename std::enable_if< std::is_same< TSingleComponentImage, ImageType >::value, ImageType >::type::Pointer
199 
200  template < typename ImageType >
201  typename std::enable_if< !std::is_same< TSingleComponentImage, ImageType >::value, ImageType >::type::Pointer
202  ConnectGradientRegularization();
203 
207 
216  typename itk::ImageToImageFilter<TOutputImage, TOutputImage>::Pointer m_LaplacianFilter;
217  typename MultiplyWithWeightsFilterType::Pointer m_MultiplyWithWeightsFilter;
218 
221  float m_Gamma{0}; //Strength of the laplacian regularization
222  float m_Tikhonov{0}; //Strength of the Tikhonov regularization
223 
225  typename TOutputImage::Pointer m_FloatingInputPointer, m_FloatingOutputPointer;
226 
229 #if ITK_VERSION_MAJOR<5
230  void VerifyInputInformation() override {}
231 #else
232  void VerifyInputInformation() const override {}
233 #endif
234 
235 
237  void GenerateInputRequestedRegion() override;
238  void GenerateOutputInformation() override;
240 
242  typename TOutputImage::ConstPointer GetInputVolume();
243  typename TOutputImage::ConstPointer GetInputProjectionStack();
244  typename TWeightsImage::ConstPointer GetInputWeights();
246 
247 };
248 } //namespace RTK
249 
250 
251 #ifndef ITK_MANUAL_INSTANTIATION
252 #include "rtkReconstructionConjugateGradientOperator.hxx"
253 #endif
254 
255 #endif
TOutputImage::ConstPointer GetInputVolume()
typename OutputImageType::Pointer OutputImagePointer
Implements the operator A used in conjugate gradient reconstruction.
void SetBackProjectionFilter(const BackProjectionFilterPointer _arg)
rtk::ThreeDCircularProjectionGeometry::ConstPointer m_Geometry
Generate an n-dimensional image with constant pixel values.
TWeightsImage::ConstPointer GetInputWeights()
void SetForwardProjectionFilter(const ForwardProjectionFilterPointer _arg)
STL namespace.
TSingleComponentImage::ConstPointer GetSupportMask()
void SetInputWeights(const TWeightsImage *weights)
TOutputImage::ConstPointer GetInputProjectionStack()
typename BackProjectionFilterType::Pointer BackProjectionFilterPointer
Projection geometry for a source and a 2-D flat panel.
typename ForwardProjectionFilterType::Pointer ForwardProjectionFilterPointer
void SetSupportMask(const TSingleComponentImage *SupportMask)
std::enable_if< std::is_same< TSingleComponentImage, ImageType >::value, ImageType >::type::Pointer ConnectGradientRegularization()
itk::ImageToImageFilter< TOutputImage, TOutputImage >::Pointer m_LaplacianFilter
void SetInputProjectionStack(const TOutputImage *projs)
void SetInputVolume(const TOutputImage *vol)
std::conditional< std::is_same< TSingleComponentImage, TOutputImage >::value, PlainMultiplyFilterType, MatrixVectorMultiplyFilterType >::type MultiplyWithWeightsFilterType
#define itkSetMacro(name, type)