RTK  2.6.0
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  * 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 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 namespace rtk
31 {
32 namespace Functor
33 {
42 template <class TInput, class TCoordRepType, class TOutput = TInput>
44 {
45 public:
47  {
48  for (std::size_t i = 0; i < itk::ITK_MAX_THREADS; i++)
49  {
50  m_AttenuationRay[i] = 0;
51  m_AttenuationPixel[i] = 0;
52  m_Ex1[i] = 1;
53  }
54  }
55 
57  bool
59  {
60  return false;
61  }
62 
63  bool
65  {
66  return !(*this != other);
67  }
68 
69  inline TOutput
70  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 
78  m_AttenuationRay[threadId] += w * (p + m_AttenuationMinusEmissionMapsPtrDiff)[i];
79  m_AttenuationPixel[threadId] += w * (p + m_AttenuationMinusEmissionMapsPtrDiff)[i];
80  return weight * p[i];
81  }
82 
83  void
85  {
86  m_AttenuationMinusEmissionMapsPtrDiff = pd;
87  }
88  TOutput *
90  {
91  return m_AttenuationRay;
92  }
93  TOutput *
95  {
96  return m_AttenuationPixel;
97  }
98  TOutput *
100  {
101  return m_Ex1;
102  }
103 
104 private:
106  TInput m_AttenuationRay[itk::ITK_MAX_THREADS];
107  TInput m_AttenuationPixel[itk::ITK_MAX_THREADS];
108  TInput m_Ex1[itk::ITK_MAX_THREADS];
109 };
110 
118 template <class TInput, class TOutput>
119 class ITK_TEMPLATE_EXPORT ComputeAttenuationCorrection
120 {
121 public:
123 
124  ComputeAttenuationCorrection() = default;
125  ~ComputeAttenuationCorrection() = default;
126  bool
128  {
129  return false;
130  }
131 
132  bool
134  {
135  return !(*this != other);
136  }
137 
138  inline void
139  operator()(const ThreadIdType threadId, TOutput & sumValue, const TInput volumeValue, const VectorType & stepInMM)
140  {
141  TInput ex2 = exp(-m_AttenuationRay[threadId] * stepInMM.GetNorm());
142  TInput wf;
143 
144  if (m_AttenuationPixel[threadId] > 0)
145  {
146  wf = (m_Ex1[threadId] - ex2) / m_AttenuationPixel[threadId];
147  }
148  else
149  {
150  wf = m_Ex1[threadId] * stepInMM.GetNorm();
151  }
152 
153  m_Ex1[threadId] = ex2;
154  m_AttenuationPixel[threadId] = 0;
155  sumValue += wf * volumeValue;
156  }
157 
158  void
159  SetAttenuationRayVector(TInput * attenuationRayVector)
160  {
161  m_AttenuationRay = attenuationRayVector;
162  }
163  void
164  SetAttenuationPixelVector(TInput * attenuationPixelVector)
165  {
166  m_AttenuationPixel = attenuationPixelVector;
167  }
168  void
169  SetEx1(TInput * ex1)
170  {
171  m_Ex1 = ex1;
172  }
173 
174 private:
177  TInput * m_Ex1;
178 };
179 
187 template <class TInput, class TOutput>
188 class ITK_TEMPLATE_EXPORT ProjectedValueAccumulationAttenuated
189 {
190 public:
192 
195  bool
197  {
198  return false;
199  }
200 
201  bool
203  {
204  return !(*this != other);
205  }
206 
207  inline void
208  operator()(const ThreadIdType threadId,
209  const TInput & input,
210  TOutput & output,
211  const TOutput & rayCastValue,
212  const VectorType & /*stepInMM*/,
213  const VectorType & itkNotUsed(source),
214  const VectorType & itkNotUsed(sourceToPixel),
215  const VectorType & itkNotUsed(nearestPoint),
216  const VectorType & itkNotUsed(farthestPoint))
217  {
218  output = input + rayCastValue;
219  m_Attenuation[threadId] = 0;
220  m_Ex1[threadId] = 1;
221  }
222 
223  void
224  SetAttenuationVector(TInput * attenuationVector)
225  {
226  m_Attenuation = attenuationVector;
227  }
228  void
229  SetEx1(TInput * ex1)
230  {
231  m_Ex1 = ex1;
232  }
233 
234 private:
235  TInput * m_Attenuation;
236  TInput * m_Ex1;
237 };
238 } // end namespace Functor
239 
255 template <
256  class TInputImage,
257  class TOutputImage,
258  class TInterpolationWeightMultiplication = Functor::InterpolationWeightMultiplicationAttenuated<
259  typename TInputImage::PixelType,
261  class TProjectedValueAccumulation =
263  class TSumAlongRay =
266  : public JosephForwardProjectionImageFilter<TInputImage,
267  TOutputImage,
268  TInterpolationWeightMultiplication,
269  TProjectedValueAccumulation,
270  TSumAlongRay>
271 {
272 public:
273  ITK_DISALLOW_COPY_AND_MOVE(JosephForwardAttenuatedProjectionImageFilter);
274 
278  TOutputImage,
279  TInterpolationWeightMultiplication,
280  TProjectedValueAccumulation,
281  TSumAlongRay>;
284  using InputPixelType = typename TInputImage::PixelType;
285  using OutputPixelType = typename TOutputImage::PixelType;
286  using OutputImageRegionType = typename TOutputImage::RegionType;
287  using CoordRepType = double;
289 
291  static constexpr unsigned int InputImageDimension = TInputImage::ImageDimension;
292 
294  itkNewMacro(Self);
295 
297 #ifdef itkOverrideGetNameOfClassMacro
298  itkOverrideGetNameOfClassMacro(JosephForwardAttenuatedProjectionImageFilter);
299 #else
301 #endif
302 
303 
304 protected:
306  ~JosephForwardAttenuatedProjectionImageFilter() override = default;
307 
309  void
310  GenerateInputRequestedRegion() override;
311 
312  void
313  BeforeThreadedGenerateData() override;
314 
317  void
318  VerifyInputInformation() const override;
319 };
320 } // end namespace rtk
321 
322 #ifndef ITK_MANUAL_INSTANTIATION
323 # include "rtkJosephForwardAttenuatedProjectionImageFilter.hxx"
324 #endif
325 
326 #endif
constexpr vcl_size_t ITK_MAX_THREADS
Function to accumulate the ray casting on the projection.
bool operator!=(const ProjectedValueAccumulationAttenuated &) const
typename TPixelType::ValueType ValueType
typename OutputImageType::RegionType OutputImageRegionType
bool operator==(const ProjectedValueAccumulationAttenuated &other) const
bool operator!=(const InterpolationWeightMultiplicationAttenuated &) const
void operator()(const ThreadIdType threadId, TOutput &sumValue, const TInput volumeValue, const VectorType &stepInMM)
bool operator==(const ComputeAttenuationCorrection &other) const
unsigned int ThreadIdType
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 &)
bool operator==(const InterpolationWeightMultiplicationAttenuated &other) const
Function to compute the attenuation correction on the projection.