5 #include <itkSymmetricEigenAnalysis.h> 6 #include <itkImageRegionConstIterator.h> 11 template <
class ImageScalarType,
class DataImageType>
12 NLMeansVectorPatchSearcher <ImageScalarType, DataImageType>
15 m_BetaParameter = 1.0;
17 m_MeanThreshold = 0.95;
18 m_VarianceThreshold = 0.5;
21 template <
class ImageScalarType,
class DataImageType>
26 unsigned int ndim = this->GetInputImage()->GetVectorLength();
27 if (m_RefPatchMean.GetSize() != ndim)
29 m_RefPatchMean.SetSize(ndim);
30 m_RefPatchCovariance.set_size(ndim,ndim);
31 m_LogRefPatchCovariance.set_size(ndim,ndim);
35 m_RefPatchMean,m_RefPatchCovariance);
38 vnl_diag_matrix <double> eVals(ndim);
40 itk::SymmetricEigenAnalysis <vnl_matrix <double>, vnl_diag_matrix <double>, vnl_matrix <double> > EigenAnalysis(ndim);
41 EigenAnalysis.ComputeEigenValuesAndVectors(m_RefPatchCovariance, eVals, eVec);
43 for (
unsigned int i = 0;i < ndim;++i)
44 eVals[i] = std::log(eVals[i]);
46 m_LogRefPatchCovariance = eVec.transpose() * eVals * eVec;
49 for (
unsigned int i = 0;i < DataImageType::ImageDimension;++i)
51 regionLocalVariance.SetIndex(i,std::max(0,(
int)refIndex[i] - 2));
52 regionLocalVariance.SetSize(i,std::min((
unsigned int)(this->GetInputImage()->GetLargestPossibleRegion().GetSize()[i] - 1),(
unsigned int)(refIndex[i] + 2)) - regionLocalVariance.GetIndex(i) + 1);
57 vnl_matrix_inverse <double> covInverse(m_NoiseCovariance);
58 m_NoiseSigma = covInverse.inverse();
61 template <
class ImageScalarType,
class DataImageType>
66 unsigned int ndim = this->GetInputImage()->GetVectorLength();
67 if (m_MovingPatchMean.GetSize() != ndim)
69 m_MovingPatchMean.SetSize(ndim);
70 m_MovingPatchCovariance.set_size(ndim,ndim);
74 m_MovingPatchMean,m_MovingPatchCovariance);
77 template <
class ImageScalarType,
class DataImageType>
80 ::TestPatchConformity(
unsigned int index,
const IndexType &refIndex,
const IndexType &movingIndex)
82 double covStdValue = m_DatabaseCovarianceDistanceStd->GetPixel(refIndex);
83 double covMeanValue = m_DatabaseCovarianceDistanceAverage->GetPixel(refIndex);
84 double meanStdValue = m_DatabaseMeanDistanceStd->GetPixel(refIndex);
85 double meanMeanValue = m_DatabaseMeanDistanceAverage->GetPixel(refIndex);
89 if (varTest > covMeanValue + m_VarianceThreshold * covStdValue)
93 m_MovingPatchNumElements,m_RefPatchCovariance,m_MovingPatchCovariance);
95 if (meanTestValue > meanMeanValue + m_MeanThreshold * meanStdValue)
101 template <
class ImageScalarType,
class DataImageType>
106 typedef itk::ImageRegionConstIterator <RefImageType> InIteratorType;
108 InIteratorType tmpIt (this->GetInputImage(), refPatch);
109 InIteratorType tmpMovingIt (this->GetComparisonImage(index), movingPatch);
113 double weightValue = 0.0;
114 unsigned int numVoxels = 0;
115 unsigned int ndim = this->GetInputImage()->GetVectorLength();
117 while (!tmpIt.IsAtEnd())
119 tmpDiffValue = tmpIt.Get() - tmpMovingIt.Get();
120 for (
unsigned int i = 0;i < ndim;++i)
121 for (
unsigned int j = i;j < ndim;++j)
124 weightValue += 2.0 * m_NoiseSigma(i,j) * tmpDiffValue[i] * tmpDiffValue[j];
126 weightValue += m_NoiseSigma(i,j) * tmpDiffValue[i] * tmpDiffValue[j];
134 weightValue = std::exp(- weightValue / (2.0 * ndim * m_BetaParameter * numVoxels));
Superclass::IndexType IndexType
Superclass::ImageRegionType ImageRegionType
unsigned int computePatchMeanAndCovariance(const itk::VectorImage< T1, Dimension > *inputImage, const itk::ImageRegion< Dimension > &patchRegion, itk::VariableLengthVector< T2 > &patchMean, vnl_matrix< T2 > &patchCov)
Computes the average and covariance matrix from a patch of a vector image.
double VectorCovarianceTest(vnl_matrix< T > &logRefPatchCov, vnl_matrix< T > &movingPatchCov)
Test if covariance matrices are different (returns distance)
double VectorMeansTest(itk::VariableLengthVector< T > &refPatchMean, itk::VariableLengthVector< T > &movingPatchMean, const unsigned int &refPatchNumElts, const unsigned int &movingPatchNumElts, vnl_matrix< T > &refPatchCov, vnl_matrix< T > &movingPatchCov)
Test if vector means are different (returns distance)
RefImageType::PixelType VectorType
void computeAverageLocalCovariance(vnl_matrix< T2 > &resVariance, itk::VectorImage< T1, Dimension > *inputImage, itk::Image< unsigned char, Dimension > *maskImage, const itk::ImageRegion< Dimension > &averagingRegion, int localNeighborhood)
Noise estimation for a patch of a vector image.
vnl_matrix< double > CovarianceType