diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index c27ad0ae8c..95c23977e8 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,2010 +1,2016 @@ /*=================================================================== 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 #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include 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), + m_HistogramBinSize(1.0), m_UseDefaultBinSize(true), 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) : NULL) { Reset(); } ImageStatisticsCalculator::Statistics::Statistics(const Statistics& other) :m_HotspotStatistics( NULL) { 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->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 != NULL; } 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 ); vnl_vector 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 != NULL) { 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() ); delete this->m_HotspotStatistics; this->m_HotspotStatistics = NULL; 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(unsigned int size) +void ImageStatisticsCalculator::SetHistogramBinSize(double size) { this->m_HistogramBinSize = size; } unsigned int 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; } 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 = NULL; } // 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; } const ImageStatisticsCalculator::HistogramType * ImageStatisticsCalculator::GetHistogram( unsigned int timeStep, unsigned int label ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return NULL; } 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 = NULL; m_InternalImageMask3D = NULL; 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 = NULL; 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 == NULL ) { throw std::runtime_error( "Image geometry invalid!" ); } const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); if ( planarFigurePlaneGeometry == NULL ) { throw std::runtime_error( "Planar-Figure not yet initialized!" ); } const PlaneGeometry *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigurePlaneGeometry ); if ( planarFigureGeometry == NULL ) { 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::StatisticsImageFilter< ImageType > StatisticsFilterType; typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); statisticsFilter->SetInput( image ); unsigned long observerTag = statisticsFilter->AddObserver( itk::ProgressEvent(), progressListener ); statisticsFilter->Update(); 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(0.0); statistics.SetVariance(statisticsFilter->GetVariance()); 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 tmpMaxIndex; vnl_vector tmpMinIndex; tmpMaxIndex.set_size(image->GetImageDimension() ); tmpMinIndex.set_size(image->GetImageDimension() ); for (unsigned int i=0; iGetIndexOfMaximum()[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 unsigned int numberOfBins = 200; if (m_UseDefaultBinSize) - m_HistogramBinSize = std::ceil( (statistics.GetMax() - statistics.GetMin() + 1)/numberOfBins ); + m_HistogramBinSize = std::ceil( (statistics.GetMax() - statistics.GetMin() + 1)/numberOfBins ); else - numberOfBins = std::floor( ( (statistics.GetMax() - statistics.GetMin() + 1) / m_HistogramBinSize) + 0.5 ); + 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 itmask(maskImage, maskImage->GetLargestPossibleRegion()); itk::ImageRegionConstIterator 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 itk::LabelStatisticsImageFilter< 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 == NULL ) { 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 ) { 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 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 ) { itkExceptionMacro( << "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->Update(); typename MaskImageType::Pointer adaptedMaskImage = adaptMaskFilter->GetOutput(); // Make sure that mask region is contained within image region if ( !image->GetLargestPossibleRegion().IsInside( adaptedMaskImage->GetLargestPossibleRegion() ) ) { itkExceptionMacro( << "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->Update(); adaptedImage = extractImageFilter->GetOutput(); } else { adaptedImage = image; } // Initialize Filter typedef itk::StatisticsImageFilter< ImageType > StatisticsFilterType; typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); statisticsFilter->SetInput( adaptedImage ); statisticsFilter->Update(); - int numberOfBins = std::floor( ( (statisticsFilter->GetMaximum() - statisticsFilter->GetMinimum() + 1) / m_HistogramBinSize) + 0.5 ); + int numberOfBins = calcNumberOfBins(statisticsFilter->GetMinimum(), statisticsFilter->GetMaximum()); typename LabelStatisticsFilterType::Pointer labelStatisticsFilter; labelStatisticsFilter = LabelStatisticsFilterType::New(); labelStatisticsFilter->SetInput( adaptedImage ); labelStatisticsFilter->SetLabelInput( adaptedMaskImage ); labelStatisticsFilter->UseHistogramsOn(); labelStatisticsFilter->SetHistogramParameters( numberOfBins, 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 labelStatisticsFilter->Update(); this->InvokeEvent( itk::EndEvent() ); labelStatisticsFilter->RemoveObserver( observerTag ); // Find all relevant labels of mask (other than 0) std::list< int > relevantLabels; bool maskNonEmpty = false; unsigned int i; for ( i = 1; i < 4096; ++i ) { if ( labelStatisticsFilter->HasLabel( i ) ) { relevantLabels.push_back( i ); maskNonEmpty = true; } } if ( maskNonEmpty ) { std::list< int >::iterator it; for ( it = relevantLabels.begin(), i = 0; it != relevantLabels.end(); ++it, ++i ) { Statistics statistics; // restore previous code 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.SetVariance(labelStatisticsFilter->GetVariance( *it )); statistics.SetSigma(labelStatisticsFilter->GetSigma( *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->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 maxIndex; vnl_vector 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; ipush_back( statistics ); } } else { histogramContainer->push_back( HistogramType::ConstPointer( m_EmptyHistogram ) ); statisticsContainer->push_back( Statistics() ); } } template ImageStatisticsCalculator::ImageExtrema ImageStatisticsCalculator::CalculateExtremaWorld( const itk::Image *inputImage, itk::Image *maskImage, double neccessaryDistanceToImageBorderInMM, unsigned int label) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef itk::ImageRegionConstIteratorWithIndex MaskImageIteratorType; typedef itk::ImageRegionConstIteratorWithIndex 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::min(); float minValue = itk::NumericTraits::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 != NULL) { 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 itk::Size 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( 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 itk::SmartPointer< itk::Image > 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(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 MaskIteratorType; MaskIteratorType maskIt(convolutionKernel,maskRegion); int numberOfSubVoxelsPerDimension = 2; // per dimension! int numberOfSubVoxels = ::pow( static_cast(numberOfSubVoxelsPerDimension), static_cast(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 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 itk::SmartPointer > ImageStatisticsCalculator::GenerateConvolutionImage( const itk::Image* 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(mmPerPixel, m_HotspotRadiusInMM); // update convolution image typedef itk::Image< TPixel, VImageDimension > InputImageType; typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; typedef itk::FFTConvolutionImageFilter ConvolutionFilterType; typename ConvolutionFilterType::Pointer convolutionFilter = ConvolutionFilterType::New(); typedef itk::ConstantBoundaryCondition 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* maskImage, itk::Point sphereCenter, double sphereRadiusInMM) { typedef itk::Image< TPixel, VImageDimension > MaskImageType; typedef itk::ImageRegionIteratorWithIndex 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* inputImage, itk::Image* 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(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::LabelStatisticsImageFilter< InputImageType, MaskImageType> LabelStatisticsFilterType; typename LabelStatisticsFilterType::Pointer labelStatisticsFilter; labelStatisticsFilter = LabelStatisticsFilterType::New(); labelStatisticsFilter->SetInput( inputImage ); labelStatisticsFilter->SetLabelInput( hotspotMaskITK ); 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 points = vtkSmartPointer::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 holePoints = NULL; if (!planarFigureHolePolyline.empty()) { holePoints = vtkSmartPointer::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 lassoStencil = vtkSmartPointer::New(); lassoStencil->SetShapeToPolygon(); lassoStencil->SetPoints( points ); vtkSmartPointer holeLassoStencil = NULL; if (holePoints.GetPointer() != NULL) { holeLassoStencil = vtkSmartPointer::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 vtkImporter = vtkSmartPointer::New(); this->ConnectPipelines( itkExporter, vtkImporter ); // Apply the generated image stencil to the input image vtkSmartPointer imageStencilFilter = vtkSmartPointer::New(); imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() ); imageStencilFilter->SetStencilConnection(lassoStencil->GetOutputPort()); imageStencilFilter->ReverseStencilOff(); imageStencilFilter->SetBackgroundValue( 0 ); imageStencilFilter->Update(); vtkSmartPointer holeStencilFilter = NULL; if (holeLassoStencil.GetPointer() != NULL) { holeStencilFilter = vtkSmartPointer::New(); holeStencilFilter->SetInputConnection(imageStencilFilter->GetOutputPort()); holeStencilFilter->SetStencilConnection(holeLassoStencil->GetOutputPort()); holeStencilFilter->ReverseStencilOn(); holeStencilFilter->SetBackgroundValue(0); holeStencilFilter->Update(); } // Export from VTK back to ITK vtkSmartPointer vtkExporter = vtkSmartPointer::New(); vtkExporter->SetInputConnection( holeStencilFilter.GetPointer() == NULL ? 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() ); } } diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.h b/Modules/ImageStatistics/mitkImageStatisticsCalculator.h index b72422193b..5d1d526574 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.h +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.h @@ -1,573 +1,579 @@ /*=================================================================== 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. ===================================================================*/ #ifndef _MITK_IMAGESTATISTICSCALCULATOR_H #define _MITK_IMAGESTATISTICSCALCULATOR_H #include #include "MitkImageStatisticsExports.h" #include #include #include "mitkImage.h" #include "mitkPlanarFigure.h" #ifndef __itkHistogram_h #include #endif #include #include // just a helper to unclutter our code // to be replaced with references to m_Member (when deprecated public members in Statistics are removed) #define mitkSetGetConstMacro(name, type) \ virtual type Get##name() const \ { \ return this->name; \ } \ \ virtual void Set##name(const type _arg) \ { \ if ( this->name != _arg ) \ { \ this->name = _arg; \ } \ } namespace mitk { /** * \brief Class for calculating statistics and histogram for an (optionally * masked) image. * * Images can be masked by either a label image (of the same dimensions as * the original image) or by a closed mitk::PlanarFigure, e.g. a circle or * polygon. When masking with a planar figure, the slice corresponding to the * plane containing the figure is extracted and then clipped with contour * defined by the figure. Planar figures need to be aligned along the main axes * of the image (axial, sagittal, coronal). Planar figures on arbitrary * rotated planes are not supported. * * For each operating mode (no masking, masking by image, masking by planar * figure), the calculated statistics and histogram are cached so that, when * switching back and forth between operation modes without modifying mask or * image, the information doesn't need to be recalculated. * * The class also has the possibility to calculate the location and separate * statistics for a region called "hotspot". The hotspot is a sphere of * user-defined size and its location is chosen in a way that the average * pixel value within the sphere is maximized. * * \warning Hotspot calculation does not work in case of 2D-images! * * Note: currently time-resolved and multi-channel pictures are not properly * supported. * * \section HotspotStatistics_caption Calculation of hotspot statistics * * Since calculation of hotspot location and statistics is not * straight-forward, the following paragraphs will describe it in more detail. * * Note: Calculation of hotspot statistics is optional and set to off by default. * Multilabel-masks are supported. * * \subsection HotspotStatistics_description Hotspot Definition * * The hotspot of an image is motivated from PET readings. It is defined * as a spherical region of fixed size which maximizes the average pixel value * within the region. The following image illustrates the concept: the * colored areas are different image intensities and the hotspot is located * in the hottest region of the image. * * Note: Only hotspots are calculated for which the whole hotspot-sphere is * inside the image by default. This behaviour can be changed by * by calling SetHotspotMustBeCompletlyInsideImage(). * \warning Note that SetHotspotMustBeCompletlyInsideImage(false) may overrate * "hot" regions at image borders, because they have a stronger influence on the * mean value! Think clearly about this fact and make sure this is what you * want/need in your application, before calling * SetHotspotMustBeCompletlyInsideImage(false)! * * * \image html hotspotexample.JPG * * \subsection HotspotStatistics_calculation Hotspot Calculation * * Since only the size of the hotspot is known initially, we need to calculate * two aspects (both implemented in CalculateHotspotStatistics() ): * - the hotspot location * - statistics of the pixels within the hotspot. * * Finding the hotspot location requires to calculate the average value at each * position. This is done by convolution of the image with a sperical kernel * image which reflects partial volumes (important in the case of low-resolution * PET images). * * Once the hotspot location is known, calculating the actual statistics is a * simple task which is implemented in CalculateHotspotStatistics() using a second * instance of the ImageStatisticsCalculator. * * Step 1: Finding the hotspot by image convolution * * As described above, we use image convolution with a rasterized sphere to * average the image at each position. To handle coarse resolutions, which would * normally force us to decide for partially contained voxels whether to count * them or not, we supersample the kernel image and use non-integer kernel values * (see GenerateHotspotSearchConvolutionKernel()), which reflect the volume part that is contained in the * sphere. For example, if three subvoxels are inside the sphere, the corresponding * kernel voxel gets a value of 0.75 (3 out of 4 subvoxels, see 2D example below). * * \image html convolutionkernelsupersampling.jpg * * Convolution itself is done by means of the itkFFTConvolutionImageFilter. * To find the hotspot location, we simply iterate the averaged image and find a * maximum location (see CalculateExtremaWorld()). In case of images with multiple * maxima the method returns value and corresponding index of the extrema that is * found by the iterator first. * * Step 2: Computation of hotspot statistics * * Once the hotspot location is found, statistics for the region are calculated * by simply iterating the input image and regarding all pixel centers inside the * hotspot-sphere for statistics. * \warning Index positions of maximum/minimum are not provided, because they are not necessarily unique * \todo If index positions of maximum/minimum are required, output needs to be changed to multiple positions / regions, etc. * * \subsection HotspotStatistics_tests Tests * * To check the correctness of the hotspot calculation, a special class * (\ref hotspottestdoc) has been created, which generates images with * known hotspot location and statistics. A number of unit tests use this class * to first generate an image of known properites and then verify that * ImageStatisticsCalculator is able to reproduce the known statistics. * */ class MITKIMAGESTATISTICS_EXPORT ImageStatisticsCalculator : public itk::Object { public: /** \brief Enum for possible masking modi. */ enum { MASKING_MODE_NONE = 0, MASKING_MODE_IMAGE = 1, MASKING_MODE_PLANARFIGURE = 2 }; typedef itk::Statistics::Histogram HistogramType; typedef HistogramType::ConstIterator HistogramConstIteratorType; /** \brief Class for common statistics, includig hotspot properties. */ class MITKIMAGESTATISTICS_EXPORT Statistics { public: Statistics(bool withHotspotStatistics = true); Statistics(const Statistics& other); virtual ~Statistics(); Statistics& operator=(Statistics const& stats); const Statistics& GetHotspotStatistics() const; // real statistics Statistics& GetHotspotStatistics(); // real statistics bool HasHotspotStatistics() const; void SetHasHotspotStatistics(bool hasHotspotStatistics); // set a flag. if set, return empty hotspotstatistics object void Reset(unsigned int dimension = 2); mitkSetGetConstMacro(Label, unsigned int) mitkSetGetConstMacro(N, unsigned int) mitkSetGetConstMacro(Min, double) mitkSetGetConstMacro(Max, double) mitkSetGetConstMacro(Mean, double) mitkSetGetConstMacro(Median, double) double GetVariance() const; /** \brief Set variance * * This method checks whether the variance is negative: * The reason that the variance may be negative is that the underlying itk::LabelStatisticsImageFilter uses a naïve algorithm * for calculating the variance ( http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance ) which can lead to negative values * due to rounding errors. * * If the variance is negative the value will be set to 0.0, else the given value will be set. */ void SetVariance( const double ); double GetSigma() const; /** \brief Set standard deviation (sigma) * * This method checks if the given standard deviation is a positive value. This is done because the underlying itk::LabelStatisticsImageFilter uses * a naïve algorithm to calculate the variance. This may lead to a negative variance and because the square root of the variance is taken it also * leads to NaN for sigma. * * If the given value is not reasonable the value will be set to 0.0, else the given value will be set. * * \see SetVariance() */ void SetSigma( const double ); mitkSetGetConstMacro(RMS, double) mitkSetGetConstMacro(MinIndex, vnl_vector) mitkSetGetConstMacro(MaxIndex, vnl_vector) mitkSetGetConstMacro(HotspotIndex, vnl_vector) private: unsigned int Label; unsigned int N; double Min; double Max; double Mean; double Median; double Variance; double Sigma; double RMS; vnl_vector MinIndex; vnl_vector MaxIndex; Statistics* m_HotspotStatistics; bool m_HasHotspotStatistics; vnl_vector HotspotIndex; //< index of hotspotsphere origin }; typedef std::vector< HistogramType::ConstPointer > HistogramContainer; typedef std::vector< Statistics > StatisticsContainer; mitkClassMacro( ImageStatisticsCalculator, itk::Object ); itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** \brief Automatically calculate bin size to obtain 200 bins. */ void SetUseDefaultBinSize(bool useDefault); /** \brief Set image from which to compute statistics. */ void SetImage( const mitk::Image *image ); /** \brief Set image for masking. */ void SetImageMask( const mitk::Image *imageMask ); /** \brief Set planar figure for masking. */ void SetPlanarFigure( mitk::PlanarFigure *planarFigure ); /** \brief Set/Get operation mode for masking */ void SetMaskingMode( unsigned int mode ); /** \brief Set/Get operation mode for masking */ itkGetMacro( MaskingMode, unsigned int ); /** \brief Set/Get operation mode for masking */ void SetMaskingModeToNone(); /** \brief Set/Get operation mode for masking */ void SetMaskingModeToImage(); /** \brief Set/Get operation mode for masking */ void SetMaskingModeToPlanarFigure(); /** \brief Set a pixel value for pixels that will be ignored in the statistics */ void SetIgnorePixelValue(double value); /** \brief Get the pixel value for pixels that will be ignored in the statistics */ double GetIgnorePixelValue(); /** \brief Set whether a pixel value should be ignored in the statistics */ void SetDoIgnorePixelValue(bool doit); /** \brief Get whether a pixel value will be ignored in the statistics */ bool GetDoIgnorePixelValue(); /** \brief Set bin size for histogram resolution.*/ - void SetHistogramBinSize( unsigned int size); + void SetHistogramBinSize( double size); /** \brief Get bin size for histogram resolution.*/ unsigned int GetHistogramBinSize(); /** \brief Sets the radius for the hotspot */ void SetHotspotRadiusInMM (double hotspotRadiusInMM); /** \brief Returns the radius of the hotspot */ double GetHotspotRadiusInMM(); /** \brief Sets whether the hotspot should be calculated */ void SetCalculateHotspot(bool calculateHotspot); /** \brief Returns true whether the hotspot should be calculated, otherwise false */ bool IsHotspotCalculated(); /** \brief Sets flag whether hotspot is completly inside the image. Please note that if set to false it can be possible that statistics are calculated for which the whole hotspot is not inside the image! \warning regarding positions at the image centers may produce unexpected hotspot locations, please see \ref HotspotStatistics_description */ void SetHotspotMustBeCompletlyInsideImage(bool hotspotIsCompletlyInsideImage, bool warn = true); /** \brief Returns true if hotspot has to be completly inside the image. */ bool GetHotspotMustBeCompletlyInsideImage() const; /** \brief Compute statistics (together with histogram) for the current * masking mode. * * Computation is not executed if statistics is already up to date. In this * case, false is returned; otherwise, true.*/ virtual bool ComputeStatistics( unsigned int timeStep = 0 ); /** \brief Retrieve the histogram depending on the current masking mode. * * \param label The label for which to retrieve the histogram in multi-label situations (ascending order). */ const HistogramType *GetHistogram( unsigned int timeStep = 0, unsigned int label = 0 ) const; /** \brief Retrieve the histogram depending on the current masking mode (for all image labels. */ const HistogramContainer &GetHistogramVector( unsigned int timeStep = 0 ) const; /** \brief Retrieve statistics depending on the current masking mode. * * \param label The label for which to retrieve the statistics in multi-label situations (ascending order). */ const Statistics &GetStatistics( unsigned int timeStep = 0, unsigned int label = 0 ) const; /** \brief Retrieve statistics depending on the current masking mode (for all image labels). */ const StatisticsContainer &GetStatisticsVector( unsigned int timeStep = 0 ) const; protected: typedef std::vector< HistogramContainer > HistogramVector; typedef std::vector< StatisticsContainer > StatisticsVector; typedef std::vector< itk::TimeStamp > TimeStampVectorType; typedef std::vector< bool > BoolVectorType; typedef itk::Image< unsigned short, 3 > MaskImage3DType; typedef itk::Image< unsigned short, 2 > MaskImage2DType; ImageStatisticsCalculator(); virtual ~ImageStatisticsCalculator(); /** \brief Depending on the masking mode, the image and mask from which to * calculate statistics is extracted from the original input image and mask * data. * * For example, a when using a PlanarFigure as mask, the 2D image slice * corresponding to the PlanarFigure will be extracted from the original * image. If masking is disabled, the original image is simply passed * through. */ void ExtractImageAndMask( unsigned int timeStep = 0 ); /** \brief If the passed vector matches any of the three principal axes * of the passed geometry, the ínteger value corresponding to the axis * is set and true is returned. */ bool GetPrincipalAxis( const BaseGeometry *geometry, Vector3D vector, unsigned int &axis ); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsUnmasked( const itk::Image< TPixel, VImageDimension > *image, StatisticsContainer* statisticsContainer, HistogramContainer *histogramContainer ); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsMasked( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage, StatisticsContainer* statisticsContainer, HistogramContainer* histogramContainer ); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateMaskFromPlanarFigure( const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ); template < typename TPixel, unsigned int VImageDimension > void InternalMaskIgnoredPixels( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage ); class ImageExtrema { public: bool Defined; double Max; double Min; vnl_vector MaxIndex; vnl_vector MinIndex; ImageExtrema() :Defined(false) ,Max(itk::NumericTraits::min()) ,Min(itk::NumericTraits::max()) { } }; /** \brief Calculates minimum, maximum, mean value and their * corresponding indices in a given ROI. As input the function * needs an image and a mask. Returns an ImageExtrema object. */ template ImageExtrema CalculateExtremaWorld( const itk::Image *inputImage, itk::Image *maskImage, double neccessaryDistanceToImageBorderInMM, unsigned int label); /** \brief Calculates the hotspot statistics depending on * masking mode. Hotspot statistics are calculated for a * hotspot which is completly located inside the image by default. */ template < typename TPixel, unsigned int VImageDimension> Statistics CalculateHotspotStatistics( const itk::Image *inputImage, itk::Image *maskImage, double radiusInMM, bool& isHotspotDefined, unsigned int label); /** Connection from ITK to VTK */ template void ConnectPipelines(ITK_Exporter exporter, vtkSmartPointer 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 void ConnectPipelines(vtkSmartPointer 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()); } void UnmaskedStatisticsProgressUpdate(); void MaskedStatisticsProgressUpdate(); /** \brief Returns size of convolution kernel depending on spacing and radius. */ template itk::Size CalculateConvolutionKernelSize(double spacing[VImageDimension], double radiusInMM); /** \brief Generates image of kernel which is needed for convolution. */ template itk::SmartPointer< itk::Image > GenerateHotspotSearchConvolutionKernel(double spacing[VImageDimension], double radiusInMM); /** \brief Convolves image with spherical kernel image. Used for hotspot calculation. */ template itk::SmartPointer< itk::Image > GenerateConvolutionImage( const itk::Image* inputImage ); /** \brief Fills pixels of the spherical hotspot mask. */ template < typename TPixel, unsigned int VImageDimension> void FillHotspotMaskPixels( itk::Image* maskImage, itk::Point sphereCenter, double sphereRadiusInMM); /** m_Image contains the input image (e.g. 2D, 3D, 3D+t)*/ mitk::Image::ConstPointer m_Image; mitk::Image::ConstPointer m_ImageMask; mitk::PlanarFigure::Pointer m_PlanarFigure; HistogramVector m_ImageHistogramVector; HistogramVector m_MaskedImageHistogramVector; HistogramVector m_PlanarFigureHistogramVector; HistogramType::Pointer m_EmptyHistogram; HistogramContainer m_EmptyHistogramContainer; StatisticsVector m_ImageStatisticsVector; StatisticsVector m_MaskedImageStatisticsVector; StatisticsVector m_PlanarFigureStatisticsVector; StatisticsVector m_MaskedImageHotspotStatisticsVector; Statistics m_EmptyStatistics; StatisticsContainer m_EmptyStatisticsContainer; unsigned int m_MaskingMode; bool m_MaskingModeChanged; /** m_InternalImage contains a image volume at one time step (e.g. 2D, 3D)*/ mitk::Image::ConstPointer m_InternalImage; MaskImage3DType::Pointer m_InternalImageMask3D; MaskImage2DType::Pointer m_InternalImageMask2D; TimeStampVectorType m_ImageStatisticsTimeStampVector; TimeStampVectorType m_MaskedImageStatisticsTimeStampVector; TimeStampVectorType m_PlanarFigureStatisticsTimeStampVector; BoolVectorType m_ImageStatisticsCalculationTriggerVector; BoolVectorType m_MaskedImageStatisticsCalculationTriggerVector; BoolVectorType m_PlanarFigureStatisticsCalculationTriggerVector; double m_IgnorePixelValue; bool m_DoIgnorePixelValue; bool m_IgnorePixelValueChanged; unsigned int m_PlanarFigureAxis; // Normal axis for PlanarFigure unsigned int m_PlanarFigureSlice; // Slice which contains PlanarFigure int m_PlanarFigureCoordinate0; // First plane-axis for PlanarFigure int m_PlanarFigureCoordinate1; // Second plane-axis for PlanarFigure - unsigned int m_HistogramBinSize; ///Bin size for histogram resoluion. + double m_HistogramBinSize; ///Bin size for histogram resoluion. bool m_UseDefaultBinSize; double m_HotspotRadiusInMM; bool m_CalculateHotspot; bool m_HotspotRadiusInMMChanged; bool m_HotspotMustBeCompletelyInsideImage; + +private: + + unsigned int calcNumberOfBins(mitk::ScalarType min, mitk::ScalarType max); + + }; } // namespace #endif diff --git a/Modules/QtWidgetsExt/include/QmitkHistogramJSWidget.h b/Modules/QtWidgetsExt/include/QmitkHistogramJSWidget.h index d881ac4f8c..7cf55969e4 100644 --- a/Modules/QtWidgetsExt/include/QmitkHistogramJSWidget.h +++ b/Modules/QtWidgetsExt/include/QmitkHistogramJSWidget.h @@ -1,278 +1,278 @@ /*=================================================================== 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. ===================================================================*/ #ifndef QMITKHISTOGRAMJSWIDGET_H #define QMITKHISTOGRAMJSWIDGET_H #include #include #include #include "MitkQtWidgetsExtExports.h" #include #include "mitkImage.h" #include "mitkPlanarFigure.h" #include #include /** * \brief Widget which shows a histogram using JavaScript. * * This class is a QWebView. It shows the histogram for a selected image * or segmentation. It also can display an intensity profile for * path elements, which lais over an image. */ class MITKQTWIDGETSEXT_EXPORT QmitkHistogramJSWidget : public QWebView { Q_OBJECT /** * \brief Measurement property. * * This property is used in JavaScript as member of the current object. * It holds a QList, containing the measurements of the current histogram. * @see GetMeasurement() */ Q_PROPERTY(QList measurement READ GetMeasurement) /** * \brief Frequency property. * * This property is used in JavaScript as member of the current object. * It holds a QList, containing the frequencies of the current histogram. * @see GetFrequency() */ Q_PROPERTY(QList frequency READ GetFrequency) /** * \brief Line graph property. * * This property is used in JavaScript as member of the current object. * It holds a boolean, which sais wether to use a line or not. * @see GetUseLineGraph() */ Q_PROPERTY(bool useLineGraph READ GetUseLineGraph) /** * @brief intensity profile property. * * This property is used in JavaScript as member of the current object. - * It holds a boolean, which sais wether to use an intensity profile or not. + * It holds a boolean, which says whether to use an intensity profile or not. * @see GetIntensityProfile() */ Q_PROPERTY(bool intensityProfile READ GetIntensityProfile) public: typedef mitk::Image::HistogramType HistogramType; typedef mitk::Image::HistogramType::ConstIterator HistogramConstIteratorType; typedef itk::PolyLineParametricPath< 3 > ParametricPathType; typedef itk::ParametricPath< 3 >::Superclass PathType; typedef mitk::PlanarFigure::PolyLineType VertexContainerType; explicit QmitkHistogramJSWidget(QWidget *parent = 0); ~QmitkHistogramJSWidget(); /** * \brief Event which notifies a change of the widget size. * * Reimplemented from QWebView::resizeEvent(), * reloads the webframe */ void resizeEvent(QResizeEvent* resizeEvent); /** * \brief Calculates the histogram. * * This function removes all frequencies of 0 until the first bin and behind the last bin. * It writes the measurement and frequency, which are given from the HistogramType, into * m_Measurement and m_Frequency. * The SignalDataChanged is called, to update the information, which is displayed in the webframe. */ void ComputeHistogram(HistogramType* histogram); /** * \brief Calculates the intensityprofile. * * If an image and a pathelement are set, this function * calculates an intensity profile for a pathelement which lies over an image. * Sets m_IntensityProfile and m_UseLineGraph to true. * The SignalDataChanged is called, to update the information, which is displayed in the webframe. */ void ComputeIntensityProfile(unsigned int timeStep = 0); /** * \brief Clears the Histogram. * * This function clears the data and calls SignalDataChanged to update * the displayed information in the webframe. */ void ClearHistogram(); /** * \brief Getter for measurement. * * @return List of measurements. */ QList GetMeasurement(); /** * \brief Getter for frequency. * * @return List of frequencies. */ QList GetFrequency(); /** * \brief Getter for uselineGraph. * * @return True if a linegraph should be used. */ bool GetUseLineGraph(); /** * \brief Getter for intensity profile. * * @return True if current histogram is an intensityprofile */ bool GetIntensityProfile(); /** * \brief Setter for reference image. * * @param image The corresponding image for an intensity profile. */ void SetImage(mitk::Image* image); /** * \brief Setter for planarFigure. * * @param planarFigure The pathelement for an intensity profile. */ void SetPlanarFigure(const mitk::PlanarFigure* planarFigure); private: /** * \brief List of frequencies. * * A QList which holds the frequencies of the current histogram * or holds the intesities of current intensity profile. */ QList m_Frequency; /** * \brief List of measurements. * * A QList which holds the measurements of the current histogram * or holds the distances of current intensity profile. */ QList m_Measurement; /** * \brief Reference image. * * Holds the image to calculate an intensity profile. */ mitk::Image::Pointer m_Image; /** * \brief Pathelement. * * Holds a not closed planar figure to calculate an intensity profile. */ mitk::PlanarFigure::ConstPointer m_PlanarFigure; bool m_UseLineGraph; bool m_IntensityProfile; /** * Holds the current histogram */ HistogramType::ConstPointer m_Histogram; /** * Path derived either form user-specified path or from PlanarFigure-generated * path */ PathType::ConstPointer m_DerivedPath; /** * Parametric path as generated from PlanarFigure */ ParametricPathType::Pointer m_ParametricPath; /** * \brief Clears data. * * Clears the QLists m_Measurement and m_Frequency */ void ClearData(); QmitkJSWebPage* m_Page; private slots: /** * \brief Adds an object to JavaScript. * * Adds an object of the widget to JavaScript. * By using this object JavaScript can react to the signals of the widget * and can access the QProperties as members of the object. */ void AddJSObject(); public slots: /** * \brief Slot for radiobutton m_barRadioButton. * * Sets m_UseLineGraph to false. * Calls signal GraphChanged to update the graph in the webframe. */ void OnBarRadioButtonSelected(); /** * \brief Slot for radiobutton m_lineRadioButton. * * Sets m_UseLineGraph to true. * Calls signal GraphChanged to update the graph in the webframe. */ void OnLineRadioButtonSelected(); signals: /** * \brief Signal data has changed. * * It has to be called when the data of the histogram or intensity profile has changed. */ void SignalDataChanged(); /** * \brief Signal graph has changed. * * It has to be called when the graph changed from barchart to linegraph. Vice versa. */ void SignalGraphChanged(); }; #endif diff --git a/Modules/QtWidgetsExt/resource/Histogram.js b/Modules/QtWidgetsExt/resource/Histogram.js index 07187f7f3c..07d1751f54 100644 --- a/Modules/QtWidgetsExt/resource/Histogram.js +++ b/Modules/QtWidgetsExt/resource/Histogram.js @@ -1,546 +1,554 @@ /*=================================================================== 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. ===================================================================*/ var margin = { top : 10, bottom : 50, left : 45, right : 20, }; var height = histogramData.height - margin.top - margin.bottom; var width = histogramData.width - margin.left - margin.right; var tension = 0.8; var connected = false; var dur = 1000; var binSize = 10; var min; var max; /* * Connecting signals from qt side with JavaScript methods. */ if (!connected) { connected = true; histogramData.SignalDataChanged.connect(updateHistogram); histogramData.SignalGraphChanged.connect(updateHistogram); } function disconnectSignals() { histogramData.SignalDataChanged.disconnect(updateHistogram); histogramData.SignalGraphChanged.disconnect(updateHistogram); delete histogramData; } /* * Predefinition of scales. */ var xScale = d3.scale.linear() .domain([d3.min(histogramData.measurement)-binSize/2,d3.max(histogramData.measurement)+binSize/2]) .range([0,width]); var yScale = d3.scale.linear() .domain([d3.min(histogramData.frequency),d3.max(histogramData.frequency)]) .range([height,margin.top]); /* * Predefinition of axis elements. */ var xAxis = d3.svg.axis() .scale(xScale) .orient("bottom") .tickFormat(d3.format("s")); var yAxis = d3.svg.axis() .scale(yScale) .orient("left") .tickFormat(d3.format("s")); /* * Predefinition of the zoom. */ var zoombie = d3.behavior.zoom().x(xScale).scaleExtent([1, 50]).on("zoom", zoom); /* * Creation of the svg element, which holds the complete histogram. */ var svg = d3.select("body") .append("svg") .attr("class", "svg") .attr("width", width + margin.right + margin.left) .attr("height", height + margin.top + margin.bottom) .append("g") .attr("transform", "translate (" + margin.left + "," + margin.top + ")") .call(zoombie) .on("mousemove", myMouseMove); /* * Appending a rectangle to the svg, to guarantee the possibility * of zooming on the whole histogram. */ svg.append("rect") .attr("width", width + margin.left + margin.right) .attr("height", height + margin.top + margin.bottom) .attr("opacity", 0); /* * Appending a second svg to main svg, which holds only the graph. */ var vis = svg.append("svg") .attr("width", width) .attr("height", height); /* * Predefinition of the lines. */ var line = d3.svg.line() .interpolate("linear") .x(function(d,i) { return xScale(histogramData.measurement[i]-binSize/2); }) .y(function(d) { return yScale(d); }); var linenull = d3.svg.line() .interpolate("linear") .x(function(d,i) { return xScale(histogramData.measurement[i]-binSize/2); }) .y(function(d) { return yScale(0); }); updateHistogram(); /* * Method to update the histogram data * and to change the displayed graph. */ function updateHistogram() { calcBinSize(); if (!histogramData.useLineGraph) { barChart(); } else if (histogramData.useLineGraph) { linePlot() } } /* * Calculation of the bin size. */ function calcBinSize() { - min = d3.min(histogramData.measurement); - max = d3.max(histogramData.measurement); - binSize = ((max - min) / (histogramData.measurement.length)); + if (1 < histogramData.measurement.length) + { + min = d3.min(histogramData.measurement); + max = d3.max(histogramData.measurement); + binSize = ((max - min) / (histogramData.measurement.length-1)); + } + else + { + binSize = 10; + } } /* * Method to display histogram as a barchart. */ function barChart() { definition(); /* * Change zoom to a fixed y-axis. */ zoombie = d3.behavior.zoom().x(xScale).scaleExtent([1, 50]).on("zoom", zoom); svg.call(zoombie); /* * Element to animate transition from linegraph to barchart. */ vis.selectAll("path.line").remove(); vis.selectAll("circle").remove(); /* * Definition of the bar elements. */ var bar = vis.selectAll("rect.bar").data(histogramData.frequency); /* * Definition how to handle new bar elements. */ bar.enter().append("rect") .attr("class", "bar") .on("mouseover", myMouseOver) .on("mouseout", myMouseOut) .attr("x", function(d,i) { return xScale(histogramData.measurement[i]-binSize/2); }) .attr("y", height) .attr("height", 0) .attr("width", barWidth) /* * Definition how to handle changed bar elements. */ bar.transition() .duration(dur) .attr("x", function(d,i) { return xScale(histogramData.measurement[i]-binSize/2); }) .attr("y", myYPostion) .attr("height", barHeight) .attr("width", barWidth); /* * Definition how to handle bar elements which doesn't exist anymore.' */ bar.exit() .transition() .duration(dur) .attr("y", height) .attr("height", 0) .remove(); /* * Update of axis elements. * First delete old ones, then generate new. */ svg.selectAll("g") .remove(); svg.append("g") .attr("class", "x axis") .attr("transform", "translate(0," + height + ")") .call(xAxis); svg.append("g") .attr("class", "y axis") .call(yAxis); } /* * Method to display histogram as a linegraph. */ function linePlot() { definition(); /* * Change zoom to a zoomable y-axis. */ zoombie = d3.behavior.zoom().x(xScale).y(yScale).scaleExtent([1, 50]).on("zoom", zoom); svg.call(zoombie); /* * Elements to animate transitions from barchart to linegraph. * Different transition when an intensity profile is generated. */ if(!histogramData.intensityProfile) { vis.selectAll("rect.bar") .transition() .duration(dur) .attr("height", 0) .remove(); } else { vis.selectAll("rect.bar") .transition() .duration(dur) .attr("y", height) // <-- .attr("height", 0) .remove(); } /* * Creating circle elements, when an intensity profile is generated to show tooltips. * Due performance losses tooltips are not supported for line histograms. */ if(histogramData.intensityProfile) { var circles = vis.selectAll("circle").data(histogramData.frequency); /* * Definition how to handle new circle elements. */ circles.enter() .append("circle") .on("mouseover", myMouseOverLine) .on("mouseout", myMouseOutLine) .attr("cx", function(d,i) { return xScale(histogramData.measurement[i]-binSize/2); }) .attr("cy", function (d) { return yScale(d) }) .attr("r", 5) .attr("opacity", 0) .style("stroke", "red") .style("stroke-width", 1) .style("fill-opacity", 0); /* * Definition how to handle bar elements which doesn't exist anymore. */ circles.exit().remove(); } else { /* * Removing of all circle elements if a line histogram is generated. */ vis.selectAll("circle").remove(); } /* * Creating a new path element. */ var graph = vis.selectAll("path.line") .data([histogramData.frequency]); /* * Definition how to handle a new path element, using predefined lines. */ graph.enter() .append("path") .attr("class", "line") .transition() .duration(dur) .attr("d", line); /* * Definition how to handle change points in an existing path element. */ graph.transition() .duration(dur) .attr("d", line); /* * Update of axis elements. * First delete old ones, then generate new. */ svg.selectAll("g") .remove(); svg.append("g") .attr("class", "x axis") .attr("transform", "translate(0," + height + ")") .call(xAxis); svg.append("g") .attr("class", "y axis") .call(yAxis); } function definition() { /* * Match scale to current data. */ xScale = d3.scale.linear() .domain([d3.min(histogramData.measurement)-binSize/2,d3.max(histogramData.measurement)+binSize/2]) .range([0,width]); yScale = d3.scale.linear() .domain([d3.min(histogramData.frequency),d3.max(histogramData.frequency)]) .range([height,margin.top]); /* * Match axes to current scale */ xAxis = d3.svg.axis() .scale(xScale) .orient("bottom") .tickFormat(d3.format("s")); + yAxis = d3.svg.axis() .scale(yScale) .orient("left") .tickFormat(d3.format("s")); } /* * Method to calculate barwidth in px. */ function barWidth(d, i) { var bw; if (i != (histogramData.measurement.length-1)) { bw =(xScale(histogramData.measurement[i + 1]) - xScale(histogramData.measurement[i])) * (histogramData.frequency.length / (histogramData.frequency.length + 1)) - 1; } else { bw =(xScale(histogramData.measurement[i]) - xScale(histogramData.measurement[i - 1])) * (histogramData.frequency.length / (histogramData.frequency.length + 1)) - 1; } /* * Ensure barwidth is not smaller than 1px. */ bw = bw > 1 ? bw : 1; return bw; } /* * Method to calculate barheight in px. * Ensure barheight is not smaller than 1px. */ function barHeight(d) { var bh; bh = height - yScale(d); bh = bh >=2 ? bh : 2; return bh; } /* * Method to calculate dynamical y positions. */ function myYPostion(d) { var myy = yScale(d); myy = (height-myy) > 2 ? myy : (height-2); if (d == 0) { return height; } return myy; } /* * Method to fit parameters of bars/line when zoomed. * Update axes elements. * Resets the view if scale is 1. */ function zoom() { if (zoombie.scale() == 1) { zoombie.translate([0,0]); xScale.domain([d3.min(histogramData.measurement)-binSize/2,d3.max(histogramData.measurement)+binSize/2]); yScale.domain([d3.min(histogramData.frequency),d3.max(histogramData.frequency)]); } if (!histogramData.useLineGraph) { svg.select(".x.axis").call(xAxis); vis.selectAll(".bar") .attr("width", barWidth) .attr("x", function(d, i) { return xScale(histogramData.measurement[i]-binSize/2); }); } else { svg.select(".x.axis").call(xAxis); svg.select(".y.axis").call(yAxis); vis.selectAll("path.line") .attr("transform", "translate(" + zoombie.translate() + ")scale(" + zoombie.scale() + ")") .style("stroke-width", 1 / zoombie.scale()); vis.selectAll("circle") .attr("cx", function(d, i) { return xScale(histogramData.measurement[i]-binSize/2); }) .attr("cy", function(d) { return yScale(d); }); } } /* * Method to show infobox, while mouse is over a bin. */ function myMouseOver() { var myBar = d3.select(this); var reScale = d3.scale.linear() .domain(xScale.range()) .range(xScale.domain()); var y = myBar.data(); var x = reScale(myBar.attr("x")); myBar.style("fill", "red"); d3.select(".infobox").style("display", "block"); if ((min >= 0) && (max <= 2)) //tooltip for float images { d3.select(".measurement").text("Greayvalue: " + (Math.round(x*1000)/1000)); } else { d3.select(".measurement").text("Greyvalue: " + (Math.round(x*10)/10) + " ... " + (Math.round((x+binSize)*10)/10)); } d3.select(".frequency").text("Frequency: " + y); } /* * Hide infobox, when mouse not over a bin. */ function myMouseOut() { var myBar = d3.select(this); myBar.style("fill", d3.rgb(0,71,185)); d3.select(".infobox").style("display", "none"); } /* * Show tooltip, while mouse is over a circle in an intensity profile. */ function myMouseOverLine() { var myCircle = d3.select(this) var reScale = d3.scale.linear() .domain(xScale.range()) .range(xScale.domain()); var y = myCircle.data(); var x = reScale(myCircle.attr("cx")); x = x >= 0 ? x : 0; myCircle.attr("opacity", 1); d3.select(".infobox").style("display", "block"); d3.select(".measurement").text("Distance: " + (Math.round(x*100)/100) + " mm"); d3.select(".frequency").text("Intensity: " + y); } /* * Hide infobox, when mouse not over a circle. */ function myMouseOutLine() { var myCircle = d3.select(this); myCircle.attr("opacity", 0); d3.select(".infobox").style("display", "none"); } /* * Update mousecoordinates when mouse is moved. * Tooltip is shown on the right side of the mouse until it reaches the right boundary, * the it switches to the left side. */ function myMouseMove() { var infobox = d3.select(".infobox"); var coords = d3.mouse(this); if ((coords[0]+120)<(width-margin.right)) { infobox.style("left", coords[0] + 75 + "px"); infobox.style("top", coords[1] + "px"); } else { infobox.style("left", coords[0] - 90 + "px"); infobox.style("top",coords[1] + "px"); } } diff --git a/Modules/QtWidgetsExt/src/QmitkHistogramJSWidget.cpp b/Modules/QtWidgetsExt/src/QmitkHistogramJSWidget.cpp index 76f0ccd1a9..a03468b6f3 100644 --- a/Modules/QtWidgetsExt/src/QmitkHistogramJSWidget.cpp +++ b/Modules/QtWidgetsExt/src/QmitkHistogramJSWidget.cpp @@ -1,240 +1,242 @@ /*=================================================================== 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 "QmitkHistogramJSWidget.h" #include "mitkPixelTypeMultiplex.h" #include #include #include "mitkRenderingManager.h" #include "mitkBaseRenderer.h" #include "mitkImageTimeSelector.h" #include "mitkExtractSliceFilter.h" #include QmitkHistogramJSWidget::QmitkHistogramJSWidget(QWidget *parent) : QWebView(parent) { // set histogram type to barchart in first instance m_UseLineGraph = false; m_Page = new QmitkJSWebPage(this); setPage(m_Page); // set html from source connect(page()->mainFrame(), SIGNAL(javaScriptWindowObjectCleared()), this, SLOT(AddJSObject())); QUrl myUrl = QUrl("qrc:///QtWidgetsExt/Histogram.html"); setUrl(myUrl); // set Scrollbars to be always disabled page()->mainFrame()->setScrollBarPolicy(Qt::Horizontal, Qt::ScrollBarAlwaysOff); page()->mainFrame()->setScrollBarPolicy(Qt::Vertical, Qt::ScrollBarAlwaysOff); m_ParametricPath = ParametricPathType::New(); } QmitkHistogramJSWidget::~QmitkHistogramJSWidget() { } // adds an Object of Type QmitkHistogramJSWidget to the JavaScript, using QtWebkitBridge void QmitkHistogramJSWidget::AddJSObject() { page()->mainFrame()->addToJavaScriptWindowObject(QString("histogramData"), this); } // reloads WebView, everytime its size has been changed, so the size of the Histogram fits to the size of the widget void QmitkHistogramJSWidget::resizeEvent(QResizeEvent* resizeEvent) { QWebView::resizeEvent(resizeEvent); // workaround for Qt Bug: https://bugs.webkit.org/show_bug.cgi?id=75984 page()->mainFrame()->evaluateJavaScript("disconnectSignals()"); this->reload(); } // method to expose data to JavaScript by using properties void QmitkHistogramJSWidget::ComputeHistogram(HistogramType* histogram) { m_Histogram = histogram; HistogramConstIteratorType startIt = m_Histogram->End(); HistogramConstIteratorType endIt = m_Histogram->End(); HistogramConstIteratorType it = m_Histogram->Begin(); ClearData(); unsigned int i = 0; bool firstValue = false; // removes frequencies of 0, which are outside the first and last bin for (; it != m_Histogram->End(); ++it) { if (it.GetFrequency() > 0.0) { endIt = it; if (!firstValue) { firstValue = true; startIt = it; } } } ++endIt; // generating Lists of measurement and frequencies for (it = startIt ; it != endIt; ++it, ++i) { QVariant frequency = QVariant::fromValue(it.GetFrequency()); QVariant measurement = it.GetMeasurementVector()[0]; m_Frequency.insert(i, frequency); m_Measurement.insert(i, measurement); } + m_IntensityProfile = false; this->SignalDataChanged(); } void QmitkHistogramJSWidget::ClearData() { m_Frequency.clear(); m_Measurement.clear(); } void QmitkHistogramJSWidget::ClearHistogram() { this->ClearData(); this->SignalDataChanged(); } QList QmitkHistogramJSWidget::GetFrequency() { return m_Frequency; } QList QmitkHistogramJSWidget::GetMeasurement() { return m_Measurement; } bool QmitkHistogramJSWidget::GetUseLineGraph() { return m_UseLineGraph; } void QmitkHistogramJSWidget::OnBarRadioButtonSelected() { if (m_UseLineGraph) { m_UseLineGraph = false; this->SignalGraphChanged(); } } void QmitkHistogramJSWidget::OnLineRadioButtonSelected() { if (!m_UseLineGraph) { m_UseLineGraph = true; this->SignalGraphChanged(); } } void QmitkHistogramJSWidget::SetImage(mitk::Image* image) { m_Image = image; } void QmitkHistogramJSWidget::SetPlanarFigure(const mitk::PlanarFigure* planarFigure) { m_PlanarFigure = planarFigure; } template void ReadPixel(mitk::PixelType, mitk::Image::Pointer image, itk::Index<3> indexPoint, double& value) { if (image->GetDimension() == 2) { mitk::ImagePixelReadAccessor readAccess(image, image->GetSliceData(0)); itk::Index<2> idx; idx[0] = indexPoint[0]; idx[1] = indexPoint[1]; value = readAccess.GetPixelByIndex(idx); } else if (image->GetDimension() == 3) { mitk::ImagePixelReadAccessor readAccess(image, image->GetVolumeData(0)); itk::Index<3> idx; idx[0] = indexPoint[0]; idx[1] = indexPoint[1]; idx[2] = indexPoint[2]; value = readAccess.GetPixelByIndex(idx); } else { //unhandled } } void QmitkHistogramJSWidget::ComputeIntensityProfile(unsigned int timeStep) { this->ClearData(); m_ParametricPath->Initialize(); if (m_PlanarFigure.IsNull()) { mitkThrow() << "PlanarFigure not set!"; } if (m_Image.IsNull()) { mitkThrow() << "Image not set!"; } mitk::Image::Pointer image; if (m_Image->GetDimension() == 4) { mitk::ImageTimeSelector::Pointer timeSelector = mitk::ImageTimeSelector::New(); timeSelector->SetInput(m_Image); timeSelector->SetTimeNr(timeStep); timeSelector->Update(); image = timeSelector->GetOutput(); } else { image = m_Image; } mitk::IntensityProfile::Pointer intensityProfile = mitk::ComputeIntensityProfile(image, const_cast(m_PlanarFigure.GetPointer())); m_Frequency.clear(); m_Measurement.clear(); int i = -1; mitk::IntensityProfile::ConstIterator end = intensityProfile->End(); for (mitk::IntensityProfile::ConstIterator it = intensityProfile->Begin(); it != end; ++it) { m_Frequency.push_back(it.GetMeasurementVector()[0]); m_Measurement.push_back(++i); } m_IntensityProfile = true; m_UseLineGraph = true; this->SignalDataChanged(); } bool QmitkHistogramJSWidget::GetIntensityProfile() { return m_IntensityProfile; } + diff --git a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsCalculationThread.cpp b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsCalculationThread.cpp index c9cf07bade..d3454f7e98 100644 --- a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsCalculationThread.cpp +++ b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsCalculationThread.cpp @@ -1,209 +1,209 @@ /*=================================================================== 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 "QmitkImageStatisticsCalculationThread.h" //QT headers #include #include QmitkImageStatisticsCalculationThread::QmitkImageStatisticsCalculationThread():QThread(), m_StatisticsImage(NULL), m_BinaryMask(NULL), m_PlanarFigureMask(NULL), m_TimeStep(0), - m_IgnoreZeros(false), m_CalculationSuccessful(false), m_StatisticChanged(false), m_HistogramBinSize(1), m_UseDefaultBinSize(true) + m_IgnoreZeros(false), m_CalculationSuccessful(false), m_StatisticChanged(false), m_HistogramBinSize(1.0), m_UseDefaultBinSize(true) { } QmitkImageStatisticsCalculationThread::~QmitkImageStatisticsCalculationThread() { } void QmitkImageStatisticsCalculationThread::Initialize( mitk::Image::Pointer image, mitk::Image::Pointer binaryImage, mitk::PlanarFigure::Pointer planarFig ) { // reset old values if( this->m_StatisticsImage.IsNotNull() ) this->m_StatisticsImage = 0; if( this->m_BinaryMask.IsNotNull() ) this->m_BinaryMask = 0; if( this->m_PlanarFigureMask.IsNotNull()) this->m_PlanarFigureMask = 0; // set new values if passed in if(image.IsNotNull()) this->m_StatisticsImage = image->Clone(); if(binaryImage.IsNotNull()) this->m_BinaryMask = binaryImage->Clone(); if(planarFig.IsNotNull()) this->m_PlanarFigureMask = planarFig->Clone(); } void QmitkImageStatisticsCalculationThread::SetUseDefaultBinSize(bool useDefault) { m_UseDefaultBinSize = useDefault; } void QmitkImageStatisticsCalculationThread::SetTimeStep( int times ) { this->m_TimeStep = times; } int QmitkImageStatisticsCalculationThread::GetTimeStep() { return this->m_TimeStep; } std::vector QmitkImageStatisticsCalculationThread::GetStatisticsData() { return this->m_StatisticsVector; } mitk::Image::Pointer QmitkImageStatisticsCalculationThread::GetStatisticsImage() { return this->m_StatisticsImage; } void QmitkImageStatisticsCalculationThread::SetIgnoreZeroValueVoxel(bool _arg) { this->m_IgnoreZeros = _arg; } bool QmitkImageStatisticsCalculationThread::GetIgnoreZeroValueVoxel() { return this->m_IgnoreZeros; } -void QmitkImageStatisticsCalculationThread::SetHistogramBinSize(unsigned int size) +void QmitkImageStatisticsCalculationThread::SetHistogramBinSize(double size) { this->m_HistogramBinSize = size; } unsigned int QmitkImageStatisticsCalculationThread::GetHistogramBinSize() { return this->m_HistogramBinSize; } std::string QmitkImageStatisticsCalculationThread::GetLastErrorMessage() { return m_message; } QmitkImageStatisticsCalculationThread::HistogramType::Pointer QmitkImageStatisticsCalculationThread::GetTimeStepHistogram(unsigned int t) { if (t >= this->m_HistogramVector.size()) return NULL; return this->m_HistogramVector[t]; } bool QmitkImageStatisticsCalculationThread::GetStatisticsChangedFlag() { return m_StatisticChanged; } bool QmitkImageStatisticsCalculationThread::GetStatisticsUpdateSuccessFlag() { return m_CalculationSuccessful; } void QmitkImageStatisticsCalculationThread::run() { bool statisticCalculationSuccessful = true; mitk::ImageStatisticsCalculator::Pointer calculator = mitk::ImageStatisticsCalculator::New(); if(this->m_StatisticsImage.IsNotNull()) { calculator->SetImage(m_StatisticsImage); calculator->SetMaskingModeToNone(); } else { statisticCalculationSuccessful = false; } // Bug 13416 : The ImageStatistics::SetImageMask() method can throw exceptions, i.e. when the dimensionality // of the masked and input image differ, we need to catch them and mark the calculation as failed // the same holds for the ::SetPlanarFigure() try { if(this->m_BinaryMask.IsNotNull()) { calculator->SetImageMask(m_BinaryMask); calculator->SetMaskingModeToImage(); } if(this->m_PlanarFigureMask.IsNotNull()) { calculator->SetPlanarFigure(m_PlanarFigureMask); calculator->SetMaskingModeToPlanarFigure(); } } catch( const itk::ExceptionObject& e) { MITK_ERROR << "ITK Exception:" << e.what(); statisticCalculationSuccessful = false; } bool statisticChanged = false; calculator->SetDoIgnorePixelValue(this->m_IgnoreZeros); calculator->SetIgnorePixelValue(0); calculator->SetHistogramBinSize( m_HistogramBinSize ); calculator->SetUseDefaultBinSize( m_UseDefaultBinSize ); for (unsigned int i = 0; i < m_StatisticsImage->GetTimeSteps(); i++) { try { statisticChanged = calculator->ComputeStatistics(i); } catch ( mitk::Exception& e) { //m_message = e.GetDescription(); MITK_ERROR<< "MITK Exception: " << e.what(); statisticCalculationSuccessful = false; } catch ( const std::runtime_error &e ) { //m_message = "Failure: " + std::string(e.what()); MITK_ERROR<< "Runtime Exception: " << e.what(); statisticCalculationSuccessful = false; } catch ( const std::exception &e ) { //m_message = "Failure: " + std::string(e.what()); MITK_ERROR<< "Standard Exception: " << e.what(); statisticCalculationSuccessful = false; } } this->m_StatisticChanged = statisticChanged; this->m_CalculationSuccessful = statisticCalculationSuccessful; if(statisticCalculationSuccessful) { this->m_StatisticsVector.clear(); this->m_HistogramVector.clear(); for (unsigned int i = 0; i < m_StatisticsImage->GetTimeSteps(); i++) { this->m_StatisticsVector.push_back(calculator->GetStatistics(i)); this->m_HistogramVector.push_back((HistogramType*)calculator->GetHistogram(i)); } } m_HistogramBinSize = calculator->GetHistogramBinSize(); } diff --git a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsCalculationThread.h b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsCalculationThread.h index d58f06da14..e842bfa7b8 100644 --- a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsCalculationThread.h +++ b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsCalculationThread.h @@ -1,115 +1,115 @@ /*=================================================================== 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. ===================================================================*/ #ifndef QMITKIMAGESTATISTICSCALCULATIONTHREAD_H_INCLUDED #define QMITKIMAGESTATISTICSCALCULATIONTHREAD_H_INCLUDED //QT headers #include #include //mitk headers #include "mitkImage.h" #include "mitkPlanarFigure.h" #include "mitkImageStatisticsCalculator.h" // itk headers #ifndef __itkHistogram_h #include #endif /** /brief This class is executed as background thread for image statistics calculation. * Documentation: This class is derived from QThread and is intended to be used by QmitkImageStatisticsView to run the image statistics calculation in a background thread keepung the gui usable. * \ingroup Plugins/MeasurementToolbox */ class QmitkImageStatisticsCalculationThread : public QThread { Q_OBJECT public: typedef itk::Statistics::Histogram HistogramType; /*! /brief standard constructor. */ QmitkImageStatisticsCalculationThread(); /*! /brief standard destructor. */ ~QmitkImageStatisticsCalculationThread(); /*! *\brief Automatically calculate bin size to obtain 200 bins. */ void SetUseDefaultBinSize(bool useDefault); /*! /brief Initializes the object with necessary data. */ void Initialize( mitk::Image::Pointer image, mitk::Image::Pointer binaryImage, mitk::PlanarFigure::Pointer planarFig ); /*! /brief returns the calculated image statistics. */ std::vector GetStatisticsData(); /*! /brief */ mitk::Image::Pointer GetStatisticsImage(); /*! /brief Set the time step of the image you want to process. */ void SetTimeStep( int times ); /*! /brief Get the time step of the image you want to process. */ int GetTimeStep(); /*! /brief Set flag to ignore zero valued voxels */ void SetIgnoreZeroValueVoxel( bool _arg ); /*! /brief Get status of zero value voxel ignoring. */ bool GetIgnoreZeroValueVoxel(); /*! /brief Set bin size for histogram resolution.*/ - void SetHistogramBinSize( unsigned int size); + void SetHistogramBinSize( double size); /*! /brief Get bin size for histogram resolution.*/ unsigned int GetHistogramBinSize(); /*! /brief Returns the histogram of the currently selected time step. */ HistogramType::Pointer GetTimeStepHistogram(unsigned int t = 0); /*! /brief Returns a flag indicating if the statistics have changed during calculation */ bool GetStatisticsChangedFlag(); /*! /brief Returns a flag the indicates if the statistics are updated successfully */ bool GetStatisticsUpdateSuccessFlag(); /*! /brief Method called once the thread is executed. */ void run(); std::string GetLastErrorMessage(); private: //member declaration mitk::Image::Pointer m_StatisticsImage; ///< member variable holds the input image for which the statistics need to be calculated. mitk::Image::Pointer m_BinaryMask; ///< member variable holds the binary mask image for segmentation image statistics calculation. mitk::PlanarFigure::Pointer m_PlanarFigureMask; ///< member variable holds the planar figure for segmentation image statistics calculation. std::vector m_StatisticsVector; ///< member variable holds the result structs. int m_TimeStep; ///< member variable holds the time step for statistics calculation bool m_IgnoreZeros; ///< member variable holds flag to indicate if zero valued voxel should be suppressed - unsigned int m_HistogramBinSize; ///< member variable holds the bin size for histogram resolution. + double m_HistogramBinSize; ///< member variable holds the bin size for histogram resolution. bool m_StatisticChanged; ///< flag set if statistics have changed bool m_CalculationSuccessful; ///< flag set if statistics calculation was successful std::vector m_HistogramVector; ///< member holds the histograms of all time steps. std::string m_message; bool m_UseDefaultBinSize; }; #endif // QMITKIMAGESTATISTICSCALCULATIONTHREAD_H_INCLUDED diff --git a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsView.cpp b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsView.cpp index 879986591e..a4b3193ae8 100644 --- a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsView.cpp +++ b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsView.cpp @@ -1,998 +1,998 @@ /*=================================================================== 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 "QmitkImageStatisticsView.h" // Qt includes #include #include // berry includes #include // mitk includes #include "mitkNodePredicateDataType.h" #include "mitkPlanarFigureInteractor.h" // itk includes #include "itksys/SystemTools.hxx" #include #include const std::string QmitkImageStatisticsView::VIEW_ID = "org.mitk.views.imagestatistics"; const int QmitkImageStatisticsView::STAT_TABLE_BASE_HEIGHT = 180; QmitkImageStatisticsView::QmitkImageStatisticsView(QObject* /*parent*/, const char* /*name*/) : m_Controls( NULL ), m_TimeStepperAdapter( NULL ), m_SelectedImage( NULL ), m_SelectedImageMask( NULL ), m_SelectedPlanarFigure( NULL ), m_ImageObserverTag( -1 ), m_ImageMaskObserverTag( -1 ), m_PlanarFigureObserverTag( -1 ), m_TimeObserverTag( -1 ), m_CurrentStatisticsValid( false ), m_StatisticsUpdatePending( false ), m_DataNodeSelectionChanged ( false ), m_Visible(false) { this->m_CalculationThread = new QmitkImageStatisticsCalculationThread; } QmitkImageStatisticsView::~QmitkImageStatisticsView() { if ( m_SelectedImage != NULL ) m_SelectedImage->RemoveObserver( m_ImageObserverTag ); if ( m_SelectedImageMask != NULL ) m_SelectedImageMask->RemoveObserver( m_ImageMaskObserverTag ); if ( m_SelectedPlanarFigure != NULL ) m_SelectedPlanarFigure->RemoveObserver( m_PlanarFigureObserverTag ); while(this->m_CalculationThread->isRunning()) // wait until thread has finished { itksys::SystemTools::Delay(100); } delete this->m_CalculationThread; } void QmitkImageStatisticsView::CreateQtPartControl(QWidget *parent) { if (m_Controls == NULL) { m_Controls = new Ui::QmitkImageStatisticsViewControls; m_Controls->setupUi(parent); CreateConnections(); m_Controls->m_ErrorMessageLabel->hide(); m_Controls->m_StatisticsWidgetStack->setCurrentIndex( 0 ); m_Controls->m_BinSizeFrame->setVisible(false); } } void QmitkImageStatisticsView::CreateConnections() { if ( m_Controls ) { connect( (QObject*)(this->m_Controls->m_ButtonCopyHistogramToClipboard), SIGNAL(clicked()),(QObject*) this, SLOT(OnClipboardHistogramButtonClicked()) ); connect( (QObject*)(this->m_Controls->m_ButtonCopyStatisticsToClipboard), SIGNAL(clicked()),(QObject*) this, SLOT(OnClipboardStatisticsButtonClicked()) ); connect( (QObject*)(this->m_Controls->m_IgnoreZerosCheckbox), SIGNAL(clicked()),(QObject*) this, SLOT(OnIgnoreZerosCheckboxClicked()) ); connect( (QObject*) this->m_CalculationThread, SIGNAL(finished()),this, SLOT( OnThreadedStatisticsCalculationEnds()),Qt::QueuedConnection); connect( (QObject*) this, SIGNAL(StatisticsUpdate()),this, SLOT( RequestStatisticsUpdate()), Qt::QueuedConnection); connect( (QObject*) this->m_Controls->m_StatisticsTable, SIGNAL(cellDoubleClicked(int,int)),this, SLOT( JumpToCoordinates(int,int)) ); connect( (QObject*) (this->m_Controls->m_barRadioButton), SIGNAL(clicked()), (QObject*) (this->m_Controls->m_JSHistogram), SLOT(OnBarRadioButtonSelected())); connect( (QObject*) (this->m_Controls->m_lineRadioButton), SIGNAL(clicked()), (QObject*) (this->m_Controls->m_JSHistogram), SLOT(OnLineRadioButtonSelected())); connect( (QObject*) (this->m_Controls->m_HistogramBinSizeSpinbox), SIGNAL(editingFinished()), this, SLOT(OnHistogramBinSizeBoxValueChanged())); connect( (QObject*)(this->m_Controls->m_UseDefaultBinSizeBox), SIGNAL(clicked()),(QObject*) this, SLOT(OnDefaultBinSizeBoxChanged()) ); } } void QmitkImageStatisticsView::OnDefaultBinSizeBoxChanged() { if (m_CalculationThread!=NULL) m_Controls->m_HistogramBinSizeSpinbox->setValue(m_CalculationThread->GetHistogramBinSize()); if (m_Controls->m_UseDefaultBinSizeBox->isChecked()) m_Controls->m_BinSizeFrame->setVisible(false); else m_Controls->m_BinSizeFrame->setVisible(true); } void QmitkImageStatisticsView::PartClosed( berry::IWorkbenchPartReference::Pointer ) { } void QmitkImageStatisticsView::OnTimeChanged(const itk::EventObject& e) { if (this->m_SelectedDataNodes.isEmpty() || this->m_SelectedImage == NULL) return; const mitk::SliceNavigationController::GeometryTimeEvent* timeEvent = dynamic_cast(&e); assert(timeEvent != NULL); unsigned int timestep = timeEvent->GetPos(); if (this->m_SelectedImage->GetTimeSteps() > 1) { for (unsigned int x = 0; x < this->m_Controls->m_StatisticsTable->columnCount(); x++) { for (unsigned int y = 0; y < this->m_Controls->m_StatisticsTable->rowCount(); y++) { QTableWidgetItem* item = this->m_Controls->m_StatisticsTable->item(y, x); if (item == NULL) break; if (x == timestep) { item->setBackgroundColor(Qt::yellow); } else { if (y % 2 == 0) item->setBackground(this->m_Controls->m_StatisticsTable->palette().base()); else item->setBackground(this->m_Controls->m_StatisticsTable->palette().alternateBase()); } } } this->m_Controls->m_StatisticsTable->viewport()->update(); } if ((this->m_SelectedImage->GetTimeSteps() == 1 && timestep == 0) || this->m_SelectedImage->GetTimeSteps() > 1) { // display histogram for selected timestep this->m_Controls->m_JSHistogram->ClearHistogram(); QmitkImageStatisticsCalculationThread::HistogramType::Pointer histogram = this->m_CalculationThread->GetTimeStepHistogram(timestep); if (histogram.IsNotNull()) { this->m_Controls->m_JSHistogram->ComputeHistogram(histogram.GetPointer()); // this->m_Controls->m_JSHistogram->SignalGraphChanged(); // hacky way to make sure the protected SignalGraphChanged() is called if (this->m_Controls->m_JSHistogram->GetUseLineGraph()) { this->m_Controls->m_JSHistogram->OnBarRadioButtonSelected(); this->m_Controls->m_JSHistogram->OnLineRadioButtonSelected(); } else { this->m_Controls->m_JSHistogram->OnLineRadioButtonSelected(); this->m_Controls->m_JSHistogram->OnBarRadioButtonSelected(); } } } } void QmitkImageStatisticsView::JumpToCoordinates(int row ,int col) { if(m_SelectedDataNodes.isEmpty()) { MITK_WARN("QmitkImageStatisticsView") << "No data node selected for statistics calculation." ; return; } mitk::Point3D world; if (row==4 && !m_WorldMinList.empty()) world = m_WorldMinList[col]; else if (row==3 && !m_WorldMaxList.empty()) world = m_WorldMaxList[col]; else return; mitk::IRenderWindowPart* part = this->GetRenderWindowPart(); if (part) { part->GetQmitkRenderWindow("axial")->GetSliceNavigationController()->SelectSliceByPoint(world); part->GetQmitkRenderWindow("sagittal")->GetSliceNavigationController()->SelectSliceByPoint(world); part->GetQmitkRenderWindow("coronal")->GetSliceNavigationController()->SelectSliceByPoint(world); mitk::SliceNavigationController::GeometryTimeEvent timeEvent(this->m_SelectedImage->GetTimeGeometry(), col); part->GetQmitkRenderWindow("axial")->GetSliceNavigationController()->SetGeometryTime(timeEvent); } } void QmitkImageStatisticsView::OnIgnoreZerosCheckboxClicked() { emit StatisticsUpdate(); } void QmitkImageStatisticsView::OnClipboardHistogramButtonClicked() { if ( m_CurrentStatisticsValid ) { const unsigned int t = this->GetRenderWindowPart()->GetTimeNavigationController()->GetTime()->GetPos(); typedef mitk::ImageStatisticsCalculator::HistogramType HistogramType; const HistogramType *histogram = this->m_CalculationThread->GetTimeStepHistogram(t).GetPointer(); QString clipboard( "Measurement \t Frequency\n" ); for ( HistogramType::ConstIterator it = histogram->Begin(); it != histogram->End(); ++it ) { - if( m_Controls->m_HistogramBinSizeSpinbox->value() == 1) + if( m_Controls->m_HistogramBinSizeSpinbox->value() == 1.0) { clipboard = clipboard.append( "%L1 \t %L2\n" ) .arg( it.GetMeasurementVector()[0], 0, 'f', 0 ) .arg( it.GetFrequency() ); } else { clipboard = clipboard.append( "%L1 \t %L2\n" ) .arg( it.GetMeasurementVector()[0], 0, 'f', 2 ) .arg( it.GetFrequency() ); } } QApplication::clipboard()->setText( clipboard, QClipboard::Clipboard ); } else { QApplication::clipboard()->clear(); } } void QmitkImageStatisticsView::OnClipboardStatisticsButtonClicked() { QLocale tempLocal; QLocale::setDefault(QLocale(QLocale::English, QLocale::UnitedStates)); if ( this->m_CurrentStatisticsValid ) { const std::vector &statistics = this->m_CalculationThread->GetStatisticsData(); const unsigned int t = this->GetRenderWindowPart()->GetTimeNavigationController()->GetTime()-> GetPos(); // Copy statistics to clipboard ("%Ln" will use the default locale for // number formatting) QString clipboard( "Mean \t StdDev \t RMS \t Max \t Min \t N \t V (mm³)\n" ); clipboard = clipboard.append( "%L1 \t %L2 \t %L3 \t %L4 \t %L5 \t %L6 \t %L7" ) .arg( statistics[t].GetMean(), 0, 'f', 10 ) .arg( statistics[t].GetSigma(), 0, 'f', 10 ) .arg( statistics[t].GetRMS(), 0, 'f', 10 ) .arg( statistics[t].GetMax(), 0, 'f', 10 ) .arg( statistics[t].GetMin(), 0, 'f', 10 ) .arg( statistics[t].GetN() ) .arg( m_Controls->m_StatisticsTable->item( 0, 6 )->text().toDouble(), 0, 'f', 10 ); QApplication::clipboard()->setText( clipboard, QClipboard::Clipboard ); } else { QApplication::clipboard()->clear(); } QLocale::setDefault(tempLocal); } void QmitkImageStatisticsView::OnSelectionChanged( berry::IWorkbenchPart::Pointer /*part*/, const QList &selectedNodes ) { if (this->m_Visible) { this->SelectionChanged( selectedNodes ); } else { this->m_DataNodeSelectionChanged = true; } } void QmitkImageStatisticsView::SelectionChanged(const QList &selectedNodes) { if( this->m_StatisticsUpdatePending ) { this->m_DataNodeSelectionChanged = true; return; // not ready for new data now! } if (selectedNodes.size() == this->m_SelectedDataNodes.size()) { int i = 0; for (; i < selectedNodes.size(); ++i) { if (selectedNodes.at(i) != this->m_SelectedDataNodes.at(i)) { break; } } // node selection did not change if (i == selectedNodes.size()) return; } //reset the feature image and image mask field m_Controls->m_SelectedFeatureImageLabel->setText("None"); m_Controls->m_SelectedMaskLabel->setText("None"); this->ReinitData(); if (selectedNodes.isEmpty()) { m_Controls->m_JSHistogram->ClearHistogram(); m_Controls->m_lineRadioButton->setEnabled(true); m_Controls->m_barRadioButton->setEnabled(true); m_Controls->m_HistogramBinSizeSpinbox->setEnabled(true); m_Controls->m_HistogramBinSizeCaptionLabel->setEnabled(true); // m_Controls->m_HistogramBinSizeLabel->setEnabled(true); m_Controls->m_InfoLabel->setText(QString("")); // m_Controls->horizontalLayout_3->setEnabled(false); m_Controls->groupBox->setEnabled(false); m_Controls->groupBox_3->setEnabled(false); } else { // m_Controls->horizontalLayout_3->setEnabled(true); m_Controls->groupBox->setEnabled(true); m_Controls->groupBox_3->setEnabled(true); } if(selectedNodes.size() == 1 || selectedNodes.size() == 2) { bool isBinary = false; selectedNodes.value(0)->GetBoolProperty("binary",isBinary); if(isBinary) { m_Controls->m_JSHistogram->ClearHistogram(); m_Controls->m_lineRadioButton->setEnabled(true); m_Controls->m_barRadioButton->setEnabled(true); m_Controls->m_HistogramBinSizeSpinbox->setEnabled(true); m_Controls->m_HistogramBinSizeCaptionLabel->setEnabled(true); // m_Controls->m_HistogramBinSizeLabel->setEnabled(true); m_Controls->m_InfoLabel->setText(QString("")); } for (int i= 0; i< selectedNodes.size(); ++i) { this->m_SelectedDataNodes.push_back(selectedNodes.at(i)); } this->m_DataNodeSelectionChanged = false; this->m_Controls->m_ErrorMessageLabel->setText( "" ); this->m_Controls->m_ErrorMessageLabel->hide(); emit StatisticsUpdate(); } else { this->m_DataNodeSelectionChanged = false; } } void QmitkImageStatisticsView::ReinitData() { while( this->m_CalculationThread->isRunning()) // wait until thread has finished { itksys::SystemTools::Delay(100); } if(this->m_SelectedImage != NULL) { this->m_SelectedImage->RemoveObserver( this->m_ImageObserverTag); this->m_SelectedImage = NULL; } if(this->m_SelectedImageMask != NULL) { this->m_SelectedImageMask->RemoveObserver( this->m_ImageMaskObserverTag); this->m_SelectedImageMask = NULL; } if(this->m_SelectedPlanarFigure != NULL) { this->m_SelectedPlanarFigure->RemoveObserver( this->m_PlanarFigureObserverTag); this->m_SelectedPlanarFigure = NULL; } this->m_SelectedDataNodes.clear(); this->m_StatisticsUpdatePending = false; m_Controls->m_ErrorMessageLabel->setText( "" ); m_Controls->m_ErrorMessageLabel->hide(); this->InvalidateStatisticsTableView(); m_Controls->m_JSHistogram->ClearHistogram(); m_Controls->m_StatisticsWidgetStack->setCurrentIndex( 0 ); } void QmitkImageStatisticsView::OnThreadedStatisticsCalculationEnds() { std::stringstream message; message << ""; m_Controls->m_ErrorMessageLabel->setText( message.str().c_str() ); m_Controls->m_ErrorMessageLabel->hide(); this->WriteStatisticsToGUI(); } void QmitkImageStatisticsView::UpdateStatistics() { mitk::IRenderWindowPart* renderPart = this->GetRenderWindowPart(); if ( renderPart == NULL ) { this->m_StatisticsUpdatePending = false; return; } m_WorldMinList.clear(); m_WorldMaxList.clear(); // classify selected nodes mitk::NodePredicateDataType::Pointer imagePredicate = mitk::NodePredicateDataType::New("Image"); std::string maskName = std::string(); std::string maskType = std::string(); std::string featureImageName = std::string(); unsigned int maskDimension = 0; // reset data from last run ITKCommandType::Pointer changeListener = ITKCommandType::New(); changeListener->SetCallbackFunction( this, &QmitkImageStatisticsView::SelectedDataModified ); mitk::DataNode::Pointer planarFigureNode; for( int i= 0 ; i < this->m_SelectedDataNodes.size(); ++i) { mitk::PlanarFigure::Pointer planarFig = dynamic_cast(this->m_SelectedDataNodes.at(i)->GetData()); if( imagePredicate->CheckNode(this->m_SelectedDataNodes.at(i)) ) { bool isMask = false; this->m_SelectedDataNodes.at(i)->GetPropertyValue("binary", isMask); if( this->m_SelectedImageMask == NULL && isMask) { this->m_SelectedImageMask = dynamic_cast(this->m_SelectedDataNodes.at(i)->GetData()); this->m_ImageMaskObserverTag = this->m_SelectedImageMask->AddObserver(itk::ModifiedEvent(), changeListener); maskName = this->m_SelectedDataNodes.at(i)->GetName(); maskType = m_SelectedImageMask->GetNameOfClass(); maskDimension = 3; } else if( !isMask ) { if(this->m_SelectedImage == NULL) { this->m_SelectedImage = static_cast(this->m_SelectedDataNodes.at(i)->GetData()); this->m_ImageObserverTag = this->m_SelectedImage->AddObserver(itk::ModifiedEvent(), changeListener); } featureImageName = this->m_SelectedDataNodes.at(i)->GetName(); } } else if (planarFig.IsNotNull()) { if(this->m_SelectedPlanarFigure == NULL) { this->m_SelectedPlanarFigure = planarFig; this->m_PlanarFigureObserverTag = this->m_SelectedPlanarFigure->AddObserver(mitk::EndInteractionPlanarFigureEvent(), changeListener); maskName = this->m_SelectedDataNodes.at(i)->GetName(); maskType = this->m_SelectedPlanarFigure->GetNameOfClass(); maskDimension = 2; planarFigureNode = m_SelectedDataNodes.at(i); } } else { std::stringstream message; message << "" << "Invalid data node type!" << ""; m_Controls->m_ErrorMessageLabel->setText( message.str().c_str() ); m_Controls->m_ErrorMessageLabel->show(); } } if(maskName == "") { maskName = "None"; maskType = ""; maskDimension = 0; } if(featureImageName == "") { featureImageName = "None"; } if (m_SelectedPlanarFigure != NULL && m_SelectedImage == NULL) { mitk::DataStorage::SetOfObjects::ConstPointer parentSet = this->GetDataStorage()->GetSources(planarFigureNode); for (int i=0; iSize(); i++) { mitk::DataNode::Pointer node = parentSet->ElementAt(i); if( imagePredicate->CheckNode(node) ) { bool isMask = false; node->GetPropertyValue("binary", isMask); if( !isMask ) { if(this->m_SelectedImage == NULL) { this->m_SelectedImage = static_cast(node->GetData()); this->m_ImageObserverTag = this->m_SelectedImage->AddObserver(itk::ModifiedEvent(), changeListener); } } } } } unsigned int timeStep = renderPart->GetTimeNavigationController()->GetTime()->GetPos(); if ( m_SelectedImage != NULL && m_SelectedImage->IsInitialized()) { // Check if a the selected image is a multi-channel image. If yes, statistics // cannot be calculated currently. if ( m_SelectedImage->GetPixelType().GetNumberOfComponents() > 1 ) { std::stringstream message; message << "Multi-component images not supported."; m_Controls->m_ErrorMessageLabel->setText( message.str().c_str() ); m_Controls->m_ErrorMessageLabel->show(); this->InvalidateStatisticsTableView(); m_Controls->m_StatisticsWidgetStack->setCurrentIndex( 0 ); m_Controls->m_JSHistogram->ClearHistogram(); m_CurrentStatisticsValid = false; this->m_StatisticsUpdatePending = false; m_Controls->m_lineRadioButton->setEnabled(true); m_Controls->m_barRadioButton->setEnabled(true); m_Controls->m_HistogramBinSizeSpinbox->setEnabled(true); m_Controls->m_HistogramBinSizeCaptionLabel->setEnabled(true); // m_Controls->m_HistogramBinSizeLabel->setEnabled(true); m_Controls->m_InfoLabel->setText(QString("")); return; } std::stringstream maskLabel; maskLabel << maskName; if ( maskDimension > 0 ) { maskLabel << " [" << maskDimension << "D " << maskType << "]"; } m_Controls->m_SelectedMaskLabel->setText( maskLabel.str().c_str() ); m_Controls->m_SelectedFeatureImageLabel->setText(featureImageName.c_str()); // check time step validity if(m_SelectedImage->GetDimension() <= 3 && timeStep > m_SelectedImage->GetDimension(3)-1) { timeStep = m_SelectedImage->GetDimension(3)-1; } // Add the used mask time step to the mask label so the user knows which mask time step was used // if the image time step is bigger than the total number of mask time steps (see // ImageStatisticsCalculator::ExtractImageAndMask) if (m_SelectedImageMask != NULL) { unsigned int maskTimeStep = timeStep; if (maskTimeStep >= m_SelectedImageMask->GetTimeSteps()) { maskTimeStep = m_SelectedImageMask->GetTimeSteps() - 1; } m_Controls->m_SelectedMaskLabel->setText(m_Controls->m_SelectedMaskLabel->text() + QString(" (t=") + QString::number(maskTimeStep) + QString(")")); } //// initialize thread and trigger it this->m_CalculationThread->SetIgnoreZeroValueVoxel( m_Controls->m_IgnoreZerosCheckbox->isChecked() ); this->m_CalculationThread->Initialize( m_SelectedImage, m_SelectedImageMask, m_SelectedPlanarFigure ); this->m_CalculationThread->SetTimeStep( timeStep ); this->m_CalculationThread->SetHistogramBinSize(m_Controls->m_HistogramBinSizeSpinbox->value()); std::stringstream message; message << "Calculating statistics..."; m_Controls->m_ErrorMessageLabel->setText( message.str().c_str() ); m_Controls->m_ErrorMessageLabel->show(); try { // Compute statistics this->m_CalculationThread->SetUseDefaultBinSize(m_Controls->m_UseDefaultBinSizeBox->isChecked()); this->m_CalculationThread->start(); } catch ( const mitk::Exception& e) { std::stringstream message; message << "" << e.GetDescription() << ""; m_Controls->m_ErrorMessageLabel->setText( message.str().c_str() ); m_Controls->m_ErrorMessageLabel->show(); this->m_StatisticsUpdatePending = false; } catch ( const std::runtime_error &e ) { // In case of exception, print error message on GUI std::stringstream message; message << "" << e.what() << ""; m_Controls->m_ErrorMessageLabel->setText( message.str().c_str() ); m_Controls->m_ErrorMessageLabel->show(); this->m_StatisticsUpdatePending = false; } catch ( const std::exception &e ) { MITK_ERROR << "Caught exception: " << e.what(); // In case of exception, print error message on GUI std::stringstream message; message << "Error! Unequal Dimensions of Image and Segmentation. No recompute possible "; m_Controls->m_ErrorMessageLabel->setText( message.str().c_str() ); m_Controls->m_ErrorMessageLabel->show(); this->m_StatisticsUpdatePending = false; } } else { this->m_StatisticsUpdatePending = false; } } void QmitkImageStatisticsView::SelectedDataModified() { if( !m_StatisticsUpdatePending ) { emit StatisticsUpdate(); } } void QmitkImageStatisticsView::NodeRemoved(const mitk::DataNode *node) { while(this->m_CalculationThread->isRunning()) // wait until thread has finished { itksys::SystemTools::Delay(100); } if (node->GetData() == m_SelectedImage) { m_SelectedImage = NULL; } } void QmitkImageStatisticsView::RequestStatisticsUpdate() { if ( !m_StatisticsUpdatePending ) { if(this->m_DataNodeSelectionChanged) { this->SelectionChanged(this->GetCurrentSelection()); } else { this->m_StatisticsUpdatePending = true; this->UpdateStatistics(); } } if (this->GetRenderWindowPart()) this->GetRenderWindowPart()->RequestUpdate(); } void QmitkImageStatisticsView::OnHistogramBinSizeBoxValueChanged() { this->UpdateStatistics(); } void QmitkImageStatisticsView::WriteStatisticsToGUI() { m_Controls->m_lineRadioButton->setEnabled(true); m_Controls->m_barRadioButton->setEnabled(true); m_Controls->m_HistogramBinSizeSpinbox->setEnabled(true); m_Controls->m_HistogramBinSizeCaptionLabel->setEnabled(true); // m_Controls->m_HistogramBinSizeLabel->setEnabled(true); m_Controls->m_InfoLabel->setText(QString("")); if(m_DataNodeSelectionChanged) { this->m_StatisticsUpdatePending = false; this->RequestStatisticsUpdate(); return; // stop visualization of results and calculate statistics of new selection } if ( this->m_CalculationThread->GetStatisticsUpdateSuccessFlag()) { if ( this->m_CalculationThread->GetStatisticsChangedFlag() ) { // Do not show any error messages m_Controls->m_ErrorMessageLabel->hide(); m_CurrentStatisticsValid = true; } if (m_Controls->m_barRadioButton->isChecked()) { m_Controls->m_JSHistogram->OnBarRadioButtonSelected(); } m_Controls->m_StatisticsWidgetStack->setCurrentIndex( 0 ); m_Controls->m_HistogramBinSizeSpinbox->setValue( this->m_CalculationThread->GetHistogramBinSize() ); //m_Controls->m_JSHistogram->ComputeHistogram( this->m_CalculationThread->GetTimeStepHistogram(this->m_CalculationThread->GetTimeStep()).GetPointer() ); this->FillStatisticsTableView( this->m_CalculationThread->GetStatisticsData(), this->m_CalculationThread->GetStatisticsImage()); } else { m_Controls->m_SelectedMaskLabel->setText( "None" ); m_Controls->m_ErrorMessageLabel->setText( m_CalculationThread->GetLastErrorMessage().c_str() ); m_Controls->m_ErrorMessageLabel->show(); // Clear statistics and histogram this->InvalidateStatisticsTableView(); m_Controls->m_StatisticsWidgetStack->setCurrentIndex( 0 ); //m_Controls->m_JSHistogram->clearHistogram(); m_CurrentStatisticsValid = false; // If a (non-closed) PlanarFigure is selected, display a line profile widget if ( m_SelectedPlanarFigure != NULL ) { // Check if the (closed) planar figure is out of bounds and so no image mask could be calculated--> Intensity Profile can not be calculated bool outOfBounds = false; if ( m_SelectedPlanarFigure->IsClosed() && m_SelectedImageMask == NULL) { outOfBounds = true; std::stringstream message; message << "Planar figure is on a rotated image plane or outside the image bounds."; m_Controls->m_InfoLabel->setText(message.str().c_str()); } // check whether PlanarFigure is initialized const mitk::PlaneGeometry *planarFigurePlaneGeometry = m_SelectedPlanarFigure->GetPlaneGeometry(); if ( planarFigurePlaneGeometry == NULL || outOfBounds) { // Clear statistics, histogram, and GUI this->InvalidateStatisticsTableView(); m_Controls->m_StatisticsWidgetStack->setCurrentIndex( 0 ); m_Controls->m_JSHistogram->ClearHistogram(); m_CurrentStatisticsValid = false; m_Controls->m_ErrorMessageLabel->hide(); m_Controls->m_SelectedMaskLabel->setText( "None" ); this->m_StatisticsUpdatePending = false; m_Controls->m_lineRadioButton->setEnabled(true); m_Controls->m_barRadioButton->setEnabled(true); m_Controls->m_HistogramBinSizeSpinbox->setEnabled(true); m_Controls->m_HistogramBinSizeCaptionLabel->setEnabled(true); // m_Controls->m_HistogramBinSizeLabel->setEnabled(true); if (!outOfBounds) m_Controls->m_InfoLabel->setText(QString("")); return; } unsigned int timeStep = this->GetRenderWindowPart()->GetTimeNavigationController()->GetTime()->GetPos(); m_Controls->m_JSHistogram->SetImage(this->m_CalculationThread->GetStatisticsImage()); m_Controls->m_JSHistogram->SetPlanarFigure(m_SelectedPlanarFigure); m_Controls->m_JSHistogram->ComputeIntensityProfile(timeStep); m_Controls->m_lineRadioButton->setEnabled(false); m_Controls->m_barRadioButton->setEnabled(false); m_Controls->m_HistogramBinSizeSpinbox->setEnabled(false); m_Controls->m_HistogramBinSizeCaptionLabel->setEnabled(false); // m_Controls->m_HistogramBinSizeLabel->setEnabled(false); std::stringstream message; message << "Only linegraph available for an intensity profile!"; m_Controls->m_InfoLabel->setText(message.str().c_str()); } } this->m_StatisticsUpdatePending = false; } void QmitkImageStatisticsView::FillStatisticsTableView( const std::vector &s, const mitk::Image *image ) { this->m_Controls->m_StatisticsTable->setColumnCount(image->GetTimeSteps()); this->m_Controls->m_StatisticsTable->horizontalHeader()->setVisible(image->GetTimeSteps() > 1); int decimals = 2; mitk::PixelType doublePix = mitk::MakeScalarPixelType< double >(); mitk::PixelType floatPix = mitk::MakeScalarPixelType< float >(); if (image->GetPixelType()==doublePix || image->GetPixelType()==floatPix) { decimals = 5; } for (unsigned int t = 0; t < image->GetTimeSteps(); t++) { this->m_Controls->m_StatisticsTable->setHorizontalHeaderItem(t, new QTableWidgetItem(QString::number(t))); if (s[t].GetMaxIndex().size()==3) { mitk::Point3D index, max, min; index[0] = s[t].GetMaxIndex()[0]; index[1] = s[t].GetMaxIndex()[1]; index[2] = s[t].GetMaxIndex()[2]; m_SelectedImage->GetGeometry()->IndexToWorld(index, max); this->m_WorldMaxList.push_back(max); index[0] = s[t].GetMinIndex()[0]; index[1] = s[t].GetMinIndex()[1]; index[2] = s[t].GetMinIndex()[2]; m_SelectedImage->GetGeometry()->IndexToWorld(index, min); this->m_WorldMinList.push_back(min); } this->m_Controls->m_StatisticsTable->setItem( 0, t, new QTableWidgetItem( QString("%1").arg(s[t].GetMean(), 0, 'f', decimals) ) ); this->m_Controls->m_StatisticsTable->setItem( 1, t, new QTableWidgetItem( QString("%1").arg(s[t].GetSigma(), 0, 'f', decimals) ) ); this->m_Controls->m_StatisticsTable->setItem( 2, t, new QTableWidgetItem( QString("%1").arg(s[t].GetRMS(), 0, 'f', decimals) ) ); QString max; max.append(QString("%1").arg(s[t].GetMax(), 0, 'f', decimals)); max += " ("; for (int i=0; im_Controls->m_StatisticsTable->setItem( 3, t, new QTableWidgetItem( max ) ); QString min; min.append(QString("%1").arg(s[t].GetMin(), 0, 'f', decimals)); min += " ("; for (int i=0; im_Controls->m_StatisticsTable->setItem( 4, t, new QTableWidgetItem( min ) ); this->m_Controls->m_StatisticsTable->setItem( 5, t, new QTableWidgetItem( QString("%1").arg(s[t].GetN()) ) ); const mitk::BaseGeometry *geometry = image->GetGeometry(); if ( geometry != NULL ) { const mitk::Vector3D &spacing = image->GetGeometry()->GetSpacing(); double volume = spacing[0] * spacing[1] * spacing[2] * (double) s[t].GetN(); this->m_Controls->m_StatisticsTable->setItem( 6, t, new QTableWidgetItem( QString("%1").arg(volume, 0, 'f', decimals) ) ); } else { this->m_Controls->m_StatisticsTable->setItem( 6, t, new QTableWidgetItem( "NA" ) ); } } this->m_Controls->m_StatisticsTable->resizeColumnsToContents(); int height = STAT_TABLE_BASE_HEIGHT; if (this->m_Controls->m_StatisticsTable->horizontalHeader()->isVisible()) height += this->m_Controls->m_StatisticsTable->horizontalHeader()->height(); if (this->m_Controls->m_StatisticsTable->horizontalScrollBar()->isVisible()) height += this->m_Controls->m_StatisticsTable->horizontalScrollBar()->height(); this->m_Controls->m_StatisticsTable->setMinimumHeight(height); // make sure the current timestep's column is highlighted (and the correct histogram is displayed) unsigned int t = this->GetRenderWindowPart()->GetTimeNavigationController()->GetTime()-> GetPos(); mitk::SliceNavigationController::GeometryTimeEvent timeEvent(this->m_SelectedImage->GetTimeGeometry(), t); this->OnTimeChanged(timeEvent); t = std::min(image->GetTimeSteps() - 1, t); QString hotspotMean; hotspotMean.append(QString("%1").arg(s[t].GetHotspotStatistics().GetMean(), 0, 'f', decimals)); hotspotMean += " ("; for (int i=0; im_Controls->m_StatisticsTable->setItem( 7, t, new QTableWidgetItem( hotspotMean ) ); QString hotspotMax; hotspotMax.append(QString("%1").arg(s[t].GetHotspotStatistics().GetMax(), 0, 'f', decimals)); hotspotMax += " ("; for (int i=0; im_Controls->m_StatisticsTable->setItem( 8, t, new QTableWidgetItem( hotspotMax ) ); QString hotspotMin; hotspotMin.append(QString("%1").arg(s[t].GetHotspotStatistics().GetMin(), 0, 'f', decimals)); hotspotMin += " ("; for (int i=0; im_Controls->m_StatisticsTable->setItem( 9, t, new QTableWidgetItem( hotspotMin ) ); } void QmitkImageStatisticsView::InvalidateStatisticsTableView() { this->m_Controls->m_StatisticsTable->horizontalHeader()->setVisible(false); this->m_Controls->m_StatisticsTable->setColumnCount(1); for ( unsigned int i = 0; i < this->m_Controls->m_StatisticsTable->rowCount(); ++i ) { { this->m_Controls->m_StatisticsTable->setItem( i, 0, new QTableWidgetItem( "NA" ) ); } } this->m_Controls->m_StatisticsTable->setMinimumHeight(STAT_TABLE_BASE_HEIGHT); } void QmitkImageStatisticsView::Activated() { } void QmitkImageStatisticsView::Deactivated() { } void QmitkImageStatisticsView::Visible() { m_Visible = true; mitk::IRenderWindowPart* renderWindow = GetRenderWindowPart(); if (renderWindow) { itk::ReceptorMemberCommand::Pointer cmdTimeEvent = itk::ReceptorMemberCommand::New(); cmdTimeEvent->SetCallbackFunction(this, &QmitkImageStatisticsView::OnTimeChanged); // It is sufficient to add the observer to the axial render window since the GeometryTimeEvent // is always triggered by all views. m_TimeObserverTag = renderWindow->GetQmitkRenderWindow("axial")-> GetSliceNavigationController()-> AddObserver(mitk::SliceNavigationController::GeometryTimeEvent(NULL, 0), cmdTimeEvent); } if (m_DataNodeSelectionChanged) { if (this->IsCurrentSelectionValid()) { this->SelectionChanged(this->GetCurrentSelection()); } else { this->SelectionChanged(this->GetDataManagerSelection()); } m_DataNodeSelectionChanged = false; } } void QmitkImageStatisticsView::Hidden() { m_Visible = false; // The slice navigation controller observer is removed here instead of in the destructor. // If it was called in the destructor, the application would freeze because the view's // destructor gets called after the render windows have been destructed. if ( m_TimeObserverTag != NULL ) { mitk::IRenderWindowPart* renderWindow = GetRenderWindowPart(); if (renderWindow) { renderWindow->GetQmitkRenderWindow("axial")->GetSliceNavigationController()-> RemoveObserver( m_TimeObserverTag ); } m_TimeObserverTag = NULL; } } void QmitkImageStatisticsView::SetFocus() { } diff --git a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsViewControls.ui b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsViewControls.ui index f9597506b9..059afa50e2 100644 --- a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsViewControls.ui +++ b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsViewControls.ui @@ -1,546 +1,552 @@ QmitkImageStatisticsViewControls true 0 0 - 281 + 548 800 Form 0 0 Qt::LeftToRight Feature Image: 0 0 None 0 0 Mask: 0 0 None -1 0 0 color: rgb(255, 0, 0); Error Message Qt::AutoText Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter Ignore zero-valued voxels false Statistics 9 9 9 100 180 16777215 16777215 Qt::ScrollBarAlwaysOff Qt::ScrollBarAsNeeded true QAbstractItemView::NoEditTriggers true true Qt::DotLine false 10 false false 80 true 80 false true true false 25 25 false false Mean StdDev RMS Max Min N V (mm³) Hotspot Peak Hotspot Max Hotspot Min 0 0 0 Copy to Clipboard Qt::Horizontal 40 20 false 150 160 Histogram false 0 0 0 0 16777215 16777215 Plot 0 0 Barchart true 0 0 0 0 Linegraph Qt::Horizontal 40 20 Use default bin size true QFrame::NoFrame QFrame::Raised 0 0 - - - - Press enter to recalculate statistics with new bin size. - - - true - - - 1 - - - 10 - - - 60 0 - 40 + 100 16777215 Bin size: + + + + Press enter to recalculate statistics with new bin size. + + + true + + + 5 + + + 0.000010000000000 + + + 1000000.000000000000000 + + + 10.000000000000000 + + + - - - - - - - 0 0 0 + + + + + + + 0 0 Copy to Clipboard Qt::Horizontal 40 20 Qt::Vertical 20 40 QmitkHistogramJSWidget QWidget
QmitkHistogramJSWidget.h
1