-
Notifications
You must be signed in to change notification settings - Fork 143
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Implementation of variance reconstruction #404
base: master
Are you sure you want to change the base?
Changes from all commits
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -128,6 +128,10 @@ class ITK_TEMPLATE_EXPORT FFTRampImageFilter | |
itkGetConstMacro(SheppLoganCutFrequency, double); | ||
itkSetMacro(SheppLoganCutFrequency, double); | ||
|
||
/** Set/Get the previous kernel update size */ | ||
itkGetConstMacro(PreviousKernelUpdateSize, SizeType); | ||
itkSetMacro(PreviousKernelUpdateSize, SizeType); | ||
Comment on lines
+131
to
+133
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Should it be protected? Any reason to expose this parameter to the user? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Setting the variable |
||
|
||
protected: | ||
FFTRampImageFilter(); | ||
~FFTRampImageFilter() override = default; | ||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,84 @@ | ||
/*========================================================================= | ||
* | ||
* Copyright RTK Consortium | ||
* | ||
* Licensed under the Apache License, Version 2.0 (the "License"); | ||
* you may not use this file except in compliance with the License. | ||
* You may obtain a copy of the License at | ||
* | ||
* http://www.apache.org/licenses/LICENSE-2.0.txt | ||
* | ||
* Unless required by applicable law or agreed to in writing, software | ||
* distributed under the License is distributed on an "AS IS" BASIS, | ||
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. | ||
* See the License for the specific language governing permissions and | ||
* limitations under the License. | ||
* | ||
*=========================================================================*/ | ||
|
||
#ifndef rtkFFTVarianceRampImageFilter_h | ||
#define rtkFFTVarianceRampImageFilter_h | ||
|
||
#include "rtkFFTRampImageFilter.h" | ||
#include "rtkFFTProjectionsConvolutionImageFilter.h" | ||
#include "itkHalfHermitianToRealInverseFFTImageFilter.h" | ||
|
||
namespace rtk | ||
{ | ||
|
||
/** \class FFTVarianceRampImageFilter | ||
* \brief Implements the variance image filter of the filtered backprojection algorithm. | ||
*/ | ||
|
||
template<class TInputImage, class TOutputImage=TInputImage, class TFFTPrecision=double> | ||
class ITK_EXPORT FFTVarianceRampImageFilter : | ||
public rtk::FFTRampImageFilter<TInputImage, TOutputImage, TFFTPrecision> | ||
{ | ||
public: | ||
typedef rtk::FFTRampImageFilter<TInputImage, TOutputImage, TFFTPrecision> Baseclass; | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This is Superclass in ITK, should be removed |
||
|
||
/** Standard class typedefs. */ | ||
typedef FFTVarianceRampImageFilter Self; | ||
typedef rtk::FFTRampImageFilter< TInputImage, | ||
TOutputImage, | ||
TFFTPrecision> Superclass; | ||
typedef itk::SmartPointer<Self> Pointer; | ||
typedef itk::SmartPointer<const Self> ConstPointer; | ||
|
||
/** Some convenient typedefs. */ | ||
typedef typename Baseclass::InputImageType InputImageType; | ||
typedef typename Baseclass::OutputImageType OutputImageType; | ||
typedef typename Baseclass::FFTPrecisionType FFTPrecisionType; | ||
typedef typename Baseclass::IndexType IndexType; | ||
typedef typename Baseclass::SizeType SizeType; | ||
|
||
typedef typename Baseclass::FFTInputImageType FFTInputImageType; | ||
typedef typename Baseclass::FFTInputImagePointer FFTInputImagePointer; | ||
typedef typename Baseclass::FFTOutputImageType FFTOutputImageType; | ||
typedef typename Baseclass::FFTOutputImagePointer FFTOutputImagePointer; | ||
|
||
/** Standard New method. */ | ||
itkNewMacro(Self); | ||
|
||
/** Runtime information support. */ | ||
itkTypeMacro(FFTVarianceRampImageFilter, FFTConvolutionImageFilter); | ||
|
||
protected: | ||
FFTVarianceRampImageFilter(); | ||
~FFTVarianceRampImageFilter() {} | ||
|
||
/** Creates and return a pointer to one line of the variance kernel in Fourier space. | ||
* Used in generate data functions. */ | ||
void UpdateFFTProjectionsConvolutionKernel(const SizeType size) ITK_OVERRIDE; | ||
|
||
FFTVarianceRampImageFilter(const Self&); //purposely not implemented | ||
void operator=(const Self&); //purposely not implemented | ||
}; // end of class | ||
|
||
} // end namespace rtk | ||
|
||
#ifndef ITK_MANUAL_INSTANTIATION | ||
#include "rtkFFTVarianceRampImageFilter.hxx" | ||
#endif | ||
|
||
#endif |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,221 @@ | ||
/*========================================================================= | ||
* | ||
* Copyright RTK Consortium | ||
* | ||
* Licensed under the Apache License, Version 2.0 (the "License"); | ||
* you may not use this file except in compliance with the License. | ||
* You may obtain a copy of the License at | ||
* | ||
* http://www.apache.org/licenses/LICENSE-2.0.txt | ||
* | ||
* Unless required by applicable law or agreed to in writing, software | ||
* distributed under the License is distributed on an "AS IS" BASIS, | ||
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. | ||
* See the License for the specific language governing permissions and | ||
* limitations under the License. | ||
* | ||
*=========================================================================*/ | ||
|
||
#ifndef rtkFFTVarianceRampImageFilter_hxx | ||
#define rtkFFTVarianceRampImageFilter_hxx | ||
|
||
namespace rtk | ||
{ | ||
|
||
template <class TInputImage, class TOutputImage, class TFFTPrecision> | ||
FFTVarianceRampImageFilter<TInputImage, TOutputImage, TFFTPrecision>::FFTVarianceRampImageFilter() = default; | ||
|
||
template<class TInputImage, class TOutputImage, class TFFTPrecision> | ||
void | ||
FFTVarianceRampImageFilter<TInputImage, TOutputImage, TFFTPrecision> | ||
::UpdateFFTProjectionsConvolutionKernel(const SizeType s) | ||
{ | ||
if(this->m_KernelFFT.GetPointer() != nullptr && s == this->GetPreviousKernelUpdateSize()) | ||
{ | ||
return; | ||
} | ||
this->SetPreviousKernelUpdateSize(s); | ||
|
||
const int width = s[0]; | ||
const int height = s[1]; | ||
|
||
// Allocate kernel | ||
SizeType size; | ||
size.Fill(1); | ||
size[0] = width; | ||
FFTInputImagePointer kernel = FFTInputImageType::New(); | ||
kernel->SetRegions( size ); | ||
kernel->Allocate(); | ||
kernel->FillBuffer(0.); | ||
|
||
// Compute kernel in space domain (see Kak & Slaney, chapter 3 equation 61 | ||
// page 72) although spacing is not squared according to equation 69 page 75 | ||
double spacing = this->GetInput()->GetSpacing()[0]; | ||
IndexType ix,jx; | ||
ix.Fill(0); | ||
jx.Fill(0); | ||
kernel->SetPixel(ix, 1./(4.*spacing) ); | ||
for(ix[0]=1, jx[0]=size[0]-1; ix[0] < typename IndexType::IndexValueType(size[0]/2); ix[0] += 2, jx[0] -= 2) | ||
{ | ||
double v = ix[0] * itk::Math::pi; | ||
v = -1. / (v * v * spacing); | ||
kernel->SetPixel(ix, v); | ||
kernel->SetPixel(jx, v); | ||
} | ||
|
||
// FFT kernel | ||
using FFTType = itk::RealToHalfHermitianForwardFFTImageFilter< FFTInputImageType, FFTOutputImageType >; | ||
typename FFTType::Pointer fftK = FFTType::New(); | ||
fftK->SetInput( kernel ); | ||
fftK->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() ); | ||
fftK->Update(); | ||
this->m_KernelFFT = fftK->GetOutput(); | ||
|
||
// Windowing (if enabled) | ||
using IteratorType = itk::ImageRegionIteratorWithIndex<typename FFTType::OutputImageType>; | ||
IteratorType itK(this->m_KernelFFT, this->m_KernelFFT->GetLargestPossibleRegion() ); | ||
|
||
unsigned int n = this->m_KernelFFT->GetLargestPossibleRegion().GetSize(0); | ||
|
||
itK.GoToBegin(); | ||
if (this->GetHannCutFrequency() > 0.) | ||
{ | ||
const unsigned int ncut = itk::Math::Round<double>(n * std::min(1.0, this->GetHannCutFrequency())); | ||
for (unsigned int i = 0; i < ncut; i++, ++itK) | ||
itK.Set(itK.Get() * TFFTPrecision(0.5 * (1 + std::cos(itk::Math::pi * i / ncut)))); | ||
} | ||
else if (this->GetCosineCutFrequency() > 0.) | ||
{ | ||
const unsigned int ncut = itk::Math::Round<double>(n * std::min(1.0, this->GetCosineCutFrequency())); | ||
for (unsigned int i = 0; i < ncut; i++, ++itK) | ||
itK.Set(itK.Get() * TFFTPrecision(std::cos(0.5 * itk::Math::pi * i / ncut))); | ||
} | ||
else if (this->GetHammingFrequency() > 0.) | ||
{ | ||
const unsigned int ncut = itk::Math::Round<double>(n * std::min(1.0, this->GetHammingFrequency())); | ||
for (unsigned int i = 0; i < ncut; i++, ++itK) | ||
itK.Set(itK.Get() * TFFTPrecision(0.54 + 0.46 * (std::cos(itk::Math::pi * i / ncut)))); | ||
} | ||
else if (this->GetRamLakCutFrequency() > 0.) | ||
{ | ||
const unsigned int ncut = itk::Math::Round<double>(n * std::min(1.0, this->GetRamLakCutFrequency())); | ||
for (unsigned int i = 0; i < ncut; i++, ++itK) | ||
{ | ||
} | ||
} | ||
else if (this->GetSheppLoganCutFrequency() > 0.) | ||
{ | ||
const unsigned int ncut = itk::Math::Round<double>(n * std::min(1.0, this->GetSheppLoganCutFrequency())); | ||
// sinc(0) --> is 1 | ||
++itK; | ||
for (unsigned int i = 1; i < ncut; i++, ++itK) | ||
{ | ||
double x = 0.5 * itk::Math::pi * i / ncut; | ||
itK.Set(itK.Get() * TFFTPrecision(std::sin(x) / x)); | ||
} | ||
} | ||
else | ||
{ | ||
itK.GoToReverseBegin(); | ||
++itK; | ||
} | ||
|
||
for(; !itK.IsAtEnd(); ++itK) | ||
{ | ||
itK.Set( itK.Get() * TFFTPrecision(0.) ); | ||
} | ||
|
||
// iFFT kernel | ||
using IFFTType = itk::HalfHermitianToRealInverseFFTImageFilter< FFTOutputImageType, FFTInputImageType >; | ||
typename IFFTType::Pointer ifftK = IFFTType::New(); | ||
ifftK->SetInput( this->m_KernelFFT ); | ||
ifftK->SetNumberOfThreads( this->GetNumberOfThreads() ); | ||
ifftK->Update(); | ||
typename FFTType::InputImageType::Pointer KernelIFFT = ifftK->GetOutput(); | ||
|
||
// calculate ratio g_C/g^2 | ||
using InverseIteratorType = itk::ImageRegionIteratorWithIndex<typename FFTType::InputImageType>; | ||
InverseIteratorType itiK(KernelIFFT, KernelIFFT->GetLargestPossibleRegion()); | ||
typename FFTType::InputImageType::PixelType sumgc = 0.; | ||
typename FFTType::InputImageType::PixelType sumg2 = 0.; | ||
typename FFTType::InputImageType::IndexType idxshifted; | ||
unsigned int const widthFFT = KernelIFFT->GetLargestPossibleRegion().GetSize()[0]; | ||
for(; !itiK.IsAtEnd(); ++itiK) | ||
{ | ||
idxshifted = itiK.GetIndex(); | ||
if(idxshifted[0] == 0) | ||
idxshifted[0] = widthFFT-1; | ||
else | ||
idxshifted[0] -= 1; | ||
|
||
sumgc += itiK.Get() * KernelIFFT->GetPixel(idxshifted); | ||
sumg2 += itiK.Get() * itiK.Get(); | ||
} | ||
typename FFTType::InputImageType::PixelType const ratiogcg2 = sumgc/sumg2; | ||
|
||
// numerical integration to calculate f_interp | ||
double const aprecision = 0.00001; | ||
double finterp = 0.; | ||
for(unsigned int i=0; i<int(1./aprecision); i++) | ||
{ | ||
double const a = double(i)*aprecision; | ||
finterp += (1-a)*(1-a) + 2*ratiogcg2*(1-a)*a + a*a; | ||
} | ||
finterp *= aprecision; | ||
|
||
// square kernel and multiply with finterp | ||
itiK.GoToBegin(); | ||
for(; !itiK.IsAtEnd(); ++itiK) | ||
{ | ||
itiK.Set(itiK.Get() * itiK.Get() * finterp); | ||
} | ||
|
||
// FFT kernel | ||
typename FFTType::Pointer fftK2 = FFTType::New(); | ||
fftK2->SetInput( ifftK->GetOutput() ); | ||
fftK2->SetNumberOfThreads( this->GetNumberOfThreads() ); | ||
fftK2->Update(); | ||
this->m_KernelFFT = fftK2->GetOutput(); | ||
|
||
Comment on lines
+128
to
+179
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This part is the only difference. Maybe we should consider putting this in a function called by the parent, which would be empty for the parent and redefined for this class? Then the reimplementation of I also wonder, could this have been done after the application of the windowing and finterp kept to the original value without the numerical integration? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. For the first part, I totally agree. It's nicer to just implement the part that actually changes. Splitting I'm not sure if I understood the second part right, but I'll try to motivate what I did here. To calculate Right now the code does the following:
I am not sure if the many times going back and forth between real space and Fourier space can cause problems (like the filter size change I fixed in this commit). Maybe it can be avoided by re-defining the squared ramp in real space? |
||
// Replicate and window if required | ||
if(this->GetHannCutFrequencyY()>0.) | ||
{ | ||
std::cerr << "HannY not implemented for variance image filter." << std::endl; | ||
|
||
size.Fill(1); | ||
size[0] = this->m_KernelFFT->GetLargestPossibleRegion().GetSize(0); | ||
size[1] = height; | ||
|
||
const unsigned int ncut = itk::Math::Round<double>( (height/2+1) * std::min(1.0, this->GetHannCutFrequencyY() ) ); | ||
|
||
this->m_KernelFFT = FFTOutputImageType::New(); | ||
this->m_KernelFFT->SetRegions( size ); | ||
this->m_KernelFFT->Allocate(); | ||
this->m_KernelFFT->FillBuffer(0.); | ||
|
||
IteratorType itTwoDK(this->m_KernelFFT, this->m_KernelFFT->GetLargestPossibleRegion() ); | ||
for(unsigned int j=0; j<ncut; j++) | ||
{ | ||
itK.GoToBegin(); | ||
const TFFTPrecision win( 0.5*( 1+std::cos(itk::Math::pi*j/ncut) ) ); | ||
for(unsigned int i=0; i<size[0]; ++itK, ++itTwoDK, i++) | ||
{ | ||
itTwoDK.Set( win * itK.Get() ); | ||
} | ||
} | ||
itTwoDK.GoToReverseBegin(); | ||
for(unsigned int j=1; j<ncut; j++) | ||
{ | ||
itK.GoToReverseBegin(); | ||
const TFFTPrecision win( 0.5*( 1+std::cos(itk::Math::pi*j/ncut) ) ); | ||
for(unsigned int i=0; i<size[0]; --itK, --itTwoDK, i++) | ||
{ | ||
itTwoDK.Set( win * itK.Get() ); | ||
} | ||
} | ||
} | ||
this->m_KernelFFT->DisconnectPipeline(); | ||
} | ||
|
||
} // end namespace rtk | ||
#endif |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm surprised you only square the constant part of the weight. Has this been tested with a divergent geometry? For simplicity, my suggestion would be to apply this filter twice instead of adding a parameters.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Indeed we have almost always used parallel geometries. We have also tested it with a divergent geometry and errors were few percent, but those could of course increase if the SID is not 2m as for pCT, but smaller. I like the idea of applying the filter twice. This would also avoid the special case with
m_VarianceMode
. However, there should be some note in thertkFFTVarianceRampImageFilter.h
that the weights filter needs to be applied twice.