Implementation of variance reconstruction #404
Conversation
…mageFilter and modifications to FFTRampImageFilter and FDKWeightProjectionFilter.
|
Thanks Jannis. I would like to add a test and I can try to implement a functional one. If I want to reconstruct the variance, should I compute variance projections, then FDKWeightProjectionFilter (VarianceMode) -> FFTVarianceRampImageFilter -> FDKBackProjectionImageFilter ? |
|
Hi Simon. Yes that's correct. In the FFTVarianceRampImageFilter we should probably also test for the correct calculation of finterp, which should be 0.46 without Hann filter, and I can provide expected values for other Hann values. The HannY option is currently not considered and I was just printing an error message, which probably is not ideal. |
| { | ||
| idxshifted = itiK.GetIndex(); | ||
| if(idxshifted[0] == 0) | ||
| idxshifted[0] = width-1; |
There was a problem hiding this comment.
Using width can actually cause problems, since the inverse Fourier of the Fourier transform is not guaranteed to have the size (see https://itk.org/Doxygen/html/classitk_1_1HalfHermitianToRealInverseFFTImageFilter.html: "the input is assumed to consist of roughly half the full complex image resulting from a real-to-complex discrete Fourier transform").
One solution would be to query KernelIFFT directly for its length. The other solution may be to use SetActualDimensionIsOdd() on the FFT filter (https://itk.org/Doxygen/html/classitk_1_1HalfHermitianToRealInverseFFTImageFilter.html#af072dcf520dee228202a11b6e99597f5) although I would have to test how this works.
Do you have any preference?
There was a problem hiding this comment.
No preference, no. Like you, I don't know what are the cons of using SetActualDimensionIsOdd(). Pick the solution you prefer!
There was a problem hiding this comment.
I added a new commit where I determine the size of the Fourier-transformed filter again instead of relying on the width variable. This seems to me the safest option.
…nce reconstruction filter
| if(this->m_VarianceMode) | ||
| { | ||
| // square if in variance mode | ||
| m_ConstantProjectionFactor[k] = m_ConstantProjectionFactor[k] * m_ConstantProjectionFactor[k]; | ||
| } |
There was a problem hiding this comment.
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.
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 the rtkFFTVarianceRampImageFilter.h that the weights filter needs to be applied twice.
| /** Set/Get the previous kernel update size */ | ||
| itkGetConstMacro(PreviousKernelUpdateSize, SizeType); | ||
| itkSetMacro(PreviousKernelUpdateSize, SizeType); |
There was a problem hiding this comment.
Should it be protected? Any reason to expose this parameter to the user?
There was a problem hiding this comment.
Setting the variable protected is better. I don't see a reason why it needs to be exposed. It just needs to be accessible here.
| public rtk::FFTRampImageFilter<TInputImage, TOutputImage, TFFTPrecision> | ||
| { | ||
| public: | ||
| typedef rtk::FFTRampImageFilter<TInputImage, TOutputImage, TFFTPrecision> Baseclass; |
There was a problem hiding this comment.
This is Superclass in ITK, should be removed
| // 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(); | ||
|
|
There was a problem hiding this comment.
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 UpdateFFTProjectionsConvolutionKernel can be removed as well as the getter and setter for PreviousKernelUpdateSize.
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.
For the first part, I totally agree. It's nicer to just implement the part that actually changes. Splitting UpdateFFTProjectionsConvolutionKernel of rtkFFTRampImageFilter.h into several steps would also make it more flexible for other RTK users that just want to modfiy one component of the filter.
I'm not sure if I understood the second part right, but I'll try to motivate what I did here.
To calculate finterp one needs access to the filter including the windowing function in real space. finterp is then an integral that depends on the filters g^2 and g_C, which in turn depend on the windowed filter in real space. The implementation that I propose is computationally most expensive, but also most flexible, since it can take any windowed ramp filter. Alternatively I could fit an function finterp(Hann) and avoid the numerical integration and going back and forth between Fourier and real space.
Right now the code does the following:
- define ramp in real space
- go to Fourier space
- apply window function
- go back to real space
- calculate
finterp - square filter and multiply by
finterp - go to Fourier space
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?
Implementation of variance reconstruction with a new FFTVarianceRampImageFilter and modifications to FFTRampImageFilter and FDKWeightProjectionFilter.