RTK  1.4.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  * 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 <math.h>
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_MAX_THREADS; i++)
49  {
50  m_AttenuationRay[i] = 0;
51  m_AttenuationPixel[i] = 0;
52  m_ex1[i] = 1;
53  }
54  }
55 
58  return false;
59  }
60 
62  {
63  return !( *this != other );
64  }
65 
66  inline TOutput operator()( const ThreadIdType threadId,
67  const double stepLengthInVoxel,
68  const TCoordRepType weight,
69  const TInput *p,
70  const int i )
71  {
72  const double w = weight*stepLengthInVoxel;
73 
76  return weight*p[i];
77  }
78 
80  TOutput * GetAttenuationRay() {return m_AttenuationRay;}
82  TOutput * GetEx1() {return m_ex1;}
83 
84 private:
89 };
90 
98 template< class TInput, class TOutput>
100 {
101 public:
103 
107  {
108  return false;
109  }
110 
111  bool operator==(const ComputeAttenuationCorrection & other) const
112  {
113  return !( *this != other );
114  }
115 
116  inline TOutput operator()(const ThreadIdType threadId,
117  const TInput volumeValue,
118  const VectorType &stepInMM)
119  {
120  TInput ex2 = exp(-m_AttenuationRay[threadId]*stepInMM.GetNorm() );
121  TInput wf;
122 
123  if(m_AttenuationPixel[threadId] > 0)
124  {
125  wf = (m_ex1[threadId]-ex2)/m_AttenuationPixel[threadId];
126  }
127  else
128  {
129  wf = m_ex1[threadId]*stepInMM.GetNorm();
130  }
131 
132  m_ex1[threadId] = ex2 ;
133  m_AttenuationPixel[threadId] = 0;
134  return wf *volumeValue;
135  }
136 
137  void SetAttenuationRayVector( TInput *attenuationRayVector) {m_AttenuationRay = attenuationRayVector;}
138  void SetAttenuationPixelVector( TInput *attenuationPixelVector) {m_AttenuationPixel = attenuationPixelVector;}
139  void SetEx1( TInput *ex1) {m_ex1 = ex1;}
140 
141 private:
144  TInput* m_ex1;
145 };
146 
154 template< class TInput, class TOutput >
156 {
157 public:
159 
163  {
164  return false;
165  }
166 
168  {
169  return !( *this != other );
170  }
171 
172  inline void operator()( const ThreadIdType threadId,
173  const TInput &input,
174  TOutput &output,
175  const TOutput &rayCastValue,
176  const VectorType & /*stepInMM*/,
177  const VectorType &itkNotUsed(source),
178  const VectorType &itkNotUsed(sourceToPixel),
179  const VectorType &itkNotUsed(nearestPoint),
180  const VectorType &itkNotUsed(farthestPoint) )
181  {
182  output = input + rayCastValue ;
183  m_Attenuation[threadId] = 0;
184  m_ex1[threadId] = 1;
185  }
186 
187  void SetAttenuationVector( TInput *attenuationVector) {m_Attenuation = attenuationVector;}
188  void SetEx1( TInput *ex1) {m_ex1 = ex1;}
189 
190 private:
191  TInput* m_Attenuation;
192  TInput* m_ex1;
193 };
194 } // end namespace Functor
195 
211 template <class TInputImage,
212  class TOutputImage,
216  >
218  public JosephForwardProjectionImageFilter<TInputImage,TOutputImage,TInterpolationWeightMultiplication, TProjectedValueAccumulation, TSumAlongRay>
219 {
220 public:
226  typedef typename TInputImage::PixelType InputPixelType;
227  typedef typename TOutputImage::PixelType OutputPixelType;
228  typedef typename TOutputImage::RegionType OutputImageRegionType;
229  typedef double CoordRepType;
231 
233  static constexpr unsigned int InputImageDimension = TInputImage::ImageDimension;
234 
236  itkNewMacro(Self);
237 
240 
241 protected:
244 
246  void GenerateInputRequestedRegion() ITK_OVERRIDE;
247 
248  void BeforeThreadedGenerateData() ITK_OVERRIDE;
249 
252  #if ITK_VERSION_MAJOR<5
253  void VerifyInputInformation() ITK_OVERRIDE;
254 #else
255  void VerifyInputInformation() const ITK_OVERRIDE;
256 #endif
257 
258 
259 private:
260  JosephForwardAttenuatedProjectionImageFilter(const Self&); //purposely not implemented
261  void operator=(const Self&); //purposely not implemented
262 };
263 } // end namespace rtk
264 
265 #ifndef ITK_MANUAL_INSTANTIATION
266 #include "rtkJosephForwardAttenuatedProjectionImageFilter.hxx"
267 #endif
268 
269 #endif
bool operator==(const ProjectedValueAccumulationAttenuated &other) const
bool operator!=(const ComputeAttenuationCorrection &) const
bool operator!=(const InterpolationWeightMultiplicationAttenuated &) const
Function to accumulate the ray casting on the projection.
itk::Vector< CoordRepType, TInputImage::ImageDimension > VectorType
constexpr std::vcl_size_t ITK_MAX_THREADS
JosephForwardProjectionImageFilter< TInputImage, TOutputImage, TInterpolationWeightMultiplication, TProjectedValueAccumulation, TSumAlongRay > Superclass
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.