diff --git a/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp b/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp
index da07560989..80efa6ab26 100644
--- a/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp
+++ b/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp
@@ -1,595 +1,603 @@
#include <mitkHotspotMaskGenerator.h>
#include <mitkImageTimeSelector.h>
#include <mitkImageCast.h>
#include <mitkPoint.h>
#include <itkImageRegionIterator.h>
#include "mitkImageAccessByItk.h"
#include <itkImageDuplicator.h>
#include <itkFFTConvolutionImageFilter.h>
#include <mitkITKImageImport.h>
namespace mitk
{
HotspotMaskGenerator::HotspotMaskGenerator():
m_HotspotRadiusinMM(6.2035049089940), // radius of a 1cm3 sphere in mm
m_HotspotMustBeCompletelyInsideImage(true),
m_Label(1)
{
m_TimeStep = 0;
m_InternalMask = mitk::Image::New();
m_InternalMaskUpdateTime = 0;
}
void HotspotMaskGenerator::SetInputImage(mitk::Image::Pointer inputImage)
{
if (inputImage != m_inputImage)
{
m_inputImage = inputImage;
m_ConvolutionImageMaxIndex.set_size(inputImage->GetDimension());
m_ConvolutionImageMinIndex.set_size(inputImage->GetDimension());
this->Modified();
}
}
void HotspotMaskGenerator::SetMask(MaskGenerator::Pointer mask)
{
if (mask != m_Mask)
{
m_Mask = mask;
this->Modified();
}
}
HotspotMaskGenerator::~HotspotMaskGenerator()
{
}
void HotspotMaskGenerator::SetHotspotRadiusInMM(double radiusInMillimeter)
{
if(radiusInMillimeter != m_HotspotRadiusinMM)
{
m_HotspotRadiusinMM = radiusInMillimeter;
this->Modified();
}
}
const double& HotspotMaskGenerator::GetHotspotRadiusinMM() const
{
return m_HotspotRadiusinMM;
}
bool HotspotMaskGenerator::GetHotspotMustBeCompletelyInsideImage() const
{
return m_HotspotMustBeCompletelyInsideImage;
}
void HotspotMaskGenerator::SetHotspotMustBeCompletelyInsideImage(bool mustBeCompletelyInImage)
{
if (m_HotspotMustBeCompletelyInsideImage != mustBeCompletelyInImage)
{
m_HotspotMustBeCompletelyInsideImage = mustBeCompletelyInImage;
this->Modified();
}
}
mitk::Image::Pointer HotspotMaskGenerator::GetMask()
{
if (IsUpdateRequired())
{
if ( m_inputImage.IsNull() )
{
throw std::runtime_error( "Error: image empty!" );
}
if ( m_TimeStep >= m_inputImage->GetTimeSteps() )
{
throw std::runtime_error( "Error: invalid time step!" );
}
mitk::ImageTimeSelector::Pointer imageTimeSelector = mitk::ImageTimeSelector::New();
imageTimeSelector->SetInput( m_inputImage );
imageTimeSelector->SetTimeNr( m_TimeStep );
imageTimeSelector->UpdateLargestPossibleRegion();
mitk::Image::Pointer timeSliceImage = imageTimeSelector->GetOutput();
m_internalImage = timeSliceImage;
m_internalMask2D = nullptr; // is this correct when this variable holds a smart pointer?
m_internalMask3D = nullptr;
if ( m_Mask != nullptr )
{
m_Mask->SetTimeStep(m_TimeStep);
mitk::Image::Pointer timeSliceMask = m_Mask->GetMask();
if ( m_internalImage->GetDimension() == 3 )
{
CastToItkImage(timeSliceMask, m_internalMask3D);
AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 3, m_internalMask3D, m_Label);
}
else if ( m_internalImage->GetDimension() == 2 )
{
CastToItkImage(timeSliceMask, m_internalMask2D);
AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 2, m_internalMask2D, m_Label);
}
else
{
throw std::runtime_error( "Error: invalid image dimension" );
}
}
else
{
if ( m_internalImage->GetDimension() == 3 )
{
AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 3, m_internalMask3D, m_Label);
}
else if ( m_internalImage->GetDimension() == 2 )
{
AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 2, m_internalMask2D, m_Label);
}
else
{
throw std::runtime_error( "Error: invalid image dimension" );
}
}
this->Modified();
}
m_InternalMaskUpdateTime = m_InternalMask->GetMTime();
return m_InternalMask;
}
+ void HotspotMaskGenerator::SetTimeStep(unsigned int timeStep)
+ {
+ if (m_TimeStep != timeStep)
+ {
+ m_TimeStep = timeStep;
+ }
+ }
+
void HotspotMaskGenerator::SetLabel(unsigned short label)
{
if (label != m_Label)
{
m_Label = label;
this->Modified();
}
}
vnl_vector<int> HotspotMaskGenerator::GetConvolutionImageMinIndex()
{
this->GetMask(); // make sure we are up to date
return m_ConvolutionImageMinIndex;
}
vnl_vector<int> HotspotMaskGenerator::GetHotspotIndex()
{
this->GetMask(); // make sure we are up to date
return m_ConvolutionImageMaxIndex;
}
template <typename TPixel, unsigned int VImageDimension >
HotspotMaskGenerator::ImageExtrema
HotspotMaskGenerator::CalculateExtremaWorld( const itk::Image<TPixel, VImageDimension>* inputImage,
typename itk::Image<unsigned short, VImageDimension>::Pointer 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>
HotspotMaskGenerator::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> >
HotspotMaskGenerator::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);
mitk::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;
mitk::Point3D subVoxelIndexPosition;
double distanceSquared = 0.0;
typedef itk::ContinuousIndex<double, VImageDimension> ContinuousIndexType;
for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt)
{
ContinuousIndexType indexPoint(maskIt.GetIndex());
mitk::Point3D voxelPosition;
for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension)
{
voxelPosition[dimension] = indexPoint[dimension];
}
maskValue = 0.0;
mitk::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> >
HotspotMaskGenerator::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 (m_HotspotMustBeCompletelyInsideImage)
{
// 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
return convolutionImage;
}
template < typename TPixel, unsigned int VImageDimension>
void
HotspotMaskGenerator::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;
// this is not very smart. I would rather use a 0 initialized mask (not the case here -> blame CalculateHotspotMask) and find the region where I need to iterate over, then iterate only over the small region
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>
void
HotspotMaskGenerator::CalculateHotspotMask(itk::Image<TPixel, VImageDimension>* inputImage,
typename itk::Image<unsigned short, VImageDimension>::Pointer maskImage,
unsigned int label)
{
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 = 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()");
}
// if mask image is not defined, create an image of the same size as inputImage and fill it with 1's
// there is maybe a better way to do this!?
if (maskImage == nullptr)
{
maskImage = MaskImageType::New();
typename MaskImageType::RegionType maskRegion = inputImage->GetLargestPossibleRegion();
typename MaskImageType::SpacingType maskSpacing = inputImage->GetSpacing();
typename MaskImageType::PointType maskOrigin = inputImage->GetOrigin();
typename MaskImageType::DirectionType maskDirection = inputImage->GetDirection();
maskImage->SetRegions(maskRegion);
maskImage->Allocate();
maskImage->SetOrigin(maskOrigin);
maskImage->SetSpacing(maskSpacing);
maskImage->SetDirection(maskDirection);
maskImage->FillBuffer(1);
label = 1;
}
// 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);
bool isHotspotDefined = convolutionImageInformation.Defined;
if (!isHotspotDefined)
{
MITK_ERROR << "No origin of hotspot-sphere was calculated!";
m_InternalMask = nullptr;
}
else
{
// create a binary mask around the "hotspot" region, fill the shape of a sphere around our hotspot center
// typename DuplicatorType::Pointer copyMachine = DuplicatorType::New();
// copyMachine->SetInputImage(inputImage);
// copyMachine->Update();
// typename CastFilterType::Pointer caster = CastFilterType::New();
// caster->SetInput( copyMachine->GetOutput() );
// caster->Update();
typename MaskImageType::Pointer hotspotMaskITK = MaskImageType::New();
hotspotMaskITK->SetOrigin(inputImage->GetOrigin());
hotspotMaskITK->SetSpacing(inputImage->GetSpacing());
hotspotMaskITK->SetLargestPossibleRegion(inputImage->GetLargestPossibleRegion());
hotspotMaskITK->SetBufferedRegion(inputImage->GetBufferedRegion());
hotspotMaskITK->SetDirection(inputImage->GetDirection());
hotspotMaskITK->SetNumberOfComponentsPerPixel(inputImage->GetNumberOfComponentsPerPixel());
hotspotMaskITK->Allocate();
hotspotMaskITK->FillBuffer(1);
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);
FillHotspotMaskPixels(hotspotMaskITK.GetPointer(), maskCenter, m_HotspotRadiusinMM);
//obtain mitk::Image::Pointer from itk::Image
mitk::Image::Pointer hotspotMaskAsMITKImage = mitk::GrabItkImageMemory(hotspotMaskITK);
m_InternalMask = hotspotMaskAsMITKImage;
m_ConvolutionImageMaxIndex = convolutionImageInformation.MaxIndex;
m_ConvolutionImageMinIndex = convolutionImageInformation.MinIndex;
}
}
bool HotspotMaskGenerator::IsUpdateRequired() const
{
unsigned long thisClassTimeStamp = this->GetMTime();
unsigned long internalMaskTimeStamp = m_InternalMask->GetMTime();
unsigned long maskGeneratorTimeStamp = m_Mask->GetMTime();
unsigned long inputImageTimeStamp = m_inputImage->GetMTime();
if (thisClassTimeStamp > m_InternalMaskUpdateTime) // inputs have changed
{
return true;
}
if (m_InternalMaskUpdateTime < maskGeneratorTimeStamp || m_InternalMaskUpdateTime < inputImageTimeStamp) // mask image has changed outside of this class
{
return true;
}
if (internalMaskTimeStamp > m_InternalMaskUpdateTime) // internal mask has been changed outside of this class
{
return true;
}
return false;
}
}
diff --git a/Modules/ImageStatistics/mitkHotspotMaskGenerator.h b/Modules/ImageStatistics/mitkHotspotMaskGenerator.h
index 909fb6b1a6..baf9ab4059 100644
--- a/Modules/ImageStatistics/mitkHotspotMaskGenerator.h
+++ b/Modules/ImageStatistics/mitkHotspotMaskGenerator.h
@@ -1,164 +1,170 @@
#ifndef MITKHOTSPOTCALCULATOR_H
#define MITKHOTSPOTCALCULATOR_H
#include <itkObject.h>
#include <mitkImage.h>
#include <itkImage.h>
#include <itkTimeStamp.h>
#include <stdexcept>
#include <MitkImageStatisticsExports.h>
#include <mitkImageTimeSelector.h>
#include <mitkMaskGenerator.h>
namespace mitk
{
/**
* @brief The HotspotMaskGenerator class is used when a hotspot has to be found in an image. A hotspot is
* the region of the image where the mean intensity is maximal (=brightest spot). It is usually used in PET scans.
* The identification of the hotspot is done as follows: First a cubic (or circular, if image is 2d)
* mask of predefined size is generated. This mask is then convolved with the input image (in fourier domain).
* The maximum value of the convolved image then corresponds to the hotspot.
* If a maskGenerator is set, only the pixels of the convolved image where the corresponding mask is == @a label
* are searched for the maximum value.
*/
class MITKIMAGESTATISTICS_EXPORT HotspotMaskGenerator: public MaskGenerator
{
public:
/** Standard Self typedef */
typedef HotspotMaskGenerator Self;
typedef MaskGenerator Superclass;
typedef itk::SmartPointer< Self > Pointer;
typedef itk::SmartPointer< const Self > ConstPointer;
/** Method for creation through the object factory. */
itkNewMacro(Self)
/** Runtime information support. */
itkTypeMacro(HotspotMaskGenerator, MaskGenerator)
/**
@brief Set the input image. Required for this class
*/
void SetInputImage(mitk::Image::Pointer inputImage);
/**
@brief Set a mask (can be nullptr if no mask is desired)
*/
void SetMask(MaskGenerator::Pointer mask);
/**
@brief Set the radius of the hotspot (in MM)
*/
void SetHotspotRadiusInMM(double radiusInMillimeter);
const double& GetHotspotRadiusinMM() const;
/**
@brief Define whether the hotspot must be completely inside the image. Default is true
*/
void SetHotspotMustBeCompletelyInsideImage(bool hotspotCompletelyInsideImage);
bool GetHotspotMustBeCompletelyInsideImage() const;
/**
@brief If a maskGenerator is set, this detemines which mask value is used
*/
void SetLabel(unsigned short label);
/**
@brief Computes and returns the hotspot mask. The hotspot mask has the same size as the input image. The hopspot has value 1, the remaining pixels are set to 0
*/
mitk::Image::Pointer GetMask();
/**
@brief Returns the image index where the hotspot is located
*/
vnl_vector<int> GetHotspotIndex();
/**
@brief Returns the index where the convolution image is minimal (darkest spot in image)
*/
vnl_vector<int> GetConvolutionImageMinIndex();
+ /**
+ * @brief SetTimeStep is used to set the time step for which the mask is to be generated
+ * @param timeStep
+ */
+ void SetTimeStep(unsigned int timeStep);
+
protected:
HotspotMaskGenerator();
~HotspotMaskGenerator();
class ImageExtrema
{
public:
bool Defined;
double Max;
double Min;
vnl_vector<int> MaxIndex;
vnl_vector<int> MinIndex;
ImageExtrema()
:Defined(false)
,Max(itk::NumericTraits<double>::min())
,Min(itk::NumericTraits<double>::max())
{
}
};
private:
/** \brief Returns size of convolution kernel depending on spacing and radius. */
template <unsigned int VImageDimension>
itk::Size<VImageDimension>
CalculateConvolutionKernelSize(double spacing[VImageDimension], double radiusInMM);
/** \brief Generates image of kernel which is needed for convolution. */
template <unsigned int VImageDimension>
itk::SmartPointer< itk::Image<float, VImageDimension> >
GenerateHotspotSearchConvolutionKernel(double spacing[VImageDimension], double radiusInMM);
/** \brief Convolves image with spherical kernel image. Used for hotspot calculation. */
template <typename TPixel, unsigned int VImageDimension>
itk::SmartPointer< itk::Image<TPixel, VImageDimension> >
GenerateConvolutionImage( const itk::Image<TPixel, VImageDimension>* inputImage );
/** \brief Fills pixels of the spherical hotspot mask. */
template < typename TPixel, unsigned int VImageDimension>
void
FillHotspotMaskPixels( itk::Image<TPixel, VImageDimension>* maskImage,
itk::Point<double, VImageDimension> sphereCenter,
double sphereRadiusInMM);
/** \brief */
template <typename TPixel, unsigned int VImageDimension>
void
CalculateHotspotMask(itk::Image<TPixel, VImageDimension>* inputImage,
typename itk::Image<unsigned short, VImageDimension>::Pointer maskImage,
unsigned int label);
template <typename TPixel, unsigned int VImageDimension >
ImageExtrema CalculateExtremaWorld( const itk::Image<TPixel, VImageDimension>* inputImage,
typename itk::Image<unsigned short, VImageDimension>::Pointer maskImage,
double neccessaryDistanceToImageBorderInMM,
unsigned int label);
bool IsUpdateRequired() const;
HotspotMaskGenerator(const HotspotMaskGenerator &);
HotspotMaskGenerator & operator=(const HotspotMaskGenerator &);
MaskGenerator::Pointer m_Mask;
mitk::Image::Pointer m_internalImage;
itk::Image<unsigned short, 2>::Pointer m_internalMask2D;
itk::Image<unsigned short, 3>::Pointer m_internalMask3D;
double m_HotspotRadiusinMM;
bool m_HotspotMustBeCompletelyInsideImage;
bool m_HotspotParamsChanged;
unsigned short m_Label;
vnl_vector<int> m_ConvolutionImageMinIndex, m_ConvolutionImageMaxIndex;
unsigned long m_InternalMaskUpdateTime;
};
}
#endif // MITKHOTSPOTCALCULATOR
diff --git a/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.cpp b/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.cpp
index ba068667dd..4e41949751 100644
--- a/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.cpp
+++ b/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.cpp
@@ -1,112 +1,120 @@
#include <mitkIgnorePixelMaskGenerator.h>
#include <mitkImageTimeSelector.h>
#include <mitkImageAccessByItk.h>
#include <itkImageIterator.h>
#include <itkImageConstIterator.h>
#include <mitkITKImageImport.h>
namespace mitk
{
void IgnorePixelMaskGenerator::SetIgnoredPixelValue(RealType pixelValue)
{
if (pixelValue != m_IgnoredPixelValue)
{
m_IgnoredPixelValue = pixelValue;
this->Modified();
}
}
+void IgnorePixelMaskGenerator::SetTimeStep(unsigned int timeStep)
+{
+ if (m_TimeStep != timeStep)
+ {
+ m_TimeStep = timeStep;
+ }
+}
+
mitk::Image::Pointer IgnorePixelMaskGenerator::GetMask()
{
if (IsUpdateRequired())
{
if (m_inputImage.IsNull())
{
MITK_ERROR << "Image not set!";
}
if (m_IgnoredPixelValue == std::numeric_limits<RealType>::min())
{
MITK_ERROR << "IgnotePixelValue not set!";
}
if (m_TimeStep > (m_inputImage->GetTimeSteps() - 1))
{
MITK_ERROR << "Invalid time step: " << m_TimeStep << ". The image has " << m_inputImage->GetTimeSteps() << " timeSteps!";
}
// extractimage time slice
ImageTimeSelector::Pointer imgTimeSel = ImageTimeSelector::New();
imgTimeSel->SetInput(m_inputImage);
imgTimeSel->SetTimeNr(m_TimeStep);
imgTimeSel->UpdateLargestPossibleRegion();
mitk::Image::Pointer timeSliceImage = imgTimeSel->GetOutput();
// update m_InternalMask
AccessByItk(timeSliceImage, InternalCalculateMask);
m_InternalMask->SetGeometry(timeSliceImage->GetGeometry());
this->Modified();
}
m_InternalMaskUpdateTime = m_InternalMask->GetMTime();
return m_InternalMask;
}
template <typename TPixel, unsigned int VImageDimension>
void IgnorePixelMaskGenerator::InternalCalculateMask(typename itk::Image<TPixel, VImageDimension>* image)
{
typedef itk::Image<TPixel, VImageDimension> ImageType;
typedef itk::Image<unsigned short, VImageDimension> MaskType;
typename MaskType::Pointer mask = MaskType::New();
mask->SetOrigin(image->GetOrigin());
mask->SetSpacing(image->GetSpacing());
mask->SetLargestPossibleRegion(image->GetLargestPossibleRegion());
mask->SetBufferedRegion(image->GetBufferedRegion());
mask->SetDirection(image->GetDirection());
mask->SetNumberOfComponentsPerPixel(image->GetNumberOfComponentsPerPixel());
mask->Allocate();
mask->FillBuffer(1);
// iterate over image and mask and set mask=1 if image=m_IgnorePixelValue
itk::ImageRegionConstIterator<ImageType> imageIterator(image, image->GetLargestPossibleRegion());
itk::ImageRegionIterator<MaskType> maskIterator(mask, mask->GetLargestPossibleRegion());
for (imageIterator.GoToBegin(); !imageIterator.IsAtEnd(); ++imageIterator, ++maskIterator)
{
if (imageIterator.Value() == static_cast<TPixel>(m_IgnoredPixelValue))
{
maskIterator.Set(0);
}
}
m_InternalMask = GrabItkImageMemory(mask);
}
bool IgnorePixelMaskGenerator::IsUpdateRequired() const
{
unsigned long thisClassTimeStamp = this->GetMTime();
unsigned long internalMaskTimeStamp = m_InternalMask->GetMTime();
unsigned long inputImageTimeStamp = m_inputImage->GetMTime();
if (thisClassTimeStamp > m_InternalMaskUpdateTime) // inputs have changed
{
return true;
}
if (m_InternalMaskUpdateTime < inputImageTimeStamp) // mask image has changed outside of this class
{
return true;
}
if (internalMaskTimeStamp > m_InternalMaskUpdateTime) // internal mask has been changed outside of this class
{
return true;
}
return false;
}
} // end namespace
diff --git a/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.h b/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.h
index d8eecde982..21887446af 100644
--- a/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.h
+++ b/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.h
@@ -1,67 +1,73 @@
#ifndef MITKIGNOREPIXELMASKGEN_
#define MITKIGNOREPIXELMASKGEN_
#include <mitkImage.h>
#include <MitkImageStatisticsExports.h>
#include <mitkMaskGenerator.h>
#include <limits>
#include <itkImage.h>
namespace mitk
{
/**
* @brief The IgnorePixelMaskGenerator class is used to generate a mask that is zero for specific pixel values in the input image. This class requires an input image.
*/
class MITKIMAGESTATISTICS_EXPORT IgnorePixelMaskGenerator: public MaskGenerator
{
public:
/** Standard Self typedef */
typedef IgnorePixelMaskGenerator Self;
typedef MaskGenerator Superclass;
typedef itk::SmartPointer< Self > Pointer;
typedef itk::SmartPointer< const Self > ConstPointer;
typedef double RealType;
/** Method for creation through the object factory. */
itkNewMacro(Self)
/** Runtime information support. */
itkTypeMacro(IgnorePixelMaskGenerator, MaskGenerator)
/**
* @brief The mask will be 0 there inputImage==pixelValue and 1 otherwise
*/
void SetIgnoredPixelValue(RealType pixelValue);
/**
* @brief Computes and returns the mask
*/
mitk::Image::Pointer GetMask();
+ /**
+ * @brief SetTimeStep is used to set the time step for which the mask is to be generated
+ * @param timeStep
+ */
+ void SetTimeStep(unsigned int timeStep);
+
protected:
IgnorePixelMaskGenerator():
m_IgnoredPixelValue(std::numeric_limits<RealType>::min())
{
m_TimeStep = 0;
m_InternalMaskUpdateTime = 0;
m_InternalMask = mitk::Image::New();
}
~IgnorePixelMaskGenerator(){}
template <typename TPixel, unsigned int VImageDimension>
void InternalCalculateMask(typename itk::Image<TPixel, VImageDimension>* image);
private:
bool IsUpdateRequired() const;
RealType m_IgnoredPixelValue;
unsigned long m_InternalMaskUpdateTime;
};
}
#endif
diff --git a/Modules/ImageStatistics/mitkImageMaskGenerator.cpp b/Modules/ImageStatistics/mitkImageMaskGenerator.cpp
index a96916aebb..d1ebd9366e 100644
--- a/Modules/ImageStatistics/mitkImageMaskGenerator.cpp
+++ b/Modules/ImageStatistics/mitkImageMaskGenerator.cpp
@@ -1,82 +1,88 @@
#include <mitkImageMaskGenerator.h>
#include <mitkImageTimeSelector.h>
#include <stdexcept>
namespace mitk {
void ImageMaskGenerator::SetImageMask(Image::Pointer maskImage)
{
if (m_internalMaskImage != maskImage)
{
m_internalMaskImage = maskImage;
this->Modified();
}
}
+void ImageMaskGenerator::SetTimeStep(unsigned int timeStep)
+{
+ if (timeStep != m_TimeStep)
+ {
+ m_TimeStep = timeStep;
+ UpdateInternalMask();
+ }
+}
+
+void ImageMaskGenerator::UpdateInternalMask()
+{
+ unsigned int timeStepForExtraction;
+
+ if (m_TimeStep >= m_internalMaskImage->GetTimeSteps())
+ {
+ MITK_WARN << "Warning: time step > number of time steps in mask image, using last time step";
+ timeStepForExtraction = m_internalMaskImage->GetTimeSteps() - 1;
+ }
+ else
+ {
+ timeStepForExtraction = m_TimeStep;
+ }
+ ImageTimeSelector::Pointer imageTimeSelector = ImageTimeSelector::New();
+ imageTimeSelector->SetInput(m_internalMaskImage);
+ imageTimeSelector->SetTimeNr(timeStepForExtraction);
+ imageTimeSelector->UpdateLargestPossibleRegion();
+
+ m_InternalMask = mitk::Image::New();
+ m_InternalMask = imageTimeSelector->GetOutput();
+}
+
+
mitk::Image::Pointer ImageMaskGenerator::GetMask()
{
if (m_internalMaskImage.IsNull())
{
MITK_ERROR << "Mask Image is nullptr";
}
-
- if (this->IsUpdateRequired())
+ if (IsUpdateRequired())
{
- unsigned int timeStepForExtraction;
-
- if (m_TimeStep >= m_internalMaskImage->GetTimeSteps())
- {
- MITK_WARN << "Warning: time step > number of time steps in mask image, using last time step";
- timeStepForExtraction = m_internalMaskImage->GetTimeSteps() - 1;
- }
- else
- {
- timeStepForExtraction = m_TimeStep;
- }
- ImageTimeSelector::Pointer imageTimeSelector = ImageTimeSelector::New();
- imageTimeSelector->SetInput(m_internalMaskImage);
- imageTimeSelector->SetTimeNr(timeStepForExtraction);
- imageTimeSelector->UpdateLargestPossibleRegion();
-
- m_InternalMask = mitk::Image::New();
- m_InternalMask = imageTimeSelector->GetOutput();
- this->Modified();
+ UpdateInternalMask();
}
- m_InternalMaskUpdateTime = m_InternalMask->GetMTime();
return m_InternalMask;
}
bool ImageMaskGenerator::IsUpdateRequired() const
{
- unsigned long thisClassTimeStamp = this->GetMTime();
unsigned long internalMaskTimeStamp = m_InternalMask->GetMTime();
unsigned long maskImageTimeStamp = m_internalMaskImage->GetMTime();
- if (thisClassTimeStamp > m_InternalMaskUpdateTime) // inputs have changed
- {
- return true;
- }
-
- if (m_InternalMaskUpdateTime < maskImageTimeStamp) // mask image has changed outside of this class
+ if (maskImageTimeStamp > internalMaskTimeStamp) // inputs have changed
{
return true;
}
- if (internalMaskTimeStamp > m_InternalMaskUpdateTime) // internal mask has been changed outside of this class
+ if (this->GetMTime() > maskImageTimeStamp) // input has changed
{
return true;
}
return false;
}
}
diff --git a/Modules/ImageStatistics/mitkImageMaskGenerator.h b/Modules/ImageStatistics/mitkImageMaskGenerator.h
index 0d6d12c826..bf8918fb9d 100644
--- a/Modules/ImageStatistics/mitkImageMaskGenerator.h
+++ b/Modules/ImageStatistics/mitkImageMaskGenerator.h
@@ -1,48 +1,51 @@
#ifndef mitkBinaryMaskGenerator
#define mitkBinaryMaskGenerator
#include <mitkImage.h>
#include <MitkImageStatisticsExports.h>
#include <mitkMaskGenerator.h>
#include <itkObject.h>
#include <itkSmartPointer.h>
namespace mitk
{
class MITKIMAGESTATISTICS_EXPORT ImageMaskGenerator: public MaskGenerator
{
public:
/** Standard Self typedef */
typedef ImageMaskGenerator Self;
typedef MaskGenerator Superclass;
typedef itk::SmartPointer< Self > Pointer;
typedef itk::SmartPointer< const Self > ConstPointer;
/** Method for creation through the object factory. */
itkNewMacro(Self)
/** Runtime information support. */
itkTypeMacro(BinaryImageMaskGenerator, MaskGenerator)
mitk::Image::Pointer GetMask();
+ void SetTimeStep(unsigned int timeStep);
+
void SetImageMask(mitk::Image::Pointer maskImage);
protected:
ImageMaskGenerator():Superclass(){
m_InternalMaskUpdateTime = 0;
m_InternalMask = mitk::Image::New();
}
private:
bool IsUpdateRequired() const;
+ void UpdateInternalMask();
mitk::Image::Pointer m_internalMaskImage;
unsigned long m_InternalMaskUpdateTime;
};
}
#endif
diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp
index 4551c82546..04f71a6dcc 100644
--- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp
+++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp
@@ -1,642 +1,642 @@
#include <limits>
#include <math.h>
#include <itkImageToHistogramFilter.h>
#include <itkMaskedImageToHistogramFilter.h>
#include <itkHistogram.h>
#include <itkExtractImageFilter.h>
#include <itkImageRegionConstIterator.h>
#include <itkChangeInformationImageFilter.h>
#include <itkMinimumMaximumImageCalculator.h>
#include <itkMinimumMaximumImageFilter.h>
#include <itkMaskImageFilter.h>
#include <itkExceptionObject.h>
#include <mitkImageStatisticsCalculator.h>
#include <mitkImage.h>
#include <mitkHistogramStatisticsCalculator.h>
#include <mitkImageAccessByItk.h>
#include <mitkImageToItk.h>
#include <mitkExtendedStatisticsImageFilter.h>
#include <mitkExtendedLabelStatisticsImageFilter.h>
#include <mitkImageTimeSelector.h>
#include <mitkMinMaxImageFilterWithIndex.h>
#include <mitkMinMaxLabelmageFilterWithIndex.h>
#include <mitkitkMaskImageFilter.h>
#include <mitkImageCast.h>
#include <mitkMaskUtilities.h>
#include "itkImageFileWriter.h"
namespace mitk
{
void ImageStatisticsCalculator::SetInputImage(mitk::Image::Pointer image)
{
if (image != m_Image)
{
m_Image = image;
m_StatisticsByTimeStep.resize(m_Image->GetTimeSteps());
m_StatisticsUpdateTimePerTimeStep.resize(m_Image->GetTimeSteps());
std::fill(m_StatisticsUpdateTimePerTimeStep.begin(), m_StatisticsUpdateTimePerTimeStep.end(), 0);
this->Modified();
}
}
void ImageStatisticsCalculator::SetMask(mitk::MaskGenerator::Pointer mask)
{
if (mask != m_MaskGenerator)
{
m_MaskGenerator = mask;
this->Modified();
}
}
void ImageStatisticsCalculator::SetSecondaryMask(mitk::MaskGenerator::Pointer mask)
{
if (mask != m_SecondaryMaskGenerator)
{
m_SecondaryMaskGenerator = mask;
this->Modified();
}
}
void ImageStatisticsCalculator::SetNBinsForHistogramStatistics(unsigned int nBins)
{
if (nBins != m_nBinsForHistogramStatistics)
{
m_nBinsForHistogramStatistics = nBins;
this->Modified();
this->m_UseBinSizeOverNBins = false;
}
if (m_UseBinSizeOverNBins)
{
this->Modified();
this->m_UseBinSizeOverNBins = false;
}
}
unsigned int ImageStatisticsCalculator::GetNBinsForHistogramStatistics() const
{
return m_nBinsForHistogramStatistics;
}
void ImageStatisticsCalculator::SetBinSizeForHistogramStatistics(double binSize)
{
if (binSize != m_binSizeForHistogramStatistics)
{
m_binSizeForHistogramStatistics = binSize;
this->Modified();
this->m_UseBinSizeOverNBins = true;
}
if (!m_UseBinSizeOverNBins)
{
this->Modified();
this->m_UseBinSizeOverNBins = true;
}
}
double ImageStatisticsCalculator::GetBinSizeForHistogramStatistics() const
{
return m_binSizeForHistogramStatistics;
}
ImageStatisticsCalculator::StatisticsContainer::Pointer ImageStatisticsCalculator::GetStatistics(unsigned int timeStep, unsigned int label)
{
if (timeStep >= m_StatisticsByTimeStep.size())
{
mitkThrow() << "invalid timeStep in ImageStatisticsCalculator_v2::GetStatistics";
}
if (m_Image.IsNull())
{
mitkThrow() << "no image";
}
if (!m_Image->IsInitialized())
{
mitkThrow() << "Image not initialized!";
}
if (IsUpdateRequired(timeStep))
{
if (m_MaskGenerator.IsNotNull())
{
m_MaskGenerator->SetTimeStep(timeStep);
m_InternalMask = m_MaskGenerator->GetMask();
if (m_MaskGenerator->GetReferenceImage().IsNotNull())
{
m_InternalImageForStatistics = m_MaskGenerator->GetReferenceImage();
}
else
{
m_InternalImageForStatistics = m_Image;
}
}
else
{
m_InternalImageForStatistics = m_Image;
}
if (m_SecondaryMaskGenerator.IsNotNull())
{
m_SecondaryMaskGenerator->SetTimeStep(timeStep);
m_SecondaryMask = m_SecondaryMaskGenerator->GetMask();
}
ImageTimeSelector::Pointer imgTimeSel = ImageTimeSelector::New();
imgTimeSel->SetInput(m_InternalImageForStatistics);
imgTimeSel->SetTimeNr(timeStep);
imgTimeSel->UpdateLargestPossibleRegion();
m_ImageTimeSlice = imgTimeSel->GetOutput();
// Calculate statistics with/without mask
if (m_MaskGenerator.IsNull() && m_SecondaryMaskGenerator.IsNull())
{
// 1) calculate statistics unmasked:
AccessByItk_1(m_ImageTimeSlice, InternalCalculateStatisticsUnmasked, timeStep)
}
else
{
// 2) calculate statistics masked
AccessByItk_1(m_ImageTimeSlice, InternalCalculateStatisticsMasked, timeStep)
}
- this->Modified();
+ //this->Modified();
}
m_StatisticsUpdateTimePerTimeStep[timeStep] = m_StatisticsByTimeStep[timeStep][m_StatisticsByTimeStep[timeStep].size()-1]->GetMTime();
for (std::vector<StatisticsContainer::Pointer>::iterator it = m_StatisticsByTimeStep[timeStep].begin(); it != m_StatisticsByTimeStep[timeStep].end(); ++it)
{
StatisticsContainer::Pointer statCont = *it;
if (statCont->GetLabel() == label)
{
return statCont->Clone();
}
}
// these lines will ony be executed if the requested label could not be found!
MITK_WARN << "Invalid label: " << label << " in time step: " << timeStep;
return StatisticsContainer::New();
}
template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsUnmasked(
typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep)
{
typedef typename itk::Image< TPixel, VImageDimension > ImageType;
typedef typename itk::ExtendedStatisticsImageFilter<ImageType> ImageStatisticsFilterType;
typedef typename itk::MinMaxImageFilterWithIndex<ImageType> MinMaxFilterType;
StatisticsContainer::Pointer statisticsResult = StatisticsContainer::New();
typename ImageStatisticsFilterType::Pointer statisticsFilter = ImageStatisticsFilterType::New();
statisticsFilter->SetInput(image);
statisticsFilter->SetCoordinateTolerance(0.001);
statisticsFilter->SetDirectionTolerance(0.001);
// TODO: this is single threaded. Implement our own image filter that does this multi threaded
// typename itk::MinimumMaximumImageCalculator<ImageType>::Pointer imgMinMaxFilter = itk::MinimumMaximumImageCalculator<ImageType>::New();
// imgMinMaxFilter->SetImage(image);
// imgMinMaxFilter->Compute();
vnl_vector<int> minIndex, maxIndex;
typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New();
minMaxFilter->SetInput(image);
minMaxFilter->UpdateLargestPossibleRegion();
typename ImageType::PixelType minval = minMaxFilter->GetMin();
typename ImageType::PixelType maxval = minMaxFilter->GetMax();
typename ImageType::IndexType tmpMinIndex = minMaxFilter->GetMinIndex();
typename ImageType::IndexType tmpMaxIndex = minMaxFilter->GetMaxIndex();
// typename ImageType::IndexType tmpMinIndex = imgMinMaxFilter->GetIndexOfMinimum();
// typename ImageType::IndexType tmpMaxIndex = imgMinMaxFilter->GetIndexOfMaximum();
minIndex.set_size(tmpMaxIndex.GetIndexDimension());
maxIndex.set_size(tmpMaxIndex.GetIndexDimension());
for (unsigned int i=0; i < tmpMaxIndex.GetIndexDimension(); i++)
{
minIndex[i] = tmpMinIndex[i];
maxIndex[i] = tmpMaxIndex[i];
}
statisticsResult->SetMinIndex(minIndex);
statisticsResult->SetMaxIndex(maxIndex);
//convert m_binSize in m_nBins if necessary
unsigned int nBinsForHistogram;
if (m_UseBinSizeOverNBins)
{
nBinsForHistogram = std::max(static_cast<double>(std::ceil(maxval - minval)) / m_binSizeForHistogramStatistics, 10.); // do not allow less than 10 bins
}
else
{
nBinsForHistogram = m_nBinsForHistogramStatistics;
}
statisticsFilter->SetHistogramParameters(nBinsForHistogram, minval, maxval);
try
{
statisticsFilter->Update();
}
catch (const itk::ExceptionObject& e)
{
mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what();
}
// no mask, therefore just one label = the whole image
m_StatisticsByTimeStep[timeStep].resize(1);
statisticsResult->SetLabel(1);
statisticsResult->SetN(image->GetLargestPossibleRegion().GetNumberOfPixels());
statisticsResult->SetMean(statisticsFilter->GetMean());
statisticsResult->SetMin(statisticsFilter->GetMinimum());
statisticsResult->SetMax(statisticsFilter->GetMaximum());
statisticsResult->SetVariance(statisticsFilter->GetVariance());
statisticsResult->SetStd(statisticsFilter->GetSigma());
statisticsResult->SetSkewness(statisticsFilter->GetSkewness());
statisticsResult->SetKurtosis(statisticsFilter->GetKurtosis());
statisticsResult->SetRMS(std::sqrt(std::pow(statisticsFilter->GetMean(), 2.) + statisticsFilter->GetVariance())); // variance = sigma^2
statisticsResult->SetMPP(statisticsFilter->GetMPP());
statisticsResult->SetEntropy(statisticsFilter->GetEntropy());
statisticsResult->SetMedian(statisticsFilter->GetMedian());
statisticsResult->SetUniformity(statisticsFilter->GetUniformity());
statisticsResult->SetUPP(statisticsFilter->GetUPP());
statisticsResult->SetHistogram(statisticsFilter->GetHistogram());
m_StatisticsByTimeStep[timeStep][0] = statisticsResult;
}
template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsMasked(
typename itk::Image< TPixel, VImageDimension >* image,
unsigned int timeStep)
{
typedef itk::Image< TPixel, VImageDimension > ImageType;
typedef itk::Image< MaskPixelType, VImageDimension > MaskType;
typedef typename MaskType::PixelType LabelPixelType;
typedef itk::ExtendedLabelStatisticsImageFilter< ImageType, MaskType > ImageStatisticsFilterType;
typedef MaskUtilities< TPixel, VImageDimension > MaskUtilType;
typedef typename itk::MinMaxLabelImageFilterWithIndex<ImageType, MaskType> MinMaxLabelFilterType;
typedef typename ImageType::PixelType InputImgPixelType;
// workaround: if m_SecondaryMaskGenerator ist not null but m_MaskGenerator is! (this is the case if we request a 'ignore zuero valued pixels'
// mask in the gui but do not define a primary mask)
bool swapMasks = false;
if (m_SecondaryMask.IsNotNull() && m_InternalMask.IsNull())
{
m_InternalMask = m_SecondaryMask;
m_SecondaryMask = nullptr;
swapMasks = true;
}
// maskImage has to have the same dimension as image
typename MaskType::Pointer maskImage = MaskType::New();
try {
// try to access the pixel values directly (no copying or casting). Only works if mask pixels are of pixelType unsigned short
maskImage = ImageToItkImage< MaskPixelType, VImageDimension >(m_InternalMask);
}
catch (itk::ExceptionObject & e)
{
// if the pixel type of the mask is not short, then we have to make a copy of m_InternalMask (and cast the values)
CastToItkImage(m_InternalMask, maskImage);
}
// if we have a secondary mask (say a ignoreZeroPixelMask) we need to combine the masks (corresponds to AND)
if (m_SecondaryMask.IsNotNull())
{
// dirty workaround for a bug when pf mask + any other mask is used in conjunction. We need a proper fix for this (Fabian Isensee is responsible and probably working on it!)
if (m_InternalMask->GetDimension() == 2 && (m_SecondaryMask->GetDimension() == 3 || m_SecondaryMask->GetDimension() == 4))
{
mitk::Image::Pointer old_img = m_SecondaryMaskGenerator->GetReferenceImage();
m_SecondaryMaskGenerator->SetInputImage(m_MaskGenerator->GetReferenceImage());
m_SecondaryMask = m_SecondaryMaskGenerator->GetMask();
m_SecondaryMaskGenerator->SetInputImage(old_img);
}
typename MaskType::Pointer secondaryMaskImage = MaskType::New();
secondaryMaskImage = ImageToItkImage< MaskPixelType, VImageDimension >(m_SecondaryMask);
// secondary mask should be a ignore zero value pixel mask derived from image. it has to be cropped to the mask region (which may be planar or simply smaller)
typename MaskUtilities<MaskPixelType, VImageDimension>::Pointer secondaryMaskMaskUtil = MaskUtilities<MaskPixelType, VImageDimension>::New();
secondaryMaskMaskUtil->SetImage(secondaryMaskImage.GetPointer());
secondaryMaskMaskUtil->SetMask(maskImage.GetPointer());
typename MaskType::Pointer adaptedSecondaryMaskImage = secondaryMaskMaskUtil->ExtractMaskImageRegion();
typename itk::MaskImageFilter2<MaskType, MaskType, MaskType>::Pointer maskFilter = itk::MaskImageFilter2<MaskType, MaskType, MaskType>::New();
maskFilter->SetInput1(maskImage);
maskFilter->SetInput2(adaptedSecondaryMaskImage);
maskFilter->SetMaskingValue(1); // all pixels of maskImage where secondaryMaskImage==1 will be kept, all the others are set to 0
maskFilter->UpdateLargestPossibleRegion();
maskImage = maskFilter->GetOutput();
}
typename MaskUtilType::Pointer maskUtil = MaskUtilType::New();
maskUtil->SetImage(image);
maskUtil->SetMask(maskImage.GetPointer());
// if mask is smaller than image, extract the image region where the mask is
typename ImageType::Pointer adaptedImage = ImageType::New();
adaptedImage = maskUtil->ExtractMaskImageRegion(); // this also checks mask sanity
// find min, max, minindex and maxindex
typename MinMaxLabelFilterType::Pointer minMaxFilter = MinMaxLabelFilterType::New();
minMaxFilter->SetInput(adaptedImage);
minMaxFilter->SetLabelInput(maskImage);
minMaxFilter->UpdateLargestPossibleRegion();
// set histogram parameters for each label individually (min/max may be different for each label)
typedef typename std::map<LabelPixelType, InputImgPixelType> MapType;
typedef typename std::pair<LabelPixelType, InputImgPixelType> PairType;
std::vector<LabelPixelType> relevantLabels = minMaxFilter->GetRelevantLabels();
MapType minVals;
MapType maxVals;
std::map<LabelPixelType, unsigned int> nBins;
for (LabelPixelType label:relevantLabels)
{
minVals.insert(PairType(label, minMaxFilter->GetMin(label)));
maxVals.insert(PairType(label, minMaxFilter->GetMax(label)));
unsigned int nBinsForHistogram;
if (m_UseBinSizeOverNBins)
{
nBinsForHistogram = std::max(static_cast<double>(std::ceil(minMaxFilter->GetMax(label) - minMaxFilter->GetMin(label))) / m_binSizeForHistogramStatistics, 10.); // do not allow less than 10 bins
}
else
{
nBinsForHistogram = m_nBinsForHistogramStatistics;
}
nBins.insert(typename std::pair<LabelPixelType, unsigned int>(label, nBinsForHistogram));
}
typename ImageStatisticsFilterType::Pointer imageStatisticsFilter = ImageStatisticsFilterType::New();
imageStatisticsFilter->SetDirectionTolerance(0.001);
imageStatisticsFilter->SetCoordinateTolerance(0.001);
imageStatisticsFilter->SetInput(adaptedImage);
imageStatisticsFilter->SetLabelInput(maskImage);
imageStatisticsFilter->SetHistogramParametersForLabels(nBins, minVals, maxVals);
imageStatisticsFilter->Update();
std::list<int> labels = imageStatisticsFilter->GetRelevantLabels();
std::list<int>::iterator it = labels.begin();
m_StatisticsByTimeStep[timeStep].resize(0);
while(it != labels.end())
{
StatisticsContainer::Pointer statisticsResult = StatisticsContainer::New();
// find min, max, minindex and maxindex
// make sure to only look in the masked region, use a masker for this
vnl_vector<int> minIndex, maxIndex;
mitk::Point3D worldCoordinateMin;
mitk::Point3D worldCoordinateMax;
mitk::Point3D indexCoordinateMin;
mitk::Point3D indexCoordinateMax;
m_InternalImageForStatistics->GetGeometry()->IndexToWorld(minMaxFilter->GetMinIndex(*it), worldCoordinateMin);
m_InternalImageForStatistics->GetGeometry()->IndexToWorld(minMaxFilter->GetMaxIndex(*it), worldCoordinateMax);
m_Image->GetGeometry()->WorldToIndex(worldCoordinateMin, indexCoordinateMin);
m_Image->GetGeometry()->WorldToIndex(worldCoordinateMax, indexCoordinateMax);
//typename ImageType::IndexType tmpMinIndex = minMaxFilter->GetMinIndex(*it);
//typename ImageType::IndexType tmpMaxIndex = minMaxFilter->GetMaxIndex(*it);
//minIndex.set_size(tmpMaxIndex.GetIndexDimension());
//maxIndex.set_size(tmpMaxIndex.GetIndexDimension());
minIndex.set_size(3);
maxIndex.set_size(3);
//for (unsigned int i=0; i < tmpMaxIndex.GetIndexDimension(); i++)
for (unsigned int i=0; i < 3; i++)
{
//minIndex[i] = tmpMinIndex[i] + (maskImage->GetOrigin()[i] - image->GetOrigin()[i]) / (double) maskImage->GetSpacing()[i];
//maxIndex[i] = tmpMaxIndex[i] + (maskImage->GetOrigin()[i] - image->GetOrigin()[i]) / (double) maskImage->GetSpacing()[i];
minIndex[i] = indexCoordinateMin[i];
maxIndex[i] = indexCoordinateMax[i];
}
statisticsResult->SetMinIndex(minIndex);
statisticsResult->SetMaxIndex(maxIndex);
// just debug
TPixel min_Filter = minMaxFilter->GetMin(*it);
TPixel max_Filter = minMaxFilter->GetMax(*it);
TPixel min_Itk = imageStatisticsFilter->GetMinimum(*it);
TPixel max_Itk = imageStatisticsFilter->GetMaximum(*it);
assert(abs(minMaxFilter->GetMax(*it) - imageStatisticsFilter->GetMaximum(*it)) < mitk::eps);
assert(abs(minMaxFilter->GetMin(*it) - imageStatisticsFilter->GetMinimum(*it)) < mitk::eps);
statisticsResult->SetN(imageStatisticsFilter->GetSum(*it) / (double) imageStatisticsFilter->GetMean(*it));
statisticsResult->SetMean(imageStatisticsFilter->GetMean(*it));
statisticsResult->SetMin(imageStatisticsFilter->GetMinimum(*it));
statisticsResult->SetMax(imageStatisticsFilter->GetMaximum(*it));
statisticsResult->SetVariance(imageStatisticsFilter->GetVariance(*it));
statisticsResult->SetStd(imageStatisticsFilter->GetSigma(*it));
statisticsResult->SetSkewness(imageStatisticsFilter->GetSkewness(*it));
statisticsResult->SetKurtosis(imageStatisticsFilter->GetKurtosis(*it));
statisticsResult->SetRMS(std::sqrt(std::pow(imageStatisticsFilter->GetMean(*it), 2.) + imageStatisticsFilter->GetVariance(*it))); // variance = sigma^2
statisticsResult->SetMPP(imageStatisticsFilter->GetMPP(*it));
statisticsResult->SetLabel(*it);
statisticsResult->SetEntropy(imageStatisticsFilter->GetEntropy(*it));
statisticsResult->SetMedian(imageStatisticsFilter->GetMedian(*it));
statisticsResult->SetUniformity(imageStatisticsFilter->GetUniformity(*it));
statisticsResult->SetUPP(imageStatisticsFilter->GetUPP(*it));
statisticsResult->SetHistogram(imageStatisticsFilter->GetHistogram(*it));
m_StatisticsByTimeStep[timeStep].push_back(statisticsResult);
++it;
}
// swap maskGenerators back
if (swapMasks)
{
m_SecondaryMask = m_InternalMask;
m_InternalMask = nullptr;
}
}
bool ImageStatisticsCalculator::IsUpdateRequired(unsigned int timeStep) const
{
unsigned long thisClassTimeStamp = this->GetMTime();
unsigned long inputImageTimeStamp = m_Image->GetMTime();
unsigned long statisticsTimeStamp = m_StatisticsUpdateTimePerTimeStep[timeStep];
if (thisClassTimeStamp > statisticsTimeStamp) // inputs have changed
{
return true;
}
if (inputImageTimeStamp > statisticsTimeStamp) // image has changed
{
return true;
}
if (m_MaskGenerator.IsNotNull())
{
unsigned long maskGeneratorTimeStamp = m_MaskGenerator->GetMTime();
if (maskGeneratorTimeStamp > statisticsTimeStamp) // there is a mask generator and it has changed
{
return true;
}
}
if (m_SecondaryMaskGenerator.IsNotNull())
{
unsigned long maskGeneratorTimeStamp = m_SecondaryMaskGenerator->GetMTime();
if (maskGeneratorTimeStamp > statisticsTimeStamp) // there is a secondary mask generator and it has changed
{
return true;
}
}
return false;
}
ImageStatisticsCalculator::StatisticsContainer::StatisticsContainer():
m_N(nan("")),
m_Mean(nan("")),
m_Min(nan("")),
m_Max(nan("")),
m_Std(nan("")),
m_Variance(nan("")),
m_Skewness(nan("")),
m_Kurtosis(nan("")),
m_RMS(nan("")),
m_MPP(nan("")),
m_Median(nan("")),
m_Uniformity(nan("")),
m_UPP(nan("")),
m_Entropy(nan(""))
{
m_minIndex.set_size(0);
m_maxIndex.set_size(0);
}
ImageStatisticsCalculator::statisticsMapType ImageStatisticsCalculator::StatisticsContainer::GetStatisticsAsMap()
{
ImageStatisticsCalculator::statisticsMapType statisticsAsMap;
statisticsAsMap["N"] = m_N;
statisticsAsMap["Mean"] = m_Mean;
statisticsAsMap["Min"] = m_Min;
statisticsAsMap["Max"] = m_Max;
statisticsAsMap["StandardDeviation"] = m_Std;
statisticsAsMap["Variance"] = m_Variance;
statisticsAsMap["Skewness"] = m_Skewness;
statisticsAsMap["Kurtosis"] = m_Kurtosis;
statisticsAsMap["RMS"] = m_RMS;
statisticsAsMap["MPP"] = m_MPP;
statisticsAsMap["Median"] = m_Median;
statisticsAsMap["Uniformity"] = m_Uniformity;
statisticsAsMap["UPP"] = m_UPP;
statisticsAsMap["Entropy"] = m_Entropy;
statisticsAsMap["Label"] = m_Label;
return statisticsAsMap;
}
void ImageStatisticsCalculator::StatisticsContainer::Reset()
{
m_N = nan("");
m_Mean = nan("");
m_Min = nan("");
m_Max = nan("");
m_Std = nan("");
m_Variance = nan("");
m_Skewness = nan("");
m_Kurtosis = nan("");
m_RMS = nan("");
m_MPP = nan("");
m_Median = nan("");
m_Uniformity = nan("");
m_UPP = nan("");
m_Entropy = nan("");
m_Histogram = HistogramType::New();
m_minIndex.set_size(0);
m_maxIndex.set_size(0);
m_Label = 0;
}
void ImageStatisticsCalculator::StatisticsContainer::Print()
{
ImageStatisticsCalculator::statisticsMapType statMap = this->GetStatisticsAsMap();
// print all map key value pairs
// const auto& val:statMap
for (auto it = statMap.begin(); it != statMap.end(); ++it)
{
std::cout << it->first << ": " << it->second << std::endl;
}
// print the min and max index
std::cout << "Min Index:" << std::endl;
for (auto it = this->GetMinIndex().begin(); it != this->GetMinIndex().end(); ++it)
{
std::cout << *it << " ";
}
std::cout << std::endl;
// print the min and max index
std::cout << "Max Index:" << std::endl;
for (auto it = this->GetMaxIndex().begin(); it != this->GetMaxIndex().end(); ++it)
{
std::cout << *it << " ";
}
std::cout << std::endl;
}
std::string ImageStatisticsCalculator::StatisticsContainer::GetAsString()
{
std::string res = "";
ImageStatisticsCalculator::statisticsMapType statMap = this->GetStatisticsAsMap();
// print all map key value pairs
// const auto& val:statMap
for (auto it = statMap.begin(); it != statMap.end(); ++it)
{
res += std::string(it->first) + ": " + std::to_string(it->second) + "\n";
}
// print the min and max index
res += "Min Index:" + std::string("\n");
for (auto it = this->GetMinIndex().begin(); it != this->GetMinIndex().end(); it++)
{
res += std::to_string(*it) + std::string(" ");
}
res += "\n";
// print the min and max index
res += "Max Index:" + std::string("\n");
for (auto it = this->GetMaxIndex().begin(); it != this->GetMaxIndex().end(); it++)
{
res += std::to_string(*it) + " ";
}
res += "\n";
return res;
}
}
diff --git a/Modules/ImageStatistics/mitkMaskGenerator.cpp b/Modules/ImageStatistics/mitkMaskGenerator.cpp
index 81c3e3a7aa..37374d440f 100644
--- a/Modules/ImageStatistics/mitkMaskGenerator.cpp
+++ b/Modules/ImageStatistics/mitkMaskGenerator.cpp
@@ -1,51 +1,51 @@
#include <mitkMaskGenerator.h>
namespace mitk
{
-void MaskGenerator::SetTimeStep(unsigned int timeStep)
-{
- if (m_TimeStep != timeStep)
- {
- m_TimeStep = timeStep;
- this->Modified();
- }
-}
MaskGenerator::MaskGenerator():
m_TimeStep(0)
{
m_inputImage = nullptr;
}
mitk::Image::Pointer MaskGenerator::GetMask()
{
return mitk::Image::New();
}
//typename itk::Region<3>::Pointer MaskGenerator::GetImageRegionOfMask(Image::Pointer image)
//{
// if (m_InternalMask.IsNull() || m_Modified)
// {
// MITK_ERROR << "Update MaskGenerator first!";
// }
// mitk::BaseGeometry::Pointer imageGeometry = image->GetGeometry();
// mitk::BaseGeometry::Pointer maskGeometry = m_InternalMask->GetGeometry();
//}
+void MaskGenerator::SetTimeStep(unsigned int timeStep)
+{
+ if (timeStep != m_TimeStep)
+ {
+ m_TimeStep = timeStep;
+ }
+}
+
void MaskGenerator::SetInputImage(mitk::Image::Pointer inputImg)
{
if (inputImg != m_inputImage)
{
m_inputImage = inputImg;
this->Modified();
}
}
mitk::Image::Pointer MaskGenerator::GetReferenceImage()
{
return m_inputImage;
}
}
diff --git a/Modules/ImageStatistics/mitkMaskGenerator.h b/Modules/ImageStatistics/mitkMaskGenerator.h
index 627c7bc4e1..3fced6ef46 100644
--- a/Modules/ImageStatistics/mitkMaskGenerator.h
+++ b/Modules/ImageStatistics/mitkMaskGenerator.h
@@ -1,70 +1,66 @@
#ifndef MITKMASKGENERATOR
#define MITKMASKGENERATOR
#include <MitkImageStatisticsExports.h>
#include <mitkImage.h>
#include <itkRegion.h>
#include <itkObject.h>
#include <itkSmartPointer.h>
namespace mitk
{
/**
* \class MaskGenerator
* \brief Base Class for all Mask Generators. Mask generators are classes that provide functionality for the
* creation of binary (or unsigned short) masks that can be applied to an image. See dervied classes for more
* information.
*/
class MITKIMAGESTATISTICS_EXPORT MaskGenerator: public itk::Object
{
public:
/** Standard Self typedef */
typedef MaskGenerator Self;
typedef itk::Object Superclass;
typedef itk::SmartPointer< Self > Pointer;
typedef itk::SmartPointer< const Self > ConstPointer;
/** Method for creation through the object factory. */
itkNewMacro(Self)
/** Runtime information support. */
itkTypeMacro(MaskGenerator, itk::Object)
//~MaskGenerator();
/**
* @brief GetMask must be overridden by derived classes.
* @return mitk::Image::Pointer of generated mask
*/
virtual mitk::Image::Pointer GetMask();
- /**
- * @brief SetTimeStep is used to set the time step for which the mask is to be generated
- * @param timeStep
- */
- void SetTimeStep(unsigned int timeStep);
-
/**
* @brief GetReferenceImage per default returns the inputImage (as set by SetInputImage). If no input image is set it will return a nullptr.
*/
virtual mitk::Image::Pointer GetReferenceImage();
/**
* @brief SetInputImage is used to set the input image to the mask generator. Some subclasses require an input image, others don't. See the documentation of the specific Mask Generator for more information.
*/
void SetInputImage(mitk::Image::Pointer inputImg);
+ virtual void SetTimeStep(unsigned int timeStep);
+
protected:
MaskGenerator();
unsigned int m_TimeStep;
mitk::Image::Pointer m_InternalMask;
mitk::Image::Pointer m_inputImage;
private:
};
}
#endif // MITKMASKGENERATOR
diff --git a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp
index 57ef871eb6..4f7b19990d 100644
--- a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp
+++ b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp
@@ -1,418 +1,426 @@
#include <mitkPlanarFigureMaskGenerator.h>
#include <mitkBaseGeometry.h>
#include <mitkITKImageImport.h>
#include "mitkImageAccessByItk.h"
#include <mitkExtractImageFilter.h>
#include <mitkConvert2Dto3DImageFilter.h>
#include <mitkImageTimeSelector.h>
#include <mitkIOUtil.h>
#include <itkCastImageFilter.h>
#include <itkVTKImageExport.h>
#include <itkVTKImageImport.h>
#include <itkImageDuplicator.h>
#include <itkExceptionObject.h>
#include <vtkPoints.h>
#include <vtkImageStencil.h>
#include <vtkImageImport.h>
#include <vtkImageExport.h>
#include <vtkLassoStencilSource.h>
#include <vtkSmartPointer.h>
namespace mitk
{
void PlanarFigureMaskGenerator::SetPlanarFigure(mitk::PlanarFigure::Pointer planarFigure)
{
if ( planarFigure.IsNull() )
{
throw std::runtime_error( "Error: planar figure empty!" );
}
if ( !planarFigure->IsClosed() )
{
throw std::runtime_error( "Masking not possible for non-closed figures" );
}
const PlaneGeometry *planarFigurePlaneGeometry = 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!" );
}
if (planarFigure != m_PlanarFigure)
{
this->Modified();
m_PlanarFigure = planarFigure;
}
}
mitk::Image::Pointer PlanarFigureMaskGenerator::GetReferenceImage()
{
if (IsUpdateRequired())
{
this->CalculateMask();
}
return m_ReferenceImage;
}
template < typename TPixel, unsigned int VImageDimension >
void PlanarFigureMaskGenerator::InternalCalculateMaskFromPlanarFigure(
const itk::Image< TPixel, VImageDimension > *image, unsigned int axis )
{
typedef itk::Image< TPixel, VImageDimension > ImageType;
typedef itk::Image< unsigned short, 2 > MaskImage2DType;
typename MaskImage2DType::Pointer maskImage = MaskImage2DType::New();
maskImage->SetOrigin(image->GetOrigin());
maskImage->SetSpacing(image->GetSpacing());
maskImage->SetLargestPossibleRegion(image->GetLargestPossibleRegion());
maskImage->SetBufferedRegion(image->GetBufferedRegion());
maskImage->SetDirection(image->GetDirection());
maskImage->SetNumberOfComponentsPerPixel(image->GetNumberOfComponentsPerPixel());
maskImage->Allocate();
maskImage->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_inputImage->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
// TODO use plane geometry normal to determine that automatically, then check whether the PF is aligned with one of the three 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;
}
// 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
// Fabian: From PlaneGeometry documentation:
// Converts a 2D point given in mm (pt2d_mm) relative to the upper-left corner of the geometry into the corresponding world-coordinate (a 3D point in mm, pt3d_mm).
// To convert a 2D point given in units (e.g., pixels in case of an image) into a 2D point given in mm (as required by this method), use IndexToWorld.
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)
{
// Fabian: same as above
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( maskImage );
// 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_InternalITKImageMask2D = duplicator->GetOutput();
}
bool PlanarFigureMaskGenerator::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;
}
void PlanarFigureMaskGenerator::CalculateMask()
{
if (m_inputImage.IsNull())
{
MITK_ERROR << "Image is not set.";
}
if (m_PlanarFigure.IsNull())
{
MITK_ERROR << "PlanarFigure is not set.";
}
if (m_TimeStep != 0)
{
MITK_WARN << "Multiple TimeSteps are not supported in PlanarFigureMaskGenerator (yet).";
}
const BaseGeometry *imageGeometry = m_inputImage->GetGeometry();
if ( imageGeometry == nullptr )
{
throw std::runtime_error( "Image geometry invalid!" );
}
if (m_inputImage->GetTimeSteps() > 0)
{
mitk::ImageTimeSelector::Pointer imgTimeSel = mitk::ImageTimeSelector::New();
imgTimeSel->SetInput(m_inputImage);
imgTimeSel->SetTimeNr(m_TimeStep);
imgTimeSel->UpdateLargestPossibleRegion();
m_InternalTimeSliceImage = imgTimeSel->GetOutput();
}
else
{
m_InternalTimeSliceImage = m_inputImage;
}
m_InternalITKImageMask2D = nullptr;
const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry();
const PlaneGeometry *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigurePlaneGeometry );
//const BaseGeometry *imageGeometry = m_inputImage->GetGeometry();
// 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
itk::Image< unsigned short, 3 >::IndexType index;
imageGeometry->WorldToIndex( planarFigureGeometry->GetOrigin(), index );
unsigned int slice = index[axis];
// extract image slice which corresponds to the planarFigure and store it in m_InternalImageSlice
mitk::Image::Pointer inputImageSlice = extract2DImageSlice(axis, slice);
//mitk::IOUtil::SaveImage(inputImageSlice, "/home/fabian/inputSliceImage.nrrd");
// Compute mask from PlanarFigure
AccessFixedDimensionByItk_1(inputImageSlice,
InternalCalculateMaskFromPlanarFigure,
2, axis)
//convert itk mask to mitk::Image::Pointer and return it
mitk::Image::Pointer planarFigureMaskImage;
planarFigureMaskImage = mitk::GrabItkImageMemory(m_InternalITKImageMask2D);
//mitk::IOUtil::SaveImage(planarFigureMaskImage, "/home/fabian/planarFigureMaskImage.nrrd");
//Convert2Dto3DImageFilter::Pointer sliceTo3DImageConverter = Convert2Dto3DImageFilter::New();
//sliceTo3DImageConverter->SetInput(planarFigureMaskImage);
//sliceTo3DImageConverter->Update();
//mitk::IOUtil::SaveImage(sliceTo3DImageConverter->GetOutput(), "/home/fabian/3DsliceImage.nrrd");
m_ReferenceImage = inputImageSlice;
//mitk::IOUtil::SaveImage(m_ReferenceImage, "/home/fabian/referenceImage.nrrd");
m_InternalMask = planarFigureMaskImage;
}
+void PlanarFigureMaskGenerator::SetTimeStep(unsigned int timeStep)
+{
+ if (timeStep != m_TimeStep)
+ {
+ m_TimeStep = timeStep;
+ }
+}
+
mitk::Image::Pointer PlanarFigureMaskGenerator::GetMask()
{
if (IsUpdateRequired())
{
this->CalculateMask();
this->Modified();
}
- m_InternalMaskUpdateTime = m_InternalMask->GetMTime();
+ m_InternalMaskUpdateTime = this->GetMTime();
return m_InternalMask;
}
mitk::Image::Pointer PlanarFigureMaskGenerator::extract2DImageSlice(unsigned int axis, unsigned int slice)
{
// Extract slice with given position and direction from image
unsigned int dimension = m_InternalTimeSliceImage->GetDimension();
mitk::Image::Pointer imageSlice = mitk::Image::New();
if (dimension == 3)
{
ExtractImageFilter::Pointer imageExtractor = ExtractImageFilter::New();
imageExtractor->SetInput( m_InternalTimeSliceImage );
imageExtractor->SetSliceDimension( axis );
imageExtractor->SetSliceIndex( slice );
imageExtractor->Update();
imageSlice = imageExtractor->GetOutput();
}
else if(dimension == 2)
{
imageSlice = m_InternalTimeSliceImage;
}
else
{
MITK_ERROR << "Unsupported image dimension. Dimension is: " << dimension << ". Only 2D and 3D images are supported.";
}
return imageSlice;
}
bool PlanarFigureMaskGenerator::IsUpdateRequired() const
{
unsigned long thisClassTimeStamp = this->GetMTime();
unsigned long internalMaskTimeStamp = m_InternalMask->GetMTime();
unsigned long planarFigureTimeStamp = m_PlanarFigure->GetMTime();
unsigned long inputImageTimeStamp = m_inputImage->GetMTime();
if (thisClassTimeStamp > m_InternalMaskUpdateTime) // inputs have changed
{
return true;
}
if (m_InternalMaskUpdateTime < planarFigureTimeStamp || m_InternalMaskUpdateTime < inputImageTimeStamp) // mask image has changed outside of this class
{
return true;
}
if (internalMaskTimeStamp > m_InternalMaskUpdateTime) // internal mask has been changed outside of this class
{
return true;
}
return false;
}
}
diff --git a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h
index 8585238d56..990e511833 100644
--- a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h
+++ b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h
@@ -1,119 +1,127 @@
#ifndef MITKPLANARFIGUREMASKGENERATOR
#define MITKPLANARFIGUREMASKGENERATOR
#include <MitkImageStatisticsExports.h>
#include <mitkImage.h>
#include <mitkPlanarFigure.h>
#include <itkImage.h>
#include <mitkMaskGenerator.h>
#include <vtkSmartPointer.h>
#include <itkVTKImageExport.h>
#include <itkVTKImageImport.h>
#include <vtkImageImport.h>
#include <vtkImageExport.h>
namespace mitk
{
/**
* \class PlanarFigureMaskGenerator
* \brief Derived from MaskGenerator. This class is used to convert a mitk::PlanarFigure into a binary image mask
*/
class MITKIMAGESTATISTICS_EXPORT PlanarFigureMaskGenerator: public MaskGenerator
{
public:
/** Standard Self typedef */
typedef PlanarFigureMaskGenerator Self;
typedef MaskGenerator Superclass;
typedef itk::SmartPointer< Self > Pointer;
typedef itk::SmartPointer< const Self > ConstPointer;
/** Method for creation through the object factory. */
itkNewMacro(Self)
/** Runtime information support. */
itkTypeMacro(PlanarFigureMaskGenerator, MaskGenerator)
/**
* @brief GetMask Computes and returns the mask
* @return mitk::Image::Pointer of the generated mask
*/
mitk::Image::Pointer GetMask();
void SetPlanarFigure(mitk::PlanarFigure::Pointer planarFigure);
mitk::Image::Pointer GetReferenceImage();
+ /**
+ * @brief SetTimeStep is used to set the time step for which the mask is to be generated
+ * @param timeStep
+ */
+ void SetTimeStep(unsigned int timeStep);
+
+
protected:
PlanarFigureMaskGenerator():Superclass(){
m_InternalMaskUpdateTime = 0;
m_InternalMask = mitk::Image::New();
m_ReferenceImage = nullptr;
}
+
private:
void CalculateMask();
template < typename TPixel, unsigned int VImageDimension >
void InternalCalculateMaskFromPlanarFigure(
const itk::Image< TPixel, VImageDimension > *image, unsigned int axis );
mitk::Image::Pointer extract2DImageSlice(unsigned int axis, unsigned int slice);
bool GetPrincipalAxis(const BaseGeometry *geometry, Vector3D vector,
unsigned int &axis );
/** Connection from ITK to VTK */
template <typename ITK_Exporter, typename VTK_Importer>
void ConnectPipelines(ITK_Exporter exporter, vtkSmartPointer<VTK_Importer> importer)
{
importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback());
importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback());
importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback());
importer->SetSpacingCallback(exporter->GetSpacingCallback());
importer->SetOriginCallback(exporter->GetOriginCallback());
importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback());
importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback());
importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback());
importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback());
importer->SetDataExtentCallback(exporter->GetDataExtentCallback());
importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback());
importer->SetCallbackUserData(exporter->GetCallbackUserData());
}
/** Connection from VTK to ITK */
template <typename VTK_Exporter, typename ITK_Importer>
void ConnectPipelines(vtkSmartPointer<VTK_Exporter> exporter, ITK_Importer importer)
{
importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback());
importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback());
importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback());
importer->SetSpacingCallback(exporter->GetSpacingCallback());
importer->SetOriginCallback(exporter->GetOriginCallback());
importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback());
importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback());
importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback());
importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback());
importer->SetDataExtentCallback(exporter->GetDataExtentCallback());
importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback());
importer->SetCallbackUserData(exporter->GetCallbackUserData());
}
bool IsUpdateRequired() const;
mitk::PlanarFigure::Pointer m_PlanarFigure;
itk::Image<unsigned short, 2>::Pointer m_InternalITKImageMask2D;
mitk::Image::Pointer m_InternalTimeSliceImage;
mitk::Image::Pointer m_ReferenceImage;
unsigned int m_PlanarFigureAxis;
unsigned long m_InternalMaskUpdateTime;
};
}
#endif // MITKPLANARFIGUREMASKGENERATOR