diff --git a/Modules/Segmentation/Algorithms/itkAdaptiveThresholdIterator.h b/Modules/Segmentation/Algorithms/itkAdaptiveThresholdIterator.h index e8e6e80175..5cdd9e93e7 100644 --- a/Modules/Segmentation/Algorithms/itkAdaptiveThresholdIterator.h +++ b/Modules/Segmentation/Algorithms/itkAdaptiveThresholdIterator.h @@ -1,245 +1,254 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __itkAdaptiveThresholdIterator_h #define __itkAdaptiveThresholdIterator_h #include "itkConditionalConstIterator.h" #include "itkImage.h" #include "itkIndex.h" #include "itkSize.h" #include #include #include #include namespace itk { /** * \class AdaptiveThresholdIterator * \brief Iterates over an image using a variable image function, * which threshold can be varied during the iteration process. * * \ingroup ImageIterators * */ template class ITK_EXPORT AdaptiveThresholdIterator : public ConditionalConstIterator { public: /** Standard class typedefs. */ typedef AdaptiveThresholdIterator Self; typedef ConditionalConstIterator Superclass; typedef TImage ImageType; // A temporary image used for storing info about all indices typedef TImage TTempImage; typename TTempImage::Pointer tempPtr; //[!] isn't really used?! /** Type of function */ typedef TFunction FunctionType; /** Type of vector used to store location info in the spatial function */ typedef typename TFunction::InputType FunctionInputType; /** Size typedef support. */ typedef typename TImage::SizeType SizeType; /** Region typedef support */ typedef typename TImage::RegionType RegionType; typedef typename TImage::IndexType IndexType; /** Internal Pixel Type */ typedef typename TImage::InternalPixelType InternalPixelType; /** External Pixel Type */ typedef typename TImage::PixelType PixelType; /** Queue containing indices representing a voxel position */ typedef std::queue IndexQueueType; /** Map used to generate the output result */ typedef std::map QueueMapType; /** Dimension of the image the iterator walks. This constant is needed so * that functions that are templated over image iterator type (as opposed to * being templated over pixel type and dimension) can have compile time * access to the dimension of the image that the iterator walks. */ itkStaticConstMacro(NDimensions, unsigned int, TImage::ImageDimension); /** Constructor establishes an iterator to walk a particular image and a * particular region of that image. This version of the constructor uses * an explicit seed pixel for the flood fill, the "startIndex" */ AdaptiveThresholdIterator(ImageType *imagePtr, FunctionType *fnPtr, IndexType startIndex); /** Constructor establishes an iterator to walk a particular image and a * particular region of that image. This version of the constructor uses * an explicit list of seed pixels for the flood fill, the "startIndex" */ AdaptiveThresholdIterator(ImageType *imagePtr, FunctionType *fnPtr, std::vector &startIndex); /** Constructor establishes an iterator to walk a particular image and a * particular region of that image. This version of the constructor * should be used when the seed pixel is unknown */ AdaptiveThresholdIterator(ImageType *imagePtr, FunctionType *fnPtr); /** Default Destructor. */ ~AdaptiveThresholdIterator() override{}; /** Initializes the iterator, called from constructor */ void InitializeIterator(); // makes the iterator go one step further void DoExtendedFloodStep(); // set-method for member-variable void SetExpansionDirection(bool upwards); // Init-method void InitRegionGrowingState(); void SetMinTH(int min); void SetMaxTH(int max); int GetSeedPointValue(void); /** switch between fine and raw leakage detection */ void SetFineDetectionMode(bool fine = false) { m_FineDetectionMode = fine; m_DetectionStop = false; } /** Get the index. This provides a read only reference to the index. * This causes the index to be calculated from pointer arithmetic and is * therefore an expensive operation. * \sa SetIndex */ const IndexType GetIndex() override { return (*m_QueueMap.find(m_RegionGrowingState)).second.front(); } // [!] is never called?! const PixelType Get(void) const override { return const_cast(this->m_Image.GetPointer()) ->GetPixel((*m_QueueMap.find(m_RegionGrowingState)).second.front()); } //[!] is never called?! void Set(const PixelType &value) { const_cast(this->m_Image.GetPointer()) ->GetPixel((*m_QueueMap.find(m_RegionGrowingState)).second.front()) = value; } void GoToBegin(); /** Is the iterator at the end of the region? */ bool IsAtEnd() override { return this->m_IsAtEnd; }; /** Walk forward one index */ void operator++() override { this->DoExtendedFloodStep(); } virtual SmartPointer GetFunction() const { return m_Function; } /** operator= is provided to make sure the handle to the image is properly * reference counted. */ Self &operator=(const Self &it) { this->m_Image = it.m_Image; // copy the smart pointer this->m_Region = it.m_Region; + this->m_InitializeValue = it.m_InitializeValue; + this->m_RegionGrowingState = it.m_RegionGrowingState; + this->m_MinTH = it.m_MinTH; + this->m_MaxTH = it.m_MaxTH; + this->m_SeedPointValue = it.m_SeedPointValue; + this->m_VoxelCounter = it.m_VoxelCounter; + this->m_LastVoxelNumber = it.m_LastVoxelNumber; + this->m_DetectedLeakagePoint = it.m_DetectedLeakagePoint; + this->m_CurrentLeakageRatio = it.m_CurrentLeakageRatio; return *this; } /** Compute whether the index of interest should be included in the flood */ bool IsPixelIncluded(const IndexType &index) const override; // Calculate the value the outputImage is initialized to static int CalculateInitializeValue(int lower, int upper) { return ((upper - lower) + 1) * (-1); }; int GetLeakagePoint(void) { return m_DetectedLeakagePoint; } protected: /* * @brief Pointer on the output image to which the result shall be written */ SmartPointer m_OutputImage; SmartPointer m_Function; /** A list of locations to start the recursive fill */ std::vector m_StartIndices; /** The origin of the source image */ typename ImageType::PointType m_ImageOrigin; /** The spacing of the source image */ typename ImageType::SpacingType m_ImageSpacing; /** Region of the source image */ RegionType m_ImageRegion; bool m_UpwardsExpansion; int m_InitializeValue; void ExpandThresholdUpwards(); void ExpandThresholdDownwards(); void IncrementRegionGrowingState(); // calculates how many steps the voxel is from the current step int EstimateDistance(IndexType); // calculates how many expansion steps will be taken unsigned int CalculateMaxRGS(); private: void InsertIndexTypeIntoQueueMap(unsigned int key, IndexType index); int m_RegionGrowingState; QueueMapType m_QueueMap; int m_MinTH; int m_MaxTH; int m_SeedPointValue; unsigned int m_VoxelCounter; unsigned int m_LastVoxelNumber; int m_DetectedLeakagePoint; float m_CurrentLeakageRatio; void CheckSeedPointValue(); /* flag for switching between raw leakage detection (bigger bronchial vessels) * and fine leakage detection (smaller bronchial vessels [starting from leaves]) */ bool m_FineDetectionMode; bool m_DetectionStop; }; } // end namespace itk #ifndef ITK_MANUAL_INSTANTIATION #include "itkAdaptiveThresholdIterator.txx" #endif #endif diff --git a/Modules/Segmentation/Algorithms/itkConnectedAdaptiveThresholdImageFilter.txx b/Modules/Segmentation/Algorithms/itkConnectedAdaptiveThresholdImageFilter.txx index 6d740e27f0..cb2a46ef18 100644 --- a/Modules/Segmentation/Algorithms/itkConnectedAdaptiveThresholdImageFilter.txx +++ b/Modules/Segmentation/Algorithms/itkConnectedAdaptiveThresholdImageFilter.txx @@ -1,300 +1,309 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _itkConnectedAdaptiveThresholdImageFilter_txx #define _itkConnectedAdaptiveThresholdImageFilter_txx #include "itkAdaptiveThresholdIterator.h" #include "itkBinaryThresholdImageFunction.h" #include "itkConnectedAdaptiveThresholdImageFilter.h" #include "itkMinimumMaximumImageFilter.h" #include "itkThresholdImageFilter.h" namespace itk { /** * Constructor */ template ConnectedAdaptiveThresholdImageFilter::ConnectedAdaptiveThresholdImageFilter() - : m_FineDetectionMode(false) + : m_OutoutImageMaskFineSegmentation(nullptr), + m_GrowingDirectionIsUpwards(true), + m_SeedpointValue(0), + m_DetectedLeakagePoint(0), + m_InitValue(0), + m_AdjLowerTh(0), + m_AdjUpperTh(0), + m_FineDetectionMode(false), + m_DiscardLastPreview(false), + m_SegmentationCancelled(false) { } template void ConnectedAdaptiveThresholdImageFilter::GenerateData() { typename ConnectedAdaptiveThresholdImageFilter::InputImageConstPointer inputImage = this->GetInput(); typename ConnectedAdaptiveThresholdImageFilter::OutputImagePointer outputImage = this->GetOutput(); typename Superclass::InputPixelObjectType::Pointer lowerThreshold = this->GetLowerInput(); typename Superclass::InputPixelObjectType::Pointer upperThreshold = this->GetUpperInput(); // kommt drauf, wie wir hier die Pipeline aufbauen this->SetLower(lowerThreshold->Get()); this->SetUpper(upperThreshold->Get()); typedef BinaryThresholdImageFunction FunctionType; typedef AdaptiveThresholdIterator IteratorType; int initValue = IteratorType::CalculateInitializeValue((int)(this->GetLower()), (int)(this->GetUpper())); // Initialize the output according to the segmentation (fine or raw) if (m_FineDetectionMode) { outputImage = this->m_OutoutImageMaskFineSegmentation; } typename ConnectedAdaptiveThresholdImageFilter::OutputImageRegionType region = outputImage->GetRequestedRegion(); outputImage->SetBufferedRegion(region); outputImage->Allocate(); if (!m_FineDetectionMode) { // only initalize the output image if we are using the raw segmentation mode outputImage->FillBuffer((typename ConnectedAdaptiveThresholdImageFilter::OutputImagePixelType)initValue); } typename FunctionType::Pointer function = FunctionType::New(); function->SetInputImage(inputImage); typename Superclass::SeedContainerType seeds; seeds = this->GetSeeds(); // pass parameters needed for region growing to iterator IteratorType it(outputImage, function, seeds); it.SetFineDetectionMode(m_FineDetectionMode); it.SetExpansionDirection(m_GrowingDirectionIsUpwards); it.SetMinTH((int)(this->GetLower())); it.SetMaxTH((int)(this->GetUpper())); it.GoToBegin(); this->m_SeedpointValue = it.GetSeedPointValue(); if ((this->GetLower()) > this->m_SeedpointValue || this->m_SeedpointValue > (this->GetUpper())) { // set m_SegmentationCancelled to true, so if it doesn't reach the point where it is set back to false // we can asssume that there was an error this->m_SegmentationCancelled = true; return; } // iterate through image until while (!it.IsAtEnd()) { // make iterator go one step further (calls method DoFloodStep()) ++it; } this->m_DetectedLeakagePoint = it.GetLeakagePoint(); this->m_SegmentationCancelled = false; } template TOutputImage *itk::ConnectedAdaptiveThresholdImageFilter::GetResultImage() { return m_OutoutImageMaskFineSegmentation; } template typename ConnectedAdaptiveThresholdImageFilter::IndexType itk::ConnectedAdaptiveThresholdImageFilter::CorrectSeedPointPosition( unsigned int sizeOfVolume, int lowerTh, int upperTh) { typename ConnectedAdaptiveThresholdImageFilter::InputImageConstPointer inputImage = this->GetInput(); typedef typename TInputImage::IndexType IndexType; IndexType itkIntelligentSeedIndex; int seedPixelValue = inputImage->GetPixel(m_SeedPointIndex); // set new seed index to the voxel with the darkest value and shortest distance to original seed if (seedPixelValue > upperTh || seedPixelValue < lowerTh) { // MITK_INFO << "seed pixel value [BEFORE] = " << seedPixelValue; // ToDo crop region itk::Index<3> workindex; for (int i = 0; i < 3; i++) { workindex[i] = m_SeedPointIndex[i] - sizeOfVolume / 2; if (workindex[i] < 0) workindex[i] = 0; } itk::Size<3> worksize; for (int i = 0; i < 3; i++) { worksize[i] = sizeOfVolume; } itk::ImageRegion<3> workregion(workindex, worksize); itk::ImageRegionIterator regionIt(const_cast(inputImage.GetPointer()), workregion); // int darkestGrayValue=seedPixelValue; int currentGrayValue; float distance = (float)(sizeOfVolume / 2); float relativeDistance = 1; // between 0 and 1 mitk::Vector3D seedVector, currentVector; mitk::FillVector3D(seedVector, m_SeedPointIndex[0], m_SeedPointIndex[1], m_SeedPointIndex[2]); currentVector = seedVector; float costValue = 0; // beware, Depending on seeking upper or lower value... for (regionIt.GoToBegin(); !regionIt.IsAtEnd(); ++regionIt) { // get current gray value currentGrayValue = regionIt.Value(); // get current seed index m_SeedPointIndex = regionIt.GetIndex(); // fill current vector mitk::FillVector3D(currentVector, m_SeedPointIndex[0], m_SeedPointIndex[1], m_SeedPointIndex[2]); // calculate distance from original seed to new seed mitk::Vector3D distVector = currentVector - seedVector; distance = fabs(distVector.GetSquaredNorm()); relativeDistance = distance / (sizeOfVolume / 2); // calculate "cost function" float currentCostValue = (1 - relativeDistance) * currentGrayValue; if (currentCostValue < costValue && currentGrayValue < upperTh) { itkIntelligentSeedIndex = regionIt.GetIndex(); costValue = currentCostValue; // MITK_INFO <<"cost value="<< costValue; // MITK_INFO <<"darkest and closest Voxel ="<< currentGrayValue; // MITK_INFO <<"m_UPPER="<< upperTh; } } // MITK_INFO<< "seed pixel value [AFTER] =" << inputImage->GetPixel(itkIntelligentSeedIndex) <<"\n"; } else { // no correction of the seed point is needed, just pass the original seed itkIntelligentSeedIndex = m_SeedPointIndex; } return itkIntelligentSeedIndex; } template void itk::ConnectedAdaptiveThresholdImageFilter::CropMask(unsigned int croppingSize) { // initialize center point of the working region itk::Index<3> workindex; for (int i = 0; i < 3; i++) { workindex[i] = m_SeedPointIndex[i] - croppingSize / 2; if (workindex[i] < 0) workindex[i] = 0; } // initialize working volume itk::Size<3> worksize; for (int i = 0; i < 3; i++) { worksize[i] = croppingSize; } // set working region itk::ImageRegion<3> workregion(workindex, worksize); // check if the entire region is inside the image if (!(m_OutoutImageMaskFineSegmentation->GetLargestPossibleRegion().IsInside(workregion))) { // if not then crop to the intersection of the image (gemeinsame Schnittmenge Bild und workingRegion) if (!(workregion.Crop(m_OutoutImageMaskFineSegmentation->GetLargestPossibleRegion()))) { MITK_ERROR << "Cropping working region failed!"; return; } } // initialize region iterator itk::ImageRegionIterator regionIt(m_OutoutImageMaskFineSegmentation, workregion); for (regionIt.GoToBegin(); !regionIt.IsAtEnd(); ++regionIt) { // and set all voxel inside the working region to zero regionIt.Set(0); } } template unsigned int itk::ConnectedAdaptiveThresholdImageFilter::AdjustIteratorMask() { typedef itk::ThresholdImageFilter ThresholdFilterType; typedef itk::MinimumMaximumImageFilter MaxFilterType; typename ThresholdFilterType::Pointer thresholdFilter = ThresholdFilterType::New(); typename MaxFilterType::Pointer maxFilter = MaxFilterType::New(); unsigned int maxValue; if (!m_DiscardLastPreview) { // get the biggest value of the image maxFilter->SetInput(m_OutoutImageMaskFineSegmentation); maxFilter->UpdateLargestPossibleRegion(); maxValue = maxFilter->GetMaximum(); } else { // use the last biggest value in the preview. This was set in SetParameterForFineSegmentation(...adjLowerTh...) [] maxValue = m_AdjLowerTh; } // set all values upper to zero (thresouldOutside uses < and > NOT <= and >=) thresholdFilter->SetInput(m_OutoutImageMaskFineSegmentation); thresholdFilter->SetOutsideValue(0); thresholdFilter->ThresholdOutside(m_AdjLowerTh, maxValue); thresholdFilter->UpdateLargestPossibleRegion(); // set all values in between lower and upper (>=lower && <=upper) to the highest value in the image thresholdFilter->SetInput(thresholdFilter->GetOutput()); thresholdFilter->SetOutsideValue(maxValue); thresholdFilter->ThresholdOutside(0, m_AdjLowerTh - 1); thresholdFilter->UpdateLargestPossibleRegion(); m_OutoutImageMaskFineSegmentation = thresholdFilter->GetOutput(); return maxValue; } template void itk::ConnectedAdaptiveThresholdImageFilter::SetParameterForFineSegmentation( TOutputImage *iteratorMaskForFineSegmentation, unsigned int adjLowerTh, unsigned int adjUpperTh, itk::Index<3> seedPoint, bool discardLeafSegmentation) { // just to make sure we´re in the right mode and the mask exsits if (m_FineDetectionMode && iteratorMaskForFineSegmentation) { m_OutoutImageMaskFineSegmentation = iteratorMaskForFineSegmentation; m_AdjLowerTh = adjLowerTh; m_AdjUpperTh = adjUpperTh; // still needed? m_SeedPointIndex = seedPoint; m_DiscardLastPreview = discardLeafSegmentation; } else { if (!m_FineDetectionMode) { MITK_ERROR << "Fine-detection-segmentation mode not set!"; } else { MITK_ERROR << "Iterator-mask-image not set!"; } } } } // end namespace itk #endif diff --git a/Modules/Segmentation/Algorithms/mitkCalculateSegmentationVolume.cpp b/Modules/Segmentation/Algorithms/mitkCalculateSegmentationVolume.cpp index 8f8305f90e..7fe3898104 100644 --- a/Modules/Segmentation/Algorithms/mitkCalculateSegmentationVolume.cpp +++ b/Modules/Segmentation/Algorithms/mitkCalculateSegmentationVolume.cpp @@ -1,122 +1,122 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include "itkImageRegionConstIteratorWithIndex.h" #include "mitkCalculateSegmentationVolume.h" #include namespace mitk { - CalculateSegmentationVolume::CalculateSegmentationVolume() {} + CalculateSegmentationVolume::CalculateSegmentationVolume() : m_Volume(0){} CalculateSegmentationVolume::~CalculateSegmentationVolume() {} template void CalculateSegmentationVolume::ItkImageProcessing(itk::Image *itkImage, TPixel *itkNotUsed(dummy)) { itk::ImageRegionConstIteratorWithIndex> iterBinaryImage( itkImage, itkImage->GetLargestPossibleRegion()); typename itk::ImageRegionConstIteratorWithIndex>::IndexType currentIndex; typename itk::ImageRegionConstIteratorWithIndex>::IndexType minIndex; for (unsigned int i = 0; i < VImageDimension; ++i) minIndex[i] = std::numeric_limits::max(); typename itk::ImageRegionConstIteratorWithIndex>::IndexType maxIndex; for (unsigned int i = 0; i < VImageDimension; ++i) maxIndex[i] = std::numeric_limits::min(); m_CenterOfMass.Fill(0.0); m_Volume = 0; while (!iterBinaryImage.IsAtEnd()) { if (iterBinaryImage.Get() > static_cast(0.0)) { // update center of mass currentIndex = iterBinaryImage.GetIndex(); itk::Vector currentPoint; for (unsigned int i = 0; i < VImageDimension; ++i) currentPoint[i] = currentIndex[i]; m_CenterOfMass = (m_CenterOfMass * (static_cast(m_Volume) / static_cast(m_Volume + 1))) // e.g. 3 points: old center * 2/3 + currentPoint * 1/3; + currentPoint / static_cast(m_Volume + 1); // update number of voxels ++m_Volume; // update bounding box for (unsigned int i = 0; i < VImageDimension; ++i) { if (currentIndex[i] < minIndex[i]) minIndex[i] = currentIndex[i]; if (currentIndex[i] > maxIndex[i]) maxIndex[i] = currentIndex[i]; } } ++iterBinaryImage; } m_MinIndexOfBoundingBox[2] = 0.0; m_MaxIndexOfBoundingBox[2] = 0.0; for (unsigned int i = 0; i < VImageDimension; ++i) { m_MinIndexOfBoundingBox[i] = minIndex[i]; m_MaxIndexOfBoundingBox[i] = maxIndex[i]; } } bool CalculateSegmentationVolume::ReadyToRun() { Image::Pointer image; GetPointerParameter("Input", image); return image.IsNotNull() && GetGroupNode(); } bool CalculateSegmentationVolume::ThreadedUpdateFunction() { // get image Image::Pointer image; GetPointerParameter("Input", image); AccessFixedDimensionByItk(image.GetPointer(), ItkImageProcessing, 3); // some magic to call the correctly templated function (we only do 3D images here!) // consider single voxel volume Vector3D spacing = image->GetSlicedGeometry()->GetSpacing(); // spacing in mm float volumeML = (ScalarType)m_Volume * spacing[0] * spacing[1] * spacing[2] / 1000.0; // convert to ml DataNode *groupNode = GetGroupNode(); if (groupNode) { groupNode->SetProperty("volume", FloatProperty::New(volumeML)); groupNode->SetProperty("centerOfMass", Vector3DProperty::New(m_CenterOfMass)); groupNode->SetProperty("boundingBoxMinimum", Vector3DProperty::New(m_MinIndexOfBoundingBox)); groupNode->SetProperty("boundingBoxMaximum", Vector3DProperty::New(m_MaxIndexOfBoundingBox)); groupNode->SetProperty("showVolume", BoolProperty::New(true)); } return true; } } // namespace diff --git a/Modules/Segmentation/Algorithms/mitkDiffImageApplier.cpp b/Modules/Segmentation/Algorithms/mitkDiffImageApplier.cpp index 5ee9584ac0..64a3787caf 100644 --- a/Modules/Segmentation/Algorithms/mitkDiffImageApplier.cpp +++ b/Modules/Segmentation/Algorithms/mitkDiffImageApplier.cpp @@ -1,355 +1,363 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkDiffImageApplier.h" #include "mitkApplyDiffImageOperation.h" #include "mitkImageAccessByItk.h" #include "mitkImageCast.h" #include "mitkImageTimeSelector.h" #include "mitkRenderingManager.h" #include "mitkSegmentationInterpolationController.h" #include #include #include mitk::DiffImageApplier::DiffImageApplier() + : m_Image(nullptr), + m_SliceDifferenceImage(nullptr), + m_SliceIndex(0), + m_SliceDimension(0), + m_TimeStep(0), + m_Dimension0(0), + m_Dimension1(0), + m_Factor(1.0) { } mitk::DiffImageApplier::~DiffImageApplier() { } void mitk::DiffImageApplier::ExecuteOperation(Operation *operation) { auto *imageOperation = dynamic_cast(operation); if (imageOperation // we actually have the kind of operation that we can handle && imageOperation->IsImageStillValid()) // AND the image is not yet deleted { m_Image = imageOperation->GetImage(); Image::Pointer image3D = m_Image; // will be changed later in case of 3D+t m_SliceDifferenceImage = imageOperation->GetDiffImage(); m_TimeStep = imageOperation->GetTimeStep(); m_Factor = imageOperation->GetFactor(); if (m_SliceDifferenceImage->GetDimension() == 2) { m_SliceIndex = imageOperation->GetSliceIndex(); m_SliceDimension = imageOperation->GetSliceDimension(); switch (m_SliceDimension) { default: case 2: m_Dimension0 = 0; m_Dimension1 = 1; break; case 1: m_Dimension0 = 0; m_Dimension1 = 2; break; case 0: m_Dimension0 = 1; m_Dimension1 = 2; break; } if (m_SliceDifferenceImage->GetDimension() != 2 || (m_Image->GetDimension() < 3 || m_Image->GetDimension() > 4) || m_SliceDifferenceImage->GetDimension(0) != m_Image->GetDimension(m_Dimension0) || m_SliceDifferenceImage->GetDimension(1) != m_Image->GetDimension(m_Dimension1) || m_SliceIndex >= m_Image->GetDimension(m_SliceDimension)) { itkExceptionMacro( "Slice and image dimensions differ or slice index is too large. Sorry, cannot work like this."); return; } if (m_Image->GetDimension() == 4) { ImageTimeSelector::Pointer timeSelector = ImageTimeSelector::New(); timeSelector->SetInput(m_Image); timeSelector->SetTimeNr(m_TimeStep); timeSelector->UpdateLargestPossibleRegion(); image3D = timeSelector->GetOutput(); } // this will do a long long if/else to find out both pixel types AccessFixedDimensionByItk(image3D, ItkImageSwitch2DDiff, 3); if (m_Factor == 1 || m_Factor == -1) { if (m_Factor == -1) { // multiply diff pixels by factor and then send this diff slice AccessFixedDimensionByItk(m_SliceDifferenceImage, ItkInvertPixelValues, 2); } // just send the diff to SegmentationInterpolationController SegmentationInterpolationController *interpolator = SegmentationInterpolationController::InterpolatorForImage(m_Image); if (interpolator) { interpolator->BlockModified(true); interpolator->SetChangedSlice(m_SliceDifferenceImage, m_SliceDimension, m_SliceIndex, m_TimeStep); } m_Image->Modified(); if (interpolator) { interpolator->BlockModified(false); } if (m_Factor == -1) // return to normal values { AccessFixedDimensionByItk(m_SliceDifferenceImage, ItkInvertPixelValues, 2); } } else // no trivial case, too lazy to do something else { m_Image->Modified(); // check if interpolation is called. prefer to send diff directly } RenderingManager::GetInstance()->RequestUpdateAll(); } else if (m_SliceDifferenceImage->GetDimension() == 3) { // ... if (m_SliceDifferenceImage->GetDimension(0) != m_Image->GetDimension(0) || m_SliceDifferenceImage->GetDimension(1) != m_Image->GetDimension(1) || m_SliceDifferenceImage->GetDimension(2) != m_Image->GetDimension(2) || m_TimeStep >= m_Image->GetDimension(3)) { itkExceptionMacro("Diff image size differs from original image size. Sorry, cannot work like this."); return; } if (m_Image->GetDimension() == 4) { ImageTimeSelector::Pointer timeSelector = ImageTimeSelector::New(); timeSelector->SetInput(m_Image); timeSelector->SetTimeNr(m_TimeStep); timeSelector->UpdateLargestPossibleRegion(); image3D = timeSelector->GetOutput(); } // this will do a long long if/else to find out both pixel types AccessFixedDimensionByItk(image3D, ItkImageSwitch3DDiff, 3); if (m_Factor == 1 || m_Factor == -1) { if (m_Factor == -1) { // multiply diff pixels by factor and then send this diff slice AccessFixedDimensionByItk(m_SliceDifferenceImage, ItkInvertPixelValues, 3); } // just send the diff to SegmentationInterpolationController SegmentationInterpolationController *interpolator = SegmentationInterpolationController::InterpolatorForImage(m_Image); if (interpolator) { interpolator->BlockModified(true); interpolator->SetChangedVolume(m_SliceDifferenceImage, m_TimeStep); } m_Image->Modified(); if (interpolator) { interpolator->BlockModified(false); } if (m_Factor == -1) // return to normal values { AccessFixedDimensionByItk(m_SliceDifferenceImage, ItkInvertPixelValues, 3); } } else // no trivial case, too lazy to do something else { m_Image->Modified(); // check if interpolation is called. prefer to send diff directly } RenderingManager::GetInstance()->RequestUpdateAll(); } else { itkExceptionMacro("Diff image must be 2D or 3D. Sorry, cannot work like this."); return; } } m_Image = nullptr; m_SliceDifferenceImage = nullptr; } mitk::DiffImageApplier *mitk::DiffImageApplier::GetInstanceForUndo() { static DiffImageApplier::Pointer s_Instance = DiffImageApplier::New(); return s_Instance; } // basically copied from mitk/Core/Algorithms/mitkImageAccessByItk.h #define myMITKDiffImageApplierFilterAccessByItk(mitkImage, itkImageTypeFunction, pixeltype, dimension, itkimage2) \ if (typeId == MapPixelComponentType::value) \ \ { \ typedef itk::Image ImageType; \ typedef mitk::ImageToItk ImageToItkType; \ itk::SmartPointer imagetoitk = ImageToItkType::New(); \ const mitk::Image *constImage = mitkImage; \ mitk::Image *nonConstImage = const_cast(constImage); \ nonConstImage->Update(); \ imagetoitk->SetInput(nonConstImage); \ imagetoitk->Update(); \ itkImageTypeFunction(imagetoitk->GetOutput(), itkimage2); \ \ } #define myMITKDiffImageApplierFilterAccessAllTypesByItk(mitkImage, itkImageTypeFunction, dimension, itkimage2) \ \ { \ myMITKDiffImageApplierFilterAccessByItk(mitkImage, itkImageTypeFunction, double, dimension, itkimage2) else myMITKDiffImageApplierFilterAccessByItk( \ mitkImage, \ itkImageTypeFunction, \ float, \ dimension, \ itkimage2) else myMITKDiffImageApplierFilterAccessByItk(mitkImage, itkImageTypeFunction, int, dimension, itkimage2) else myMITKDiffImageApplierFilterAccessByItk(mitkImage, \ itkImageTypeFunction, \ unsigned int, \ dimension, \ itkimage2) else myMITKDiffImageApplierFilterAccessByItk(mitkImage, itkImageTypeFunction, short, dimension, itkimage2) else myMITKDiffImageApplierFilterAccessByItk(mitkImage, itkImageTypeFunction, unsigned short, dimension, itkimage2) else myMITKDiffImageApplierFilterAccessByItk(mitkImage, \ itkImageTypeFunction, \ char, \ dimension, \ itkimage2) else myMITKDiffImageApplierFilterAccessByItk(mitkImage, \ itkImageTypeFunction, \ unsigned char, \ dimension, \ itkimage2) \ \ } template void mitk::DiffImageApplier::ItkImageSwitch2DDiff(itk::Image *itkImage) { const int typeId = m_SliceDifferenceImage->GetPixelType().GetComponentType(); myMITKDiffImageApplierFilterAccessAllTypesByItk(m_SliceDifferenceImage, ItkImageProcessing2DDiff, 2, itkImage); } template void mitk::DiffImageApplier::ItkImageSwitch3DDiff(itk::Image *itkImage) { const int typeId = m_SliceDifferenceImage->GetPixelType().GetComponentType(); myMITKDiffImageApplierFilterAccessAllTypesByItk(m_SliceDifferenceImage, ItkImageProcessing3DDiff, 3, itkImage); } template void mitk::DiffImageApplier::ItkImageProcessing2DDiff(itk::Image *diffImage, itk::Image *outputImage) { typedef itk::Image DiffImageType; typedef itk::Image VolumeImageType; typedef itk::ImageSliceIteratorWithIndex OutputSliceIteratorType; typedef itk::ImageRegionConstIterator DiffSliceIteratorType; typename VolumeImageType::RegionType sliceInVolumeRegion; sliceInVolumeRegion = outputImage->GetLargestPossibleRegion(); sliceInVolumeRegion.SetSize(m_SliceDimension, 1); // just one slice sliceInVolumeRegion.SetIndex(m_SliceDimension, m_SliceIndex); // exactly this slice, please OutputSliceIteratorType outputIterator(outputImage, sliceInVolumeRegion); outputIterator.SetFirstDirection(m_Dimension0); outputIterator.SetSecondDirection(m_Dimension1); DiffSliceIteratorType diffIterator(diffImage, diffImage->GetLargestPossibleRegion()); // iterate over output slice (and over input slice simultaneously) outputIterator.GoToBegin(); diffIterator.GoToBegin(); while (!outputIterator.IsAtEnd()) { while (!outputIterator.IsAtEndOfSlice()) { while (!outputIterator.IsAtEndOfLine()) { TPixel2 newValue = outputIterator.Get() + (TPixel2)((double)diffIterator.Get() * m_Factor); outputIterator.Set(newValue); ++outputIterator; ++diffIterator; } outputIterator.NextLine(); } outputIterator.NextSlice(); } } template void mitk::DiffImageApplier::ItkImageProcessing3DDiff(itk::Image *diffImage, itk::Image *outputImage) { typedef itk::Image DiffImageType; typedef itk::Image VolumeImageType; typedef itk::ImageRegionIterator OutputSliceIteratorType; typedef itk::ImageRegionConstIterator DiffSliceIteratorType; OutputSliceIteratorType outputIterator(outputImage, outputImage->GetLargestPossibleRegion()); DiffSliceIteratorType diffIterator(diffImage, diffImage->GetLargestPossibleRegion()); // iterate over output slice (and over input slice simultaneously) outputIterator.GoToBegin(); diffIterator.GoToBegin(); while (!outputIterator.IsAtEnd()) { TPixel2 newValue = outputIterator.Get() + (TPixel2)((double)diffIterator.Get() * m_Factor); outputIterator.Set(newValue); ++outputIterator; ++diffIterator; } } #ifdef _MSC_VER # pragma warning(push) # pragma warning(disable:4146) // unary minus operator applied to unsigned type, result still unsigned #endif template void mitk::DiffImageApplier::ItkInvertPixelValues(itk::Image *itkImage) { typedef itk::ImageRegionIterator> IteratorType; IteratorType iter(itkImage, itkImage->GetLargestPossibleRegion()); iter.GoToBegin(); while (!iter.IsAtEnd()) { iter.Set(-(iter.Get())); ++iter; } } #ifdef _MSC_VER # pragma warning(pop) #endif diff --git a/Modules/Segmentation/Algorithms/mitkDiffSliceOperation.cpp b/Modules/Segmentation/Algorithms/mitkDiffSliceOperation.cpp index 7099e138e6..2982be76c0 100644 --- a/Modules/Segmentation/Algorithms/mitkDiffSliceOperation.cpp +++ b/Modules/Segmentation/Algorithms/mitkDiffSliceOperation.cpp @@ -1,102 +1,104 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkDiffSliceOperation.h" #include #include mitk::DiffSliceOperation::DiffSliceOperation() : Operation(1) { m_TimeStep = 0; m_zlibSliceContainer = nullptr; m_Image = nullptr; m_WorldGeometry = nullptr; m_SliceGeometry = nullptr; m_ImageIsValid = false; + m_DeleteObserverTag = 0; } mitk::DiffSliceOperation::DiffSliceOperation(mitk::Image *imageVolume, Image *slice, SlicedGeometry3D *sliceGeometry, unsigned int timestep, BaseGeometry *currentWorldGeometry) : Operation(1) { m_WorldGeometry = currentWorldGeometry->Clone(); /* Quick fix for bug 12338. Guard object - fix this when clone method of PlaneGeometry is cloning the reference geometry (see bug 13392)*/ // xxxx m_GuardReferenceGeometry = mitk::BaseGeometry::New(); m_GuardReferenceGeometry = dynamic_cast(m_WorldGeometry.GetPointer())->GetReferenceGeometry(); /*---------------------------------------------------------------------------------------------------*/ m_SliceGeometry = sliceGeometry->Clone(); m_TimeStep = timestep; m_zlibSliceContainer = CompressedImageContainer::New(); m_zlibSliceContainer->SetImage(slice); m_Image = imageVolume; + m_DeleteObserverTag = 0; if (m_Image) { /*add an observer to listen to the delete event of the image, this is necessary because the operation is then * invalid*/ itk::SimpleMemberCommand::Pointer command = itk::SimpleMemberCommand::New(); command->SetCallbackFunction(this, &DiffSliceOperation::OnImageDeleted); // get the id of the observer, used to remove it later on m_DeleteObserverTag = imageVolume->AddObserver(itk::DeleteEvent(), command); m_ImageIsValid = true; } else m_ImageIsValid = false; } mitk::DiffSliceOperation::~DiffSliceOperation() { m_WorldGeometry = nullptr; m_zlibSliceContainer = nullptr; if (m_ImageIsValid) { // if the image is still there, we have to remove the observer from it m_Image->RemoveObserver(m_DeleteObserverTag); } m_Image = nullptr; } mitk::Image::Pointer mitk::DiffSliceOperation::GetSlice() { Image::Pointer image = m_zlibSliceContainer->GetImage(); return image; } bool mitk::DiffSliceOperation::IsValid() { return m_ImageIsValid && m_zlibSliceContainer.IsNotNull() && (m_WorldGeometry.IsNotNull()); // TODO improve } void mitk::DiffSliceOperation::OnImageDeleted() { // if our imageVolume is removed e.g. from the datastorage the operation is no lnger valid m_ImageIsValid = false; } diff --git a/Modules/Segmentation/Algorithms/mitkImageToContourFilter.cpp b/Modules/Segmentation/Algorithms/mitkImageToContourFilter.cpp index ef6fe0753f..f48870f246 100644 --- a/Modules/Segmentation/Algorithms/mitkImageToContourFilter.cpp +++ b/Modules/Segmentation/Algorithms/mitkImageToContourFilter.cpp @@ -1,157 +1,158 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkImageToContourFilter.h" #include "mitkImageAccessByItk.h" #include "mitkImageCast.h" #include "mitkVtkRepresentationProperty.h" #include "vtkLinearTransform.h" #include "vtkMatrix4x4.h" #include "vtkProperty.h" #include "vtkSmartPointer.h" #include mitk::ImageToContourFilter::ImageToContourFilter() { this->m_UseProgressBar = false; this->m_ProgressStepSize = 1; + this->m_SliceGeometry = nullptr; } mitk::ImageToContourFilter::~ImageToContourFilter() { } void mitk::ImageToContourFilter::GenerateData() { mitk::Image::ConstPointer sliceImage = ImageToSurfaceFilter::GetInput(); if (!sliceImage) { MITK_ERROR << "mitk::ImageToContourFilter: No input available. Please set the input!" << std::endl; itkExceptionMacro("mitk::ImageToContourFilter: No input available. Please set the input!"); return; } if (sliceImage->GetDimension() > 2 || sliceImage->GetDimension() < 2) { MITK_ERROR << "mitk::ImageToImageFilter::GenerateData() works only with 2D images. Please assure that your input " "image is 2D!" << std::endl; itkExceptionMacro( "mitk::ImageToImageFilter::GenerateData() works only with 2D images. Please assure that your input image is 2D!"); return; } m_SliceGeometry = sliceImage->GetGeometry(); AccessFixedDimensionByItk(sliceImage, Itk2DContourExtraction, 2); // Setting progressbar if (this->m_UseProgressBar) mitk::ProgressBar::GetInstance()->Progress(this->m_ProgressStepSize); } template void mitk::ImageToContourFilter::Itk2DContourExtraction(const itk::Image *sliceImage) { typedef itk::Image ImageType; typedef itk::ContourExtractor2DImageFilter ContourExtractor; typedef itk::ConstantPadImageFilter PadFilterType; typename PadFilterType::Pointer padFilter = PadFilterType::New(); typename ImageType::SizeType lowerExtendRegion; lowerExtendRegion[0] = 1; lowerExtendRegion[1] = 1; typename ImageType::SizeType upperExtendRegion; upperExtendRegion[0] = 1; upperExtendRegion[1] = 1; /* * We need to pad here, since the ITK contour extractor fails if the * segmentation touches more than one image edge. * By padding the image for one row at each edge we overcome this issue */ padFilter->SetInput(sliceImage); padFilter->SetConstant(0); padFilter->SetPadLowerBound(lowerExtendRegion); padFilter->SetPadUpperBound(upperExtendRegion); typename ContourExtractor::Pointer contourExtractor = ContourExtractor::New(); contourExtractor->SetInput(padFilter->GetOutput()); contourExtractor->SetContourValue(0.5); contourExtractor->Update(); unsigned int foundPaths = contourExtractor->GetNumberOfOutputs(); vtkSmartPointer contourSurface = vtkSmartPointer::New(); vtkSmartPointer points = vtkSmartPointer::New(); vtkSmartPointer polygons = vtkSmartPointer::New(); unsigned int pointId(0); for (unsigned int i = 0; i < foundPaths; i++) { const ContourPath *currentPath = contourExtractor->GetOutput(i)->GetVertexList(); vtkSmartPointer polygon = vtkSmartPointer::New(); polygon->GetPointIds()->SetNumberOfIds(currentPath->Size()); Point3D currentPoint; Point3D currentWorldPoint; for (unsigned int j = 0; j < currentPath->Size(); j++) { currentPoint[0] = currentPath->ElementAt(j)[0]; currentPoint[1] = currentPath->ElementAt(j)[1]; currentPoint[2] = 0; m_SliceGeometry->IndexToWorld(currentPoint, currentWorldPoint); points->InsertPoint(pointId, currentWorldPoint[0], currentWorldPoint[1], currentWorldPoint[2]); polygon->GetPointIds()->SetId(j, pointId); pointId++; } // for2 polygons->InsertNextCell(polygon); } // for1 contourSurface->SetPoints(points); contourSurface->SetPolys(polygons); contourSurface->BuildLinks(); Surface::Pointer finalSurface = this->GetOutput(); finalSurface->SetVtkPolyData(contourSurface); } void mitk::ImageToContourFilter::GenerateOutputInformation() { Superclass::GenerateOutputInformation(); } void mitk::ImageToContourFilter::SetUseProgressBar(bool status) { this->m_UseProgressBar = status; } void mitk::ImageToContourFilter::SetProgressStepSize(unsigned int stepSize) { this->m_ProgressStepSize = stepSize; } diff --git a/Modules/Segmentation/Controllers/mitkSegmentationInterpolationController.cpp b/Modules/Segmentation/Controllers/mitkSegmentationInterpolationController.cpp index 1031296c03..87040ca083 100644 --- a/Modules/Segmentation/Controllers/mitkSegmentationInterpolationController.cpp +++ b/Modules/Segmentation/Controllers/mitkSegmentationInterpolationController.cpp @@ -1,589 +1,590 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkSegmentationInterpolationController.h" #include "mitkImageCast.h" #include "mitkImageReadAccessor.h" #include "mitkImageTimeSelector.h" #include #include //#include #include "mitkShapeBasedInterpolationAlgorithm.h" #include #include #include mitk::SegmentationInterpolationController::InterpolatorMapType mitk::SegmentationInterpolationController::s_InterpolatorForImage; // static member initialization mitk::SegmentationInterpolationController *mitk::SegmentationInterpolationController::InterpolatorForImage( const Image *image) { auto iter = s_InterpolatorForImage.find(image); if (iter != s_InterpolatorForImage.end()) { return iter->second; } else { return nullptr; } } -mitk::SegmentationInterpolationController::SegmentationInterpolationController() : m_BlockModified(false) +mitk::SegmentationInterpolationController::SegmentationInterpolationController() + : m_BlockModified(false), m_2DInterpolationActivated(false) { } void mitk::SegmentationInterpolationController::Activate2DInterpolation(bool status) { m_2DInterpolationActivated = status; } mitk::SegmentationInterpolationController *mitk::SegmentationInterpolationController::GetInstance() { static mitk::SegmentationInterpolationController::Pointer m_Instance; if (m_Instance.IsNull()) { m_Instance = SegmentationInterpolationController::New(); } return m_Instance; } mitk::SegmentationInterpolationController::~SegmentationInterpolationController() { // remove this from the list of interpolators for (auto iter = s_InterpolatorForImage.begin(); iter != s_InterpolatorForImage.end(); ++iter) { if (iter->second == this) { s_InterpolatorForImage.erase(iter); break; } } } void mitk::SegmentationInterpolationController::OnImageModified(const itk::EventObject &) { if (!m_BlockModified && m_Segmentation.IsNotNull() && m_2DInterpolationActivated) { SetSegmentationVolume(m_Segmentation); } } void mitk::SegmentationInterpolationController::BlockModified(bool block) { m_BlockModified = block; } void mitk::SegmentationInterpolationController::SetSegmentationVolume(const Image *segmentation) { // clear old information (remove all time steps m_SegmentationCountInSlice.clear(); // delete this from the list of interpolators auto iter = s_InterpolatorForImage.find(segmentation); if (iter != s_InterpolatorForImage.end()) { s_InterpolatorForImage.erase(iter); } if (!segmentation) return; if (segmentation->GetDimension() > 4 || segmentation->GetDimension() < 3) { itkExceptionMacro("SegmentationInterpolationController needs a 3D-segmentation or 3D+t, not 2D."); } if (m_Segmentation != segmentation) { // observe Modified() event of image itk::ReceptorMemberCommand::Pointer command = itk::ReceptorMemberCommand::New(); command->SetCallbackFunction(this, &SegmentationInterpolationController::OnImageModified); segmentation->AddObserver(itk::ModifiedEvent(), command); } m_Segmentation = segmentation; m_SegmentationCountInSlice.resize(m_Segmentation->GetTimeSteps()); for (unsigned int timeStep = 0; timeStep < m_Segmentation->GetTimeSteps(); ++timeStep) { m_SegmentationCountInSlice[timeStep].resize(3); for (unsigned int dim = 0; dim < 3; ++dim) { m_SegmentationCountInSlice[timeStep][dim].clear(); m_SegmentationCountInSlice[timeStep][dim].resize(m_Segmentation->GetDimension(dim)); m_SegmentationCountInSlice[timeStep][dim].assign(m_Segmentation->GetDimension(dim), 0); } } s_InterpolatorForImage.insert(std::make_pair(m_Segmentation, this)); // for all timesteps // scan whole image for (unsigned int timeStep = 0; timeStep < m_Segmentation->GetTimeSteps(); ++timeStep) { ImageTimeSelector::Pointer timeSelector = ImageTimeSelector::New(); timeSelector->SetInput(m_Segmentation); timeSelector->SetTimeNr(timeStep); timeSelector->UpdateLargestPossibleRegion(); Image::Pointer segmentation3D = timeSelector->GetOutput(); AccessFixedDimensionByItk_2(segmentation3D, ScanWholeVolume, 3, m_Segmentation, timeStep); } // PrintStatus(); SetReferenceVolume(m_ReferenceImage); Modified(); } void mitk::SegmentationInterpolationController::SetReferenceVolume(const Image *referenceImage) { m_ReferenceImage = referenceImage; if (m_ReferenceImage.IsNull()) return; // no image set - ignore it then assert(m_Segmentation.IsNotNull()); // should never happen // ensure the reference image has the same dimensionality and extents as the segmentation image if (m_ReferenceImage.IsNull() || m_Segmentation.IsNull() || m_ReferenceImage->GetDimension() != m_Segmentation->GetDimension() || m_ReferenceImage->GetPixelType().GetNumberOfComponents() != 1 || m_Segmentation->GetPixelType().GetNumberOfComponents() != 1) { MITK_WARN << "Segmentation image has different image characteristics than reference image." << std::endl; m_ReferenceImage = nullptr; return; } for (unsigned int dim = 0; dim < m_Segmentation->GetDimension(); ++dim) if (m_ReferenceImage->GetDimension(dim) != m_Segmentation->GetDimension(dim)) { MITK_WARN << "original patient image does not match segmentation (different extent in dimension " << dim << "), ignoring patient image" << std::endl; m_ReferenceImage = nullptr; return; } } void mitk::SegmentationInterpolationController::SetChangedVolume(const Image *sliceDiff, unsigned int timeStep) { if (!sliceDiff) return; if (sliceDiff->GetDimension() != 3) return; AccessFixedDimensionByItk_1(sliceDiff, ScanChangedVolume, 3, timeStep); // PrintStatus(); Modified(); } void mitk::SegmentationInterpolationController::SetChangedSlice(const Image *sliceDiff, unsigned int sliceDimension, unsigned int sliceIndex, unsigned int timeStep) { if (!sliceDiff) return; if (sliceDimension > 2) return; if (timeStep >= m_SegmentationCountInSlice.size()) return; if (sliceIndex >= m_SegmentationCountInSlice[timeStep][sliceDimension].size()) return; unsigned int dim0(0); unsigned int dim1(1); // determine the other two dimensions switch (sliceDimension) { default: case 2: dim0 = 0; dim1 = 1; break; case 1: dim0 = 0; dim1 = 2; break; case 0: dim0 = 1; dim1 = 2; break; } mitk::ImageReadAccessor readAccess(sliceDiff); auto *rawSlice = (unsigned char *)readAccess.GetData(); if (!rawSlice) return; AccessFixedDimensionByItk_1( sliceDiff, ScanChangedSlice, 2, SetChangedSliceOptions(sliceDimension, sliceIndex, dim0, dim1, timeStep, rawSlice)); Modified(); } template void mitk::SegmentationInterpolationController::ScanChangedSlice(const itk::Image *, const SetChangedSliceOptions &options) { auto *pixelData((DATATYPE *)options.pixelData); unsigned int timeStep(options.timeStep); unsigned int sliceDimension(options.sliceDimension); unsigned int sliceIndex(options.sliceIndex); if (sliceDimension > 2) return; if (sliceIndex >= m_SegmentationCountInSlice[timeStep][sliceDimension].size()) return; unsigned int dim0(options.dim0); unsigned int dim1(options.dim1); int numberOfPixels(0); // number of pixels in this slice that are not 0 unsigned int dim0max = m_SegmentationCountInSlice[timeStep][dim0].size(); unsigned int dim1max = m_SegmentationCountInSlice[timeStep][dim1].size(); // scan the slice from two directions // and set the flags for the two dimensions of the slice for (unsigned int v = 0; v < dim1max; ++v) { for (unsigned int u = 0; u < dim0max; ++u) { DATATYPE value = *(pixelData + u + v * dim0max); assert((signed)m_SegmentationCountInSlice[timeStep][dim0][u] + (signed)value >= 0); // just for debugging. This must always be true, otherwise some counting is going wrong assert((signed)m_SegmentationCountInSlice[timeStep][dim1][v] + (signed)value >= 0); m_SegmentationCountInSlice[timeStep][dim0][u] = static_cast(m_SegmentationCountInSlice[timeStep][dim0][u] + value); m_SegmentationCountInSlice[timeStep][dim1][v] = static_cast(m_SegmentationCountInSlice[timeStep][dim1][v] + value); numberOfPixels += static_cast(value); } } // flag for the dimension of the slice itself assert((signed)m_SegmentationCountInSlice[timeStep][sliceDimension][sliceIndex] + numberOfPixels >= 0); m_SegmentationCountInSlice[timeStep][sliceDimension][sliceIndex] += numberOfPixels; // MITK_INFO << "scan t=" << timeStep << " from (0,0) to (" << dim0max << "," << dim1max << ") (" << pixelData << "-" // << pixelData+dim0max*dim1max-1 << ") in slice " << sliceIndex << " found " << numberOfPixels << " pixels" << // std::endl; } template void mitk::SegmentationInterpolationController::ScanChangedVolume(const itk::Image *diffImage, unsigned int timeStep) { typedef itk::ImageSliceConstIteratorWithIndex> IteratorType; IteratorType iter(diffImage, diffImage->GetLargestPossibleRegion()); iter.SetFirstDirection(0); iter.SetSecondDirection(1); int numberOfPixels(0); // number of pixels in this slice that are not 0 typename IteratorType::IndexType index; unsigned int x = 0; unsigned int y = 0; unsigned int z = 0; iter.GoToBegin(); while (!iter.IsAtEnd()) { while (!iter.IsAtEndOfSlice()) { while (!iter.IsAtEndOfLine()) { index = iter.GetIndex(); x = index[0]; y = index[1]; z = index[2]; TPixel value = iter.Get(); assert((signed)m_SegmentationCountInSlice[timeStep][0][x] + (signed)value >= 0); // just for debugging. This must always be true, otherwise some counting is going wrong assert((signed)m_SegmentationCountInSlice[timeStep][1][y] + (signed)value >= 0); m_SegmentationCountInSlice[timeStep][0][x] = static_cast(m_SegmentationCountInSlice[timeStep][0][x] + value); m_SegmentationCountInSlice[timeStep][1][y] = static_cast(m_SegmentationCountInSlice[timeStep][1][y] + value); numberOfPixels += static_cast(value); ++iter; } iter.NextLine(); } assert((signed)m_SegmentationCountInSlice[timeStep][2][z] + numberOfPixels >= 0); m_SegmentationCountInSlice[timeStep][2][z] += numberOfPixels; numberOfPixels = 0; iter.NextSlice(); } } template void mitk::SegmentationInterpolationController::ScanWholeVolume(const itk::Image *, const Image *volume, unsigned int timeStep) { if (!volume) return; if (timeStep >= m_SegmentationCountInSlice.size()) return; ImageReadAccessor readAccess(volume, volume->GetVolumeData(timeStep)); for (unsigned int slice = 0; slice < volume->GetDimension(2); ++slice) { const auto *rawVolume = static_cast(readAccess.GetData()); // we again promise not to change anything, we'll just count const DATATYPE *rawSlice = rawVolume + (volume->GetDimension(0) * volume->GetDimension(1) * slice); ScanChangedSlice(nullptr, SetChangedSliceOptions(2, slice, 0, 1, timeStep, rawSlice)); } } void mitk::SegmentationInterpolationController::PrintStatus() { unsigned int timeStep(0); // if needed, put a loop over time steps around everyting, but beware, output will be long MITK_INFO << "Interpolator status (timestep 0): dimensions " << m_SegmentationCountInSlice[timeStep][0].size() << " " << m_SegmentationCountInSlice[timeStep][1].size() << " " << m_SegmentationCountInSlice[timeStep][2].size() << std::endl; MITK_INFO << "Slice 0: " << m_SegmentationCountInSlice[timeStep][2][0] << std::endl; // row "x" for (unsigned int index = 0; index < m_SegmentationCountInSlice[timeStep][0].size(); ++index) { if (m_SegmentationCountInSlice[timeStep][0][index] > 0) MITK_INFO << "O"; else MITK_INFO << "."; } MITK_INFO << std::endl; // rows "y" and "z" (diagonal) for (unsigned int index = 1; index < m_SegmentationCountInSlice[timeStep][1].size(); ++index) { if (m_SegmentationCountInSlice[timeStep][1][index] > 0) MITK_INFO << "O"; else MITK_INFO << "."; if (m_SegmentationCountInSlice[timeStep][2].size() > index) // if we also have a z value here, then print it, too { for (unsigned int indent = 1; indent < index; ++indent) MITK_INFO << " "; if (m_SegmentationCountInSlice[timeStep][2][index] > 0) MITK_INFO << m_SegmentationCountInSlice[timeStep][2][index]; //"O"; else MITK_INFO << "."; } MITK_INFO << std::endl; } // z indices that are larger than the biggest y index for (unsigned int index = m_SegmentationCountInSlice[timeStep][1].size(); index < m_SegmentationCountInSlice[timeStep][2].size(); ++index) { for (unsigned int indent = 0; indent < index; ++indent) MITK_INFO << " "; if (m_SegmentationCountInSlice[timeStep][2][index] > 0) MITK_INFO << m_SegmentationCountInSlice[timeStep][2][index]; //"O"; else MITK_INFO << "."; MITK_INFO << std::endl; } } mitk::Image::Pointer mitk::SegmentationInterpolationController::Interpolate(unsigned int sliceDimension, unsigned int sliceIndex, const mitk::PlaneGeometry *currentPlane, unsigned int timeStep) { if (m_Segmentation.IsNull()) return nullptr; if (!currentPlane) { return nullptr; } if (timeStep >= m_SegmentationCountInSlice.size()) return nullptr; if (sliceDimension > 2) return nullptr; unsigned int upperLimit = m_SegmentationCountInSlice[timeStep][sliceDimension].size(); if (sliceIndex >= upperLimit - 1) return nullptr; // can't interpolate first and last slice if (sliceIndex < 1) return nullptr; if (m_SegmentationCountInSlice[timeStep][sliceDimension][sliceIndex] > 0) return nullptr; // slice contains a segmentation, won't interpolate anything then unsigned int lowerBound(0); unsigned int upperBound(0); bool bounds(false); for (lowerBound = sliceIndex - 1; /*lowerBound >= 0*/; --lowerBound) { if (m_SegmentationCountInSlice[timeStep][sliceDimension][lowerBound] > 0) { bounds = true; break; } if (lowerBound == 0) break; // otherwise overflow and start at something like 4294967295 } if (!bounds) return nullptr; bounds = false; for (upperBound = sliceIndex + 1; upperBound < upperLimit; ++upperBound) { if (m_SegmentationCountInSlice[timeStep][sliceDimension][upperBound] > 0) { bounds = true; break; } } if (!bounds) return nullptr; // ok, we have found two neighboring slices with segmentations (and we made sure that the current slice does NOT // contain anything // MITK_INFO << "Interpolate in timestep " << timeStep << ", dimension " << sliceDimension << ": estimate slice " << // sliceIndex << " from slices " << lowerBound << " and " << upperBound << std::endl; mitk::Image::Pointer lowerMITKSlice; mitk::Image::Pointer upperMITKSlice; mitk::Image::Pointer resultImage; try { // Setting up the ExtractSliceFilter mitk::ExtractSliceFilter::Pointer extractor = ExtractSliceFilter::New(); extractor->SetInput(m_Segmentation); extractor->SetTimeStep(timeStep); extractor->SetResliceTransformByGeometry(m_Segmentation->GetTimeGeometry()->GetGeometryForTimeStep(timeStep)); extractor->SetVtkOutputRequest(false); // Reslicing the current plane extractor->SetWorldGeometry(currentPlane); extractor->Modified(); extractor->Update(); resultImage = extractor->GetOutput(); resultImage->DisconnectPipeline(); // Creating PlaneGeometry for lower slice mitk::PlaneGeometry::Pointer reslicePlane = currentPlane->Clone(); // Transforming the current origin so that it matches the lower slice mitk::Point3D origin = currentPlane->GetOrigin(); m_Segmentation->GetSlicedGeometry(timeStep)->WorldToIndex(origin, origin); origin[sliceDimension] = lowerBound; m_Segmentation->GetSlicedGeometry(timeStep)->IndexToWorld(origin, origin); reslicePlane->SetOrigin(origin); // Extract the lower slice extractor = ExtractSliceFilter::New(); extractor->SetInput(m_Segmentation); extractor->SetTimeStep(timeStep); extractor->SetResliceTransformByGeometry(m_Segmentation->GetTimeGeometry()->GetGeometryForTimeStep(timeStep)); extractor->SetVtkOutputRequest(false); extractor->SetWorldGeometry(reslicePlane); extractor->Modified(); extractor->Update(); lowerMITKSlice = extractor->GetOutput(); lowerMITKSlice->DisconnectPipeline(); // Transforming the current origin so that it matches the upper slice m_Segmentation->GetSlicedGeometry(timeStep)->WorldToIndex(origin, origin); origin[sliceDimension] = upperBound; m_Segmentation->GetSlicedGeometry(timeStep)->IndexToWorld(origin, origin); reslicePlane->SetOrigin(origin); // Extract the upper slice extractor = ExtractSliceFilter::New(); extractor->SetInput(m_Segmentation); extractor->SetTimeStep(timeStep); extractor->SetResliceTransformByGeometry(m_Segmentation->GetTimeGeometry()->GetGeometryForTimeStep(timeStep)); extractor->SetVtkOutputRequest(false); extractor->SetWorldGeometry(reslicePlane); extractor->Modified(); extractor->Update(); upperMITKSlice = extractor->GetOutput(); upperMITKSlice->DisconnectPipeline(); if (lowerMITKSlice.IsNull() || upperMITKSlice.IsNull()) return nullptr; } catch (const std::exception &e) { MITK_ERROR << "Error in 2D interpolation: " << e.what(); return nullptr; } // interpolation algorithm gets some inputs // two segmentations (guaranteed to be of the same data type, but no special data type guaranteed) // orientation (sliceDimension) of the segmentations // position of the two slices (sliceIndices) // one volume image (original patient image) // // interpolation algorithm can use e.g. itk::ImageSliceConstIteratorWithIndex to // inspect the original patient image at appropriate positions mitk::SegmentationInterpolationAlgorithm::Pointer algorithm = mitk::ShapeBasedInterpolationAlgorithm::New().GetPointer(); return algorithm->Interpolate(lowerMITKSlice.GetPointer(), lowerBound, upperMITKSlice.GetPointer(), upperBound, sliceIndex, sliceDimension, resultImage, timeStep, m_ReferenceImage); } diff --git a/Modules/Segmentation/DataManagement/mitkExtrudedContour.cpp b/Modules/Segmentation/DataManagement/mitkExtrudedContour.cpp index edbcc4dd1a..565f3a275a 100644 --- a/Modules/Segmentation/DataManagement/mitkExtrudedContour.cpp +++ b/Modules/Segmentation/DataManagement/mitkExtrudedContour.cpp @@ -1,387 +1,387 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkExtrudedContour.h" #include "mitkBaseProcess.h" #include "mitkNumericTypes.h" #include "mitkProportionalTimeGeometry.h" #include #include #include #include #include #include #include #include #include #include // vtkButterflySubdivisionFilter * subdivs; #include #include #include #include #include mitk::ExtrudedContour::ExtrudedContour() - : m_Contour(nullptr), m_ClippingGeometry(nullptr), m_AutomaticVectorGeneration(false) + : m_Contour(nullptr), m_ClippingGeometry(nullptr), m_AutomaticVectorGeneration(false), m_Decimate(nullptr) { ProportionalTimeGeometry::Pointer timeGeometry = ProportionalTimeGeometry::New(); timeGeometry->Initialize(1); SetTimeGeometry(timeGeometry); FillVector3D(m_Vector, 0.0, 0.0, 1.0); m_RightVector.Fill(0.0); m_ExtrusionFilter = vtkLinearExtrusionFilter::New(); m_ExtrusionFilter->CappingOff(); m_ExtrusionFilter->SetExtrusionTypeToVectorExtrusion(); double vtkvector[3] = {0, 0, 1}; // set extrusion vector m_ExtrusionFilter->SetVector(vtkvector); m_TriangleFilter = vtkTriangleFilter::New(); m_TriangleFilter->SetInputConnection(m_ExtrusionFilter->GetOutputPort()); m_SubdivisionFilter = vtkLinearSubdivisionFilter::New(); m_SubdivisionFilter->SetInputConnection(m_TriangleFilter->GetOutputPort()); m_SubdivisionFilter->SetNumberOfSubdivisions(4); m_ClippingBox = vtkPlanes::New(); m_ClipPolyDataFilter = vtkClipPolyData::New(); m_ClipPolyDataFilter->SetInputConnection(m_SubdivisionFilter->GetOutputPort()); m_ClipPolyDataFilter->SetClipFunction(m_ClippingBox); m_ClipPolyDataFilter->InsideOutOn(); m_Polygon = vtkPolygon::New(); m_ProjectionPlane = mitk::PlaneGeometry::New(); } mitk::ExtrudedContour::~ExtrudedContour() { m_ClipPolyDataFilter->Delete(); m_ClippingBox->Delete(); m_SubdivisionFilter->Delete(); m_TriangleFilter->Delete(); m_ExtrusionFilter->Delete(); m_Polygon->Delete(); } bool mitk::ExtrudedContour::IsInside(const Point3D &worldPoint) const { static double polygonNormal[3] = {0.0, 0.0, 1.0}; // project point onto plane float xt[3]; itk2vtk(worldPoint, xt); xt[0] = worldPoint[0] - m_Origin[0]; xt[1] = worldPoint[1] - m_Origin[1]; xt[2] = worldPoint[2] - m_Origin[2]; float dist = xt[0] * m_Normal[0] + xt[1] * m_Normal[1] + xt[2] * m_Normal[2]; xt[0] -= dist * m_Normal[0]; xt[1] -= dist * m_Normal[1]; xt[2] -= dist * m_Normal[2]; double x[3]; x[0] = xt[0] * m_Right[0] + xt[1] * m_Right[1] + xt[2] * m_Right[2]; x[1] = xt[0] * m_Down[0] + xt[1] * m_Down[1] + xt[2] * m_Down[2]; x[2] = 0; // determine whether it's in the selection loop and then evaluate point // in polygon only if absolutely necessary. if (x[0] >= this->m_ProjectedContourBounds[0] && x[0] <= this->m_ProjectedContourBounds[1] && x[1] >= this->m_ProjectedContourBounds[2] && x[1] <= this->m_ProjectedContourBounds[3] && this->m_Polygon->PointInPolygon(x, m_Polygon->Points->GetNumberOfPoints(), ((vtkDoubleArray *)this->m_Polygon->Points->GetData())->GetPointer(0), (double *)this->m_ProjectedContourBounds, polygonNormal) == 1) return true; else return false; } mitk::ScalarType mitk::ExtrudedContour::GetVolume() { return -1.0; } void mitk::ExtrudedContour::UpdateOutputInformation() { if (this->GetSource()) { this->GetSource()->UpdateOutputInformation(); } if (GetMTime() > m_LastCalculateExtrusionTime) { BuildGeometry(); BuildSurface(); } // if ( ( m_CalculateBoundingBox ) && ( m_PolyDataSeries.size() > 0 ) ) // CalculateBoundingBox(); } void mitk::ExtrudedContour::BuildSurface() { if (m_Contour.IsNull()) { SetVtkPolyData(nullptr); return; } // set extrusion contour vtkPolyData *polyData = vtkPolyData::New(); vtkCellArray *polys = vtkCellArray::New(); polys->InsertNextCell(m_Polygon->GetPointIds()); polyData->SetPoints(m_Polygon->GetPoints()); // float vtkpoint[3]; // unsigned int i, numPts = m_Polygon->GetNumberOfPoints(); // for(i=0; im_Polygon->Points->GetPoint(i); // pointids[i]=loopPoints->InsertNextPoint(vtkpoint); //} // polys->InsertNextCell( i, pointids ); // delete [] pointids; // polyData->SetPoints( loopPoints ); polyData->SetPolys(polys); polys->Delete(); m_ExtrusionFilter->SetInputData(polyData); polyData->Delete(); // set extrusion scale factor m_ExtrusionFilter->SetScaleFactor(GetGeometry()->GetExtentInMM(2)); SetVtkPolyData(m_SubdivisionFilter->GetOutput()); // if(m_ClippingGeometry.IsNull()) //{ // SetVtkPolyData(m_SubdivisionFilter->GetOutput()); //} // else //{ // m_ClipPolyDataFilter->SetInput(m_SubdivisionFilter->GetOutput()); // mitk::BoundingBox::BoundsArrayType bounds=m_ClippingGeometry->GetBounds(); // m_ClippingBox->SetBounds(bounds[0], bounds[1], bounds[2], bounds[3], bounds[4], bounds[5]); // m_ClippingBox->SetTransform(GetGeometry()->GetVtkTransform()); // m_ClipPolyDataFilter->SetClipFunction(m_ClippingBox); // m_ClipPolyDataFilter->SetValue(0); // SetVtkPolyData(m_ClipPolyDataFilter->GetOutput()); //} m_LastCalculateExtrusionTime.Modified(); } void mitk::ExtrudedContour::BuildGeometry() { if (m_Contour.IsNull()) return; // Initialize(1); Vector3D nullvector; nullvector.Fill(0.0); float xProj[3]; unsigned int i; unsigned int numPts = 20; // m_Contour->GetNumberOfPoints(); mitk::Contour::PathPointer path = m_Contour->GetContourPath(); mitk::Contour::PathType::InputType cstart = path->StartOfInput(); mitk::Contour::PathType::InputType cend = path->EndOfInput(); mitk::Contour::PathType::InputType cstep = (cend - cstart) / numPts; mitk::Contour::PathType::InputType ccur; // Part I: guarantee/calculate legal vectors m_Vector.Normalize(); itk2vtk(m_Vector, m_Normal); // check m_Vector if (mitk::Equal(m_Vector, nullvector) || m_AutomaticVectorGeneration) { if (m_AutomaticVectorGeneration == false) itkWarningMacro("Extrusion vector is 0 (" << m_Vector << "); trying to use normal of polygon"); vtkPoints *loopPoints = vtkPoints::New(); // mitk::Contour::PointsContainerIterator pointsIt = m_Contour->GetPoints()->Begin(); double vtkpoint[3]; unsigned int i = 0; for (i = 0, ccur = cstart; i < numPts; ++i, ccur += cstep) { itk2vtk(path->Evaluate(ccur), vtkpoint); loopPoints->InsertNextPoint(vtkpoint); } // Make sure points define a loop with a m_Normal vtkPolygon::ComputeNormal(loopPoints, m_Normal); loopPoints->Delete(); vtk2itk(m_Normal, m_Vector); if (mitk::Equal(m_Vector, nullvector)) { itkExceptionMacro("Cannot calculate normal of polygon"); } } // check m_RightVector if ((mitk::Equal(m_RightVector, nullvector)) || (mitk::Equal(m_RightVector * m_Vector, 0.0) == false)) { if (mitk::Equal(m_RightVector, nullvector)) { itkDebugMacro("Right vector is 0. Calculating."); } else { itkWarningMacro("Right vector (" << m_RightVector << ") not perpendicular to extrusion vector " << m_Vector << ": " << m_RightVector * m_Vector); } // calculate a legal m_RightVector if (mitk::Equal(m_Vector[1], 0.0f) == false) { FillVector3D(m_RightVector, 1.0f, -m_Vector[0] / m_Vector[1], 0.0f); m_RightVector.Normalize(); } else { FillVector3D(m_RightVector, 0.0f, 1.0f, 0.0f); } } // calculate down-vector VnlVector rightDV = m_RightVector.GetVnlVector(); rightDV.normalize(); vnl2vtk(rightDV, m_Right); VnlVector downDV = vnl_cross_3d(m_Vector.GetVnlVector(), rightDV); downDV.normalize(); vnl2vtk(downDV, m_Down); // Part II: calculate plane as base for extrusion, project the contour // on this plane and store as polygon for IsInside test and BoundingBox calculation // initialize m_ProjectionPlane, yet with origin at 0 m_ProjectionPlane->InitializeStandardPlane(rightDV, downDV); // create vtkPolygon from contour and simultaneously determine 2D bounds of // contour projected on m_ProjectionPlane // mitk::Contour::PointsContainerIterator pointsIt = m_Contour->GetPoints()->Begin(); m_Polygon->Points->Reset(); m_Polygon->Points->SetNumberOfPoints(numPts); m_Polygon->PointIds->Reset(); m_Polygon->PointIds->SetNumberOfIds(numPts); mitk::Point2D pt2d; mitk::Point3D pt3d; mitk::Point2D min, max; min.Fill(itk::NumericTraits::max()); max.Fill(itk::NumericTraits::min()); xProj[2] = 0.0; for (i = 0, ccur = cstart; i < numPts; ++i, ccur += cstep) { pt3d.CastFrom(path->Evaluate(ccur)); m_ProjectionPlane->Map(pt3d, pt2d); xProj[0] = pt2d[0]; if (pt2d[0] < min[0]) min[0] = pt2d[0]; if (pt2d[0] > max[0]) max[0] = pt2d[0]; xProj[1] = pt2d[1]; if (pt2d[1] < min[1]) min[1] = pt2d[1]; if (pt2d[1] > max[1]) max[1] = pt2d[1]; m_Polygon->Points->SetPoint(i, xProj); m_Polygon->PointIds->SetId(i, i); } // shift parametric origin to (0,0) for (i = 0; i < numPts; ++i) { double *pt = this->m_Polygon->Points->GetPoint(i); pt[0] -= min[0]; pt[1] -= min[1]; itkDebugMacro(<< i << ": (" << pt[0] << "," << pt[1] << "," << pt[2] << ")"); } this->m_Polygon->GetBounds(m_ProjectedContourBounds); // m_ProjectedContourBounds[4]=-1.0; m_ProjectedContourBounds[5]=1.0; // calculate origin (except translation along the normal) and bounds // of m_ProjectionPlane: // origin is composed of the minimum x-/y-coordinates of the polygon, // bounds from the extent of the polygon, both after projecting on the plane mitk::Point3D origin; m_ProjectionPlane->Map(min, origin); ScalarType bounds[6] = {0, max[0] - min[0], 0, max[1] - min[1], 0, 1}; m_ProjectionPlane->SetBounds(bounds); m_ProjectionPlane->SetOrigin(origin); // Part III: initialize geometry if (m_ClippingGeometry.IsNotNull()) { ScalarType min_dist = itk::NumericTraits::max(), max_dist = itk::NumericTraits::min(), dist; unsigned char i; for (i = 0; i < 8; ++i) { dist = m_ProjectionPlane->SignedDistance(m_ClippingGeometry->GetCornerPoint(i)); if (dist < min_dist) min_dist = dist; if (dist > max_dist) max_dist = dist; } // incorporate translation along the normal into origin origin = origin + m_Vector * min_dist; m_ProjectionPlane->SetOrigin(origin); bounds[5] = max_dist - min_dist; } else bounds[5] = 20; itk2vtk(origin, m_Origin); mitk::BaseGeometry::Pointer g3d = GetGeometry(0); assert(g3d.IsNotNull()); g3d->SetBounds(bounds); g3d->SetIndexToWorldTransform(m_ProjectionPlane->GetIndexToWorldTransform()); ProportionalTimeGeometry::Pointer timeGeometry = ProportionalTimeGeometry::New(); timeGeometry->Initialize(g3d, 1); SetTimeGeometry(timeGeometry); } unsigned long mitk::ExtrudedContour::GetMTime() const { unsigned long latestTime = Superclass::GetMTime(); if (m_Contour.IsNotNull()) { unsigned long localTime; localTime = m_Contour->GetMTime(); if (localTime > latestTime) latestTime = localTime; } return latestTime; } diff --git a/Modules/Segmentation/Interactions/mitkFastMarchingTool3D.cpp b/Modules/Segmentation/Interactions/mitkFastMarchingTool3D.cpp index ea76fab6e4..a01ad6399a 100644 --- a/Modules/Segmentation/Interactions/mitkFastMarchingTool3D.cpp +++ b/Modules/Segmentation/Interactions/mitkFastMarchingTool3D.cpp @@ -1,462 +1,464 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkFastMarchingTool3D.h" #include "mitkToolManager.h" #include "mitkBaseRenderer.h" #include "mitkInteractionConst.h" #include "mitkRenderingManager.h" #include "itkOrImageFilter.h" #include "mitkImageCast.h" #include "mitkImageTimeSelector.h" // us #include #include #include #include namespace mitk { MITK_TOOL_MACRO(MITKSEGMENTATION_EXPORT, FastMarchingTool3D, "FastMarching3D tool"); } mitk::FastMarchingTool3D::FastMarchingTool3D() : /*FeedbackContourTool*/ AutoSegmentationTool(), m_NeedUpdate(true), m_CurrentTimeStep(0), m_LowerThreshold(0), m_UpperThreshold(200), m_StoppingValue(100), m_Sigma(1.0), m_Alpha(-0.5), - m_Beta(3.0) + m_Beta(3.0), + m_PointSetAddObserverTag(0), + m_PointSetRemoveObserverTag(0) { } mitk::FastMarchingTool3D::~FastMarchingTool3D() { } bool mitk::FastMarchingTool3D::CanHandle(BaseData *referenceData) const { if (referenceData == nullptr) return false; auto *image = dynamic_cast(referenceData); if (image == nullptr) return false; if (image->GetDimension() < 3) return false; return true; } const char **mitk::FastMarchingTool3D::GetXPM() const { return nullptr; // mitkFastMarchingTool3D_xpm; } us::ModuleResource mitk::FastMarchingTool3D::GetIconResource() const { us::Module *module = us::GetModuleContext()->GetModule(); us::ModuleResource resource = module->GetResource("FastMarching_48x48.png"); return resource; } const char *mitk::FastMarchingTool3D::GetName() const { return "Fast Marching 3D"; } void mitk::FastMarchingTool3D::SetUpperThreshold(double value) { m_UpperThreshold = value / 10.0; m_ThresholdFilter->SetUpperThreshold(m_UpperThreshold); m_NeedUpdate = true; } void mitk::FastMarchingTool3D::SetLowerThreshold(double value) { m_LowerThreshold = value / 10.0; m_ThresholdFilter->SetLowerThreshold(m_LowerThreshold); m_NeedUpdate = true; } void mitk::FastMarchingTool3D::SetBeta(double value) { if (m_Beta != value) { m_Beta = value; m_SigmoidFilter->SetBeta(m_Beta); m_NeedUpdate = true; } } void mitk::FastMarchingTool3D::SetSigma(double value) { if (m_Sigma != value) { if (value > 0.0) { m_Sigma = value; m_GradientMagnitudeFilter->SetSigma(m_Sigma); m_NeedUpdate = true; } } } void mitk::FastMarchingTool3D::SetAlpha(double value) { if (m_Alpha != value) { m_Alpha = value; m_SigmoidFilter->SetAlpha(m_Alpha); m_NeedUpdate = true; } } void mitk::FastMarchingTool3D::SetStoppingValue(double value) { if (m_StoppingValue != value) { m_StoppingValue = value; m_FastMarchingFilter->SetStoppingValue(m_StoppingValue); m_NeedUpdate = true; } } void mitk::FastMarchingTool3D::Activated() { Superclass::Activated(); m_ResultImageNode = mitk::DataNode::New(); m_ResultImageNode->SetName("FastMarching_Preview"); m_ResultImageNode->SetBoolProperty("helper object", true); m_ResultImageNode->SetColor(0.0, 1.0, 0.0); m_ResultImageNode->SetVisibility(true); m_ToolManager->GetDataStorage()->Add(this->m_ResultImageNode, m_ToolManager->GetReferenceData(0)); m_SeedsAsPointSet = mitk::PointSet::New(); m_SeedsAsPointSetNode = mitk::DataNode::New(); m_SeedsAsPointSetNode->SetData(m_SeedsAsPointSet); m_SeedsAsPointSetNode->SetName("3D_FastMarching_PointSet"); m_SeedsAsPointSetNode->SetBoolProperty("helper object", true); m_SeedsAsPointSetNode->SetColor(0.0, 1.0, 0.0); m_SeedsAsPointSetNode->SetVisibility(true); // Create PointSetData Interactor m_SeedPointInteractor = mitk::PointSetDataInteractor::New(); // Load the according state machine for regular point set interaction m_SeedPointInteractor->LoadStateMachine("PointSet.xml"); // Set the configuration file that defines the triggers for the transitions m_SeedPointInteractor->SetEventConfig("PointSetConfig.xml"); // set the DataNode (which already is added to the DataStorage m_SeedPointInteractor->SetDataNode(m_SeedsAsPointSetNode); m_ReferenceImageAsITK = InternalImageType::New(); m_ProgressCommand = mitk::ToolCommand::New(); m_ThresholdFilter = ThresholdingFilterType::New(); m_ThresholdFilter->SetLowerThreshold(m_LowerThreshold); m_ThresholdFilter->SetUpperThreshold(m_UpperThreshold); m_ThresholdFilter->SetOutsideValue(0); m_ThresholdFilter->SetInsideValue(1.0); m_SmoothFilter = SmoothingFilterType::New(); m_SmoothFilter->AddObserver(itk::ProgressEvent(), m_ProgressCommand); m_SmoothFilter->SetTimeStep(0.05); m_SmoothFilter->SetNumberOfIterations(2); m_SmoothFilter->SetConductanceParameter(9.0); m_GradientMagnitudeFilter = GradientFilterType::New(); m_GradientMagnitudeFilter->AddObserver(itk::ProgressEvent(), m_ProgressCommand); m_GradientMagnitudeFilter->SetSigma(m_Sigma); m_SigmoidFilter = SigmoidFilterType::New(); m_SigmoidFilter->AddObserver(itk::ProgressEvent(), m_ProgressCommand); m_SigmoidFilter->SetAlpha(m_Alpha); m_SigmoidFilter->SetBeta(m_Beta); m_SigmoidFilter->SetOutputMinimum(0.0); m_SigmoidFilter->SetOutputMaximum(1.0); m_FastMarchingFilter = FastMarchingFilterType::New(); m_FastMarchingFilter->AddObserver(itk::ProgressEvent(), m_ProgressCommand); m_FastMarchingFilter->SetStoppingValue(m_StoppingValue); m_SeedContainer = NodeContainer::New(); m_SeedContainer->Initialize(); m_FastMarchingFilter->SetTrialPoints(m_SeedContainer); // set up pipeline m_SmoothFilter->SetInput(m_ReferenceImageAsITK); m_GradientMagnitudeFilter->SetInput(m_SmoothFilter->GetOutput()); m_SigmoidFilter->SetInput(m_GradientMagnitudeFilter->GetOutput()); m_FastMarchingFilter->SetInput(m_SigmoidFilter->GetOutput()); m_ThresholdFilter->SetInput(m_FastMarchingFilter->GetOutput()); m_ToolManager->GetDataStorage()->Add(m_SeedsAsPointSetNode, m_ToolManager->GetWorkingData(0)); itk::SimpleMemberCommand::Pointer pointAddedCommand = itk::SimpleMemberCommand::New(); pointAddedCommand->SetCallbackFunction(this, &mitk::FastMarchingTool3D::OnAddPoint); m_PointSetAddObserverTag = m_SeedsAsPointSet->AddObserver(mitk::PointSetAddEvent(), pointAddedCommand); itk::SimpleMemberCommand::Pointer pointRemovedCommand = itk::SimpleMemberCommand::New(); pointRemovedCommand->SetCallbackFunction(this, &mitk::FastMarchingTool3D::OnDelete); m_PointSetRemoveObserverTag = m_SeedsAsPointSet->AddObserver(mitk::PointSetRemoveEvent(), pointRemovedCommand); this->Initialize(); } void mitk::FastMarchingTool3D::Deactivated() { m_ToolManager->GetDataStorage()->Remove(this->m_ResultImageNode); m_ToolManager->GetDataStorage()->Remove(this->m_SeedsAsPointSetNode); this->ClearSeeds(); this->m_SmoothFilter->RemoveAllObservers(); this->m_SigmoidFilter->RemoveAllObservers(); this->m_GradientMagnitudeFilter->RemoveAllObservers(); this->m_FastMarchingFilter->RemoveAllObservers(); m_ResultImageNode = nullptr; mitk::RenderingManager::GetInstance()->RequestUpdateAll(); unsigned int numberOfPoints = m_SeedsAsPointSet->GetSize(); for (unsigned int i = 0; i < numberOfPoints; ++i) { mitk::Point3D point = m_SeedsAsPointSet->GetPoint(i); auto *doOp = new mitk::PointOperation(mitk::OpREMOVE, point, 0); m_SeedsAsPointSet->ExecuteOperation(doOp); } // Deactivate Interaction m_SeedPointInteractor->SetDataNode(nullptr); m_ToolManager->GetDataStorage()->Remove(m_SeedsAsPointSetNode); m_SeedsAsPointSetNode = nullptr; m_SeedsAsPointSet->RemoveObserver(m_PointSetAddObserverTag); m_SeedsAsPointSet->RemoveObserver(m_PointSetRemoveObserverTag); Superclass::Deactivated(); } void mitk::FastMarchingTool3D::Initialize() { m_ReferenceImage = dynamic_cast(m_ToolManager->GetReferenceData(0)->GetData()); if (m_ReferenceImage->GetTimeGeometry()->CountTimeSteps() > 1) { mitk::ImageTimeSelector::Pointer timeSelector = ImageTimeSelector::New(); timeSelector->SetInput(m_ReferenceImage); timeSelector->SetTimeNr(m_CurrentTimeStep); timeSelector->UpdateLargestPossibleRegion(); m_ReferenceImage = timeSelector->GetOutput(); } CastToItkImage(m_ReferenceImage, m_ReferenceImageAsITK); m_SmoothFilter->SetInput(m_ReferenceImageAsITK); m_NeedUpdate = true; } void mitk::FastMarchingTool3D::ConfirmSegmentation() { // combine preview image with current working segmentation if (dynamic_cast(m_ResultImageNode->GetData())) { // logical or combination of preview and segmentation slice OutputImageType::Pointer segmentationImageInITK = OutputImageType::New(); mitk::Image::Pointer workingImage = dynamic_cast(GetTargetSegmentationNode()->GetData()); if (workingImage->GetTimeGeometry()->CountTimeSteps() > 1) { mitk::ImageTimeSelector::Pointer timeSelector = mitk::ImageTimeSelector::New(); timeSelector->SetInput(workingImage); timeSelector->SetTimeNr(m_CurrentTimeStep); timeSelector->UpdateLargestPossibleRegion(); CastToItkImage(timeSelector->GetOutput(), segmentationImageInITK); } else { CastToItkImage(workingImage, segmentationImageInITK); } typedef itk::OrImageFilter OrImageFilterType; OrImageFilterType::Pointer orFilter = OrImageFilterType::New(); orFilter->SetInput(0, m_ThresholdFilter->GetOutput()); orFilter->SetInput(1, segmentationImageInITK); orFilter->Update(); // set image volume in current time step from itk image workingImage->SetVolume((void *)(m_ThresholdFilter->GetOutput()->GetPixelContainer()->GetBufferPointer()), m_CurrentTimeStep); this->m_ResultImageNode->SetVisibility(false); this->ClearSeeds(); workingImage->Modified(); } mitk::RenderingManager::GetInstance()->RequestUpdateAll(); m_ToolManager->ActivateTool(-1); } void mitk::FastMarchingTool3D::OnAddPoint() { // Add a new seed point for FastMarching algorithm mitk::Point3D clickInIndex; m_ReferenceImage->GetGeometry()->WorldToIndex(m_SeedsAsPointSet->GetPoint(m_SeedsAsPointSet->GetSize() - 1), clickInIndex); itk::Index<3> seedPosition; seedPosition[0] = clickInIndex[0]; seedPosition[1] = clickInIndex[1]; seedPosition[2] = clickInIndex[2]; NodeType node; const double seedValue = 0.0; node.SetValue(seedValue); node.SetIndex(seedPosition); this->m_SeedContainer->InsertElement(this->m_SeedContainer->Size(), node); m_FastMarchingFilter->Modified(); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); m_NeedUpdate = true; this->Update(); m_ReadyMessage.Send(); } void mitk::FastMarchingTool3D::OnDelete() { // delete last seed point if (!(this->m_SeedContainer->empty())) { // delete last element of seeds container this->m_SeedContainer->pop_back(); m_FastMarchingFilter->Modified(); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); m_NeedUpdate = true; this->Update(); } } void mitk::FastMarchingTool3D::Update() { const unsigned int progress_steps = 200; if (m_NeedUpdate) { m_ProgressCommand->AddStepsToDo(progress_steps); // remove interaction with poinset while updating m_SeedPointInteractor->SetDataNode(nullptr); CurrentlyBusy.Send(true); try { m_ThresholdFilter->Update(); } catch (itk::ExceptionObject &excep) { MITK_ERROR << "Exception caught: " << excep.GetDescription(); m_ProgressCommand->SetProgress(progress_steps); CurrentlyBusy.Send(false); std::string msg = excep.GetDescription(); ErrorMessage.Send(msg); return; } m_ProgressCommand->SetProgress(progress_steps); CurrentlyBusy.Send(false); // make output visible mitk::Image::Pointer result = mitk::Image::New(); CastToMitkImage(m_ThresholdFilter->GetOutput(), result); result->GetGeometry()->SetOrigin(m_ReferenceImage->GetGeometry()->GetOrigin()); result->GetGeometry()->SetIndexToWorldTransform(m_ReferenceImage->GetGeometry()->GetIndexToWorldTransform()); m_ResultImageNode->SetData(result); m_ResultImageNode->SetVisibility(true); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); // add interaction with poinset again m_SeedPointInteractor->SetDataNode(m_SeedsAsPointSetNode); } } void mitk::FastMarchingTool3D::ClearSeeds() { // clear seeds for FastMarching as well as the PointSet for visualization if (this->m_SeedContainer.IsNotNull()) this->m_SeedContainer->Initialize(); if (this->m_SeedsAsPointSet.IsNotNull()) { // remove observers from current pointset m_SeedsAsPointSet->RemoveObserver(m_PointSetAddObserverTag); m_SeedsAsPointSet->RemoveObserver(m_PointSetRemoveObserverTag); // renew pointset this->m_SeedsAsPointSet = mitk::PointSet::New(); this->m_SeedsAsPointSetNode->SetData(this->m_SeedsAsPointSet); m_SeedsAsPointSetNode->SetName("Seeds_Preview"); m_SeedsAsPointSetNode->SetBoolProperty("helper object", true); m_SeedsAsPointSetNode->SetColor(0.0, 1.0, 0.0); m_SeedsAsPointSetNode->SetVisibility(true); // add callback function for adding and removing points itk::SimpleMemberCommand::Pointer pointAddedCommand = itk::SimpleMemberCommand::New(); pointAddedCommand->SetCallbackFunction(this, &mitk::FastMarchingTool3D::OnAddPoint); m_PointSetAddObserverTag = m_SeedsAsPointSet->AddObserver(mitk::PointSetAddEvent(), pointAddedCommand); itk::SimpleMemberCommand::Pointer pointRemovedCommand = itk::SimpleMemberCommand::New(); pointRemovedCommand->SetCallbackFunction(this, &mitk::FastMarchingTool3D::OnDelete); m_PointSetRemoveObserverTag = m_SeedsAsPointSet->AddObserver(mitk::PointSetRemoveEvent(), pointRemovedCommand); } if (this->m_FastMarchingFilter.IsNotNull()) m_FastMarchingFilter->Modified(); this->m_NeedUpdate = true; } void mitk::FastMarchingTool3D::Reset() { // clear all seeds and preview empty result this->ClearSeeds(); m_ResultImageNode->SetVisibility(false); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void mitk::FastMarchingTool3D::SetCurrentTimeStep(int t) { if (m_CurrentTimeStep != t) { m_CurrentTimeStep = t; this->Initialize(); } } diff --git a/Modules/Segmentation/Interactions/mitkLiveWireTool2D.cpp b/Modules/Segmentation/Interactions/mitkLiveWireTool2D.cpp index b5965a9746..4aabc3e86e 100644 --- a/Modules/Segmentation/Interactions/mitkLiveWireTool2D.cpp +++ b/Modules/Segmentation/Interactions/mitkLiveWireTool2D.cpp @@ -1,669 +1,670 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include #include #include #include "mitkLiveWireTool2D.h" #include "mitkLiveWireTool2D.xpm" namespace mitk { MITK_TOOL_MACRO(MITKSEGMENTATION_EXPORT, LiveWireTool2D, "LiveWire tool"); } class RemoveFromDataStorage { public: RemoveFromDataStorage(mitk::DataStorage::Pointer dataStorage) : m_DataStorage(dataStorage) {} void operator()(mitk::DataNode *dataNode) { m_DataStorage->Remove(dataNode); } void operator()(const std::pair &dataNode) { m_DataStorage->Remove(dataNode.first); } private: mitk::DataStorage::Pointer m_DataStorage; }; -mitk::LiveWireTool2D::LiveWireTool2D() : SegTool2D("LiveWireTool"), m_PlaneGeometry(nullptr) +mitk::LiveWireTool2D::LiveWireTool2D() + : SegTool2D("LiveWireTool"), m_PlaneGeometry(nullptr), m_CreateAndUseDynamicCosts(false) { } mitk::LiveWireTool2D::~LiveWireTool2D() { this->ClearSegmentation(); } void mitk::LiveWireTool2D::RemoveHelperObjects() { DataStorage *dataStorage = m_ToolManager->GetDataStorage(); if (!m_EditingContours.empty()) std::for_each(m_EditingContours.begin(), m_EditingContours.end(), RemoveFromDataStorage(dataStorage)); if (!m_WorkingContours.empty()) std::for_each(m_WorkingContours.begin(), m_WorkingContours.end(), RemoveFromDataStorage(dataStorage)); if (m_EditingContourNode.IsNotNull()) dataStorage->Remove(m_EditingContourNode); if (m_LiveWireContourNode.IsNotNull()) dataStorage->Remove(m_LiveWireContourNode); if (m_ContourModelNode.IsNotNull()) dataStorage->Remove(m_ContourModelNode); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void mitk::LiveWireTool2D::ReleaseHelperObjects() { this->RemoveHelperObjects(); if (!m_EditingContours.empty()) m_EditingContours.clear(); if (!m_WorkingContours.empty()) m_WorkingContours.clear(); m_EditingContourNode = nullptr; m_EditingContour = nullptr; m_LiveWireContourNode = nullptr; m_LiveWireContour = nullptr; m_ContourModelNode = nullptr; m_Contour = nullptr; } void mitk::LiveWireTool2D::ReleaseInteractors() { this->EnableContourLiveWireInteraction(false); m_LiveWireNodes.clear(); } void mitk::LiveWireTool2D::ConnectActionsAndFunctions() { CONNECT_CONDITION("CheckContourClosed", OnCheckPoint); CONNECT_FUNCTION("InitObject", OnInitLiveWire); CONNECT_FUNCTION("AddPoint", OnAddPoint); CONNECT_FUNCTION("CtrlAddPoint", OnAddPoint); CONNECT_FUNCTION("MovePoint", OnMouseMoveNoDynamicCosts); CONNECT_FUNCTION("FinishContour", OnFinish); CONNECT_FUNCTION("DeletePoint", OnLastSegmentDelete); CONNECT_FUNCTION("CtrlMovePoint", OnMouseMoved); } const char **mitk::LiveWireTool2D::GetXPM() const { return mitkLiveWireTool2D_xpm; } us::ModuleResource mitk::LiveWireTool2D::GetIconResource() const { return us::GetModuleContext()->GetModule()->GetResource("LiveWire_48x48.png"); } us::ModuleResource mitk::LiveWireTool2D::GetCursorIconResource() const { return us::GetModuleContext()->GetModule()->GetResource("LiveWire_Cursor_32x32.png"); } const char *mitk::LiveWireTool2D::GetName() const { return "Live Wire"; } void mitk::LiveWireTool2D::Activated() { Superclass::Activated(); this->ResetToStartState(); this->EnableContourLiveWireInteraction(true); } void mitk::LiveWireTool2D::Deactivated() { this->ConfirmSegmentation(); Superclass::Deactivated(); } void mitk::LiveWireTool2D::EnableContourLiveWireInteraction(bool on) { for (auto interactor : m_LiveWireNodes) { if (on) interactor->EnableInteraction(true); else interactor->EnableInteraction(false); } } void mitk::LiveWireTool2D::ConfirmSegmentation() { DataNode *workingNode(m_ToolManager->GetWorkingData(0)); if (!workingNode) return; auto *workingImage = dynamic_cast(workingNode->GetData()); if (!workingImage) return; // for all contours in list (currently created by tool) auto itWorkingContours = this->m_WorkingContours.begin(); std::vector sliceList; sliceList.reserve(m_WorkingContours.size()); while (itWorkingContours != this->m_WorkingContours.end()) { // if node contains data if (itWorkingContours->first->GetData()) { // if this is a contourModel auto *contourModel = dynamic_cast(itWorkingContours->first->GetData()); if (contourModel) { // for each timestep of this contourModel for (TimeStepType currentTimestep = 0; currentTimestep < contourModel->GetTimeGeometry()->CountTimeSteps(); ++currentTimestep) { // get the segmentation image slice at current timestep mitk::Image::Pointer workingSlice = this->GetAffectedImageSliceAs2DImage(itWorkingContours->second, workingImage, currentTimestep); mitk::ContourModel::Pointer projectedContour = mitk::ContourModelUtils::ProjectContourTo2DSlice(workingSlice, contourModel, true, false); mitk::ContourModelUtils::FillContourInSlice(projectedContour, workingSlice, workingImage, 1.0); // write back to image volume SliceInformation sliceInfo(workingSlice, itWorkingContours->second, currentTimestep); sliceList.push_back(sliceInfo); this->WriteSliceToVolume(sliceInfo); } } } ++itWorkingContours; } this->WriteBackSegmentationResult(sliceList, false); this->ClearSegmentation(); } void mitk::LiveWireTool2D::ClearSegmentation() { this->ReleaseHelperObjects(); this->ReleaseInteractors(); this->ResetToStartState(); } bool mitk::LiveWireTool2D::IsPositionEventInsideImageRegion(mitk::InteractionPositionEvent *positionEvent, mitk::BaseData *data) { bool IsPositionEventInsideImageRegion = data != nullptr && data->GetGeometry()->IsInside(positionEvent->GetPositionInWorld()); if (!IsPositionEventInsideImageRegion) { MITK_WARN("LiveWireTool2D") << "PositionEvent is outside ImageRegion!"; return false; } return true; } void mitk::LiveWireTool2D::OnInitLiveWire(StateMachineAction *, InteractionEvent *interactionEvent) { auto *positionEvent = dynamic_cast(interactionEvent); if (!positionEvent) return; mitk::DataNode *workingDataNode = m_ToolManager->GetWorkingData(0); if (!IsPositionEventInsideImageRegion(positionEvent, workingDataNode->GetData())) { this->ResetToStartState(); return; } m_LastEventSender = positionEvent->GetSender(); m_LastEventSlice = m_LastEventSender->GetSlice(); int timestep = positionEvent->GetSender()->GetTimeStep(); m_Contour = mitk::ContourModel::New(); m_Contour->Expand(timestep + 1); m_ContourModelNode = mitk::DataNode::New(); m_ContourModelNode->SetData(m_Contour); m_ContourModelNode->SetName("working contour node"); m_ContourModelNode->SetProperty("layer", IntProperty::New(100)); m_ContourModelNode->AddProperty("fixedLayer", BoolProperty::New(true)); m_ContourModelNode->SetProperty("helper object", mitk::BoolProperty::New(true)); m_ContourModelNode->AddProperty("contour.color", ColorProperty::New(1, 1, 0), nullptr, true); m_ContourModelNode->AddProperty("contour.points.color", ColorProperty::New(1.0, 0.0, 0.1), nullptr, true); m_ContourModelNode->AddProperty("contour.controlpoints.show", BoolProperty::New(true), nullptr, true); m_LiveWireContour = mitk::ContourModel::New(); m_LiveWireContour->Expand(timestep + 1); m_LiveWireContourNode = mitk::DataNode::New(); m_LiveWireContourNode->SetData(m_LiveWireContour); m_LiveWireContourNode->SetName("active livewire node"); m_LiveWireContourNode->SetProperty("layer", IntProperty::New(101)); m_LiveWireContourNode->AddProperty("fixedLayer", BoolProperty::New(true)); m_LiveWireContourNode->SetProperty("helper object", mitk::BoolProperty::New(true)); m_LiveWireContourNode->AddProperty("contour.color", ColorProperty::New(0.1, 1.0, 0.1), nullptr, true); m_LiveWireContourNode->AddProperty("contour.width", mitk::FloatProperty::New(4.0), nullptr, true); m_EditingContour = mitk::ContourModel::New(); m_EditingContour->Expand(timestep + 1); m_EditingContourNode = mitk::DataNode::New(); m_EditingContourNode->SetData(m_EditingContour); m_EditingContourNode->SetName("editing node"); m_EditingContourNode->SetProperty("layer", IntProperty::New(102)); m_EditingContourNode->AddProperty("fixedLayer", BoolProperty::New(true)); m_EditingContourNode->SetProperty("helper object", mitk::BoolProperty::New(true)); m_EditingContourNode->AddProperty("contour.color", ColorProperty::New(0.1, 1.0, 0.1), nullptr, true); m_EditingContourNode->AddProperty("contour.points.color", ColorProperty::New(0.0, 0.0, 1.0), nullptr, true); m_EditingContourNode->AddProperty("contour.width", mitk::FloatProperty::New(4.0), nullptr, true); m_ToolManager->GetDataStorage()->Add(m_ContourModelNode, workingDataNode); m_ToolManager->GetDataStorage()->Add(m_LiveWireContourNode, workingDataNode); m_ToolManager->GetDataStorage()->Add(m_EditingContourNode, workingDataNode); // set current slice as input for ImageToLiveWireContourFilter m_WorkingSlice = this->GetAffectedReferenceSlice(positionEvent); mitk::Point3D newOrigin = m_WorkingSlice->GetSlicedGeometry()->GetOrigin(); m_WorkingSlice->GetSlicedGeometry()->WorldToIndex(newOrigin, newOrigin); m_WorkingSlice->GetSlicedGeometry()->IndexToWorld(newOrigin, newOrigin); m_WorkingSlice->GetSlicedGeometry()->SetOrigin(newOrigin); m_LiveWireFilter = mitk::ImageLiveWireContourModelFilter::New(); m_LiveWireFilter->SetInput(m_WorkingSlice); // map click to pixel coordinates mitk::Point3D click = positionEvent->GetPositionInWorld(); itk::Index<3> idx; m_WorkingSlice->GetGeometry()->WorldToIndex(click, idx); // get the pixel the gradient in region of 5x5 itk::Index<3> indexWithHighestGradient; AccessFixedDimensionByItk_2(m_WorkingSlice, FindHighestGradientMagnitudeByITK, 2, idx, indexWithHighestGradient); // itk::Index to mitk::Point3D click[0] = indexWithHighestGradient[0]; click[1] = indexWithHighestGradient[1]; click[2] = indexWithHighestGradient[2]; m_WorkingSlice->GetGeometry()->IndexToWorld(click, click); // set initial start point m_Contour->AddVertex(click, true, timestep); m_LiveWireFilter->SetStartPoint(click); // remember plane geometry to determine if events were triggered in same plane m_PlaneGeometry = interactionEvent->GetSender()->GetCurrentWorldPlaneGeometry(); m_CreateAndUseDynamicCosts = true; // render assert(positionEvent->GetSender()->GetRenderWindow()); mitk::RenderingManager::GetInstance()->RequestUpdate(positionEvent->GetSender()->GetRenderWindow()); } void mitk::LiveWireTool2D::OnAddPoint(StateMachineAction *, InteractionEvent *interactionEvent) { // complete LiveWire interaction for last segment // add current LiveWire contour to the finished contour and reset // to start new segment and computation auto *positionEvent = dynamic_cast(interactionEvent); if (!positionEvent) return; if (m_PlaneGeometry != nullptr) { // this checks that the point is in the correct slice if (m_PlaneGeometry->DistanceFromPlane(positionEvent->GetPositionInWorld()) > mitk::sqrteps) { return; } } int timestep = positionEvent->GetSender()->GetTimeStep(); // add repulsive points to avoid to get the same path again typedef mitk::ImageLiveWireContourModelFilter::InternalImageType::IndexType IndexType; mitk::ContourModel::ConstVertexIterator iter = m_LiveWireContour->IteratorBegin(timestep); for (; iter != m_LiveWireContour->IteratorEnd(timestep); iter++) { IndexType idx; this->m_WorkingSlice->GetGeometry()->WorldToIndex((*iter)->Coordinates, idx); this->m_LiveWireFilter->AddRepulsivePoint(idx); } // remove duplicate first vertex, it's already contained in m_Contour m_LiveWireContour->RemoveVertexAt(0, timestep); // set last added point as control point m_LiveWireContour->SetControlVertexAt(m_LiveWireContour->GetNumberOfVertices(timestep) - 1, timestep); // merge contours m_Contour->Concatenate(m_LiveWireContour, timestep); // clear the livewire contour and reset the corresponding datanode m_LiveWireContour->Clear(timestep); // set new start point m_LiveWireFilter->SetStartPoint(positionEvent->GetPositionInWorld()); if (m_CreateAndUseDynamicCosts) { // use dynamic cost map for next update m_LiveWireFilter->CreateDynamicCostMap(m_Contour); m_LiveWireFilter->SetUseDynamicCostMap(true); // m_CreateAndUseDynamicCosts = false; } // render assert(positionEvent->GetSender()->GetRenderWindow()); mitk::RenderingManager::GetInstance()->RequestUpdate(positionEvent->GetSender()->GetRenderWindow()); } void mitk::LiveWireTool2D::OnMouseMoved(StateMachineAction *, InteractionEvent *interactionEvent) { // compute LiveWire segment from last control point to current mouse position auto *positionEvent = dynamic_cast(interactionEvent); if (!positionEvent) return; // actual LiveWire computation int timestep = positionEvent->GetSender()->GetTimeStep(); m_LiveWireFilter->SetEndPoint(positionEvent->GetPositionInWorld()); m_LiveWireFilter->SetTimeStep(timestep); m_LiveWireFilter->Update(); m_LiveWireContour = this->m_LiveWireFilter->GetOutput(); m_LiveWireContourNode->SetData(this->m_LiveWireContour); // render assert(positionEvent->GetSender()->GetRenderWindow()); positionEvent->GetSender()->GetRenderingManager()->RequestUpdate(positionEvent->GetSender()->GetRenderWindow()); } void mitk::LiveWireTool2D::OnMouseMoveNoDynamicCosts(StateMachineAction *, InteractionEvent *interactionEvent) { // do not use dynamic cost map m_LiveWireFilter->SetUseDynamicCostMap(false); OnMouseMoved(nullptr, interactionEvent); m_LiveWireFilter->SetUseDynamicCostMap(true); } bool mitk::LiveWireTool2D::OnCheckPoint(const InteractionEvent *interactionEvent) { // check double click on first control point to finish the LiveWire tool // // Check distance to first point. // Transition YES if click close to first control point // const auto *positionEvent = dynamic_cast(interactionEvent); if (positionEvent) { int timestep = positionEvent->GetSender()->GetTimeStep(); mitk::Point3D click = positionEvent->GetPositionInWorld(); mitk::Point3D first = this->m_Contour->GetVertexAt(0, timestep)->Coordinates; if (first.EuclideanDistanceTo(click) < 4.5) { // allow to finish return true; } else { return false; } } return false; } void mitk::LiveWireTool2D::OnFinish(StateMachineAction *, InteractionEvent *interactionEvent) { // finish livewire tool interaction auto *positionEvent = dynamic_cast(interactionEvent); if (!positionEvent) return; // Have to do that here so that the m_LastEventSender is set correctly mitk::SegTool2D::AddContourmarker(); // actual timestep int timestep = positionEvent->GetSender()->GetTimeStep(); // remove last control point being added by double click m_Contour->RemoveVertexAt(m_Contour->GetNumberOfVertices(timestep) - 1, timestep); // save contour and corresponding plane geometry to list std::pair cp( m_ContourModelNode, (positionEvent->GetSender()->GetCurrentWorldPlaneGeometry()->Clone().GetPointer())); this->m_WorkingContours.push_back(cp); std::pair ecp( m_EditingContourNode, (positionEvent->GetSender()->GetCurrentWorldPlaneGeometry()->Clone().GetPointer())); this->m_EditingContours.push_back(ecp); m_LiveWireFilter->SetUseDynamicCostMap(false); this->FinishTool(); } void mitk::LiveWireTool2D::FinishTool() { auto numberOfTimesteps = static_cast(m_Contour->GetTimeGeometry()->CountTimeSteps()); // close contour in each timestep for (int i = 0; i <= numberOfTimesteps; i++) m_Contour->Close(i); m_ToolManager->GetDataStorage()->Remove(m_LiveWireContourNode); // clear live wire contour node m_LiveWireContourNode = nullptr; m_LiveWireContour = nullptr; // A new ContourModelLiveWireInteractor is created that will listen to new events // set the livewire interactor to edit control points m_ContourInteractor = mitk::ContourModelLiveWireInteractor::New(); m_ContourInteractor->SetDataNode(m_ContourModelNode); // TODO load statemachine and config m_ContourInteractor->LoadStateMachine("ContourModelModificationInteractor.xml", us::GetModuleContext()->GetModule()); // Set the configuration file that defines the triggers for the transitions m_ContourInteractor->SetEventConfig("ContourModelModificationConfig.xml", us::GetModuleContext()->GetModule()); m_ContourInteractor->SetWorkingImage(this->m_WorkingSlice); m_ContourInteractor->SetEditingContourModelNode(this->m_EditingContourNode); m_ContourModelNode->SetDataInteractor(m_ContourInteractor.GetPointer()); this->m_LiveWireNodes.push_back(m_ContourInteractor); } void mitk::LiveWireTool2D::OnLastSegmentDelete(StateMachineAction *, InteractionEvent *interactionEvent) { int timestep = interactionEvent->GetSender()->GetTimeStep(); // if last point of current contour will be removed go to start state and remove nodes if (m_Contour->GetNumberOfVertices(timestep) <= 1) { m_ToolManager->GetDataStorage()->Remove(m_LiveWireContourNode); m_ToolManager->GetDataStorage()->Remove(m_ContourModelNode); m_ToolManager->GetDataStorage()->Remove(m_EditingContourNode); m_LiveWireContour = mitk::ContourModel::New(); m_Contour = mitk::ContourModel::New(); m_ContourModelNode->SetData(m_Contour); m_LiveWireContourNode->SetData(m_LiveWireContour); this->ResetToStartState(); // go to start state } else // remove last segment from contour and reset livewire contour { m_LiveWireContour = mitk::ContourModel::New(); m_LiveWireContourNode->SetData(m_LiveWireContour); mitk::ContourModel::Pointer newContour = mitk::ContourModel::New(); newContour->Expand(m_Contour->GetTimeSteps()); auto begin = m_Contour->IteratorBegin(); // iterate from last point to next active point auto newLast = m_Contour->IteratorBegin() + (m_Contour->GetNumberOfVertices() - 1); // go at least one down if (newLast != begin) { newLast--; } // search next active control point while (newLast != begin && !((*newLast)->IsControlPoint)) { newLast--; } // set position of start point for livewire filter to coordinates of the new last point m_LiveWireFilter->SetStartPoint((*newLast)->Coordinates); auto it = m_Contour->IteratorBegin(); // fill new Contour while (it <= newLast) { newContour->AddVertex((*it)->Coordinates, (*it)->IsControlPoint, timestep); it++; } newContour->SetClosed(m_Contour->IsClosed()); // set new contour visible m_ContourModelNode->SetData(newContour); m_Contour = newContour; assert(interactionEvent->GetSender()->GetRenderWindow()); mitk::RenderingManager::GetInstance()->RequestUpdate(interactionEvent->GetSender()->GetRenderWindow()); } } template void mitk::LiveWireTool2D::FindHighestGradientMagnitudeByITK(itk::Image *inputImage, itk::Index<3> &index, itk::Index<3> &returnIndex) { typedef itk::Image InputImageType; typedef typename InputImageType::IndexType IndexType; unsigned long xMAX = inputImage->GetLargestPossibleRegion().GetSize()[0]; unsigned long yMAX = inputImage->GetLargestPossibleRegion().GetSize()[1]; returnIndex[0] = index[0]; returnIndex[1] = index[1]; returnIndex[2] = 0.0; double gradientMagnitude = 0.0; double maxGradientMagnitude = 0.0; /* the size and thus the region of 7x7 is only used to calculate the gradient magnitude in that region not for searching the maximum value */ // maximum value in each direction for size typename InputImageType::SizeType size; size[0] = 7; size[1] = 7; // minimum value in each direction for startRegion IndexType startRegion; startRegion[0] = index[0] - 3; startRegion[1] = index[1] - 3; if (startRegion[0] < 0) startRegion[0] = 0; if (startRegion[1] < 0) startRegion[1] = 0; if (xMAX - index[0] < 7) startRegion[0] = xMAX - 7; if (yMAX - index[1] < 7) startRegion[1] = yMAX - 7; index[0] = startRegion[0] + 3; index[1] = startRegion[1] + 3; typename InputImageType::RegionType region; region.SetSize(size); region.SetIndex(startRegion); typedef typename itk::GradientMagnitudeImageFilter GradientMagnitudeFilterType; typename GradientMagnitudeFilterType::Pointer gradientFilter = GradientMagnitudeFilterType::New(); gradientFilter->SetInput(inputImage); gradientFilter->GetOutput()->SetRequestedRegion(region); gradientFilter->Update(); typename InputImageType::Pointer gradientMagnImage; gradientMagnImage = gradientFilter->GetOutput(); IndexType currentIndex; currentIndex[0] = 0; currentIndex[1] = 0; // search max (approximate) gradient magnitude for (int x = -1; x <= 1; ++x) { currentIndex[0] = index[0] + x; for (int y = -1; y <= 1; ++y) { currentIndex[1] = index[1] + y; gradientMagnitude = gradientMagnImage->GetPixel(currentIndex); // check for new max if (maxGradientMagnitude < gradientMagnitude) { maxGradientMagnitude = gradientMagnitude; returnIndex[0] = currentIndex[0]; returnIndex[1] = currentIndex[1]; returnIndex[2] = 0.0; } // end if } // end for y currentIndex[1] = index[1]; } // end for x } diff --git a/Modules/Segmentation/Interactions/mitkRegionGrowingTool.cpp b/Modules/Segmentation/Interactions/mitkRegionGrowingTool.cpp index 82174afa8b..28e42ef5a7 100644 --- a/Modules/Segmentation/Interactions/mitkRegionGrowingTool.cpp +++ b/Modules/Segmentation/Interactions/mitkRegionGrowingTool.cpp @@ -1,654 +1,655 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkRegionGrowingTool.h" #include "mitkApplicationCursor.h" #include "mitkBaseRenderer.h" #include "mitkImageDataItem.h" #include "mitkImageToContourModelFilter.h" #include "mitkOverwriteSliceImageFilter.h" #include "mitkRegionGrowingTool.xpm" #include "mitkRenderingManager.h" #include "mitkToolManager.h" #include "mitkExtractDirectedPlaneImageFilterNew.h" #include "mitkLabelSetImage.h" #include "mitkOverwriteDirectedPlaneImageFilter.h" // us #include #include #include #include // ITK #include "mitkITKImageImport.h" #include "mitkImageAccessByItk.h" #include #include #include #include #include namespace mitk { MITK_TOOL_MACRO(MITKSEGMENTATION_EXPORT, RegionGrowingTool, "Region growing tool"); } #define ROUND(a) ((a) > 0 ? (int)((a) + 0.5) : -(int)(0.5 - (a))) mitk::RegionGrowingTool::RegionGrowingTool() : FeedbackContourTool("PressMoveRelease"), m_SeedValue(0), m_ScreenYDifference(0), m_ScreenXDifference(0), m_MouseDistanceScaleFactor(0.5), m_FillFeedbackContour(true), - m_ConnectedComponentValue(1) + m_ConnectedComponentValue(1), + m_PaintingPixelValue(0) { } mitk::RegionGrowingTool::~RegionGrowingTool() { } void mitk::RegionGrowingTool::ConnectActionsAndFunctions() { CONNECT_FUNCTION("PrimaryButtonPressed", OnMousePressed); CONNECT_FUNCTION("Move", OnMouseMoved); CONNECT_FUNCTION("Release", OnMouseReleased); } const char **mitk::RegionGrowingTool::GetXPM() const { return mitkRegionGrowingTool_xpm; } us::ModuleResource mitk::RegionGrowingTool::GetIconResource() const { us::Module *module = us::GetModuleContext()->GetModule(); us::ModuleResource resource = module->GetResource("RegionGrowing_48x48.png"); return resource; } us::ModuleResource mitk::RegionGrowingTool::GetCursorIconResource() const { us::Module *module = us::GetModuleContext()->GetModule(); us::ModuleResource resource = module->GetResource("RegionGrowing_Cursor_32x32.png"); return resource; } const char *mitk::RegionGrowingTool::GetName() const { return "Region Growing"; } void mitk::RegionGrowingTool::Activated() { Superclass::Activated(); } void mitk::RegionGrowingTool::Deactivated() { Superclass::Deactivated(); } // Get the average pixel value of square/cube with radius=neighborhood around index template void mitk::RegionGrowingTool::GetNeighborhoodAverage(itk::Image *itkImage, itk::Index index, ScalarType *result, unsigned int neighborhood) { // maybe assert that image dimension is only 2 or 3? auto neighborhoodInt = (int)neighborhood; TPixel averageValue(0); unsigned int numberOfPixels = (2 * neighborhood + 1) * (2 * neighborhood + 1); if (imageDimension == 3) { numberOfPixels *= (2 * neighborhood + 1); } MITK_DEBUG << "Getting neighborhood of " << numberOfPixels << " pixels around " << index; itk::Index currentIndex; for (int i = (0 - neighborhoodInt); i <= neighborhoodInt; ++i) { currentIndex[0] = index[0] + i; for (int j = (0 - neighborhoodInt); j <= neighborhoodInt; ++j) { currentIndex[1] = index[1] + j; if (imageDimension == 3) { for (int k = (0 - neighborhoodInt); k <= neighborhoodInt; ++k) { currentIndex[2] = index[2] + k; if (itkImage->GetLargestPossibleRegion().IsInside(currentIndex)) { averageValue += itkImage->GetPixel(currentIndex); } else { numberOfPixels -= 1; } } } else { if (itkImage->GetLargestPossibleRegion().IsInside(currentIndex)) { averageValue += itkImage->GetPixel(currentIndex); } else { numberOfPixels -= 1; } } } } *result = (ScalarType)averageValue; *result /= numberOfPixels; } // Check whether index lies inside a segmentation template void mitk::RegionGrowingTool::IsInsideSegmentation(itk::Image *itkImage, itk::Index index, bool *result) { if (itkImage->GetPixel(index) > 0) { *result = true; } else { *result = false; } } // Do the region growing (i.e. call an ITK filter that does it) template void mitk::RegionGrowingTool::StartRegionGrowing(itk::Image *inputImage, itk::Index seedIndex, std::array thresholds, mitk::Image::Pointer &outputImage) { MITK_DEBUG << "Starting region growing at index " << seedIndex << " with lower threshold " << thresholds[0] << " and upper threshold " << thresholds[1]; typedef itk::Image InputImageType; typedef itk::Image OutputImageType; typedef itk::ConnectedThresholdImageFilter RegionGrowingFilterType; typename RegionGrowingFilterType::Pointer regionGrower = RegionGrowingFilterType::New(); // perform region growing in desired segmented region regionGrower->SetInput(inputImage); regionGrower->AddSeed(seedIndex); regionGrower->SetLower(thresholds[0]); regionGrower->SetUpper(thresholds[1]); try { regionGrower->Update(); } catch (...) { return; // Should we do something? } typename OutputImageType::Pointer resultImage = regionGrower->GetOutput(); // Smooth result: Every pixel is replaced by the majority of the neighborhood typedef itk::NeighborhoodIterator NeighborhoodIteratorType; typedef itk::ImageRegionIterator ImageIteratorType; typename NeighborhoodIteratorType::RadiusType radius; radius.Fill(2); // for now, maybe make this something the user can adjust in the preferences? typedef itk::ImageDuplicator< OutputImageType > DuplicatorType; typename DuplicatorType::Pointer duplicator = DuplicatorType::New(); duplicator->SetInputImage(resultImage); duplicator->Update(); typename OutputImageType::Pointer resultDup = duplicator->GetOutput(); NeighborhoodIteratorType neighborhoodIterator(radius, resultDup, resultDup->GetRequestedRegion()); ImageIteratorType imageIterator(resultImage, resultImage->GetRequestedRegion()); for (neighborhoodIterator.GoToBegin(), imageIterator.GoToBegin(); !neighborhoodIterator.IsAtEnd(); ++neighborhoodIterator, ++imageIterator) { DefaultSegmentationDataType voteYes(0); DefaultSegmentationDataType voteNo(0); for (unsigned int i = 0; i < neighborhoodIterator.Size(); ++i) { if (neighborhoodIterator.GetPixel(i) > 0) { voteYes += 1; } else { voteNo += 1; } } if (voteYes > voteNo) { imageIterator.Set(1); } else { imageIterator.Set(0); } } if (resultImage.IsNull()) { MITK_DEBUG << "Region growing result is empty."; } // Can potentially have multiple regions, use connected component image filter to label disjunct regions typedef itk::ConnectedComponentImageFilter ConnectedComponentImageFilterType; typename ConnectedComponentImageFilterType::Pointer connectedComponentFilter = ConnectedComponentImageFilterType::New(); connectedComponentFilter->SetInput(resultImage); connectedComponentFilter->Update(); typename OutputImageType::Pointer resultImageCC = connectedComponentFilter->GetOutput(); m_ConnectedComponentValue = resultImageCC->GetPixel(seedIndex); outputImage = mitk::GrabItkImageMemory(resultImageCC); } void mitk::RegionGrowingTool::OnMousePressed(StateMachineAction *, InteractionEvent *interactionEvent) { auto *positionEvent = dynamic_cast(interactionEvent); if (!positionEvent) return; MITK_DEBUG << "OnMousePressed"; m_LastEventSender = positionEvent->GetSender(); m_LastEventSlice = m_LastEventSender->GetSlice(); m_LastScreenPosition = positionEvent->GetPointerPositionOnScreen(); // ReferenceSlice is from the underlying image, WorkingSlice from the active segmentation (can be empty) m_ReferenceSlice = FeedbackContourTool::GetAffectedReferenceSlice(positionEvent); m_WorkingSlice = FeedbackContourTool::GetAffectedWorkingSlice(positionEvent); if (m_WorkingSlice.IsNotNull()) // can't do anything without a working slice (i.e. a possibly empty segmentation) { MITK_DEBUG << "OnMousePressed: got working slice"; // 2. Determine if the user clicked inside or outside of the segmentation/working slice (i.e. the whole volume) mitk::BaseGeometry::Pointer workingSliceGeometry; workingSliceGeometry = m_WorkingSlice->GetGeometry(); workingSliceGeometry->WorldToIndex(positionEvent->GetPositionInWorld(), m_SeedPoint); itk::Index<2> indexInWorkingSlice2D; indexInWorkingSlice2D[0] = m_SeedPoint[0]; indexInWorkingSlice2D[1] = m_SeedPoint[1]; if (workingSliceGeometry->IsIndexInside(m_SeedPoint)) { MITK_DEBUG << "OnMousePressed: point " << positionEvent->GetPositionInWorld() << " (index coordinates " << m_SeedPoint << ") is inside working slice"; // 3. determine the pixel value under the last click to determine what to do bool inside(true); AccessFixedDimensionByItk_2(m_WorkingSlice, IsInsideSegmentation, 2, indexInWorkingSlice2D, &inside); m_PaintingPixelValue = inside ? 0 : 1; if (inside) { MITK_DEBUG << "Clicked inside segmentation"; // For now, we're doing nothing when the user clicks inside the segmentation. Behaviour can be implemented via // OnMousePressedInside() // When you do, be sure to remove the m_PaintingPixelValue check in OnMouseMoved() and OnMouseReleased() return; } else { MITK_DEBUG << "Clicked outside of segmentation"; OnMousePressedOutside(nullptr, interactionEvent); } } } } // Use this to implement a behaviour for when the user clicks inside a segmentation (for example remove something) // Old IpPic code is kept as comment for reference void mitk::RegionGrowingTool::OnMousePressedInside() { // mitk::InteractionPositionEvent* positionEvent = dynamic_cast( interactionEvent // ); // //const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent()); // checked in // OnMousePressed // // 3.1.1. Create a skeletonization of the segmentation and try to find a nice cut // // apply the skeletonization-and-cut algorithm // // generate contour to remove // // set m_ReferenceSlice = nullptr so nothing will happen during mouse move // // remember to fill the contour with 0 in mouserelease // mitkIpPicDescriptor* segmentationHistory = ipMITKSegmentationCreateGrowerHistory( workingPicSlice, // m_LastWorkingSeed, nullptr ); // free again // if (segmentationHistory) // { // tCutResult cutContour = ipMITKSegmentationGetCutPoints( workingPicSlice, segmentationHistory, // initialWorkingOffset ); // tCutResult is a ipSegmentation type // mitkIpPicFree( segmentationHistory ); // if (cutContour.cutIt) // { // int timestep = positionEvent->GetSender()->GetTimeStep(); // // 3.1.2 copy point from float* to mitk::Contour // ContourModel::Pointer contourInImageIndexCoordinates = ContourModel::New(); // contourInImageIndexCoordinates->Expand(timestep + 1); // contourInImageIndexCoordinates->SetClosed(true, timestep); // Point3D newPoint; // for (int index = 0; index < cutContour.deleteSize; ++index) // { // newPoint[0] = cutContour.deleteCurve[ 2 * index + 0 ] - 0.5;//correction is needed because the // output of the algorithm is center based // newPoint[1] = cutContour.deleteCurve[ 2 * index + 1 ] - 0.5;//and we want our contour displayed // corner based. // newPoint[2] = 0.0; // contourInImageIndexCoordinates->AddVertex( newPoint, timestep ); // } // free(cutContour.traceline); // free(cutContour.deleteCurve); // perhaps visualize this for fun? // free(cutContour.onGradient); // ContourModel::Pointer contourInWorldCoordinates = FeedbackContourTool::BackProjectContourFrom2DSlice( // m_WorkingSlice->GetGeometry(), contourInImageIndexCoordinates, true ); // true: sub 0.5 for // ipSegmentation correction // FeedbackContourTool::SetFeedbackContour( contourInWorldCoordinates ); // FeedbackContourTool::SetFeedbackContourVisible(true); // mitk::RenderingManager::GetInstance()->RequestUpdate( positionEvent->GetSender()->GetRenderWindow() ); // m_FillFeedbackContour = true; // } // else // { // m_FillFeedbackContour = false; // } // } // else // { // m_FillFeedbackContour = false; // } // m_ReferenceSlice = nullptr; // return true; } void mitk::RegionGrowingTool::OnMousePressedOutside(StateMachineAction *, InteractionEvent *interactionEvent) { auto *positionEvent = dynamic_cast(interactionEvent); if (positionEvent) { // Get geometry and indices mitk::BaseGeometry::Pointer workingSliceGeometry; workingSliceGeometry = m_WorkingSlice->GetGeometry(); itk::Index<2> indexInWorkingSlice2D; indexInWorkingSlice2D[0] = m_SeedPoint[0]; indexInWorkingSlice2D[1] = m_SeedPoint[1]; mitk::BaseGeometry::Pointer referenceSliceGeometry; referenceSliceGeometry = m_ReferenceSlice->GetGeometry(); itk::Index<3> indexInReferenceSlice; itk::Index<2> indexInReferenceSlice2D; referenceSliceGeometry->WorldToIndex(positionEvent->GetPositionInWorld(), indexInReferenceSlice); indexInReferenceSlice2D[0] = indexInReferenceSlice[0]; indexInReferenceSlice2D[1] = indexInReferenceSlice[1]; // Get seed neighborhood ScalarType averageValue(0); AccessFixedDimensionByItk_3(m_ReferenceSlice, GetNeighborhoodAverage, 2, indexInReferenceSlice2D, &averageValue, 1); m_SeedValue = averageValue; MITK_DEBUG << "Seed value is " << m_SeedValue; // Get level window settings LevelWindow lw(0, 500); // default window 0 to 500, can we do something smarter here? m_ToolManager->GetReferenceData(0)->GetLevelWindow( lw); // will fill lw if levelwindow property is present, otherwise won't touch it. ScalarType currentVisibleWindow = lw.GetWindow(); MITK_DEBUG << "Level window width is " << currentVisibleWindow; m_InitialThresholds[0] = m_SeedValue - currentVisibleWindow / 20.0; // 20 is arbitrary (though works reasonably // well), is there a better alternative (maybe // option in preferences)? m_InitialThresholds[1] = m_SeedValue + currentVisibleWindow / 20.0; m_Thresholds[0] = m_InitialThresholds[0]; m_Thresholds[1] = m_InitialThresholds[1]; // Perform region growing mitk::Image::Pointer resultImage = mitk::Image::New(); AccessFixedDimensionByItk_3( m_ReferenceSlice, StartRegionGrowing, 2, indexInWorkingSlice2D, m_Thresholds, resultImage); resultImage->SetGeometry(workingSliceGeometry); // Extract contour if (resultImage.IsNotNull() && m_ConnectedComponentValue >= 1) { float isoOffset = 0.33; mitk::ImageToContourModelFilter::Pointer contourExtractor = mitk::ImageToContourModelFilter::New(); contourExtractor->SetInput(resultImage); contourExtractor->SetContourValue(m_ConnectedComponentValue - isoOffset); contourExtractor->Update(); ContourModel::Pointer resultContour = ContourModel::New(); resultContour = contourExtractor->GetOutput(); // Show contour if (resultContour.IsNotNull()) { ContourModel::Pointer resultContourWorld = FeedbackContourTool::BackProjectContourFrom2DSlice( workingSliceGeometry, FeedbackContourTool::ProjectContourTo2DSlice(m_WorkingSlice, resultContour)); // this is not a beautiful solution, just one that works, check T22412 for details int timestep = positionEvent->GetSender()->GetTimeStep(); if (0 != timestep) { int size = resultContourWorld->GetNumberOfVertices(0); auto resultContourTimeWorld = mitk::ContourModel::New(); resultContourTimeWorld->Expand(timestep + 1); for (int loop = 0; loop < size; ++loop) { resultContourTimeWorld->AddVertex(resultContourWorld->GetVertexAt(loop, 0), timestep); } FeedbackContourTool::SetFeedbackContour(resultContourTimeWorld); } else { FeedbackContourTool::SetFeedbackContour(resultContourWorld); } FeedbackContourTool::SetFeedbackContourVisible(true); mitk::RenderingManager::GetInstance()->RequestUpdate(m_LastEventSender->GetRenderWindow()); } } } } void mitk::RegionGrowingTool::OnMouseMoved(StateMachineAction *, InteractionEvent *interactionEvent) { // Until OnMousePressedInside() implements a behaviour, we're just returning here whenever m_PaintingPixelValue is 0, // i.e. when the user clicked inside the segmentation if (m_PaintingPixelValue == 0) { return; } auto *positionEvent = dynamic_cast(interactionEvent); if (m_ReferenceSlice.IsNotNull() && positionEvent) { // Get geometry and indices mitk::BaseGeometry::Pointer workingSliceGeometry; workingSliceGeometry = m_WorkingSlice->GetGeometry(); itk::Index<2> indexInWorkingSlice2D; indexInWorkingSlice2D[0] = m_SeedPoint[0]; indexInWorkingSlice2D[1] = m_SeedPoint[1]; m_ScreenYDifference += positionEvent->GetPointerPositionOnScreen()[1] - m_LastScreenPosition[1]; m_ScreenXDifference += positionEvent->GetPointerPositionOnScreen()[0] - m_LastScreenPosition[0]; m_LastScreenPosition = positionEvent->GetPointerPositionOnScreen(); // Moving the mouse up and down adjusts the width of the threshold window, moving it left and right shifts the // threshold window m_Thresholds[0] = std::min( m_SeedValue, m_InitialThresholds[0] - (m_ScreenYDifference - m_ScreenXDifference) * m_MouseDistanceScaleFactor); m_Thresholds[1] = std::max( m_SeedValue, m_InitialThresholds[1] + (m_ScreenYDifference + m_ScreenXDifference) * m_MouseDistanceScaleFactor); MITK_DEBUG << "Screen difference X: " << m_ScreenXDifference; // Perform region growing again and show the result mitk::Image::Pointer resultImage = mitk::Image::New(); AccessFixedDimensionByItk_3( m_ReferenceSlice, StartRegionGrowing, 2, indexInWorkingSlice2D, m_Thresholds, resultImage); resultImage->SetGeometry(workingSliceGeometry); // Update the contour if (resultImage.IsNotNull() && m_ConnectedComponentValue >= 1) { float isoOffset = 0.33; mitk::ImageToContourModelFilter::Pointer contourExtractor = mitk::ImageToContourModelFilter::New(); contourExtractor->SetInput(resultImage); contourExtractor->SetContourValue(m_ConnectedComponentValue - isoOffset); contourExtractor->Update(); ContourModel::Pointer resultContour = ContourModel::New(); resultContour = contourExtractor->GetOutput(); // Show contour if (resultContour.IsNotNull()) { ContourModel::Pointer resultContourWorld = FeedbackContourTool::BackProjectContourFrom2DSlice( workingSliceGeometry, FeedbackContourTool::ProjectContourTo2DSlice(m_WorkingSlice, resultContour)); // this is not a beautiful solution, just one that works, check T22412 for details int timestep = positionEvent->GetSender()->GetTimeStep(); if (0 != timestep) { int size = resultContourWorld->GetNumberOfVertices(0); auto resultContourTimeWorld = mitk::ContourModel::New(); resultContourTimeWorld->Expand(timestep + 1); for (int loop = 0; loop < size; ++loop) { resultContourTimeWorld->AddVertex(resultContourWorld->GetVertexAt(loop, 0), timestep); } FeedbackContourTool::SetFeedbackContour(resultContourTimeWorld); } else { FeedbackContourTool::SetFeedbackContour(resultContourWorld); } FeedbackContourTool::SetFeedbackContourVisible(true); mitk::RenderingManager::GetInstance()->ForceImmediateUpdate(positionEvent->GetSender()->GetRenderWindow()); } } } } void mitk::RegionGrowingTool::OnMouseReleased(StateMachineAction *, InteractionEvent *interactionEvent) { // Until OnMousePressedInside() implements a behaviour, we're just returning here whenever m_PaintingPixelValue is 0, // i.e. when the user clicked inside the segmentation if (m_PaintingPixelValue == 0) { return; } auto *positionEvent = dynamic_cast(interactionEvent); if (m_WorkingSlice.IsNotNull() && m_FillFeedbackContour && positionEvent) { // Project contour into working slice ContourModel *feedbackContour(FeedbackContourTool::GetFeedbackContour()); ContourModel::Pointer projectedContour; // this is not a beautiful solution, just one that works, check T22412 for details int timestep = positionEvent->GetSender()->GetTimeStep(); if (0 != timestep) { int size = feedbackContour->GetNumberOfVertices(timestep); auto feedbackContourTime = mitk::ContourModel::New(); feedbackContourTime->Expand(timestep + 1); for (int loop = 0; loop < size; ++loop) { feedbackContourTime->AddVertex(feedbackContour->GetVertexAt(loop, timestep), 0); } projectedContour = FeedbackContourTool::ProjectContourTo2DSlice(m_WorkingSlice, feedbackContourTime, false, false); } else { projectedContour = FeedbackContourTool::ProjectContourTo2DSlice(m_WorkingSlice, feedbackContour, false, false); } // If there is a projected contour, fill it if (projectedContour.IsNotNull()) { // Get working data to pass to following method so we don't overwrite locked labels in a LabelSetImage mitk::DataNode *workingNode(m_ToolManager->GetWorkingData(0)); mitk::LabelSetImage *labelImage = workingNode != nullptr ? dynamic_cast(workingNode->GetData()) : nullptr; MITK_DEBUG << "Filling Segmentation"; if (labelImage != nullptr) { // m_PaintingPixelValue only decides whether to paint or not // For LabelSetImages we want to paint with the active label value auto activeLabel = labelImage->GetActiveLabel(labelImage->GetActiveLayer())->GetValue(); mitk::ContourModelUtils::FillContourInSlice(projectedContour, 0, m_WorkingSlice, labelImage, m_PaintingPixelValue * activeLabel); } else { mitk::ContourModelUtils::FillContourInSlice(projectedContour, 0, m_WorkingSlice, m_WorkingSlice, m_PaintingPixelValue); } this->WriteBackSegmentationResult(positionEvent, m_WorkingSlice); FeedbackContourTool::SetFeedbackContourVisible(false); } m_ScreenYDifference = 0; m_ScreenXDifference = 0; } } diff --git a/Modules/SurfaceInterpolation/mitkCreateDistanceImageFromSurfaceFilter.cpp b/Modules/SurfaceInterpolation/mitkCreateDistanceImageFromSurfaceFilter.cpp index aebf5e8798..b7bb5d067c 100644 --- a/Modules/SurfaceInterpolation/mitkCreateDistanceImageFromSurfaceFilter.cpp +++ b/Modules/SurfaceInterpolation/mitkCreateDistanceImageFromSurfaceFilter.cpp @@ -1,620 +1,621 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkCreateDistanceImageFromSurfaceFilter.h" #include "mitkImageCast.h" #include "vtkCellArray.h" #include "vtkCellData.h" #include "vtkDoubleArray.h" #include "vtkPolyData.h" #include "vtkSmartPointer.h" #include "itkImageRegionIteratorWithIndex.h" #include "itkNeighborhoodIterator.h" #include void mitk::CreateDistanceImageFromSurfaceFilter::CreateEmptyDistanceImage() { // Determine the bounds of the input points in index- and world-coordinates DistanceImageType::PointType minPointInWorldCoordinates, maxPointInWorldCoordinates; DistanceImageType::IndexType minPointInIndexCoordinates, maxPointInIndexCoordinates; DetermineBounds( minPointInWorldCoordinates, maxPointInWorldCoordinates, minPointInIndexCoordinates, maxPointInIndexCoordinates); // Calculate the extent of the region that contains all given points in MM. // To do this, we take the difference between the maximal and minimal // index-coordinates (must not be less than 1) and multiply it with the // spacing of the reference-image. Vector3D extentMM; for (unsigned int dim = 0; dim < 3; ++dim) { extentMM[dim] = (std::abs(maxPointInIndexCoordinates[dim] - minPointInIndexCoordinates[dim])) * m_ReferenceImage->GetSpacing()[dim]; } /* * Now create an empty distance image. The created image will always have the same number of pixels, independent from * the original image (e.g. always consists of 500000 pixels) and will have an isotropic spacing. * The spacing is calculated like the following: * The image's volume = 500000 Pixels = extentX*spacing*extentY*spacing*extentZ*spacing * So the spacing is: spacing = ( extentX*extentY*extentZ / 500000 )^(1/3) */ double basis = (extentMM[0] * extentMM[1] * extentMM[2]) / m_DistanceImageVolume; double exponent = 1.0 / 3.0; m_DistanceImageSpacing = pow(basis, exponent); // calculate the number of pixels of the distance image for each direction unsigned int numberOfXPixel = extentMM[0] / m_DistanceImageSpacing; unsigned int numberOfYPixel = extentMM[1] / m_DistanceImageSpacing; unsigned int numberOfZPixel = extentMM[2] / m_DistanceImageSpacing; // We increase the sizeOfRegion by 4 as we decrease the origin by 2 later. // This expansion of the region is necessary to achieve a complete // interpolation. DistanceImageType::SizeType sizeOfRegion; sizeOfRegion[0] = numberOfXPixel + 8; sizeOfRegion[1] = numberOfYPixel + 8; sizeOfRegion[2] = numberOfZPixel + 8; // The region starts at index 0,0,0 DistanceImageType::IndexType initialOriginAsIndex; initialOriginAsIndex.Fill(0); DistanceImageType::PointType originAsWorld = minPointInWorldCoordinates; DistanceImageType::RegionType lpRegion; lpRegion.SetSize(sizeOfRegion); lpRegion.SetIndex(initialOriginAsIndex); // We initialize the itk::Image with // * origin and direction to have it correctly placed and rotated in the world // * the largest possible region to set the extent to be calculated // * the isotropic spacing that we have calculated above m_DistanceImageITK = DistanceImageType::New(); m_DistanceImageITK->SetOrigin(originAsWorld); m_DistanceImageITK->SetDirection(m_ReferenceImage->GetDirection()); m_DistanceImageITK->SetRegions(lpRegion); m_DistanceImageITK->SetSpacing(m_DistanceImageSpacing); m_DistanceImageITK->Allocate(); // First of all the image is initialized with the value 10*m_DistanceImageSpacing for each pixel m_DistanceImageDefaultBufferValue = 10 * m_DistanceImageSpacing; m_DistanceImageITK->FillBuffer(m_DistanceImageDefaultBufferValue); // Now we move the origin of the distanceImage 2 index-Coordinates // in all directions DistanceImageType::IndexType originAsIndex; m_DistanceImageITK->TransformPhysicalPointToIndex(originAsWorld, originAsIndex); originAsIndex[0] -= 2; originAsIndex[1] -= 2; originAsIndex[2] -= 2; m_DistanceImageITK->TransformIndexToPhysicalPoint(originAsIndex, originAsWorld); m_DistanceImageITK->SetOrigin(originAsWorld); } mitk::CreateDistanceImageFromSurfaceFilter::CreateDistanceImageFromSurfaceFilter() + : m_DistanceImageSpacing(0.0), m_DistanceImageDefaultBufferValue(0.0) { m_DistanceImageVolume = 50000; this->m_UseProgressBar = false; this->m_ProgressStepSize = 5; mitk::Image::Pointer output = mitk::Image::New(); this->SetNthOutput(0, output.GetPointer()); } mitk::CreateDistanceImageFromSurfaceFilter::~CreateDistanceImageFromSurfaceFilter() { } void mitk::CreateDistanceImageFromSurfaceFilter::GenerateData() { this->PreprocessContourPoints(); this->CreateEmptyDistanceImage(); // First of all we have to build the equation-system from the existing contour-edge-points this->CreateSolutionMatrixAndFunctionValues(); if (this->m_UseProgressBar) mitk::ProgressBar::GetInstance()->Progress(1); m_Weights = m_SolutionMatrix.partialPivLu().solve(m_FunctionValues); if (this->m_UseProgressBar) mitk::ProgressBar::GetInstance()->Progress(2); // The last step is to create the distance map with the interpolated distance function this->FillDistanceImage(); if (this->m_UseProgressBar) mitk::ProgressBar::GetInstance()->Progress(2); m_Centers.clear(); m_Normals.clear(); } void mitk::CreateDistanceImageFromSurfaceFilter::PreprocessContourPoints() { unsigned int numberOfInputs = this->GetNumberOfIndexedInputs(); if (numberOfInputs == 0) { MITK_ERROR << "mitk::CreateDistanceImageFromSurfaceFilter: No input available. Please set an input!" << std::endl; itkExceptionMacro("mitk::CreateDistanceImageFromSurfaceFilter: No input available. Please set an input!"); return; } // First of all we have to extract the nomals and the surface points. // Duplicated points can be eliminated vtkSmartPointer polyData; vtkSmartPointer currentCellNormals; vtkSmartPointer existingPolys; vtkSmartPointer existingPoints; double p[3]; PointType currentPoint; PointType normal; for (unsigned int i = 0; i < numberOfInputs; i++) { auto currentSurface = this->GetInput(i); polyData = currentSurface->GetVtkPolyData(); if (polyData->GetNumberOfPolys() == 0) { MITK_INFO << "mitk::CreateDistanceImageFromSurfaceFilter: No input-polygons available. Please be sure the input " "surface consists of polygons!" << std::endl; } currentCellNormals = vtkDoubleArray::SafeDownCast(polyData->GetCellData()->GetNormals()); existingPolys = polyData->GetPolys(); existingPoints = polyData->GetPoints(); existingPolys->InitTraversal(); vtkIdType *cell(nullptr); vtkIdType cellSize(0); for (existingPolys->InitTraversal(); existingPolys->GetNextCell(cellSize, cell);) { for (vtkIdType j = 0; j < cellSize; j++) { existingPoints->GetPoint(cell[j], p); currentPoint.copy_in(p); int count = std::count(m_Centers.begin(), m_Centers.end(), currentPoint); if (count == 0) { double currentNormal[3]; currentCellNormals->GetTuple(cell[j], currentNormal); normal.copy_in(currentNormal); m_Normals.push_back(normal); m_Centers.push_back(currentPoint); } } // end for all points } // end for all cells } // end for all outputs } void mitk::CreateDistanceImageFromSurfaceFilter::CreateSolutionMatrixAndFunctionValues() { // For we can now calculate the exact size of the centers we initialize the data structures unsigned int numberOfCenters = m_Centers.size(); m_Centers.reserve(numberOfCenters * 3); m_FunctionValues.resize(numberOfCenters * 3); m_FunctionValues.fill(0); PointType currentPoint; PointType normal; // Create inner points for (unsigned int i = 0; i < numberOfCenters; i++) { currentPoint = m_Centers.at(i); normal = m_Normals.at(i); currentPoint[0] = currentPoint[0] - normal[0] * m_DistanceImageSpacing; currentPoint[1] = currentPoint[1] - normal[1] * m_DistanceImageSpacing; currentPoint[2] = currentPoint[2] - normal[2] * m_DistanceImageSpacing; m_Centers.push_back(currentPoint); m_FunctionValues[numberOfCenters + i] = -m_DistanceImageSpacing; } // Create outer points for (unsigned int i = 0; i < numberOfCenters; i++) { currentPoint = m_Centers.at(i); normal = m_Normals.at(i); currentPoint[0] = currentPoint[0] + normal[0] * m_DistanceImageSpacing; currentPoint[1] = currentPoint[1] + normal[1] * m_DistanceImageSpacing; currentPoint[2] = currentPoint[2] + normal[2] * m_DistanceImageSpacing; m_Centers.push_back(currentPoint); m_FunctionValues[numberOfCenters * 2 + i] = m_DistanceImageSpacing; } // Now we have created all centers and all function values. Next step is to create the solution matrix numberOfCenters = m_Centers.size(); m_SolutionMatrix.resize(numberOfCenters, numberOfCenters); m_Weights.resize(numberOfCenters); PointType p1; PointType p2; double norm; for (unsigned int i = 0; i < numberOfCenters; i++) { for (unsigned int j = 0; j < numberOfCenters; j++) { // Calculate the RBF value. Currently using Phi(r) = r with r is the euclidian distance between two points p1 = m_Centers.at(i); p2 = m_Centers.at(j); p1 = p1 - p2; norm = p1.two_norm(); m_SolutionMatrix(i, j) = norm; } } } void mitk::CreateDistanceImageFromSurfaceFilter::FillDistanceImage() { /* * Now we must calculate the distance for each pixel. But instead of calculating the distance value * for all of the image's pixels we proceed similar to the region growing algorithm: * * 1. Take the first pixel from the narrowband_point_list and calculate the distance for each neighbor (6er) * 2. If the current index's distance value is below a certain threshold push it into the list * 3. Next iteration take the next index from the list and originAsIndex with 1. again * * This is done until the narrowband_point_list is empty. */ typedef itk::ImageRegionIteratorWithIndex ImageIterator; typedef itk::NeighborhoodIterator NeighborhoodImageIterator; std::queue narrowbandPoints; PointType currentPoint = m_Centers.at(0); double distance = this->CalculateDistanceValue(currentPoint); // create itk::Point from vnl_vector DistanceImageType::PointType currentPointAsPoint; currentPointAsPoint[0] = currentPoint[0]; currentPointAsPoint[1] = currentPoint[1]; currentPointAsPoint[2] = currentPoint[2]; // Transform the input point in world-coordinates to index-coordinates DistanceImageType::IndexType currentIndex; m_DistanceImageITK->TransformPhysicalPointToIndex(currentPointAsPoint, currentIndex); assert( m_DistanceImageITK->GetLargestPossibleRegion().IsInside(currentIndex)); // we are quite certain this should hold narrowbandPoints.push(currentIndex); m_DistanceImageITK->SetPixel(currentIndex, distance); NeighborhoodImageIterator::RadiusType radius; radius.Fill(1); NeighborhoodImageIterator nIt(radius, m_DistanceImageITK, m_DistanceImageITK->GetLargestPossibleRegion()); unsigned int relativeNbIdx[] = {4, 10, 12, 14, 16, 22}; bool isInBounds = false; while (!narrowbandPoints.empty()) { nIt.SetLocation(narrowbandPoints.front()); narrowbandPoints.pop(); unsigned int *relativeNb = &relativeNbIdx[0]; for (int i = 0; i < 6; i++) { nIt.GetPixel(*relativeNb, isInBounds); if (isInBounds && nIt.GetPixel(*relativeNb) == m_DistanceImageDefaultBufferValue) { currentIndex = nIt.GetIndex(*relativeNb); // Transform the currently checked point from index-coordinates to // world-coordinates m_DistanceImageITK->TransformIndexToPhysicalPoint(currentIndex, currentPointAsPoint); // create a vnl_vector currentPoint[0] = currentPointAsPoint[0]; currentPoint[1] = currentPointAsPoint[1]; currentPoint[2] = currentPointAsPoint[2]; // and check the distance distance = this->CalculateDistanceValue(currentPoint); if (std::fabs(distance) <= m_DistanceImageSpacing * 2) { nIt.SetPixel(*relativeNb, distance); narrowbandPoints.push(currentIndex); } } relativeNb++; } } ImageIterator imgRegionIterator(m_DistanceImageITK, m_DistanceImageITK->GetLargestPossibleRegion()); imgRegionIterator.GoToBegin(); double prevPixelVal = 1; DistanceImageType::IndexType _size; _size.Fill(-1); _size += m_DistanceImageITK->GetLargestPossibleRegion().GetSize(); // Set every pixel inside the surface to -m_DistanceImageDefaultBufferValue except the edge point (so that the // received surface is closed) while (!imgRegionIterator.IsAtEnd()) { if (imgRegionIterator.Get() == m_DistanceImageDefaultBufferValue && prevPixelVal < 0) { while (imgRegionIterator.Get() == m_DistanceImageDefaultBufferValue) { if (imgRegionIterator.GetIndex()[0] == _size[0] || imgRegionIterator.GetIndex()[1] == _size[1] || imgRegionIterator.GetIndex()[2] == _size[2] || imgRegionIterator.GetIndex()[0] == 0U || imgRegionIterator.GetIndex()[1] == 0U || imgRegionIterator.GetIndex()[2] == 0U) { imgRegionIterator.Set(m_DistanceImageDefaultBufferValue); prevPixelVal = m_DistanceImageDefaultBufferValue; ++imgRegionIterator; break; } else { imgRegionIterator.Set((-1) * m_DistanceImageDefaultBufferValue); ++imgRegionIterator; prevPixelVal = (-1) * m_DistanceImageDefaultBufferValue; } } } else if (imgRegionIterator.GetIndex()[0] == _size[0] || imgRegionIterator.GetIndex()[1] == _size[1] || imgRegionIterator.GetIndex()[2] == _size[2] || imgRegionIterator.GetIndex()[0] == 0U || imgRegionIterator.GetIndex()[1] == 0U || imgRegionIterator.GetIndex()[2] == 0U) { imgRegionIterator.Set(m_DistanceImageDefaultBufferValue); prevPixelVal = m_DistanceImageDefaultBufferValue; ++imgRegionIterator; } else { prevPixelVal = imgRegionIterator.Get(); ++imgRegionIterator; } } Image::Pointer resultImage = this->GetOutput(); // Cast the created distance-Image from itk::Image to the mitk::Image // that is our output. CastToMitkImage(m_DistanceImageITK, resultImage); } double mitk::CreateDistanceImageFromSurfaceFilter::CalculateDistanceValue(PointType p) { double distanceValue(0); PointType p1; PointType p2; double norm; CenterList::iterator centerIter; unsigned int count(0); for (centerIter = m_Centers.begin(); centerIter != m_Centers.end(); centerIter++) { p1 = *centerIter; p2 = p - p1; norm = p2.two_norm(); distanceValue = distanceValue + (norm * m_Weights[count]); ++count; } return distanceValue; } void mitk::CreateDistanceImageFromSurfaceFilter::GenerateOutputInformation() { } void mitk::CreateDistanceImageFromSurfaceFilter::PrintEquationSystem() { std::stringstream out; out << "Nummber of rows: " << m_SolutionMatrix.rows() << " ****** Number of columns: " << m_SolutionMatrix.cols() << endl; out << "[ "; for (int i = 0; i < m_SolutionMatrix.rows(); i++) { for (int j = 0; j < m_SolutionMatrix.cols(); j++) { out << m_SolutionMatrix(i, j) << " "; } out << ";" << endl; } out << " ]\n\n\n"; for (unsigned int i = 0; i < m_Centers.size(); i++) { out << m_Centers.at(i) << ";" << endl; } std::cout << "Equation system: \n\n\n" << out.str(); } void mitk::CreateDistanceImageFromSurfaceFilter::SetInput(const mitk::Surface *surface) { this->SetInput(0, surface); } void mitk::CreateDistanceImageFromSurfaceFilter::SetInput(unsigned int idx, const mitk::Surface *surface) { if (this->GetInput(idx) != surface) { this->SetNthInput(idx, const_cast(surface)); this->Modified(); } } const mitk::Surface *mitk::CreateDistanceImageFromSurfaceFilter::GetInput() { if (this->GetNumberOfIndexedInputs() < 1) return nullptr; return static_cast(this->ProcessObject::GetInput(0)); } const mitk::Surface *mitk::CreateDistanceImageFromSurfaceFilter::GetInput(unsigned int idx) { if (this->GetNumberOfIndexedInputs() < 1) return nullptr; return static_cast(this->ProcessObject::GetInput(idx)); } void mitk::CreateDistanceImageFromSurfaceFilter::RemoveInputs(mitk::Surface *input) { DataObjectPointerArraySizeType nb = this->GetNumberOfIndexedInputs(); for (DataObjectPointerArraySizeType i = 0; i < nb; i++) { if (this->GetInput(i) == input) { this->RemoveInput(i); return; } } } void mitk::CreateDistanceImageFromSurfaceFilter::Reset() { for (unsigned int i = 0; i < this->GetNumberOfIndexedInputs(); i++) { this->PopBackInput(); } this->SetNumberOfIndexedInputs(0); this->SetNumberOfIndexedOutputs(1); mitk::Image::Pointer output = mitk::Image::New(); this->SetNthOutput(0, output.GetPointer()); } void mitk::CreateDistanceImageFromSurfaceFilter::SetUseProgressBar(bool status) { this->m_UseProgressBar = status; } void mitk::CreateDistanceImageFromSurfaceFilter::SetProgressStepSize(unsigned int stepSize) { this->m_ProgressStepSize = stepSize; } void mitk::CreateDistanceImageFromSurfaceFilter::SetReferenceImage(itk::ImageBase<3>::Pointer referenceImage) { m_ReferenceImage = referenceImage; } void mitk::CreateDistanceImageFromSurfaceFilter::DetermineBounds( DistanceImageType::PointType &minPointInWorldCoordinates, DistanceImageType::PointType &maxPointInWorldCoordinates, DistanceImageType::IndexType &minPointInIndexCoordinates, DistanceImageType::IndexType &maxPointInIndexCoordinates) { PointType firstCenter = m_Centers.at(0); DistanceImageType::PointType tmpPoint; tmpPoint[0] = firstCenter[0]; tmpPoint[1] = firstCenter[1]; tmpPoint[2] = firstCenter[2]; // transform the first point from world-coordinates to index-coordinates itk::ContinuousIndex tmpIndex; m_ReferenceImage->TransformPhysicalPointToContinuousIndex(tmpPoint, tmpIndex); // initialize the variables with this first point DistanceImageType::IndexValueType xmin = tmpIndex[0]; DistanceImageType::IndexValueType ymin = tmpIndex[1]; DistanceImageType::IndexValueType zmin = tmpIndex[2]; DistanceImageType::IndexValueType xmax = tmpIndex[0]; DistanceImageType::IndexValueType ymax = tmpIndex[1]; DistanceImageType::IndexValueType zmax = tmpIndex[2]; // iterate over the rest of the points auto centerIter = m_Centers.begin(); for (++centerIter; centerIter != m_Centers.end(); centerIter++) { tmpPoint[0] = (*centerIter)[0]; tmpPoint[1] = (*centerIter)[1]; tmpPoint[2] = (*centerIter)[2]; // transform each point from world-coordinates to index-coordinates m_ReferenceImage->TransformPhysicalPointToContinuousIndex(tmpPoint, tmpIndex); // and set the variables accordingly to find the minimum // and maximum in all directions in index-coordinates if (xmin > tmpIndex[0]) { xmin = tmpIndex[0]; } if (ymin > tmpIndex[1]) { ymin = tmpIndex[1]; } if (zmin > tmpIndex[2]) { zmin = tmpIndex[2]; } if (xmax < tmpIndex[0]) { xmax = tmpIndex[0]; } if (ymax < tmpIndex[1]) { ymax = tmpIndex[1]; } if (zmax < tmpIndex[2]) { zmax = tmpIndex[2]; } } // put the found coordinates into Index-Points minPointInIndexCoordinates[0] = xmin; minPointInIndexCoordinates[1] = ymin; minPointInIndexCoordinates[2] = zmin; maxPointInIndexCoordinates[0] = xmax; maxPointInIndexCoordinates[1] = ymax; maxPointInIndexCoordinates[2] = zmax; // and transform them into world-coordinates m_ReferenceImage->TransformIndexToPhysicalPoint(minPointInIndexCoordinates, minPointInWorldCoordinates); m_ReferenceImage->TransformIndexToPhysicalPoint(maxPointInIndexCoordinates, maxPointInWorldCoordinates); } diff --git a/Modules/SurfaceInterpolation/mitkImageToPointCloudFilter.cpp b/Modules/SurfaceInterpolation/mitkImageToPointCloudFilter.cpp index 766e08fb81..8497ab49a5 100644 --- a/Modules/SurfaceInterpolation/mitkImageToPointCloudFilter.cpp +++ b/Modules/SurfaceInterpolation/mitkImageToPointCloudFilter.cpp @@ -1,166 +1,166 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkImageToPointCloudFilter.h" #include #include #include #include #include #include #include #include #include -mitk::ImageToPointCloudFilter::ImageToPointCloudFilter() +mitk::ImageToPointCloudFilter::ImageToPointCloudFilter() : m_Geometry(nullptr) { m_Method = DetectionMethod(0); this->SetNumberOfRequiredInputs(1); this->SetNumberOfIndexedOutputs(1); } mitk::ImageToPointCloudFilter::~ImageToPointCloudFilter() { } void mitk::ImageToPointCloudFilter::GenerateData() { mitk::Image::ConstPointer image = ImageToUnstructuredGridFilter::GetInput(); m_Geometry = image->GetGeometry(); if (image.IsNull()) { MITK_ERROR << "mitk::ImageToContourFilter: No input available. " "Please set the input!" << std::endl; return; } mitk::Image::Pointer notConstImage = const_cast(image.GetPointer()); switch (m_Method) { case 0: AccessByItk_1(notConstImage.GetPointer(), StdDeviations, 2) break; case 1: AccessByItk_1(notConstImage.GetPointer(), StdDeviations, 3) break; case 2: AccessByItk_1(notConstImage.GetPointer(), StdDeviations, 4) break; default: AccessByItk_1(notConstImage.GetPointer(), StdDeviations, 2) break; } } template void mitk::ImageToPointCloudFilter::StdDeviations(itk::Image *image, int amount) { typedef itk::Image InputImageType; typedef itk::CastImageFilter ImagePTypeToFloatPTypeCasterType; typedef itk::LaplacianImageFilter LaplacianFilterType; typename LaplacianFilterType::Pointer lapFilter = LaplacianFilterType::New(); typename ImagePTypeToFloatPTypeCasterType::Pointer caster = ImagePTypeToFloatPTypeCasterType::New(); caster->SetInput(image); caster->Update(); FloatImageType::Pointer fImage = caster->GetOutput(); lapFilter->SetInput(fImage); lapFilter->UpdateLargestPossibleRegion(); mitk::Image::Pointer edgeImage = mitk::ImportItkImage(lapFilter->GetOutput()); mitk::ImageStatisticsCalculator::Pointer statCalc = mitk::ImageStatisticsCalculator::New(); statCalc->SetInputImage(edgeImage); mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer stats = statCalc->GetStatistics(); double mean = stats->GetMean(); double stdDev = stats->GetStd(); double upperThreshold = mean + stdDev * amount; double lowerThreshold = mean - stdDev * amount; typename itk::ImageRegionIterator it(lapFilter->GetOutput(), lapFilter->GetOutput()->GetRequestedRegion()); vtkSmartPointer points = vtkSmartPointer::New(); double greatX = 0, greatY = 0, greatZ = 0; it.GoToBegin(); while (!it.IsAtEnd()) { if (it.Get() > lowerThreshold && it.Get() < upperThreshold) { it.Set(0); } else { it.Set(1); mitk::Point3D imagePoint; mitk::Point3D worldPoint; imagePoint[0] = it.GetIndex()[0]; imagePoint[1] = it.GetIndex()[1]; imagePoint[2] = it.GetIndex()[2]; m_Geometry->IndexToWorld(imagePoint, worldPoint); if (worldPoint[0] > greatX) greatX = worldPoint[0]; if (worldPoint[1] > greatY) greatY = worldPoint[1]; if (worldPoint[2] > greatZ) greatZ = worldPoint[2]; points->InsertNextPoint(worldPoint[0], worldPoint[1], worldPoint[2]); m_NumberOfExtractedPoints++; } ++it; } /*need to build the UnstructuredGrid with at least one vertex otherwise its not visible*/ vtkSmartPointer verts = vtkSmartPointer::New(); verts->GetPointIds()->SetNumberOfIds(m_NumberOfExtractedPoints); for (int i = 0; i < m_NumberOfExtractedPoints; i++) { verts->GetPointIds()->SetId(i, i); } vtkSmartPointer uGrid = vtkSmartPointer::New(); uGrid->Allocate(1); uGrid->InsertNextCell(verts->GetCellType(), verts->GetPointIds()); uGrid->SetPoints(points); mitk::UnstructuredGrid::Pointer outputGrid = mitk::UnstructuredGrid::New(); outputGrid->SetVtkUnstructuredGrid(uGrid); this->SetNthOutput(0, outputGrid); } void mitk::ImageToPointCloudFilter::GenerateOutputInformation() { Superclass::GenerateOutputInformation(); }