diff --git a/Modules/MitkExt/Algorithms/mitkImageStatisticsCalculator.cpp b/Modules/MitkExt/Algorithms/mitkImageStatisticsCalculator.cpp index f1ab8103ce..9a33945c9c 100644 --- a/Modules/MitkExt/Algorithms/mitkImageStatisticsCalculator.cpp +++ b/Modules/MitkExt/Algorithms/mitkImageStatisticsCalculator.cpp @@ -1,1068 +1,1028 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date: 2009-05-12 19:56:03 +0200 (Di, 12 Mai 2009) $ Version: $Revision: 17179 $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #include "mitkImageStatisticsCalculator.h" #include "mitkImageAccessByItk.h" #include "mitkExtractImageFilter.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include +#if ( ( VTK_MAJOR_VERSION <= 5 ) && ( VTK_MINOR_VERSION<=8) ) + #include "mitkvtkLassoStencilSource.h" +#else + #include "vtkLassoStencilSource.h" +#endif + + #include #include - namespace mitk { ImageStatisticsCalculator::ImageStatisticsCalculator() : m_MaskingMode( MASKING_MODE_NONE ), m_MaskingModeChanged( false ), m_DoIgnorePixelValue(false), m_IgnorePixelValue(0.0), m_IgnorePixelValueChanged(false) { m_EmptyHistogram = HistogramType::New(); HistogramType::SizeType histogramSize; 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_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(); for ( unsigned int t = 0; t < m_Image->GetTimeSteps(); ++t ) { m_MaskedImageStatisticsTimeStampVector[t].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[t] = true; } } } void ImageStatisticsCalculator::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(); 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() ) { itkExceptionMacro( << "Image not set!" ); } 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 ); Statistics *statistics; HistogramType::ConstPointer *histogram; switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(!m_DoIgnorePixelValue) { statistics = &m_ImageStatisticsVector[timeStep]; histogram = &m_ImageHistogramVector[timeStep]; m_ImageStatisticsTimeStampVector[timeStep].Modified(); m_ImageStatisticsCalculationTriggerVector[timeStep] = false; } else { statistics = &m_MaskedImageStatisticsVector[timeStep]; histogram = &m_MaskedImageHistogramVector[timeStep]; m_MaskedImageStatisticsTimeStampVector[timeStep].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false; } break; } case MASKING_MODE_IMAGE: statistics = &m_MaskedImageStatisticsVector[timeStep]; histogram = &m_MaskedImageHistogramVector[timeStep]; m_MaskedImageStatisticsTimeStampVector[timeStep].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false; break; case MASKING_MODE_PLANARFIGURE: statistics = &m_PlanarFigureStatisticsVector[timeStep]; histogram = &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, *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 && !m_DoIgnorePixelValue ) { 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::ConstPointer(); m_InternalImageMask3D = MaskImage3DType::Pointer(); m_InternalImageMask2D = MaskImage2DType::Pointer(); return true; } const ImageStatisticsCalculator::HistogramType * ImageStatisticsCalculator::GetHistogram( unsigned int timeStep ) 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]; 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 ) 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]; 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) { CastToItkImage( timeSliceImage, m_InternalImageMask3D ); m_InternalImageMask3D->FillBuffer(1); } break; } case MASKING_MODE_IMAGE: { if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() > 1) ) { if ( timeStep < m_ImageMask->GetTimeSteps() ) { 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 has not enough time steps!" ); } } 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!" ); } // Find slice number corresponding to PlanarFigure in input image MaskImage3DType::IndexType index; imageGeometry->WorldToIndex( planarFigureGeometry->GetOrigin(), index ); unsigned int slice = index[axis]; // Extract slice with given position and direction from image ExtractImageFilter::Pointer imageExtractor = ExtractImageFilter::New(); imageExtractor->SetInput( timeSliceImage ); imageExtractor->SetSliceDimension( axis ); imageExtractor->SetSliceIndex( slice ); imageExtractor->Update(); m_InternalImage = imageExtractor->GetOutput(); // 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, Statistics &statistics, typename HistogramType::ConstPointer *histogram ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef typename ImageType::IndexType IndexType; typedef itk::Statistics::ScalarImageToHistogramGenerator< ImageType > HistogramGeneratorType; // 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 histogram typename HistogramGeneratorType::Pointer histogramGenerator = HistogramGeneratorType::New(); histogramGenerator->SetInput( image ); histogramGenerator->SetMarginalScale( 100 ); // Defines y-margin width of histogram histogramGenerator->SetNumberOfBins( 768 ); // CT range [-1024, +2048] --> bin size 4 values histogramGenerator->SetHistogramMin( -1024.0 ); histogramGenerator->SetHistogramMax( 2048.0 ); histogramGenerator->Compute(); *histogram = histogramGenerator->GetOutput(); // 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 ); } 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 = itmask.Begin(); itimage = itimage.Begin(); 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, Statistics &statistics, typename HistogramType::ConstPointer *histogram ) { 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; // 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 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]; for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) { double indexCoordDistance = (maskOrigin[i] - imageOrigin[i]) / imageSpacing[i]; double misalignment = indexCoordDistance - floor( indexCoordDistance + 0.5 ); if ( fabs( misalignment ) > imageSpacing[i] / 20.0 ) { itkExceptionMacro( << "Pixels/voxels of mask and image are not sufficiently aligned! (Misalignment: " << misalignment << ")" ); } 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 typename LabelStatisticsFilterType::Pointer labelStatisticsFilter; labelStatisticsFilter = LabelStatisticsFilterType::New(); labelStatisticsFilter->SetInput( adaptedImage ); labelStatisticsFilter->SetLabelInput( adaptedMaskImage ); labelStatisticsFilter->UseHistogramsOn(); labelStatisticsFilter->SetHistogramParameters( 384, -1024.0, 2048.0); // 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 label of mask (other than 0) bool maskNonEmpty = false; unsigned int i; for ( i = 1; i < 4096; ++i ) { if ( labelStatisticsFilter->HasLabel( i ) ) { maskNonEmpty = true; break; } } if ( maskNonEmpty ) { *histogram = labelStatisticsFilter->GetHistogram( i ); statistics.N = labelStatisticsFilter->GetCount( i ); statistics.Min = labelStatisticsFilter->GetMinimum( i ); statistics.Max = labelStatisticsFilter->GetMaximum( i ); statistics.Mean = labelStatisticsFilter->GetMean( i ); statistics.Median = labelStatisticsFilter->GetMedian( i ); statistics.Sigma = labelStatisticsFilter->GetSigma( i ); statistics.RMS = sqrt( statistics.Mean * statistics.Mean + statistics.Sigma * statistics.Sigma ); } else { *histogram = m_EmptyHistogram; statistics.Reset(); } } 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 ); - // 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) + // 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 ); const mitk::Geometry3D *imageGeometry3D = m_Image->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 + // store the polyline contour as vtkPoints object bool outOfBounds = false; - vtkPoints *points = vtkPoints::New(); + vtkSmartPointer points = vtkPoints::New(); typename PlanarFigure::PolyLineType::const_iterator it; const Vector3D& imageSpacing3D = imageGeometry3D->GetSpacing(); 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 ); - // Ensure correct pixel center - point3D[0] -= 0.5 / imageSpacing3D[0]; - point3D[1] -= 0.5 / imageSpacing3D[1]; - // 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 ); - // Add point to polyline array - points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 ); + points->InsertNextPoint( point3D[i0], point3D[i1], 0 ); } - polyline->SetPoints( points ); - points->Delete(); if ( outOfBounds ) { - polyline->Delete(); throw std::runtime_error( "Figure at least partially outside of image bounds!" ); } - unsigned int numberOfPoints = planarFigurePolyline.size(); - vtkIdType *ptIds = new vtkIdType[numberOfPoints]; - for ( vtkIdType 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->SetInput( 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->SetInput( extrudeFilter->GetOutput() ); - - + // create a vtkLassoStencilSource and set the points of the Polygon + vtkSmartPointer lassoStencil = vtkLassoStencilSource::New(); + lassoStencil->SetShapeToPolygon(); + lassoStencil->SetPoints( points ); // 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() ); - vtkImageImport *vtkImporter = vtkImageImport::New(); + vtkSmartPointer vtkImporter = vtkImageImport::New(); this->ConnectPipelines( itkExporter, vtkImporter ); - vtkImporter->Update(); - // Apply the generated image stencil to the input image - vtkImageStencil *imageStencilFilter = vtkImageStencil::New(); - imageStencilFilter->SetInput( vtkImporter->GetOutput() ); - imageStencilFilter->SetStencil( polyDataToImageStencil->GetOutput() ); + vtkSmartPointer imageStencilFilter = vtkImageStencil::New(); + imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() ); + imageStencilFilter->SetStencil( lassoStencil->GetOutput() ); imageStencilFilter->ReverseStencilOff(); imageStencilFilter->SetBackgroundValue( 0 ); imageStencilFilter->Update(); - // Export from VTK back to ITK - vtkImageExport *vtkExporter = vtkImageExport::New(); - vtkExporter->SetInput( imageStencilFilter->GetOutput() ); + vtkSmartPointer vtkExporter = vtkImageExport::New(); + vtkExporter->SetInputConnection( imageStencilFilter->GetOutputPort() ); vtkExporter->Update(); typename ImageImportType::Pointer itkImporter = ImageImportType::New(); this->ConnectPipelines( vtkExporter, itkImporter ); itkImporter->Update(); - - //typedef itk::ImageFileWriter< MaskImage2DType > FileWriterType; - //FileWriterType::Pointer writer = FileWriterType::New(); - //writer->SetInput( itkImporter->GetOutput() ); - //writer->SetFileName( "c:/test.mhd" ); - //writer->Update(); - - // Store mask m_InternalImageMask2D = itkImporter->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 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/MitkExt/Algorithms/mitkImageStatisticsCalculator.h b/Modules/MitkExt/Algorithms/mitkImageStatisticsCalculator.h index 0cfa7a3a8e..15c36a257d 100644 --- a/Modules/MitkExt/Algorithms/mitkImageStatisticsCalculator.h +++ b/Modules/MitkExt/Algorithms/mitkImageStatisticsCalculator.h @@ -1,298 +1,298 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date: 2009-05-12 19:56:03 +0200 (Di, 12 Mai 2009) $ Version: $Revision: 17179 $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #ifndef _MITK_IMAGESTATISTICSCALCULATOR_H #define _MITK_IMAGESTATISTICSCALCULATOR_H #include #include "MitkExtExports.h" #include #include #include "mitkImage.h" #include "mitkImageTimeSelector.h" #include "mitkPlanarFigure.h" - +#include namespace mitk { /** * \brief Class for calculating statistics and histogram for an (optionally * masked) image. * * Images can be masked by either a (binary) image (of the same dimensions as * the original image) or by a closed mitk::PlanarFigure, e.g. a circle or * polygon. When masking with a planar figure, the slice corresponding to the * plane containing the figure is extracted and then clipped with contour * defined by the figure. Planar figures need to be aligned along the main axes * of the image (transversal, sagittal, coronal). Planar figures on arbitrary * rotated planes are not supported. * * For each operating mode (no masking, masking by image, masking by planar * figure), the calculated statistics and histogram are cached so that, when * switching back and forth between operation modes without modifying mask or * image, the information doesn't need to be recalculated. * * Note: currently time-resolved and multi-channel pictures are not properly * supported. */ class MitkExt_EXPORT ImageStatisticsCalculator : public itk::Object { public: enum { MASKING_MODE_NONE = 0, MASKING_MODE_IMAGE, MASKING_MODE_PLANARFIGURE }; typedef mitk::Image::HistogramType HistogramType; typedef mitk::Image::HistogramType::ConstIterator HistogramConstIteratorType; struct Statistics { unsigned int N; double Min; double Max; double Mean; double Median; double Variance; double Sigma; double RMS; void Reset() { N = 0; Min = 0.0; Max = 0.0; Mean = 0.0; Median = 0.0; Variance = 0.0; Sigma = 0.0; RMS = 0.0; } }; mitkClassMacro( ImageStatisticsCalculator, itk::Object ); itkNewMacro( ImageStatisticsCalculator ); /** \brief Set image from which to compute statistics. */ void SetImage( const mitk::Image *image ); /** \brief Set binary image for masking. */ void SetImageMask( const mitk::Image *imageMask ); /** \brief Set planar figure for masking. */ void SetPlanarFigure( const mitk::PlanarFigure *planarFigure ); /** \brief Set/Get operation mode for masking */ void SetMaskingMode( unsigned int mode ); /** \brief Set/Get operation mode for masking */ itkGetMacro( MaskingMode, unsigned int ); /** \brief Set/Get operation mode for masking */ void SetMaskingModeToNone(); /** \brief Set/Get operation mode for masking */ void SetMaskingModeToImage(); /** \brief Set/Get operation mode for masking */ void SetMaskingModeToPlanarFigure(); /** \brief Set a pixel value for pixels that will be ignored in the statistics */ void SetIgnorePixelValue(double value); /** \brief Get the pixel value for pixels that will be ignored in the statistics */ double GetIgnorePixelValue(); /** \brief Set wether a pixel value should be ignored in the statistics */ void SetDoIgnorePixelValue(bool doit); /** \brief Get wether a pixel value will be ignored in the statistics */ bool GetDoIgnorePixelValue(); /** \brief Compute statistics (together with histogram) for the current * masking mode. * * Computation is not executed if statistics is already up to date. In this * case, false is returned; otherwise, true.*/ virtual bool ComputeStatistics( unsigned int timeStep = 0 ); /** \brief Retrieve the histogram depending on the current masking mode. */ const HistogramType *GetHistogram( unsigned int timeStep = 0 ) const; /** \brief Retrieve statistics depending on the current masking mode. */ const Statistics &GetStatistics( unsigned int timeStep = 0 ) const; protected: typedef std::vector< HistogramType::ConstPointer > HistogramVectorType; typedef std::vector< Statistics > StatisticsVectorType; typedef std::vector< itk::TimeStamp > TimeStampVectorType; typedef std::vector< bool > BoolVectorType; typedef itk::Image< unsigned short, 3 > MaskImage3DType; typedef itk::Image< unsigned short, 2 > MaskImage2DType; ImageStatisticsCalculator(); virtual ~ImageStatisticsCalculator(); /** \brief Depending on the masking mode, the image and mask from which to * calculate statistics is extracted from the original input image and mask * data. * * For example, a when using a PlanarFigure as mask, the 2D image slice * corresponding to the PlanarFigure will be extracted from the original * image. If masking is disabled, the original image is simply passed * through. */ void ExtractImageAndMask( unsigned int timeStep = 0 ); /** \brief If the passed vector matches any of the three principal axes * of the passed geometry, the ínteger value corresponding to the axis * is set and true is returned. */ bool GetPrincipalAxis( const Geometry3D *geometry, Vector3D vector, unsigned int &axis ); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsUnmasked( const itk::Image< TPixel, VImageDimension > *image, Statistics &statistics, typename HistogramType::ConstPointer *histogram ); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsMasked( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage, Statistics &statistics, typename HistogramType::ConstPointer *histogram ); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateMaskFromPlanarFigure( const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ); template < typename TPixel, unsigned int VImageDimension > void InternalMaskIgnoredPixels( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage ); /** Connection from ITK to VTK */ template - void ConnectPipelines(ITK_Exporter exporter, VTK_Importer* importer) + void ConnectPipelines(ITK_Exporter exporter, vtkSmartPointer importer) { importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); importer->SetSpacingCallback(exporter->GetSpacingCallback()); importer->SetOriginCallback(exporter->GetOriginCallback()); importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); importer->SetCallbackUserData(exporter->GetCallbackUserData()); } /** Connection from VTK to ITK */ template - void ConnectPipelines(VTK_Exporter* exporter, ITK_Importer importer) + void ConnectPipelines(vtkSmartPointer exporter, ITK_Importer importer) { importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); importer->SetSpacingCallback(exporter->GetSpacingCallback()); importer->SetOriginCallback(exporter->GetOriginCallback()); importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); importer->SetCallbackUserData(exporter->GetCallbackUserData()); } void UnmaskedStatisticsProgressUpdate(); void MaskedStatisticsProgressUpdate(); mitk::Image::ConstPointer m_Image; mitk::Image::ConstPointer m_ImageMask; mitk::PlanarFigure::ConstPointer m_PlanarFigure; HistogramVectorType m_ImageHistogramVector; HistogramVectorType m_MaskedImageHistogramVector; HistogramVectorType m_PlanarFigureHistogramVector; HistogramType::Pointer m_EmptyHistogram; StatisticsVectorType m_ImageStatisticsVector; StatisticsVectorType m_MaskedImageStatisticsVector; StatisticsVectorType m_PlanarFigureStatisticsVector; Statistics m_EmptyStatistics; unsigned int m_MaskingMode; bool m_MaskingModeChanged; mitk::Image::ConstPointer m_InternalImage; MaskImage3DType::Pointer m_InternalImageMask3D; MaskImage2DType::Pointer m_InternalImageMask2D; TimeStampVectorType m_ImageStatisticsTimeStampVector; TimeStampVectorType m_MaskedImageStatisticsTimeStampVector; TimeStampVectorType m_PlanarFigureStatisticsTimeStampVector; BoolVectorType m_ImageStatisticsCalculationTriggerVector; BoolVectorType m_MaskedImageStatisticsCalculationTriggerVector; BoolVectorType m_PlanarFigureStatisticsCalculationTriggerVector; double m_IgnorePixelValue; bool m_DoIgnorePixelValue; bool m_IgnorePixelValueChanged; }; } #endif // #define _MITK_IMAGESTATISTICSCALCULATOR_H diff --git a/Modules/MitkExt/Algorithms/mitkvtkImageStencilRaster.h b/Modules/MitkExt/Algorithms/mitkvtkImageStencilRaster.h new file mode 100644 index 0000000000..1815d054fa --- /dev/null +++ b/Modules/MitkExt/Algorithms/mitkvtkImageStencilRaster.h @@ -0,0 +1,98 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkImageStencilData.h + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +// .NAME vtkImageStencilData - efficient description of an image stencil +// .SECTION Description +// vtkImageStencilData describes an image stencil in a manner which is +// efficient both in terms of speed and storage space. The stencil extents +// are stored for each x-row across the image (multiple extents per row if +// necessary) and can be retrieved via the GetNextExtent() method. +// .SECTION see also +// vtkImageStencilSource vtkImageStencil + +#ifndef __vtkImageStencilRaster_h +#define __vtkImageStencilRaster_h + +#include "vtkImageStencilData.h" +#include "vtkDataObject.h" +#include "MitkExtExports.h" + + +//BTX +// Description: +// This is a helper class for stencil creation. It is a raster with +// infinite resolution in the X direction (approximately, since it uses +// double precision). Lines that represent polygon edges can be drawn +// into this raster, and then filled given a tolerance. +class MitkExt_EXPORT vtkImageStencilRaster +{ +public: + // Description: + // Create a raster with the specified whole y extent. + vtkImageStencilRaster(const int wholeExtent[2]); + + // Description: + // Destructor. + ~vtkImageStencilRaster(); + + // Description: + // Reset the raster to its original state, but keep the same whole + // extent. Pre-allocate the specified 1D allocateExtent, which must be + // within the whole extent. + void PrepareForNewData(const int allocateExtent[2] = 0); + + // Description: + // Insert a line into the raster, given the two end points. + // The "inflect1" and "inflect2" should be set if you want + // to add a small vertical tolerance to either endpoints. + void InsertLine(const double p1[2], const double p2[2], + bool inflect1, bool inflect2); + + // Description: + // Fill the specified extent of a vtkImageStencilData with the raster, + // after permuting the raster according to xj and yj. + void FillStencilData(vtkImageStencilData *data, const int extent[6], + int xj = 0, int yj = 1); + + // Description: + // The tolerance for float-to-int conversions. + void SetTolerance(double tol) { this->Tolerance = tol; } + double GetTolerance() { return this->Tolerance; } + +protected: + // Description: + // Ensure that the raster is initialized for the specified range + // of y values, which must be within the Extent. + void PrepareExtent(int ymin, int ymax); + + // Description: + // Insert an x point into the raster. If the y value is larger + // than the y extent, the extent will grow automatically. + void InsertPoint(int y, double x); + + int Extent[2]; + int UsedExtent[2]; + double **Raster; + double Tolerance; + +private: + vtkImageStencilRaster(const vtkImageStencilRaster&); // Not implemented. + void operator=(const vtkImageStencilRaster&); // Not implemented. +}; +//ETX + +#endif + + + diff --git a/Modules/MitkExt/Algorithms/mitkvtkLassoStencilSource.h b/Modules/MitkExt/Algorithms/mitkvtkLassoStencilSource.h new file mode 100644 index 0000000000..86ada7a3b7 --- /dev/null +++ b/Modules/MitkExt/Algorithms/mitkvtkLassoStencilSource.h @@ -0,0 +1,107 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkLassoStencilSource.h + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +// .NAME vtkLassoStencilSource - Create a stencil from a contour +// .SECTION Description +// vtkLassoStencilSource will create an image stencil from a +// set of points that define a contour. Its output can be +// used with vtkImageStecil or other vtk classes that apply +// a stencil to an image. +// .SECTION See Also +// vtkROIStencilSource vtkPolyDataToImageStencil +// .SECTION Thanks +// Thanks to David Gobbi for contributing this class to VTK. + +#ifndef __vtkLassoStencilSource_h +#define __vtkLassoStencilSource_h + + +#include "vtkImageStencilSource.h" +#include "MitkExtExports.h" + +class vtkPoints; +class vtkSpline; +class vtkLSSPointMap; + +class MitkExt_EXPORT vtkLassoStencilSource : public vtkImageStencilSource +{ +public: + static vtkLassoStencilSource *New(); + vtkTypeMacro(vtkLassoStencilSource, vtkImageStencilSource); + void PrintSelf(ostream& os, vtkIndent indent); + +//BTX + enum { + POLYGON = 0, + SPLINE = 1, + }; +//ETX + + // Description: + // The shape to use, default is "Polygon". The spline is a + // cardinal spline. Bezier splines are not yet supported. + vtkGetMacro(Shape, int); + vtkSetClampMacro(Shape, int, POLYGON, SPLINE); + void SetShapeToPolygon() { this->SetShape(POLYGON); }; + void SetShapeToSpline() { this->SetShape(SPLINE); }; + virtual const char *GetShapeAsString(); + + // Description: + // The points that make up the lassoo. The loop does not + // have to be closed, the last point will automatically be + // connected to the first point by a straight line segment. + virtual void SetPoints(vtkPoints *points); + vtkGetObjectMacro(Points, vtkPoints); + + // Description: + // The slice orientation. The default is 2, which is XY. + // Other values are 0, which is YZ, and 1, which is XZ. + vtkGetMacro(SliceOrientation, int); + vtkSetClampMacro(SliceOrientation, int, 0, 2); + + // Description: + // The points for a particular slice. This will override the + // points that were set by calling SetPoints() for the slice. + // To clear the setting, call SetSlicePoints(slice, NULL). + virtual void SetSlicePoints(int i, vtkPoints *points); + virtual vtkPoints *GetSlicePoints(int i); + + // Description: + // Remove points from all slices. + virtual void RemoveAllSlicePoints(); + + // Description: + // Overload GetMTime() to include the timestamp on the points. + unsigned long GetMTime(); + +protected: + vtkLassoStencilSource(); + ~vtkLassoStencilSource(); + + virtual int RequestData(vtkInformation *, vtkInformationVector **, + vtkInformationVector *); + + int Shape; + int SliceOrientation; + vtkPoints *Points; + vtkSpline *SplineX; + vtkSpline *SplineY; + vtkLSSPointMap *PointMap; + +private: + vtkLassoStencilSource(const vtkLassoStencilSource&); // Not implemented. + void operator=(const vtkLassoStencilSource&); // Not implemented. +}; + +#endif diff --git a/Modules/MitkExt/Algorithms/vtkImageStencilRaster.cxx b/Modules/MitkExt/Algorithms/vtkImageStencilRaster.cxx new file mode 100644 index 0000000000..6b00e1f26f --- /dev/null +++ b/Modules/MitkExt/Algorithms/vtkImageStencilRaster.cxx @@ -0,0 +1,452 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkImageStencilData.cxx + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +#include "mitkvtkImageStencilRaster.h" + +#include "vtkImageStencilSource.h" +#include "vtkInformation.h" +#include "vtkInformationVector.h" +#include "vtkDemandDrivenPipeline.h" +#include "vtkDataSetAttributes.h" +#include "vtkDataArray.h" +#include "vtkObjectFactory.h" +#include "vtkMath.h" + +#include +#include + +//---------------------------------------------------------------------------- +// tolerance for float-to-int conversion in stencil operations + +#define VTK_STENCIL_TOL 7.62939453125e-06 + +//---------------------------------------------------------------------------- +vtkImageStencilRaster::vtkImageStencilRaster(const int extent[2]) +{ + int rsize = extent[1] - extent[0] + 1; + + // the "raster" is a sequence of pointer-pairs, where the first pointer + // points to the first value in the raster line, and the second pointer + // points to one location past the last vale in the raster line. The + // difference is the number of x values stored in the raster line. + this->Raster = new double*[2*static_cast(rsize)]; + + // the extent is the range of y values (one line per y) + this->Extent[0] = extent[0]; + this->Extent[1] = extent[1]; + + // tolerance should be larger than expected roundoff errors + this->Tolerance = VTK_STENCIL_TOL; + + // no extent is used initially + this->UsedExtent[0] = 0; + this->UsedExtent[1] = -1; +} + +//---------------------------------------------------------------------------- +vtkImageStencilRaster::~vtkImageStencilRaster() +{ + if (this->UsedExtent[1] >= this->UsedExtent[0]) + { + size_t i = 2*static_cast(this->UsedExtent[0] - this->Extent[0]); + size_t imax = 2*static_cast(this->UsedExtent[1] - this->Extent[0]); + + do + { + if (this->Raster[i]) + { + delete [] this->Raster[i]; + } + i += 2; + } + while (i <= imax); + } + delete [] this->Raster; +} + +//---------------------------------------------------------------------------- +void vtkImageStencilRaster::PrepareForNewData(const int allocateExtent[2]) +{ + if (this->UsedExtent[1] >= this->UsedExtent[0]) + { + // reset and re-use the allocated raster lines + size_t i = 2*static_cast(this->UsedExtent[0]-this->Extent[0]); + size_t imax=2*static_cast(this->UsedExtent[1]-this->Extent[0]); + do + { + this->Raster[i+1] = this->Raster[i]; + i += 2; + } + while (i <= imax); + } + + if (allocateExtent && allocateExtent[1] >= allocateExtent[1]) + { + this->PrepareExtent(allocateExtent[0], allocateExtent[1]); + } +} + +//---------------------------------------------------------------------------- +void vtkImageStencilRaster::PrepareExtent(int ymin, int ymax) +{ + // this does not do any allocation, it just initializes any + // raster lines are not already part of the UsedExtent, and + // then expands the UsedExtent to include [ymin, ymax] + + if (this->UsedExtent[1] < this->UsedExtent[0]) + { + size_t i = 2*static_cast(ymin - this->Extent[0]); + size_t imax = 2*static_cast(ymax - this->Extent[0]); + + do + { + this->Raster[i] = 0; + } + while (++i <= imax); + + this->UsedExtent[0] = ymin; + this->UsedExtent[1] = ymax; + + return; + } + + if (ymin < this->UsedExtent[0]) + { + size_t i = 2*static_cast(ymin - this->Extent[0]); + size_t imax = 2*static_cast(this->UsedExtent[0]-this->Extent[0]-1); + + do + { + this->Raster[i] = 0; + } + while (++i <= imax); + + this->UsedExtent[0] = ymin; + } + + if (ymax > this->UsedExtent[1]) + { + size_t i = 2*static_cast(this->UsedExtent[1]+1 - this->Extent[0]); + size_t imax = 2*static_cast(ymax - this->Extent[0]); + + do + { + this->Raster[i] = 0; + } + while (++i <= imax); + + this->UsedExtent[1] = ymax; + } +} + +//---------------------------------------------------------------------------- +void vtkImageStencilRaster::InsertPoint(int y, double x) +{ + size_t pos = 2*static_cast(y - this->Extent[0]); + double* &rhead = this->Raster[pos]; + double* &rtail = this->Raster[pos+1]; + + // current size is the diff between the tail and the head + size_t n = rtail - rhead; + + // no allocation on this raster line yet + if (rhead == 0) + { + rhead = new double[2]; + rtail = rhead; + } + // grow whenever size reaches a power of two + else if (n > 1 && (n & (n-1)) == 0) + { + double *ptr = new double[2*n]; + for (size_t i = 0; i < n; i++) + { + ptr[i] = rhead[i]; + } + delete [] rhead; + rhead = ptr; + rtail = ptr + n; + } + + // insert the value + *rtail++ = x; +} + +//---------------------------------------------------------------------------- +void vtkImageStencilRaster::InsertLine( + const double pt1[2], const double pt2[2], + bool inflection1, bool inflection2) +{ + double x1 = pt1[0]; + double x2 = pt2[0]; + double y1 = pt1[1]; + double y2 = pt2[1]; + + // swap end points if necessary + if (y1 > y2) + { + x1 = pt2[0]; + x2 = pt1[0]; + y1 = pt2[1]; + y2 = pt1[1]; + bool tmp = inflection1; + inflection1 = inflection2; + inflection2 = tmp; + } + + // find min and max of x values + double xmin = x1; + double xmax = x2; + if (x1 > x2) + { + xmin = x2; + xmax = x1; + } + + // check for parallel to the x-axis + if (y1 == y2) + { + return; + } + + double ymin = y1; + double ymax = y2; + + // if an end is an inflection point, include a tolerance + ymin -= inflection1*this->Tolerance; + ymax += inflection2*this->Tolerance; + + // Integer y values for start and end of line + int iy1, iy2; + iy1 = this->Extent[0]; + iy2 = this->Extent[1]; + + // Check for out of bounds + if (ymax < iy1 || ymin >= iy2) + { + return; + } + + // Guard against extentY + if (ymin >= iy1) + { + iy1 = vtkMath::Floor(ymin) + 1; + } + if (ymax < iy2) + { + iy2 = vtkMath::Floor(ymax); + } + + // Expand allocated extent if necessary + if (iy1 < this->UsedExtent[0] || + iy2 > this->UsedExtent[1]) + { + this->PrepareExtent(iy1, iy2); + } + + // Precompute values for a Bresenham-like line algorithm + double grad = (x2 - x1)/(y2 - y1); + double delta = (iy1 - y1)*grad; + + // Go along y and place each x in the proper raster line + for (int y = iy1; y <= iy2; y++) + { + double x = x1 + delta; + // incrementing delta has less roundoff error than incrementing x, + // since delta will typically be smaller than x + delta += grad; + + // clamp x (because of tolerance, it might not be in range) + if (x < xmin) + { + x = xmin; + } + else if (x > xmax) + { + x = xmax; + } + + this->InsertPoint(y, x); + } +} + +//---------------------------------------------------------------------------- +void vtkImageStencilRaster::FillStencilData( + vtkImageStencilData *data, const int extent[6], int xj, int yj) +{ + if (xj != 0) + { + // slices are stacked in the x direction + int xmin = extent[2*xj]; + int xmax = extent[2*xj+1]; + int ymin = this->UsedExtent[0]; + int ymax = this->UsedExtent[1]; + int zmin = extent[0]; + int zmax = extent[1]; + + for (int idY = ymin; idY <= ymax; idY++) + { + size_t pos = 2*static_cast(idY - this->Extent[0]); + double *rline = this->Raster[pos]; + double *rlineEnd = this->Raster[pos+1]; + + if (rline == 0) + { + continue; + } + + vtkstd::sort(rline, rlineEnd); + + int xy[2]; + xy[2-xj] = idY; + + int lastr = VTK_INT_MIN; + + size_t l = rlineEnd - rline; + l = l - (l & 1); // force l to be an even number + for (size_t k = 0; k < l; k += 2) + { + double x1 = rline[k] - this->Tolerance; + double x2 = rline[k+1] + this->Tolerance; + + // make sure one of the ends is in bounds + if (x2 < xmin || x1 >= xmax) + { + continue; + } + + // clip the line segment with the bounds + int r1 = xmin; + int r2 = xmax; + + if (x1 >= xmin) + { + r1 = vtkMath::Floor(x1) + 1; + } + if (x2 < xmax) + { + r2 = vtkMath::Floor(x2); + } + + // ensure no overlap occurs with previous + if (r1 <= lastr) + { + r1 = lastr + 1; + } + lastr = r2; + + for (int idX = r1; idX <= r2; idX++) + { + xy[xj-1] = idX; + data->InsertNextExtent(zmin, zmax, xy[0], xy[1]); + } + } + } + } + else + { + // slices stacked in the y or z direction + int zj = 3 - yj; + int xmin = extent[0]; + int xmax = extent[1]; + int ymin = this->UsedExtent[0]; + int ymax = this->UsedExtent[1]; + int zmin = extent[2*zj]; + int zmax = extent[2*zj+1]; + + // convert each raster line into extents for the stencil + for (int idY = ymin; idY <= ymax; idY++) + { + size_t pos = 2*static_cast(idY - this->Extent[0]); + double *rline = this->Raster[pos]; + double *rlineEnd = this->Raster[pos+1]; + + if (rline == 0) + { + continue; + } + + vtkstd::sort(rline, rlineEnd); + + int lastr = VTK_INT_MIN; + + // go through each raster line and fill the stencil + size_t l = rlineEnd - rline; + l = l - (l & 1); // force l to be an even number + for (size_t k = 0; k < l; k += 2) + { + int yz[2]; + + yz[yj-1] = idY; + yz[2-yj] = zmin; + + double x1 = rline[k] - this->Tolerance; + double x2 = rline[k+1] + this->Tolerance; + + if (x2 < xmin || x1 >= xmax) + { + continue; + } + + int r1 = xmin; + int r2 = xmax; + + if (x1 >= xmin) + { + r1 = vtkMath::Floor(x1) + 1; + } + if (x2 < xmax) + { + r2 = vtkMath::Floor(x2); + } + + // ensure no overlap occurs between extents + if (r1 <= lastr) + { + r1 = lastr + 1; + } + lastr = r2; + + if (r2 >= r1) + { + data->InsertNextExtent(r1, r2, yz[0], yz[1]); + } + } + } + + // copy the result to all other slices + if (zmin < zmax) + { + for (int idY = ymin; idY <= ymax; idY++) + { + int r1, r2; + int yz[2]; + + yz[yj-1] = idY; + yz[2-yj] = zmin; + + int iter = 0; + while (data->GetNextExtent(r1, r2, xmin, xmax, yz[0], yz[1], iter)) + { + for (int idZ = zmin + 1; idZ <= zmax; idZ++) + { + yz[2-yj] = idZ; + data->InsertNextExtent(r1, r2, yz[0], yz[1]); + } + yz[2-yj] = zmin; + } + } + } + } +} diff --git a/Modules/MitkExt/Algorithms/vtkLassoStencilSource.cxx b/Modules/MitkExt/Algorithms/vtkLassoStencilSource.cxx new file mode 100644 index 0000000000..dbe66d903d --- /dev/null +++ b/Modules/MitkExt/Algorithms/vtkLassoStencilSource.cxx @@ -0,0 +1,609 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkLassoStencilSource.cxx + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +#include "mitkvtkLassoStencilSource.h" + +#include "vtkMath.h" +#include "vtkPoints.h" +#include "vtkCardinalSpline.h" +#include "vtkDataArray.h" +#include "vtkImageData.h" +#include "vtkImageStencilData.h" +#include "vtkInformation.h" +#include "vtkInformationVector.h" +#include "vtkObjectFactory.h" +#include "vtkStreamingDemandDrivenPipeline.h" +#include "vtkSmartPointer.h" + +#include "mitkvtkImageStencilRaster.h" + +#include +#include + +vtkStandardNewMacro(vtkLassoStencilSource); +vtkCxxSetObjectMacro(vtkLassoStencilSource, Points, vtkPoints); + +//---------------------------------------------------------------------------- +class vtkLSSPointMap : public vtkstd::map > +{ +}; + +//---------------------------------------------------------------------------- +vtkLassoStencilSource::vtkLassoStencilSource() +{ + this->SetNumberOfInputPorts(0); + + this->Shape = vtkLassoStencilSource::POLYGON; + this->SliceOrientation = 2; + this->Points = NULL; + this->SplineX = vtkCardinalSpline::New(); + this->SplineY = vtkCardinalSpline::New(); + + this->PointMap = new vtkLSSPointMap(); +} + +//---------------------------------------------------------------------------- +vtkLassoStencilSource::~vtkLassoStencilSource() +{ + this->SetPoints(NULL); + if (this->SplineX) + { + this->SplineX->Delete(); + this->SplineX = NULL; + } + if (this->SplineY) + { + this->SplineY->Delete(); + this->SplineY = NULL; + } + if (this->PointMap) + { + delete this->PointMap; + this->PointMap = NULL; + } +} + +//---------------------------------------------------------------------------- +void vtkLassoStencilSource::PrintSelf(ostream& os, vtkIndent indent) +{ + this->Superclass::PrintSelf(os,indent); + + os << indent << "Shape: " << this->GetShapeAsString() << "\n"; + os << indent << "Points: " << this->Points << "\n"; + os << indent << "SliceOrientation: " << this->GetSliceOrientation() << "\n"; + os << indent << "SlicePoints: " << this->PointMap->size() << "\n"; +} + +//---------------------------------------------------------------------------- +const char *vtkLassoStencilSource::GetShapeAsString() +{ + switch (this->Shape) + { + case vtkLassoStencilSource::POLYGON: + return "Polygon"; + case vtkLassoStencilSource::SPLINE: + return "Spline"; + } + return ""; +} + +//---------------------------------------------------------------------------- +unsigned long int vtkLassoStencilSource::GetMTime() +{ + unsigned long mTime = this->vtkImageStencilSource::GetMTime(); + + if ( this->Points != NULL ) + { + unsigned long t = this->Points->GetMTime(); + if (t > mTime) + { + mTime = t; + } + } + + if ( !this->PointMap->empty() ) + { + vtkLSSPointMap::iterator iter = this->PointMap->begin(); + while ( iter != this->PointMap->end() ) + { + unsigned long t = iter->second->GetMTime(); + if (t > mTime) + { + mTime = t; + } + iter++; + } + } + + return mTime; +} + +//---------------------------------------------------------------------------- +void vtkLassoStencilSource::SetSlicePoints(int i, vtkPoints *points) +{ + vtkLSSPointMap::iterator iter = this->PointMap->find(i); + if (iter != this->PointMap->end()) + { + if (iter->second == points) + { + return; + } + else if (points == 0) + { + this->PointMap->erase(iter); + } + else + { + iter->second = points; + } + } + else + { + if (points == NULL) + { + return; + } + else + { + this->PointMap->insert(iter, vtkLSSPointMap::value_type(i, points)); + } + } + + this->Modified(); +} + +//---------------------------------------------------------------------------- +void vtkLassoStencilSource::RemoveAllSlicePoints() +{ + this->PointMap->clear(); +} + +//---------------------------------------------------------------------------- +vtkPoints *vtkLassoStencilSource::GetSlicePoints(int i) +{ + vtkLSSPointMap::iterator iter = this->PointMap->find(i); + if (iter != this->PointMap->end()) + { + return iter->second; + } + return NULL; +} + +//---------------------------------------------------------------------------- +// tolerance for stencil operations + +#define VTK_STENCIL_TOL 7.62939453125e-06 + +//---------------------------------------------------------------------------- +// Compute a reduced extent based on the bounds of the shape. +static void vtkLassoStencilSourceSubExtent( + vtkPoints *points, + const double origin[3], const double spacing[3], + const int extent[6], int subextent[6]) +{ + double bounds[6]; + points->GetBounds(bounds); + + for (int i = 0; i < 3; i++) + { + double emin = (bounds[2*i] - origin[i])/spacing[i] - VTK_STENCIL_TOL; + double emax = (bounds[2*i+1] - origin[i])/spacing[i] + VTK_STENCIL_TOL; + + subextent[2*i] = extent[2*i]; + subextent[2*i+1] = extent[2*i+1]; + + if (extent[2*i] < emin) + { + subextent[2*i] = VTK_INT_MAX; + if (extent[2*i+1] >= emin) + { + subextent[2*i] = vtkMath::Floor(emin) + 1; + } + } + + if (extent[2*i+1] > emax) + { + subextent[2*i+1] = VTK_INT_MIN; + if (extent[2*i] <= emax) + { + subextent[2*i+1] = vtkMath::Floor(emax); + } + } + + } +} + +//---------------------------------------------------------------------------- +// Rasterize a polygon into the stencil +static int vtkLassoStencilSourcePolygon( + vtkPoints *points, vtkImageStencilData *data, vtkImageStencilRaster *raster, + const int extent[6], const double origin[3], const double spacing[3], + int xj, int yj) +{ + // get the bounds of the polygon + int subextent[6]; + vtkLassoStencilSourceSubExtent(points, origin, spacing, extent, subextent); + + // allocate the raster + raster->PrepareForNewData(&subextent[2*yj]); + + // rasterize each line + vtkIdType n = points->GetNumberOfPoints(); + double p[3]; + double p0[2], p1[2], p2[2], p3[2]; + + points->GetPoint(n-1, p); + p0[0] = (p[xj] - origin[xj])/spacing[xj]; + p0[1] = (p[yj] - origin[yj])/spacing[yj]; + + points->GetPoint(0, p); + p1[0] = (p[xj] - origin[xj])/spacing[xj]; + p1[1] = (p[yj] - origin[yj])/spacing[yj]; + + double dx = p1[0] - p0[0]; + double dy = p1[1] - p0[1]; + if (dx*dx + dy*dy <= VTK_STENCIL_TOL*VTK_STENCIL_TOL) + { + n -= 1; + points->GetPoint(n-1, p); + p0[0] = (p[xj] - origin[xj])/spacing[xj]; + p0[1] = (p[yj] - origin[yj])/spacing[yj]; + } + + points->GetPoint(1, p); + p2[0] = (p[xj] - origin[xj])/spacing[xj]; + p2[1] = (p[yj] - origin[yj])/spacing[yj]; + + // inflection means the line changes vertical direction + bool inflection1, inflection2; + inflection1 = ( (p1[1] - p0[1])*(p2[1] - p1[1]) <= 0 ); + + for (vtkIdType i = 0; i < n; i++) + { + points->GetPoint((i+2)%n, p); + p3[0] = (p[xj] - origin[xj])/spacing[xj]; + p3[1] = (p[yj] - origin[yj])/spacing[yj]; + + inflection2 = ( (p2[1] - p1[1])*(p3[1] - p2[1]) <= 0 ); + + raster->InsertLine(p1, p2, inflection1, inflection2); + + p0[0] = p1[0]; p0[1] = p1[1]; + p1[0] = p2[0]; p1[1] = p2[1]; + p2[0] = p3[0]; p2[1] = p3[1]; + inflection1 = inflection2; + } + + raster->FillStencilData(data, extent, xj, yj); + + return 1; +} + + +//---------------------------------------------------------------------------- +// Generate the splines for the given set of points. The splines +// will be closed if the final point is equal to the first point. +// The parametric value for the resulting spline will be valid over +// the range [0, tmax] where the tmax value is returned by reference. +static void vtkLassoStencilSourceCreateSpline(vtkPoints *points, + const double origin[3], const double spacing[3], + int xj, int yj, vtkSpline *xspline, vtkSpline *yspline, + double &tmax, double &dmax) +{ + // initialize the spline + xspline->RemoveAllPoints(); + yspline->RemoveAllPoints(); + xspline->ClosedOff(); + yspline->ClosedOff(); + + // get the number of points and the first/last point + vtkIdType n = points->GetNumberOfPoints(); + double p[3]; + double p0[2], p1[2]; + + points->GetPoint(n-1, p); + p0[0] = (p[xj] - origin[xj])/spacing[xj]; + p0[1] = (p[yj] - origin[yj])/spacing[yj]; + + points->GetPoint(0, p); + p1[0] = (p[xj] - origin[xj])/spacing[xj]; + p1[1] = (p[yj] - origin[yj])/spacing[yj]; + + // factor between real distance and parametric distance + double f = 1.0; + // the length of the implicit segment for closed loops + double lastd = 0; + + // aspect ratio + double xf = 1.0; + double yf = 1.0; + if (spacing[xj] > spacing[yj]) + { + xf = spacing[xj]/spacing[yj]; + } + else + { + yf = spacing[yj]/spacing[xj]; + } + + // if first and last point are same, spline is closed + double dx = (p1[0] - p0[0])*xf; + double dy = (p1[1] - p0[1])*yf; + double d2 = dx*dx + dy*dy; + while (d2 <= VTK_STENCIL_TOL*VTK_STENCIL_TOL && n > 1) + { + n -= 1; + points->GetPoint(n-1, p); + p0[0] = (p[xj] - origin[xj])/spacing[xj]; + p0[1] = (p[yj] - origin[yj])/spacing[yj]; + + xspline->ClosedOn(); + yspline->ClosedOn(); + + // vtkSpline considers the parametric length of the implicit + // segment of closed loops to be unity, so set "f" so that + // multiplying the real length of that segment by "f" gives unity. + dx = (p1[0] - p0[0])*xf; + dy = (p1[1] - p0[1])*yf; + d2 = dx*dx + dy*dy; + lastd = sqrt(d2); + if (lastd > 0) + { + f = 1.0/lastd; + } + } + + // Add all the points to the spline. + double d = 0.0; + for (vtkIdType i = 0; i < n; i++) + { + p0[0] = p1[0]; p0[1] = p1[1]; + + points->GetPoint(i, p); + p1[0] = (p[xj] - origin[xj])/spacing[xj]; + p1[1] = (p[yj] - origin[yj])/spacing[yj]; + + dx = (p1[0] - p0[0])*xf; + dy = (p1[1] - p0[1])*yf; + + d += sqrt(dx*dx + dy*dy); + + double t = f*d; + + xspline->AddPoint(t, p1[0]); + yspline->AddPoint(t, p1[1]); + } + + // Do the spline precomputations + xspline->Compute(); + yspline->Compute(); + + // The spline is valid over t = [0, tmax] + d += lastd; + tmax = f*d; + dmax = d; +} + +//---------------------------------------------------------------------------- +// Rasterize a spline contour into the stencil +static int vtkLassoStencilSourceSpline( + vtkPoints *points, vtkImageStencilData *data, vtkImageStencilRaster *raster, + const int extent[6], const double origin[3], const double spacing[3], + int xj, int yj, vtkSpline *xspline, vtkSpline *yspline) +{ + // create the splines + double tmax, dmax; + vtkLassoStencilSourceCreateSpline( + points, origin, spacing, xj, yj, xspline, yspline, tmax, dmax); + + if (dmax <= VTK_STENCIL_TOL) + { + return 1; + } + + // get the bounds of the polygon as a first guess of the spline bounds + int subextent[6]; + vtkLassoStencilSourceSubExtent(points, origin, spacing, extent, subextent); + + // allocate the raster + raster->PrepareForNewData(&subextent[2*yj]); + + // go around the spline + vtkIdType n = vtkMath::Floor(dmax)+1; + double delta = tmax/n; + + double p0[2], p1[2], p2[2], p3[2]; + + double t = tmax; + if (xspline->GetClosed()) + { + t = (n-1)*tmax/n; + } + else + { + n = n + 1; + } + + p0[0] = xspline->Evaluate(t); + p0[1] = yspline->Evaluate(t); + + t = 0; + p1[0] = xspline->Evaluate(t); + p1[1] = yspline->Evaluate(t); + + t = delta; + p2[0] = xspline->Evaluate(t); + p2[1] = yspline->Evaluate(t); + + // inflection means the line changes vertical direction + bool inflection1, inflection2; + inflection1 = ( (p1[1] - p0[1])*(p2[1] - p1[1]) <= 0 ); + + for (vtkIdType i = 0; i < n; i++) + { + t += delta; + if (i == n-2) + { + t = 0; + } + + p3[0] = xspline->Evaluate(t); + p3[1] = yspline->Evaluate(t); + + inflection2 = ( (p2[1] - p1[1])*(p3[1] - p2[1]) <= 0 ); + + raster->InsertLine(p1, p2, inflection1, inflection2); + + p0[0] = p1[0]; p0[1] = p1[1]; + p1[0] = p2[0]; p1[1] = p2[1]; + p2[0] = p3[0]; p2[1] = p3[1]; + inflection1 = inflection2; + } + + raster->FillStencilData(data, extent, xj, yj); + + return 1; +} + +//---------------------------------------------------------------------------- +static int vtkLassoStencilSourceExecute( + vtkPoints *points, vtkImageStencilData *data, vtkImageStencilRaster *raster, + int extent[6], double origin[3], double spacing[3], int shape, + int xj, int yj, vtkSpline *xspline, vtkSpline *yspline) +{ + int result = 1; + + if (points == 0 || points->GetNumberOfPoints() < 3) + { + return 1; + } + + switch (shape) + { + case vtkLassoStencilSource::POLYGON: + result = vtkLassoStencilSourcePolygon( + points, data, raster, extent, origin, spacing, xj, yj); + break; + case vtkLassoStencilSource::SPLINE: + result = vtkLassoStencilSourceSpline( + points, data, raster, extent, origin, spacing, xj, yj, + xspline, yspline); + break; + } + + return result; +} + + +//---------------------------------------------------------------------------- +int vtkLassoStencilSource::RequestData( + vtkInformation *request, + vtkInformationVector **inputVector, + vtkInformationVector *outputVector) +{ + int extent[6]; + double origin[3]; + double spacing[3]; + int result = 1; + + this->Superclass::RequestData(request, inputVector, outputVector); + + vtkInformation *outInfo = outputVector->GetInformationObject(0); + vtkImageStencilData *data = vtkImageStencilData::SafeDownCast( + outInfo->Get(vtkDataObject::DATA_OBJECT())); + + outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), extent); + outInfo->Get(vtkDataObject::ORIGIN(), origin); + outInfo->Get(vtkDataObject::SPACING(), spacing); + + int slabExtent[6]; + slabExtent[0] = extent[0]; slabExtent[1] = extent[1]; + slabExtent[2] = extent[2]; slabExtent[3] = extent[3]; + slabExtent[4] = extent[4]; slabExtent[5] = extent[5]; + + int xj = 0; + int yj = 1; + int zj = 2; + + if (this->SliceOrientation == 0) + { + xj = 1; + yj = 2; + zj = 0; + } + else if (this->SliceOrientation == 1) + { + xj = 0; + yj = 2; + zj = 1; + } + + vtkImageStencilRaster raster(&extent[2*yj]); + raster.SetTolerance(VTK_STENCIL_TOL); + + int zmin = extent[2*zj]; + int zmax = extent[2*zj+1]; + + vtkLSSPointMap::iterator iter = this->PointMap->lower_bound(zmin); + vtkLSSPointMap::iterator maxiter = this->PointMap->upper_bound(zmax); + + while (iter != maxiter && result != 0) + { + this->SetProgress((slabExtent[2*zj] - zmin)*1.0/(zmax - zmin + 1)); + + int i = iter->first; + vtkPoints *points = iter->second; + + // fill in the slices with no SlicePoints + if (this->Points && i > slabExtent[2*zj]) + { + slabExtent[2*zj+1] = i-1; + + result = vtkLassoStencilSourceExecute( + this->Points, data, &raster, slabExtent, origin, spacing, + this->Shape, xj, yj, this->SplineX, this->SplineY); + } + + // do the slice with its SlicePoints + if (result) + { + slabExtent[2*zj] = i; + slabExtent[2*zj+1] = i; + + result = vtkLassoStencilSourceExecute( + points, data, &raster, slabExtent, origin, spacing, + this->Shape, xj, yj, this->SplineX, this->SplineY); + + slabExtent[2*zj] = slabExtent[2*zj+1] + 1; + } + + ++iter; + } + + this->SetProgress((slabExtent[2*zj] - zmin)*1.0/(zmax - zmin + 1)); + + // fill in the rest + if (result && slabExtent[2*zj] <= zmax) + { + slabExtent[2*zj+1] = zmax; + + result = vtkLassoStencilSourceExecute( + this->Points, data, &raster, slabExtent, origin, spacing, + this->Shape, xj, yj, this->SplineX, this->SplineY); + + this->SetProgress(1.0); + } + + return result; +} diff --git a/Modules/MitkExt/files.cmake b/Modules/MitkExt/files.cmake index 8316b38ec7..f0f21d93a3 100644 --- a/Modules/MitkExt/files.cmake +++ b/Modules/MitkExt/files.cmake @@ -1,196 +1,206 @@ SET(CPP_FILES Algorithms/mitkMaskAndCutRoiImageFilter.cpp Algorithms/mitkBoundingObjectToSegmentationFilter.cpp Algorithms/vtkPointSetSlicer.cxx Algorithms/mitkCoreExtObjectFactory.cpp Algorithms/mitkImageStatisticsCalculator.cpp Algorithms/mitkAngleCorrectByPointFilter.cpp Algorithms/mitkAutoCropImageFilter.cpp Algorithms/mitkBoundingObjectCutter.cpp Algorithms/mitkCalculateSegmentationVolume.cpp Algorithms/mitkContourSetToPointSetFilter.cpp Algorithms/mitkContourUtils.cpp Algorithms/mitkCorrectorAlgorithm.cpp Algorithms/mitkCylindricToCartesianFilter.cpp Algorithms/mitkDiffImageApplier.cpp Algorithms/mitkDopplerToStrainRateFilter.cpp Algorithms/mitkExtractDirectedPlaneImageFilter.cpp Algorithms/mitkExtractDirectedPlaneImageFilterNew.cpp Algorithms/mitkExtractImageFilter.cpp Algorithms/mitkGeometryClipImageFilter.cpp Algorithms/mitkGeometryDataSource.cpp Algorithms/mitkHeightFieldSurfaceClipImageFilter.cpp Algorithms/mitkImageToLookupTableFilter.cpp Algorithms/mitkImageToSurfaceFilter.cpp Algorithms/mitkInterpolateLinesFilter.cpp Algorithms/mitkLabeledImageToSurfaceFilter.cpp Algorithms/mitkLabeledImageVolumeCalculator.cpp Algorithms/mitkLookupTableSource.cpp Algorithms/mitkManualSegmentationToSurfaceFilter.cpp Algorithms/mitkMaskImageFilter.cpp Algorithms/mitkMeshSource.cpp Algorithms/mitkNonBlockingAlgorithm.cpp Algorithms/mitkOverwriteSliceImageFilter.cpp Algorithms/mitkOverwriteDirectedPlaneImageFilter.cpp Algorithms/mitkPadImageFilter.cpp Algorithms/mitkPlaneCutFilter.cpp Algorithms/mitkPlaneFit.cpp Algorithms/mitkPlanesPerpendicularToLinesFilter.cpp Algorithms/mitkPointLocator.cpp Algorithms/mitkPointSetToCurvedGeometryFilter.cpp Algorithms/mitkPointSetToGeometryDataFilter.cpp Algorithms/mitkPointSetIndexToWorldTransformFilter.cpp Algorithms/mitkSurfaceIndexToWorldTransformFilter.cpp Algorithms/mitkPolygonToRingFilter.cpp Algorithms/mitkProbeFilter.cpp Algorithms/mitkSegmentationSink.cpp Algorithms/mitkShapeBasedInterpolationAlgorithm.cpp Algorithms/mitkShowSegmentationAsSurface.cpp Algorithms/mitkSimpleHistogram.cpp Algorithms/mitkSimpleUnstructuredGridHistogram.cpp Algorithms/mitkSurfaceToImageFilter.cpp Algorithms/mitkUnstructuredGridHistogram.cpp Algorithms/mitkUnstructuredGridSource.cpp Algorithms/mitkVolumeVisualizationImagePreprocessor.cpp Controllers/mitkMovieGenerator.cpp Controllers/mitkMultiStepper.cpp Controllers/mitkSegmentationInterpolationController.cpp Controllers/mitkToolManager.cpp DataManagement/mitkAffineTransformationOperation.cpp DataManagement/mitkApplyDiffImageOperation.cpp DataManagement/mitkBoundingObject.cpp DataManagement/mitkBoundingObjectGroup.cpp DataManagement/mitkCellOperation.cpp DataManagement/mitkColorConversions.cpp DataManagement/mitkColorSequence.cpp DataManagement/mitkColorSequenceCycleH.cpp DataManagement/mitkColorSequenceHalfTones.cpp DataManagement/mitkColorSequenceRainbow.cpp DataManagement/mitkCompressedImageContainer.cpp DataManagement/mitkCone.cpp DataManagement/mitkContour.cpp DataManagement/mitkContourSet.cpp DataManagement/mitkCuboid.cpp DataManagement/mitkCylinder.cpp DataManagement/mitkDataStorageSelection.cpp DataManagement/mitkDelegateManager.cpp DataManagement/mitkDrawOperation.cpp DataManagement/mitkEllipsoid.cpp DataManagement/mitkExternAbstractTransformGeometry.cpp DataManagement/mitkExtrudedContour.cpp DataManagement/mitkFrameOfReferenceUIDManager.cpp DataManagement/mitkGridRepresentationProperty.cpp DataManagement/mitkGridVolumeMapperProperty.cpp DataManagement/mitkItkBaseDataAdapter.cpp DataManagement/mitkLabeledImageLookupTable.cpp DataManagement/mitkLineOperation.cpp DataManagement/mitkMesh.cpp DataManagement/mitkObjectSet.cpp DataManagement/mitkOrganTypeProperty.cpp DataManagement/mitkPlaneLandmarkProjector.cpp DataManagement/mitkPlane.cpp DataManagement/mitkPropertyManager.cpp DataManagement/mitkPropertyObserver.cpp DataManagement/mitkSeedsImage.cpp DataManagement/mitkSeedsImageLookupTableSource.cpp DataManagement/mitkSphereLandmarkProjector.cpp # DataManagement/mitkUSLookupTableSource.cpp DataManagement/mitkUnstructuredGrid.cpp DataManagement/mitkVideoSource.cpp DataManagement/vtkObjectSet.cpp IO/mitkObjFileIOFactory.cpp IO/mitkObjFileReader.cpp IO/mitkPACSPlugin.cpp IO/mitkParRecFileIOFactory.cpp IO/mitkParRecFileReader.cpp IO/mitkStlVolumeTimeSeriesIOFactory.cpp IO/mitkStlVolumeTimeSeriesReader.cpp IO/mitkUnstructuredGridVtkWriter.cpp IO/mitkUnstructuredGridVtkWriterFactory.cpp IO/mitkVtkUnstructuredGridIOFactory.cpp IO/mitkVtkUnstructuredGridReader.cpp IO/mitkVtkVolumeTimeSeriesIOFactory.cpp IO/mitkVtkVolumeTimeSeriesReader.cpp Interactions/mitkAutoSegmentationTool.cpp Interactions/mitkConferenceEventMapper.cpp Interactions/mitkConnectPointsInteractor.cpp Interactions/mitkContourInteractor.cpp Interactions/mitkContourTool.cpp #Interactions/mitkCoordinateSupplier.cpp #Interactions/mitkDisplayCoordinateOperation.cpp #Interactions/mitkDisplayInteractor.cpp Interactions/mitkAffineInteractor3D.cpp Interactions/mitkDisplayPointSetInteractor.cpp #Interactions/mitkDisplayVectorInteractor.cpp Interactions/mitkExtrudedContourInteractor.cpp Interactions/mitkFeedbackContourTool.cpp Interactions/mitkInteractionDebug.cpp Interactions/mitkInteractionDebugger.cpp Interactions/mitkPaintbrushTool.cpp Interactions/mitkPointInteractor.cpp Interactions/mitkPointSelectorInteractor.cpp #Interactions/mitkPositionTracker.cpp Interactions/mitkSeedsInteractor.cpp Interactions/mitkSegTool2D.cpp Interactions/mitkSegmentationsProcessingTool.cpp Interactions/mitkSetRegionTool.cpp Interactions/mitkSocketClient.cpp Interactions/mitkSurfaceDeformationInteractor3D.cpp Interactions/mitkSurfaceInteractor.cpp Interactions/mitkTool.cpp Interactions/mitkAddContourTool.cpp Interactions/mitkAutoCropTool.cpp Interactions/mitkBinaryThresholdTool.cpp Interactions/mitkCalculateGrayValueStatisticsTool.cpp Interactions/mitkCalculateVolumetryTool.cpp Interactions/mitkCorrectorTool2D.cpp Interactions/mitkCreateSurfaceTool.cpp Interactions/mitkEraseRegionTool.cpp Interactions/mitkFillRegionTool.cpp Interactions/mitkRegionGrowingTool.cpp Interactions/mitkSubtractContourTool.cpp Interactions/mitkDrawPaintbrushTool.cpp Interactions/mitkErasePaintbrushTool.cpp Interactions/mitkMorphologicTool.cpp Interactions/mitkErodeTool.cpp Interactions/mitkDilateTool.cpp Interactions/mitkOpeningTool.cpp Interactions/mitkClosingTool.cpp Interactions/mitkBinaryThresholdULTool.cpp Interactions/mitkPixelManipulationTool.cpp Interactions/mitkRegionGrow3DTool.cpp Rendering/mitkContourMapper2D.cpp Rendering/mitkContourSetMapper2D.cpp Rendering/mitkContourSetVtkMapper3D.cpp Rendering/mitkContourVtkMapper3D.cpp Rendering/mitkEnhancedPointSetVtkMapper3D.cpp Rendering/mitkImageBackground2D.cpp Rendering/mitkLineMapper2D.cpp # Rendering/mitkLineVtkMapper3D.cpp Rendering/mitkMeshMapper2D.cpp Rendering/mitkMeshVtkMapper3D.cpp Rendering/mitkNativeRenderWindowInteractor.cpp Rendering/mitkSplineMapper2D.cpp Rendering/mitkSplineVtkMapper3D.cpp Rendering/mitkUnstructuredGridMapper2D.cpp Rendering/mitkUnstructuredGridVtkMapper3D.cpp Rendering/mitkVectorImageMapper2D.cpp Rendering/vtkUnstructuredGridMapper.cpp Rendering/vtkMaskedGlyph2D.cpp Rendering/vtkMaskedGlyph3D.cpp Rendering/vtkMitkVolumeTextureMapper3D.cpp Rendering/vtkMitkOpenGLVolumeTextureMapper3D.cpp Rendering/mitkGPUVolumeMapper3D.cpp Rendering/vtkMitkGPUVolumeRayCastMapper.cpp Rendering/vtkMitkOpenGLGPUVolumeRayCastMapper.cpp Rendering/vtkMitkOpenGLGPUVolumeRayCastMapperShaders.cpp ) IF(WIN32 AND NOT MINGW) SET(CPP_FILES Controllers/mitkMovieGeneratorWin32.cpp ${CPP_FILES} ) ENDIF(WIN32 AND NOT MINGW) +IF ( ${VTK_MAJOR_VERSION}.${VTK_MINOR_VERSION} VERSION_LESS 5.8 ) + MESSAGE(STATUS "Using VTK 5.8 classes from MITK respository") + SET(CPP_FILES ${CPP_FILES} + Algorithms/vtkImageStencilRaster.cxx + Algorithms/vtkLassoStencilSource.cxx + ) +ENDIF ( ${VTK_MAJOR_VERSION}.${VTK_MINOR_VERSION} VERSION_LESS 5.8 ) + + +