Skip to content

Commit

Permalink
Forgotten files...
Browse files Browse the repository at this point in the history
  • Loading branch information
@astamm @ocommowi
astamm authored and ocommowi committed Oct 24, 2018
1 parent 2a68dc6 commit 4243abd
Show file tree
Hide file tree
Showing 4 changed files with 712 additions and 0 deletions.
37 changes: 37 additions & 0 deletions Anima/filtering/rician_to_gaussian_noise/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,37 @@
if(BUILD_TOOLS AND USE_ARB)

project(animaRicianToGaussianNoise)

## #############################################################################
## List Sources
## #############################################################################

list_source_files(${PROJECT_NAME}
${CMAKE_CURRENT_SOURCE_DIR}
)

## #############################################################################
## add executable
## #############################################################################

add_executable(${PROJECT_NAME}
${${PROJECT_NAME}_CFILES}
)


## #############################################################################
## Link
## #############################################################################

target_link_libraries(${PROJECT_NAME}
${ITKIO_LIBRARIES}
AnimaSpecialFunctions
)

## #############################################################################
## install
## #############################################################################

set_exe_install_rules(${PROJECT_NAME})

endif()
104 changes: 104 additions & 0 deletions Anima/filtering/rician_to_gaussian_noise/animaRicianToGaussianImageFilter.h
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,104 @@
#pragma once

#include <animaMaskedImageToImageFilter.h>

#include <boost/math/distributions/normal.hpp>
#include <boost/math/distributions/rayleigh.hpp>

namespace anima
{
template <class ComponentType,unsigned int ImageDimension>
class RicianToGaussianImageFilter :
public anima::MaskedImageToImageFilter<itk::Image<ComponentType,ImageDimension>,itk::Image<ComponentType,ImageDimension> >
{
public:
/** Standard class typedefs. */
typedef RicianToGaussianImageFilter Self;

typedef itk::Image<ComponentType,ImageDimension> InputImageType;
typedef typename InputImageType::PixelType InputPixelType;

typedef itk::Image<ComponentType,ImageDimension> OutputImageType;
typedef typename OutputImageType::PixelType OutputPixelType;

typedef anima::MaskedImageToImageFilter<InputImageType, OutputImageType> Superclass;

typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;

/** Method for creation through the object factory. */
itkNewMacro(Self)

/** Run-time type information (and related methods) */
itkTypeMacro(RicianToGaussianImageFilter, MaskedImageToImageFilter)

/** Superclass typedefs. */
typedef typename Superclass::InputImageRegionType InputImageRegionType;
typedef typename Superclass::OutputImageRegionType OutputImageRegionType;
typedef typename Superclass::MaskImageType MaskImageType;
typedef typename MaskImageType::SizeType SizeType;

typedef typename InputImageType::Pointer InputPointerType;
typedef typename OutputImageType::Pointer OutputPointerType;

itkSetMacro(MaximumNumberOfIterations, unsigned int)
itkSetMacro(Epsilon, double)
itkSetMacro(Sigma, double)
itkSetMacro(MeanImage, InputPointerType)
itkSetMacro(VarianceImage, InputPointerType)

itkSetMacro(Scale, double)
itkGetConstMacro(Scale, double)

itkSetMacro(Alpha, double)
itkGetConstMacro(Alpha, double)

OutputPointerType GetLocationImage() {return this->GetOutput(0);}
OutputPointerType GetScaleImage() {return this->GetOutput(1);}
OutputPointerType GetGaussianImage() {return this->GetOutput(2);}

protected:
RicianToGaussianImageFilter()
: Superclass()
{
m_MaximumNumberOfIterations = 100;
m_Epsilon = 1.0e-8;
m_Sigma = 1.0;
m_Scale = 0.0;
m_Alpha = 0.005;
m_ThreadScaleSamples.clear();
m_NeighborWeights.clear();
m_Radius.Fill(0);
m_MeanImage = NULL;
m_VarianceImage = NULL;

this->SetNumberOfRequiredOutputs(3);
this->SetNthOutput(0, this->MakeOutput(0));
this->SetNthOutput(1, this->MakeOutput(1));
this->SetNthOutput(2, this->MakeOutput(2));
}

virtual ~RicianToGaussianImageFilter() {}

void BeforeThreadedGenerateData(void) ITK_OVERRIDE;
void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
itk::ThreadIdType threadId) ITK_OVERRIDE;
void AfterThreadedGenerateData(void) ITK_OVERRIDE;

private:
ITK_DISALLOW_COPY_AND_ASSIGN(RicianToGaussianImageFilter);

unsigned int m_MaximumNumberOfIterations;
double m_Epsilon, m_Sigma, m_Scale, m_Alpha;
std::vector<std::vector<double> > m_ThreadScaleSamples;
SizeType m_Radius;
std::vector<double> m_NeighborWeights;
InputPointerType m_MeanImage, m_VarianceImage;

boost::math::normal_distribution<> m_NormalDistribution;
boost::math::rayleigh_distribution<> m_RayleighDistribution;
};

} // end of namespace anima

#include "animaRicianToGaussianImageFilter.hxx"
247 changes: 247 additions & 0 deletions Anima/filtering/rician_to_gaussian_noise/animaRicianToGaussianImageFilter.hxx
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,247 @@
#pragma once
#include "animaRicianToGaussianImageFilter.h"
#include <animaChiDistribution.h>
#include <animaVectorOperations.h>

#include <itkConstNeighborhoodIterator.h>
#include <itkImageRegionConstIterator.h>
#include <itkImageRegionIterator.h>
#include <itkGaussianOperator.h>

namespace anima
{

template <class ComponentType,unsigned int ImageDimension>
void
RicianToGaussianImageFilter<ComponentType,ImageDimension>
::BeforeThreadedGenerateData(void)
{
// Compute spatial weights of neighbors beforehand
typename InputImageType::SpacingType spacing = this->GetInput()->GetSpacing();
typedef itk::GaussianOperator<double,ImageDimension> GaussianOperatorType;
std::vector<GaussianOperatorType> gaussianKernels(ImageDimension);

for (unsigned int i = 0;i < ImageDimension;++i)
{
unsigned int reverse_i = ImageDimension - i - 1;
double stddev = m_Sigma / spacing[i];
gaussianKernels[reverse_i].SetDirection(i);
gaussianKernels[reverse_i].SetVariance(stddev * stddev);
gaussianKernels[reverse_i].SetMaximumError(1e-3);
gaussianKernels[reverse_i].CreateDirectional();
gaussianKernels[reverse_i].ScaleCoefficients(1.0e4);
m_Radius[i] = gaussianKernels[reverse_i].GetRadius(i);
}

m_NeighborWeights.clear();
for (unsigned int i = 0;i < gaussianKernels[0].Size();++i)
{
for (unsigned int j = 0;j < gaussianKernels[1].Size();++j)
{
for (unsigned int k = 0;k < gaussianKernels[2].Size();++k)
{
double weight = gaussianKernels[0][i] * gaussianKernels[1][j] * gaussianKernels[2][k];
m_NeighborWeights.push_back(weight);
}
}
}

// Initialize thread containers for global scale estimation
unsigned int numThreads = this->GetNumberOfThreads();
m_ThreadScaleSamples.resize(numThreads);
this->Superclass::BeforeThreadedGenerateData();
}

template <class ComponentType,unsigned int ImageDimension>
void
RicianToGaussianImageFilter<ComponentType,ImageDimension>
::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
itk::ThreadIdType threadId)
{
typedef itk::ConstNeighborhoodIterator<InputImageType> InputIteratorType;
typedef itk::ConstNeighborhoodIterator<MaskImageType> MaskIteratorType;
typedef itk::ImageRegionConstIterator<InputImageType> InputSimpleIteratorType;
typedef itk::ImageRegionIterator<OutputImageType> OutputIteratorType;

InputIteratorType inputItr(m_Radius, this->GetInput(), outputRegionForThread);
unsigned int neighborhoodSize = inputItr.Size();

MaskIteratorType maskItr(m_Radius, this->GetComputationMask(), outputRegionForThread);

InputSimpleIteratorType meanItr, varItr;
if (m_MeanImage)
meanItr = InputSimpleIteratorType(m_MeanImage, outputRegionForThread);
if (m_VarianceImage)
varItr = InputSimpleIteratorType(m_VarianceImage, outputRegionForThread);

OutputIteratorType locationItr(this->GetOutput(0), outputRegionForThread);
OutputIteratorType scaleItr(this->GetOutput(1), outputRegionForThread);
OutputIteratorType signalItr(this->GetOutput(2), outputRegionForThread);

std::vector<double> samples, weights;
bool isInBounds;
typename InputImageType::IndexType currentIndex, neighborIndex;
typename InputImageType::PointType currentPoint, neighborPoint;

while (!maskItr.IsAtEnd())
{
// Discard voxels outside of the brain
if (maskItr.GetCenterPixel() == 0)
{
locationItr.Set(0);
scaleItr.Set(0);
signalItr.Set(0);

++inputItr;
++maskItr;
if (m_MeanImage)
++meanItr;
if (m_VarianceImage)
++varItr;
++locationItr;
++scaleItr;
++signalItr;

continue;
}

// Get signal at current central voxel
double inputSignal = inputItr.GetCenterPixel();

// Rice-corrupted signals should all be positive
if (inputSignal <= 0)
{
locationItr.Set(0);
scaleItr.Set(0);
signalItr.Set(0);

++inputItr;
++maskItr;
if (m_MeanImage)
++meanItr;
if (m_VarianceImage)
++varItr;
++locationItr;
++scaleItr;
++signalItr;

continue;
}

// Estimation of location and scale
double location = 0;
double scale = m_Scale;

if (m_MeanImage && m_VarianceImage)
{
// If mean and variance images are available, use them instead of neighborhood
double sigmaValue = std::sqrt(varItr.Get());
double rValue = meanItr.Get() / sigmaValue;
double k1Value = 0;
double thetaValue = anima::FixedPointFinder(rValue, 1, k1Value);
k1Value = anima::KummerFunction(-thetaValue * thetaValue / 2.0, -0.5, 1);

scale = sigmaValue / std::sqrt(anima::XiFunction(thetaValue * thetaValue, 1, k1Value, M_PI / 2.0));
location = thetaValue * scale;
}
else
{
// Use neighbors to create samples
samples.clear();
weights.clear();

for (unsigned int i = 0; i < neighborhoodSize; ++i)
{
double tmpVal = static_cast<double>(inputItr.GetPixel(i, isInBounds));

if (isInBounds && !std::isnan(tmpVal) && std::isfinite(tmpVal))
{
if (maskItr.GetPixel(i) != maskItr.GetCenterPixel())
continue;

double weight = m_NeighborWeights[i];

if (weight < m_Epsilon)
continue;

samples.push_back(tmpVal);
weights.push_back(weight);
}
}

if (samples.size() == 1)
location = inputSignal;
else
anima::GetRiceParameters(samples, weights, location, scale);
}

// Transform Rice signal in Gaussian signal
double outputSignal = inputSignal;
double snrValue = location / scale;

if (!std::isfinite(snrValue) || std::isnan(snrValue))
{
std::cout << snrValue << " " << location << " " << scale << std::endl;
itkExceptionMacro("Estimated SNR is invalid");
}

if (snrValue <= 0)
outputSignal = 0.0;
else if (snrValue <= m_Epsilon)
{
double unifSignal = boost::math::cdf(m_RayleighDistribution, inputSignal / scale);

if (unifSignal >= 1.0 - m_Alpha || unifSignal <= m_Alpha)
unifSignal = boost::math::cdf(m_RayleighDistribution, snrValue);

outputSignal = scale * boost::math::quantile(m_NormalDistribution, unifSignal);
}
else if (snrValue <= 600) // if SNR if > 600 keep signal as is, else...
{
double unifSignal = anima::GetRiceCDF(inputSignal, location, scale);

if (unifSignal >= 1.0 - m_Alpha || unifSignal <= m_Alpha)
unifSignal = anima::GetRiceCDF(location, location, scale);

outputSignal = location + scale * boost::math::quantile(m_NormalDistribution, unifSignal);
}

m_ThreadScaleSamples[threadId].push_back(scale);

locationItr.Set(static_cast<OutputPixelType>(location));
scaleItr.Set(static_cast<OutputPixelType>(scale));
signalItr.Set(static_cast<OutputPixelType>(outputSignal));

++inputItr;
++maskItr;
if (m_MeanImage)
++meanItr;
if (m_VarianceImage)
++varItr;
++locationItr;
++scaleItr;
++signalItr;
}
}

template <class ComponentType,unsigned int ImageDimension>
void
RicianToGaussianImageFilter<ComponentType,ImageDimension>
::AfterThreadedGenerateData(void)
{
if (m_Scale == 0)
{
std::vector<double> scaleSamples;
for (unsigned int i = 0;i < this->GetNumberOfThreads();++i)
scaleSamples.insert(scaleSamples.end(), m_ThreadScaleSamples[i].begin(), m_ThreadScaleSamples[i].end());

m_Scale = anima::GetMedian(scaleSamples);
}

m_ThreadScaleSamples.clear();
m_NeighborWeights.clear();

this->Superclass::AfterThreadedGenerateData();
}

} // end of namespace anima
Loading

0 comments on commit 4243abd

Please sign in to comment.