Skip to content

Commit e6cb103

Browse files
astammastamm
committed
Forgotten files...
1 parent 39bc7c5 commit e6cb103

4 files changed

Lines changed: 712 additions & 0 deletions

File tree

‎Anima/filtering/rician_to_gaussian_noise/CMakeLists.txt‎

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,37 @@
1+
if(BUILD_TOOLS AND USE_ARB)
2+
3+
project(animaRicianToGaussianNoise)
4+
5+
## #############################################################################
6+
## List Sources
7+
## #############################################################################
8+
9+
list_source_files(${PROJECT_NAME}
10+
${CMAKE_CURRENT_SOURCE_DIR}
11+
)
12+
13+
## #############################################################################
14+
## add executable
15+
## #############################################################################
16+
17+
add_executable(${PROJECT_NAME}
18+
${${PROJECT_NAME}_CFILES}
19+
)
20+
21+
22+
## #############################################################################
23+
## Link
24+
## #############################################################################
25+
26+
target_link_libraries(${PROJECT_NAME}
27+
${ITKIO_LIBRARIES}
28+
AnimaSpecialFunctions
29+
)
30+
31+
## #############################################################################
32+
## install
33+
## #############################################################################
34+
35+
set_exe_install_rules(${PROJECT_NAME})
36+
37+
endif()

‎Anima/filtering/rician_to_gaussian_noise/animaRicianToGaussianImageFilter.h‎

Lines changed: 104 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,104 @@
1+
#pragma once
2+
3+
#include <animaMaskedImageToImageFilter.h>
4+
5+
#include <boost/math/distributions/normal.hpp>
6+
#include <boost/math/distributions/rayleigh.hpp>
7+
8+
namespace anima
9+
{
10+
template <class ComponentType,unsigned int ImageDimension>
11+
class RicianToGaussianImageFilter :
12+
public anima::MaskedImageToImageFilter<itk::Image<ComponentType,ImageDimension>,itk::Image<ComponentType,ImageDimension> >
13+
{
14+
public:
15+
/** Standard class typedefs. */
16+
typedef RicianToGaussianImageFilter Self;
17+
18+
typedef itk::Image<ComponentType,ImageDimension> InputImageType;
19+
typedef typename InputImageType::PixelType InputPixelType;
20+
21+
typedef itk::Image<ComponentType,ImageDimension> OutputImageType;
22+
typedef typename OutputImageType::PixelType OutputPixelType;
23+
24+
typedef anima::MaskedImageToImageFilter<InputImageType, OutputImageType> Superclass;
25+
26+
typedef itk::SmartPointer<Self> Pointer;
27+
typedef itk::SmartPointer<const Self> ConstPointer;
28+
29+
/** Method for creation through the object factory. */
30+
itkNewMacro(Self)
31+
32+
/** Run-time type information (and related methods) */
33+
itkTypeMacro(RicianToGaussianImageFilter, MaskedImageToImageFilter)
34+
35+
/** Superclass typedefs. */
36+
typedef typename Superclass::InputImageRegionType InputImageRegionType;
37+
typedef typename Superclass::OutputImageRegionType OutputImageRegionType;
38+
typedef typename Superclass::MaskImageType MaskImageType;
39+
typedef typename MaskImageType::SizeType SizeType;
40+
41+
typedef typename InputImageType::Pointer InputPointerType;
42+
typedef typename OutputImageType::Pointer OutputPointerType;
43+
44+
itkSetMacro(MaximumNumberOfIterations, unsigned int)
45+
itkSetMacro(Epsilon, double)
46+
itkSetMacro(Sigma, double)
47+
itkSetMacro(MeanImage, InputPointerType)
48+
itkSetMacro(VarianceImage, InputPointerType)
49+
50+
itkSetMacro(Scale, double)
51+
itkGetConstMacro(Scale, double)
52+
53+
itkSetMacro(Alpha, double)
54+
itkGetConstMacro(Alpha, double)
55+
56+
OutputPointerType GetLocationImage() {return this->GetOutput(0);}
57+
OutputPointerType GetScaleImage() {return this->GetOutput(1);}
58+
OutputPointerType GetGaussianImage() {return this->GetOutput(2);}
59+
60+
protected:
61+
RicianToGaussianImageFilter()
62+
: Superclass()
63+
{
64+
m_MaximumNumberOfIterations = 100;
65+
m_Epsilon = 1.0e-8;
66+
m_Sigma = 1.0;
67+
m_Scale = 0.0;
68+
m_Alpha = 0.005;
69+
m_ThreadScaleSamples.clear();
70+
m_NeighborWeights.clear();
71+
m_Radius.Fill(0);
72+
m_MeanImage = NULL;
73+
m_VarianceImage = NULL;
74+
75+
this->SetNumberOfRequiredOutputs(3);
76+
this->SetNthOutput(0, this->MakeOutput(0));
77+
this->SetNthOutput(1, this->MakeOutput(1));
78+
this->SetNthOutput(2, this->MakeOutput(2));
79+
}
80+
81+
virtual ~RicianToGaussianImageFilter() {}
82+
83+
void BeforeThreadedGenerateData(void) ITK_OVERRIDE;
84+
void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
85+
itk::ThreadIdType threadId) ITK_OVERRIDE;
86+
void AfterThreadedGenerateData(void) ITK_OVERRIDE;
87+
88+
private:
89+
ITK_DISALLOW_COPY_AND_ASSIGN(RicianToGaussianImageFilter);
90+
91+
unsigned int m_MaximumNumberOfIterations;
92+
double m_Epsilon, m_Sigma, m_Scale, m_Alpha;
93+
std::vector<std::vector<double> > m_ThreadScaleSamples;
94+
SizeType m_Radius;
95+
std::vector<double> m_NeighborWeights;
96+
InputPointerType m_MeanImage, m_VarianceImage;
97+
98+
boost::math::normal_distribution<> m_NormalDistribution;
99+
boost::math::rayleigh_distribution<> m_RayleighDistribution;
100+
};
101+
102+
} // end of namespace anima
103+
104+
#include "animaRicianToGaussianImageFilter.hxx"

‎Anima/filtering/rician_to_gaussian_noise/animaRicianToGaussianImageFilter.hxx‎

Lines changed: 247 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,247 @@
1+
#pragma once
2+
#include "animaRicianToGaussianImageFilter.h"
3+
#include <animaChiDistribution.h>
4+
#include <animaVectorOperations.h>
5+
6+
#include <itkConstNeighborhoodIterator.h>
7+
#include <itkImageRegionConstIterator.h>
8+
#include <itkImageRegionIterator.h>
9+
#include <itkGaussianOperator.h>
10+
11+
namespace anima
12+
{
13+
14+
template <class ComponentType,unsigned int ImageDimension>
15+
void
16+
RicianToGaussianImageFilter<ComponentType,ImageDimension>
17+
::BeforeThreadedGenerateData(void)
18+
{
19+
// Compute spatial weights of neighbors beforehand
20+
typename InputImageType::SpacingType spacing = this->GetInput()->GetSpacing();
21+
typedef itk::GaussianOperator<double,ImageDimension> GaussianOperatorType;
22+
std::vector<GaussianOperatorType> gaussianKernels(ImageDimension);
23+
24+
for (unsigned int i = 0;i < ImageDimension;++i)
25+
{
26+
unsigned int reverse_i = ImageDimension - i - 1;
27+
double stddev = m_Sigma / spacing[i];
28+
gaussianKernels[reverse_i].SetDirection(i);
29+
gaussianKernels[reverse_i].SetVariance(stddev * stddev);
30+
gaussianKernels[reverse_i].SetMaximumError(1e-3);
31+
gaussianKernels[reverse_i].CreateDirectional();
32+
gaussianKernels[reverse_i].ScaleCoefficients(1.0e4);
33+
m_Radius[i] = gaussianKernels[reverse_i].GetRadius(i);
34+
}
35+
36+
m_NeighborWeights.clear();
37+
for (unsigned int i = 0;i < gaussianKernels[0].Size();++i)
38+
{
39+
for (unsigned int j = 0;j < gaussianKernels[1].Size();++j)
40+
{
41+
for (unsigned int k = 0;k < gaussianKernels[2].Size();++k)
42+
{
43+
double weight = gaussianKernels[0][i] * gaussianKernels[1][j] * gaussianKernels[2][k];
44+
m_NeighborWeights.push_back(weight);
45+
}
46+
}
47+
}
48+
49+
// Initialize thread containers for global scale estimation
50+
unsigned int numThreads = this->GetNumberOfThreads();
51+
m_ThreadScaleSamples.resize(numThreads);
52+
this->Superclass::BeforeThreadedGenerateData();
53+
}
54+
55+
template <class ComponentType,unsigned int ImageDimension>
56+
void
57+
RicianToGaussianImageFilter<ComponentType,ImageDimension>
58+
::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
59+
itk::ThreadIdType threadId)
60+
{
61+
typedef itk::ConstNeighborhoodIterator<InputImageType> InputIteratorType;
62+
typedef itk::ConstNeighborhoodIterator<MaskImageType> MaskIteratorType;
63+
typedef itk::ImageRegionConstIterator<InputImageType> InputSimpleIteratorType;
64+
typedef itk::ImageRegionIterator<OutputImageType> OutputIteratorType;
65+
66+
InputIteratorType inputItr(m_Radius, this->GetInput(), outputRegionForThread);
67+
unsigned int neighborhoodSize = inputItr.Size();
68+
69+
MaskIteratorType maskItr(m_Radius, this->GetComputationMask(), outputRegionForThread);
70+
71+
InputSimpleIteratorType meanItr, varItr;
72+
if (m_MeanImage)
73+
meanItr = InputSimpleIteratorType(m_MeanImage, outputRegionForThread);
74+
if (m_VarianceImage)
75+
varItr = InputSimpleIteratorType(m_VarianceImage, outputRegionForThread);
76+
77+
OutputIteratorType locationItr(this->GetOutput(0), outputRegionForThread);
78+
OutputIteratorType scaleItr(this->GetOutput(1), outputRegionForThread);
79+
OutputIteratorType signalItr(this->GetOutput(2), outputRegionForThread);
80+
81+
std::vector<double> samples, weights;
82+
bool isInBounds;
83+
typename InputImageType::IndexType currentIndex, neighborIndex;
84+
typename InputImageType::PointType currentPoint, neighborPoint;
85+
86+
while (!maskItr.IsAtEnd())
87+
{
88+
// Discard voxels outside of the brain
89+
if (maskItr.GetCenterPixel() == 0)
90+
{
91+
locationItr.Set(0);
92+
scaleItr.Set(0);
93+
signalItr.Set(0);
94+
95+
++inputItr;
96+
++maskItr;
97+
if (m_MeanImage)
98+
++meanItr;
99+
if (m_VarianceImage)
100+
++varItr;
101+
++locationItr;
102+
++scaleItr;
103+
++signalItr;
104+
105+
continue;
106+
}
107+
108+
// Get signal at current central voxel
109+
double inputSignal = inputItr.GetCenterPixel();
110+
111+
// Rice-corrupted signals should all be positive
112+
if (inputSignal <= 0)
113+
{
114+
locationItr.Set(0);
115+
scaleItr.Set(0);
116+
signalItr.Set(0);
117+
118+
++inputItr;
119+
++maskItr;
120+
if (m_MeanImage)
121+
++meanItr;
122+
if (m_VarianceImage)
123+
++varItr;
124+
++locationItr;
125+
++scaleItr;
126+
++signalItr;
127+
128+
continue;
129+
}
130+
131+
// Estimation of location and scale
132+
double location = 0;
133+
double scale = m_Scale;
134+
135+
if (m_MeanImage && m_VarianceImage)
136+
{
137+
// If mean and variance images are available, use them instead of neighborhood
138+
double sigmaValue = std::sqrt(varItr.Get());
139+
double rValue = meanItr.Get() / sigmaValue;
140+
double k1Value = 0;
141+
double thetaValue = anima::FixedPointFinder(rValue, 1, k1Value);
142+
k1Value = anima::KummerFunction(-thetaValue * thetaValue / 2.0, -0.5, 1);
143+
144+
scale = sigmaValue / std::sqrt(anima::XiFunction(thetaValue * thetaValue, 1, k1Value, M_PI / 2.0));
145+
location = thetaValue * scale;
146+
}
147+
else
148+
{
149+
// Use neighbors to create samples
150+
samples.clear();
151+
weights.clear();
152+
153+
for (unsigned int i = 0; i < neighborhoodSize; ++i)
154+
{
155+
double tmpVal = static_cast<double>(inputItr.GetPixel(i, isInBounds));
156+
157+
if (isInBounds && !std::isnan(tmpVal) && std::isfinite(tmpVal))
158+
{
159+
if (maskItr.GetPixel(i) != maskItr.GetCenterPixel())
160+
continue;
161+
162+
double weight = m_NeighborWeights[i];
163+
164+
if (weight < m_Epsilon)
165+
continue;
166+
167+
samples.push_back(tmpVal);
168+
weights.push_back(weight);
169+
}
170+
}
171+
172+
if (samples.size() == 1)
173+
location = inputSignal;
174+
else
175+
anima::GetRiceParameters(samples, weights, location, scale);
176+
}
177+
178+
// Transform Rice signal in Gaussian signal
179+
double outputSignal = inputSignal;
180+
double snrValue = location / scale;
181+
182+
if (!std::isfinite(snrValue) || std::isnan(snrValue))
183+
{
184+
std::cout << snrValue << " " << location << " " << scale << std::endl;
185+
itkExceptionMacro("Estimated SNR is invalid");
186+
}
187+
188+
if (snrValue <= 0)
189+
outputSignal = 0.0;
190+
else if (snrValue <= m_Epsilon)
191+
{
192+
double unifSignal = boost::math::cdf(m_RayleighDistribution, inputSignal / scale);
193+
194+
if (unifSignal >= 1.0 - m_Alpha || unifSignal <= m_Alpha)
195+
unifSignal = boost::math::cdf(m_RayleighDistribution, snrValue);
196+
197+
outputSignal = scale * boost::math::quantile(m_NormalDistribution, unifSignal);
198+
}
199+
else if (snrValue <= 600) // if SNR if > 600 keep signal as is, else...
200+
{
201+
double unifSignal = anima::GetRiceCDF(inputSignal, location, scale);
202+
203+
if (unifSignal >= 1.0 - m_Alpha || unifSignal <= m_Alpha)
204+
unifSignal = anima::GetRiceCDF(location, location, scale);
205+
206+
outputSignal = location + scale * boost::math::quantile(m_NormalDistribution, unifSignal);
207+
}
208+
209+
m_ThreadScaleSamples[threadId].push_back(scale);
210+
211+
locationItr.Set(static_cast<OutputPixelType>(location));
212+
scaleItr.Set(static_cast<OutputPixelType>(scale));
213+
signalItr.Set(static_cast<OutputPixelType>(outputSignal));
214+
215+
++inputItr;
216+
++maskItr;
217+
if (m_MeanImage)
218+
++meanItr;
219+
if (m_VarianceImage)
220+
++varItr;
221+
++locationItr;
222+
++scaleItr;
223+
++signalItr;
224+
}
225+
}
226+
227+
template <class ComponentType,unsigned int ImageDimension>
228+
void
229+
RicianToGaussianImageFilter<ComponentType,ImageDimension>
230+
::AfterThreadedGenerateData(void)
231+
{
232+
if (m_Scale == 0)
233+
{
234+
std::vector<double> scaleSamples;
235+
for (unsigned int i = 0;i < this->GetNumberOfThreads();++i)
236+
scaleSamples.insert(scaleSamples.end(), m_ThreadScaleSamples[i].begin(), m_ThreadScaleSamples[i].end());
237+
238+
m_Scale = anima::GetMedian(scaleSamples);
239+
}
240+
241+
m_ThreadScaleSamples.clear();
242+
m_NeighborWeights.clear();
243+
244+
this->Superclass::AfterThreadedGenerateData();
245+
}
246+
247+
} // end of namespace anima

0 commit comments

Comments
 (0)