diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp index a0975aa7e5..0162c0858d 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp @@ -1,1256 +1,1256 @@ /*=================================================================== 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 "mitkPartialVolumeAnalysisHistogramCalculator.h" #include "mitkImageAccessByItk.h" #include "mitkExtractImageFilter.h" #include #include #include #include #include #include "itkResampleImageFilter.h" #include "itkGaussianInterpolateImageFunction.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "itkGaussianInterpolateImageFunction.h" #include "itkBSplineInterpolateImageFunction.h" #include "itkNearestNeighborInterpolateImageFunction.h" #include "itkImageMaskSpatialObject.h" #include "itkRegionOfInterestImageFilter.h" #include "itkListSample.h" #include #include namespace mitk { PartialVolumeAnalysisHistogramCalculator::PartialVolumeAnalysisHistogramCalculator() : m_MaskingMode( MASKING_MODE_NONE ), m_MaskingModeChanged( false ), m_NumberOfBins(256), m_UpsamplingFactor(1), m_GaussianSigma(0), m_ForceUpdate(false), m_PlanarFigureThickness(0) { m_EmptyHistogram = HistogramType::New(); m_EmptyHistogram->SetMeasurementVectorSize(1); HistogramType::SizeType histogramSize(1); histogramSize.Fill( 256 ); m_EmptyHistogram->Initialize( histogramSize ); m_EmptyStatistics.Reset(); } PartialVolumeAnalysisHistogramCalculator::~PartialVolumeAnalysisHistogramCalculator() { } void PartialVolumeAnalysisHistogramCalculator::SetImage( const mitk::Image *image ) { if ( m_Image != image ) { m_Image = image; this->Modified(); m_ImageStatisticsTimeStamp.Modified(); m_ImageStatisticsCalculationTriggerBool = true; } } void PartialVolumeAnalysisHistogramCalculator::AddAdditionalResamplingImage( const mitk::Image *image ) { m_AdditionalResamplingImages.push_back(image); this->Modified(); m_ImageStatisticsTimeStamp.Modified(); m_ImageStatisticsCalculationTriggerBool = true; } void PartialVolumeAnalysisHistogramCalculator::SetModified( ) { this->Modified(); m_ImageStatisticsTimeStamp.Modified(); m_ImageStatisticsCalculationTriggerBool = true; m_MaskedImageStatisticsTimeStamp.Modified(); m_MaskedImageStatisticsCalculationTriggerBool = true; m_PlanarFigureStatisticsTimeStamp.Modified(); m_PlanarFigureStatisticsCalculationTriggerBool = true; } void PartialVolumeAnalysisHistogramCalculator::SetImageMask( const mitk::Image *imageMask ) { if ( m_Image.IsNull() ) { itkExceptionMacro( << "Image needs to be set first!" ); } if ( m_Image->GetTimeSteps() != imageMask->GetTimeSteps() ) { itkExceptionMacro( << "Image and image mask need to have equal number of time steps!" ); } if ( m_ImageMask != imageMask ) { m_ImageMask = imageMask; this->Modified(); m_MaskedImageStatisticsTimeStamp.Modified(); m_MaskedImageStatisticsCalculationTriggerBool = true; } } void PartialVolumeAnalysisHistogramCalculator::SetPlanarFigure( const mitk::PlanarFigure *planarFigure ) { if ( m_Image.IsNull() ) { itkExceptionMacro( << "Image needs to be set first!" ); } if ( m_PlanarFigure != planarFigure ) { m_PlanarFigure = planarFigure; this->Modified(); m_PlanarFigureStatisticsTimeStamp.Modified(); m_PlanarFigureStatisticsCalculationTriggerBool = true; } } void PartialVolumeAnalysisHistogramCalculator::SetMaskingMode( unsigned int mode ) { if ( m_MaskingMode != mode ) { m_MaskingMode = mode; m_MaskingModeChanged = true; this->Modified(); } } void PartialVolumeAnalysisHistogramCalculator::SetMaskingModeToNone() { if ( m_MaskingMode != MASKING_MODE_NONE ) { m_MaskingMode = MASKING_MODE_NONE; m_MaskingModeChanged = true; this->Modified(); } } void PartialVolumeAnalysisHistogramCalculator::SetMaskingModeToImage() { if ( m_MaskingMode != MASKING_MODE_IMAGE ) { m_MaskingMode = MASKING_MODE_IMAGE; m_MaskingModeChanged = true; this->Modified(); } } void PartialVolumeAnalysisHistogramCalculator::SetMaskingModeToPlanarFigure() { if ( m_MaskingMode != MASKING_MODE_PLANARFIGURE ) { m_MaskingMode = MASKING_MODE_PLANARFIGURE; m_MaskingModeChanged = true; this->Modified(); } } bool PartialVolumeAnalysisHistogramCalculator::ComputeStatistics() { MITK_DEBUG << "ComputeStatistics() start"; if ( m_Image.IsNull() ) { itkExceptionMacro( << "Image not set!" ); } if ( m_Image->GetReferenceCount() == 1 ) { MITK_DEBUG << "No Stats calculated; no one else holds a reference on it"; return false; } // 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_ImageStatisticsTimeStamp.GetMTime(); unsigned long maskedImageMTime = m_MaskedImageStatisticsTimeStamp.GetMTime(); unsigned long planarFigureMTime = m_PlanarFigureStatisticsTimeStamp.GetMTime(); bool imageStatisticsCalculationTrigger = m_ImageStatisticsCalculationTriggerBool; bool maskedImageStatisticsCalculationTrigger = m_MaskedImageStatisticsCalculationTriggerBool; bool planarFigureStatisticsCalculationTrigger = m_PlanarFigureStatisticsCalculationTriggerBool; if ( /*prevent calculation without mask*/!m_ForceUpdate &&( m_MaskingMode == MASKING_MODE_NONE || ( ((m_MaskingMode != MASKING_MODE_NONE) || (imageMTime > m_Image->GetMTime() && !imageStatisticsCalculationTrigger)) && ((m_MaskingMode != MASKING_MODE_IMAGE) || (m_ImageMask.IsNotNull() && maskedImageMTime > m_ImageMask->GetMTime() && !maskedImageStatisticsCalculationTrigger)) && ((m_MaskingMode != MASKING_MODE_PLANARFIGURE) || (m_PlanarFigure.IsNotNull() && planarFigureMTime > m_PlanarFigure->GetMTime() && !planarFigureStatisticsCalculationTrigger)) ) ) ) { MITK_DEBUG << "Returning, statistics already up to date!"; if ( m_MaskingModeChanged ) { m_MaskingModeChanged = false; return true; } else { return false; } } // Reset state changed flag m_MaskingModeChanged = false; // Depending on masking mode, extract and/or generate the required image // and mask data from the user input this->ExtractImageAndMask( ); Statistics *statistics; HistogramType::ConstPointer *histogram; switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: statistics = &m_ImageStatistics; histogram = &m_ImageHistogram; m_ImageStatisticsTimeStamp.Modified(); m_ImageStatisticsCalculationTriggerBool = false; break; case MASKING_MODE_IMAGE: statistics = &m_MaskedImageStatistics; histogram = &m_MaskedImageHistogram; m_MaskedImageStatisticsTimeStamp.Modified(); m_MaskedImageStatisticsCalculationTriggerBool = false; break; case MASKING_MODE_PLANARFIGURE: statistics = &m_PlanarFigureStatistics; histogram = &m_PlanarFigureHistogram; m_PlanarFigureStatisticsTimeStamp.Modified(); m_PlanarFigureStatisticsCalculationTriggerBool = false; break; } // Calculate statistics and histogram(s) if ( m_InternalImage->GetDimension() == 3 ) { if ( m_MaskingMode == MASKING_MODE_NONE ) { // Reset state changed flag AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateStatisticsUnmasked, 3, *statistics, histogram ); } else { AccessFixedDimensionByItk_3( m_InternalImage, InternalCalculateStatisticsMasked, 3, m_InternalImageMask3D.GetPointer(), *statistics, histogram ); } } else if ( m_InternalImage->GetDimension() == 2 ) { if ( m_MaskingMode == MASKING_MODE_NONE ) { AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateStatisticsUnmasked, 2, *statistics, histogram ); } else { AccessFixedDimensionByItk_3( m_InternalImage, InternalCalculateStatisticsMasked, 2, m_InternalImageMask2D.GetPointer(), *statistics, histogram ); } } else { MITK_ERROR << "ImageStatistics: Image dimension not supported!"; } // Release unused image smart pointers to free memory // m_InternalImage = mitk::Image::Pointer(); m_InternalImageMask3D = MaskImage3DType::Pointer(); m_InternalImageMask2D = MaskImage2DType::Pointer(); return true; } const PartialVolumeAnalysisHistogramCalculator::HistogramType * PartialVolumeAnalysisHistogramCalculator::GetHistogram( ) const { if ( m_Image.IsNull() ) { return NULL; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: return m_ImageHistogram; case MASKING_MODE_IMAGE: return m_MaskedImageHistogram; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureHistogram; } } const PartialVolumeAnalysisHistogramCalculator::Statistics & PartialVolumeAnalysisHistogramCalculator::GetStatistics( ) const { if ( m_Image.IsNull() ) { return m_EmptyStatistics; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: return m_ImageStatistics; case MASKING_MODE_IMAGE: return m_MaskedImageStatistics; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureStatistics; } } void PartialVolumeAnalysisHistogramCalculator::ExtractImageAndMask( ) { MITK_DEBUG << "ExtractImageAndMask( ) start"; if ( m_Image.IsNull() ) { throw std::runtime_error( "Error: image empty!" ); } mitk::Image *timeSliceImage = const_cast(m_Image.GetPointer());//imageTimeSelector->GetOutput(); switch ( m_MaskingMode ) { case MASKING_MODE_NONE: { m_InternalImage = timeSliceImage; int s = m_AdditionalResamplingImages.size(); m_InternalAdditionalResamplingImages.resize(s); for(int i=0; i(m_AdditionalResamplingImages[i].GetPointer()); } m_InternalImageMask2D = NULL; m_InternalImageMask3D = NULL; break; } case MASKING_MODE_IMAGE: { if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() > 1) ) { ImageTimeSelector::Pointer maskedImageTimeSelector = ImageTimeSelector::New(); maskedImageTimeSelector->SetInput( m_ImageMask ); maskedImageTimeSelector->SetTimeNr( 0 ); maskedImageTimeSelector->UpdateLargestPossibleRegion(); mitk::Image *timeSliceMaskedImage = maskedImageTimeSelector->GetOutput(); InternalMaskImage(timeSliceMaskedImage); if(m_UpsamplingFactor != 1) { InternalResampleImage(m_InternalImageMask3D); } AccessFixedDimensionByItk_1( timeSliceImage, InternalResampleImageFromMask, 3, -1 ); int s = m_AdditionalResamplingImages.size(); m_InternalAdditionalResamplingImages.resize(s); for(int i=0; iIsClosed() ) { throw std::runtime_error( "Masking not possible for non-closed figures" ); } const Geometry3D *imageGeometry = timeSliceImage->GetUpdatedGeometry(); if ( imageGeometry == NULL ) { throw std::runtime_error( "Image geometry invalid!" ); } const Geometry2D *planarFigureGeometry2D = m_PlanarFigure->GetGeometry2D(); if ( planarFigureGeometry2D == NULL ) { throw std::runtime_error( "Planar-Figure not yet initialized!" ); } const PlaneGeometry *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigureGeometry2D ); if ( planarFigureGeometry == NULL ) { throw std::runtime_error( "Non-planar planar figures not supported!" ); } // unsigned int axis = 2; // unsigned int slice = 0; AccessFixedDimensionByItk_3( timeSliceImage, InternalReorientImagePlane, 3, timeSliceImage->GetGeometry(), m_PlanarFigure->GetGeometry(), -1 ); AccessFixedDimensionByItk_1( m_InternalImage, InternalCalculateMaskFromPlanarFigure, 3, 2 ); int s = m_AdditionalResamplingImages.size(); for(int i=0; iGetGeometry(), m_PlanarFigure->GetGeometry(), i ); AccessFixedDimensionByItk_1( m_InternalAdditionalResamplingImages[i], InternalCropAdditionalImage, 3, i ); } } } } bool PartialVolumeAnalysisHistogramCalculator::GetPrincipalAxis( const Geometry3D *geometry, Vector3D vector, unsigned int &axis ) { vector.Normalize(); for ( unsigned int i = 0; i < 3; ++i ) { Vector3D axisVector = geometry->GetAxisVector( i ); axisVector.Normalize(); if ( fabs( fabs( axisVector * vector ) - 1.0) < mitk::eps ) { axis = i; return true; } } return false; } void PartialVolumeAnalysisHistogramCalculator::InternalMaskImage( mitk::Image *image ) { typedef itk::ImageMaskSpatialObject<3> ImageMaskSpatialObject; typedef itk::Image< unsigned char, 3 > ImageType; typedef itk::ImageRegion<3> RegionType; typedef mitk::ImageToItk CastType; CastType::Pointer caster = CastType::New(); caster->SetInput(image); caster->Update(); ImageMaskSpatialObject::Pointer maskSO = ImageMaskSpatialObject::New(); maskSO->SetImage ( caster->GetOutput() ); m_InternalMask3D = maskSO->GetAxisAlignedBoundingBoxRegion(); // check if bounding box is empty, if so set it to 1,1,1 // to prevent empty mask image if (m_InternalMask3D.GetSize()[0] == 0 ) { m_InternalMask3D.SetSize(0,1); m_InternalMask3D.SetSize(1,1); m_InternalMask3D.SetSize(2,1); } MITK_DEBUG << "Bounding Box Region: " << m_InternalMask3D; typedef itk::RegionOfInterestImageFilter< ImageType, MaskImage3DType > ROIFilterType; ROIFilterType::Pointer roi = ROIFilterType::New(); roi->SetRegionOfInterest(m_InternalMask3D); roi->SetInput(caster->GetOutput()); roi->Update(); m_InternalImageMask3D = roi->GetOutput(); MITK_DEBUG << "Created m_InternalImageMask3D"; } template < typename TPixel, unsigned int VImageDimension > void PartialVolumeAnalysisHistogramCalculator::InternalReorientImagePlane( const itk::Image< TPixel, VImageDimension > *image, mitk::Geometry3D* , mitk::Geometry3D* planegeo3D, int additionalIndex ) { MITK_DEBUG << "InternalReorientImagePlane() start"; typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< float, VImageDimension > FloatImageType; typedef itk::ResampleImageFilter ResamplerType; typename ResamplerType::Pointer resampler = ResamplerType::New(); mitk::PlaneGeometry* planegeo = dynamic_cast(planegeo3D); float upsamp = m_UpsamplingFactor; float gausssigma = m_GaussianSigma; // Spacing typename ResamplerType::SpacingType spacing = planegeo->GetSpacing(); spacing[0] = image->GetSpacing()[0] / upsamp; spacing[1] = image->GetSpacing()[1] / upsamp; spacing[2] = image->GetSpacing()[2]; if(m_PlanarFigureThickness) { spacing[2] = image->GetSpacing()[2] / upsamp; } resampler->SetOutputSpacing( spacing ); // Size typename ResamplerType::SizeType size; size[0] = planegeo->GetParametricExtentInMM(0) / spacing[0]; size[1] = planegeo->GetParametricExtentInMM(1) / spacing[1]; size[2] = 1+2*m_PlanarFigureThickness; // klaus add +2*m_PlanarFigureThickness MITK_DEBUG << "setting size2:="<SetSize( size ); // Origin typename mitk::Point3D orig = planegeo->GetOrigin(); typename mitk::Point3D corrorig; planegeo3D->WorldToIndex(orig,corrorig); corrorig[0] += 0.5/upsamp; corrorig[1] += 0.5/upsamp; if(m_PlanarFigureThickness) { float thickyyy = m_PlanarFigureThickness; thickyyy/=upsamp; corrorig[2] -= thickyyy; // klaus add -= (float)m_PlanarFigureThickness/upsamp statt += 0 } planegeo3D->IndexToWorld(corrorig,corrorig); resampler->SetOutputOrigin(corrorig ); // Direction typename ResamplerType::DirectionType direction; typename mitk::AffineTransform3D::MatrixType matrix = planegeo->GetIndexToWorldTransform()->GetMatrix(); for(int c=0; cSetOutputDirection( direction ); // Gaussian interpolation if(gausssigma != 0) { double sigma[3]; for( unsigned int d = 0; d < 3; d++ ) { sigma[d] = gausssigma * image->GetSpacing()[d]; } double alpha = 2.0; typedef itk::GaussianInterpolateImageFunction GaussianInterpolatorType; typename GaussianInterpolatorType::Pointer interpolator = GaussianInterpolatorType::New(); interpolator->SetInputImage( image ); interpolator->SetParameters( sigma, alpha ); resampler->SetInterpolator( interpolator ); } else { // typedef typename itk::BSplineInterpolateImageFunction // InterpolatorType; typedef typename itk::LinearInterpolateImageFunction InterpolatorType; typename InterpolatorType::Pointer interpolator = InterpolatorType::New(); interpolator->SetInputImage( image ); resampler->SetInterpolator( interpolator ); } // Other resampling options resampler->SetInput( image ); resampler->SetDefaultPixelValue(0); MITK_DEBUG << "Resampling requested image plane ... "; resampler->Update(); MITK_DEBUG << " ... done"; if(additionalIndex < 0) { this->m_InternalImage = mitk::Image::New(); this->m_InternalImage->InitializeByItk( resampler->GetOutput() ); this->m_InternalImage->SetVolume( resampler->GetOutput()->GetBufferPointer() ); } else { unsigned int myIndex = additionalIndex; this->m_InternalAdditionalResamplingImages.push_back(mitk::Image::New()); this->m_InternalAdditionalResamplingImages[myIndex]->InitializeByItk( resampler->GetOutput() ); this->m_InternalAdditionalResamplingImages[myIndex]->SetVolume( resampler->GetOutput()->GetBufferPointer() ); } } template < typename TPixel, unsigned int VImageDimension > void PartialVolumeAnalysisHistogramCalculator::InternalResampleImageFromMask( const itk::Image< TPixel, VImageDimension > *image, int additionalIndex ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typename ImageType::Pointer outImage = ImageType::New(); outImage->SetSpacing( m_InternalImageMask3D->GetSpacing() ); // Set the image spacing outImage->SetOrigin( m_InternalImageMask3D->GetOrigin() ); // Set the image origin outImage->SetDirection( m_InternalImageMask3D->GetDirection() ); // Set the image direction outImage->SetRegions( m_InternalImageMask3D->GetLargestPossibleRegion() ); outImage->Allocate(); outImage->FillBuffer(0); typedef itk::InterpolateImageFunction BaseInterpType; typedef itk::GaussianInterpolateImageFunction GaussianInterpolatorType; typedef itk::LinearInterpolateImageFunction LinearInterpolatorType; typename BaseInterpType::Pointer interpolator; if(m_GaussianSigma != 0) { double sigma[3]; for( unsigned int d = 0; d < 3; d++ ) { sigma[d] = m_GaussianSigma * image->GetSpacing()[d]; } double alpha = 2.0; interpolator = GaussianInterpolatorType::New(); dynamic_cast(interpolator.GetPointer())->SetParameters( sigma, alpha ); } else { interpolator = LinearInterpolatorType::New(); } interpolator->SetInputImage( image ); itk::ImageRegionConstIterator itmask(m_InternalImageMask3D, m_InternalImageMask3D->GetLargestPossibleRegion()); itk::ImageRegionIterator itimage(outImage, outImage->GetLargestPossibleRegion()); itmask.GoToBegin(); itimage.GoToBegin(); itk::Point< double, 3 > point; itk::ContinuousIndex< double, 3 > index; while( !itmask.IsAtEnd() ) { if(itmask.Get() != 0) { outImage->TransformIndexToPhysicalPoint (itimage.GetIndex(), point); image->TransformPhysicalPointToContinuousIndex(point, index); itimage.Set(interpolator->EvaluateAtContinuousIndex(index)); } ++itmask; ++itimage; } if(additionalIndex < 0) { this->m_InternalImage = mitk::Image::New(); this->m_InternalImage->InitializeByItk( outImage.GetPointer() ); this->m_InternalImage->SetVolume( outImage->GetBufferPointer() ); } else { this->m_InternalAdditionalResamplingImages[additionalIndex] = mitk::Image::New(); this->m_InternalAdditionalResamplingImages[additionalIndex]->InitializeByItk( outImage.GetPointer() ); this->m_InternalAdditionalResamplingImages[additionalIndex]->SetVolume( outImage->GetBufferPointer() ); } } void PartialVolumeAnalysisHistogramCalculator::InternalResampleImage( const MaskImage3DType *image ) { typedef itk::ResampleImageFilter ResamplerType; ResamplerType::Pointer resampler = ResamplerType::New(); // Size ResamplerType::SizeType size; size[0] = m_UpsamplingFactor * image->GetLargestPossibleRegion().GetSize()[0]; size[1] = m_UpsamplingFactor * image->GetLargestPossibleRegion().GetSize()[1]; size[2] = m_UpsamplingFactor * image->GetLargestPossibleRegion().GetSize()[2];; resampler->SetSize( size ); // Origin mitk::Point3D orig = image->GetOrigin(); resampler->SetOutputOrigin(orig ); // Spacing ResamplerType::SpacingType spacing; spacing[0] = image->GetSpacing()[0] / m_UpsamplingFactor; spacing[1] = image->GetSpacing()[1] / m_UpsamplingFactor; spacing[2] = image->GetSpacing()[2] / m_UpsamplingFactor; resampler->SetOutputSpacing( spacing ); resampler->SetOutputDirection( image->GetDirection() ); typedef itk::NearestNeighborInterpolateImageFunction InterpolatorType; InterpolatorType::Pointer interpolator = InterpolatorType::New(); interpolator->SetInputImage( image ); resampler->SetInterpolator( interpolator ); // Other resampling options resampler->SetInput( image ); resampler->SetDefaultPixelValue(0); resampler->Update(); m_InternalImageMask3D = resampler->GetOutput(); } template < typename TPixel, unsigned int VImageDimension > void PartialVolumeAnalysisHistogramCalculator::InternalCalculateStatisticsUnmasked( const itk::Image< TPixel, VImageDimension > *image, Statistics &statistics, typename HistogramType::ConstPointer *histogram ) { MITK_DEBUG << "InternalCalculateStatisticsUnmasked()"; typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned char, VImageDimension > MaskImageType; typedef typename ImageType::IndexType IndexType; // Progress listening... typedef itk::SimpleMemberCommand< PartialVolumeAnalysisHistogramCalculator > ITKCommandType; ITKCommandType::Pointer progressListener; progressListener = ITKCommandType::New(); progressListener->SetCallbackFunction( this, &PartialVolumeAnalysisHistogramCalculator::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() ); statistics.N = image->GetBufferedRegion().GetNumberOfPixels(); statistics.Min = statisticsFilter->GetMinimum(); statistics.Max = statisticsFilter->GetMaximum(); statistics.Mean = statisticsFilter->GetMean(); statistics.Median = 0.0; statistics.Sigma = statisticsFilter->GetSigma(); statistics.RMS = sqrt( statistics.Mean * statistics.Mean + statistics.Sigma * statistics.Sigma ); typename ImageType::Pointer inImage = const_cast(image); // Calculate histogram typedef itk::Statistics::ScalarImageToHistogramGenerator< ImageType > HistogramGeneratorType; typename HistogramGeneratorType::Pointer histogramGenerator = HistogramGeneratorType::New(); histogramGenerator->SetInput( inImage ); histogramGenerator->SetMarginalScale( 10 ); // Defines y-margin width of histogram histogramGenerator->SetNumberOfBins( m_NumberOfBins ); // CT range [-1024, +2048] --> bin size 4 values histogramGenerator->SetHistogramMin( statistics.Min ); histogramGenerator->SetHistogramMax( statistics.Max ); histogramGenerator->Compute(); *histogram = histogramGenerator->GetOutput(); } template < typename TPixel, unsigned int VImageDimension > void PartialVolumeAnalysisHistogramCalculator::InternalCalculateStatisticsMasked( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned char, VImageDimension > *maskImage, Statistics &, typename HistogramType::ConstPointer *histogram ) { MITK_DEBUG << "InternalCalculateStatisticsMasked() start"; typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned char, VImageDimension > MaskImageType; typedef typename ImageType::IndexType IndexType; // generate a list sample of angles at positions // where the fiber-prob is higher than .2*maxprob typedef TPixel MeasurementType; const unsigned int MeasurementVectorLength = 1; typedef itk::Vector< MeasurementType , MeasurementVectorLength > MeasurementVectorType; typedef itk::Statistics::ListSample< MeasurementVectorType > ListSampleType; typename ListSampleType::Pointer listSample = ListSampleType::New(); listSample->SetMeasurementVectorSize( MeasurementVectorLength ); itk::ImageRegionConstIterator itmask(maskImage, maskImage->GetLargestPossibleRegion()); itk::ImageRegionConstIterator itimage(image, image->GetLargestPossibleRegion()); itmask.GoToBegin(); itimage.GoToBegin(); while( !itmask.IsAtEnd() ) { if(itmask.Get() != 0) { // apend to list MeasurementVectorType mv; mv[0] = ( MeasurementType ) itimage.Get(); listSample->PushBack(mv); } ++itmask; ++itimage; } // generate a histogram from the list sample typedef double HistogramMeasurementType; typedef itk::Statistics::Histogram< HistogramMeasurementType, itk::Statistics::DenseFrequencyContainer2 > HistogramType; typedef itk::Statistics::SampleToHistogramFilter< ListSampleType, HistogramType > GeneratorType; typename GeneratorType::Pointer generator = GeneratorType::New(); typename GeneratorType::HistogramType::SizeType size(MeasurementVectorLength); size.Fill(m_NumberOfBins); generator->SetHistogramSize( size ); generator->SetInput( listSample ); generator->SetMarginalScale( 10.0 ); generator->Update(); *histogram = generator->GetOutput(); } template < typename TPixel, unsigned int VImageDimension > void PartialVolumeAnalysisHistogramCalculator::InternalCropAdditionalImage( itk::Image< TPixel, VImageDimension > *image, int additionalIndex ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::RegionOfInterestImageFilter< ImageType, ImageType > ROIFilterType; typename ROIFilterType::Pointer roi = ROIFilterType::New(); roi->SetRegionOfInterest(m_CropRegion); roi->SetInput(image); roi->Update(); m_InternalAdditionalResamplingImages[additionalIndex] = mitk::Image::New(); m_InternalAdditionalResamplingImages[additionalIndex]->InitializeByItk(roi->GetOutput()); m_InternalAdditionalResamplingImages[additionalIndex]->SetVolume(roi->GetOutput()->GetBufferPointer()); } template < typename TPixel, unsigned int VImageDimension > void PartialVolumeAnalysisHistogramCalculator::InternalCalculateMaskFromPlanarFigure( itk::Image< TPixel, VImageDimension > *image, unsigned int axis ) { MITK_DEBUG << "InternalCalculateMaskFromPlanarFigure() start"; typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::CastImageFilter< ImageType, MaskImage3DType > CastFilterType; // Generate mask image as new image with same header as input image and // initialize with "1". MaskImage3DType::Pointer newMaskImage = MaskImage3DType::New(); newMaskImage->SetSpacing( image->GetSpacing() ); // Set the image spacing newMaskImage->SetOrigin( image->GetOrigin() ); // Set the image origin newMaskImage->SetDirection( image->GetDirection() ); // Set the image direction newMaskImage->SetRegions( image->GetLargestPossibleRegion() ); newMaskImage->Allocate(); newMaskImage->FillBuffer( 1 ); // Generate VTK polygon from (closed) PlanarFigure polyline // (The polyline points are shifted by -0.5 in z-direction to make sure // that the extrusion filter, which afterwards elevates all points by +0.5 // in z-direction, creates a 3D object which is cut by the the plane z=0) const mitk::Geometry2D *planarFigureGeometry2D = m_PlanarFigure->GetGeometry2D(); const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); const mitk::Geometry3D *imageGeometry3D = m_InternalImage->GetGeometry( 0 ); vtkPolyData *polyline = vtkPolyData::New(); polyline->Allocate( 1, 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; } // Create VTK polydata object of polyline contour bool outOfBounds = false; vtkPoints *points = vtkPoints::New(); typename PlanarFigure::PolyLineType::const_iterator it; for ( it = planarFigurePolyline.begin(); it != planarFigurePolyline.end(); ++it ) { Point3D point3D; // Convert 2D point back to the local index coordinates of the selected // image - mitk::Point2D point2D = it->Point; + mitk::Point2D point2D = *it; planarFigureGeometry2D->WorldToIndex(point2D, point2D); point2D[0] -= 0.5/m_UpsamplingFactor; point2D[1] -= 0.5/m_UpsamplingFactor; planarFigureGeometry2D->IndexToWorld(point2D, point2D); planarFigureGeometry2D->Map( point2D, 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 ); point3D[i0] += 0.5; point3D[i1] += 0.5; // Add point to polyline array points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 ); } polyline->SetPoints( points ); points->Delete(); if ( outOfBounds ) { polyline->Delete(); throw std::runtime_error( "Figure at least partially outside of image bounds!" ); } std::size_t numberOfPoints = planarFigurePolyline.size(); vtkIdType *ptIds = new vtkIdType[numberOfPoints]; for ( std::size_t i = 0; i < numberOfPoints; ++i ) { ptIds[i] = i; } polyline->InsertNextCell( VTK_POLY_LINE, numberOfPoints, ptIds ); // Extrude the generated contour polygon vtkLinearExtrusionFilter *extrudeFilter = vtkLinearExtrusionFilter::New(); extrudeFilter->SetInputData( polyline ); extrudeFilter->SetScaleFactor( 1 ); extrudeFilter->SetExtrusionTypeToNormalExtrusion(); extrudeFilter->SetVector( 0.0, 0.0, 1.0 ); // Make a stencil from the extruded polygon vtkPolyDataToImageStencil *polyDataToImageStencil = vtkPolyDataToImageStencil::New(); polyDataToImageStencil->SetInputConnection( extrudeFilter->GetOutputPort() ); // Export from ITK to VTK (to use a VTK filter) typedef itk::VTKImageImport< MaskImage3DType > ImageImportType; typedef itk::VTKImageExport< MaskImage3DType > ImageExportType; typename ImageExportType::Pointer itkExporter = ImageExportType::New(); itkExporter->SetInput( newMaskImage ); vtkImageImport *vtkImporter = vtkImageImport::New(); this->ConnectPipelines( itkExporter, vtkImporter ); vtkImporter->Update(); // Apply the generated image stencil to the input image vtkImageStencil *imageStencilFilter = vtkImageStencil::New(); imageStencilFilter->SetInputData( vtkImporter->GetOutput() ); imageStencilFilter->SetStencilConnection(polyDataToImageStencil->GetOutputPort() ); imageStencilFilter->ReverseStencilOff(); imageStencilFilter->SetBackgroundValue( 0 ); imageStencilFilter->Update(); // Export from VTK back to ITK vtkImageExport *vtkExporter = vtkImageExport::New(); vtkExporter->SetInputData( imageStencilFilter->GetOutput() ); vtkExporter->Update(); typename ImageImportType::Pointer itkImporter = ImageImportType::New(); this->ConnectPipelines( vtkExporter, itkImporter ); itkImporter->Update(); // calculate cropping bounding box m_InternalImageMask3D = itkImporter->GetOutput(); m_InternalImageMask3D->SetDirection(image->GetDirection()); itk::ImageRegionIterator itmask(m_InternalImageMask3D, m_InternalImageMask3D->GetLargestPossibleRegion()); itmask.GoToBegin(); while( !itmask.IsAtEnd() ) { if(itmask.Get() != 0) { typename ImageType::IndexType index = itmask.GetIndex(); for(unsigned int thick=0; thick<2*m_PlanarFigureThickness+1; thick++) { index[axis] = thick; m_InternalImageMask3D->SetPixel(index, itmask.Get()); } } ++itmask; } itmask.GoToBegin(); itk::ImageRegionIterator itimage(image, image->GetLargestPossibleRegion()); itimage.GoToBegin(); typename ImageType::SizeType lowersize; lowersize.Fill(std::numeric_limits::max()); typename ImageType::SizeType uppersize; uppersize.Fill(std::numeric_limits::min()); while( !itmask.IsAtEnd() ) { if(itmask.Get() == 0) { itimage.Set(0); } else { typename ImageType::IndexType index = itimage.GetIndex(); typename ImageType::SizeType signedindex; signedindex[0] = index[0]; signedindex[1] = index[1]; signedindex[2] = index[2]; lowersize[0] = signedindex[0] < lowersize[0] ? signedindex[0] : lowersize[0]; lowersize[1] = signedindex[1] < lowersize[1] ? signedindex[1] : lowersize[1]; lowersize[2] = signedindex[2] < lowersize[2] ? signedindex[2] : lowersize[2]; uppersize[0] = signedindex[0] > uppersize[0] ? signedindex[0] : uppersize[0]; uppersize[1] = signedindex[1] > uppersize[1] ? signedindex[1] : uppersize[1]; uppersize[2] = signedindex[2] > uppersize[2] ? signedindex[2] : uppersize[2]; } ++itmask; ++itimage; } typename ImageType::IndexType index; index[0] = lowersize[0]; index[1] = lowersize[1]; index[2] = lowersize[2]; typename ImageType::SizeType size; size[0] = uppersize[0] - lowersize[0] + 1; size[1] = uppersize[1] - lowersize[1] + 1; size[2] = uppersize[2] - lowersize[2] + 1; m_CropRegion = itk::ImageRegion<3>(index, size); // crop internal image typedef itk::RegionOfInterestImageFilter< ImageType, ImageType > ROIFilterType; typename ROIFilterType::Pointer roi = ROIFilterType::New(); roi->SetRegionOfInterest(m_CropRegion); roi->SetInput(image); roi->Update(); m_InternalImage = mitk::Image::New(); m_InternalImage->InitializeByItk(roi->GetOutput()); m_InternalImage->SetVolume(roi->GetOutput()->GetBufferPointer()); // crop internal mask typedef itk::RegionOfInterestImageFilter< MaskImage3DType, MaskImage3DType > ROIMaskFilterType; typename ROIMaskFilterType::Pointer roi2 = ROIMaskFilterType::New(); roi2->SetRegionOfInterest(m_CropRegion); roi2->SetInput(m_InternalImageMask3D); roi2->Update(); m_InternalImageMask3D = roi2->GetOutput(); // Clean up VTK objects polyline->Delete(); extrudeFilter->Delete(); polyDataToImageStencil->Delete(); vtkImporter->Delete(); imageStencilFilter->Delete(); //vtkExporter->Delete(); // TODO: crashes when outcommented; memory leak?? delete[] ptIds; } void PartialVolumeAnalysisHistogramCalculator::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 PartialVolumeAnalysisHistogramCalculator::MaskedStatisticsProgressUpdate() { this->InvokeEvent( itk::ProgressEvent() ); } } diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index 9a4869d897..074abe1407 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,1272 +1,1272 @@ /*=================================================================== 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 #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_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::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; } 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_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; return true; } 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 ) { CastToItkImage( timeSliceImage, m_InternalImageMask3D ); m_InternalImageMask3D->FillBuffer(1); } if( m_InternalImage->GetDimension() == 2 ) { CastToItkImage( timeSliceImage, 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 Geometry3D *imageGeometry = timeSliceImage->GetGeometry(); if ( imageGeometry == NULL ) { throw std::runtime_error( "Image geometry invalid!" ); } const Geometry2D *planarFigureGeometry2D = m_PlanarFigure->GetGeometry2D(); if ( planarFigureGeometry2D == NULL ) { throw std::runtime_error( "Planar-Figure not yet initialized!" ); } const PlaneGeometry *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigureGeometry2D ); 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 Geometry3D *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; } 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 itk::Image< unsigned short, VImageDimension > MaskImageType; 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.Label = 1; statistics.N = image->GetBufferedRegion().GetNumberOfPixels(); statistics.Min = statisticsFilter->GetMinimum(); statistics.Max = statisticsFilter->GetMaximum(); statistics.Mean = statisticsFilter->GetMean(); statistics.Median = 0.0; statistics.Sigma = statisticsFilter->GetSigma(); statistics.RMS = sqrt( statistics.Mean * statistics.Mean + statistics.Sigma * statistics.Sigma ); statistics.MinIndex.set_size(image->GetImageDimension()); statistics.MaxIndex.set_size(image->GetImageDimension()); for (unsigned int i=0; iGetIndexOfMaximum()[i]; statistics.MinIndex[i] = minMaxFilter->GetIndexOfMinimum()[i]; } statisticsContainer->push_back( statistics ); // Calculate histogram typename HistogramGeneratorType::Pointer histogramGenerator = HistogramGeneratorType::New(); histogramGenerator->SetInput( image ); histogramGenerator->SetMarginalScale( 100 ); histogramGenerator->SetNumberOfBins( 768 ); histogramGenerator->SetHistogramMin( statistics.Min ); histogramGenerator->SetHistogramMax( statistics.Max ); 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]; } // 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 = ( m_DoIgnorePixelValue && (m_MaskingMode == MASKING_MODE_NONE) ) ? 768 : 384; 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 ) { histogramContainer->push_back( HistogramType::ConstPointer( labelStatisticsFilter->GetHistogram( (*it) ) ) ); Statistics statistics; statistics.Label = (*it); statistics.N = labelStatisticsFilter->GetCount( *it ); statistics.Min = labelStatisticsFilter->GetMinimum( *it ); statistics.Max = labelStatisticsFilter->GetMaximum( *it ); statistics.Mean = labelStatisticsFilter->GetMean( *it ); statistics.Median = labelStatisticsFilter->GetMedian( *it ); statistics.Sigma = labelStatisticsFilter->GetSigma( *it ); statistics.RMS = sqrt( statistics.Mean * statistics.Mean + statistics.Sigma * statistics.Sigma ); // restrict image to mask area for min/max index calculation typedef itk::MaskImageFilter< ImageType, MaskImageType, ImageType > MaskImageFilterType; typename MaskImageFilterType::Pointer masker = MaskImageFilterType::New(); masker->SetOutsideValue( (statistics.Min+statistics.Max)/2 ); 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() ); statistics.MinIndex.set_size(adaptedImage->GetImageDimension()); statistics.MaxIndex.set_size(adaptedImage->GetImageDimension()); typename MinMaxFilterType::IndexType tempMaxIndex = minMaxFilter->GetIndexOfMaximum(); typename MinMaxFilterType::IndexType tempMinIndex = minMaxFilter->GetIndexOfMinimum(); // FIX BUG 14644 //If a PlanarFigure is used for segmentation the //adaptedImage is a single slice (2D). Adding the // 3. dimension. if (m_MaskingMode == MASKING_MODE_PLANARFIGURE && m_Image->GetDimension()==3) { statistics.MaxIndex.set_size(m_Image->GetDimension()); statistics.MaxIndex[m_PlanarFigureCoordinate0]=tempMaxIndex[0]; statistics.MaxIndex[m_PlanarFigureCoordinate1]=tempMaxIndex[1]; statistics.MaxIndex[m_PlanarFigureAxis]=m_PlanarFigureSlice; statistics.MinIndex.set_size(m_Image->GetDimension()); statistics.MinIndex[m_PlanarFigureCoordinate0]=tempMinIndex[0]; statistics.MinIndex[m_PlanarFigureCoordinate1]=tempMinIndex[1]; statistics.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 < 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::Geometry2D *planarFigureGeometry2D = m_PlanarFigure->GetGeometry2D(); const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 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); const mitk::Geometry3D *imageGeometry3D = m_Image->GetGeometry( 0 ); // 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 - planarFigureGeometry2D->Map( it->Point, point3D ); + planarFigureGeometry2D->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) { - planarFigureGeometry2D->Map(it->Point, point3D); + planarFigureGeometry2D->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/mitkIntensityProfile.cpp b/Modules/ImageStatistics/mitkIntensityProfile.cpp index f899c6cd26..b5fa201dda 100644 --- a/Modules/ImageStatistics/mitkIntensityProfile.cpp +++ b/Modules/ImageStatistics/mitkIntensityProfile.cpp @@ -1,337 +1,337 @@ /*=================================================================== 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 #include #include #include #include #include #include #include "mitkIntensityProfile.h" using namespace mitk; template static void ReadPixel(const PixelType&, Image::Pointer image, const Index3D& index, ScalarType* returnValue) { switch (image->GetDimension()) { case 2: { ImagePixelReadAccessor readAccess(image, image->GetSliceData(0)); *returnValue = readAccess.GetPixelByIndex(reinterpret_cast&>(index)); break; } case 3: { ImagePixelReadAccessor readAccess(image, image->GetVolumeData(0)); *returnValue = readAccess.GetPixelByIndex(index); break; } default: *returnValue = 0; break; } } static IntensityProfile::Pointer ComputeIntensityProfile(Image::Pointer image, itk::PolyLineParametricPath<3>::Pointer path) { IntensityProfile::Pointer intensityProfile = IntensityProfile::New(); itk::PolyLineParametricPath<3>::InputType input = path->StartOfInput(); Geometry3D* imageGeometry = image->GetGeometry(); const PixelType pixelType = image->GetPixelType(); IntensityProfile::MeasurementVectorType measurementVector; itk::PolyLineParametricPath<3>::OffsetType offset; Point3D worldPoint; Index3D index; do { imageGeometry->IndexToWorld(path->Evaluate(input), worldPoint); imageGeometry->WorldToIndex(worldPoint, index); mitkPixelTypeMultiplex3(ReadPixel, pixelType, image, index, measurementVector.GetDataPointer()); intensityProfile->PushBack(measurementVector); offset = path->IncrementInput(input); } while ((offset[0] | offset[1] | offset[2]) != 0); return intensityProfile; } template static typename itk::InterpolateImageFunction::Pointer CreateInterpolateImageFunction(InterpolateImageFunction::Enum interpolator) { switch (interpolator) { case InterpolateImageFunction::NearestNeighbor: return itk::NearestNeighborInterpolateImageFunction::New().GetPointer(); case InterpolateImageFunction::Linear: return itk::LinearInterpolateImageFunction::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Blackman_3: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Blackman_4: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Blackman_5: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Cosine_3: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Cosine_4: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Cosine_5: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Hamming_3: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Hamming_4: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Hamming_5: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Lanczos_3: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Lanczos_4: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Lanczos_5: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Welch_3: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Welch_4: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Welch_5: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); default: return itk::NearestNeighborInterpolateImageFunction::New().GetPointer(); } } template static void ComputeIntensityProfile(itk::Image* image, itk::PolyLineParametricPath<3>::Pointer path, unsigned int numSamples, InterpolateImageFunction::Enum interpolator, IntensityProfile::Pointer intensityProfile) { typename itk::InterpolateImageFunction >::Pointer interpolateImageFunction = CreateInterpolateImageFunction >(interpolator); interpolateImageFunction->SetInputImage(image); const itk::PolyLineParametricPath<3>::InputType startOfInput = path->StartOfInput(); const itk::PolyLineParametricPath<3>::InputType delta = 1.0 / (numSamples - 1); IntensityProfile::MeasurementVectorType measurementVector; for (unsigned int i = 0; i < numSamples; ++i) { measurementVector[0] = interpolateImageFunction->EvaluateAtContinuousIndex(path->Evaluate(startOfInput + i * delta)); intensityProfile->PushBack(measurementVector); } } static IntensityProfile::Pointer ComputeIntensityProfile(Image::Pointer image, itk::PolyLineParametricPath<3>::Pointer path, unsigned int numSamples, InterpolateImageFunction::Enum interpolator) { IntensityProfile::Pointer intensityProfile = IntensityProfile::New(); AccessFixedDimensionByItk_n(image, ComputeIntensityProfile, 3, (path, numSamples, interpolator, intensityProfile)); return intensityProfile; } class AddPolyLineElementToPath { public: AddPolyLineElementToPath(const Geometry2D* planarFigureGeometry, const Geometry3D* imageGeometry, itk::PolyLineParametricPath<3>::Pointer path) : m_PlanarFigureGeometry(planarFigureGeometry), m_ImageGeometry(imageGeometry), m_Path(path) { } void operator()(const PlanarFigure::PolyLineElement& polyLineElement) { - m_PlanarFigureGeometry->Map(polyLineElement.Point, m_WorldPoint); + m_PlanarFigureGeometry->Map(polyLineElement, m_WorldPoint); m_ImageGeometry->WorldToIndex(m_WorldPoint, m_ContinuousIndexPoint); m_Vertex.CastFrom(m_ContinuousIndexPoint); m_Path->AddVertex(m_Vertex); } private: const Geometry2D* m_PlanarFigureGeometry; const Geometry3D* m_ImageGeometry; itk::PolyLineParametricPath<3>::Pointer m_Path; Point3D m_WorldPoint; Point3D m_ContinuousIndexPoint; itk::PolyLineParametricPath<3>::ContinuousIndexType m_Vertex; }; static itk::PolyLineParametricPath<3>::Pointer CreatePathFromPlanarFigure(Geometry3D* imageGeometry, PlanarFigure* planarFigure) { itk::PolyLineParametricPath<3>::Pointer path = itk::PolyLineParametricPath<3>::New(); const PlanarFigure::PolyLineType polyLine = planarFigure->GetPolyLine(0); std::for_each(polyLine.begin(), polyLine.end(), AddPolyLineElementToPath(planarFigure->GetGeometry2D(), imageGeometry, path)); return path; } static void AddPointToPath(const Geometry3D* imageGeometry, const Point3D& point, itk::PolyLineParametricPath<3>::Pointer path) { Point3D continuousIndexPoint; imageGeometry->WorldToIndex(point, continuousIndexPoint); itk::PolyLineParametricPath<3>::ContinuousIndexType vertex; vertex.CastFrom(continuousIndexPoint); path->AddVertex(vertex); } static itk::PolyLineParametricPath<3>::Pointer CreatePathFromPoints(Geometry3D* imageGeometry, const Point3D& startPoint, const Point3D& endPoint) { itk::PolyLineParametricPath<3>::Pointer path = itk::PolyLineParametricPath<3>::New(); AddPointToPath(imageGeometry, startPoint, path); AddPointToPath(imageGeometry, endPoint, path); return path; } IntensityProfile::Pointer mitk::ComputeIntensityProfile(Image::Pointer image, PlanarFigure::Pointer planarFigure) { return ::ComputeIntensityProfile(image, CreatePathFromPlanarFigure(image->GetGeometry(), planarFigure)); } IntensityProfile::Pointer mitk::ComputeIntensityProfile(Image::Pointer image, PlanarLine::Pointer planarLine, unsigned int numSamples, InterpolateImageFunction::Enum interpolator) { return ::ComputeIntensityProfile(image, CreatePathFromPlanarFigure(image->GetGeometry(), planarLine.GetPointer()), numSamples, interpolator); } IntensityProfile::Pointer mitk::ComputeIntensityProfile(Image::Pointer image, const Point3D& startPoint, const Point3D& endPoint, unsigned int numSamples, InterpolateImageFunction::Enum interpolator) { return ::ComputeIntensityProfile(image, CreatePathFromPoints(image->GetGeometry(), startPoint, endPoint), numSamples, interpolator); } IntensityProfile::InstanceIdentifier mitk::ComputeGlobalMaximum(IntensityProfile::Pointer intensityProfile) { IntensityProfile::MeasurementType max = -vcl_numeric_limits::max(); IntensityProfile::InstanceIdentifier maxIndex = 0; IntensityProfile::ConstIterator end = intensityProfile->End(); IntensityProfile::MeasurementType measurement; for (IntensityProfile::ConstIterator it = intensityProfile->Begin(); it != end; ++it) { measurement = it.GetMeasurementVector()[0]; if (measurement > max) { max = measurement; maxIndex = it.GetInstanceIdentifier(); } } return maxIndex; } IntensityProfile::InstanceIdentifier mitk::ComputeGlobalMinimum(IntensityProfile::Pointer intensityProfile) { IntensityProfile::MeasurementType min = vcl_numeric_limits::max(); IntensityProfile::InstanceIdentifier minIndex = 0; IntensityProfile::ConstIterator end = intensityProfile->End(); IntensityProfile::MeasurementType measurement; for (IntensityProfile::ConstIterator it = intensityProfile->Begin(); it != end; ++it) { measurement = it.GetMeasurementVector()[0]; if (measurement < min) { min = measurement; minIndex = it.GetInstanceIdentifier(); } } return minIndex; } IntensityProfile::InstanceIdentifier mitk::ComputeCenterOfMaximumArea(IntensityProfile::Pointer intensityProfile, IntensityProfile::InstanceIdentifier radius) { const IntensityProfile::MeasurementType min = intensityProfile->GetMeasurementVector(ComputeGlobalMinimum(intensityProfile))[0]; const IntensityProfile::InstanceIdentifier areaWidth = 1 + 2 * radius; IntensityProfile::MeasurementType maxArea = 0; for (IntensityProfile::InstanceIdentifier i = 0; i < areaWidth; ++i) maxArea += intensityProfile->GetMeasurementVector(i)[0] - min; const IntensityProfile::InstanceIdentifier lastIndex = intensityProfile->Size() - areaWidth; IntensityProfile::InstanceIdentifier centerOfMaxArea = radius; IntensityProfile::MeasurementType area = maxArea; for (IntensityProfile::InstanceIdentifier i = 1; i <= lastIndex; ++i) { area += intensityProfile->GetMeasurementVector(i + areaWidth - 1)[0] - min; area -= intensityProfile->GetMeasurementVector(i - 1)[0] - min; if (area > maxArea) { maxArea = area; centerOfMaxArea = i + radius; // TODO: If multiple areas in the neighborhood have the same intensity chose the middle one instead of the first one. } } return centerOfMaxArea; } std::vector mitk::CreateVectorFromIntensityProfile(IntensityProfile::Pointer intensityProfile) { std::vector result; result.reserve(intensityProfile->Size()); IntensityProfile::ConstIterator end = intensityProfile->End(); for (IntensityProfile::ConstIterator it = intensityProfile->Begin(); it != end; ++it) result.push_back(it.GetMeasurementVector()[0]); return result; } IntensityProfile::Pointer mitk::CreateIntensityProfileFromVector(const std::vector& vector) { const IntensityProfile::InstanceIdentifier size = vector.size(); IntensityProfile::Pointer result = IntensityProfile::New(); result->Resize(size); for (IntensityProfile::InstanceIdentifier i = 0; i < size; ++i) result->SetMeasurement(i, 0, vector[i]); return result; } diff --git a/Modules/PlanarFigure/CMakeLists.txt b/Modules/PlanarFigure/CMakeLists.txt index cf6f567b1b..3f37aaf0f2 100644 --- a/Modules/PlanarFigure/CMakeLists.txt +++ b/Modules/PlanarFigure/CMakeLists.txt @@ -1,9 +1,9 @@ MITK_CREATE_MODULE( INCLUDE_DIRS Algorithms DataManagement Interactions IO Rendering DEPENDS MitkSceneSerializationBase - WARNINGS_AS_ERRORS + # WARNINGS_AS_ERRORS ) IF( BUILD_TESTING ) add_subdirectory(Testing) ENDIF() diff --git a/Modules/PlanarFigure/DataManagement/mitkPlanarAngle.cpp b/Modules/PlanarFigure/DataManagement/mitkPlanarAngle.cpp index dd95f19306..4ad9310b33 100644 --- a/Modules/PlanarFigure/DataManagement/mitkPlanarAngle.cpp +++ b/Modules/PlanarFigure/DataManagement/mitkPlanarAngle.cpp @@ -1,175 +1,171 @@ /*=================================================================== 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 "mitkPlanarAngle.h" #include "mitkGeometry2D.h" mitk::PlanarAngle::PlanarAngle() : FEATURE_ID_ANGLE( this->AddFeature( "Angle", "deg" ) ) { // Start with two control points this->ResetNumberOfControlPoints( 2 ); this->SetNumberOfPolyLines(1); this->SetNumberOfHelperPolyLines(1); m_HelperPolyLinesToBePainted->InsertElement( 0, false ); } mitk::PlanarAngle::~PlanarAngle() { } void mitk::PlanarAngle::GeneratePolyLine() { this->ClearPolyLines(); - // Generate poly-line for angle for ( unsigned int i=0; iGetNumberOfControlPoints(); i++ ) - { - mitk::PlanarFigure::PolyLineElement element( this->GetControlPoint( i ), i ); - this->AppendPointToPolyLine( 0, element ); - } + this->AppendPointToPolyLine(0, this->GetControlPoint(i)); } void mitk::PlanarAngle::GenerateHelperPolyLine(double mmPerDisplayUnit, unsigned int displayHeight) { // Generate helper-poly-line for angle if ( this->GetNumberOfControlPoints() < 3) { m_HelperPolyLinesToBePainted->SetElement(0, false); return; //We do not need to draw an angle as there are no two arms yet } this->ClearHelperPolyLines(); const Point2D centerPoint = this->GetControlPoint( 1 ); const Point2D boundaryPointOne = this->GetControlPoint( 0 ); const Point2D boundaryPointTwo = this->GetControlPoint( 2 ); double radius = centerPoint.EuclideanDistanceTo( boundaryPointOne ); if ( radius > centerPoint.EuclideanDistanceTo( boundaryPointTwo ) ) { radius = centerPoint.EuclideanDistanceTo( boundaryPointTwo ); } //Fixed size radius depending on screen size for the angle double nonScalingRadius = displayHeight * mmPerDisplayUnit * 0.05; if (nonScalingRadius > radius) { m_HelperPolyLinesToBePainted->SetElement(0, false); return; //if the arc has a radius that is longer than the shortest arm it should not be painted } m_HelperPolyLinesToBePainted->SetElement(0, true); radius = nonScalingRadius; double angle = this->GetQuantity( FEATURE_ID_ANGLE ); //Determine from which arm the angle should be drawn Vector2D v0 = boundaryPointOne - centerPoint; Vector2D v1 = boundaryPointTwo - centerPoint; Vector2D v2; v2[0] = 1.0; v2[1] = 0.0; v0[0] = v0[0] * cos( 0.001 ) - v0[1] * sin( 0.001 ); //rotate one arm a bit v0[1] = v0[0] * sin( 0.001 ) + v0[1] * cos( 0.001 ); v0.Normalize(); v1.Normalize(); double testAngle = acos( v0 * v1 ); //if the rotated arm is closer to the other arm than before it is the one from which we start drawing //else we start drawing from the other arm (we want to draw in the mathematically positive direction) if( angle > testAngle ) { v1[0] = v0[0] * cos( -0.001 ) - v0[1] * sin( -0.001 ); v1[1] = v0[0] * sin( -0.001 ) + v0[1] * cos( -0.001 ); //We determine if the arm is mathematically forward or backward //assuming we rotate between -pi and pi if ( acos( v0 * v2 ) > acos ( v1 * v2 )) { testAngle = acos( v1 * v2 ); } else { testAngle = -acos( v1 * v2 ); } } else { v0[0] = v1[0] * cos( -0.001 ) - v1[1] * sin( -0.001 ); v0[1] = v1[0] * sin( -0.001 ) + v1[1] * cos( -0.001 ); //We determine if the arm is mathematically forward or backward //assuming we rotate between -pi and pi if ( acos( v0 * v2 ) < acos ( v1 * v2 )) { testAngle = acos( v1 * v2 ); } else { testAngle = -acos( v1 * v2 ); } } // Generate poly-line with 16 segments for ( int t = 0; t < 16; ++t ) { double alpha = (double) t * angle / 15.0 + testAngle; Point2D polyLinePoint; polyLinePoint[0] = centerPoint[0] + radius * cos( alpha ); polyLinePoint[1] = centerPoint[1] + radius * sin( alpha ); - AppendPointToHelperPolyLine( 0, PolyLineElement( polyLinePoint, t ) ); + this->AppendPointToHelperPolyLine(0, polyLinePoint); } } void mitk::PlanarAngle::EvaluateFeaturesInternal() { if ( this->GetNumberOfControlPoints() < 3 ) { // Angle not yet complete. return; } // Calculate angle between lines const Point2D &p0 = this->GetControlPoint( 0 ); const Point2D &p1 = this->GetControlPoint( 1 ); const Point2D &p2 = this->GetControlPoint( 2 ); Vector2D v0 = p1 - p0; Vector2D v1 = p1 - p2; v0.Normalize(); v1.Normalize(); double angle = acos( v0 * v1 ); this->SetQuantity( FEATURE_ID_ANGLE, angle ); } void mitk::PlanarAngle::PrintSelf( std::ostream& os, itk::Indent indent) const { Superclass::PrintSelf( os, indent ); } diff --git a/Modules/PlanarFigure/DataManagement/mitkPlanarArrow.cpp b/Modules/PlanarFigure/DataManagement/mitkPlanarArrow.cpp index 8c62768eba..896807d286 100644 --- a/Modules/PlanarFigure/DataManagement/mitkPlanarArrow.cpp +++ b/Modules/PlanarFigure/DataManagement/mitkPlanarArrow.cpp @@ -1,137 +1,117 @@ /*=================================================================== 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 "mitkPlanarArrow.h" #include "mitkGeometry2D.h" mitk::PlanarArrow::PlanarArrow() : FEATURE_ID_LENGTH( this->AddFeature( "Length", "mm" ) ) { // Directed arrow has two control points this->ResetNumberOfControlPoints( 2 ); m_ArrowTipScaleFactor = -1.0; this->SetNumberOfPolyLines( 1 ); this->SetNumberOfHelperPolyLines( 2 ); - //m_PolyLines->InsertElement( 0, VertexContainerType::New()); - // Create helper polyline object (for drawing the orthogonal orientation line) - //m_HelperPolyLines->InsertElement( 0, VertexContainerType::New()); - //m_HelperPolyLines->InsertElement( 1, VertexContainerType::New()); - //m_HelperPolyLines->ElementAt( 0 )->Reserve( 2 ); - //m_HelperPolyLines->ElementAt( 1 )->Reserve( 2 ); m_HelperPolyLinesToBePainted->InsertElement( 0, false ); m_HelperPolyLinesToBePainted->InsertElement( 1, false ); } mitk::PlanarArrow::~PlanarArrow() { } void mitk::PlanarArrow::GeneratePolyLine() { this->ClearPolyLines(); - this->AppendPointToPolyLine( 0, PolyLineElement( this->GetControlPoint( 0 ), 0 )); - this->AppendPointToPolyLine( 0, PolyLineElement( this->GetControlPoint( 1 ), 0 )); - - // TODO: start line at specified start point... - // Generate poly-line - //m_PolyLines->ElementAt( 0 )->Reserve( 2 ); - //m_PolyLines->ElementAt( 0 )->ElementAt( 0 ) = m_ControlPoints->ElementAt( 0 ); - //m_PolyLines->ElementAt( 0 )->ElementAt( 1 ) = m_ControlPoints->ElementAt( 1 ); + this->AppendPointToPolyLine(0, this->GetControlPoint(0)); + this->AppendPointToPolyLine(0, this->GetControlPoint(1)); } void mitk::PlanarArrow::GenerateHelperPolyLine(double mmPerDisplayUnit, unsigned int displayHeight) { // Generate helper polyline (orientation line orthogonal to first line) // if the third control point is currently being set if ( this->GetNumberOfControlPoints() != 2 ) { m_HelperPolyLinesToBePainted->SetElement( 0, false ); m_HelperPolyLinesToBePainted->SetElement( 1, false ); return; } this->ClearHelperPolyLines(); m_HelperPolyLinesToBePainted->SetElement( 0, true ); m_HelperPolyLinesToBePainted->SetElement( 1, true ); //Fixed size depending on screen size for the angle float scaleFactor = 0.015; if ( m_ArrowTipScaleFactor > 0.0 ) { scaleFactor = m_ArrowTipScaleFactor; } double nonScalingLength = displayHeight * mmPerDisplayUnit * scaleFactor; // Calculate arrow peak const Point2D p1 = this->GetControlPoint( 0 ); const Point2D p2 = this->GetControlPoint( 1 ); - //const Point2D& p1 = m_ControlPoints->ElementAt( 0 ); - //const Point2D& p2 = m_ControlPoints->ElementAt( 1 ); Vector2D n1 = p1 - p2; n1.Normalize(); double degrees = 100.0; Vector2D temp; temp[0] = n1[0] * cos(degrees) - n1[1] * sin(degrees); temp[1] = n1[0] * sin(degrees) + n1[1] * cos(degrees); Vector2D temp2; temp2[0] = n1[0] * cos(-degrees) - n1[1] * sin(-degrees); temp2[1] = n1[0] * sin(-degrees) + n1[1] * cos(-degrees); - this->AppendPointToHelperPolyLine( 0, PolyLineElement( p1, 0 )); - this->AppendPointToHelperPolyLine( 0, PolyLineElement( p1 - temp * nonScalingLength, 0 )); - this->AppendPointToHelperPolyLine( 1, PolyLineElement( p1, 0 )); - this->AppendPointToHelperPolyLine( 1, PolyLineElement( p1 - temp2 * nonScalingLength, 0 )); - - //m_HelperPolyLines->ElementAt( 0 )->ElementAt( 0 ) = p1; - //m_HelperPolyLines->ElementAt( 0 )->ElementAt( 1 ) = p1 - temp * nonScalingLength; - //m_HelperPolyLines->ElementAt( 1 )->ElementAt( 0 ) = p1; - //m_HelperPolyLines->ElementAt( 1 )->ElementAt( 1 ) = p1 - temp2 * nonScalingLength; - + this->AppendPointToHelperPolyLine(0, p1); + this->AppendPointToHelperPolyLine(0, p1 - temp * nonScalingLength); + this->AppendPointToHelperPolyLine(1, p1); + this->AppendPointToHelperPolyLine(1, p1 - temp2 * nonScalingLength); } void mitk::PlanarArrow::EvaluateFeaturesInternal() { // Calculate line length const Point3D &p0 = this->GetWorldControlPoint( 0 ); const Point3D &p1 = this->GetWorldControlPoint( 1 ); double length = p0.EuclideanDistanceTo( p1 ); this->SetQuantity( FEATURE_ID_LENGTH, length ); } void mitk::PlanarArrow::PrintSelf( std::ostream& os, itk::Indent indent) const { Superclass::PrintSelf( os, indent ); } void mitk::PlanarArrow::SetArrowTipScaleFactor( float scale ) { m_ArrowTipScaleFactor = scale; } diff --git a/Modules/PlanarFigure/DataManagement/mitkPlanarBezierCurve.cpp b/Modules/PlanarFigure/DataManagement/mitkPlanarBezierCurve.cpp index d2293d3a67..3690f4abda 100644 --- a/Modules/PlanarFigure/DataManagement/mitkPlanarBezierCurve.cpp +++ b/Modules/PlanarFigure/DataManagement/mitkPlanarBezierCurve.cpp @@ -1,128 +1,128 @@ /*=================================================================== 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 "mitkPlanarBezierCurve.h" #include #include mitk::PlanarBezierCurve::PlanarBezierCurve() : FEATURE_ID_LENGTH(Superclass::AddFeature("Length", "mm")), m_NumberOfSegments(100) { this->ResetNumberOfControlPoints(2); this->SetNumberOfPolyLines(1); this->SetNumberOfHelperPolyLines(1); } mitk::PlanarBezierCurve::~PlanarBezierCurve() { } void mitk::PlanarBezierCurve::EvaluateFeaturesInternal() { std::vector points; points.reserve(m_NumberOfSegments + 1); PolyLineType::const_iterator end = m_PolyLines[0].end(); for (PolyLineType::const_iterator it = m_PolyLines[0].begin(); it != end; ++it) points.push_back(it->Point); double length = 0.0; for (unsigned int i = 0; i < m_NumberOfSegments; ++i) length += points[i].EuclideanDistanceTo(points[i + 1]); this->SetQuantity(FEATURE_ID_LENGTH, length); } unsigned int mitk::PlanarBezierCurve::GetNumberOfSegments() const { return m_NumberOfSegments; } void mitk::PlanarBezierCurve::SetNumberOfSegments(unsigned int numSegments) { m_NumberOfSegments = std::max(1U, numSegments); if (this->IsPlaced()) { this->GeneratePolyLine(); this->Modified(); } } void mitk::PlanarBezierCurve::GenerateHelperPolyLine(double, unsigned int) { this->ClearHelperPolyLines(); unsigned int numHelperPolyLinePoints = m_ControlPoints.size(); for (unsigned int i = 0; i < numHelperPolyLinePoints; ++i) - this->AppendPointToHelperPolyLine(0, PolyLineElement(m_ControlPoints[i], i)); + this->AppendPointToHelperPolyLine(0, m_ControlPoints[i]); } void mitk::PlanarBezierCurve::GeneratePolyLine() { this->ClearPolyLines(); const unsigned int numPolyLinePoints = m_NumberOfSegments + 1; for (unsigned int i = 0; i < numPolyLinePoints; ++i) - this->AppendPointToPolyLine(0, PolyLineElement(this->ComputeDeCasteljauPoint(i / static_cast(m_NumberOfSegments)), i)); + this->AppendPointToPolyLine(0, this->ComputeDeCasteljauPoint(i / static_cast(m_NumberOfSegments))); } mitk::Point2D mitk::PlanarBezierCurve::ComputeDeCasteljauPoint(mitk::ScalarType t) { unsigned int n = m_ControlPoints.size() - 1; if (m_DeCasteljauPoints.size() != n) m_DeCasteljauPoints.resize(n); for (unsigned int i = 0; i < n; ++i) { m_DeCasteljauPoints[i][0] = (1 - t) * m_ControlPoints[i][0] + t * m_ControlPoints[i + 1][0]; m_DeCasteljauPoints[i][1] = (1 - t) * m_ControlPoints[i][1] + t * m_ControlPoints[i + 1][1]; } for (--n; n > 0; --n) { for (unsigned int i = 0; i < n; ++i) { m_DeCasteljauPoints[i][0] = (1 - t) * m_DeCasteljauPoints[i][0] + t * m_DeCasteljauPoints[i + 1][0]; m_DeCasteljauPoints[i][1] = (1 - t) * m_DeCasteljauPoints[i][1] + t * m_DeCasteljauPoints[i + 1][1]; } } return m_DeCasteljauPoints[0]; } unsigned int mitk::PlanarBezierCurve::GetMaximumNumberOfControlPoints() const { return std::numeric_limits::max(); } unsigned int mitk::PlanarBezierCurve::GetMinimumNumberOfControlPoints() const { return 2; } bool mitk::PlanarBezierCurve::IsHelperToBePainted(unsigned int index) { return index == 0 && m_ControlPoints.size() > 2; } diff --git a/Modules/PlanarFigure/DataManagement/mitkPlanarCircle.cpp b/Modules/PlanarFigure/DataManagement/mitkPlanarCircle.cpp index 19b1f9542e..ae44fc627d 100644 --- a/Modules/PlanarFigure/DataManagement/mitkPlanarCircle.cpp +++ b/Modules/PlanarFigure/DataManagement/mitkPlanarCircle.cpp @@ -1,165 +1,165 @@ /*=================================================================== 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 "mitkPlanarCircle.h" #include "mitkGeometry2D.h" #include "mitkProperties.h" mitk::PlanarCircle::PlanarCircle() : FEATURE_ID_RADIUS( this->AddFeature( "Radius", "mm" ) ), FEATURE_ID_DIAMETER( this->AddFeature( "Diameter", "mm" ) ), FEATURE_ID_AREA( this->AddFeature( "Area", "mm2" ) ), m_MinRadius(0), m_MaxRadius(100), m_MinMaxRadiusContraintsActive(false) { // Circle has two control points this->ResetNumberOfControlPoints( 2 ); this->SetNumberOfPolyLines( 1 ); this->SetProperty( "closed", mitk::BoolProperty::New(true) ); } mitk::PlanarCircle::~PlanarCircle() { } bool mitk::PlanarCircle::SetControlPoint( unsigned int index, const Point2D &point, bool /*createIfDoesNotExist*/ ) { // moving center point if(index == 0) { const Point2D ¢erPoint = GetControlPoint( 0 ); Point2D boundaryPoint = GetControlPoint( 1 ); vnl_vector vec = (point.GetVnlVector() - centerPoint.GetVnlVector()); boundaryPoint[0] += vec[0]; boundaryPoint[1] += vec[1]; PlanarFigure::SetControlPoint( 0, point ); PlanarFigure::SetControlPoint( 1, boundaryPoint ); return true; } else if ( index == 1 ) { PlanarFigure::SetControlPoint( index, point ); return true; } return false; } mitk::Point2D mitk::PlanarCircle::ApplyControlPointConstraints(unsigned int index, const Point2D &point) { if ( this->GetGeometry2D() == NULL ) { return point; } Point2D indexPoint; this->GetGeometry2D()->WorldToIndex( point, indexPoint ); BoundingBox::BoundsArrayType bounds = this->GetGeometry2D()->GetBounds(); if ( indexPoint[0] < bounds[0] ) { indexPoint[0] = bounds[0]; } if ( indexPoint[0] > bounds[1] ) { indexPoint[0] = bounds[1]; } if ( indexPoint[1] < bounds[2] ) { indexPoint[1] = bounds[2]; } if ( indexPoint[1] > bounds[3] ) { indexPoint[1] = bounds[3]; } Point2D constrainedPoint; this->GetGeometry2D()->IndexToWorld( indexPoint, constrainedPoint ); if(m_MinMaxRadiusContraintsActive) { if( index != 0) { const Point2D ¢erPoint = this->GetControlPoint(0); double euclideanDinstanceFromCenterToPoint1 = centerPoint.EuclideanDistanceTo(point); Vector2D vectorProjectedPoint; vectorProjectedPoint = point - centerPoint; vectorProjectedPoint.Normalize(); if( euclideanDinstanceFromCenterToPoint1 > m_MaxRadius ) { vectorProjectedPoint *= m_MaxRadius; constrainedPoint = centerPoint; constrainedPoint += vectorProjectedPoint; } else if( euclideanDinstanceFromCenterToPoint1 < m_MinRadius ) { vectorProjectedPoint *= m_MinRadius; constrainedPoint = centerPoint; constrainedPoint += vectorProjectedPoint; } } } return constrainedPoint; } void mitk::PlanarCircle::GeneratePolyLine() { // TODO: start circle at specified boundary point... // clear the PolyLine-Contrainer, it will be reconstructed soon enough... this->ClearPolyLines(); const Point2D ¢erPoint = GetControlPoint( 0 ); const Point2D &boundaryPoint = GetControlPoint( 1 ); double radius = centerPoint.EuclideanDistanceTo( boundaryPoint ); // Generate poly-line with 64 segments for ( int t = 0; t < 64; ++t ) { double alpha = (double) t * vnl_math::pi / 32.0; // construct the new polyline point ... Point2D polyLinePoint; polyLinePoint[0] = centerPoint[0] + radius * cos( alpha ); polyLinePoint[1] = centerPoint[1] + radius * sin( alpha ); // ... and append it to the PolyLine. // No extending supported here, so we can set the index of the PolyLineElement to '0' - AppendPointToPolyLine( 0, PolyLineElement( polyLinePoint, 0 ) ); + this->AppendPointToPolyLine(0, polyLinePoint); } } void mitk::PlanarCircle::GenerateHelperPolyLine(double /*mmPerDisplayUnit*/, unsigned int /*displayHeight*/) { // A circle does not require a helper object } void mitk::PlanarCircle::EvaluateFeaturesInternal() { // Calculate circle radius and area const Point3D &p0 = this->GetWorldControlPoint( 0 ); const Point3D &p1 = this->GetWorldControlPoint( 1 ); double radius = p0.EuclideanDistanceTo( p1 ); double area = vnl_math::pi * radius * radius; this->SetQuantity( FEATURE_ID_RADIUS, radius ); this->SetQuantity( FEATURE_ID_DIAMETER, 2*radius ); this->SetQuantity( FEATURE_ID_AREA, area ); } void mitk::PlanarCircle::PrintSelf( std::ostream& os, itk::Indent indent) const { Superclass::PrintSelf( os, indent ); } diff --git a/Modules/PlanarFigure/DataManagement/mitkPlanarCross.cpp b/Modules/PlanarFigure/DataManagement/mitkPlanarCross.cpp index d08afba80e..5f202fc289 100644 --- a/Modules/PlanarFigure/DataManagement/mitkPlanarCross.cpp +++ b/Modules/PlanarFigure/DataManagement/mitkPlanarCross.cpp @@ -1,350 +1,344 @@ /*=================================================================== 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 "mitkPlanarCross.h" #include "mitkGeometry2D.h" #include "mitkProperties.h" mitk::PlanarCross::PlanarCross() : FEATURE_ID_LONGESTDIAMETER( this->AddFeature( "Longest Axis", "mm" ) ), FEATURE_ID_SHORTAXISDIAMETER( this->AddFeature( "Short Axis", "mm" ) ) { // Cross has two control points at the beginning this->ResetNumberOfControlPoints( 2 ); // Create property for SingleLineMode (default: false) this->SetProperty( "SingleLineMode", mitk::BoolProperty::New( false ) ); // Create helper polyline object (for drawing the orthogonal orientation line) this->SetNumberOfHelperPolyLines( 1 ); m_HelperPolyLinesToBePainted->InsertElement( 0, false ); } mitk::PlanarCross::~PlanarCross() { } void mitk::PlanarCross::SetSingleLineMode( bool singleLineMode ) { this->SetProperty( "SingleLineMode", mitk::BoolProperty::New( singleLineMode ) ); this->Modified(); } bool mitk::PlanarCross::GetSingleLineMode() const { mitk::BoolProperty* singleLineMode = dynamic_cast< mitk::BoolProperty* >( this->GetProperty( "SingleLineMode" ).GetPointer() ); if ( singleLineMode != NULL ) { return singleLineMode->GetValue(); } return false; } bool mitk::PlanarCross::ResetOnPointSelect() { if ( this->GetSingleLineMode() ) { // In single line mode --> nothing to reset return false; } switch ( m_SelectedControlPoint ) { default: // Nothing selected --> nothing to reset return false; case 0: { // Control point 0 selected: exchange points 0 and 1 Point2D tmpPoint = this->GetControlPoint( 0 ); this->SetControlPoint( 0, this->GetControlPoint( 1 ) ); this->SetControlPoint( 1, tmpPoint ); // FALLS THROUGH! } case 1: { // Control point 0 or 1 selected: reset number of control points to two this->ResetNumberOfControlPoints( 2 ); this->SelectControlPoint( 1 ); return true; } case 2: { // Control point 2 selected: replace point 0 with point 3 and point 1 with point 2 this->SetControlPoint( 0, this->GetControlPoint( 3 ) ); this->SetControlPoint( 1, this->GetControlPoint( 2 ) ); // Adjust selected control point, reset number of control points to two this->ResetNumberOfControlPoints( 2 ); this->SelectControlPoint( 1 ); return true; } case 3: { // Control point 3 selected: replace point 0 with point 2 and point 1 with point 3 this->SetControlPoint( 0, this->GetControlPoint( 2 ) ); this->SetControlPoint( 1, this->GetControlPoint( 3 ) ); // Adjust selected control point, reset number of control points to two this->ResetNumberOfControlPoints( 2 ); this->SelectControlPoint( 1 ); return true; } } } unsigned int mitk::PlanarCross::GetNumberOfFeatures() const { if ( this->GetSingleLineMode() || (this->GetNumberOfControlPoints() < 4) ) { return 1; } else { return 2; } } mitk::Point2D mitk::PlanarCross::ApplyControlPointConstraints( unsigned int index, const Point2D& point ) { // Apply spatial constraints from superclass and from this class until the resulting constrained // point converges. Although not an optimal implementation, this iterative approach // helps to respect both constraints from the superclass and from this class. Without this, // situations may occur where control points are constrained by the superclass, but again // moved out of the superclass bounds by the subclass, or vice versa. unsigned int count = 0; // ensures stop of approach if point does not converge in reasonable time Point2D confinedPoint = point; Point2D superclassConfinedPoint; do { superclassConfinedPoint = Superclass::ApplyControlPointConstraints( index, confinedPoint ); confinedPoint = this->InternalApplyControlPointConstraints( index, superclassConfinedPoint ); ++count; } while ( (confinedPoint.EuclideanDistanceTo( superclassConfinedPoint ) > mitk::eps) && (count < 32) ); return confinedPoint; } mitk::Point2D mitk::PlanarCross::InternalApplyControlPointConstraints( unsigned int index, const Point2D& point ) { // Apply constraints depending on current interaction state switch ( index ) { case 2: { // Check if 3rd control point is outside of the range (2D area) defined by the first // line (via the first two control points); if it is outside, clip it to the bounds const Point2D p1 = this->GetControlPoint( 0 ); const Point2D p2 = this->GetControlPoint( 1 ); Vector2D n1 = p2 - p1; n1.Normalize(); Vector2D v1 = point - p1; double dotProduct = n1 * v1; Point2D crossPoint = p1 + n1 * dotProduct;; Vector2D crossVector = point - crossPoint; if ( dotProduct < 0.0 ) { // Out-of-bounds on the left: clip point to left boundary return (p1 + crossVector); } else if ( dotProduct > p2.EuclideanDistanceTo( p1 ) ) { // Out-of-bounds on the right: clip point to right boundary return (p2 + crossVector); } else { // Pass back original point return point; } } case 3: { // Constrain 4th control point so that with the 3rd control point it forms // a line orthogonal to the first line (constraint 1); the 4th control point // must lie on the opposite side of the line defined by the first two control // points than the 3rd control point (constraint 2) const Point2D p1 = this->GetControlPoint( 0 ); const Point2D p2 = this->GetControlPoint( 1 ); const Point2D p3 = this->GetControlPoint( 2 ); // Calculate distance of original point from orthogonal line the corrected // point should lie on to project the point onto this line Vector2D n1 = p2 - p1; n1.Normalize(); Vector2D v1 = point - p3; double dotProduct1 = n1 * v1; Point2D pointOnLine = point - n1 * dotProduct1; // Project new point onto line [p1, p2] Vector2D v2 = pointOnLine - p1; double dotProduct2 = n1 * v2; Point2D crossingPoint = p1 + n1 * dotProduct2; // Determine whether the projected point on the line, or the crossing point should be // used (according to the second constraint in the comment above) if ( (pointOnLine.SquaredEuclideanDistanceTo( p3 ) > crossingPoint.SquaredEuclideanDistanceTo( p3 )) && (pointOnLine.SquaredEuclideanDistanceTo( p3 ) > pointOnLine.SquaredEuclideanDistanceTo( crossingPoint )) ) { return pointOnLine; } else { return crossingPoint; } } default: return point; } } void mitk::PlanarCross::GeneratePolyLine() { - this->SetNumberOfPolyLines( 1 ); - + this->SetNumberOfPolyLines(1); this->ClearPolyLines(); - if ( this->GetNumberOfControlPoints() > 2) - { + if (this->GetNumberOfControlPoints() > 2) this->SetNumberOfPolyLines( 2 ); - } - for ( unsigned int i = 0; i < this->GetNumberOfControlPoints(); ++i ) + for (unsigned int i = 0; i < this->GetNumberOfControlPoints(); ++i) { if (i < 2) - { - this->AppendPointToPolyLine( 0, mitk::PlanarFigure::PolyLineElement( this->GetControlPoint( i ), i ) ); - } + this->AppendPointToPolyLine(0, this->GetControlPoint(i)); + if (i > 1) - { - this->AppendPointToPolyLine( 1, mitk::PlanarFigure::PolyLineElement( this->GetControlPoint( i ), i ) ); - } + this->AppendPointToPolyLine(1, this->GetControlPoint(i)); } } void mitk::PlanarCross::GenerateHelperPolyLine(double /*mmPerDisplayUnit*/, unsigned int /*displayHeight*/) { // Generate helper polyline (orientation line orthogonal to first line) // if the third control point is currently being set if ( this->GetNumberOfControlPoints() != 3 ) { m_HelperPolyLinesToBePainted->SetElement( 0, false ); return; } m_HelperPolyLinesToBePainted->SetElement( 0, true ); this->ClearHelperPolyLines(); // Calculate cross point of first line (p1 to p2) and orthogonal line through // the third control point (p3) const Point2D p1 = this->GetControlPoint( 0 ); const Point2D p2 = this->GetControlPoint( 1 ); const Point2D p3 = this->GetControlPoint( 2 ); Vector2D n1 = p2 - p1; n1.Normalize(); Vector2D v1 = p3 - p1; Point2D crossPoint = p1 + n1 * (n1 * v1); Vector2D v2 = crossPoint - p3; if ( v2.GetNorm() < 1.0 ) { // If third point is on the first line, draw orthogonal "infinite" line // through cross point on line Vector2D v0; v0[0] = n1[1]; v0[1] = -n1[0]; - this->AppendPointToHelperPolyLine( 0, mitk::PlanarFigure::PolyLineElement( p3 - v0 * 10000.0, 0 ) ) ; - this->AppendPointToHelperPolyLine( 0, mitk::PlanarFigure::PolyLineElement( p3 + v0 * 10000.0, 0 ) ) ; + this->AppendPointToHelperPolyLine(0, p3 - v0 * 10000.0); + this->AppendPointToHelperPolyLine(0, p3 + v0 * 10000.0); } else { // Else, draw orthogonal line starting from third point and crossing the // first line, open-ended only on the other side - this->AppendPointToHelperPolyLine( 0, mitk::PlanarFigure::PolyLineElement( p3, 0 ) ) ; - this->AppendPointToHelperPolyLine( 0, mitk::PlanarFigure::PolyLineElement( p3 + v2 * 10000.0, 0 ) ) ; + this->AppendPointToHelperPolyLine(0, p3); + this->AppendPointToHelperPolyLine(0, p3 + v2 * 10000.0); } } void mitk::PlanarCross::EvaluateFeaturesInternal() { // Calculate length of first line const Point3D &p0 = this->GetWorldControlPoint( 0 ); const Point3D &p1 = this->GetWorldControlPoint( 1 ); double l1 = p0.EuclideanDistanceTo( p1 ); // Calculate length of second line double l2 = 0.0; if ( !this->GetSingleLineMode() && (this->GetNumberOfControlPoints() > 3) ) { const Point3D &p2 = this->GetWorldControlPoint( 2 ); const Point3D &p3 = this->GetWorldControlPoint( 3 ); l2 = p2.EuclideanDistanceTo( p3 ); } double longestDiameter; double shortAxisDiameter; if ( l1 > l2 ) { longestDiameter = l1; shortAxisDiameter = l2; } else { longestDiameter = l2; shortAxisDiameter = l1; } this->SetQuantity( FEATURE_ID_LONGESTDIAMETER, longestDiameter ); this->SetQuantity( FEATURE_ID_SHORTAXISDIAMETER, shortAxisDiameter ); } void mitk::PlanarCross::PrintSelf( std::ostream& os, itk::Indent indent) const { Superclass::PrintSelf( os, indent ); } diff --git a/Modules/PlanarFigure/DataManagement/mitkPlanarDoubleEllipse.cpp b/Modules/PlanarFigure/DataManagement/mitkPlanarDoubleEllipse.cpp index 22fb6754d0..68dfab03db 100644 --- a/Modules/PlanarFigure/DataManagement/mitkPlanarDoubleEllipse.cpp +++ b/Modules/PlanarFigure/DataManagement/mitkPlanarDoubleEllipse.cpp @@ -1,251 +1,251 @@ /*=================================================================== 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 "mitkPlanarDoubleEllipse.h" #include #include mitk::PlanarDoubleEllipse::PlanarDoubleEllipse() : FEATURE_ID_MAJOR_AXIS(Superclass::AddFeature("Major Axis", "mm")), FEATURE_ID_MINOR_AXIS(Superclass::AddFeature("Minor Axis", "mm")), FEATURE_ID_THICKNESS(Superclass::AddFeature("Thickness", "mm")), m_NumberOfSegments(64), m_ConstrainCircle(true), m_ConstrainThickness(true) { this->ResetNumberOfControlPoints(4); this->SetNumberOfPolyLines(2); this->SetProperty("closed", mitk::BoolProperty::New(true)); } mitk::PlanarDoubleEllipse::~PlanarDoubleEllipse() { } mitk::Point2D mitk::PlanarDoubleEllipse::ApplyControlPointConstraints(unsigned int index, const Point2D& point) { if (index == 2 && !m_ConstrainCircle) { Point2D centerPoint = this->GetControlPoint(0); Vector2D outerMajorVector = this->GetControlPoint(1) - centerPoint; Vector2D minorDirection; minorDirection[0] = outerMajorVector[1]; minorDirection[1] = -outerMajorVector[0]; minorDirection.Normalize(); double outerMajorRadius = outerMajorVector.GetNorm(); double innerMajorRadius = (this->GetControlPoint(3) - centerPoint).GetNorm(); ScalarType radius = std::max(outerMajorRadius - innerMajorRadius, std::min(centerPoint.EuclideanDistanceTo(point), outerMajorRadius)); return centerPoint + minorDirection * radius; } else if (index == 3 && !m_ConstrainThickness) { Point2D centerPoint = this->GetControlPoint(0); Vector2D outerMajorVector = this->GetControlPoint(1) - centerPoint; double outerMajorRadius = outerMajorVector.GetNorm(); double outerMinorRadius = (this->GetControlPoint(2) - centerPoint).GetNorm(); ScalarType radius = std::max(outerMajorRadius - outerMinorRadius, std::min(centerPoint.EuclideanDistanceTo(point), outerMajorRadius)); outerMajorVector.Normalize(); return centerPoint - outerMajorVector * radius; } return point; } void mitk::PlanarDoubleEllipse::EvaluateFeaturesInternal() { Point2D centerPoint = this->GetControlPoint(0); ScalarType outerMajorRadius = centerPoint.EuclideanDistanceTo(this->GetControlPoint(1)); this->SetQuantity(FEATURE_ID_MAJOR_AXIS, 2 * outerMajorRadius); this->SetQuantity(FEATURE_ID_MINOR_AXIS, 2 * centerPoint.EuclideanDistanceTo(this->GetControlPoint(2))); this->SetQuantity(FEATURE_ID_THICKNESS, outerMajorRadius - centerPoint.EuclideanDistanceTo(this->GetControlPoint(3))); } void mitk::PlanarDoubleEllipse::GenerateHelperPolyLine(double, unsigned int) { } void mitk::PlanarDoubleEllipse::GeneratePolyLine() { this->ClearPolyLines(); Point2D centerPoint = this->GetControlPoint(0); Point2D outerMajorPoint = this->GetControlPoint(1); Vector2D direction = outerMajorPoint - centerPoint; direction.Normalize(); const ScalarType deltaAngle = vnl_math::pi / (m_NumberOfSegments / 2); int start = 0; int end = m_NumberOfSegments; if (direction[1] < 0.0) { direction[0] = -direction[0]; end = m_NumberOfSegments / 2; start = -end; } vnl_matrix_fixed rotation; rotation[1][0] = std::sin(std::acos(direction[0])); rotation[0][0] = direction[0]; rotation[1][1] = direction[0]; rotation[0][1] = -rotation[1][0]; ScalarType outerMajorRadius = centerPoint.EuclideanDistanceTo(outerMajorPoint); ScalarType outerMinorRadius = centerPoint.EuclideanDistanceTo(this->GetControlPoint(2)); ScalarType innerMajorRadius = centerPoint.EuclideanDistanceTo(this->GetControlPoint(3)); ScalarType innerMinorRadius = innerMajorRadius - (outerMajorRadius - outerMinorRadius); ScalarType angle; ScalarType cosAngle; ScalarType sinAngle; vnl_vector_fixed vector; Point2D point; for (int i = start; i < end; ++i) { angle = i * deltaAngle; cosAngle = std::cos(angle); sinAngle = std::sin(angle); vector[0] = outerMajorRadius * cosAngle; vector[1] = outerMinorRadius * sinAngle; vector = rotation * vector; point[0] = centerPoint[0] + vector[0]; point[1] = centerPoint[1] + vector[1]; - this->AppendPointToPolyLine(0, PolyLineElement(point, 0)); + this->AppendPointToPolyLine(0, point); vector[0] = innerMajorRadius * cosAngle; vector[1] = innerMinorRadius * sinAngle; vector = rotation * vector; point[0] = centerPoint[0] + vector[0]; point[1] = centerPoint[1] + vector[1]; - this->AppendPointToPolyLine(1, PolyLineElement(point, 0)); + this->AppendPointToPolyLine(1, point); } } unsigned int mitk::PlanarDoubleEllipse::GetNumberOfSegments() const { return m_NumberOfSegments; } void mitk::PlanarDoubleEllipse::SetNumberOfSegments(unsigned int numSegments) { m_NumberOfSegments = std::max(4U, numSegments); if (this->IsPlaced()) { this->GeneratePolyLine(); this->Modified(); } } unsigned int mitk::PlanarDoubleEllipse::GetMaximumNumberOfControlPoints() const { return 4; } unsigned int mitk::PlanarDoubleEllipse::GetMinimumNumberOfControlPoints() const { return 4; } bool mitk::PlanarDoubleEllipse::SetControlPoint(unsigned int index, const Point2D& point, bool createIfDoesNotExist) { switch (index) { case 0: { Point2D centerPoint = this->GetControlPoint(0); Vector2D vector = point - centerPoint; Superclass::SetControlPoint(0, point, createIfDoesNotExist); Superclass::SetControlPoint(1, this->GetControlPoint(1) + vector, createIfDoesNotExist); Superclass::SetControlPoint(2, this->GetControlPoint(2) + vector, createIfDoesNotExist); Superclass::SetControlPoint(3, this->GetControlPoint(3) + vector, createIfDoesNotExist); break; } case 1: { Vector2D vector = point - this->GetControlPoint(1); Superclass::SetControlPoint(1, point, createIfDoesNotExist); Point2D centerPoint = this->GetControlPoint(0); Vector2D outerMajorVector = point - centerPoint; Vector2D outerMinorVector; outerMinorVector[0] = outerMajorVector[1]; outerMinorVector[1] = -outerMajorVector[0]; if (!m_ConstrainCircle) { outerMinorVector.Normalize(); outerMinorVector *= centerPoint.EuclideanDistanceTo(this->GetControlPoint(2)); } Superclass::SetControlPoint(2, centerPoint + outerMinorVector, createIfDoesNotExist); Vector2D innerMajorVector = outerMajorVector; if (!m_ConstrainThickness) { innerMajorVector.Normalize(); innerMajorVector *= centerPoint.EuclideanDistanceTo(this->GetControlPoint(3) - vector); } Superclass::SetControlPoint(3, centerPoint - innerMajorVector, createIfDoesNotExist); break; } case 2: { m_ConstrainCircle = false; Superclass::SetControlPoint(2, point, createIfDoesNotExist); break; } case 3: { m_ConstrainThickness = false; Superclass::SetControlPoint(3, point, createIfDoesNotExist); break; } default: return false; } return true; } diff --git a/Modules/PlanarFigure/DataManagement/mitkPlanarEllipse.cpp b/Modules/PlanarFigure/DataManagement/mitkPlanarEllipse.cpp index ac13384947..8b1ffda3a7 100644 --- a/Modules/PlanarFigure/DataManagement/mitkPlanarEllipse.cpp +++ b/Modules/PlanarFigure/DataManagement/mitkPlanarEllipse.cpp @@ -1,275 +1,275 @@ /*=================================================================== 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 #include "mitkPlanarEllipse.h" #include "mitkGeometry2D.h" #include "mitkProperties.h" #include mitk::PlanarEllipse::PlanarEllipse() : m_MinRadius(0), m_MaxRadius(100), m_MinMaxRadiusContraintsActive(false), m_TreatAsCircle(true) { // Ellipse has three control points this->ResetNumberOfControlPoints( 4 ); this->SetNumberOfPolyLines( 2 ); this->SetProperty( "closed", mitk::BoolProperty::New(true) ); } mitk::PlanarEllipse::~PlanarEllipse() { } bool mitk::PlanarEllipse::SetControlPoint( unsigned int index, const Point2D &point, bool createIfDoesNotExist ) { if(index == 0) // moving center point and control points accordingly { const Point2D ¢erPoint = GetControlPoint( 0 ); Point2D boundaryPoint1 = GetControlPoint( 1 ); Point2D boundaryPoint2 = GetControlPoint( 2 ); Point2D boundaryPoint3 = GetControlPoint( 3 ); vnl_vector vec = (point.GetVnlVector() - centerPoint.GetVnlVector()); boundaryPoint1[0] += vec[0]; boundaryPoint1[1] += vec[1]; boundaryPoint2[0] += vec[0]; boundaryPoint2[1] += vec[1]; boundaryPoint3[0] += vec[0]; boundaryPoint3[1] += vec[1]; PlanarFigure::SetControlPoint( 0, point, createIfDoesNotExist ); PlanarFigure::SetControlPoint( 1, boundaryPoint1, createIfDoesNotExist ); PlanarFigure::SetControlPoint( 2, boundaryPoint2, createIfDoesNotExist ); PlanarFigure::SetControlPoint( 3, boundaryPoint3, createIfDoesNotExist ); return true; } else if (index < 3) { PlanarFigure::SetControlPoint( index, point, createIfDoesNotExist ); int otherIndex = index+1; if (otherIndex > 2) otherIndex = 1; const Point2D ¢erPoint = GetControlPoint( 0 ); Point2D otherPoint = GetControlPoint( otherIndex ); Point2D point3 = GetControlPoint( 3 ); Vector2D vec1 = point - centerPoint; Vector2D vec2; if (index == 1 && m_TreatAsCircle ) { float x = vec1[0]; vec2[0] = vec1[1]; vec2[1] = x; if (index==1) vec2[0] *= -1; else vec2[1] *= -1; otherPoint = centerPoint+vec2; PlanarFigure::SetControlPoint( otherIndex, otherPoint, createIfDoesNotExist ); float r = centerPoint.EuclideanDistanceTo(otherPoint); // adjust additional third control point Point2D p3 = this->GetControlPoint(3); Vector2D vec3; vec3[0] = p3[0]-centerPoint[0]; vec3[1] = p3[1]-centerPoint[1]; if (vec3[0]!=0 || vec3[1]!=0) { vec3.Normalize(); vec3 *= r; } else { vec3[0] = r; vec3[1] = 0; } point3 = centerPoint + vec3; PlanarFigure::SetControlPoint( 3, point3, createIfDoesNotExist ); } else if ( vec1.GetNorm() > 0 ) { float r = centerPoint.EuclideanDistanceTo(otherPoint); float x = vec1[0]; vec2[0] = vec1[1]; vec2[1] = x; if (index==1) vec2[0] *= -1; else vec2[1] *= -1; vec2.Normalize(); vec2 *= r; if ( vec2.GetNorm() > 0 ) { otherPoint = centerPoint+vec2; PlanarFigure::SetControlPoint( otherIndex, otherPoint, createIfDoesNotExist ); } // adjust third control point Vector2D vec3 = point3 - centerPoint; vec3.Normalize(); double r1 = centerPoint.EuclideanDistanceTo( GetControlPoint( 1 ) ); double r2 = centerPoint.EuclideanDistanceTo( GetControlPoint( 2 ) ); Point2D newPoint = centerPoint + vec3*std::max(r1, r2); PlanarFigure::SetControlPoint( 3, newPoint, createIfDoesNotExist ); m_TreatAsCircle = false; } return true; } else if (index == 3) { Point2D centerPoint = GetControlPoint( 0 ); Vector2D vec3 = point - centerPoint; vec3.Normalize(); double r1 = centerPoint.EuclideanDistanceTo( GetControlPoint( 1 ) ); double r2 = centerPoint.EuclideanDistanceTo( GetControlPoint( 2 ) ); Point2D newPoint = centerPoint + vec3*std::max(r1, r2); PlanarFigure::SetControlPoint( index, newPoint, createIfDoesNotExist ); m_TreatAsCircle = false; return true; } return false; } void mitk::PlanarEllipse::PlaceFigure( const mitk::Point2D &point ) { PlanarFigure::PlaceFigure( point ); m_SelectedControlPoint = 1; } mitk::Point2D mitk::PlanarEllipse::ApplyControlPointConstraints(unsigned int index, const Point2D &point) { return point; Point2D indexPoint; this->GetGeometry2D()->WorldToIndex( point, indexPoint ); BoundingBox::BoundsArrayType bounds = this->GetGeometry2D()->GetBounds(); if ( indexPoint[0] < bounds[0] ) { indexPoint[0] = bounds[0]; } if ( indexPoint[0] > bounds[1] ) { indexPoint[0] = bounds[1]; } if ( indexPoint[1] < bounds[2] ) { indexPoint[1] = bounds[2]; } if ( indexPoint[1] > bounds[3] ) { indexPoint[1] = bounds[3]; } Point2D constrainedPoint; this->GetGeometry2D()->IndexToWorld( indexPoint, constrainedPoint ); if(m_MinMaxRadiusContraintsActive) { if( index != 0) { const Point2D ¢erPoint = this->GetControlPoint(0); double euclideanDinstanceFromCenterToPoint1 = centerPoint.EuclideanDistanceTo(point); Vector2D vectorProjectedPoint; vectorProjectedPoint = point - centerPoint; vectorProjectedPoint.Normalize(); if( euclideanDinstanceFromCenterToPoint1 > m_MaxRadius ) { vectorProjectedPoint *= m_MaxRadius; constrainedPoint = centerPoint; constrainedPoint += vectorProjectedPoint; } else if( euclideanDinstanceFromCenterToPoint1 < m_MinRadius ) { vectorProjectedPoint *= m_MinRadius; constrainedPoint = centerPoint; constrainedPoint += vectorProjectedPoint; } } } return constrainedPoint; } void mitk::PlanarEllipse::GeneratePolyLine() { // clear the PolyLine-Contrainer, it will be reconstructed soon enough... this->ClearPolyLines(); const Point2D ¢erPoint = GetControlPoint( 0 ); const Point2D &boundaryPoint1 = GetControlPoint( 1 ); const Point2D &boundaryPoint2 = GetControlPoint( 2 ); Vector2D dir = boundaryPoint1 - centerPoint; dir.Normalize(); vnl_matrix_fixed rot; // differentiate between clockwise and counterclockwise rotation int start = 0; int end = 64; if (dir[1]<0) { dir[0] = -dir[0]; start = -32; end = 32; } // construct rotation matrix to align ellipse with control point vector rot[0][0] = dir[0]; rot[1][1] = rot[0][0]; rot[1][0] = sin(acos(rot[0][0])); rot[0][1] = -rot[1][0]; double radius1 = centerPoint.EuclideanDistanceTo( boundaryPoint1 ); double radius2 = centerPoint.EuclideanDistanceTo( boundaryPoint2 ); // Generate poly-line with 64 segments for ( int t = start; t < end; ++t ) { double alpha = (double) t * vnl_math::pi / 32.0; // construct the new polyline point ... vnl_vector_fixed< float, 2 > vec; vec[0] = radius1 * cos( alpha ); vec[1] = radius2 * sin( alpha ); vec = rot*vec; Point2D polyLinePoint; polyLinePoint[0] = centerPoint[0] + vec[0]; polyLinePoint[1] = centerPoint[1] + vec[1]; // ... and append it to the PolyLine. // No extending supported here, so we can set the index of the PolyLineElement to '0' - AppendPointToPolyLine( 0, PolyLineElement( polyLinePoint, 0 ) ); + this->AppendPointToPolyLine(0, polyLinePoint); } - AppendPointToPolyLine( 1, PolyLineElement( centerPoint, 0 ) ); - AppendPointToPolyLine( 1, PolyLineElement( GetControlPoint( 3 ), 0 ) ); + this->AppendPointToPolyLine(1, centerPoint); + this->AppendPointToPolyLine(1, this->GetControlPoint(3)); } void mitk::PlanarEllipse::GenerateHelperPolyLine(double /*mmPerDisplayUnit*/, unsigned int /*displayHeight*/) { // A circle does not require a helper object } void mitk::PlanarEllipse::EvaluateFeaturesInternal() { } void mitk::PlanarEllipse::PrintSelf( std::ostream& os, itk::Indent indent) const { Superclass::PrintSelf( os, indent ); } diff --git a/Modules/PlanarFigure/DataManagement/mitkPlanarFigure.cpp b/Modules/PlanarFigure/DataManagement/mitkPlanarFigure.cpp index ca8fe70fc5..264a7973af 100644 --- a/Modules/PlanarFigure/DataManagement/mitkPlanarFigure.cpp +++ b/Modules/PlanarFigure/DataManagement/mitkPlanarFigure.cpp @@ -1,730 +1,757 @@ /*=================================================================== 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 "mitkPlanarFigure.h" #include "mitkGeometry2D.h" #include "mitkProperties.h" #include #include "algorithm" +mitk::PlanarFigure::PolyLineElement::PolyLineElement(Point2D point, int index) + : Point(point), + Index(index) +{ +} + +mitk::PlanarFigure::PolyLineElement::PolyLineElement(const Point2D& point) + : Point(point), + Index(-1) +{ +} + +mitk::PlanarFigure::PolyLineElement::operator mitk::Point2D&() +{ + return Point; +} + +mitk::PlanarFigure::PolyLineElement::operator const mitk::Point2D&() const +{ + return Point; +} mitk::PlanarFigure::PlanarFigure() : m_SelectedControlPoint( -1 ), m_PreviewControlPointVisible( false ), m_FigurePlaced( false ), m_Geometry2D( NULL ), m_PolyLineUpToDate(false), m_HelperLinesUpToDate(false), m_FeaturesUpToDate(false), m_FeaturesMTime( 0 ) { m_HelperPolyLinesToBePainted = BoolContainerType::New(); m_DisplaySize.first = 0.0; m_DisplaySize.second = 0; this->SetProperty( "closed", mitk::BoolProperty::New( false ) ); // Currently only single-time-step geometries are supported this->InitializeTimeGeometry( 1 ); } mitk::PlanarFigure::~PlanarFigure() { } mitk::PlanarFigure::PlanarFigure(const Self& other) : BaseData(other), m_ControlPoints(other.m_ControlPoints), m_NumberOfControlPoints(other.m_NumberOfControlPoints), m_SelectedControlPoint(other.m_SelectedControlPoint), m_PolyLines(other.m_PolyLines), m_HelperPolyLines(other.m_HelperPolyLines), m_HelperPolyLinesToBePainted(other.m_HelperPolyLinesToBePainted->Clone()), m_PreviewControlPoint(other.m_PreviewControlPoint), m_PreviewControlPointVisible(other.m_PreviewControlPointVisible), m_FigurePlaced(other.m_FigurePlaced), m_Geometry2D(other.m_Geometry2D), // do not clone since SetGeometry2D() doesn't clone either m_PolyLineUpToDate(other.m_PolyLineUpToDate), m_HelperLinesUpToDate(other.m_HelperLinesUpToDate), m_FeaturesUpToDate(other.m_FeaturesUpToDate), m_Features(other.m_Features), m_FeaturesMTime(other.m_FeaturesMTime), m_DisplaySize(other.m_DisplaySize) { } void mitk::PlanarFigure::SetGeometry2D( mitk::Geometry2D *geometry ) { this->SetGeometry( geometry ); m_Geometry2D = dynamic_cast(GetGeometry(0));//geometry; } const mitk::Geometry2D *mitk::PlanarFigure::GetGeometry2D() const { return m_Geometry2D; } bool mitk::PlanarFigure::IsClosed() const { mitk::BoolProperty* closed = dynamic_cast< mitk::BoolProperty* >( this->GetProperty( "closed" ).GetPointer() ); if ( closed != NULL ) { return closed->GetValue(); } return false; } void mitk::PlanarFigure::PlaceFigure( const mitk::Point2D& point ) { for ( unsigned int i = 0; i < this->GetNumberOfControlPoints(); ++i ) { m_ControlPoints.push_back( this->ApplyControlPointConstraints( i, point ) ); } m_FigurePlaced = true; m_SelectedControlPoint = 1; } bool mitk::PlanarFigure::AddControlPoint( const mitk::Point2D& point, int position ) { // if we already have the maximum number of control points, do nothing if ( m_NumberOfControlPoints < this->GetMaximumNumberOfControlPoints() ) { // if position has not been defined or position would be the last control point, just append the new one // we also append a new point if we click onto the line between the first two control-points if the second control-point is selected // -> special case for PlanarCross if ( position == -1 || position > (int)m_NumberOfControlPoints-1 || (position == 1 && m_SelectedControlPoint == 2) ) { if ( m_ControlPoints.size() > this->GetMaximumNumberOfControlPoints()-1 ) { // get rid of deprecated control points in the list. This is necessary // as ::ResetNumberOfControlPoints() only sets the member, does not resize the list! m_ControlPoints.resize( this->GetNumberOfControlPoints() ); } m_ControlPoints.push_back( this->ApplyControlPointConstraints( m_NumberOfControlPoints, point ) ); m_SelectedControlPoint = m_NumberOfControlPoints; } else { // insert the point at the given position and set it as selected point ControlPointListType::iterator iter = m_ControlPoints.begin() + position; m_ControlPoints.insert( iter, this->ApplyControlPointConstraints( position, point ) ); for( unsigned int i = 0; i < m_ControlPoints.size(); ++i ) { if( point == m_ControlPoints.at(i) ) { m_SelectedControlPoint = i; } } } // polylines & helperpolylines need to be repainted m_PolyLineUpToDate = false; m_HelperLinesUpToDate = false; m_FeaturesUpToDate = false; // one control point more ++m_NumberOfControlPoints; return true; } else { return false; } } bool mitk::PlanarFigure::SetControlPoint( unsigned int index, const Point2D& point, bool createIfDoesNotExist ) { bool controlPointSetCorrectly = false; if (createIfDoesNotExist) { if ( m_NumberOfControlPoints <= index ) { m_ControlPoints.push_back( this->ApplyControlPointConstraints( index, point ) ); m_NumberOfControlPoints++; } else { m_ControlPoints.at( index ) = this->ApplyControlPointConstraints( index, point ); } controlPointSetCorrectly = true; } else if ( index < m_NumberOfControlPoints ) { m_ControlPoints.at( index ) = this->ApplyControlPointConstraints( index, point ); controlPointSetCorrectly = true; } else { return false; } if ( controlPointSetCorrectly ) { m_PolyLineUpToDate = false; m_HelperLinesUpToDate = false; m_FeaturesUpToDate = false; } return controlPointSetCorrectly; } bool mitk::PlanarFigure::SetCurrentControlPoint( const Point2D& point ) { if ( (m_SelectedControlPoint < 0) || (m_SelectedControlPoint >= (int)m_NumberOfControlPoints) ) { return false; } return this->SetControlPoint(m_SelectedControlPoint, point, false); } unsigned int mitk::PlanarFigure::GetNumberOfControlPoints() const { return m_NumberOfControlPoints; } bool mitk::PlanarFigure::SelectControlPoint( unsigned int index ) { if ( index < this->GetNumberOfControlPoints() ) { m_SelectedControlPoint = index; return true; } else { return false; } } bool mitk::PlanarFigure::DeselectControlPoint() { bool wasSelected = ( m_SelectedControlPoint != -1); m_SelectedControlPoint = -1; return wasSelected; } void mitk::PlanarFigure::SetPreviewControlPoint( const Point2D& point ) { m_PreviewControlPoint = point; m_PreviewControlPointVisible = true; } void mitk::PlanarFigure::ResetPreviewContolPoint() { m_PreviewControlPointVisible = false; } mitk::Point2D mitk::PlanarFigure::GetPreviewControlPoint() { return m_PreviewControlPoint; } bool mitk::PlanarFigure::IsPreviewControlPointVisible() { return m_PreviewControlPointVisible; } mitk::Point2D mitk::PlanarFigure::GetControlPoint( unsigned int index ) const { if ( index < m_NumberOfControlPoints ) { return m_ControlPoints.at( index ); } itkExceptionMacro( << "GetControlPoint(): Invalid index!" ); } mitk::Point3D mitk::PlanarFigure::GetWorldControlPoint( unsigned int index ) const { Point3D point3D; if ( (m_Geometry2D != NULL) && (index < m_NumberOfControlPoints) ) { m_Geometry2D->Map( m_ControlPoints.at( index ), point3D ); return point3D; } itkExceptionMacro( << "GetWorldControlPoint(): Invalid index!" ); } const mitk::PlanarFigure::PolyLineType mitk::PlanarFigure::GetPolyLine(unsigned int index) { mitk::PlanarFigure::PolyLineType polyLine; if ( index > m_PolyLines.size() || !m_PolyLineUpToDate ) { this->GeneratePolyLine(); m_PolyLineUpToDate = true; } return m_PolyLines.at( index );; } const mitk::PlanarFigure::PolyLineType mitk::PlanarFigure::GetPolyLine(unsigned int index) const { return m_PolyLines.at( index ); } void mitk::PlanarFigure::ClearPolyLines() { for ( std::vector::size_type i=0; iGenerateHelperPolyLine(mmPerDisplayUnit, displayHeight); m_HelperLinesUpToDate = true; // store these parameters to be able to check next time if somebody zoomed in or out m_DisplaySize.first = mmPerDisplayUnit; m_DisplaySize.second = displayHeight; } helperPolyLine = m_HelperPolyLines.at(index); } return helperPolyLine; } void mitk::PlanarFigure::ClearHelperPolyLines() { for ( std::vector::size_type i=0; iGeneratePolyLine(); } this->EvaluateFeaturesInternal(); m_FeaturesUpToDate = true; } } void mitk::PlanarFigure::UpdateOutputInformation() { // Bounds are NOT calculated here, since the Geometry2D defines a fixed // frame (= bounds) for the planar figure. Superclass::UpdateOutputInformation(); this->GetTimeGeometry()->Update(); } void mitk::PlanarFigure::SetRequestedRegionToLargestPossibleRegion() { } bool mitk::PlanarFigure::RequestedRegionIsOutsideOfTheBufferedRegion() { return false; } bool mitk::PlanarFigure::VerifyRequestedRegion() { return true; } void mitk::PlanarFigure::SetRequestedRegion(const itk::DataObject * /*data*/ ) { } void mitk::PlanarFigure::ResetNumberOfControlPoints( int numberOfControlPoints ) { // DO NOT resize the list here, will cause crash!! m_NumberOfControlPoints = numberOfControlPoints; } mitk::Point2D mitk::PlanarFigure::ApplyControlPointConstraints( unsigned int /*index*/, const Point2D& point ) { if ( m_Geometry2D == NULL ) { return point; } Point2D indexPoint; m_Geometry2D->WorldToIndex( point, indexPoint ); BoundingBox::BoundsArrayType bounds = m_Geometry2D->GetBounds(); if ( indexPoint[0] < bounds[0] ) { indexPoint[0] = bounds[0]; } if ( indexPoint[0] > bounds[1] ) { indexPoint[0] = bounds[1]; } if ( indexPoint[1] < bounds[2] ) { indexPoint[1] = bounds[2]; } if ( indexPoint[1] > bounds[3] ) { indexPoint[1] = bounds[3]; } Point2D constrainedPoint; m_Geometry2D->IndexToWorld( indexPoint, constrainedPoint ); return constrainedPoint; } unsigned int mitk::PlanarFigure::AddFeature( const char *featureName, const char *unitName ) { unsigned int index = m_Features.size(); Feature newFeature( featureName, unitName ); m_Features.push_back( newFeature ); return index; } void mitk::PlanarFigure::SetFeatureName( unsigned int index, const char *featureName ) { if ( index < m_Features.size() ) { m_Features[index].Name = featureName; } } void mitk::PlanarFigure::SetFeatureUnit( unsigned int index, const char *unitName ) { if ( index < m_Features.size() ) { m_Features[index].Unit = unitName; } } void mitk::PlanarFigure::SetQuantity( unsigned int index, double quantity ) { if ( index < m_Features.size() ) { m_Features[index].Quantity = quantity; } } void mitk::PlanarFigure::ActivateFeature( unsigned int index ) { if ( index < m_Features.size() ) { m_Features[index].Active = true; } } void mitk::PlanarFigure::DeactivateFeature( unsigned int index ) { if ( index < m_Features.size() ) { m_Features[index].Active = false; } } void mitk::PlanarFigure::InitializeTimeGeometry( unsigned int timeSteps ) { mitk::Geometry2D::Pointer geometry2D = mitk::Geometry2D::New(); geometry2D->Initialize(); if ( timeSteps > 1 ) { mitk::ScalarType timeBounds[] = {0.0, 1.0}; geometry2D->SetTimeBounds( timeBounds ); } // The geometry is propagated automatically to all time steps, // if EvenlyTimed is true... ProportionalTimeGeometry::Pointer timeGeometry = ProportionalTimeGeometry::New(); timeGeometry->Initialize(geometry2D, timeSteps); SetTimeGeometry(timeGeometry); } void mitk::PlanarFigure::PrintSelf( std::ostream& os, itk::Indent indent) const { Superclass::PrintSelf( os, indent ); os << indent << this->GetNameOfClass() << ":\n"; if (this->IsClosed()) os << indent << "This figure is closed\n"; else os << indent << "This figure is not closed\n"; os << indent << "Minimum number of control points: " << this->GetMinimumNumberOfControlPoints() << std::endl; os << indent << "Maximum number of control points: " << this->GetMaximumNumberOfControlPoints() << std::endl; os << indent << "Current number of control points: " << this->GetNumberOfControlPoints() << std::endl; os << indent << "Control points:" << std::endl; for ( unsigned int i = 0; i < this->GetNumberOfControlPoints(); ++i ) { //os << indent.GetNextIndent() << i << ": " << m_ControlPoints->ElementAt( i ) << std::endl; os << indent.GetNextIndent() << i << ": " << m_ControlPoints.at( i ) << std::endl; } os << indent << "Geometry:\n"; this->GetGeometry2D()->Print(os, indent.GetNextIndent()); } unsigned short mitk::PlanarFigure::GetPolyLinesSize() { if ( !m_PolyLineUpToDate ) { this->GeneratePolyLine(); m_PolyLineUpToDate = true; } return m_PolyLines.size(); } unsigned short mitk::PlanarFigure::GetHelperPolyLinesSize() { return m_HelperPolyLines.size(); } bool mitk::PlanarFigure::IsHelperToBePainted(unsigned int index) { return m_HelperPolyLinesToBePainted->GetElement( index ); } bool mitk::PlanarFigure::ResetOnPointSelect() { return false; } void mitk::PlanarFigure::RemoveControlPoint( unsigned int index ) { if ( index > m_ControlPoints.size() ) return; if ( (m_ControlPoints.size() -1) < this->GetMinimumNumberOfControlPoints() ) return; ControlPointListType::iterator iter; iter = m_ControlPoints.begin() + index; m_ControlPoints.erase( iter ); m_PolyLineUpToDate = false; m_HelperLinesUpToDate = false; m_FeaturesUpToDate = false; --m_NumberOfControlPoints; } void mitk::PlanarFigure::RemoveLastControlPoint() { RemoveControlPoint( m_ControlPoints.size()-1 ); } void mitk::PlanarFigure::DeepCopy(Self::Pointer oldFigure) { //DeepCopy only same types of planar figures //Notice to get typeid polymorph you have to use the *operator if(typeid(*oldFigure) != typeid(*this)) { itkExceptionMacro( << "DeepCopy(): Inconsistent type of source (" << typeid(*oldFigure).name() << ") and destination figure (" << typeid(*this).name() << ")!" ); return; } m_ControlPoints.clear(); this->ClearPolyLines(); this->ClearHelperPolyLines(); // clone base data members SetPropertyList(oldFigure->GetPropertyList()->Clone()); /// deep copy members m_FigurePlaced = oldFigure->m_FigurePlaced; m_SelectedControlPoint = oldFigure->m_SelectedControlPoint; m_FeaturesMTime = oldFigure->m_FeaturesMTime; m_Features = oldFigure->m_Features; m_NumberOfControlPoints = oldFigure->m_NumberOfControlPoints; //copy geometry 2D of planar figure Geometry2D::Pointer affineGeometry = oldFigure->m_Geometry2D->Clone(); SetGeometry2D(affineGeometry.GetPointer()); for(unsigned long index=0; index < oldFigure->GetNumberOfControlPoints(); index++) { m_ControlPoints.push_back( oldFigure->GetControlPoint( index )); } //After setting the control points we can generate the polylines this->GeneratePolyLine(); } void mitk::PlanarFigure::SetNumberOfPolyLines( unsigned int numberOfPolyLines ) { m_PolyLines.resize(numberOfPolyLines); } void mitk::PlanarFigure::SetNumberOfHelperPolyLines( unsigned int numberOfHerlperPolyLines ) { m_HelperPolyLines.resize(numberOfHerlperPolyLines); } void mitk::PlanarFigure::AppendPointToPolyLine( unsigned int index, PolyLineElement element ) { if ( index < m_PolyLines.size() ) { - m_PolyLines.at( index ).push_back( element ); + if(element.Index == -1) + element.Index = m_PolyLines[index].size(); + + m_PolyLines[index].push_back(element); m_PolyLineUpToDate = false; } else { MITK_ERROR << "Tried to add point to PolyLine " << index+1 << ", although only " << m_PolyLines.size() << " exists"; } } void mitk::PlanarFigure::AppendPointToHelperPolyLine( unsigned int index, PolyLineElement element ) { if ( index < m_HelperPolyLines.size() ) { - m_HelperPolyLines.at( index ).push_back( element ); + if(element.Index == -1) + element.Index = m_HelperPolyLines[index].size(); + + m_HelperPolyLines[index].push_back(element); m_HelperLinesUpToDate = false; } else { MITK_ERROR << "Tried to add point to HelperPolyLine " << index+1 << ", although only " << m_HelperPolyLines.size() << " exists"; } } diff --git a/Modules/PlanarFigure/DataManagement/mitkPlanarFigure.h b/Modules/PlanarFigure/DataManagement/mitkPlanarFigure.h index d13bbcf44f..5d54a93869 100644 --- a/Modules/PlanarFigure/DataManagement/mitkPlanarFigure.h +++ b/Modules/PlanarFigure/DataManagement/mitkPlanarFigure.h @@ -1,418 +1,423 @@ /*=================================================================== 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_PLANAR_FIGURE_H_ #define _MITK_PLANAR_FIGURE_H_ #include #include "mitkBaseData.h" #include "mitkCommon.h" #include namespace mitk { class Geometry2D; /** * \brief Base-class for geometric planar (2D) figures, such as * lines, circles, rectangles, polygons, etc. * * \warning Currently does not support time-resolved data handling * * Behavior and appearance of PlanarFigures are controlled by various properties; for a detailed * list of appearance properties see mitk::PlanarFigureMapper2D * * The following properties control general PlanarFigure behavior: * *
    *
  • "selected": true if the planar figure is selected *
  • "planarfigure.ishovering": true if the mouse "hovers" over the planar figure *
  • "planarfigure.iseditable": true if the planar figure can be edited (otherwise, * it can only be picked/selected, but its control points cannot be edited); default is true *
  • "planarfigure.isextendable": true if new control points can be inserted into the list of control points; * default is false *
* * * TODO: Implement local 2D transform (including center of rotation...) * */ class MitkPlanarFigure_EXPORT PlanarFigure : public BaseData { public: mitkClassMacro( PlanarFigure, BaseData ) itkCloneMacro( Self ) + /** \brief Treat as Point2D by implicitly using conversion operators. + * + * \deprecatedSince{2014_06} "struct PolyLineElement {...};" will be changed to "typedef Point2D PolyLineElement;". + */ struct PolyLineElement { - PolyLineElement( Point2D point, int index ) - : Point( point ), Index( index ) - { - }; + DEPRECATED(PolyLineElement(Point2D point, int index)); + PolyLineElement(const Point2D& point); + + operator Point2D&(); + operator const Point2D&() const; - Point2D Point; - int Index; + DEPRECATED(Point2D Point); + DEPRECATED(int Index); }; typedef itk::VectorContainer< unsigned long, bool> BoolContainerType; typedef std::deque< Point2D > ControlPointListType; - typedef std::list< PolyLineElement > PolyLineType; + typedef std::vector< PolyLineElement > PolyLineType; /** \brief Sets the 2D geometry on which this figure will be placed. * * In most cases, this is a Geometry already owned by another object, e.g. * describing the slice of the image on which measurements will be * performed. */ virtual void SetGeometry2D( mitk::Geometry2D *geometry ); /** \brief Returns (previously set) 2D geometry of this figure. */ virtual const Geometry2D *GetGeometry2D() const; /** \brief True if the planar figure is closed. * * Default is false. The "closed" boolean property must be set in sub-classes. */ virtual bool IsClosed() const; /** \brief True if the planar figure has been placed (and can be * displayed/interacted with). */ virtual bool IsPlaced() const { return m_FigurePlaced; }; /** \brief Place figure at the given point (in 2D index coordinates) onto * the given 2D geometry. * * By default, the first two control points of the figure are set to the * passed point. Further points can be set via AddControlPoint(), if the * current number of control points is below the maximum number of control * points. * * Can be re-implemented in sub-classes as needed. */ virtual void PlaceFigure( const Point2D& point ); /** * \brief Adds / inserts new control-points * * This method adds a new control-point with the coordinates defined by point at the given index. * If 'index' == -1 or index is greater than the number of control-points the new point is appended * to the back of the list of control points. * If a control-point already exists for 'index', an additional point is inserted at that position. * It is not possible to add more points if the maximum number of control-points (GetMaximumNumberOfControlPoints()) * has been reached. */ virtual bool AddControlPoint( const Point2D& point, int index = -1 ); virtual bool SetControlPoint( unsigned int index, const Point2D& point, bool createIfDoesNotExist = false); virtual bool SetCurrentControlPoint( const Point2D& point ); /** \brief Returns the current number of 2D control points defining this figure. */ unsigned int GetNumberOfControlPoints() const; /** \brief Returns the minimum number of control points needed to represent * this figure. * * Must be implemented in sub-classes. */ virtual unsigned int GetMinimumNumberOfControlPoints() const = 0; /** \brief Returns the maximum number of control points allowed for * this figure (e.g. 3 for triangles). * * Must be implemented in sub-classes. */ virtual unsigned int GetMaximumNumberOfControlPoints() const = 0; /** \brief Selects currently active control points. */ virtual bool SelectControlPoint( unsigned int index ); /** \brief Deselect control point; no control point active. */ virtual bool DeselectControlPoint(); /** \brief Return currently selected control point. */ virtual int GetSelectedControlPoint() const { return m_SelectedControlPoint; } /** \brief Returns specified control point in 2D world coordinates. */ Point2D GetControlPoint( unsigned int index ) const; /** \brief Returns specified control point in world coordinates. */ Point3D GetWorldControlPoint( unsigned int index ) const; /** \brief Returns the polyline representing the planar figure * (for rendering, measurements, etc.). */ const PolyLineType GetPolyLine(unsigned int index); /** \brief Returns the polyline representing the planar figure * (for rendering, measurments, etc.). */ const PolyLineType GetPolyLine(unsigned int index) const; /** \brief Returns the polyline that should be drawn the same size at every scale * (for text, angles, etc.). */ const PolyLineType GetHelperPolyLine( unsigned int index, double mmPerDisplayUnit, unsigned int displayHeight ); /** \brief Sets the position of the PreviewControlPoint. Automatically sets it visible.*/ void SetPreviewControlPoint( const Point2D& point ); /** \brief Marks the PreviewControlPoint as invisible.*/ void ResetPreviewContolPoint(); /** \brief Returns whether or not the PreviewControlPoint is visible.*/ bool IsPreviewControlPointVisible(); /** \brief Returns the coordinates of the PreviewControlPoint. */ Point2D GetPreviewControlPoint(); /** \brief Returns the number of features available for this PlanarFigure * (such as, radius, area, ...). */ virtual unsigned int GetNumberOfFeatures() const; /** \brief Returns the name (identifier) of the specified features. */ const char *GetFeatureName( unsigned int index ) const; /** \brief Returns the physical unit of the specified features. */ const char *GetFeatureUnit( unsigned int index ) const; /** Returns quantity of the specified feature (e.g., length, radius, * area, ... ) */ double GetQuantity( unsigned int index ) const; /** \brief Returns true if the feature with the specified index exists and * is active (an inactive feature may e.g. be the area of a non-closed * polygon. */ bool IsFeatureActive( unsigned int index ) const; /** \brief Returns true if the feature with the specified index exists and is set visible */ bool IsFeatureVisible( unsigned int index ) const; /** \brief Defines if the feature with the specified index will be shown as an * overlay in the RenderWindow */ void SetFeatureVisible( unsigned int index, bool visible ); /** \brief Calculates quantities of all features of this planar figure. */ virtual void EvaluateFeatures(); /** \brief Intherited from parent */ virtual void UpdateOutputInformation(); /** \brief Intherited from parent */ virtual void SetRequestedRegionToLargestPossibleRegion(); /** \brief Intherited from parent */ virtual bool RequestedRegionIsOutsideOfTheBufferedRegion(); /** \brief Intherited from parent */ virtual bool VerifyRequestedRegion(); /** \brief Intherited from parent */ virtual void SetRequestedRegion( const itk::DataObject *data); /** \brief Returns the current number of polylines */ virtual unsigned short GetPolyLinesSize(); /** \brief Returns the current number of helperpolylines */ virtual unsigned short GetHelperPolyLinesSize(); /** \brief Returns whether a helper polyline should be painted or not */ virtual bool IsHelperToBePainted(unsigned int index); /** \brief Returns true if the planar figure is reset to "add points" mode * when a point is selected. * * Default return value is false. Subclasses can overwrite this method and * execute any reset / initialization statements required. */ virtual bool ResetOnPointSelect(); /** \brief removes the point with the given index from the list of controlpoints. */ virtual void RemoveControlPoint( unsigned int index ); /** \brief Removes last control point */ virtual void RemoveLastControlPoint(); /** \brief Copies contents and state of a figre provided as parameter to the current object. * * Requires a matching type of both figures. * * \note Deprecated, use Clone() instead. */ DEPRECATED(void DeepCopy(Self::Pointer oldFigure)); /** \brief Allow sub-classes to apply constraints on control points. * * Sub-classes can define spatial constraints to certain control points by * overwriting this method and returning a constrained point. By default, * the points are constrained by the image bounds. */ virtual Point2D ApplyControlPointConstraints( unsigned int /*index*/, const Point2D& point ); protected: PlanarFigure(); virtual ~PlanarFigure(); PlanarFigure(const Self& other); /** \brief Set the initial number of control points of the planar figure */ void ResetNumberOfControlPoints( int numberOfControlPoints ); /** Adds feature (e.g., circumference, radius, angle, ...) to feature vector * of a planar figure object and returns integer ID for the feature element. * Should be called in sub-class constructors. */ virtual unsigned int AddFeature( const char *featureName, const char *unitName ); /** Sets the name of the specified feature. INTERNAL METHOD. */ void SetFeatureName( unsigned int index, const char *featureName ); /** Sets the physical unit of the specified feature. INTERNAL METHOD. */ void SetFeatureUnit( unsigned int index, const char *unitName ); /** Sets quantity of the specified feature. INTERNAL METHOD. */ void SetQuantity( unsigned int index, double quantity ); /** Sets the specified feature as active. INTERAL METHOD. */ void ActivateFeature( unsigned int index ); /** Sets the specified feature as active. INTERAL METHOD. */ void DeactivateFeature( unsigned int index ); /** \brief Generates the poly-line representation of the planar figure. * Must be implemented in sub-classes. */ virtual void GeneratePolyLine() = 0; /** \brief Generates the poly-lines that should be drawn the same size regardless of zoom. * Must be implemented in sub-classes. */ virtual void GenerateHelperPolyLine(double mmPerDisplayUnit, unsigned int displayHeight) = 0; /** \brief Calculates quantities of all features of this planar figure. * Must be implemented in sub-classes. */ virtual void EvaluateFeaturesInternal() = 0; /** \brief Initializes the TimeGeometry describing the (time-resolved) * geometry of this figure. Note that each time step holds one Geometry2D. */ virtual void InitializeTimeGeometry( unsigned int timeSteps = 1 ); /** \brief defines the number of PolyLines that will be available */ void SetNumberOfPolyLines( unsigned int numberOfPolyLines ); /** \brief Append a point to the PolyLine # index */ void AppendPointToPolyLine( unsigned int index, PolyLineElement element ); /** \brief clears the list of PolyLines. Call before re-calculating a new Polyline. */ void ClearPolyLines(); /** \brief defines the number of HelperPolyLines that will be available */ void SetNumberOfHelperPolyLines( unsigned int numberOfHelperPolyLines ); /** \brief Append a point to the HelperPolyLine # index */ void AppendPointToHelperPolyLine( unsigned int index, PolyLineElement element ); /** \brief clears the list of HelperPolyLines. Call before re-calculating a new HelperPolyline. */ void ClearHelperPolyLines(); virtual void PrintSelf( std::ostream& os, itk::Indent indent ) const; ControlPointListType m_ControlPoints; unsigned int m_NumberOfControlPoints; // Currently selected control point; -1 means no point selected int m_SelectedControlPoint; std::vector m_PolyLines; std::vector m_HelperPolyLines; BoolContainerType::Pointer m_HelperPolyLinesToBePainted; // this point is used to store the coordiantes an additional 'ControlPoint' that is rendered // when the mouse cursor is above the figure (and not a control-point) and when the // property 'planarfigure.isextendable' is set to true Point2D m_PreviewControlPoint; bool m_PreviewControlPointVisible; bool m_FigurePlaced; private: // not implemented to prevent PlanarFigure::New() calls which would create an itk::Object. static Pointer New(); struct Feature { Feature( const char *name, const char *unit ) : Name( name ), Unit( unit ), Quantity( 0.0 ), Active( true ), Visible( true ) { } std::string Name; std::string Unit; double Quantity; bool Active; bool Visible; }; virtual itk::LightObject::Pointer InternalClone() const = 0; Geometry2D *m_Geometry2D; bool m_PolyLineUpToDate; bool m_HelperLinesUpToDate; bool m_FeaturesUpToDate; // Vector of features available for this geometric figure typedef std::vector< Feature > FeatureVectorType; FeatureVectorType m_Features; unsigned long m_FeaturesMTime; // this pair is used to store the mmInDisplayUnits (m_DisplaySize.first) and the displayHeight (m_DisplaySize.second) // that the helperPolyLines have been calculated for. // It's used to determine whether or not GetHelperPolyLine() needs to recalculate the HelperPolyLines. std::pair m_DisplaySize; }; } // namespace mitk #endif //_MITK_PLANAR_FIGURE_H_ diff --git a/Modules/PlanarFigure/DataManagement/mitkPlanarFourPointAngle.cpp b/Modules/PlanarFigure/DataManagement/mitkPlanarFourPointAngle.cpp index fb7ec96ab7..0aa9aee2e1 100644 --- a/Modules/PlanarFigure/DataManagement/mitkPlanarFourPointAngle.cpp +++ b/Modules/PlanarFigure/DataManagement/mitkPlanarFourPointAngle.cpp @@ -1,84 +1,79 @@ /*=================================================================== 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 "mitkPlanarFourPointAngle.h" #include "mitkGeometry2D.h" mitk::PlanarFourPointAngle::PlanarFourPointAngle() : FEATURE_ID_ANGLE( this->AddFeature( "Angle", "deg" ) ) { // Four point angle has two control points this->ResetNumberOfControlPoints( 2 ); this->SetNumberOfPolyLines( 2 ); } mitk::PlanarFourPointAngle::~PlanarFourPointAngle() { } void mitk::PlanarFourPointAngle::GeneratePolyLine() { this->ClearPolyLines(); - // TODO: start line at specified start point... - // Generate poly-line - for ( unsigned int i = 0; i < this->GetNumberOfControlPoints(); ++i ) - { - int index = i/2; - this->AppendPointToPolyLine( index, PolyLineElement( GetControlPoint( i ), i ) ); - } + for (unsigned int i = 0; i < this->GetNumberOfControlPoints(); ++i) + this->AppendPointToPolyLine(i / 2, this->GetControlPoint(i)); } void mitk::PlanarFourPointAngle::GenerateHelperPolyLine(double /*mmPerDisplayUnit*/, unsigned int /*displayHeight*/) { // Generate helper-poly-line for an four point angle // Need to discuss a sensible implementation } void mitk::PlanarFourPointAngle::EvaluateFeaturesInternal() { if ( this->GetNumberOfControlPoints() < 4 ) { // Angle not yet complete. return; } // Calculate angle between lines const Point2D &p0 = this->GetControlPoint( 0 ); const Point2D &p1 = this->GetControlPoint( 1 ); const Point2D &p2 = this->GetControlPoint( 2 ); const Point2D &p3 = this->GetControlPoint( 3 ); Vector2D v0 = p1 - p0; Vector2D v1 = p3 - p2; v0.Normalize(); v1.Normalize(); double angle = acos( v0 * v1 ); this->SetQuantity( FEATURE_ID_ANGLE, angle ); } void mitk::PlanarFourPointAngle::PrintSelf( std::ostream& os, itk::Indent indent) const { Superclass::PrintSelf( os, indent ); } diff --git a/Modules/PlanarFigure/DataManagement/mitkPlanarLine.cpp b/Modules/PlanarFigure/DataManagement/mitkPlanarLine.cpp index 288c98643c..75f3c3d7e6 100644 --- a/Modules/PlanarFigure/DataManagement/mitkPlanarLine.cpp +++ b/Modules/PlanarFigure/DataManagement/mitkPlanarLine.cpp @@ -1,66 +1,65 @@ /*=================================================================== 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 "mitkPlanarLine.h" #include "mitkGeometry2D.h" mitk::PlanarLine::PlanarLine() : FEATURE_ID_LENGTH( this->AddFeature( "Length", "mm" ) ) { // Line has two control points this->ResetNumberOfControlPoints( 2 ); this->SetNumberOfPolyLines( 1 ); } mitk::PlanarLine::~PlanarLine() { } void mitk::PlanarLine::GeneratePolyLine() { this->ClearPolyLines(); - // TODO: start line at specified start point... - // Generate poly-line - this->AppendPointToPolyLine( 0 , mitk::PlanarFigure::PolyLineElement( this->GetControlPoint(0), 0) ); - this->AppendPointToPolyLine( 0 , mitk::PlanarFigure::PolyLineElement( this->GetControlPoint(1), 0) ); + + this->AppendPointToPolyLine(0, this->GetControlPoint(0)); + this->AppendPointToPolyLine(0, this->GetControlPoint(1)); } void mitk::PlanarLine::GenerateHelperPolyLine(double /*mmPerDisplayUnit*/, unsigned int /*displayHeight*/) { // A line does not require a helper object } void mitk::PlanarLine::EvaluateFeaturesInternal() { // Calculate line length const Point3D &p0 = this->GetWorldControlPoint( 0 ); const Point3D &p1 = this->GetWorldControlPoint( 1 ); double length = p0.EuclideanDistanceTo( p1 ); this->SetQuantity( FEATURE_ID_LENGTH, length ); } void mitk::PlanarLine::PrintSelf( std::ostream& os, itk::Indent indent) const { Superclass::PrintSelf( os, indent ); } diff --git a/Modules/PlanarFigure/DataManagement/mitkPlanarPolygon.cpp b/Modules/PlanarFigure/DataManagement/mitkPlanarPolygon.cpp index 5ac118bd61..26c6ba3d1d 100644 --- a/Modules/PlanarFigure/DataManagement/mitkPlanarPolygon.cpp +++ b/Modules/PlanarFigure/DataManagement/mitkPlanarPolygon.cpp @@ -1,284 +1,280 @@ /*=================================================================== 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 "mitkPlanarPolygon.h" #include "mitkGeometry2D.h" #include "mitkProperties.h" // stl related includes #include mitk::PlanarPolygon::PlanarPolygon() : FEATURE_ID_CIRCUMFERENCE( this->AddFeature( "Circumference", "mm" ) ), FEATURE_ID_AREA( this->AddFeature( "Area", "mm2" ) ) { // Polygon has at least two control points this->ResetNumberOfControlPoints( 2 ); this->SetNumberOfPolyLines( 1 ); // Polygon is closed by default this->SetProperty( "closed", mitk::BoolProperty::New( true ) ); this->SetProperty( "subdivision", mitk::BoolProperty::New( false ) ); } mitk::PlanarPolygon::~PlanarPolygon() { } void mitk::PlanarPolygon::SetClosed( bool closed ) { this->SetProperty( "closed", mitk::BoolProperty::New( closed ) ); if ( !closed ) { // For non-closed polygons: use "Length" as feature name; disable area this->SetFeatureName( FEATURE_ID_CIRCUMFERENCE, "Length" ); this->DeactivateFeature( FEATURE_ID_AREA ); } else { // For closed polygons: use "Circumference" as feature name; enable area this->SetFeatureName( FEATURE_ID_CIRCUMFERENCE, "Circumference" ); this->ActivateFeature( FEATURE_ID_AREA ); } this->Modified(); } void mitk::PlanarPolygon::GeneratePolyLine() { this->ClearPolyLines(); - for ( ControlPointListType::size_type i=0; iAppendPointToPolyLine( 0, elem ); - } + for (ControlPointListType::size_type i = 0; i < m_ControlPoints.size(); ++i) + this->AppendPointToPolyLine(0, this->GetControlPoint(i)); } void mitk::PlanarPolygon::GenerateHelperPolyLine(double /*mmPerDisplayUnit*/, unsigned int /*displayHeight*/) { // A polygon does not require helper objects } void mitk::PlanarPolygon::EvaluateFeaturesInternal() { // Calculate circumference double circumference = 0.0; unsigned int i,j; ControlPointListType polyLinePoints; polyLinePoints.clear(); PolyLineType::iterator iter; for( iter = m_PolyLines[0].begin(); iter != m_PolyLines[0].end(); ++iter ) { polyLinePoints.push_back((*iter).Point); } if(polyLinePoints.empty()) return; for ( i = 0; i <(polyLinePoints.size()-1); ++i ) { circumference += polyLinePoints[i].EuclideanDistanceTo( polyLinePoints[i + 1] ); } if ( this->IsClosed() ) { circumference += polyLinePoints[i].EuclideanDistanceTo( polyLinePoints.front() ); } this->SetQuantity( FEATURE_ID_CIRCUMFERENCE, circumference ); // Calculate polygon area (if closed) double area = 0.0; bool intersection = false; if ( this->IsClosed() && (this->GetGeometry2D() != NULL) ) { // does PlanarPolygon overlap/intersect itself? unsigned int numberOfPoints = polyLinePoints.size(); if( numberOfPoints >= 4) { for ( i = 0; i < (numberOfPoints - 1); ++i ) { // line 1 Point2D p0 = polyLinePoints[i]; Point2D p1 = polyLinePoints[i + 1]; // check for intersection with all other lines for (j = i+1; j < (numberOfPoints - 1); ++j ) { Point2D p2 = polyLinePoints[j]; Point2D p3 = polyLinePoints[j + 1]; intersection = CheckForLineIntersection(p0,p1,p2,p3); if (intersection) break; } if (intersection) break; // only because the inner loop might have changed "intersection" // last line from p_x to p_0 Point2D p2 = polyLinePoints.front(); Point2D p3 = polyLinePoints.back(); intersection = CheckForLineIntersection(p0,p1,p2,p3); if (intersection) break; } } // calculate area for ( i = 0; i < polyLinePoints.size(); ++i ) { Point2D p0 = polyLinePoints[i]; Point2D p1 = polyLinePoints[ (i + 1) % polyLinePoints.size() ]; area += p0[0] * p1[1] - p1[0] * p0[1]; } area /= 2.0; } // set area if appropiate (i.e. closed and not intersected) if(this->IsClosed() && !intersection) { SetQuantity( FEATURE_ID_AREA, fabs( area ) ); this->ActivateFeature( FEATURE_ID_AREA ); } else { SetQuantity( FEATURE_ID_AREA, 0 ); this->DeactivateFeature( FEATURE_ID_AREA ); } } void mitk::PlanarPolygon::PrintSelf( std::ostream& os, itk::Indent indent) const { Superclass::PrintSelf( os, indent ); if (this->IsClosed()) os << indent << "Polygon is closed\n"; else os << indent << "Polygon is not closed\n"; } // based on // http://flassari.is/2008/11/line-line-intersection-in-cplusplus/ bool mitk::PlanarPolygon::CheckForLineIntersection( const mitk::Point2D& p1, const mitk::Point2D& p2, const mitk::Point2D& p3, const mitk::Point2D& p4, Point2D& intersection ) const { // do not check for intersections with control points if(p1 == p2 || p1 == p3 || p1 == p4 || p2 == p3 || p2 == p4 || p3 == p4) return false; // Store the values for fast access and easy // equations-to-code conversion double x1 = p1[0], x2 = p2[0], x3 = p3[0], x4 = p4[0]; double y1 = p1[1], y2 = p2[1], y3 = p3[1], y4 = p4[1]; double d = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4); // If d is zero, there is no intersection //if (d < mitk::eps) return false; if (d == 0) return false; // Get the x and y double pre = (x1*y2 - y1*x2); double post = (x3*y4 - y3*x4); double x = ( pre * (x3 - x4) - (x1 - x2) * post ) / d; double y = ( pre * (y3 - y4) - (y1 - y2) * post ) / d; double tolerance = 0.001; // Check if the x coordinates are within both lines, including tolerance if ( x < ( std::min(x1, x2) - tolerance ) || x > ( std::max(x1, x2) + tolerance ) || x < ( std::min(x3, x4) - tolerance ) || x > ( std::max(x3, x4) + tolerance ) ) { return false; } // Check if the y coordinates are within both lines, including tolerance if ( y < ( std::min(y1, y2) - tolerance ) || y > ( std::max(y1, y2) + tolerance ) || y < ( std::min(y3, y4) - tolerance ) || y > ( std::max(y3, y4) + tolerance ) ) { return false; } // point of intersection Point2D ret; ret[0] = x; ret[1] = y; intersection = ret; return true; } bool mitk::PlanarPolygon::CheckForLineIntersection( const mitk::Point2D& p1, const mitk::Point2D& p2, const mitk::Point2D& p3, const mitk::Point2D& p4 ) const { mitk::Point2D intersection; return mitk::PlanarPolygon::CheckForLineIntersection( p1, p2, p3, p4, intersection ); } std::vector mitk::PlanarPolygon::CheckForLineIntersection( const mitk::Point2D& p1, const mitk::Point2D& p2 ) const { std::vector intersectionList; ControlPointListType polyLinePoints; PolyLineType tempList = m_PolyLines[0]; PolyLineType::iterator iter; for( iter = tempList.begin(); iter != tempList.end(); ++iter ) { polyLinePoints.push_back((*iter).Point); } for ( ControlPointListType::size_type i=0; iIsClosed() ) { mitk::Point2D intersection, lastControlPoint, firstControlPoint; lastControlPoint = polyLinePoints.back(); firstControlPoint = polyLinePoints.front(); if ( mitk::PlanarPolygon::CheckForLineIntersection( lastControlPoint, firstControlPoint, p1, p2, intersection ) ) { intersectionList.push_back( intersection ); } } return intersectionList; } diff --git a/Modules/PlanarFigure/DataManagement/mitkPlanarRectangle.cpp b/Modules/PlanarFigure/DataManagement/mitkPlanarRectangle.cpp index 6b6b7e8611..cfafc5ad35 100644 --- a/Modules/PlanarFigure/DataManagement/mitkPlanarRectangle.cpp +++ b/Modules/PlanarFigure/DataManagement/mitkPlanarRectangle.cpp @@ -1,150 +1,146 @@ /*=================================================================== 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 "mitkProperties.h" #include "mitkPlanarRectangle.h" #include "mitkGeometry2D.h" mitk::PlanarRectangle::PlanarRectangle() : FEATURE_ID_CIRCUMFERENCE( this->AddFeature( "Circumference", "mm" ) ), FEATURE_ID_AREA( this->AddFeature( "Area", "mm2" ) ) { // Rectangle has four control points this->ResetNumberOfControlPoints( 4 ); this->SetProperty( "closed", mitk::BoolProperty::New(true) ); this->SetNumberOfPolyLines( 1 ); } mitk::PlanarRectangle::~PlanarRectangle() { } bool mitk::PlanarRectangle::SetControlPoint( unsigned int index, const Point2D &point, bool createIfDoesNotExist ) { // heres the deal with the rectangle: // when a point is moved all corresponding corner points are moved with him // e.g. if the lower right point (index=3) is moved the upper right point (index=1) // is moved in the same x direction // and the lower left point (index=2) is moved in the same y direction // the upper left point (index=0) is left untouched bool set = PlanarFigure::SetControlPoint( index, point, createIfDoesNotExist ); if(set) { // can be made better ... unsigned int horizontalCorrespondingPointIndex = 1; unsigned int verticalCorrespondingPointIndex = 3; if(index == 1) { horizontalCorrespondingPointIndex = 0; verticalCorrespondingPointIndex = 2; } else if(index == 2) { horizontalCorrespondingPointIndex = 3; verticalCorrespondingPointIndex = 1; } else if(index == 3) { horizontalCorrespondingPointIndex = 2; verticalCorrespondingPointIndex = 0; } Point2D verticalCorrespondingPoint = GetControlPoint( verticalCorrespondingPointIndex ); verticalCorrespondingPoint[0] = point[0]; PlanarFigure::SetControlPoint( verticalCorrespondingPointIndex, verticalCorrespondingPoint ); Point2D horizontalCorrespondingPoint = GetControlPoint( horizontalCorrespondingPointIndex ); horizontalCorrespondingPoint[1] = point[1]; PlanarFigure::SetControlPoint( horizontalCorrespondingPointIndex, horizontalCorrespondingPoint ); } return set; } void mitk::PlanarRectangle::PlaceFigure( const mitk::Point2D &point ) { PlanarFigure::PlaceFigure( point ); m_SelectedControlPoint = 3; } void mitk::PlanarRectangle::GeneratePolyLine() { - // TODO: start polygon at specified initalize point... + this->ClearPolyLines(); - ClearPolyLines(); - - for ( unsigned int i = 0; i < this->GetNumberOfControlPoints(); ++i ) - { - AppendPointToPolyLine( 0, PolyLineElement( GetControlPoint(i), i ) ); - } + for (unsigned int i = 0; i < this->GetNumberOfControlPoints(); ++i) + this->AppendPointToPolyLine(0, this->GetControlPoint(i)); } void mitk::PlanarRectangle::GenerateHelperPolyLine(double /*mmPerDisplayUnit*/, unsigned int /*displayHeight*/) { // A polygon does not require helper objects } void mitk::PlanarRectangle::EvaluateFeaturesInternal() { // Calculate circumference double circumference = 0.0; unsigned int i; for ( i = 0; i < this->GetNumberOfControlPoints(); ++i ) { circumference += this->GetWorldControlPoint( i ).EuclideanDistanceTo( this->GetWorldControlPoint( (i + 1) % this->GetNumberOfControlPoints() ) ); } this->SetQuantity( FEATURE_ID_CIRCUMFERENCE, circumference ); // Calculate rectangle area (well, done a bit clumsy...) double area = 0.0; if ( this->GetGeometry2D() != NULL ) { for ( i = 0; i < this->GetNumberOfControlPoints(); ++i ) { Point2D p0 = this->GetControlPoint( i ); Point2D p1 = this->GetControlPoint( (i + 1) % this->GetNumberOfControlPoints() ); area += p0[0] * p1[1] - p1[0] * p0[1]; } area /= 2.0; } this->SetQuantity( FEATURE_ID_AREA, fabs(area) ); } void mitk::PlanarRectangle::PrintSelf( std::ostream& os, itk::Indent indent) const { Superclass::PrintSelf( os, indent ); os << indent << "Number of control points: " << this->GetNumberOfControlPoints() << std::endl; os << indent << "Control points:" << std::endl; for ( unsigned int i = 0; i < this->GetNumberOfControlPoints(); ++i ) { os << indent << indent << i << ": " < #include mitk::PlanarFigureReader::PlanarFigureReader() : PlanarFigureSource(), FileReader(), m_FileName(""), m_FilePrefix(""), m_FilePattern(""), m_Success(false) { this->SetNumberOfRequiredOutputs(1); this->SetNumberOfIndexedOutputs(1); this->SetNthOutput(0, this->MakeOutput(0)); m_CanReadFromMemory = true; //this->Modified(); //this->GetOutput()->Modified(); //this->GetOutput()->ReleaseData(); } mitk::PlanarFigureReader::~PlanarFigureReader() {} void mitk::PlanarFigureReader::GenerateData() { +MITK_INFO << "PlanarFigureReader::GenerateData()"; m_Success = false; this->SetNumberOfIndexedOutputs(0); // reset all outputs, we add new ones depending on the file content TiXmlDocument document; if(m_ReadFromMemory) { if(m_MemoryBuffer == NULL || m_MemorySize == 0) { //check itkWarningMacro( << "Sorry, memory buffer has not been set!" ); return; } if(m_MemoryBuffer[ m_MemorySize - 1 ] == '\0') { document.Parse(m_MemoryBuffer); } else { char * tmpArray = new char[(int)m_MemorySize+1]; tmpArray[m_MemorySize] = '\0'; memcpy(tmpArray,m_MemoryBuffer,m_MemorySize); document.Parse(m_MemoryBuffer); delete [] tmpArray; } } else { if (m_FileName.empty()) { itkWarningMacro( << "Sorry, filename has not been set!" ); return; } if (this->CanReadFile( m_FileName.c_str()) == false) { itkWarningMacro( << "Sorry, can't read file " << m_FileName << "!" ); return; } if (!document.LoadFile(m_FileName)) { MITK_ERROR << "Could not open/read/parse " << m_FileName << ". TinyXML reports: '" << document.ErrorDesc() << "'. " << "The error occurred in row " << document.ErrorRow() << ", column " << document.ErrorCol() << "."; return; } } int fileVersion = 1; TiXmlElement* versionObject = document.FirstChildElement("Version"); if (versionObject != NULL) { if ( versionObject->QueryIntAttribute( "FileVersion", &fileVersion ) != TIXML_SUCCESS ) { MITK_WARN << m_FileName << " does not contain version information! Trying version 1 format." << std::endl; } } else { MITK_WARN << m_FileName << " does not contain version information! Trying version 1 format." << std::endl; } if (fileVersion != 1) // add file version selection and version specific file parsing here, if newer file versions are created { MITK_WARN << "File version > 1 is not supported by this reader."; return; } /* file version 1 reader code */ for( TiXmlElement* pfElement = document.FirstChildElement("PlanarFigure"); pfElement != NULL; pfElement = pfElement->NextSiblingElement("PlanarFigure") ) { if (pfElement == NULL) continue; std::string type = pfElement->Attribute("type"); mitk::PlanarFigure::Pointer planarFigure = NULL; if (type == "PlanarAngle") { planarFigure = mitk::PlanarAngle::New(); } else if (type == "PlanarCircle") { planarFigure = mitk::PlanarCircle::New(); } else if (type == "PlanarEllipse") { planarFigure = mitk::PlanarEllipse::New(); } else if (type == "PlanarCross") { planarFigure = mitk::PlanarCross::New(); } else if (type == "PlanarFourPointAngle") { planarFigure = mitk::PlanarFourPointAngle::New(); } else if (type == "PlanarLine") { planarFigure = mitk::PlanarLine::New(); } else if (type == "PlanarPolygon") { planarFigure = mitk::PlanarPolygon::New(); } else if (type == "PlanarSubdivisionPolygon") { planarFigure = mitk::PlanarSubdivisionPolygon::New(); } else if (type == "PlanarRectangle") { planarFigure = mitk::PlanarRectangle::New(); } else if (type == "PlanarArrow") { planarFigure = mitk::PlanarArrow::New(); } else if (type == "PlanarDoubleEllipse") { planarFigure = mitk::PlanarDoubleEllipse::New(); } else if (type == "PlanarBezierCurve") { planarFigure = mitk::PlanarBezierCurve::New(); } else { // unknown type MITK_WARN << "encountered unknown planar figure type '" << type << "'. Skipping this element."; continue; } // Read properties of the planar figure for( TiXmlElement* propertyElement = pfElement->FirstChildElement("property"); propertyElement != NULL; propertyElement = propertyElement->NextSiblingElement("property") ) { const char* keya = propertyElement->Attribute("key"); std::string key( keya ? keya : ""); const char* typea = propertyElement->Attribute("type"); std::string type( typea ? typea : ""); // hand propertyElement to specific reader std::stringstream propertyDeserializerClassName; propertyDeserializerClassName << type << "Serializer"; std::list readers = itk::ObjectFactoryBase::CreateAllInstance(propertyDeserializerClassName.str().c_str()); if (readers.size() < 1) { MITK_ERROR << "No property reader found for " << type; } if (readers.size() > 1) { MITK_WARN << "Multiple property readers found for " << type << ". Using arbitrary first one."; } for ( std::list::iterator iter = readers.begin(); iter != readers.end(); ++iter ) { if (BasePropertySerializer* reader = dynamic_cast( iter->GetPointer() ) ) { BaseProperty::Pointer property = reader->Deserialize( propertyElement->FirstChildElement() ); if (property.IsNotNull()) { planarFigure->GetPropertyList()->ReplaceProperty(key, property); } else { MITK_ERROR << "There were errors while loading property '" << key << "' of type " << type << ". Your data may be corrupted"; } break; } } } // If we load a planarFigure, it has definitely been placed correctly. // If we do not set this property here, we cannot load old planarFigures // without messing up the interaction (PF-Interactor needs this property. planarFigure->GetPropertyList()->SetBoolProperty( "initiallyplaced", true ); // Read geometry of containing plane TiXmlElement* geoElement = pfElement->FirstChildElement("Geometry"); if (geoElement != NULL) { try { // Create plane geometry mitk::PlaneGeometry::Pointer planeGeo = mitk::PlaneGeometry::New(); // Extract and set plane transform parameters DoubleList transformList = this->GetDoubleAttributeListFromXMLNode( geoElement->FirstChildElement( "transformParam" ), "param", 12 ); typedef mitk::Geometry3D::TransformType TransformType; TransformType::ParametersType parameters; parameters.SetSize( 12 ); unsigned int i; DoubleList::iterator it; for ( it = transformList.begin(), i = 0; it != transformList.end(); ++it, ++i ) { parameters.SetElement( i, *it ); } typedef mitk::Geometry3D::TransformType TransformType; TransformType::Pointer affineGeometry = TransformType::New(); affineGeometry->SetParameters( parameters ); planeGeo->SetIndexToWorldTransform( affineGeometry ); // Extract and set plane bounds DoubleList boundsList = this->GetDoubleAttributeListFromXMLNode( geoElement->FirstChildElement( "boundsParam" ), "bound", 6 ); typedef mitk::Geometry3D::BoundsArrayType BoundsArrayType; BoundsArrayType bounds; for ( it = boundsList.begin(), i = 0; it != boundsList.end(); ++it, ++i ) { bounds[i] = *it; } planeGeo->SetBounds( bounds ); // Extract and set spacing and origin Vector3D spacing = this->GetVectorFromXMLNode(geoElement->FirstChildElement("Spacing")); planeGeo->SetSpacing( spacing ); Point3D origin = this->GetPointFromXMLNode(geoElement->FirstChildElement("Origin")); planeGeo->SetOrigin( origin ); planarFigure->SetGeometry2D(planeGeo); } catch (...) { } } TiXmlElement* cpElement = pfElement->FirstChildElement("ControlPoints"); bool first = true; if (cpElement != NULL) for( TiXmlElement* vertElement = cpElement->FirstChildElement("Vertex"); vertElement != NULL; vertElement = vertElement->NextSiblingElement("Vertex")) { if (vertElement == NULL) continue; int id = 0; mitk::Point2D::ValueType x = 0.0; mitk::Point2D::ValueType y = 0.0; if (vertElement->QueryIntAttribute("id", &id) == TIXML_WRONG_TYPE) return; // TODO: can we do a better error handling? if (vertElement->QueryDoubleAttribute("x", &x) == TIXML_WRONG_TYPE) return; // TODO: can we do a better error handling? if (vertElement->QueryDoubleAttribute("y", &y) == TIXML_WRONG_TYPE) return; // TODO: can we do a better error handling? Point2D p; p.SetElement(0, x); p.SetElement(1, y); +MITK_INFO << "Vertex " << id << ": (" << p[0] << ", " << p[1] << ")"; if (first == true) // needed to set m_FigurePlaced to true { planarFigure->PlaceFigure(p); first = false; } planarFigure->SetControlPoint(id, p, true); } // Calculate feature quantities of this PlanarFigure planarFigure->EvaluateFeatures(); // Make sure that no control point is currently selected planarFigure->DeselectControlPoint(); // \TODO: what about m_FigurePlaced and m_SelectedControlPoint ?? this->SetNthOutput( this->GetNumberOfOutputs(), planarFigure ); // add planarFigure as new output of this filter } m_Success = true; } mitk::Point3D mitk::PlanarFigureReader::GetPointFromXMLNode(TiXmlElement* e) { if (e == NULL) throw std::invalid_argument("node invalid"); // TODO: can we do a better error handling? mitk::Point3D point; mitk::ScalarType p(-1.0); if (e->QueryDoubleAttribute("x", &p) == TIXML_WRONG_TYPE) throw std::invalid_argument("node malformatted"); // TODO: can we do a better error handling? point.SetElement(0, p); if (e->QueryDoubleAttribute("y", &p) == TIXML_WRONG_TYPE) throw std::invalid_argument("node malformatted"); // TODO: can we do a better error handling? point.SetElement(1, p); if (e->QueryDoubleAttribute("z", &p) == TIXML_WRONG_TYPE) throw std::invalid_argument("node malformatted"); // TODO: can we do a better error handling? point.SetElement(2, p); return point; } mitk::Vector3D mitk::PlanarFigureReader::GetVectorFromXMLNode(TiXmlElement* e) { if (e == NULL) throw std::invalid_argument("node invalid"); // TODO: can we do a better error handling? mitk::Vector3D vector; mitk::ScalarType p(-1.0); if (e->QueryDoubleAttribute("x", &p) == TIXML_WRONG_TYPE) throw std::invalid_argument("node malformatted"); // TODO: can we do a better error handling? vector.SetElement(0, p); if (e->QueryDoubleAttribute("y", &p) == TIXML_WRONG_TYPE) throw std::invalid_argument("node malformatted"); // TODO: can we do a better error handling? vector.SetElement(1, p); if (e->QueryDoubleAttribute("z", &p) == TIXML_WRONG_TYPE) throw std::invalid_argument("node malformatted"); // TODO: can we do a better error handling? vector.SetElement(2, p); return vector; } mitk::PlanarFigureReader::DoubleList mitk::PlanarFigureReader::GetDoubleAttributeListFromXMLNode(TiXmlElement* e, const char *attributeNameBase, unsigned int count) { DoubleList list; if (e == NULL) throw std::invalid_argument("node invalid"); // TODO: can we do a better error handling? for ( unsigned int i = 0; i < count; ++i ) { mitk::ScalarType p(-1.0); std::stringstream attributeName; attributeName << attributeNameBase << i; if (e->QueryDoubleAttribute( attributeName.str().c_str(), &p ) == TIXML_WRONG_TYPE) throw std::invalid_argument("node malformatted"); // TODO: can we do a better error handling? list.push_back( p ); } return list; } void mitk::PlanarFigureReader::GenerateOutputInformation() { } int mitk::PlanarFigureReader::CanReadFile ( const char *name ) { if (std::string(name).empty()) return false; return (itksys::SystemTools::LowerCase(itksys::SystemTools::GetFilenameLastExtension(name)) == ".pf"); //assume, we can read all .pf files //TiXmlDocument document(name); //if (document.LoadFile() == false) // return false; //return (document.FirstChildElement("PlanarFigure") != NULL); } bool mitk::PlanarFigureReader::CanReadFile(const std::string filename, const std::string, const std::string) { if (filename.empty()) return false; return (itksys::SystemTools::LowerCase(itksys::SystemTools::GetFilenameLastExtension(filename)) == ".pf"); //assume, we can read all .pf files //TiXmlDocument document(filename); //if (document.LoadFile() == false) // return false; //return (document.FirstChildElement("PlanarFigure") != NULL); } void mitk::PlanarFigureReader::ResizeOutputs( const unsigned int& num ) { unsigned int prevNum = this->GetNumberOfOutputs(); this->SetNumberOfIndexedOutputs( num ); for ( unsigned int i = prevNum; i < num; ++i ) { this->SetNthOutput( i, this->MakeOutput( i ).GetPointer() ); } } diff --git a/Modules/PlanarFigure/Interactions/mitkPlanarFigureInteractor.cpp b/Modules/PlanarFigure/Interactions/mitkPlanarFigureInteractor.cpp index bc7ecae241..fd29413b95 100644 --- a/Modules/PlanarFigure/Interactions/mitkPlanarFigureInteractor.cpp +++ b/Modules/PlanarFigure/Interactions/mitkPlanarFigureInteractor.cpp @@ -1,947 +1,947 @@ /*=================================================================== 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. ===================================================================*/ #define PLANARFIGUREINTERACTOR_DBG MITK_DEBUG("PlanarFigureInteractor") << __LINE__ << ": " #include "mitkPlanarFigureInteractor.h" #include "mitkPlanarFigure.h" #include "mitkPlanarPolygon.h" #include "mitkInteractionPositionEvent.h" #include "mitkInternalEvent.h" #include "mitkBaseRenderer.h" #include "mitkRenderingManager.h" //how precise must the user pick the point //default value mitk::PlanarFigureInteractor::PlanarFigureInteractor() : DataInteractor() , m_Precision( 6.5 ) , m_MinimumPointDistance( 25.0 ) , m_IsHovering( false ) , m_LastPointWasValid( false ) { } mitk::PlanarFigureInteractor::~PlanarFigureInteractor() { } void mitk::PlanarFigureInteractor::ConnectActionsAndFunctions() { CONNECT_CONDITION("figure_is_on_current_slice", CheckFigureOnRenderingGeometry); CONNECT_CONDITION("figure_is_placed", CheckFigurePlaced); CONNECT_CONDITION("minimal_figure_is_finished", CheckMinimalFigureFinished); CONNECT_CONDITION("hovering_above_figure", CheckFigureHovering); CONNECT_CONDITION("hovering_above_point", CheckControlPointHovering); CONNECT_CONDITION("figure_is_selected", CheckSelection); CONNECT_CONDITION("point_is_valid", CheckPointValidity); CONNECT_CONDITION("figure_is_finished", CheckFigureFinished); CONNECT_CONDITION("reset_on_point_select_needed", CheckResetOnPointSelect); CONNECT_CONDITION("points_can_be_added_or_removed", CheckFigureIsExtendable); CONNECT_FUNCTION( "finalize_figure", FinalizeFigure); CONNECT_FUNCTION( "hide_preview_point", HidePreviewPoint ) CONNECT_FUNCTION( "hide_control_points", HideControlPoints ) CONNECT_FUNCTION( "set_preview_point_position", SetPreviewPointPosition ) CONNECT_FUNCTION( "move_current_point", MoveCurrentPoint); CONNECT_FUNCTION( "deselect_point", DeselectPoint); CONNECT_FUNCTION( "add_new_point", AddPoint); CONNECT_FUNCTION( "add_initial_point", AddInitialPoint); CONNECT_FUNCTION( "remove_selected_point", RemoveSelectedPoint); CONNECT_FUNCTION( "request_context_menu", RequestContextMenu); CONNECT_FUNCTION( "select_figure", SelectFigure ); CONNECT_FUNCTION( "select_point", SelectPoint ); CONNECT_FUNCTION( "end_interaction", EndInteraction ); CONNECT_FUNCTION( "start_hovering", StartHovering ) CONNECT_FUNCTION( "end_hovering", EndHovering ); } bool mitk::PlanarFigureInteractor::CheckFigurePlaced( const InteractionEvent* /*interactionEvent*/ ) { mitk::PlanarFigure *planarFigure = dynamic_cast( GetDataNode()->GetData() ); bool isFigureFinished = false; planarFigure->GetPropertyList()->GetBoolProperty( "initiallyplaced", isFigureFinished ); return planarFigure->IsPlaced() && isFigureFinished; } bool mitk::PlanarFigureInteractor::MoveCurrentPoint(StateMachineAction*, InteractionEvent* interactionEvent) { mitk::InteractionPositionEvent* positionEvent = dynamic_cast( interactionEvent ); if ( positionEvent == NULL ) return false; bool isEditable = true; GetDataNode()->GetBoolProperty( "planarfigure.iseditable", isEditable ); mitk::PlanarFigure *planarFigure = dynamic_cast( GetDataNode()->GetData() ); mitk::Geometry2D *planarFigureGeometry = dynamic_cast< Geometry2D * >( planarFigure->GetGeometry( 0 ) ); // Extract point in 2D world coordinates (relative to Geometry2D of // PlanarFigure) Point2D point2D; if ( !this->TransformPositionEventToPoint2D( positionEvent, planarFigureGeometry, point2D ) || !isEditable ) { return false; } // check if the control points shall be hidden during interaction bool hidecontrolpointsduringinteraction = false; GetDataNode()->GetBoolProperty( "planarfigure.hidecontrolpointsduringinteraction", hidecontrolpointsduringinteraction ); // hide the control points if necessary //interactionEvent->GetSender()->GetDataStorage()->BlockNodeModifiedEvents( true ); GetDataNode()->SetBoolProperty( "planarfigure.drawcontrolpoints", !hidecontrolpointsduringinteraction ); //interactionEvent->GetSender()->GetDataStorage()->BlockNodeModifiedEvents( false ); // Move current control point to this point planarFigure->SetCurrentControlPoint( point2D ); // Re-evaluate features planarFigure->EvaluateFeatures(); // Update rendered scene interactionEvent->GetSender()->GetRenderingManager()->RequestUpdateAll(); return true; } bool mitk::PlanarFigureInteractor::FinalizeFigure( StateMachineAction*, InteractionEvent* interactionEvent ) { mitk::PlanarFigure *planarFigure = dynamic_cast( GetDataNode()->GetData() ); planarFigure->Modified(); planarFigure->DeselectControlPoint(); planarFigure->RemoveLastControlPoint(); planarFigure->SetProperty( "initiallyplaced", mitk::BoolProperty::New( true ) ); GetDataNode()->SetBoolProperty( "planarfigure.drawcontrolpoints", true ); GetDataNode()->Modified(); planarFigure->InvokeEvent( EndPlacementPlanarFigureEvent() ); planarFigure->InvokeEvent( EndInteractionPlanarFigureEvent() ); interactionEvent->GetSender()->GetRenderingManager()->RequestUpdateAll(); return false; } bool mitk::PlanarFigureInteractor::EndInteraction( StateMachineAction*, InteractionEvent* interactionEvent ) { mitk::PlanarFigure *planarFigure = dynamic_cast( GetDataNode()->GetData() ); GetDataNode()->SetBoolProperty( "planarfigure.drawcontrolpoints", true ); planarFigure->Modified(); planarFigure->InvokeEvent( EndInteractionPlanarFigureEvent() ); interactionEvent->GetSender()->GetRenderingManager()->RequestUpdateAll(); return false; } bool mitk::PlanarFigureInteractor::EndHovering( StateMachineAction*, InteractionEvent* interactionEvent ) { mitk::PlanarFigure *planarFigure = dynamic_cast( GetDataNode()->GetData() ); planarFigure->ResetPreviewContolPoint(); // Invoke end-hover event once the mouse is exiting the figure area m_IsHovering = false; planarFigure->InvokeEvent( EndHoverPlanarFigureEvent() ); // Set bool property to indicate that planar figure is no longer in "hovering" mode GetDataNode()->SetBoolProperty( "planarfigure.ishovering", false ); interactionEvent->GetSender()->GetRenderingManager()->RequestUpdateAll(); return false; } bool mitk::PlanarFigureInteractor::CheckMinimalFigureFinished( const InteractionEvent* /*interactionEvent*/ ) { mitk::PlanarFigure *planarFigure = dynamic_cast( GetDataNode()->GetData() ); return ( planarFigure->GetNumberOfControlPoints() >= planarFigure->GetMinimumNumberOfControlPoints() ); } bool mitk::PlanarFigureInteractor::CheckFigureFinished( const InteractionEvent* /*interactionEvent*/ ) { mitk::PlanarFigure *planarFigure = dynamic_cast( GetDataNode()->GetData() ); return ( planarFigure->GetNumberOfControlPoints() >= planarFigure->GetMaximumNumberOfControlPoints() ); } bool mitk::PlanarFigureInteractor::CheckFigureIsExtendable( const InteractionEvent* /*interactionEvent*/ ) { bool isExtendable = false; GetDataNode()->GetBoolProperty("planarfigure.isextendable", isExtendable); return isExtendable; } bool mitk::PlanarFigureInteractor::DeselectPoint(StateMachineAction*, InteractionEvent* /*interactionEvent*/) { mitk::PlanarFigure *planarFigure = dynamic_cast( GetDataNode()->GetData() ); bool wasSelected = planarFigure->DeselectControlPoint(); if ( wasSelected ) { // Issue event so that listeners may update themselves planarFigure->Modified(); planarFigure->InvokeEvent( EndInteractionPlanarFigureEvent() ); GetDataNode()->SetBoolProperty( "planarfigure.drawcontrolpoints", true ); // GetDataNode()->SetBoolProperty( "planarfigure.ishovering", false ); GetDataNode()->Modified(); } return true; } bool mitk::PlanarFigureInteractor::AddPoint(StateMachineAction*, InteractionEvent* interactionEvent) { mitk::InteractionPositionEvent* positionEvent = dynamic_cast( interactionEvent ); if ( positionEvent == NULL ) return false; bool selected = false; bool isEditable = true; GetDataNode()->GetBoolProperty("selected", selected); GetDataNode()->GetBoolProperty( "planarfigure.iseditable", isEditable ); if ( !selected || !isEditable ) { return false; } mitk::PlanarFigure *planarFigure = dynamic_cast( GetDataNode()->GetData() ); mitk::Geometry2D *planarFigureGeometry = dynamic_cast< Geometry2D * >( planarFigure->GetGeometry( 0 ) ); // If the planarFigure already has reached the maximum number if ( planarFigure->GetNumberOfControlPoints() >= planarFigure->GetMaximumNumberOfControlPoints() ) { return false; } // Extract point in 2D world coordinates (relative to Geometry2D of // PlanarFigure) Point2D point2D, projectedPoint; if ( !this->TransformPositionEventToPoint2D( positionEvent, planarFigureGeometry, point2D ) ) { return false; } // TODO: check segment of polyline we clicked in int nextIndex = -1; // We only need to check which position to insert the control point // when interacting with a PlanarPolygon. For all other types // new control points will always be appended /* * Added check for "initiallyplaced" due to bug 13097: * * There are two possible cases in which a point can be inserted into a PlanarPolygon: * * 1. The figure is currently drawn -> the point will be appended at the end of the figure * 2. A point is inserted at a userdefined position after the initial placement of the figure is finished * * In the second case we need to determine the proper insertion index. In the first case the index always has * to be -1 so that the point is appended to the end. * * These changes are necessary because of a mac os x specific issue: If a users draws a PlanarPolygon then the * next point to be added moves according to the mouse position. If then the user left clicks in order to add * a point one would assume the last move position is identical to the left click position. This is actually the * case for windows and linux but somehow NOT for mac. Because of the insertion logic of a new point in the * PlanarFigure then for mac the wrong current selected point is determined. * * With this check here this problem can be avoided. However a redesign of the insertion logic should be considered */ bool isFigureFinished = false; planarFigure->GetPropertyList()->GetBoolProperty( "initiallyplaced", isFigureFinished ); mitk::BaseRenderer *renderer = interactionEvent->GetSender(); const Geometry2D *projectionPlane = renderer->GetCurrentWorldGeometry2D(); if ( dynamic_cast( planarFigure ) && isFigureFinished) { nextIndex = this->IsPositionOverFigure( positionEvent, planarFigure, planarFigureGeometry, projectionPlane, renderer->GetDisplayGeometry(), projectedPoint ); } // Add point as new control point renderer->GetDisplayGeometry()->DisplayToWorld( projectedPoint, projectedPoint ); if ( planarFigure->IsPreviewControlPointVisible() ) { point2D = planarFigure->GetPreviewControlPoint(); } planarFigure->AddControlPoint( point2D, nextIndex ); if ( planarFigure->IsPreviewControlPointVisible() ) { planarFigure->SelectControlPoint( nextIndex ); planarFigure->ResetPreviewContolPoint(); } // Re-evaluate features planarFigure->EvaluateFeatures(); //this->LogPrintPlanarFigureQuantities( planarFigure ); // Update rendered scene renderer->GetRenderingManager()->RequestUpdateAll(); return true; } bool mitk::PlanarFigureInteractor::AddInitialPoint(StateMachineAction*, InteractionEvent* interactionEvent) { mitk::InteractionPositionEvent* positionEvent = dynamic_cast( interactionEvent ); if ( positionEvent == NULL ) return false; mitk::PlanarFigure *planarFigure = dynamic_cast( GetDataNode()->GetData() ); mitk::BaseRenderer *renderer = interactionEvent->GetSender(); mitk::Geometry2D *planarFigureGeometry = dynamic_cast< Geometry2D * >( planarFigure->GetGeometry( 0 ) ); // Invoke event to notify listeners that placement of this PF starts now planarFigure->InvokeEvent( StartPlacementPlanarFigureEvent() ); // Use Geometry2D of the renderer clicked on for this PlanarFigure mitk::PlaneGeometry *planeGeometry = const_cast< mitk::PlaneGeometry * >( dynamic_cast< const mitk::PlaneGeometry * >( renderer->GetSliceNavigationController()->GetCurrentPlaneGeometry() ) ); if ( planeGeometry != NULL ) { planarFigureGeometry = planeGeometry; planarFigure->SetGeometry2D( planeGeometry ); } else { return false; } // Extract point in 2D world coordinates (relative to Geometry2D of // PlanarFigure) Point2D point2D; if ( !this->TransformPositionEventToPoint2D( positionEvent, planarFigureGeometry, point2D ) ) { return false; } // Place PlanarFigure at this point planarFigure->PlaceFigure( point2D ); // Re-evaluate features planarFigure->EvaluateFeatures(); //this->LogPrintPlanarFigureQuantities( planarFigure ); // Set a bool property indicating that the figure has been placed in // the current RenderWindow. This is required so that the same render // window can be re-aligned to the Geometry2D of the PlanarFigure later // on in an application. GetDataNode()->SetBoolProperty( "PlanarFigureInitializedWindow", true, renderer ); // Update rendered scene renderer->GetRenderingManager()->RequestUpdateAll(); return true; } bool mitk::PlanarFigureInteractor::StartHovering( StateMachineAction*, InteractionEvent* interactionEvent ) { mitk::InteractionPositionEvent* positionEvent = dynamic_cast( interactionEvent ); if ( positionEvent == NULL ) return false; mitk::PlanarFigure *planarFigure = dynamic_cast( GetDataNode()->GetData() ); mitk::BaseRenderer *renderer = interactionEvent->GetSender(); if ( !m_IsHovering ) { // Invoke hover event once when the mouse is entering the figure area m_IsHovering = true; planarFigure->InvokeEvent( StartHoverPlanarFigureEvent() ); // Set bool property to indicate that planar figure is currently in "hovering" mode GetDataNode()->SetBoolProperty( "planarfigure.ishovering", true ); renderer->GetRenderingManager()->RequestUpdateAll(); } return true; } bool mitk::PlanarFigureInteractor::SetPreviewPointPosition( StateMachineAction*, InteractionEvent* interactionEvent ) { mitk::InteractionPositionEvent* positionEvent = dynamic_cast( interactionEvent ); if ( positionEvent == NULL ) return false; mitk::PlanarFigure *planarFigure = dynamic_cast( GetDataNode()->GetData() ); mitk::BaseRenderer *renderer = interactionEvent->GetSender(); planarFigure->DeselectControlPoint(); mitk::Point2D pointProjectedOntoLine = positionEvent->GetPointerPositionOnScreen(); bool selected(false); bool isExtendable(false); bool isEditable(true); GetDataNode()->GetBoolProperty("selected", selected); GetDataNode()->GetBoolProperty("planarfigure.isextendable", isExtendable); GetDataNode()->GetBoolProperty("planarfigure.iseditable", isEditable ); if ( selected && isExtendable && isEditable ) { renderer->GetDisplayGeometry()->DisplayToWorld( pointProjectedOntoLine, pointProjectedOntoLine ); planarFigure->SetPreviewControlPoint( pointProjectedOntoLine ); } renderer->GetRenderingManager()->RequestUpdateAll(); return true; } bool mitk::PlanarFigureInteractor::HideControlPoints( StateMachineAction*, InteractionEvent* /*interactionEvent*/ ) { GetDataNode()->SetBoolProperty( "planarfigure.drawcontrolpoints", false ); return true; } bool mitk::PlanarFigureInteractor::HidePreviewPoint( StateMachineAction*, InteractionEvent* interactionEvent ) { mitk::PlanarFigure *planarFigure = dynamic_cast( GetDataNode()->GetData() ); planarFigure->ResetPreviewContolPoint(); mitk::BaseRenderer *renderer = interactionEvent->GetSender(); renderer->GetRenderingManager()->RequestUpdateAll(); return true; } bool mitk::PlanarFigureInteractor::CheckFigureHovering( const InteractionEvent* interactionEvent ) { const mitk::InteractionPositionEvent* positionEvent = dynamic_cast( interactionEvent ); if ( positionEvent == NULL ) return false; mitk::PlanarFigure *planarFigure = dynamic_cast( GetDataNode()->GetData() ); mitk::BaseRenderer *renderer = interactionEvent->GetSender(); mitk::Geometry2D *planarFigureGeometry = dynamic_cast< Geometry2D * >( planarFigure->GetGeometry( 0 ) ); const Geometry2D *projectionPlane = renderer->GetCurrentWorldGeometry2D(); mitk::Point2D pointProjectedOntoLine; int previousControlPoint = this->IsPositionOverFigure( positionEvent, planarFigure, planarFigureGeometry, projectionPlane, renderer->GetDisplayGeometry(), pointProjectedOntoLine ); bool isHovering = (previousControlPoint != -1); if ( isHovering ) { return true; } else { return false; } return false; } bool mitk::PlanarFigureInteractor::CheckControlPointHovering( const InteractionEvent* interactionEvent ) { const mitk::InteractionPositionEvent* positionEvent = dynamic_cast( interactionEvent ); if ( positionEvent == NULL ) return false; mitk::PlanarFigure *planarFigure = dynamic_cast( GetDataNode()->GetData() ); mitk::BaseRenderer *renderer = interactionEvent->GetSender(); mitk::Geometry2D *planarFigureGeometry = dynamic_cast< Geometry2D * >( planarFigure->GetGeometry( 0 ) ); const Geometry2D *projectionPlane = renderer->GetCurrentWorldGeometry2D(); int pointIndex = -1; pointIndex = mitk::PlanarFigureInteractor::IsPositionInsideMarker( positionEvent, planarFigure, planarFigureGeometry, projectionPlane, renderer->GetDisplayGeometry() ); if ( pointIndex >= 0 ) { return true; } else { return false; } } bool mitk::PlanarFigureInteractor::CheckSelection( const InteractionEvent* /*interactionEvent*/ ) { bool selected = false; GetDataNode()->GetBoolProperty("selected", selected); return selected; } bool mitk::PlanarFigureInteractor::SelectFigure( StateMachineAction*, InteractionEvent* /*interactionEvent*/ ) { mitk::PlanarFigure *planarFigure = dynamic_cast( GetDataNode()->GetData() ); planarFigure->InvokeEvent( SelectPlanarFigureEvent() ); return false; } bool mitk::PlanarFigureInteractor::SelectPoint( StateMachineAction*, InteractionEvent* interactionEvent ) { mitk::InteractionPositionEvent* positionEvent = dynamic_cast( interactionEvent ); if ( positionEvent == NULL ) return false; mitk::PlanarFigure *planarFigure = dynamic_cast( GetDataNode()->GetData() ); mitk::BaseRenderer *renderer = interactionEvent->GetSender(); mitk::Geometry2D *planarFigureGeometry = dynamic_cast< Geometry2D * >( planarFigure->GetGeometry( 0 ) ); const Geometry2D *projectionPlane = renderer->GetCurrentWorldGeometry2D(); int pointIndex = -1; pointIndex = mitk::PlanarFigureInteractor::IsPositionInsideMarker( positionEvent, planarFigure, planarFigureGeometry, projectionPlane, renderer->GetDisplayGeometry() ); if ( pointIndex >= 0 ) { // If mouse is above control point, mark it as selected planarFigure->SelectControlPoint( pointIndex ); } else { planarFigure->DeselectControlPoint(); } return false; } bool mitk::PlanarFigureInteractor::CheckPointValidity( const InteractionEvent* interactionEvent ) { // Check if the distance of the current point to the previously set point in display coordinates // is sufficient (if a previous point exists) // Extract display position const mitk::InteractionPositionEvent* positionEvent = dynamic_cast( interactionEvent ); if ( positionEvent == NULL ) return false; mitk::PlanarFigure *planarFigure = dynamic_cast( GetDataNode()->GetData() ); m_LastPointWasValid = IsMousePositionAcceptableAsNewControlPoint( positionEvent, planarFigure ); return m_LastPointWasValid; } bool mitk::PlanarFigureInteractor::RemoveSelectedPoint(StateMachineAction*, InteractionEvent* interactionEvent) { mitk::PlanarFigure *planarFigure = dynamic_cast( GetDataNode()->GetData() ); mitk::BaseRenderer *renderer = interactionEvent->GetSender(); int selectedControlPoint = planarFigure->GetSelectedControlPoint(); planarFigure->RemoveControlPoint( selectedControlPoint ); // Re-evaluate features planarFigure->EvaluateFeatures(); planarFigure->Modified(); GetDataNode()->SetBoolProperty( "planarfigure.drawcontrolpoints", true ); planarFigure->InvokeEvent( EndInteractionPlanarFigureEvent() ); renderer->GetRenderingManager()->RequestUpdateAll(); HandleEvent( mitk::InternalEvent::New( renderer, this, "Dummy-Event" ), GetDataNode() ); return true; } bool mitk::PlanarFigureInteractor::RequestContextMenu(StateMachineAction*, InteractionEvent* /*interactionEvent*/) { mitk::PlanarFigure *planarFigure = dynamic_cast( GetDataNode()->GetData() ); bool selected = false; GetDataNode()->GetBoolProperty("selected", selected); // no need to invoke this if the figure is already selected if ( !selected ) { planarFigure->InvokeEvent( SelectPlanarFigureEvent() ); } planarFigure->InvokeEvent( ContextMenuPlanarFigureEvent() ); return true; } bool mitk::PlanarFigureInteractor::CheckResetOnPointSelect( const InteractionEvent* /*interactionEvent*/ ) { mitk::PlanarFigure *planarFigure = dynamic_cast( GetDataNode()->GetData() ); // Invoke tmpEvent to notify listeners that interaction with this PF starts now planarFigure->InvokeEvent( StartInteractionPlanarFigureEvent() ); // Reset the PlanarFigure if required return planarFigure->ResetOnPointSelect(); } bool mitk::PlanarFigureInteractor::CheckFigureOnRenderingGeometry( const InteractionEvent* interactionEvent ) { const mitk::InteractionPositionEvent* posEvent = dynamic_cast(interactionEvent); if ( posEvent == NULL ) return false; mitk::Point3D worldPoint3D = posEvent->GetPositionInWorld(); mitk::PlanarFigure *planarFigure = dynamic_cast( GetDataNode()->GetData() ); mitk::Geometry2D *planarFigureGeometry2D = dynamic_cast< Geometry2D * >( planarFigure->GetGeometry( 0 ) ); double planeThickness = planarFigureGeometry2D->GetExtentInMM( 2 ); if ( planarFigureGeometry2D->Distance( worldPoint3D ) > planeThickness ) { // don't react, when interaction is too far away return false; } return true; } void mitk::PlanarFigureInteractor::SetPrecision( mitk::ScalarType precision ) { m_Precision = precision; } void mitk::PlanarFigureInteractor::SetMinimumPointDistance( ScalarType minimumDistance ) { m_MinimumPointDistance = minimumDistance; } bool mitk::PlanarFigureInteractor::TransformPositionEventToPoint2D( const InteractionPositionEvent *positionEvent, const Geometry2D *planarFigureGeometry, Point2D &point2D ) { mitk::Point3D worldPoint3D = positionEvent->GetPositionInWorld(); // TODO: proper handling of distance tolerance if ( planarFigureGeometry->Distance( worldPoint3D ) > 0.1 ) { return false; } // Project point onto plane of this PlanarFigure planarFigureGeometry->Map( worldPoint3D, point2D ); return true; } bool mitk::PlanarFigureInteractor::TransformObjectToDisplay( const mitk::Point2D &point2D, mitk::Point2D &displayPoint, const mitk::Geometry2D *objectGeometry, const mitk::Geometry2D *rendererGeometry, const mitk::DisplayGeometry *displayGeometry ) const { mitk::Point3D point3D; // Map circle point from local 2D geometry into 3D world space objectGeometry->Map( point2D, point3D ); // TODO: proper handling of distance tolerance if ( displayGeometry->Distance( point3D ) < 0.1 ) { // Project 3D world point onto display geometry rendererGeometry->Map( point3D, displayPoint ); displayGeometry->WorldToDisplay( displayPoint, displayPoint ); return true; } return false; } bool mitk::PlanarFigureInteractor::IsPointNearLine( const mitk::Point2D& point, const mitk::Point2D& startPoint, const mitk::Point2D& endPoint, mitk::Point2D& projectedPoint ) const { mitk::Vector2D n1 = endPoint - startPoint; n1.Normalize(); // Determine dot products between line vector and startpoint-point / endpoint-point vectors double l1 = n1 * (point - startPoint); double l2 = -n1 * (point - endPoint); // Determine projection of specified point onto line defined by start / end point mitk::Point2D crossPoint = startPoint + n1 * l1; projectedPoint = crossPoint; // Point is inside encompassing rectangle IF // - its distance to its projected point is small enough // - it is not further outside of the line than the defined tolerance if (((crossPoint.SquaredEuclideanDistanceTo(point) < 20.0) && (l1 > 0.0) && (l2 > 0.0)) || endPoint.SquaredEuclideanDistanceTo(point) < 20.0 || startPoint.SquaredEuclideanDistanceTo(point) < 20.0) { return true; } return false; } int mitk::PlanarFigureInteractor::IsPositionOverFigure( const InteractionPositionEvent *positionEvent, PlanarFigure *planarFigure, const Geometry2D *planarFigureGeometry, const Geometry2D *rendererGeometry, const DisplayGeometry *displayGeometry, Point2D& pointProjectedOntoLine ) const { mitk::Point2D displayPosition = positionEvent->GetPointerPositionOnScreen(); // Iterate over all polylines of planar figure, and check if // any one is close to the current display position typedef mitk::PlanarFigure::PolyLineType VertexContainerType; Point2D polyLinePoint; Point2D firstPolyLinePoint; Point2D previousPolyLinePoint; for ( unsigned short loop=0; loopGetPolyLinesSize(); ++loop ) { const VertexContainerType polyLine = planarFigure->GetPolyLine( loop ); bool firstPoint( true ); for ( VertexContainerType::const_iterator it = polyLine.begin(); it != polyLine.end(); ++it ) { // Get plane coordinates of this point of polyline (if possible) if ( !this->TransformObjectToDisplay( it->Point, polyLinePoint, planarFigureGeometry, rendererGeometry, displayGeometry ) ) { break; // Poly line invalid (not on current 2D plane) --> skip it } if ( firstPoint ) { firstPolyLinePoint = polyLinePoint; firstPoint = false; } else if ( this->IsPointNearLine( displayPosition, previousPolyLinePoint, polyLinePoint, pointProjectedOntoLine ) ) { // Point is close enough to line segment --> Return index of the segment - return it->Index; + return std::distance(polyLine.begin(), it); } previousPolyLinePoint = polyLinePoint; } // For closed figures, also check last line segment if ( planarFigure->IsClosed() && this->IsPointNearLine( displayPosition, polyLinePoint, firstPolyLinePoint, pointProjectedOntoLine ) ) { return 0; // Return index of first control point } } return -1; } int mitk::PlanarFigureInteractor::IsPositionInsideMarker( const InteractionPositionEvent* positionEvent, const PlanarFigure *planarFigure, const Geometry2D *planarFigureGeometry, const Geometry2D *rendererGeometry, const DisplayGeometry *displayGeometry ) const { mitk::Point2D displayPosition = positionEvent->GetPointerPositionOnScreen(); // Iterate over all control points of planar figure, and check if // any one is close to the current display position mitk::Point2D displayControlPoint; int numberOfControlPoints = planarFigure->GetNumberOfControlPoints(); for ( int i=0; iTransformObjectToDisplay( planarFigure->GetControlPoint(i), displayControlPoint, planarFigureGeometry, rendererGeometry, displayGeometry ) ) { // TODO: variable size of markers if ( displayPosition.SquaredEuclideanDistanceTo( displayControlPoint ) < 20.0 ) { return i; } } } return -1; } void mitk::PlanarFigureInteractor::LogPrintPlanarFigureQuantities( const PlanarFigure *planarFigure ) { MITK_INFO << "PlanarFigure: " << planarFigure->GetNameOfClass(); for ( unsigned int i = 0; i < planarFigure->GetNumberOfFeatures(); ++i ) { MITK_INFO << "* " << planarFigure->GetFeatureName( i ) << ": " << planarFigure->GetQuantity( i ) << " " << planarFigure->GetFeatureUnit( i ); } } bool mitk::PlanarFigureInteractor::IsMousePositionAcceptableAsNewControlPoint( const mitk::InteractionPositionEvent* positionEvent, const PlanarFigure* planarFigure ) { assert(positionEvent && planarFigure); BaseRenderer* renderer = positionEvent->GetSender(); assert(renderer); // Get the timestep to support 3D+t int timeStep( renderer->GetTimeStep( planarFigure ) ); bool tooClose(false); const Geometry2D *renderingPlane = renderer->GetCurrentWorldGeometry2D(); mitk::Geometry2D *planarFigureGeometry = dynamic_cast< mitk::Geometry2D * >( planarFigure->GetGeometry( timeStep ) ); Point2D point2D, correctedPoint; // Get the point2D from the positionEvent if ( !this->TransformPositionEventToPoint2D( positionEvent, planarFigureGeometry, point2D ) ) { return false; } // apply the controlPoint constraints of the planarFigure to get the // coordinates that would actually be used. correctedPoint = const_cast( planarFigure )->ApplyControlPointConstraints( 0, point2D ); // map the 2D coordinates of the new point to world-coordinates // and transform those to display-coordinates mitk::Point3D newPoint3D; planarFigureGeometry->Map( correctedPoint, newPoint3D ); mitk::Point2D newDisplayPosition; renderingPlane->Map( newPoint3D, newDisplayPosition ); renderer->GetDisplayGeometry()->WorldToDisplay( newDisplayPosition, newDisplayPosition ); for( int i=0; i < (int)planarFigure->GetNumberOfControlPoints(); i++ ) { if ( i != planarFigure->GetSelectedControlPoint() ) { // Try to convert previous point to current display coordinates mitk::Point3D previousPoint3D; // map the 2D coordinates of the control-point to world-coordinates planarFigureGeometry->Map( planarFigure->GetControlPoint( i ), previousPoint3D ); if ( renderer->GetDisplayGeometry()->Distance( previousPoint3D ) < 0.1 ) // ugly, but assert makes this work { mitk::Point2D previousDisplayPosition; // transform the world-coordinates into display-coordinates renderingPlane->Map( previousPoint3D, previousDisplayPosition ); renderer->GetDisplayGeometry()->WorldToDisplay( previousDisplayPosition, previousDisplayPosition ); //Calculate the distance. We use display-coordinates here to make // the check independent of the zoom-level of the rendering scene. double a = newDisplayPosition[0] - previousDisplayPosition[0]; double b = newDisplayPosition[1] - previousDisplayPosition[1]; // If point is to close, do not set a new point tooClose = (a * a + b * b < m_MinimumPointDistance ); } if ( tooClose ) return false; // abort loop early } } return !tooClose; // default } void mitk::PlanarFigureInteractor::ConfigurationChanged() { mitk::PropertyList::Pointer properties = GetAttributes(); std::string precision = ""; if (properties->GetStringProperty("precision", precision)) { m_Precision = atof(precision.c_str()); } else { m_Precision = (ScalarType) 6.5; } std::string minPointDistance = ""; if (properties->GetStringProperty("minPointDistance", minPointDistance)) { m_MinimumPointDistance = atof(minPointDistance.c_str()); } else { m_MinimumPointDistance = (ScalarType) 25.0; } } diff --git a/Modules/PlanarFigureSegmentation/mitkPlanarFigureSegmentationController.cpp b/Modules/PlanarFigureSegmentation/mitkPlanarFigureSegmentationController.cpp index 64bb34165d..f242d24454 100644 --- a/Modules/PlanarFigureSegmentation/mitkPlanarFigureSegmentationController.cpp +++ b/Modules/PlanarFigureSegmentation/mitkPlanarFigureSegmentationController.cpp @@ -1,313 +1,312 @@ /*=================================================================== 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 "mitkPlanarFigureSegmentationController.h" #include "mitkSurfaceToImageFilter.h" #include #include #include #include #include #include #include "mitkImageWriter.h" #include "mitkSurfaceVtkWriter.h" #include "mitkImageToSurfaceFilter.h" #include "mitkImageAccessByItk.h" #include "mitkImageCast.h" mitk::PlanarFigureSegmentationController::PlanarFigureSegmentationController() : itk::Object() , m_ReduceFilter( NULL ) , m_NormalsFilter( NULL ) , m_DistanceImageCreator( NULL ) , m_ReferenceImage( NULL ) , m_SegmentationAsImage( NULL ) { InitializeFilters(); } mitk::PlanarFigureSegmentationController::~PlanarFigureSegmentationController() { } void mitk::PlanarFigureSegmentationController::SetReferenceImage( mitk::Image::Pointer referenceImage ) { m_ReferenceImage = referenceImage; } void mitk::PlanarFigureSegmentationController::AddPlanarFigure( mitk::PlanarFigure::Pointer planarFigure ) { if ( planarFigure.IsNull() ) return; bool newFigure = true; std::size_t indexOfFigure = 0; for( std::size_t i=0; iCreateSurfaceFromPlanarFigure( planarFigure ); m_SurfaceList.push_back( figureAsSurface ); if (!m_PlanarFigureList.empty()) { indexOfFigure = m_PlanarFigureList.size() -1 ; } } else { figureAsSurface = this->CreateSurfaceFromPlanarFigure( planarFigure ); m_SurfaceList.at(indexOfFigure) = figureAsSurface; } if ( m_ReduceFilter.IsNull() ) { InitializeFilters(); } m_ReduceFilter->SetInput( indexOfFigure, figureAsSurface ); m_NormalsFilter->SetInput( indexOfFigure, m_ReduceFilter->GetOutput( indexOfFigure ) ); m_DistanceImageCreator->SetInput( indexOfFigure, m_NormalsFilter->GetOutput( indexOfFigure ) ); } void mitk::PlanarFigureSegmentationController::RemovePlanarFigure( mitk::PlanarFigure::Pointer planarFigure ) { if ( planarFigure.IsNull() ) return; bool figureFound = false; std::size_t indexOfFigure = 0; for( std::size_t i=0; iRemoveInputs( m_NormalsFilter->GetOutput( indexOfFigure ) ); // m_NormalsFilter->RemoveInput( m_ReduceFilter->GetOutput( indexOfFigure ) ); // m_ReduceFilter->RemoveInput( const_cast(m_ReduceFilter->GetInput(indexOfFigure)) ); } else { // this is not very nice! If the figure that has been removed is NOT the last // one in the list we have to create new filters and add all remaining // inputs again. // // Has to be done as the filters do not work when removing an input // other than the last one. // create new filters InitializeFilters(); // and add all existing surfaces SurfaceListType::iterator surfaceIter = m_SurfaceList.begin(); for ( surfaceIter = m_SurfaceList.begin(); surfaceIter!=m_SurfaceList.end(); surfaceIter++ ) { m_ReduceFilter->SetInput( indexOfFigure, (*surfaceIter) ); m_NormalsFilter->SetInput( indexOfFigure, m_ReduceFilter->GetOutput( indexOfFigure ) ); m_DistanceImageCreator->SetInput( indexOfFigure, m_NormalsFilter->GetOutput( indexOfFigure ) ); } } PlanarFigureListType::iterator whereIter = m_PlanarFigureList.begin(); whereIter += indexOfFigure; m_PlanarFigureList.erase( whereIter ); SurfaceListType::iterator surfaceIter = m_SurfaceList.begin(); surfaceIter += indexOfFigure; m_SurfaceList.erase( surfaceIter ); } template void mitk::PlanarFigureSegmentationController::GetImageBase(itk::Image* input, itk::ImageBase<3>::Pointer& result) { result = input; } mitk::Image::Pointer mitk::PlanarFigureSegmentationController::GetInterpolationResult() { m_SegmentationAsImage = NULL; if ( m_PlanarFigureList.size() == 0 ) { m_SegmentationAsImage = mitk::Image::New(); m_SegmentationAsImage->Initialize(mitk::MakeScalarPixelType() , *m_ReferenceImage->GetTimeGeometry()); return m_SegmentationAsImage; } itk::ImageBase<3>::Pointer itkImage; AccessFixedDimensionByItk_1( m_ReferenceImage.GetPointer(), GetImageBase, 3, itkImage ); m_DistanceImageCreator->SetReferenceImage( itkImage.GetPointer() ); m_ReduceFilter->Update(); m_NormalsFilter->Update(); m_DistanceImageCreator->Update(); mitk::Image::Pointer distanceImage = m_DistanceImageCreator->GetOutput(); // Cleanup the pipeline distanceImage->DisconnectPipeline(); m_DistanceImageCreator = NULL; m_NormalsFilter = NULL; m_ReduceFilter = NULL; itkImage = NULL; // If this bool flag is true, the distanceImage will be written to the // filesystem as nrrd-image and as surface-representation. bool debugOutput(false); if ( debugOutput ) { mitk::ImageWriter::Pointer imageWriter = mitk::ImageWriter::New(); imageWriter->SetInput( distanceImage ); imageWriter->SetExtension( ".nrrd" ); imageWriter->SetFileName( "v:/DistanceImage" ); imageWriter->Update(); } mitk::ImageToSurfaceFilter::Pointer imageToSurfaceFilter = mitk::ImageToSurfaceFilter::New(); imageToSurfaceFilter->SetInput( distanceImage ); imageToSurfaceFilter->SetThreshold( 0 ); imageToSurfaceFilter->Update(); mitk::Surface::Pointer segmentationAsSurface = imageToSurfaceFilter->GetOutput(); // Cleanup the pipeline segmentationAsSurface->DisconnectPipeline(); imageToSurfaceFilter = NULL; if ( debugOutput ) { mitk::SurfaceVtkWriter::Pointer surfaceWriter = mitk::SurfaceVtkWriter::New(); surfaceWriter->SetInput( segmentationAsSurface ); surfaceWriter->SetExtension( ".vtk" ); surfaceWriter->SetFileName( "v:/DistanceImageAsSurface.vtk" ); surfaceWriter->Update(); } mitk::SurfaceToImageFilter::Pointer surfaceToImageFilter = mitk::SurfaceToImageFilter::New(); surfaceToImageFilter->SetInput( segmentationAsSurface ); surfaceToImageFilter->SetImage( m_ReferenceImage ); surfaceToImageFilter->SetMakeOutputBinary(true); surfaceToImageFilter->Update(); m_SegmentationAsImage = surfaceToImageFilter->GetOutput(); // Cleanup the pipeline m_SegmentationAsImage->DisconnectPipeline(); return m_SegmentationAsImage; } mitk::Surface::Pointer mitk::PlanarFigureSegmentationController::CreateSurfaceFromPlanarFigure( mitk::PlanarFigure::Pointer figure ) { if ( figure.IsNull() ) { MITK_ERROR << "Given PlanarFigure is NULL. Please provide valid PlanarFigure."; return NULL; } mitk::Surface::Pointer newSurface = mitk::Surface::New(); vtkSmartPointer points = vtkSmartPointer::New(); vtkSmartPointer polygon = vtkSmartPointer::New(); vtkSmartPointer cells = vtkSmartPointer::New(); vtkSmartPointer polyData = vtkSmartPointer::New(); const mitk::Geometry2D* figureGeometry = figure->GetGeometry2D(); // Get the polyline mitk::PlanarFigure::PolyLineType planarPolyLine = figure->GetPolyLine(0); mitk::PlanarFigure::PolyLineType::iterator iter; // iterate over the polyline, ... int pointCounter = 0; for( iter = planarPolyLine.begin(); iter != planarPolyLine.end(); iter++ ) { // ... determine the world-coordinates - mitk::Point2D polyLinePoint = iter->Point; mitk::Point3D pointInWorldCoordiantes; - figureGeometry->Map( polyLinePoint, pointInWorldCoordiantes ); + figureGeometry->Map( *iter, pointInWorldCoordiantes ); // and add them as new points to the vtkPoints points->InsertNextPoint( pointInWorldCoordiantes[0], pointInWorldCoordiantes[1], pointInWorldCoordiantes[2] ); ++pointCounter; } // create a polygon with the points of the polyline polygon->GetPointIds()->SetNumberOfIds( pointCounter ); for(int i = 0; i < pointCounter; i++) { polygon->GetPointIds()->SetId(i,i); } // initialize the vtkCellArray and vtkPolyData cells->InsertNextCell(polygon); polyData->SetPoints(points); polyData->SetPolys( cells ); // set the polydata to the surface newSurface->SetVtkPolyData( polyData ); return newSurface; } mitk::PlanarFigureSegmentationController::PlanarFigureListType mitk::PlanarFigureSegmentationController::GetAllPlanarFigures() { return m_PlanarFigureList; } void mitk::PlanarFigureSegmentationController::InitializeFilters() { m_ReduceFilter = mitk::ReduceContourSetFilter::New(); m_ReduceFilter->SetReductionType(ReduceContourSetFilter::NTH_POINT); m_ReduceFilter->SetStepSize( 10 ); m_NormalsFilter = mitk::ComputeContourSetNormalsFilter::New(); m_DistanceImageCreator = mitk::CreateDistanceImageFromSurfaceFilter::New(); }