diff --git a/Modules/ContourModel/Algorithms/mitkContourModelUtils.h b/Modules/ContourModel/Algorithms/mitkContourModelUtils.h index ec8b5fe1b9..29ccd2919f 100644 --- a/Modules/ContourModel/Algorithms/mitkContourModelUtils.h +++ b/Modules/ContourModel/Algorithms/mitkContourModelUtils.h @@ -1,155 +1,155 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef mitkContourModelUtils_h #define mitkContourModelUtils_h #include #include #include #include namespace mitk { /** * \brief Helpful methods for working with contours and images * * */ class MITKCONTOURMODEL_EXPORT ContourModelUtils : public itk::Object { public: mitkClassMacroItkParent(ContourModelUtils, itk::Object); /** \brief Projects a contour onto an image point by point. Converts from world to index coordinates. \param slice \param contourIn3D */ static ContourModel::Pointer ProjectContourTo2DSlice(const Image *slice, const ContourModel *contourIn3D); /** \brief Projects a slice index coordinates of a contour back into world coordinates. \param sliceGeometry \param contourIn2D */ static ContourModel::Pointer BackProjectContourFrom2DSlice(const BaseGeometry *sliceGeometry, const ContourModel *contourIn2D); /** \brief Fill a contour in a 2D slice with a specified pixel value. This version always uses the contour of time step 0 and fills the image. \deprecated This function is deprecated. Use FillContourInSlice2() (in conjunction e.g. with TransferLabelContent()) instead. \pre sliceImage points to a valid instance \pre projectedContour points to a valid instance */ //[[deprecated]] DEPRECATED(static void FillContourInSlice(const ContourModel *projectedContour, Image *sliceImage, const Image* workingImage, int paintingPixelValue = 1)); /** \brief Fill a contour in a 2D slice with a specified pixel value. This overloaded version uses the contour at the passed contourTimeStep to fill the passed image slice. \deprecated This function is deprecated. Use FillContourInSlice2() (in conjunction e.g. with TransferLabelContentAtTimeStep()) instead. \pre sliceImage points to a valid instance \pre projectedContour points to a valid instance */ //[[deprecated]] DEPRECATED(static void FillContourInSlice(const ContourModel *projectedContour, TimeStepType contourTimeStep, Image *sliceImage, const Image* workingImage, int paintingPixelValue = 1)); /** \brief Fill a contour in a 2D slice with a specified pixel value. This version always uses the contour of time step 0 and fills the image. \param projectedContour Pointer to the contour that should be projected. \param sliceImage Pointer to the image which content should be altered by adding the contour with the specified paintingPixelValue. \param paintingPixelValue \pre sliceImage points to a valid instance \pre projectedContour points to a valid instance */ static void FillContourInSlice2(const ContourModel* projectedContour, Image* sliceImage, int paintingPixelValue = 1); /** \brief Fill a contour in a 2D slice with a specified pixel value. This overloaded version uses the contour at the passed contourTimeStep to fill the passed image slice. \param projectedContour Pointer to the contour that should be projected. \param contourTimeStep \param sliceImage Pointer to the image which content should be altered by \param paintingPixelValue adding the contour with the specified paintingPixelValue. \pre sliceImage points to a valid instance \pre projectedContour points to a valid instance */ static void FillContourInSlice2(const ContourModel* projectedContour, TimeStepType contourTimeStep, Image* sliceImage, int paintingPixelValue = 1); /** \brief Fills the paintingPixelValue into every pixel of resultImage as indicated by filledImage. If a LableSet image is specified it also by incorporating the rules of LabelSet images when filling the content. - \param filledImage Pointer to the image content that should be checked to decided if a pixel in resultImage should + \param filledImage Pointer to the image content that should be checked to decide if a pixel in resultImage should be filled with paintingPixelValue or not. \param resultImage Pointer to the image content that should be overwritten guided by the content of filledImage. \param image Pointer to an mitk image that allows to define the LabelSet image which states steer the filling process. If an LabelSet instance is passed its states (e.g. locked labels etc...) will be used. If nullptr or an normal image is passed, then simply any pixel position indicated by filledImage will be overwritten. \param paintingPixelValue the pixelvalue/label that should be used in the result image when filling. \param fillForegroundThreshold The threshold value that decides if a pixel in the filled image counts as foreground (>=fillForegroundThreshold) or not. \deprecated This function is deprecated. Use TransferLabelContent() instead. */ [[deprecated]] static void FillSliceInSlice(vtkSmartPointer filledImage, vtkSmartPointer resultImage, const Image* image, int paintingPixelValue, double fillForegroundThreshold = 1.0); /** \brief Move the contour in time step 0 to to a new contour model at the given time step. */ static ContourModel::Pointer MoveZerothContourTimeStep(const ContourModel *contour, TimeStepType timeStep); /** \brief Retrieves the active pixel value of a (labelset) image. If the image is basic image, the pixel value 1 (one) will be returned. If the image is actually a labelset image, the pixel value of the active label of the active layer will be returned. \param workingImage The (labelset) image to retrieve the active pixel value of. */ static int GetActivePixelValue(const Image* workingImage); protected: ContourModelUtils(); ~ContourModelUtils() override; }; } #endif diff --git a/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp b/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp index 37fd7ad7a7..09a04be161 100644 --- a/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp +++ b/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp @@ -1,542 +1,542 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include #include #include #include #include #include "mitkImageAccessByItk.h" #include #include #include #include namespace mitk { HotspotMaskGenerator::HotspotMaskGenerator(): m_HotspotRadiusInMM(6.2035049089940), // radius of a 1cm3 sphere in mm m_HotspotMustBeCompletelyInsideImage(true), m_Label(1) { m_InternalMask = mitk::Image::New(); m_InternalMaskUpdateTime = 0; } HotspotMaskGenerator::~HotspotMaskGenerator() { } unsigned int HotspotMaskGenerator::GetNumberOfMasks() const { return 1; } mitk::Image::ConstPointer HotspotMaskGenerator::DoGetMask(unsigned int) { if (IsUpdateRequired()) { if ( m_InputImage.IsNull() ) { throw std::runtime_error( "Error: image empty!" ); } if ( m_InputImage->GetTimeGeometry()->IsValidTimePoint(m_TimePoint) ) { throw std::runtime_error( "Error: invalid time point!" ); } auto timeSliceImage = SelectImageByTimePoint(m_InputImage, m_TimePoint); m_internalMask2D = nullptr; // is this correct when this variable holds a smart pointer? m_internalMask3D = nullptr; if ( m_Mask != nullptr ) { m_Mask->SetTimePoint(m_TimePoint); mitk::Image::ConstPointer timeSliceMask = m_Mask->GetMask(0); if ( timeSliceImage->GetDimension() == 3 ) { itk::Image::Pointer noneConstMaskImage; //needed to work around the fact that CastToItkImage currently does not support const itk images. CastToItkImage(timeSliceMask, noneConstMaskImage); m_internalMask3D = noneConstMaskImage; AccessFixedDimensionByItk_2(timeSliceImage, CalculateHotspotMask, 3, m_internalMask3D.GetPointer(), m_Label); } else if ( timeSliceImage->GetDimension() == 2 ) { itk::Image::Pointer noneConstMaskImage; //needed to work around the fact that CastToItkImage currently does not support const itk images. CastToItkImage(timeSliceMask, noneConstMaskImage); m_internalMask2D = noneConstMaskImage; AccessFixedDimensionByItk_2(timeSliceImage, CalculateHotspotMask, 2, m_internalMask2D.GetPointer(), m_Label); } else { throw std::runtime_error( "Error: invalid image dimension" ); } } else { if ( timeSliceImage->GetDimension() == 3 ) { AccessFixedDimensionByItk_2(timeSliceImage, CalculateHotspotMask, 3, m_internalMask3D.GetPointer(), m_Label); } else if ( timeSliceImage->GetDimension() == 2 ) { AccessFixedDimensionByItk_2(timeSliceImage, CalculateHotspotMask, 2, m_internalMask2D.GetPointer(), m_Label); } else { throw std::runtime_error( "Error: invalid image dimension" ); } } this->Modified(); } m_InternalMaskUpdateTime = m_InternalMask->GetMTime(); return m_InternalMask; } template HotspotMaskGenerator::ImageExtrema HotspotMaskGenerator::CalculateExtremaWorld( const itk::Image* inputImage, const itk::Image* maskImage, - double neccessaryDistanceToImageBorderInMM, + double necessaryDistanceToImageBorderInMM, unsigned int label ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef itk::ImageRegionConstIteratorWithIndex MaskImageIteratorType; typedef itk::ImageRegionConstIteratorWithIndex InputImageIndexIteratorType; typename ImageType::SpacingType spacing = inputImage->GetSpacing(); ImageExtrema minMax; minMax.Defined = false; minMax.MaxIndex.set_size(VImageDimension); minMax.MaxIndex.set_size(VImageDimension); typename ImageType::RegionType allowedExtremaRegion = inputImage->GetLargestPossibleRegion(); - bool keepDistanceToImageBorders( neccessaryDistanceToImageBorderInMM > 0 ); + bool keepDistanceToImageBorders( necessaryDistanceToImageBorderInMM > 0 ); if (keepDistanceToImageBorders) { itk::IndexValueType distanceInPixels[VImageDimension]; for(unsigned short dimension = 0; dimension < VImageDimension; ++dimension) { // To confirm that the whole hotspot is inside the image we have to keep a specific distance to the image-borders, which is as long as // the radius. To get the amount of indices we divide the radius by spacing and add 0.5 because voxels are center based: // For example with a radius of 2.2 and a spacing of 1 two indices are enough because 2.2 / 1 + 0.5 = 2.7 => 2. // But with a radius of 2.7 we need 3 indices because 2.7 / 1 + 0.5 = 3.2 => 3 - distanceInPixels[dimension] = int( neccessaryDistanceToImageBorderInMM / spacing[dimension] + 0.5); + distanceInPixels[dimension] = int( necessaryDistanceToImageBorderInMM / spacing[dimension] + 0.5); } allowedExtremaRegion.ShrinkByRadius(distanceInPixels); } InputImageIndexIteratorType imageIndexIt(inputImage, allowedExtremaRegion); float maxValue = itk::NumericTraits::min(); float minValue = itk::NumericTraits::max(); typename ImageType::IndexType maxIndex; typename ImageType::IndexType minIndex; for(unsigned short i = 0; i < VImageDimension; ++i) { maxIndex[i] = 0; minIndex[i] = 0; } if (maskImage != nullptr) { MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); typename ImageType::IndexType imageIndex; typename ImageType::IndexType maskIndex; for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { imageIndex = maskIndex = maskIt.GetIndex(); if(maskIt.Get() == label) { if( allowedExtremaRegion.IsInside(imageIndex) ) { imageIndexIt.SetIndex( imageIndex ); double value = imageIndexIt.Get(); minMax.Defined = true; //Calculate minimum, maximum and corresponding index-values if( value > maxValue ) { maxIndex = imageIndexIt.GetIndex(); maxValue = value; } if(value < minValue ) { minIndex = imageIndexIt.GetIndex(); minValue = value; } } } } } else { for(imageIndexIt.GoToBegin(); !imageIndexIt.IsAtEnd(); ++imageIndexIt) { double value = imageIndexIt.Get(); minMax.Defined = true; //Calculate minimum, maximum and corresponding index-values if( value > maxValue ) { maxIndex = imageIndexIt.GetIndex(); maxValue = value; } if(value < minValue ) { minIndex = imageIndexIt.GetIndex(); minValue = value; } } } minMax.MaxIndex.set_size(VImageDimension); minMax.MinIndex.set_size(VImageDimension); for(unsigned int i = 0; i < minMax.MaxIndex.size(); ++i) { minMax.MaxIndex[i] = maxIndex[i]; } for(unsigned int i = 0; i < minMax.MinIndex.size(); ++i) { minMax.MinIndex[i] = minIndex[i]; } minMax.Max = maxValue; minMax.Min = minValue; return minMax; } template itk::Size HotspotMaskGenerator::CalculateConvolutionKernelSize( double spacing[VImageDimension], double radiusInMM ) { typedef itk::Image< float, VImageDimension > KernelImageType; typedef typename KernelImageType::SizeType SizeType; SizeType maskSize; for(unsigned int i = 0; i < VImageDimension; ++i) { maskSize[i] = static_cast( 2 * radiusInMM / spacing[i]); // We always want an uneven size to have a clear center point in the convolution mask if(maskSize[i] % 2 == 0 ) { ++maskSize[i]; } } return maskSize; } template itk::SmartPointer< itk::Image > HotspotMaskGenerator::GenerateHotspotSearchConvolutionKernel(double mmPerPixel[VImageDimension], double radiusInMM ) { std::stringstream ss; for (unsigned int i = 0; i < VImageDimension; ++i) { ss << mmPerPixel[i]; if (i < VImageDimension -1) ss << ","; } MITK_DEBUG << "Update convolution kernel for spacing (" << ss.str() << ") and radius " << radiusInMM << "mm"; double radiusInMMSquared = radiusInMM * radiusInMM; typedef itk::Image< float, VImageDimension > KernelImageType; typename KernelImageType::Pointer convolutionKernel = KernelImageType::New(); // Calculate size and allocate mask image typedef typename KernelImageType::SizeType SizeType; SizeType maskSize = this->CalculateConvolutionKernelSize(mmPerPixel, radiusInMM); mitk::Point3D convolutionMaskCenterIndex; convolutionMaskCenterIndex.Fill(0.0); for(unsigned int i = 0; i < VImageDimension; ++i) { convolutionMaskCenterIndex[i] = 0.5 * (double)(maskSize[i]-1); } typedef typename KernelImageType::IndexType IndexType; IndexType maskIndex; maskIndex.Fill(0); typedef typename KernelImageType::RegionType RegionType; RegionType maskRegion; maskRegion.SetSize(maskSize); maskRegion.SetIndex(maskIndex); convolutionKernel->SetRegions(maskRegion); convolutionKernel->SetSpacing(mmPerPixel); convolutionKernel->Allocate(); // Fill mask image values by subsampling the image grid typedef itk::ImageRegionIteratorWithIndex MaskIteratorType; MaskIteratorType maskIt(convolutionKernel,maskRegion); int numberOfSubVoxelsPerDimension = 2; // per dimension! int numberOfSubVoxels = ::pow( static_cast(numberOfSubVoxelsPerDimension), static_cast(VImageDimension) ); double subVoxelSizeInPixels = 1.0 / (double)numberOfSubVoxelsPerDimension; double valueOfOneSubVoxel = 1.0 / (double)numberOfSubVoxels; mitk::Point3D subVoxelIndexPosition; double distanceSquared = 0.0; typedef itk::ContinuousIndex ContinuousIndexType; for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { ContinuousIndexType indexPoint(maskIt.GetIndex()); mitk::Point3D voxelPosition; for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) { voxelPosition[dimension] = indexPoint[dimension]; } double maskValue = 0.0; mitk::Vector3D subVoxelOffset; subVoxelOffset.Fill(0.0); // iterate sub-voxels by iterating all possible offsets for (subVoxelOffset[0] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[0] < +0.5; subVoxelOffset[0] += subVoxelSizeInPixels) { for (subVoxelOffset[1] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[1] < +0.5; subVoxelOffset[1] += subVoxelSizeInPixels) { for (subVoxelOffset[2] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[2] < +0.5; subVoxelOffset[2] += subVoxelSizeInPixels) { subVoxelIndexPosition = voxelPosition + subVoxelOffset; // this COULD be integrated into the for-loops if necessary (add voxelPosition to initializer and end condition) distanceSquared = (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] * (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] + (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] * (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] + (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2] * (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2]; if (distanceSquared <= radiusInMMSquared) { maskValue += valueOfOneSubVoxel; } } } } maskIt.Set( maskValue ); } return convolutionKernel; } template itk::SmartPointer > HotspotMaskGenerator::GenerateConvolutionImage( const itk::Image* inputImage ) { double mmPerPixel[VImageDimension]; for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) { mmPerPixel[dimension] = inputImage->GetSpacing()[dimension]; } // update convolution kernel typedef itk::Image< float, VImageDimension > KernelImageType; typename KernelImageType::Pointer convolutionKernel = this->GenerateHotspotSearchConvolutionKernel(mmPerPixel, m_HotspotRadiusInMM); // update convolution image typedef itk::Image< TPixel, VImageDimension > InputImageType; typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; typedef itk::FFTConvolutionImageFilter ConvolutionFilterType; itk::VnlFFTImageFilterInitFactory::RegisterFactories(); typename ConvolutionFilterType::Pointer convolutionFilter = ConvolutionFilterType::New(); typedef itk::ConstantBoundaryCondition BoundaryConditionType; BoundaryConditionType boundaryCondition; boundaryCondition.SetConstant(0.0); if (m_HotspotMustBeCompletelyInsideImage) { // overwrite default boundary condition convolutionFilter->SetBoundaryCondition(&boundaryCondition); } convolutionFilter->SetInput(inputImage); convolutionFilter->SetKernelImage(convolutionKernel); convolutionFilter->SetNormalize(true); MITK_DEBUG << "Update Convolution image for hotspot search"; convolutionFilter->UpdateLargestPossibleRegion(); typename ConvolutionImageType::Pointer convolutionImage = convolutionFilter->GetOutput(); convolutionImage->SetSpacing( inputImage->GetSpacing() ); // only workaround because convolution filter seems to ignore spacing of input image return convolutionImage; } template < typename TPixel, unsigned int VImageDimension> void HotspotMaskGenerator::FillHotspotMaskPixels( itk::Image* maskImage, itk::Point sphereCenter, double sphereRadiusInMM ) { typedef itk::Image< TPixel, VImageDimension > MaskImageType; typedef itk::ImageRegionIteratorWithIndex MaskImageIteratorType; MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); typename MaskImageType::IndexType maskIndex; typename MaskImageType::PointType worldPosition; // this is not very smart. I would rather use a 0 initialized mask (not the case here -> blame CalculateHotspotMask) and find the region where I need to iterate over, then iterate only over the small region for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { maskIndex = maskIt.GetIndex(); maskImage->TransformIndexToPhysicalPoint(maskIndex, worldPosition); maskIt.Set( worldPosition.EuclideanDistanceTo(sphereCenter) <= sphereRadiusInMM ? 1 : 0 ); } } template void HotspotMaskGenerator::CalculateHotspotMask(const itk::Image* inputImage, const itk::Image* maskImage, unsigned int label) { typedef itk::Image< TPixel, VImageDimension > InputImageType; typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typename ConvolutionImageType::Pointer convolutionImage = this->GenerateConvolutionImage(inputImage); if (convolutionImage.IsNull()) { MITK_ERROR << "Empty convolution image in CalculateHotspotStatistics(). We should never reach this state (logic error)."; throw std::logic_error("Empty convolution image in CalculateHotspotStatistics()"); } typename MaskImageType::ConstPointer usedMask = maskImage; // if mask image is not defined, create an image of the same size as inputImage and fill it with 1's // there is maybe a better way to do this!? if (maskImage == nullptr) { auto defaultMask = MaskImageType::New(); typename MaskImageType::RegionType maskRegion = inputImage->GetLargestPossibleRegion(); typename MaskImageType::SpacingType maskSpacing = inputImage->GetSpacing(); typename MaskImageType::PointType maskOrigin = inputImage->GetOrigin(); typename MaskImageType::DirectionType maskDirection = inputImage->GetDirection(); defaultMask->SetRegions(maskRegion); defaultMask->Allocate(); defaultMask->SetOrigin(maskOrigin); defaultMask->SetSpacing(maskSpacing); defaultMask->SetDirection(maskDirection); defaultMask->FillBuffer(1); usedMask = defaultMask; label = 1; } // find maximum in convolution image, given the current mask double requiredDistanceToBorder = m_HotspotMustBeCompletelyInsideImage ? m_HotspotRadiusInMM : -1.0; ImageExtrema convolutionImageInformation = CalculateExtremaWorld(convolutionImage.GetPointer(), usedMask.GetPointer(), requiredDistanceToBorder, label); bool isHotspotDefined = convolutionImageInformation.Defined; if (!isHotspotDefined) { MITK_ERROR << "No origin of hotspot-sphere was calculated!"; m_InternalMask = nullptr; } else { // create a binary mask around the "hotspot" region, fill the shape of a sphere around our hotspot center // typename DuplicatorType::Pointer copyMachine = DuplicatorType::New(); // copyMachine->SetInputImage(inputImage); // copyMachine->Update(); // typename CastFilterType::Pointer caster = CastFilterType::New(); // caster->SetInput( copyMachine->GetOutput() ); // caster->Update(); typename MaskImageType::Pointer hotspotMaskITK = MaskImageType::New(); hotspotMaskITK->SetOrigin(inputImage->GetOrigin()); hotspotMaskITK->SetSpacing(inputImage->GetSpacing()); hotspotMaskITK->SetLargestPossibleRegion(inputImage->GetLargestPossibleRegion()); hotspotMaskITK->SetBufferedRegion(inputImage->GetBufferedRegion()); hotspotMaskITK->SetDirection(inputImage->GetDirection()); hotspotMaskITK->SetNumberOfComponentsPerPixel(inputImage->GetNumberOfComponentsPerPixel()); hotspotMaskITK->Allocate(); hotspotMaskITK->FillBuffer(1); typedef typename InputImageType::IndexType IndexType; IndexType maskCenterIndex; for (unsigned int d =0; d< VImageDimension;++d) { maskCenterIndex[d]=convolutionImageInformation.MaxIndex[d]; } typename ConvolutionImageType::PointType maskCenter; inputImage->TransformIndexToPhysicalPoint(maskCenterIndex,maskCenter); FillHotspotMaskPixels(hotspotMaskITK.GetPointer(), maskCenter, m_HotspotRadiusInMM); //obtain mitk::Image::Pointer from itk::Image mitk::Image::Pointer hotspotMaskAsMITKImage = mitk::GrabItkImageMemory(hotspotMaskITK); m_InternalMask = hotspotMaskAsMITKImage; m_ConvolutionImageMaxIndex = convolutionImageInformation.MaxIndex; m_ConvolutionImageMinIndex = convolutionImageInformation.MinIndex; } } bool HotspotMaskGenerator::IsUpdateRequired() const { unsigned long thisClassTimeStamp = this->GetMTime(); unsigned long internalMaskTimeStamp = m_InternalMask->GetMTime(); unsigned long maskGeneratorTimeStamp = m_Mask->GetMTime(); unsigned long inputImageTimeStamp = m_InputImage->GetMTime(); if (thisClassTimeStamp > m_InternalMaskUpdateTime) // inputs have changed { return true; } if (m_InternalMaskUpdateTime < maskGeneratorTimeStamp || m_InternalMaskUpdateTime < inputImageTimeStamp) // mask image has changed outside of this class { return true; } if (internalMaskTimeStamp > m_InternalMaskUpdateTime) // internal mask has been changed outside of this class { return true; } return false; } } diff --git a/Modules/Segmentation/Algorithms/itkContourExtractor2DImageFilter.txx b/Modules/Segmentation/Algorithms/itkContourExtractor2DImageFilter.txx index 36428ed91a..852c3ad118 100644 --- a/Modules/Segmentation/Algorithms/itkContourExtractor2DImageFilter.txx +++ b/Modules/Segmentation/Algorithms/itkContourExtractor2DImageFilter.txx @@ -1,536 +1,536 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ /*=================================================================== This file is based heavily on a corresponding ITK filter. ===================================================================*/ #ifndef __itkContourExtractor2DImageFilter_txx #define __itkContourExtractor2DImageFilter_txx #include "itkConstShapedNeighborhoodIterator.h" #include "itkConstShapedNeighborhoodIterator.h" #include "itkContourExtractor2DImageFilter.h" #include "itkProgressReporter.h" #include namespace itk { // Constructor template ContourExtractor2DImageFilter::ContourExtractor2DImageFilter() { this->m_ContourValue = NumericTraits::Zero; this->m_ReverseContourOrientation = false; this->m_VertexConnectHighPixels = false; this->m_UseCustomRegion = false; this->m_NumberOfContoursCreated = 0; } // Destructor template ContourExtractor2DImageFilter::~ContourExtractor2DImageFilter() { } template void ContourExtractor2DImageFilter::GenerateData() { // Make sure the structures for containing, looking up, and numbering the // growing contours are empty and ready. m_Contours.clear(); m_ContourStarts.clear(); m_ContourEnds.clear(); m_NumberOfContoursCreated = 0; // Set up an iterator to "march the squares" across the image. // We associate each 2px-by-2px square with the pixel in the upper left of // that square. We then iterate across the image, examining these 2x2 squares // and building the contour. By iterating the upper-left pixel of our // "current square" across every pixel in the image except those on the // bottom row and rightmost column, we have visited every valid square in the // image. InputRegionType region = this->GetInput()->GetRequestedRegion(); typename InputRegionType::SizeType shrunkSize = region.GetSize(); shrunkSize[0] -= 1; shrunkSize[1] -= 1; InputRegionType shrunkRegion(region.GetIndex(), shrunkSize); // Set up a progress reporter ProgressReporter progress(this, 0, shrunkRegion.GetNumberOfPixels()); // A 1-pixel radius sets up a neighborhood with the following indices: // 0 1 2 // 3 4 5 // 6 7 8 // We are interested only in the square of 4,5,7,8 which is the 2x2 square // with the center pixel at the top-left. So we only activate the // corresponding offsets, and only query pixels 4, 5, 7, and 8 with the // iterator's GetPixel method. typedef ConstShapedNeighborhoodIterator SquareIterator; typename SquareIterator::RadiusType radius = {{1, 1}}; SquareIterator it(radius, this->GetInput(), shrunkRegion); InputOffsetType none = {{0, 0}}; InputOffsetType right = {{1, 0}}; InputOffsetType down = {{0, 1}}; InputOffsetType diag = {{1, 1}}; it.ActivateOffset(none); it.ActivateOffset(right); it.ActivateOffset(down); it.ActivateOffset(diag); for (it.GoToBegin(); !it.IsAtEnd(); ++it) { // There are sixteen different possible square types, diagramed below. // A + indicates that the vertex is above the contour value, and a - // indicates that the vertex is below or equal to the contour value. // The vertices of each square are here numbered: // 01 // 23 // and treated as a binary value with the bits in that order. Thus each // square can be so numbered: // 0-- 1+- 2-+ 3++ 4-- 5+- 6-+ 7++ // -- -- -- -- +- +- +- +- // // 8-- 9+- 10-+ 11++ 12-- 13+- 14-+ 15++ // -+ -+ -+ -+ ++ ++ ++ ++ // // The position of the line segment that cuts through (or doesn't, in case // 0 and 15) each square is clear, except in cases 6 and 9. In this case, // where the segments are placed is determined by // m_VertexConnectHighPixels. If m_VertexConnectHighPixels is false, then // lines like are drawn through square 6, and lines like are drawn through // square 9. Otherwise, the situation is reversed. // Finally, recall that we draw the lines so that (moving from tail to // head) the lower-valued pixels are on the left of the line. So, for // example, case 1 entails a line slanting from the middle of the top of // the square to the middle of the left side of the square. // (1) Determine what number square we are currently inspecting. Remember // that as far as the neighborhood iterator is concerned, our square // 01 is numbered as 45 // 23 78 InputPixelType v0, v1, v2, v3; v0 = it.GetPixel(4); v1 = it.GetPixel(5); v2 = it.GetPixel(7); v3 = it.GetPixel(8); InputIndexType index = it.GetIndex(); unsigned char squareCase = 0; if (v0 > m_ContourValue) squareCase += 1; if (v1 > m_ContourValue) squareCase += 2; if (v2 > m_ContourValue) squareCase += 4; if (v3 > m_ContourValue) squareCase += 8; // Set up macros to find the ContinuousIndex where the contour intersects // one of the sides of the square. Normally macros should, of course, be // eschewed, but since this is an inner loop not calling the function four // times when two would do is probably worth while. Plus, copy-pasting // these into the switch below is even worse. InterpolateContourPosition // takes the values at two vertices, the index of the first vertex, and the // offset between the two vertices. #define TOP_ this->InterpolateContourPosition(v0, v1, index, right) #define BOTTOM_ this->InterpolateContourPosition(v2, v3, index + down, right) #define LEFT_ this->InterpolateContourPosition(v0, v2, index, down) #define RIGHT_ this->InterpolateContourPosition(v1, v3, index + right, down) // (2) Add line segments to the growing contours as defined by the cases. // AddSegment takes a "from" vertex and a "to" vertex, and adds it to the // a growing contour, creates a new contour, or merges two together. switch (squareCase) { case 0: // no line break; case 1: // top to left this->AddSegment(TOP_, LEFT_); break; case 2: // right to top this->AddSegment(RIGHT_, TOP_); break; case 3: // right to left this->AddSegment(RIGHT_, LEFT_); break; case 4: // left to bottom this->AddSegment(LEFT_, BOTTOM_); break; case 5: // top to bottom this->AddSegment(TOP_, BOTTOM_); break; case 6: if (m_VertexConnectHighPixels) { // left to top this->AddSegment(LEFT_, TOP_); // right to bottom this->AddSegment(RIGHT_, BOTTOM_); } else { // right to top this->AddSegment(RIGHT_, TOP_); // left to bottom this->AddSegment(LEFT_, BOTTOM_); } break; case 7: // right to bottom this->AddSegment(RIGHT_, BOTTOM_); break; case 8: // bottom to right this->AddSegment(BOTTOM_, RIGHT_); break; case 9: if (m_VertexConnectHighPixels) { // top to right this->AddSegment(TOP_, RIGHT_); // bottom to left this->AddSegment(BOTTOM_, LEFT_); } else { // top to left this->AddSegment(TOP_, LEFT_); // bottom to right this->AddSegment(BOTTOM_, RIGHT_); } break; case 10: // bottom to top this->AddSegment(BOTTOM_, TOP_); break; case 11: // bottom to left this->AddSegment(BOTTOM_, LEFT_); break; case 12: // left to right this->AddSegment(LEFT_, RIGHT_); break; case 13: // top to right this->AddSegment(TOP_, RIGHT_); break; case 14: // left to top this->AddSegment(LEFT_, TOP_); break; case 15: // no line break; } // switch squareCase progress.CompletedPixel(); } // iteration - // Now create the outputs paths from the dequeues we've been using. + // Now create the outputs paths from the deques we've been using. this->FillOutputs(); m_Contours.clear(); m_ContourStarts.clear(); m_ContourEnds.clear(); m_NumberOfContoursCreated = 0; } template inline typename ContourExtractor2DImageFilter::VertexType ContourExtractor2DImageFilter::InterpolateContourPosition(InputPixelType fromValue, InputPixelType toValue, InputIndexType fromIndex, InputOffsetType toOffset) { VertexType output; // Now calculate the fraction of the way from 'from' to 'to' that the contour // crosses. Interpolate linearly: y = v0 + (v1 - v0) * x, and solve for the // x that gives y = m_ContourValue: x = (m_ContourValue - v0) / (v1 - v0). // This assumes that v0 and v1 are separated by exactly ONE unit. So the to // Offset. value must have exactly one component 1 and the other component 0. // Also this assumes that fromValue and toValue are different. Otherwise we // can't interpolate anything! itkAssertOrThrowMacro((fromValue != toValue), "source and destination are the same"); itkAssertOrThrowMacro(((toOffset[0] == 0 && toOffset[1] == 1) || (toOffset[0] == 1 && toOffset[1] == 0)), "toOffset has unexpected values"); double x = (m_ContourValue - static_cast(fromValue)) / (toValue - static_cast(fromValue)); output[0] = fromIndex[0] + x * toOffset[0]; output[1] = fromIndex[1] + x * toOffset[1]; return output; } template void ContourExtractor2DImageFilter::AddSegment(VertexType from, VertexType to) { if (from == to) { // Arc is degenerate: ignore, and the from/two point will be connected // later by other squares. Degeneracy happens when (and only when) a square // has exactly one vertex that is the contour value, and the rest are above // that value. return; } // Try to find an existing contour that starts where the new segment ends. VertexMapIterator newTail = m_ContourStarts.find(to); // Try to find an existing contour that ends where the new segment starts. VertexMapIterator newHead = m_ContourEnds.find(from); if (newTail != m_ContourStarts.end() && newHead != m_ContourEnds.end()) { // We need to connect these two contours. The act of connecting them will // add the needed arc. auto tail = newTail->second; itkAssertOrThrowMacro((tail->front() == to), "End doesn't match Beginning"); auto head = newHead->second; itkAssertOrThrowMacro((head->back() == from), "Beginning doesn't match End"); if (head == tail) { // We've closed a contour. Add the end point, and remove from the maps head->push_back(to); m_ContourStarts.erase(newTail); // erase the front of tail. Because head and tail are the same contour, // don't worry about erasing the front of head! m_ContourEnds.erase(newHead); // erase the end of head/tail. } else { // We have found two distinct contours that need to be joined. Careful // here: we want to keep the first segment in the list when merging so // that contours are always returned in top-to-bottom, right-to-left // order (with regard to the image pixel found to be inside the contour). if (tail->m_ContourNumber > head->m_ContourNumber) { // if tail was created later than head... // Copy tail to the end of head and remove // tail from everything. head->insert(head->end(), tail->begin(), tail->end()); // Now remove 'tail' from the list and the maps because it has been // subsumed. m_ContourStarts.erase(newTail); int erased = m_ContourEnds.erase(tail->back()); // There should be exactly one entry in the hash for that endpoint if (erased != 1) { itkWarningMacro(<< "There should be exactly one entry in the hash for that endpoint, but there are " << erased); } m_Contours.erase(tail); // remove from the master list // Now remove the old end of 'head' from the ends map and add // the new end. m_ContourEnds.erase(newHead); m_ContourEnds.insert(VertexContourRefPair(head->back(), head)); } else { // Copy head to the beginning of tail and remove // head from everything. tail->insert(tail->begin(), head->begin(), head->end()); // Now remove 'head' from the list and the maps because // it has been subsumed. m_ContourEnds.erase(newHead); int erased = m_ContourStarts.erase(head->front()); if (erased != 1) { itkWarningMacro(<< "There should be exactly one entry in the hash for that endpoint, but there are " << erased); } m_Contours.erase(head); // remove from the master list // Now remove the old start of 'tail' from the starts map and // add the new start. m_ContourStarts.erase(newTail); m_ContourStarts.insert(VertexContourRefPair(tail->front(), tail)); } } } else if (newTail == m_ContourStarts.end() && newHead == m_ContourEnds.end()) { // No contours found: add a new one. // Make it on the heap. It will be copied into m_Contours. ContourType contour; // Add the endpoints contour.push_front(from); contour.push_back(to); contour.m_ContourNumber = m_NumberOfContoursCreated++; // Add the contour to the end of the list and get a reference to it. m_Contours.push_back(contour); // recall that end() is an iterator to one past the back! auto newContour = --m_Contours.end(); // add the endpoints and an iterator pointing to the contour // in the list to the maps. m_ContourStarts.insert(VertexContourRefPair(from, newContour)); m_ContourEnds.insert(VertexContourRefPair(to, newContour)); } else if (newTail != m_ContourStarts.end() && newHead == m_ContourEnds.end()) { // Found a single contour to which the new arc should be prepended. auto tail = newTail->second; itkAssertOrThrowMacro((tail->front() == to), "End doesn't match Beginning"); tail->push_front(from); // erase the old start of this contour m_ContourStarts.erase(newTail); // Now add the new start of this contour. m_ContourStarts.insert(VertexContourRefPair(from, tail)); } else if (newTail == m_ContourStarts.end() && newHead != m_ContourEnds.end()) { // Found a single contour to which the new arc should be appended. auto head = newHead->second; itkAssertOrThrowMacro((head->back() == from), "Beginning doesn't match End"); head->push_back(to); // erase the old end of this contour m_ContourEnds.erase(newHead); // Now add the new start of this contour. m_ContourEnds.insert(VertexContourRefPair(to, head)); } } template void ContourExtractor2DImageFilter::FillOutputs() { this->SetNumberOfIndexedOutputs(m_Contours.size()); int i = 0; for (auto it = m_Contours.begin(); it != m_Contours.end(); it++, i++) { OutputPathPointer output = this->GetOutput(i); if (output.IsNull()) { // Static cast is OK because we know PathSource will make its templated // class type output = static_cast(this->MakeOutput(i).GetPointer()); this->SetNthOutput(i, output.GetPointer()); } typename VertexListType::Pointer path = const_cast(output->GetVertexList()); path->Initialize(); path->reserve(it->size()); // use std::vector version of 'reserve()' // instead of VectorContainer::Reserve() to work around // the fact that the latter is essentially std::vector::resize(), // which is not what we want. // Now put all the points from the contour deque into the path and // mark output as modified typedef typename ContourType::const_iterator ConstIteratorType; if (m_ReverseContourOrientation) { ConstIteratorType itC = (*it).end(); do { itC--; path->push_back(*itC); } while (itC != (*it).begin()); } else { ConstIteratorType itC = (*it).begin(); while (itC != (*it).end()) { path->push_back(*itC); itC++; } } output->Modified(); } } template void ContourExtractor2DImageFilter::SetRequestedRegion(const InputRegionType region) { itkDebugMacro("setting RequestedRegion to " << region); m_UseCustomRegion = true; if (this->m_RequestedRegion != region) { this->m_RequestedRegion = region; this->Modified(); } } template void ContourExtractor2DImageFilter::ClearRequestedRegion() { itkDebugMacro("Clearing RequestedRegion."); if (this->m_UseCustomRegion == true) { this->m_UseCustomRegion = false; this->Modified(); } } template void ContourExtractor2DImageFilter::GenerateInputRequestedRegion() { InputImageType *input = const_cast(this->GetInput()); if (!input) return; if (m_UseCustomRegion) { InputRegionType requestedRegion = m_RequestedRegion; if (requestedRegion.Crop(input->GetLargestPossibleRegion())) { input->SetRequestedRegion(requestedRegion); return; } else { // Couldn't crop the region (requested region is outside the largest // possible region). Throw an exception. // store what we tried to request (prior to trying to crop) input->SetRequestedRegion(requestedRegion); // build an exception InvalidRequestedRegionError e(__FILE__, __LINE__); e.SetLocation(ITK_LOCATION); e.SetDescription("Requested region is outside the largest possible region."); e.SetDataObject(input); throw e; } } else { input->SetRequestedRegion(input->GetLargestPossibleRegion()); } } /** * Standard "PrintSelf" method */ template void ContourExtractor2DImageFilter::PrintSelf(std::ostream &os, Indent indent) const { Superclass::PrintSelf(os, indent); os << indent << "ReverseContourOrientation: " << m_ReverseContourOrientation << std::endl; os << indent << "VertexConnectHighPixels: " << m_VertexConnectHighPixels << std::endl; os << indent << "UseCustomRegion: " << m_UseCustomRegion << std::endl; os << indent << "NumericTraits: " << m_UseCustomRegion << std::endl; os << indent << "NumberOfContoursCreated: " << m_NumberOfContoursCreated << std::endl; if (m_UseCustomRegion) { os << indent << "Custom region: " << m_RequestedRegion << std::endl; } typedef typename NumericTraits::PrintType InputRealPrintType; os << indent << "Contour value: " << static_cast(m_ContourValue) << std::endl; } } // end namespace itk #endif