RTK  2.0.1
Reconstruction Toolkit
rtkJosephForwardAttenuatedProjectionImageFilter.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 rtkJosephForwardAttenuatedProjectionImageFilter_h
20 #define rtkJosephForwardAttenuatedProjectionImageFilter_h
21 
22 #include "rtkConfiguration.h"
25 #include "rtkMacro.h"
26 #include <itkPixelTraits.h>
27 #include <cmath>
28 #include <vector>
29 
30 #if ITK_VERSION_MAJOR>4 && defined ( ITK_FUTURE_LEGACY_REMOVE )
31  using itk::ITK_MAX_THREADS;
32 #endif
33 
34 namespace rtk
35 {
36 namespace Functor
37 {
46 template< class TInput, class TCoordRepType, class TOutput = TInput >
48 {
49 public:
51  {
52  for (std::size_t i = 0; i < ITK_MAX_THREADS; i++)
53  {
54  m_AttenuationRay[i] = 0;
55  m_AttenuationPixel[i] = 0;
56  m_ex1[i] = 1;
57  }
58  }
59 
62  return false;
63  }
64 
66  {
67  return !( *this != other );
68  }
69 
70  inline TOutput operator()( const ThreadIdType threadId,
71  const double stepLengthInVoxel,
72  const TCoordRepType weight,
73  const TInput *p,
74  const int i )
75  {
76  const double w = weight*stepLengthInVoxel;
77 
80  return weight*p[i];
81  }
82 
84  TOutput * GetAttenuationRay() {return m_AttenuationRay;}
86  TOutput * GetEx1() {return m_ex1;}
87 
88 private:
93 };
94 
102 template< class TInput, class TOutput>
104 {
105 public:
107 
108  ComputeAttenuationCorrection()= default;
109  ~ComputeAttenuationCorrection() = default;
111  {
112  return false;
113  }
114 
115  bool operator==(const ComputeAttenuationCorrection & other) const
116  {
117  return !( *this != other );
118  }
119 
120  inline TOutput operator()(const ThreadIdType threadId,
121  const TInput volumeValue,
122  const VectorType &stepInMM)
123  {
124  TInput ex2 = exp(-m_AttenuationRay[threadId]*stepInMM.GetNorm() );
125  TInput wf;
126 
127  if(m_AttenuationPixel[threadId] > 0)
128  {
129  wf = (m_ex1[threadId]-ex2)/m_AttenuationPixel[threadId];
130  }
131  else
132  {
133  wf = m_ex1[threadId]*stepInMM.GetNorm();
134  }
135 
136  m_ex1[threadId] = ex2;
137  m_AttenuationPixel[threadId] = 0;
138  return wf *volumeValue;
139  }
140 
141  void SetAttenuationRayVector( TInput *attenuationRayVector) {m_AttenuationRay = attenuationRayVector;}
142  void SetAttenuationPixelVector( TInput *attenuationPixelVector) {m_AttenuationPixel = attenuationPixelVector;}
143  void SetEx1( TInput *ex1) {m_ex1 = ex1;}
144 
145 private:
148  TInput* m_ex1;
149 };
150 
158 template< class TInput, class TOutput >
160 {
161 public:
163 
167  {
168  return false;
169  }
170 
172  {
173  return !( *this != other );
174  }
175 
176  inline void operator()( const ThreadIdType threadId,
177  const TInput &input,
178  TOutput &output,
179  const TOutput &rayCastValue,
180  const VectorType & /*stepInMM*/,
181  const VectorType &itkNotUsed(source),
182  const VectorType &itkNotUsed(sourceToPixel),
183  const VectorType &itkNotUsed(nearestPoint),
184  const VectorType &itkNotUsed(farthestPoint) )
185  {
186  output = input + rayCastValue;
187  m_Attenuation[threadId] = 0;
188  m_ex1[threadId] = 1;
189  }
190 
191  void SetAttenuationVector( TInput *attenuationVector) {m_Attenuation = attenuationVector;}
192  void SetEx1( TInput *ex1) {m_ex1 = ex1;}
193 
194 private:
195  TInput* m_Attenuation;
196  TInput* m_ex1;
197 };
198 } // end namespace Functor
199 
215 template <class TInputImage,
216  class TOutputImage,
220  >
222  public JosephForwardProjectionImageFilter<TInputImage,TOutputImage,TInterpolationWeightMultiplication, TProjectedValueAccumulation, TSumAlongRay>
223 {
224 public:
225  ITK_DISALLOW_COPY_AND_ASSIGN(JosephForwardAttenuatedProjectionImageFilter);
226 
232  using InputPixelType = typename TInputImage::PixelType;
233  using OutputPixelType = typename TOutputImage::PixelType;
234  using OutputImageRegionType = typename TOutputImage::RegionType;
235  using CoordRepType = double;
237 
239  static constexpr unsigned int InputImageDimension = TInputImage::ImageDimension;
240 
242  itkNewMacro(Self);
243 
246 
247 protected:
249  ~JosephForwardAttenuatedProjectionImageFilter() override = default;
250 
252  void GenerateInputRequestedRegion() override;
253 
254  void BeforeThreadedGenerateData() override;
255 
258  #if ITK_VERSION_MAJOR<5
259  void VerifyInputInformation() override;
260 #else
261  void VerifyInputInformation() const override;
262 #endif
263 
264 
265 };
266 } // end namespace rtk
267 
268 #ifndef ITK_MANUAL_INSTANTIATION
269 #include "rtkJosephForwardAttenuatedProjectionImageFilter.hxx"
270 #endif
271 
272 #endif
bool operator==(const ProjectedValueAccumulationAttenuated &other) const
Base class for forward projection, i.e. accumulation along x-ray lines.
bool operator!=(const ComputeAttenuationCorrection &) const
bool operator!=(const InterpolationWeightMultiplicationAttenuated &) const
Function to accumulate the ray casting on the projection.
constexpr std::vcl_size_t ITK_MAX_THREADS
bool operator!=(const ProjectedValueAccumulationAttenuated &) const
bool operator==(const InterpolationWeightMultiplicationAttenuated &other) const
unsigned int ThreadIdType
TOutput operator()(const ThreadIdType threadId, const TInput volumeValue, const VectorType &stepInMM)
bool operator==(const ComputeAttenuationCorrection &other) const
TOutput operator()(const ThreadIdType threadId, const double stepLengthInVoxel, const TCoordRepType weight, const TInput *p, const int i)
Function to multiply the interpolation weights with the projected volume values and attenuation map...
void operator()(const ThreadIdType threadId, const TInput &input, TOutput &output, const TOutput &rayCastValue, const VectorType &, const VectorType &, const VectorType &, const VectorType &, const VectorType &)
Function to compute the attenuation correction on the projection.