diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp
index 03bd6e601e..f55900fafc 100644
--- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp
+++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp
@@ -1,2226 +1,2240 @@
/*===================================================================
The Medical Imaging Interaction Toolkit (MITK)
Copyright (c) German Cancer Research Center,
Division of Medical and Biological Informatics.
All rights reserved.
This software is distributed WITHOUT ANY WARRANTY; without
even the implied warranty of MERCHANTABILITY or FITNESS FOR
A PARTICULAR PURPOSE.
See LICENSE.txt or http://www.mitk.org for details.
===================================================================*/
#include "mitkImageStatisticsCalculator.h"
#include "mitkImageAccessByItk.h"
#include "mitkImageCast.h"
#include "mitkExtractImageFilter.h"
#include "mitkImageTimeSelector.h"
#include "mitkITKImageImport.h"
#include <mitkExtendedStatisticsImageFilter.h>
#include <mitkExtendedLabelStatisticsImageFilter.h>
#include <itkScalarImageToHistogramGenerator.h>
#include <itkChangeInformationImageFilter.h>
#include <itkExtractImageFilter.h>
#include <itkMinimumMaximumImageCalculator.h>
#include <itkStatisticsImageFilter.h>
#include <itkLabelStatisticsImageFilter.h>
#include <itkMaskImageFilter.h>
#include <itkImageFileWriter.h>
#include <itkRescaleIntensityImageFilter.h>
#include <itkMaskImageFilter.h>
#include <itkImageRegionConstIterator.h>
#include <itkImageRegionIterator.h>
#include <itkCastImageFilter.h>
#include <itkVTKImageImport.h>
#include <itkVTKImageExport.h>
#include <itkImageDuplicator.h>
#include <vtkPoints.h>
#include <vtkImageStencil.h>
#include <vtkImageImport.h>
#include <vtkImageExport.h>
#include <vtkLassoStencilSource.h>
#include <itkFFTConvolutionImageFilter.h>
#include <itkConstantBoundaryCondition.h>
#include <itkImageDuplicator.h>
#include <itkContinuousIndex.h>
#include <itkNumericTraits.h>
#include <list>
#include <exception>
#include "itkImage.h"
//#define DEBUG_HOTSPOTSEARCH
#define _USE_MATH_DEFINES
#include <math.h>
#include "vtkLassoStencilSource.h"
namespace mitk
{
ImageStatisticsCalculator::ImageStatisticsCalculator()
: m_MaskingMode( MASKING_MODE_NONE ),
m_MaskingModeChanged( false ),
m_IgnorePixelValue(0.0),
m_DoIgnorePixelValue(false),
m_IgnorePixelValueChanged(false),
m_PlanarFigureAxis (0),
m_PlanarFigureSlice (0),
m_PlanarFigureCoordinate0 (0),
m_PlanarFigureCoordinate1 (0),
m_HistogramBinSize(1.0),
m_UseDefaultBinSize(true),
m_UseBinSizeBasedOnVOIRegion(false),
m_HotspotRadiusInMM(6.2035049089940), // radius of a 1cm3 sphere in mm
m_CalculateHotspot(false),
m_HotspotRadiusInMMChanged(false),
m_HotspotMustBeCompletelyInsideImage(true)
{
m_EmptyHistogram = HistogramType::New();
m_EmptyHistogram->SetMeasurementVectorSize(1);
HistogramType::SizeType histogramSize(1);
histogramSize.Fill( 256 );
m_EmptyHistogram->Initialize( histogramSize );
m_EmptyStatistics.Reset();
}
ImageStatisticsCalculator::~ImageStatisticsCalculator()
{
}
void ImageStatisticsCalculator::SetUseDefaultBinSize(bool useDefault)
{
m_UseDefaultBinSize = useDefault;
}
ImageStatisticsCalculator::Statistics::Statistics(bool withHotspotStatistics)
:m_HotspotStatistics(withHotspotStatistics ? new Statistics(false) : nullptr)
{
Reset();
}
ImageStatisticsCalculator::Statistics::Statistics(const Statistics& other)
:m_HotspotStatistics( nullptr)
{
this->SetLabel( other.GetLabel() );
this->SetN( other.GetN() );
this->SetMin( other.GetMin() );
this->SetMax( other.GetMax() );
this->SetMedian( other.GetMedian() );
this->SetMean( other.GetMean() );
this->SetVariance( other.GetVariance() );
this->SetKurtosis( other.GetKurtosis() );
this->SetSkewness( other.GetSkewness() );
this->SetUniformity( other.GetUniformity() );
this->SetEntropy( other.GetEntropy() );
this->SetUPP( other.GetUPP() );
this->SetMPP( other.GetMPP() );
this->SetSigma( other.GetSigma() );
this->SetRMS( other.GetRMS() );
this->SetMaxIndex( other.GetMaxIndex() );
this->SetMinIndex( other.GetMinIndex() );
this->SetHotspotIndex( other.GetHotspotIndex() );
if (other.m_HotspotStatistics)
{
this->m_HotspotStatistics = new Statistics(false);
*this->m_HotspotStatistics = *other.m_HotspotStatistics;
}
}
bool ImageStatisticsCalculator::Statistics::HasHotspotStatistics() const
{
return m_HotspotStatistics != nullptr;
}
void ImageStatisticsCalculator::Statistics::SetHasHotspotStatistics(bool hasHotspotStatistics)
{
m_HasHotspotStatistics = hasHotspotStatistics;
}
ImageStatisticsCalculator::Statistics::~Statistics()
{
delete m_HotspotStatistics;
}
double ImageStatisticsCalculator::Statistics::GetVariance() const
{
return this->Variance;
}
void ImageStatisticsCalculator::Statistics::SetVariance( const double value )
{
if( this->Variance != value )
{
if( value < 0.0 )
{
this->Variance = 0.0; // if given value is negative set variance to 0.0
}
else
{
this->Variance = value;
}
}
}
double ImageStatisticsCalculator::Statistics::GetSigma() const
{
return this->Sigma;
}
void ImageStatisticsCalculator::Statistics::SetSigma( const double value )
{
if( this->Sigma != value )
{
// for some compiler the value != value works to check for NaN but not for all
// but we can always be sure that the standard deviation is a positive value
if( value != value || value < 0.0 )
{
// if standard deviation is NaN we just assume 0.0
this->Sigma = 0.0;
}
else
{
this->Sigma = value;
}
}
}
void ImageStatisticsCalculator::Statistics::Reset(unsigned int dimension)
{
SetLabel(0);
SetN( 0 );
SetMin( 0.0 );
SetMax( 0.0 );
SetMedian( 0.0 );
SetVariance( 0.0 );
SetMean( 0.0 );
SetSigma( 0.0 );
SetRMS( 0.0 );
SetSkewness( 0.0 );
SetKurtosis( 0.0 );
SetUniformity( 0.0 );
SetEntropy( 0.0 );
SetMPP( 0.0 );
SetUPP( 0.0 );
vnl_vector<int> zero;
zero.set_size(dimension);
for(unsigned int i = 0; i < dimension; ++i)
{
zero[i] = 0;
}
SetMaxIndex(zero);
SetMinIndex(zero);
SetHotspotIndex(zero);
if (m_HotspotStatistics != nullptr)
{
m_HotspotStatistics->Reset(dimension);
}
}
const ImageStatisticsCalculator::Statistics&
ImageStatisticsCalculator::Statistics::GetHotspotStatistics() const
{
if (m_HotspotStatistics)
{
return *m_HotspotStatistics;
}
else
{
throw std::logic_error("Object has no hostspot statistics, see HasHotspotStatistics()");
}
}
ImageStatisticsCalculator::Statistics&
ImageStatisticsCalculator::Statistics::GetHotspotStatistics()
{
if (m_HotspotStatistics)
{
return *m_HotspotStatistics;
}
else
{
throw std::logic_error("Object has no hostspot statistics, see HasHotspotStatistics()");
}
}
ImageStatisticsCalculator::Statistics&
ImageStatisticsCalculator::Statistics::operator=(ImageStatisticsCalculator::Statistics const& other)
{
if (this == &other)
return *this;
this->SetLabel( other.GetLabel() );
this->SetN( other.GetN() );
this->SetMin( other.GetMin() );
this->SetMax( other.GetMax() );
this->SetMean( other.GetMean() );
this->SetMedian( other.GetMedian() );
this->SetVariance( other.GetVariance() );
this->SetSigma( other.GetSigma() );
this->SetRMS( other.GetRMS() );
this->SetMinIndex( other.GetMinIndex() );
this->SetMaxIndex( other.GetMaxIndex() );
this->SetHotspotIndex( other.GetHotspotIndex() );
this->SetSkewness( other.GetSkewness() );
this->SetKurtosis( other.GetKurtosis() );
this->SetUniformity( other.GetUniformity() );
this->SetEntropy( other.GetEntropy() );
this->SetUPP( other.GetUPP() );
this->SetMPP( other.GetMPP() );
delete this->m_HotspotStatistics;
this->m_HotspotStatistics = nullptr;
if (other.m_HotspotStatistics)
{
this->m_HotspotStatistics = new Statistics(false);
*this->m_HotspotStatistics = *other.m_HotspotStatistics;
}
return *this;
}
void ImageStatisticsCalculator::SetImage( const mitk::Image *image )
{
if ( m_Image != image )
{
m_Image = image;
this->Modified();
unsigned int numberOfTimeSteps = image->GetTimeSteps();
// Initialize vectors to time-size of this image
m_ImageHistogramVector.resize( numberOfTimeSteps );
m_MaskedImageHistogramVector.resize( numberOfTimeSteps );
m_PlanarFigureHistogramVector.resize( numberOfTimeSteps );
m_ImageStatisticsVector.resize( numberOfTimeSteps );
m_MaskedImageStatisticsVector.resize( numberOfTimeSteps );
m_PlanarFigureStatisticsVector.resize( numberOfTimeSteps );
m_ImageStatisticsTimeStampVector.resize( numberOfTimeSteps );
m_MaskedImageStatisticsTimeStampVector.resize( numberOfTimeSteps );
m_PlanarFigureStatisticsTimeStampVector.resize( numberOfTimeSteps );
m_ImageStatisticsCalculationTriggerVector.resize( numberOfTimeSteps );
m_MaskedImageStatisticsCalculationTriggerVector.resize( numberOfTimeSteps );
m_PlanarFigureStatisticsCalculationTriggerVector.resize( numberOfTimeSteps );
for ( unsigned int t = 0; t < image->GetTimeSteps(); ++t )
{
m_ImageStatisticsTimeStampVector[t].Modified();
m_ImageStatisticsCalculationTriggerVector[t] = true;
}
}
}
void ImageStatisticsCalculator::SetImageMask( const mitk::Image *imageMask )
{
if ( m_Image.IsNull() )
{
itkExceptionMacro( << "Image needs to be set first!" );
}
if ( m_ImageMask != imageMask )
{
m_ImageMask = imageMask;
this->Modified();
for ( unsigned int t = 0; t < m_Image->GetTimeSteps(); ++t )
{
m_MaskedImageStatisticsTimeStampVector[t].Modified();
m_MaskedImageStatisticsCalculationTriggerVector[t] = true;
}
}
}
void ImageStatisticsCalculator::SetPlanarFigure( mitk::PlanarFigure *planarFigure )
{
if ( m_Image.IsNull() )
{
itkExceptionMacro( << "Image needs to be set first!" );
}
if ( m_PlanarFigure != planarFigure )
{
m_PlanarFigure = planarFigure;
this->Modified();
for ( unsigned int t = 0; t < m_Image->GetTimeSteps(); ++t )
{
m_PlanarFigureStatisticsTimeStampVector[t].Modified();
m_PlanarFigureStatisticsCalculationTriggerVector[t] = true;
}
}
}
void ImageStatisticsCalculator::SetMaskingMode( unsigned int mode )
{
if ( m_MaskingMode != mode )
{
m_MaskingMode = mode;
m_MaskingModeChanged = true;
this->Modified();
}
}
void ImageStatisticsCalculator::SetMaskingModeToNone()
{
if ( m_MaskingMode != MASKING_MODE_NONE )
{
m_MaskingMode = MASKING_MODE_NONE;
m_MaskingModeChanged = true;
this->Modified();
}
}
void ImageStatisticsCalculator::SetMaskingModeToImage()
{
if ( m_MaskingMode != MASKING_MODE_IMAGE )
{
m_MaskingMode = MASKING_MODE_IMAGE;
m_MaskingModeChanged = true;
this->Modified();
}
}
void ImageStatisticsCalculator::SetMaskingModeToPlanarFigure()
{
if ( m_MaskingMode != MASKING_MODE_PLANARFIGURE )
{
m_MaskingMode = MASKING_MODE_PLANARFIGURE;
m_MaskingModeChanged = true;
this->Modified();
}
}
void ImageStatisticsCalculator::SetIgnorePixelValue(double value)
{
if ( m_IgnorePixelValue != value )
{
m_IgnorePixelValue = value;
if(m_DoIgnorePixelValue)
{
m_IgnorePixelValueChanged = true;
}
this->Modified();
}
}
double ImageStatisticsCalculator::GetIgnorePixelValue()
{
return m_IgnorePixelValue;
}
void ImageStatisticsCalculator::SetDoIgnorePixelValue(bool value)
{
if ( m_DoIgnorePixelValue != value )
{
m_DoIgnorePixelValue = value;
m_IgnorePixelValueChanged = true;
this->Modified();
}
}
bool ImageStatisticsCalculator::GetDoIgnorePixelValue()
{
return m_DoIgnorePixelValue;
}
void ImageStatisticsCalculator::SetHistogramBinSize(double size)
{
this->m_HistogramBinSize = size;
}
double ImageStatisticsCalculator::GetHistogramBinSize()
{
return this->m_HistogramBinSize;
}
void ImageStatisticsCalculator::SetHotspotRadiusInMM(double value)
{
if ( m_HotspotRadiusInMM != value )
{
m_HotspotRadiusInMM = value;
if(m_CalculateHotspot)
{
m_HotspotRadiusInMMChanged = true;
//MITK_INFO <<"Hotspot radius changed, new convolution required";
}
this->Modified();
}
}
double ImageStatisticsCalculator::GetHotspotRadiusInMM()
{
return m_HotspotRadiusInMM;
}
void ImageStatisticsCalculator::SetCalculateHotspot(bool on)
{
if ( m_CalculateHotspot != on )
{
m_CalculateHotspot = on;
m_HotspotRadiusInMMChanged = true;
//MITK_INFO <<"Hotspot calculation changed, new convolution required";
this->Modified();
}
}
bool ImageStatisticsCalculator::IsHotspotCalculated()
{
return m_CalculateHotspot;
}
void ImageStatisticsCalculator::SetHotspotMustBeCompletlyInsideImage(bool hotspotMustBeCompletelyInsideImage, bool warn)
{
m_HotspotMustBeCompletelyInsideImage = hotspotMustBeCompletelyInsideImage;
if (!m_HotspotMustBeCompletelyInsideImage && warn)
{
MITK_WARN << "Hotspot calculation will extrapolate pixels at image borders. Be aware of the consequences for the hotspot location.";
}
}
bool ImageStatisticsCalculator::GetHotspotMustBeCompletlyInsideImage() const
{
return m_HotspotMustBeCompletelyInsideImage;
}
/* Implementation of the min max values for setting the range of the histogram */
template < typename TPixel, unsigned int VImageDimension >
void ImageStatisticsCalculator::GetMinAndMaxValue( double &min,
double &max, int &counter, double &sigma,
const itk::Image< TPixel, VImageDimension > *InputImage,
itk::Image< unsigned short, VImageDimension > *MaskedImage )
{
typedef itk::Image< unsigned short, VImageDimension > MaskImageType;
typedef itk::Image< TPixel, VImageDimension > ImageType;
typedef itk::ImageRegionConstIteratorWithIndex<ImageType> Imageie;
typedef itk::ImageRegionConstIteratorWithIndex<MaskImageType> Imageie2;
Imageie2 labelIterator2( MaskedImage, MaskedImage->GetRequestedRegion() );
Imageie labelIterator3( InputImage, InputImage->GetRequestedRegion() );
max = 0;
min = 0;
counter = 0;
sigma = 0;
double SumOfSquares = 0;
double sumSquared = 0;
double actualPielValue = 0;
int counterOfPixelsInROI = 0;
for( labelIterator2.GoToBegin(); !labelIterator2.IsAtEnd(); ++labelIterator2, ++labelIterator3)
{
if( labelIterator2.Value()== 1.0)
{
counter++;
counterOfPixelsInROI++;
actualPielValue = labelIterator3.Value();
sumSquared = sumSquared + actualPielValue;
SumOfSquares = SumOfSquares + std::pow(actualPielValue,2);
if(counterOfPixelsInROI == 1)
{
max = actualPielValue;
min = actualPielValue;
}
if(actualPielValue >= max)
{
max = actualPielValue;
}
else if(actualPielValue <= min)
{
min = actualPielValue;
}
}
}
if (counter > 1)
{
sigma = ( SumOfSquares - std::pow( sumSquared, 2) / counter ) / ( counter-1 );
}
else
{
sigma = 0;
}
}
bool ImageStatisticsCalculator::ComputeStatistics( unsigned int timeStep )
{
if (m_Image.IsNull() )
{
mitkThrow() << "Image not set!";
}
if (!m_Image->IsInitialized())
{
mitkThrow() << "Image not initialized!";
}
if ( m_Image->GetReferenceCount() == 1 )
{
// Image no longer valid; we are the only ones to still hold a reference on it
return false;
}
if ( timeStep >= m_Image->GetTimeSteps() )
{
throw std::runtime_error( "Error: invalid time step!" );
}
// If a mask was set but we are the only ones to still hold a reference on
// it, delete it.
if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() == 1) )
{
m_ImageMask = nullptr;
}
// Check if statistics is already up-to-date
unsigned long imageMTime = m_ImageStatisticsTimeStampVector[timeStep].GetMTime();
unsigned long maskedImageMTime = m_MaskedImageStatisticsTimeStampVector[timeStep].GetMTime();
unsigned long planarFigureMTime = m_PlanarFigureStatisticsTimeStampVector[timeStep].GetMTime();
bool imageStatisticsCalculationTrigger = m_ImageStatisticsCalculationTriggerVector[timeStep];
bool maskedImageStatisticsCalculationTrigger = m_MaskedImageStatisticsCalculationTriggerVector[timeStep];
bool planarFigureStatisticsCalculationTrigger = m_PlanarFigureStatisticsCalculationTriggerVector[timeStep];
if ( !m_IgnorePixelValueChanged
&& !m_HotspotRadiusInMMChanged
&& ((m_MaskingMode != MASKING_MODE_NONE) || (imageMTime > m_Image->GetMTime() && !imageStatisticsCalculationTrigger))
&& ((m_MaskingMode != MASKING_MODE_IMAGE) || (maskedImageMTime > m_ImageMask->GetMTime() && !maskedImageStatisticsCalculationTrigger))
&& ((m_MaskingMode != MASKING_MODE_PLANARFIGURE) || (planarFigureMTime > m_PlanarFigure->GetMTime() && !planarFigureStatisticsCalculationTrigger)) )
{
// Statistics is up to date!
if ( m_MaskingModeChanged )
{
m_MaskingModeChanged = false;
}
else
{
return false;
}
}
// Reset state changed flag
m_MaskingModeChanged = false;
m_IgnorePixelValueChanged = false;
// Depending on masking mode, extract and/or generate the required image
// and mask data from the user input
this->ExtractImageAndMask( timeStep );
StatisticsContainer *statisticsContainer;
HistogramContainer *histogramContainer;
switch ( m_MaskingMode )
{
case MASKING_MODE_NONE:
default:
if(!m_DoIgnorePixelValue)
{
statisticsContainer = &m_ImageStatisticsVector[timeStep];
histogramContainer = &m_ImageHistogramVector[timeStep];
m_ImageStatisticsTimeStampVector[timeStep].Modified();
m_ImageStatisticsCalculationTriggerVector[timeStep] = false;
}
else
{
statisticsContainer = &m_MaskedImageStatisticsVector[timeStep];
histogramContainer = &m_MaskedImageHistogramVector[timeStep];
m_MaskedImageStatisticsTimeStampVector[timeStep].Modified();
m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false;
}
break;
case MASKING_MODE_IMAGE:
statisticsContainer = &m_MaskedImageStatisticsVector[timeStep];
histogramContainer = &m_MaskedImageHistogramVector[timeStep];
m_MaskedImageStatisticsTimeStampVector[timeStep].Modified();
m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false;
break;
case MASKING_MODE_PLANARFIGURE:
statisticsContainer = &m_PlanarFigureStatisticsVector[timeStep];
histogramContainer = &m_PlanarFigureHistogramVector[timeStep];
m_PlanarFigureStatisticsTimeStampVector[timeStep].Modified();
m_PlanarFigureStatisticsCalculationTriggerVector[timeStep] = false;
break;
}
// Calculate statistics and histogram(s)
if ( m_InternalImage->GetDimension() == 3 )
{
if ( m_MaskingMode == MASKING_MODE_NONE && !m_DoIgnorePixelValue )
{
AccessFixedDimensionByItk_2(
m_InternalImage,
InternalCalculateStatisticsUnmasked,
3,
statisticsContainer,
histogramContainer );
}
else
{
AccessFixedDimensionByItk_3(
m_InternalImage,
InternalCalculateStatisticsMasked,
3,
m_InternalImageMask3D.GetPointer(),
statisticsContainer,
histogramContainer );
}
}
else if ( m_InternalImage->GetDimension() == 2 )
{
if ( m_MaskingMode == MASKING_MODE_NONE && !m_DoIgnorePixelValue )
{
AccessFixedDimensionByItk_2(
m_InternalImage,
InternalCalculateStatisticsUnmasked,
2,
statisticsContainer,
histogramContainer );
}
else
{
AccessFixedDimensionByItk_3(
m_InternalImage,
InternalCalculateStatisticsMasked,
2,
m_InternalImageMask2D.GetPointer(),
statisticsContainer,
histogramContainer );
}
}
else
{
MITK_ERROR << "ImageStatistics: Image dimension not supported!";
}
// Release unused image smart pointers to free memory
m_InternalImage = mitk::Image::ConstPointer();
m_InternalImageMask3D = MaskImage3DType::Pointer();
m_InternalImageMask2D = MaskImage2DType::Pointer();
return true;
}
ImageStatisticsCalculator::BinFrequencyType
ImageStatisticsCalculator::GetBinsAndFreuqencyForHistograms( unsigned int timeStep , unsigned int label ) const
{
const HistogramType *binsAndFrequencyToCalculate = this->GetHistogram(0);
// ToDo: map should be created on stack not on heap
std::map<int, double> returnedHistogramMap;
unsigned int size = binsAndFrequencyToCalculate->Size();
for( unsigned int bin=0; bin < size; ++bin )
{
double frequency = binsAndFrequencyToCalculate->GetFrequency( bin, 0 );
//if( frequency > mitk::eps )
{
returnedHistogramMap.insert( std::pair<int, double>(binsAndFrequencyToCalculate->GetMeasurement( bin, 0 ), binsAndFrequencyToCalculate->GetFrequency( bin, 0 ) ) );
}
}
return returnedHistogramMap;
}
const ImageStatisticsCalculator::HistogramType *
ImageStatisticsCalculator::GetHistogram( unsigned int timeStep, unsigned int label ) const
{
if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) )
{
return nullptr;
}
switch ( m_MaskingMode )
{
case MASKING_MODE_NONE:
default:
{
if(m_DoIgnorePixelValue)
return m_MaskedImageHistogramVector[timeStep][label];
return m_ImageHistogramVector[timeStep][label];
}
case MASKING_MODE_IMAGE:
return m_MaskedImageHistogramVector[timeStep][label];
case MASKING_MODE_PLANARFIGURE:
return m_PlanarFigureHistogramVector[timeStep][label];
}
}
const ImageStatisticsCalculator::HistogramContainer &
ImageStatisticsCalculator::GetHistogramVector( unsigned int timeStep ) const
{
if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) )
{
return m_EmptyHistogramContainer;
}
switch ( m_MaskingMode )
{
case MASKING_MODE_NONE:
default:
{
if(m_DoIgnorePixelValue)
return m_MaskedImageHistogramVector[timeStep];
return m_ImageHistogramVector[timeStep];
}
case MASKING_MODE_IMAGE:
return m_MaskedImageHistogramVector[timeStep];
case MASKING_MODE_PLANARFIGURE:
return m_PlanarFigureHistogramVector[timeStep];
}
}
const ImageStatisticsCalculator::Statistics &
ImageStatisticsCalculator::GetStatistics( unsigned int timeStep, unsigned int label ) const
{
if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) )
{
return m_EmptyStatistics;
}
switch ( m_MaskingMode )
{
case MASKING_MODE_NONE:
default:
{
if(m_DoIgnorePixelValue)
return m_MaskedImageStatisticsVector[timeStep][label];
return m_ImageStatisticsVector[timeStep][label];
}
case MASKING_MODE_IMAGE:
return m_MaskedImageStatisticsVector[timeStep][label];
case MASKING_MODE_PLANARFIGURE:
return m_PlanarFigureStatisticsVector[timeStep][label];
}
}
const ImageStatisticsCalculator::StatisticsContainer &
ImageStatisticsCalculator::GetStatisticsVector( unsigned int timeStep ) const
{
if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) )
{
return m_EmptyStatisticsContainer;
}
switch ( m_MaskingMode )
{
case MASKING_MODE_NONE:
default:
{
if(m_DoIgnorePixelValue)
return m_MaskedImageStatisticsVector[timeStep];
return m_ImageStatisticsVector[timeStep];
}
case MASKING_MODE_IMAGE:
return m_MaskedImageStatisticsVector[timeStep];
case MASKING_MODE_PLANARFIGURE:
return m_PlanarFigureStatisticsVector[timeStep];
}
}
void ImageStatisticsCalculator::ExtractImageAndMask( unsigned int timeStep )
{
if ( m_Image.IsNull() )
{
throw std::runtime_error( "Error: image empty!" );
}
if ( timeStep >= m_Image->GetTimeSteps() )
{
throw std::runtime_error( "Error: invalid time step!" );
}
ImageTimeSelector::Pointer imageTimeSelector = ImageTimeSelector::New();
imageTimeSelector->SetInput( m_Image );
imageTimeSelector->SetTimeNr( timeStep );
imageTimeSelector->UpdateLargestPossibleRegion();
mitk::Image *timeSliceImage = imageTimeSelector->GetOutput();
switch ( m_MaskingMode )
{
case MASKING_MODE_NONE:
{
m_InternalImage = timeSliceImage;
m_InternalImageMask2D = nullptr;
m_InternalImageMask3D = nullptr;
if(m_DoIgnorePixelValue)
{
if( m_InternalImage->GetDimension() == 3 )
{
if(itk::ImageIOBase::USHORT != timeSliceImage->GetPixelType().GetComponentType())
CastToItkImage( timeSliceImage, m_InternalImageMask3D );
else
CastToItkImage( timeSliceImage->Clone(), m_InternalImageMask3D );
m_InternalImageMask3D->FillBuffer(1);
}
if( m_InternalImage->GetDimension() == 2 )
{
if(itk::ImageIOBase::USHORT != timeSliceImage->GetPixelType().GetComponentType())
CastToItkImage( timeSliceImage, m_InternalImageMask2D );
else
CastToItkImage( timeSliceImage->Clone(), m_InternalImageMask2D );
m_InternalImageMask2D->FillBuffer(1);
}
}
break;
}
case MASKING_MODE_IMAGE:
{
if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() > 1) )
{
if ( timeStep >= m_ImageMask->GetTimeSteps() )
{
// Use the last mask time step in case the current time step is bigger than the total
// number of mask time steps.
// It makes more sense setting this to the last mask time step than to 0.
// For instance if you have a mask with 2 time steps and an image with 5:
// If time step 0 is selected, the mask will use time step 0.
// If time step 1 is selected, the mask will use time step 1.
// If time step 2+ is selected, the mask will use time step 1.
// If you have a mask with only one time step instead, this will always default to 0.
timeStep = m_ImageMask->GetTimeSteps() - 1;
}
ImageTimeSelector::Pointer maskedImageTimeSelector = ImageTimeSelector::New();
maskedImageTimeSelector->SetInput( m_ImageMask );
maskedImageTimeSelector->SetTimeNr( timeStep );
maskedImageTimeSelector->UpdateLargestPossibleRegion();
mitk::Image *timeSliceMaskedImage = maskedImageTimeSelector->GetOutput();
m_InternalImage = timeSliceImage;
CastToItkImage( timeSliceMaskedImage, m_InternalImageMask3D );
}
else
{
throw std::runtime_error( "Error: image mask empty!" );
}
break;
}
case MASKING_MODE_PLANARFIGURE:
{
m_InternalImageMask2D = nullptr;
if ( m_PlanarFigure.IsNull() )
{
throw std::runtime_error( "Error: planar figure empty!" );
}
if ( !m_PlanarFigure->IsClosed() )
{
throw std::runtime_error( "Masking not possible for non-closed figures" );
}
const BaseGeometry *imageGeometry = timeSliceImage->GetGeometry();
if ( imageGeometry == nullptr )
{
throw std::runtime_error( "Image geometry invalid!" );
}
const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry();
if ( planarFigurePlaneGeometry == nullptr )
{
throw std::runtime_error( "Planar-Figure not yet initialized!" );
}
const PlaneGeometry *planarFigureGeometry =
dynamic_cast< const PlaneGeometry * >( planarFigurePlaneGeometry );
if ( planarFigureGeometry == nullptr )
{
throw std::runtime_error( "Non-planar planar figures not supported!" );
}
// Find principal direction of PlanarFigure in input image
unsigned int axis;
if ( !this->GetPrincipalAxis( imageGeometry,
planarFigureGeometry->GetNormal(), axis ) )
{
throw std::runtime_error( "Non-aligned planar figures not supported!" );
}
m_PlanarFigureAxis = axis;
// Find slice number corresponding to PlanarFigure in input image
MaskImage3DType::IndexType index;
imageGeometry->WorldToIndex( planarFigureGeometry->GetOrigin(), index );
unsigned int slice = index[axis];
m_PlanarFigureSlice = slice;
// Extract slice with given position and direction from image
unsigned int dimension = timeSliceImage->GetDimension();
if (dimension != 2)
{
ExtractImageFilter::Pointer imageExtractor = ExtractImageFilter::New();
imageExtractor->SetInput( timeSliceImage );
imageExtractor->SetSliceDimension( axis );
imageExtractor->SetSliceIndex( slice );
imageExtractor->Update();
m_InternalImage = imageExtractor->GetOutput();
}
else
{
m_InternalImage = timeSliceImage;
}
// Compute mask from PlanarFigure
AccessFixedDimensionByItk_1(
m_InternalImage,
InternalCalculateMaskFromPlanarFigure,
2, axis );
}
}
if(m_DoIgnorePixelValue)
{
if ( m_InternalImage->GetDimension() == 3 )
{
AccessFixedDimensionByItk_1(
m_InternalImage,
InternalMaskIgnoredPixels,
3,
m_InternalImageMask3D.GetPointer() );
}
else if ( m_InternalImage->GetDimension() == 2 )
{
AccessFixedDimensionByItk_1(
m_InternalImage,
InternalMaskIgnoredPixels,
2,
m_InternalImageMask2D.GetPointer() );
}
}
}
bool ImageStatisticsCalculator::GetPrincipalAxis(
const BaseGeometry *geometry, Vector3D vector,
unsigned int &axis )
{
vector.Normalize();
for ( unsigned int i = 0; i < 3; ++i )
{
Vector3D axisVector = geometry->GetAxisVector( i );
axisVector.Normalize();
if ( fabs( fabs( axisVector * vector ) - 1.0) < mitk::eps )
{
axis = i;
return true;
}
}
return false;
}
unsigned int ImageStatisticsCalculator::calcNumberOfBins(mitk::ScalarType min, mitk::ScalarType max)
{
return std::ceil( ( (max - min ) / m_HistogramBinSize) );
}
template < typename TPixel, unsigned int VImageDimension >
void ImageStatisticsCalculator::InternalCalculateStatisticsUnmasked(
const itk::Image< TPixel, VImageDimension > *image,
StatisticsContainer *statisticsContainer,
HistogramContainer* histogramContainer )
{
typedef itk::Image< TPixel, VImageDimension > ImageType;
typedef typename ImageType::IndexType IndexType;
typedef itk::Statistics::ScalarImageToHistogramGenerator< ImageType >
HistogramGeneratorType;
statisticsContainer->clear();
histogramContainer->clear();
// Progress listening...
typedef itk::SimpleMemberCommand< ImageStatisticsCalculator > ITKCommandType;
ITKCommandType::Pointer progressListener;
progressListener = ITKCommandType::New();
progressListener->SetCallbackFunction( this,
&ImageStatisticsCalculator::UnmaskedStatisticsProgressUpdate );
// Issue 100 artificial progress events since ScalarIMageToHistogramGenerator
// does not (yet?) support progress reporting
this->InvokeEvent( itk::StartEvent() );
for ( unsigned int i = 0; i < 100; ++i )
{
this->UnmaskedStatisticsProgressUpdate();
}
// Calculate statistics (separate filter)
typedef itk::ExtendedStatisticsImageFilter< ImageType > StatisticsFilterType;
typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New();
statisticsFilter->SetInput( image );
statisticsFilter->SetBinSize( 100 );
+ statisticsFilter->SetCoordinateTolerance( 0.001 );
+ statisticsFilter->SetDirectionTolerance( 0.001 );
unsigned long observerTag = statisticsFilter->AddObserver( itk::ProgressEvent(), progressListener );
try
{
statisticsFilter->Update();
}
catch (const itk::ExceptionObject& e)
{
mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what();
}
catch( const std::exception& e )
{
//mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what();
}
statisticsFilter->RemoveObserver( observerTag );
this->InvokeEvent( itk::EndEvent() );
// Calculate minimum and maximum
typedef itk::MinimumMaximumImageCalculator< ImageType > MinMaxFilterType;
typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New();
minMaxFilter->SetImage( image );
unsigned long observerTag2 = minMaxFilter->AddObserver( itk::ProgressEvent(), progressListener );
minMaxFilter->Compute();
minMaxFilter->RemoveObserver( observerTag2 );
this->InvokeEvent( itk::EndEvent() );
Statistics statistics;
statistics.Reset();
statistics.SetLabel(1);
statistics.SetN(image->GetBufferedRegion().GetNumberOfPixels());
statistics.SetMin(statisticsFilter->GetMinimum());
statistics.SetMax(statisticsFilter->GetMaximum());
statistics.SetMean(statisticsFilter->GetMean());
statistics.SetMedian(statisticsFilter->GetMedian());
statistics.SetVariance(statisticsFilter->GetVariance());
statistics.SetSkewness(statisticsFilter->GetSkewness());
statistics.SetKurtosis(statisticsFilter->GetKurtosis());
statistics.SetUniformity( statisticsFilter->GetUniformity());
statistics.SetEntropy( statisticsFilter->GetEntropy());
statistics.SetUPP( statisticsFilter->GetUPP());
statistics.SetMPP( statisticsFilter->GetMPP());
statistics.SetSigma(statisticsFilter->GetSigma());
statistics.SetRMS(sqrt( statistics.GetMean() * statistics.GetMean() + statistics.GetSigma() * statistics.GetSigma() ));
statistics.GetMinIndex().set_size(image->GetImageDimension());
statistics.GetMaxIndex().set_size(image->GetImageDimension());
vnl_vector<int> tmpMaxIndex;
vnl_vector<int> tmpMinIndex;
tmpMaxIndex.set_size(image->GetImageDimension() );
tmpMinIndex.set_size(image->GetImageDimension() );
for (unsigned int i=0; i<statistics.GetMaxIndex().size(); i++)
{
tmpMaxIndex[i] = minMaxFilter->GetIndexOfMaximum()[i];
tmpMinIndex[i] = minMaxFilter->GetIndexOfMinimum()[i];
}
statistics.SetMinIndex(tmpMaxIndex);
statistics.SetMinIndex(tmpMinIndex);
if( IsHotspotCalculated() && VImageDimension == 3 )
{
typedef itk::Image< unsigned short, VImageDimension > MaskImageType;
typename MaskImageType::Pointer nullMask;
bool isHotspotDefined(false);
Statistics hotspotStatistics = this->CalculateHotspotStatistics(image, nullMask.GetPointer(), m_HotspotRadiusInMM, isHotspotDefined, 0 );
if (isHotspotDefined)
{
statistics.SetHasHotspotStatistics(true);
statistics.GetHotspotStatistics() = hotspotStatistics;
}
else
{
statistics.SetHasHotspotStatistics(false);
}
if(statistics.GetHotspotStatistics().HasHotspotStatistics() )
{
MITK_DEBUG << "Hotspot statistics available";
statistics.SetHotspotIndex(hotspotStatistics.GetHotspotIndex());
}
else
{
MITK_ERROR << "No hotspot statistics available!";
}
}
statisticsContainer->push_back( statistics );
// Calculate histogram
// calculate bin size or number of bins
unsigned int numberOfBins = 200; // default number of bins
if (m_UseDefaultBinSize)
{
m_HistogramBinSize = std::ceil( (statistics.GetMax() - statistics.GetMin() + 1)/numberOfBins );
}
else
{
numberOfBins = calcNumberOfBins(statistics.GetMin(), statistics.GetMax());
}
typename HistogramGeneratorType::Pointer histogramGenerator = HistogramGeneratorType::New();
histogramGenerator->SetInput( image );
histogramGenerator->SetMarginalScale( 100 );
histogramGenerator->SetNumberOfBins( numberOfBins );
histogramGenerator->SetHistogramMin( statistics.GetMin() );
histogramGenerator->SetHistogramMax( statistics.GetMax() );
histogramGenerator->Compute();
histogramContainer->push_back( histogramGenerator->GetOutput() );
}
template < typename TPixel, unsigned int VImageDimension >
void ImageStatisticsCalculator::InternalMaskIgnoredPixels(
const itk::Image< TPixel, VImageDimension > *image,
itk::Image< unsigned short, VImageDimension > *maskImage )
{
typedef itk::Image< TPixel, VImageDimension > ImageType;
typedef itk::Image< unsigned short, VImageDimension > MaskImageType;
itk::ImageRegionIterator<MaskImageType>
itmask(maskImage, maskImage->GetLargestPossibleRegion());
itk::ImageRegionConstIterator<ImageType>
itimage(image, image->GetLargestPossibleRegion());
itmask.GoToBegin();
itimage.GoToBegin();
while( !itmask.IsAtEnd() )
{
if(m_IgnorePixelValue == itimage.Get())
{
itmask.Set(0);
}
++itmask;
++itimage;
}
}
template < typename TPixel, unsigned int VImageDimension >
void ImageStatisticsCalculator::InternalCalculateStatisticsMasked(
const itk::Image< TPixel, VImageDimension > *image,
itk::Image< unsigned short, VImageDimension > *maskImage,
StatisticsContainer* statisticsContainer,
HistogramContainer* histogramContainer )
{
typedef itk::Image< TPixel, VImageDimension > ImageType;
typedef itk::Image< unsigned short, VImageDimension > MaskImageType;
typedef typename ImageType::IndexType IndexType;
typedef typename ImageType::PointType PointType;
typedef typename ImageType::SpacingType SpacingType;
typedef typename ImageType::Pointer ImagePointer;
typedef itk::ExtendedLabelStatisticsImageFilter< ImageType, MaskImageType > LabelStatisticsFilterType;
typedef itk::ChangeInformationImageFilter< MaskImageType > ChangeInformationFilterType;
typedef itk::ExtractImageFilter< ImageType, ImageType > ExtractImageFilterType;
statisticsContainer->clear();
histogramContainer->clear();
// Make sure that mask is set
if ( maskImage == nullptr )
{
itkExceptionMacro( << "Mask image needs to be set!" );
}
// Make sure that spacing of mask and image are the same
//SpacingType imageSpacing = image->GetSpacing();
//SpacingType maskSpacing = maskImage->GetSpacing();
//PointType zeroPoint; zeroPoint.Fill( 0.0 );
//if ( (zeroPoint + imageSpacing).SquaredEuclideanDistanceTo( (zeroPoint + maskSpacing) ) > mitk::eps )
//{
// itkExceptionMacro( << "Mask needs to have same spacing as image! (Image spacing: " << imageSpacing << "; Mask spacing: " << maskSpacing << ")" );
//}
// Make sure that orientation of mask and image are the same
typedef typename ImageType::DirectionType DirectionType;
DirectionType imageDirection = image->GetDirection();
DirectionType maskDirection = maskImage->GetDirection();
for( int i = 0; i < imageDirection.ColumnDimensions; ++i )
{
for( int j = 0; j < imageDirection.ColumnDimensions; ++j )
{
double differenceDirection = imageDirection[i][j] - maskDirection[i][j];
if ( fabs( differenceDirection ) > mitk::eps )
{
double differenceDirection = imageDirection[i][j] - maskDirection[i][j];
if ( fabs( differenceDirection ) > 0.001 /*mitk::eps*/ ) // TODO: temp fix (bug 17121)
{
itkExceptionMacro( << "Mask needs to have same direction as image! (Image direction: " << imageDirection << "; Mask direction: " << maskDirection << ")" );
}
}
}
}
// Make sure that the voxels of mask and image are correctly "aligned", i.e., voxel boundaries are the same in both images
PointType imageOrigin = image->GetOrigin();
PointType maskOrigin = maskImage->GetOrigin();
long offset[ImageType::ImageDimension];
typedef itk::ContinuousIndex<double, VImageDimension> ContinousIndexType;
ContinousIndexType maskOriginContinousIndex, imageOriginContinousIndex;
image->TransformPhysicalPointToContinuousIndex(maskOrigin, maskOriginContinousIndex);
image->TransformPhysicalPointToContinuousIndex(imageOrigin, imageOriginContinousIndex);
for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i )
{
double misalignment = maskOriginContinousIndex[i] - floor( maskOriginContinousIndex[i] + 0.5 );
if ( fabs( misalignment ) > mitk::eps )
{
itkWarningMacro( << "Pixels/voxels of mask and image are not sufficiently aligned! (Misalignment: " << misalignment << ")" );
}
double indexCoordDistance = maskOriginContinousIndex[i] - imageOriginContinousIndex[i];
offset[i] = int( indexCoordDistance + image->GetBufferedRegion().GetIndex()[i] + 0.5 );
}
// Adapt the origin and region (index/size) of the mask so that the origin of both are the same
typename ChangeInformationFilterType::Pointer adaptMaskFilter;
adaptMaskFilter = ChangeInformationFilterType::New();
adaptMaskFilter->ChangeOriginOn();
adaptMaskFilter->ChangeRegionOn();
adaptMaskFilter->SetInput( maskImage );
adaptMaskFilter->SetOutputOrigin( image->GetOrigin() );
adaptMaskFilter->SetOutputOffset( offset );
+ adaptMaskFilter->SetCoordinateTolerance( 0.001 );
+ adaptMaskFilter->SetDirectionTolerance( 0.001 );
+
typename MaskImageType::Pointer adaptedMaskImage;
try
{
adaptMaskFilter->Update();
adaptedMaskImage = adaptMaskFilter->GetOutput();
}
catch( const itk::ExceptionObject &e)
{
mitkThrow() << "Attempt to adapt shifted origin of the mask image failed due to ITK Exception: \n" << e.what();
}
catch( const std::exception& e )
{
//mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what();
}
// Make sure that mask region is contained within image region
if ( adaptedMaskImage.IsNotNull() &&
!image->GetLargestPossibleRegion().IsInside( adaptedMaskImage->GetLargestPossibleRegion() ) )
{
itkWarningMacro( << "Mask region needs to be inside of image region! (Image region: "
<< image->GetLargestPossibleRegion() << "; Mask region: " << adaptedMaskImage->GetLargestPossibleRegion() << ")" );
}
// If mask region is smaller than image region, extract the sub-sampled region from the original image
typename ImageType::SizeType imageSize = image->GetBufferedRegion().GetSize();
typename ImageType::SizeType maskSize = maskImage->GetBufferedRegion().GetSize();
bool maskSmallerImage = false;
for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i )
{
if ( maskSize[i] < imageSize[i] )
{
maskSmallerImage = true;
}
}
typename ImageType::ConstPointer adaptedImage;
if ( maskSmallerImage )
{
typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New();
extractImageFilter->SetInput( image );
extractImageFilter->SetExtractionRegion( adaptedMaskImage->GetBufferedRegion() );
+ extractImageFilter->SetCoordinateTolerance( 0.001 );
+ extractImageFilter->SetDirectionTolerance( 0.001 );
extractImageFilter->Update();
adaptedImage = extractImageFilter->GetOutput();
}
else
{
adaptedImage = image;
}
// Initialize Filter
typedef itk::StatisticsImageFilter< ImageType > StatisticsFilterType;
typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New();
statisticsFilter->SetInput( adaptedImage );
try
{
statisticsFilter->Update();
}
catch( const itk::ExceptionObject& e)
{
mitkThrow() << "Image statistics initialization computation failed with ITK Exception: \n " << e.what();
}
catch( const std::exception& e )
{
//mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what();
}
// Calculate bin size or number of bins
unsigned int numberOfBins = 200; // default number of bins
double maximum = 0.0;
double minimum = 0.0;
if (m_UseBinSizeBasedOnVOIRegion)
{
maximum = statisticsFilter->GetMaximum();
minimum = statisticsFilter->GetMinimum();
if (m_UseDefaultBinSize)
{
m_HistogramBinSize = std::ceil( static_cast<double>((statisticsFilter->GetMaximum() - statisticsFilter->GetMinimum() + 1)/numberOfBins) );
}
else
{
numberOfBins = calcNumberOfBins(statisticsFilter->GetMinimum(), statisticsFilter->GetMaximum());
}
}
else
{
double sig = 0.0;
int counter = 0;
//Find the min and max values for the Roi to set the range for the histogram
GetMinAndMaxValue( minimum, maximum, counter, sig, image, maskImage);
numberOfBins = maximum - minimum;
if(maximum - minimum <= 10)
{
numberOfBins = 100;
}
}
- typename LabelStatisticsFilterType::Pointer labelStatisticsFilter = LabelStatisticsFilterType::New();
+ typename LabelStatisticsFilterType::Pointer labelStatisticsFilter = LabelStatisticsFilterType::New();
labelStatisticsFilter->SetInput( adaptedImage );
labelStatisticsFilter->SetLabelInput( adaptedMaskImage );
+ labelStatisticsFilter->SetCoordinateTolerance( 0.001 );
+ labelStatisticsFilter->SetDirectionTolerance( 0.001 );
labelStatisticsFilter->UseHistogramsOn();
labelStatisticsFilter->SetHistogramParameters( numberOfBins, floor(minimum), ceil(maximum) ); //statisticsFilter->GetMinimum() statisticsFilter->GetMaximum()
// Add progress listening
typedef itk::SimpleMemberCommand< ImageStatisticsCalculator > ITKCommandType;
ITKCommandType::Pointer progressListener;
progressListener = ITKCommandType::New();
progressListener->SetCallbackFunction( this,
&ImageStatisticsCalculator::MaskedStatisticsProgressUpdate );
unsigned long observerTag = labelStatisticsFilter->AddObserver(
itk::ProgressEvent(), progressListener );
// Execute filter
this->InvokeEvent( itk::StartEvent() );
// Make sure that only the mask region is considered (otherwise, if the mask region is smaller
// than the image region, the Update() would result in an exception).
labelStatisticsFilter->GetOutput()->SetRequestedRegion( adaptedMaskImage->GetLargestPossibleRegion() );
// Execute the filter
try
{
labelStatisticsFilter->Update();
}
catch( const itk::ExceptionObject& e)
{
mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what();
}
catch( const std::exception& e )
{
//mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what();
}
this->InvokeEvent( itk::EndEvent() );
if( observerTag )
labelStatisticsFilter->RemoveObserver( observerTag );
// Find all relevant labels of mask (other than 0)
std::list< int > relevantLabels = labelStatisticsFilter->GetRelevantLabels();
unsigned int i;
if ( labelStatisticsFilter->GetMaskingNonEmpty() )
{
std::list< int >::iterator it;
for ( it = relevantLabels.begin(), i = 0;
it != relevantLabels.end();
++it, ++i )
{
Statistics statistics; // restore previous code
labelStatisticsFilter->GetHistogram(*it) ;
histogramContainer->push_back( HistogramType::ConstPointer( labelStatisticsFilter->GetHistogram( (*it) ) ) );
statistics.SetLabel (*it);
statistics.SetN(labelStatisticsFilter->GetCount( *it ));
statistics.SetMin(labelStatisticsFilter->GetMinimum( *it ));
statistics.SetMax(labelStatisticsFilter->GetMaximum( *it ));
statistics.SetMean(labelStatisticsFilter->GetMean( *it ));
statistics.SetMedian(labelStatisticsFilter->GetMedian( *it));
statistics.SetMedian(labelStatisticsFilter->GetMedian( *it ));
statistics.SetVariance(labelStatisticsFilter->GetVariance( *it ));
statistics.SetSigma(labelStatisticsFilter->GetSigma( *it ));
statistics.SetSkewness(labelStatisticsFilter->GetSkewness( *it ));
statistics.SetKurtosis(labelStatisticsFilter->GetKurtosis( *it ));
statistics.SetUniformity( labelStatisticsFilter->GetUniformity( *it ));
statistics.SetEntropy( labelStatisticsFilter->GetEntropy( *it ));
statistics.SetUPP( labelStatisticsFilter->GetUPP( *it));
statistics.SetMPP( labelStatisticsFilter->GetMPP( *it));
statistics.SetRMS(sqrt( statistics.GetMean() * statistics.GetMean()
+ statistics.GetSigma() * statistics.GetSigma() ));
// restrict image to mask area for min/max index calculation
typedef itk::MaskImageFilter< ImageType, MaskImageType, ImageType > MaskImageFilterType;
typename MaskImageFilterType::Pointer masker = MaskImageFilterType::New();
bool isMinAndMaxSameValue = (statistics.GetMin() == statistics.GetMax());
// bug 17962: following is a workaround for the case when min and max are the same, we can probably find a nicer way here
double outsideValue = (isMinAndMaxSameValue ? (statistics.GetMax()/2) : (statistics.GetMin()+statistics.GetMax())/2);
masker->SetOutsideValue( outsideValue );
masker->SetInput1(adaptedImage);
masker->SetInput2(adaptedMaskImage);
+ masker->SetCoordinateTolerance( 0.001 );
+ masker->SetDirectionTolerance( 0.001 );
masker->Update();
// get index of minimum and maximum
typedef itk::MinimumMaximumImageCalculator< ImageType > MinMaxFilterType;
typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New();
minMaxFilter->SetImage( masker->GetOutput() );
unsigned long observerTag2 = minMaxFilter->AddObserver( itk::ProgressEvent(), progressListener );
minMaxFilter->Compute();
minMaxFilter->RemoveObserver( observerTag2 );
this->InvokeEvent( itk::EndEvent() );
typename MinMaxFilterType::IndexType tempMaxIndex = minMaxFilter->GetIndexOfMaximum();
// bug 17962: following is a workaround for the case when min and max are the same, we can probably find a nicer way here
typename MinMaxFilterType::IndexType tempMinIndex =
(isMinAndMaxSameValue ? minMaxFilter->GetIndexOfMaximum() : minMaxFilter->GetIndexOfMinimum());
// FIX BUG 14644
//If a PlanarFigure is used for segmentation the
//adaptedImage is a single slice (2D). Adding the
// 3. dimension.
vnl_vector<int> maxIndex;
vnl_vector<int> minIndex;
maxIndex.set_size(m_Image->GetDimension());
minIndex.set_size(m_Image->GetDimension());
if (m_MaskingMode == MASKING_MODE_PLANARFIGURE && m_Image->GetDimension()==3)
{
maxIndex[m_PlanarFigureCoordinate0] = tempMaxIndex[0];
maxIndex[m_PlanarFigureCoordinate1] = tempMaxIndex[1];
maxIndex[m_PlanarFigureAxis] = m_PlanarFigureSlice;
minIndex[m_PlanarFigureCoordinate0] = tempMinIndex[0] ;
minIndex[m_PlanarFigureCoordinate1] = tempMinIndex[1];
minIndex[m_PlanarFigureAxis] = m_PlanarFigureSlice;
} else
{
for (unsigned int i = 0; i<maxIndex.size(); i++)
{
maxIndex[i] = tempMaxIndex[i];
minIndex[i] = tempMinIndex[i];
}
}
// FIX END
statistics.SetMaxIndex(maxIndex);
statistics.SetMinIndex(minIndex);
/*****************************************************Calculate Hotspot Statistics**********************************************/
if(IsHotspotCalculated() && VImageDimension == 3)
{
bool isDefined(false);
Statistics hotspotStatistics = CalculateHotspotStatistics(adaptedImage.GetPointer(), adaptedMaskImage.GetPointer(),GetHotspotRadiusInMM(), isDefined, *it);
statistics.GetHotspotStatistics() = hotspotStatistics;
if(statistics.GetHotspotStatistics().HasHotspotStatistics())
{
MITK_DEBUG << "Hotspot statistics available";
statistics.SetHotspotIndex( hotspotStatistics.GetHotspotIndex() );
}
else
{
MITK_ERROR << "No hotspot statistics available!";
}
}
statisticsContainer->push_back( statistics );
}
}
else
{
histogramContainer->push_back( HistogramType::ConstPointer( m_EmptyHistogram ) );
statisticsContainer->push_back( Statistics() );
}
}
template <typename TPixel, unsigned int VImageDimension >
ImageStatisticsCalculator::ImageExtrema
ImageStatisticsCalculator::CalculateExtremaWorld(
const itk::Image<TPixel, VImageDimension> *inputImage,
itk::Image<unsigned short, VImageDimension> *maskImage,
double neccessaryDistanceToImageBorderInMM,
unsigned int label)
{
typedef itk::Image< TPixel, VImageDimension > ImageType;
typedef itk::Image< unsigned short, VImageDimension > MaskImageType;
typedef itk::ImageRegionConstIteratorWithIndex<MaskImageType> MaskImageIteratorType;
typedef itk::ImageRegionConstIteratorWithIndex<ImageType> InputImageIndexIteratorType;
typename ImageType::SpacingType spacing = inputImage->GetSpacing();
ImageExtrema minMax;
minMax.Defined = false;
minMax.MaxIndex.set_size(VImageDimension);
minMax.MaxIndex.set_size(VImageDimension);
typename ImageType::RegionType allowedExtremaRegion = inputImage->GetLargestPossibleRegion();
bool keepDistanceToImageBorders( neccessaryDistanceToImageBorderInMM > 0 );
if (keepDistanceToImageBorders)
{
long distanceInPixels[VImageDimension];
for(unsigned short dimension = 0; dimension < VImageDimension; ++dimension)
{
// To confirm that the whole hotspot is inside the image we have to keep a specific distance to the image-borders, which is as long as
// the radius. To get the amount of indices we divide the radius by spacing and add 0.5 because voxels are center based:
// For example with a radius of 2.2 and a spacing of 1 two indices are enough because 2.2 / 1 + 0.5 = 2.7 => 2.
// But with a radius of 2.7 we need 3 indices because 2.7 / 1 + 0.5 = 3.2 => 3
distanceInPixels[dimension] = int( neccessaryDistanceToImageBorderInMM / spacing[dimension] + 0.5);
}
allowedExtremaRegion.ShrinkByRadius(distanceInPixels);
}
InputImageIndexIteratorType imageIndexIt(inputImage, allowedExtremaRegion);
float maxValue = itk::NumericTraits<float>::min();
float minValue = itk::NumericTraits<float>::max();
typename ImageType::IndexType maxIndex;
typename ImageType::IndexType minIndex;
for(unsigned short i = 0; i < VImageDimension; ++i)
{
maxIndex[i] = 0;
minIndex[i] = 0;
}
if (maskImage != nullptr)
{
MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion());
typename ImageType::IndexType imageIndex;
typename ImageType::PointType worldPosition;
typename ImageType::IndexType maskIndex;
for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt)
{
imageIndex = maskIndex = maskIt.GetIndex();
if(maskIt.Get() == label)
{
if( allowedExtremaRegion.IsInside(imageIndex) )
{
imageIndexIt.SetIndex( imageIndex );
double value = imageIndexIt.Get();
minMax.Defined = true;
//Calculate minimum, maximum and corresponding index-values
if( value > maxValue )
{
maxIndex = imageIndexIt.GetIndex();
maxValue = value;
}
if(value < minValue )
{
minIndex = imageIndexIt.GetIndex();
minValue = value;
}
}
}
}
}
else
{
for(imageIndexIt.GoToBegin(); !imageIndexIt.IsAtEnd(); ++imageIndexIt)
{
double value = imageIndexIt.Get();
minMax.Defined = true;
//Calculate minimum, maximum and corresponding index-values
if( value > maxValue )
{
maxIndex = imageIndexIt.GetIndex();
maxValue = value;
}
if(value < minValue )
{
minIndex = imageIndexIt.GetIndex();
minValue = value;
}
}
}
minMax.MaxIndex.set_size(VImageDimension);
minMax.MinIndex.set_size(VImageDimension);
for(unsigned int i = 0; i < minMax.MaxIndex.size(); ++i)
{
minMax.MaxIndex[i] = maxIndex[i];
}
for(unsigned int i = 0; i < minMax.MinIndex.size(); ++i)
{
minMax.MinIndex[i] = minIndex[i];
}
minMax.Max = maxValue;
minMax.Min = minValue;
return minMax;
}
template <unsigned int VImageDimension>
itk::Size<VImageDimension>
ImageStatisticsCalculator
::CalculateConvolutionKernelSize(double spacing[VImageDimension], double radiusInMM)
{
typedef itk::Image< float, VImageDimension > KernelImageType;
typedef typename KernelImageType::SizeType SizeType;
SizeType maskSize;
for(unsigned int i = 0; i < VImageDimension; ++i)
{
maskSize[i] = static_cast<int>( 2 * radiusInMM / spacing[i]);
// We always want an uneven size to have a clear center point in the convolution mask
if(maskSize[i] % 2 == 0 )
{
++maskSize[i];
}
}
return maskSize;
}
template <unsigned int VImageDimension>
itk::SmartPointer< itk::Image<float, VImageDimension> >
ImageStatisticsCalculator
::GenerateHotspotSearchConvolutionKernel(double mmPerPixel[VImageDimension], double radiusInMM)
{
std::stringstream ss;
for (unsigned int i = 0; i < VImageDimension; ++i)
{
ss << mmPerPixel[i];
if (i < VImageDimension -1)
ss << ",";
}
MITK_DEBUG << "Update convolution kernel for spacing (" << ss.str() << ") and radius " << radiusInMM << "mm";
double radiusInMMSquared = radiusInMM * radiusInMM;
typedef itk::Image< float, VImageDimension > KernelImageType;
typename KernelImageType::Pointer convolutionKernel = KernelImageType::New();
// Calculate size and allocate mask image
typedef typename KernelImageType::SizeType SizeType;
SizeType maskSize = this->CalculateConvolutionKernelSize<VImageDimension>(mmPerPixel, radiusInMM);
Point3D convolutionMaskCenterIndex; convolutionMaskCenterIndex.Fill(0.0);
for(unsigned int i = 0; i < VImageDimension; ++i)
{
convolutionMaskCenterIndex[i] = 0.5 * (double)(maskSize[i]-1);
}
typedef typename KernelImageType::IndexType IndexType;
IndexType maskIndex;
maskIndex.Fill(0);
typedef typename KernelImageType::RegionType RegionType;
RegionType maskRegion;
maskRegion.SetSize(maskSize);
maskRegion.SetIndex(maskIndex);
convolutionKernel->SetRegions(maskRegion);
convolutionKernel->SetSpacing(mmPerPixel);
convolutionKernel->Allocate();
// Fill mask image values by subsampling the image grid
typedef itk::ImageRegionIteratorWithIndex<KernelImageType> MaskIteratorType;
MaskIteratorType maskIt(convolutionKernel,maskRegion);
int numberOfSubVoxelsPerDimension = 2; // per dimension!
int numberOfSubVoxels = ::pow( static_cast<float>(numberOfSubVoxelsPerDimension), static_cast<float>(VImageDimension) );
double subVoxelSizeInPixels = 1.0 / (double)numberOfSubVoxelsPerDimension;
double valueOfOneSubVoxel = 1.0 / (double)numberOfSubVoxels;
double maskValue = 0.0;
Point3D subVoxelIndexPosition;
double distanceSquared = 0.0;
typedef itk::ContinuousIndex<double, VImageDimension> ContinuousIndexType;
for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt)
{
ContinuousIndexType indexPoint(maskIt.GetIndex());
Point3D voxelPosition;
for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension)
{
voxelPosition[dimension] = indexPoint[dimension];
}
maskValue = 0.0;
Vector3D subVoxelOffset; subVoxelOffset.Fill(0.0);
// iterate sub-voxels by iterating all possible offsets
for (subVoxelOffset[0] = -0.5 + subVoxelSizeInPixels / 2.0;
subVoxelOffset[0] < +0.5;
subVoxelOffset[0] += subVoxelSizeInPixels)
{
for (subVoxelOffset[1] = -0.5 + subVoxelSizeInPixels / 2.0;
subVoxelOffset[1] < +0.5;
subVoxelOffset[1] += subVoxelSizeInPixels)
{
for (subVoxelOffset[2] = -0.5 + subVoxelSizeInPixels / 2.0;
subVoxelOffset[2] < +0.5;
subVoxelOffset[2] += subVoxelSizeInPixels)
{
subVoxelIndexPosition = voxelPosition + subVoxelOffset; // this COULD be integrated into the for-loops if neccessary (add voxelPosition to initializer and end condition)
distanceSquared =
(subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] * (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0]
+ (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] * (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1]
+ (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2] * (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2];
if (distanceSquared <= radiusInMMSquared)
{
maskValue += valueOfOneSubVoxel;
}
}
}
}
maskIt.Set( maskValue );
}
return convolutionKernel;
}
template <typename TPixel, unsigned int VImageDimension>
itk::SmartPointer<itk::Image<TPixel, VImageDimension> >
ImageStatisticsCalculator::GenerateConvolutionImage( const itk::Image<TPixel, VImageDimension>* inputImage )
{
double mmPerPixel[VImageDimension];
for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension)
{
mmPerPixel[dimension] = inputImage->GetSpacing()[dimension];
}
// update convolution kernel
typedef itk::Image< float, VImageDimension > KernelImageType;
typename KernelImageType::Pointer convolutionKernel = this->GenerateHotspotSearchConvolutionKernel<VImageDimension>(mmPerPixel, m_HotspotRadiusInMM);
// update convolution image
typedef itk::Image< TPixel, VImageDimension > InputImageType;
typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType;
typedef itk::FFTConvolutionImageFilter<InputImageType,
KernelImageType,
ConvolutionImageType> ConvolutionFilterType;
typename ConvolutionFilterType::Pointer convolutionFilter = ConvolutionFilterType::New();
typedef itk::ConstantBoundaryCondition<InputImageType, InputImageType> BoundaryConditionType;
BoundaryConditionType boundaryCondition;
boundaryCondition.SetConstant(0.0);
if (GetHotspotMustBeCompletlyInsideImage())
{
// overwrite default boundary condition
convolutionFilter->SetBoundaryCondition(&boundaryCondition);
}
convolutionFilter->SetInput(inputImage);
convolutionFilter->SetKernelImage(convolutionKernel);
convolutionFilter->SetNormalize(true);
MITK_DEBUG << "Update Convolution image for hotspot search";
convolutionFilter->UpdateLargestPossibleRegion();
typename ConvolutionImageType::Pointer convolutionImage = convolutionFilter->GetOutput();
convolutionImage->SetSpacing( inputImage->GetSpacing() ); // only workaround because convolution filter seems to ignore spacing of input image
m_HotspotRadiusInMMChanged = false;
return convolutionImage;
}
template < typename TPixel, unsigned int VImageDimension>
void
ImageStatisticsCalculator
::FillHotspotMaskPixels( itk::Image<TPixel, VImageDimension>* maskImage,
itk::Point<double, VImageDimension> sphereCenter,
double sphereRadiusInMM)
{
typedef itk::Image< TPixel, VImageDimension > MaskImageType;
typedef itk::ImageRegionIteratorWithIndex<MaskImageType> MaskImageIteratorType;
MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion());
typename MaskImageType::IndexType maskIndex;
typename MaskImageType::PointType worldPosition;
for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt)
{
maskIndex = maskIt.GetIndex();
maskImage->TransformIndexToPhysicalPoint(maskIndex, worldPosition);
maskIt.Set( worldPosition.EuclideanDistanceTo(sphereCenter) <= sphereRadiusInMM ? 1 : 0 );
}
}
template < typename TPixel, unsigned int VImageDimension>
ImageStatisticsCalculator::Statistics
ImageStatisticsCalculator::CalculateHotspotStatistics(
const itk::Image<TPixel, VImageDimension>* inputImage,
itk::Image<unsigned short, VImageDimension>* maskImage,
double radiusInMM,
bool& isHotspotDefined,
unsigned int label)
{
// get convolution image (updated in GenerateConvolutionImage())
typedef itk::Image< TPixel, VImageDimension > InputImageType;
typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType;
typedef itk::Image< float, VImageDimension > KernelImageType;
typedef itk::Image< unsigned short, VImageDimension > MaskImageType;
//typename ConvolutionImageType::Pointer convolutionImage = dynamic_cast<ConvolutionImageType*>(this->GenerateConvolutionImage(inputImage));
typename ConvolutionImageType::Pointer convolutionImage = this->GenerateConvolutionImage(inputImage);
if (convolutionImage.IsNull())
{
MITK_ERROR << "Empty convolution image in CalculateHotspotStatistics(). We should never reach this state (logic error).";
throw std::logic_error("Empty convolution image in CalculateHotspotStatistics()");
}
// find maximum in convolution image, given the current mask
double requiredDistanceToBorder = m_HotspotMustBeCompletelyInsideImage ? m_HotspotRadiusInMM : -1.0;
ImageExtrema convolutionImageInformation = CalculateExtremaWorld(convolutionImage.GetPointer(), maskImage, requiredDistanceToBorder, label);
isHotspotDefined = convolutionImageInformation.Defined;
if (!isHotspotDefined)
{
m_EmptyStatistics.Reset(VImageDimension);
MITK_ERROR << "No origin of hotspot-sphere was calculated! Returning empty statistics";
return m_EmptyStatistics;
}
else
{
// create a binary mask around the "hotspot" region, fill the shape of a sphere around our hotspot center
typedef itk::ImageDuplicator< InputImageType > DuplicatorType;
typename DuplicatorType::Pointer copyMachine = DuplicatorType::New();
copyMachine->SetInputImage(inputImage);
copyMachine->Update();
typedef itk::CastImageFilter< InputImageType, MaskImageType > CastFilterType;
typename CastFilterType::Pointer caster = CastFilterType::New();
caster->SetInput( copyMachine->GetOutput() );
caster->Update();
typename MaskImageType::Pointer hotspotMaskITK = caster->GetOutput();
typedef typename InputImageType::IndexType IndexType;
IndexType maskCenterIndex;
for (unsigned int d =0; d< VImageDimension;++d) maskCenterIndex[d]=convolutionImageInformation.MaxIndex[d];
typename ConvolutionImageType::PointType maskCenter;
inputImage->TransformIndexToPhysicalPoint(maskCenterIndex,maskCenter);
this->FillHotspotMaskPixels(hotspotMaskITK.GetPointer(), maskCenter, radiusInMM);
// calculate statistics within the binary mask
typedef itk::ExtendedLabelStatisticsImageFilter< InputImageType, MaskImageType> LabelStatisticsFilterType;
typename LabelStatisticsFilterType::Pointer labelStatisticsFilter;
labelStatisticsFilter = LabelStatisticsFilterType::New();
labelStatisticsFilter->SetInput( inputImage );
labelStatisticsFilter->SetLabelInput( hotspotMaskITK );
+ labelStatisticsFilter->SetCoordinateTolerance( 0.001 );
+ labelStatisticsFilter->SetDirectionTolerance( 0.001 );
+
labelStatisticsFilter->Update();
Statistics hotspotStatistics;
hotspotStatistics.SetHotspotIndex(convolutionImageInformation.MaxIndex);
hotspotStatistics.SetMean(convolutionImageInformation.Max);
if ( labelStatisticsFilter->HasLabel( 1 ) )
{
hotspotStatistics.SetLabel (1);
hotspotStatistics.SetN(labelStatisticsFilter->GetCount(1));
hotspotStatistics.SetMin(labelStatisticsFilter->GetMinimum(1));
hotspotStatistics.SetMax(labelStatisticsFilter->GetMaximum(1));
hotspotStatistics.SetMedian(labelStatisticsFilter->GetMedian(1));
hotspotStatistics.SetVariance(labelStatisticsFilter->GetVariance(1));
hotspotStatistics.SetSigma(labelStatisticsFilter->GetSigma(1));
hotspotStatistics.SetRMS(sqrt( hotspotStatistics.GetMean() * hotspotStatistics.GetMean()
+ hotspotStatistics.GetSigma() * hotspotStatistics.GetSigma() ));
MITK_DEBUG << "Statistics for inside hotspot: Mean " << hotspotStatistics.GetMean()
<< ", SD " << hotspotStatistics.GetSigma()
<< ", Max " << hotspotStatistics.GetMax()
<< ", Min " << hotspotStatistics.GetMin();
}
else
{
MITK_ERROR << "Uh oh! Unable to calculate statistics for hotspot region...";
return m_EmptyStatistics;
}
return hotspotStatistics;
}
}
template < typename TPixel, unsigned int VImageDimension >
void ImageStatisticsCalculator::InternalCalculateMaskFromPlanarFigure(
const itk::Image< TPixel, VImageDimension > *image, unsigned int axis )
{
typedef itk::Image< TPixel, VImageDimension > ImageType;
typedef itk::CastImageFilter< ImageType, MaskImage2DType > CastFilterType;
// Generate mask image as new image with same header as input image and
// initialize with 1.
typename CastFilterType::Pointer castFilter = CastFilterType::New();
castFilter->SetInput( image );
castFilter->Update();
castFilter->GetOutput()->FillBuffer( 1 );
// all PolylinePoints of the PlanarFigure are stored in a vtkPoints object.
// These points are used by the vtkLassoStencilSource to create
// a vtkImageStencil.
const mitk::PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry();
const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 );
const mitk::BaseGeometry *imageGeometry3D = m_Image->GetGeometry( 0 );
// If there is a second poly line in a closed planar figure, treat it as a hole.
PlanarFigure::PolyLineType planarFigureHolePolyline;
if (m_PlanarFigure->GetPolyLinesSize() == 2)
planarFigureHolePolyline = m_PlanarFigure->GetPolyLine(1);
// Determine x- and y-dimensions depending on principal axis
int i0, i1;
switch ( axis )
{
case 0:
i0 = 1;
i1 = 2;
break;
case 1:
i0 = 0;
i1 = 2;
break;
case 2:
default:
i0 = 0;
i1 = 1;
break;
}
m_PlanarFigureCoordinate0= i0;
m_PlanarFigureCoordinate1= i1;
// store the polyline contour as vtkPoints object
bool outOfBounds = false;
vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
typename PlanarFigure::PolyLineType::const_iterator it;
for ( it = planarFigurePolyline.begin();
it != planarFigurePolyline.end();
++it )
{
Point3D point3D;
// Convert 2D point back to the local index coordinates of the selected
// image
planarFigurePlaneGeometry->Map( *it, point3D );
// Polygons (partially) outside of the image bounds can not be processed
// further due to a bug in vtkPolyDataToImageStencil
if ( !imageGeometry3D->IsInside( point3D ) )
{
outOfBounds = true;
}
imageGeometry3D->WorldToIndex( point3D, point3D );
points->InsertNextPoint( point3D[i0], point3D[i1], 0 );
}
vtkSmartPointer<vtkPoints> holePoints = nullptr;
if (!planarFigureHolePolyline.empty())
{
holePoints = vtkSmartPointer<vtkPoints>::New();
Point3D point3D;
PlanarFigure::PolyLineType::const_iterator end = planarFigureHolePolyline.end();
for (it = planarFigureHolePolyline.begin(); it != end; ++it)
{
planarFigurePlaneGeometry->Map(*it, point3D);
imageGeometry3D->WorldToIndex(point3D, point3D);
holePoints->InsertNextPoint(point3D[i0], point3D[i1], 0);
}
}
// mark a malformed 2D planar figure ( i.e. area = 0 ) as out of bounds
// this can happen when all control points of a rectangle lie on the same line = two of the three extents are zero
double bounds[6] = {0, 0, 0, 0, 0, 0};
points->GetBounds( bounds );
bool extent_x = (fabs(bounds[0] - bounds[1])) < mitk::eps;
bool extent_y = (fabs(bounds[2] - bounds[3])) < mitk::eps;
bool extent_z = (fabs(bounds[4] - bounds[5])) < mitk::eps;
// throw an exception if a closed planar figure is deformed, i.e. has only one non-zero extent
if ( m_PlanarFigure->IsClosed() &&
((extent_x && extent_y) || (extent_x && extent_z) || (extent_y && extent_z)))
{
mitkThrow() << "Figure has a zero area and cannot be used for masking.";
}
if ( outOfBounds )
{
throw std::runtime_error( "Figure at least partially outside of image bounds!" );
}
// create a vtkLassoStencilSource and set the points of the Polygon
vtkSmartPointer<vtkLassoStencilSource> lassoStencil = vtkSmartPointer<vtkLassoStencilSource>::New();
lassoStencil->SetShapeToPolygon();
lassoStencil->SetPoints( points );
vtkSmartPointer<vtkLassoStencilSource> holeLassoStencil = nullptr;
if (holePoints.GetPointer() != nullptr)
{
holeLassoStencil = vtkSmartPointer<vtkLassoStencilSource>::New();
holeLassoStencil->SetShapeToPolygon();
holeLassoStencil->SetPoints(holePoints);
}
// Export from ITK to VTK (to use a VTK filter)
typedef itk::VTKImageImport< MaskImage2DType > ImageImportType;
typedef itk::VTKImageExport< MaskImage2DType > ImageExportType;
typename ImageExportType::Pointer itkExporter = ImageExportType::New();
itkExporter->SetInput( castFilter->GetOutput() );
vtkSmartPointer<vtkImageImport> vtkImporter = vtkSmartPointer<vtkImageImport>::New();
this->ConnectPipelines( itkExporter, vtkImporter );
// Apply the generated image stencil to the input image
vtkSmartPointer<vtkImageStencil> imageStencilFilter = vtkSmartPointer<vtkImageStencil>::New();
imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() );
imageStencilFilter->SetStencilConnection(lassoStencil->GetOutputPort());
imageStencilFilter->ReverseStencilOff();
imageStencilFilter->SetBackgroundValue( 0 );
imageStencilFilter->Update();
vtkSmartPointer<vtkImageStencil> holeStencilFilter = nullptr;
if (holeLassoStencil.GetPointer() != nullptr)
{
holeStencilFilter = vtkSmartPointer<vtkImageStencil>::New();
holeStencilFilter->SetInputConnection(imageStencilFilter->GetOutputPort());
holeStencilFilter->SetStencilConnection(holeLassoStencil->GetOutputPort());
holeStencilFilter->ReverseStencilOn();
holeStencilFilter->SetBackgroundValue(0);
holeStencilFilter->Update();
}
// Export from VTK back to ITK
vtkSmartPointer<vtkImageExport> vtkExporter = vtkSmartPointer<vtkImageExport>::New();
vtkExporter->SetInputConnection( holeStencilFilter.GetPointer() == nullptr
? imageStencilFilter->GetOutputPort()
: holeStencilFilter->GetOutputPort());
vtkExporter->Update();
typename ImageImportType::Pointer itkImporter = ImageImportType::New();
this->ConnectPipelines( vtkExporter, itkImporter );
itkImporter->Update();
typedef itk::ImageDuplicator< ImageImportType::OutputImageType > DuplicatorType;
DuplicatorType::Pointer duplicator = DuplicatorType::New();
duplicator->SetInputImage( itkImporter->GetOutput() );
duplicator->Update();
// Store mask
m_InternalImageMask2D = duplicator->GetOutput();
}
void ImageStatisticsCalculator::UnmaskedStatisticsProgressUpdate()
{
// Need to throw away every second progress event to reach a final count of
// 100 since two consecutive filters are used in this case
static int updateCounter = 0;
if ( updateCounter++ % 2 == 0 )
{
this->InvokeEvent( itk::ProgressEvent() );
}
}
void ImageStatisticsCalculator::MaskedStatisticsProgressUpdate()
{
this->InvokeEvent( itk::ProgressEvent() );
}
}