RTK  2.4.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  * 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 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, typename TSingleComponentImage = TOutputImage, typename TWeightsImage = TOutputImage>
118 class ITK_TEMPLATE_EXPORT ReconstructionConjugateGradientOperator : public ConjugateGradientOperator<TOutputImage>
119 {
120 public:
121  ITK_DISALLOW_COPY_AND_MOVE(ReconstructionConjugateGradientOperator);
122 
127 #ifdef RTK_USE_CUDA
128  using GradientImageType =
129  itk::CudaImage<itk::CovariantVector<typename TOutputImage::PixelType, TOutputImage::ImageDimension>,
130  TOutputImage::ImageDimension>;
131 #else
132  using GradientImageType =
134  TOutputImage::ImageDimension>;
135 #endif
136 
138  itkNewMacro(Self);
139 
141  void
142  SetInputVolume(const TOutputImage * vol);
143  void
144  SetInputProjectionStack(const TOutputImage * projs);
145  void
146  SetInputWeights(const TWeightsImage * weights);
148 
150  itkTypeMacro(rtkReconstructionConjugateGradientOperator, ConjugateGradientOperator);
151 
154 
157 
160  using AddFilterType = itk::AddImageFilter<TOutputImage>;
161 
162  // If TOutputImage is an itk::Image of floats or double, so are the weights, and a simple Multiply filter is required
163  // If TOutputImage is an itk::Image of itk::Vector<float (or double)>, a BlockDiagonalMatrixVectorMultiply filter
164  // is needed. Thus the meta-programming construct
167  typedef typename std::conditional<std::is_same<TSingleComponentImage, TOutputImage>::value,
170 
171  using OutputImagePointer = typename TOutputImage::Pointer;
172 
174  void
175  SetBackProjectionFilter(const BackProjectionFilterPointer _arg);
176 
178  void
179  SetForwardProjectionFilter(const ForwardProjectionFilterPointer _arg);
180 
182  void
183  SetSupportMask(const TSingleComponentImage * SupportMask);
184  typename TSingleComponentImage::ConstPointer
185  GetSupportMask();
187 
189  itkSetConstObjectMacro(Geometry, ThreeDCircularProjectionGeometry);
190 
193  itkSetMacro(Gamma, float);
194  itkGetMacro(Gamma, float);
195  itkSetMacro(Tikhonov, float);
196  itkGetMacro(Tikhonov, float);
198 
199 protected:
201  ~ReconstructionConjugateGradientOperator() override = default;
202 
204  void
205  VerifyPreconditions() ITKv5_CONST override;
206 
208  void
209  GenerateData() override;
210 
211  template <typename ImageType>
212  typename std::enable_if<std::is_same<TSingleComponentImage, ImageType>::value, ImageType>::type::Pointer
213  ConnectGradientRegularization();
214 
215  template <typename ImageType>
216  typename std::enable_if<!std::is_same<TSingleComponentImage, ImageType>::value, ImageType>::type::Pointer
217  ConnectGradientRegularization();
218 
220  BackProjectionFilterPointer m_BackProjectionFilter;
221  ForwardProjectionFilterPointer m_ForwardProjectionFilter;
222 
223  typename ConstantSourceType::Pointer m_ConstantProjectionsSource;
224  typename ConstantSourceType::Pointer m_ConstantVolumeSource;
225  typename MultiplyFilterType::Pointer m_MultiplyOutputVolumeFilter;
226  typename MultiplyFilterType::Pointer m_MultiplyInputVolumeFilter;
227  typename MultiplyFilterType::Pointer m_MultiplyLaplacianFilter;
228  typename MultiplyFilterType::Pointer m_MultiplyTikhonovFilter;
229  typename AddFilterType::Pointer m_AddLaplacianFilter;
230  typename AddFilterType::Pointer m_AddTikhonovFilter;
231  typename itk::ImageToImageFilter<TOutputImage, TOutputImage>::Pointer m_LaplacianFilter;
232  typename MultiplyWithWeightsFilterType::Pointer m_MultiplyWithWeightsFilter;
233 
236  float m_Gamma{ 0 }; // Strength of the laplacian regularization
237  float m_Tikhonov{ 0 }; // Strength of the Tikhonov regularization
238 
240  typename TOutputImage::Pointer m_FloatingInputPointer, m_FloatingOutputPointer;
241 
244  void
245  VerifyInputInformation() const override
246  {}
247 
249  void
250  GenerateInputRequestedRegion() override;
251  void
252  GenerateOutputInformation() override;
254 
256  typename TOutputImage::ConstPointer
257  GetInputVolume();
258  typename TOutputImage::ConstPointer
259  GetInputProjectionStack();
260  typename TWeightsImage::ConstPointer
261  GetInputWeights();
262 };
263 } // namespace rtk
265 
266 
267 #ifndef ITK_MANUAL_INSTANTIATION
268 # include "rtkReconstructionConjugateGradientOperator.hxx"
269 #endif
270 
271 #endif
typename OutputImageType::Pointer OutputImagePointer
Implements the operator A used in conjugate gradient reconstruction.
Generate an n-dimensional image with constant pixel values.
STL namespace.
std::conditional< std::is_same< TSingleComponentImage, TOutputImage >::value, PlainMultiplyFilterType, MatrixVectorMultiplyFilterType >::type MultiplyWithWeightsFilterType
typename BackProjectionFilterType::Pointer BackProjectionFilterPointer
Projection geometry for a source and a 2-D flat panel.
#define itkSetMacro(name, type)
typename ForwardProjectionFilterType::Pointer ForwardProjectionFilterPointer