diff --git a/Modules/Segmentation/Algorithms/mitkShowSegmentationAsSmoothedSurface.cpp b/Modules/Segmentation/Algorithms/mitkShowSegmentationAsSmoothedSurface.cpp index adf280419f..28774b66b6 100644 --- a/Modules/Segmentation/Algorithms/mitkShowSegmentationAsSmoothedSurface.cpp +++ b/Modules/Segmentation/Algorithms/mitkShowSegmentationAsSmoothedSurface.cpp @@ -1,513 +1,514 @@ /*=================================================================== 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 "mitkShowSegmentationAsSmoothedSurface.h" #include "mitkImageToItk.h" #include "mitkImageCast.h" #include "itkIntelligentBinaryClosingFilter.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include using namespace mitk; using namespace std; ShowSegmentationAsSmoothedSurface::ShowSegmentationAsSmoothedSurface() { } ShowSegmentationAsSmoothedSurface::~ShowSegmentationAsSmoothedSurface() { } void ShowSegmentationAsSmoothedSurface::Initialize(const NonBlockingAlgorithm *other) { Superclass::Initialize(other); bool syncVisibility = false; if (other != nullptr) other->GetParameter("Sync visibility", syncVisibility); SetParameter("Sync visibility", syncVisibility); SetParameter("Wireframe", false); // The Smoothing value is used as variance for a Gauss filter. // A reasonable default value equals the image spacing in mm. SetParameter("Smoothing", 1.0f); // Valid range for decimation value is [0, 1). High values // increase decimation, especially when very close to 1. // A value of 0 disables decimation. SetParameter("Decimation", 0.5f); // Valid range for closing value is [0, 1]. Higher values // increase closing. A value of 0 disables closing. SetParameter("Closing", 0.0f); } bool ShowSegmentationAsSmoothedSurface::ReadyToRun() { try { mitk::Image::Pointer image; GetPointerParameter("Input", image); return image.IsNotNull() && GetGroupNode(); } catch (const invalid_argument &) { return false; } } bool ShowSegmentationAsSmoothedSurface::ThreadedUpdateFunction() { Image::Pointer image; GetPointerParameter("Input", image); float smoothing; GetParameter("Smoothing", smoothing); float decimation; GetParameter("Decimation", decimation); float closing; GetParameter("Closing", closing); int timeNr = 0; GetParameter("TimeNr", timeNr); if (image->GetDimension() == 4) MITK_INFO << "CREATING SMOOTHED POLYGON MODEL (t = " << timeNr << ')'; else MITK_INFO << "CREATING SMOOTHED POLYGON MODEL"; MITK_INFO << " Smoothing = " << smoothing; MITK_INFO << " Decimation = " << decimation; MITK_INFO << " Closing = " << closing; Geometry3D::Pointer geometry = dynamic_cast(image->GetGeometry()->Clone().GetPointer()); // Make ITK image out of MITK image typedef itk::Image CharImageType; typedef itk::Image ShortImageType; typedef itk::Image FloatImageType; if (image->GetDimension() == 4) { ImageTimeSelector::Pointer imageTimeSelector = ImageTimeSelector::New(); imageTimeSelector->SetInput(image); imageTimeSelector->SetTimeNr(timeNr); imageTimeSelector->UpdateLargestPossibleRegion(); image = imageTimeSelector->GetOutput(0); } ImageToItk::Pointer imageToItkFilter = ImageToItk::New(); try { imageToItkFilter->SetInput(image); } catch (const itk::ExceptionObject &e) { // Most probably the input image type is wrong. Binary images are expected to be // >unsigned< char images. MITK_ERROR << e.GetDescription() << endl; return false; } imageToItkFilter->Update(); CharImageType::Pointer itkImage = imageToItkFilter->GetOutput(); // Get bounding box and relabel MITK_INFO << "Extracting VOI..."; int imageLabel = 1; bool roiFound = false; CharImageType::IndexType minIndex; minIndex.Fill(numeric_limits::max()); CharImageType::IndexType maxIndex; maxIndex.Fill(numeric_limits::min()); itk::ImageRegionIteratorWithIndex iter(itkImage, itkImage->GetLargestPossibleRegion()); for (iter.GoToBegin(); !iter.IsAtEnd(); ++iter) { if (iter.Get() == imageLabel) { roiFound = true; iter.Set(1); CharImageType::IndexType currentIndex = iter.GetIndex(); for (unsigned int dim = 0; dim < 3; ++dim) { minIndex[dim] = min(currentIndex[dim], minIndex[dim]); maxIndex[dim] = max(currentIndex[dim], maxIndex[dim]); } } else { iter.Set(0); } } if (!roiFound) { ProgressBar::GetInstance()->Progress(8); MITK_ERROR << "Didn't found segmentation labeled with " << imageLabel << "!" << endl; return false; } ProgressBar::GetInstance()->Progress(1); // Extract and pad bounding box typedef itk::RegionOfInterestImageFilter ROIFilterType; ROIFilterType::Pointer roiFilter = ROIFilterType::New(); CharImageType::RegionType region; CharImageType::SizeType size; for (unsigned int dim = 0; dim < 3; ++dim) { size[dim] = maxIndex[dim] - minIndex[dim] + 1; } region.SetIndex(minIndex); region.SetSize(size); roiFilter->SetInput(itkImage); roiFilter->SetRegionOfInterest(region); roiFilter->ReleaseDataFlagOn(); roiFilter->ReleaseDataBeforeUpdateFlagOn(); typedef itk::ConstantPadImageFilter PadFilterType; PadFilterType::Pointer padFilter = PadFilterType::New(); const PadFilterType::SizeValueType pad[3] = { 10, 10, 10 }; padFilter->SetInput(roiFilter->GetOutput()); padFilter->SetConstant(0); padFilter->SetPadLowerBound(pad); padFilter->SetPadUpperBound(pad); padFilter->ReleaseDataFlagOn(); padFilter->ReleaseDataBeforeUpdateFlagOn(); padFilter->Update(); CharImageType::Pointer roiImage = padFilter->GetOutput(); roiImage->DisconnectPipeline(); roiFilter = nullptr; padFilter = nullptr; // Correct origin of real geometry (changed by cropping and padding) typedef Geometry3D::TransformType TransformType; TransformType::Pointer transform = TransformType::New(); TransformType::OutputVectorType translation; for (unsigned int dim = 0; dim < 3; ++dim) translation[dim] = (int)minIndex[dim] - (int)pad[dim]; transform->SetIdentity(); transform->Translate(translation); geometry->Compose(transform, true); ProgressBar::GetInstance()->Progress(1); // Median MITK_INFO << "Median..."; typedef itk::BinaryMedianImageFilter MedianFilterType; MedianFilterType::Pointer medianFilter = MedianFilterType::New(); CharImageType::SizeType radius = { 0 }; medianFilter->SetRadius(radius); medianFilter->SetBackgroundValue(0); medianFilter->SetForegroundValue(1); medianFilter->SetInput(roiImage); medianFilter->ReleaseDataFlagOn(); medianFilter->ReleaseDataBeforeUpdateFlagOn(); medianFilter->Update(); ProgressBar::GetInstance()->Progress(1); // Intelligent closing MITK_INFO << "Intelligent closing..."; unsigned int surfaceRatio = (unsigned int)((1.0f - closing) * 100.0f); typedef itk::IntelligentBinaryClosingFilter ClosingFilterType; ClosingFilterType::Pointer closingFilter = ClosingFilterType::New(); closingFilter->SetInput(medianFilter->GetOutput()); closingFilter->ReleaseDataFlagOn(); closingFilter->ReleaseDataBeforeUpdateFlagOn(); closingFilter->SetSurfaceRatio(surfaceRatio); closingFilter->Update(); ShortImageType::Pointer closedImage = closingFilter->GetOutput(); closedImage->DisconnectPipeline(); roiImage = nullptr; medianFilter = nullptr; closingFilter = nullptr; ProgressBar::GetInstance()->Progress(1); // Gaussian blur MITK_INFO << "Gauss..."; typedef itk::BinaryThresholdImageFilter BinaryThresholdToFloatFilterType; BinaryThresholdToFloatFilterType::Pointer binThresToFloatFilter = BinaryThresholdToFloatFilterType::New(); binThresToFloatFilter->SetInput(closedImage); binThresToFloatFilter->SetLowerThreshold(1); binThresToFloatFilter->SetUpperThreshold(1); binThresToFloatFilter->SetInsideValue(100); binThresToFloatFilter->SetOutsideValue(0); binThresToFloatFilter->ReleaseDataFlagOn(); binThresToFloatFilter->ReleaseDataBeforeUpdateFlagOn(); typedef itk::DiscreteGaussianImageFilter GaussianFilterType; // From the following line on, IntelliSense (VS 2008) is broken. Any idea how to fix it? GaussianFilterType::Pointer gaussFilter = GaussianFilterType::New(); gaussFilter->SetInput(binThresToFloatFilter->GetOutput()); gaussFilter->SetUseImageSpacing(true); gaussFilter->SetVariance(smoothing); gaussFilter->ReleaseDataFlagOn(); gaussFilter->ReleaseDataBeforeUpdateFlagOn(); typedef itk::BinaryThresholdImageFilter BinaryThresholdFromFloatFilterType; BinaryThresholdFromFloatFilterType::Pointer binThresFromFloatFilter = BinaryThresholdFromFloatFilterType::New(); binThresFromFloatFilter->SetInput(gaussFilter->GetOutput()); binThresFromFloatFilter->SetLowerThreshold(50); binThresFromFloatFilter->SetUpperThreshold(255); binThresFromFloatFilter->SetInsideValue(1); binThresFromFloatFilter->SetOutsideValue(0); binThresFromFloatFilter->ReleaseDataFlagOn(); binThresFromFloatFilter->ReleaseDataBeforeUpdateFlagOn(); binThresFromFloatFilter->Update(); CharImageType::Pointer blurredImage = binThresFromFloatFilter->GetOutput(); blurredImage->DisconnectPipeline(); closedImage = nullptr; binThresToFloatFilter = nullptr; gaussFilter = nullptr; ProgressBar::GetInstance()->Progress(1); // Fill holes MITK_INFO << "Filling cavities..."; typedef itk::ConnectedThresholdImageFilter ConnectedThresholdFilterType; ConnectedThresholdFilterType::Pointer connectedThresFilter = ConnectedThresholdFilterType::New(); CharImageType::IndexType corner; corner[0] = 0; corner[1] = 0; corner[2] = 0; connectedThresFilter->SetInput(blurredImage); connectedThresFilter->SetSeed(corner); connectedThresFilter->SetLower(0); connectedThresFilter->SetUpper(0); connectedThresFilter->SetReplaceValue(2); connectedThresFilter->ReleaseDataFlagOn(); connectedThresFilter->ReleaseDataBeforeUpdateFlagOn(); typedef itk::BinaryThresholdImageFilter BinaryThresholdFilterType; BinaryThresholdFilterType::Pointer binThresFilter = BinaryThresholdFilterType::New(); binThresFilter->SetInput(connectedThresFilter->GetOutput()); binThresFilter->SetLowerThreshold(0); binThresFilter->SetUpperThreshold(0); binThresFilter->SetInsideValue(50); binThresFilter->SetOutsideValue(0); binThresFilter->ReleaseDataFlagOn(); binThresFilter->ReleaseDataBeforeUpdateFlagOn(); typedef itk::AddImageFilter AddFilterType; AddFilterType::Pointer addFilter = AddFilterType::New(); addFilter->SetInput1(blurredImage); addFilter->SetInput2(binThresFilter->GetOutput()); addFilter->ReleaseDataFlagOn(); addFilter->ReleaseDataBeforeUpdateFlagOn(); addFilter->Update(); ProgressBar::GetInstance()->Progress(1); // Surface extraction MITK_INFO << "Surface extraction..."; Image::Pointer filteredImage = Image::New(); CastToMitkImage(addFilter->GetOutput(), filteredImage); filteredImage->SetGeometry(geometry); ImageToSurfaceFilter::Pointer imageToSurfaceFilter = ImageToSurfaceFilter::New(); imageToSurfaceFilter->SetInput(filteredImage); imageToSurfaceFilter->SetThreshold(50); imageToSurfaceFilter->SmoothOn(); imageToSurfaceFilter->SetDecimate(ImageToSurfaceFilter::NoDecimation); m_Surface = imageToSurfaceFilter->GetOutput(0); ProgressBar::GetInstance()->Progress(1); // Mesh decimation if (decimation > 0.0f && decimation < 1.0f) { MITK_INFO << "Quadric mesh decimation..."; vtkQuadricDecimation *quadricDecimation = vtkQuadricDecimation::New(); quadricDecimation->SetInputData(m_Surface->GetVtkPolyData()); quadricDecimation->SetTargetReduction(decimation); quadricDecimation->AttributeErrorMetricOn(); quadricDecimation->GlobalWarningDisplayOff(); quadricDecimation->Update(); vtkCleanPolyData* cleaner = vtkCleanPolyData::New(); cleaner->SetInputConnection(quadricDecimation->GetOutputPort()); cleaner->PieceInvariantOn(); cleaner->ConvertLinesToPointsOn(); cleaner->ConvertStripsToPolysOn(); cleaner->PointMergingOn(); cleaner->Update(); m_Surface->SetVtkPolyData(cleaner->GetOutput()); } ProgressBar::GetInstance()->Progress(1); // Compute Normals vtkPolyDataNormals* computeNormals = vtkPolyDataNormals::New(); computeNormals->SetInputData(m_Surface->GetVtkPolyData()); computeNormals->SetFeatureAngle(360.0f); + computeNormals->AutoOrientNormalsOn(); computeNormals->FlipNormalsOff(); computeNormals->Update(); m_Surface->SetVtkPolyData(computeNormals->GetOutput()); return true; } void ShowSegmentationAsSmoothedSurface::ThreadedUpdateSuccessful() { DataNode::Pointer node = DataNode::New(); bool wireframe = false; GetParameter("Wireframe", wireframe); if (wireframe) { VtkRepresentationProperty *representation = dynamic_cast( node->GetProperty("material.representation")); if (representation != nullptr) representation->SetRepresentationToWireframe(); } node->SetProperty("opacity", FloatProperty::New(1.0)); node->SetProperty("line width", IntProperty::New(1)); node->SetProperty("scalar visibility", BoolProperty::New(false)); std::string groupNodeName = "surface"; DataNode *groupNode = GetGroupNode(); if (groupNode != nullptr) groupNode->GetName(groupNodeName); node->SetProperty("name", StringProperty::New(groupNodeName)); node->SetData(m_Surface); BaseProperty *colorProperty = groupNode->GetProperty("color"); if (colorProperty != nullptr) node->ReplaceProperty("color", colorProperty->Clone()); else node->SetProperty("color", ColorProperty::New(1.0f, 0.0f, 0.0f)); bool showResult = true; GetParameter("Show result", showResult); bool syncVisibility = false; GetParameter("Sync visibility", syncVisibility); Image::Pointer image; GetPointerParameter("Input", image); BaseProperty *organTypeProperty = image->GetProperty("organ type"); if (organTypeProperty != nullptr) m_Surface->SetProperty("organ type", organTypeProperty); BaseProperty *visibleProperty = groupNode->GetProperty("visible"); if (visibleProperty != nullptr && syncVisibility) node->ReplaceProperty("visible", visibleProperty->Clone()); else node->SetProperty("visible", BoolProperty::New(showResult)); InsertBelowGroupNode(node); Superclass::ThreadedUpdateSuccessful(); } diff --git a/Modules/Segmentation/Algorithms/mitkShowSegmentationAsSurface.cpp b/Modules/Segmentation/Algorithms/mitkShowSegmentationAsSurface.cpp index 6f645835fc..33352af7f5 100644 --- a/Modules/Segmentation/Algorithms/mitkShowSegmentationAsSurface.cpp +++ b/Modules/Segmentation/Algorithms/mitkShowSegmentationAsSurface.cpp @@ -1,233 +1,235 @@ /*=================================================================== 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 "mitkShowSegmentationAsSurface.h" #include "mitkManualSegmentationToSurfaceFilter.h" #include "mitkVtkRepresentationProperty.h" #include #include namespace mitk { ShowSegmentationAsSurface::ShowSegmentationAsSurface() :m_UIDGeneratorSurfaces("Surface_"), m_AddToTree(false) { } ShowSegmentationAsSurface::~ShowSegmentationAsSurface() { } void ShowSegmentationAsSurface::Initialize(const NonBlockingAlgorithm* other) { Superclass::Initialize(other); bool syncVisibility(false); if (other) { other->GetParameter("Sync visibility", syncVisibility); } SetParameter("Sync visibility", syncVisibility ); SetParameter("Median kernel size", 3u); SetParameter("Apply median", true ); SetParameter("Smooth", true ); SetParameter("Gaussian SD", 1.5f ); SetParameter("Decimate mesh", true ); SetParameter("Decimation rate", 0.8f ); SetParameter("Wireframe", false ); } bool ShowSegmentationAsSurface::ReadyToRun() { try { Image::Pointer image; GetPointerParameter("Input", image); return image.IsNotNull() && GetGroupNode(); } catch (std::invalid_argument&) { return false; } } bool ShowSegmentationAsSurface::ThreadedUpdateFunction() { Image::Pointer image; GetPointerParameter("Input", image); bool smooth(true); GetParameter("Smooth", smooth); bool applyMedian(true); GetParameter("Apply median", applyMedian); bool decimateMesh(true); GetParameter("Decimate mesh", decimateMesh); unsigned int medianKernelSize(3); GetParameter("Median kernel size", medianKernelSize); float gaussianSD(1.5); GetParameter("Gaussian SD", gaussianSD ); float reductionRate(0.8); GetParameter("Decimation rate", reductionRate ); MITK_INFO << "Creating polygon model with smoothing " << smooth << " gaussianSD " << gaussianSD << " median " << applyMedian << " median kernel " << medianKernelSize << " mesh reduction " << decimateMesh << " reductionRate " << reductionRate; ManualSegmentationToSurfaceFilter::Pointer surfaceFilter = ManualSegmentationToSurfaceFilter::New(); surfaceFilter->SetInput( image ); surfaceFilter->SetThreshold( 0.5 ); //expects binary image with zeros and ones surfaceFilter->SetUseGaussianImageSmooth(smooth); // apply gaussian to thresholded image ? surfaceFilter->SetSmooth(smooth); if (smooth) { surfaceFilter->InterpolationOn(); surfaceFilter->SetGaussianStandardDeviation( gaussianSD ); } surfaceFilter->SetMedianFilter3D(applyMedian); // apply median to segmentation before marching cubes ? if (applyMedian) { surfaceFilter->SetMedianKernelSize(medianKernelSize, medianKernelSize, medianKernelSize); // apply median to segmentation before marching cubes } //fix to avoid vtk warnings see bug #5390 if ( image->GetDimension() > 3 ) decimateMesh = false; if (decimateMesh) { surfaceFilter->SetDecimate( ImageToSurfaceFilter::QuadricDecimation ); surfaceFilter->SetTargetReduction( reductionRate ); } else { surfaceFilter->SetDecimate( ImageToSurfaceFilter::NoDecimation ); } surfaceFilter->UpdateLargestPossibleRegion(); // calculate normals for nicer display m_Surface = surfaceFilter->GetOutput(); vtkPolyData* polyData = m_Surface->GetVtkPolyData(); if (!polyData) throw std::logic_error("Could not create polygon model"); polyData->SetVerts(0); polyData->SetLines(0); if ( smooth || applyMedian || decimateMesh) { vtkPolyDataNormals* normalsGen = vtkPolyDataNormals::New(); + normalsGen->AutoOrientNormalsOn(); + normalsGen->FlipNormalsOff(); normalsGen->SetInputData( polyData ); normalsGen->Update(); m_Surface->SetVtkPolyData( normalsGen->GetOutput() ); normalsGen->Delete(); } else { m_Surface->SetVtkPolyData( polyData ); } return true; } void ShowSegmentationAsSurface::ThreadedUpdateSuccessful() { m_Node = DataNode::New(); bool wireframe(false); GetParameter("Wireframe", wireframe ); if (wireframe) { VtkRepresentationProperty *np = dynamic_cast(m_Node->GetProperty("material.representation")); if (np) np->SetRepresentationToWireframe(); } m_Node->SetProperty("opacity", FloatProperty::New(0.3) ); m_Node->SetProperty("line width", IntProperty::New(1) ); m_Node->SetProperty("scalar visibility", BoolProperty::New(false) ); std::string groupNodesName ("surface"); DataNode* groupNode = GetGroupNode(); if (groupNode) { groupNode->GetName( groupNodesName ); //if parameter smooth is set add extension to node name bool smooth(true); GetParameter("Smooth", smooth); if(smooth) groupNodesName.append("_smoothed"); } m_Node->SetProperty( "name", StringProperty::New(groupNodesName) ); // synchronize this object's color with the parent's color //surfaceNode->SetProperty( "color", parentNode->GetProperty( "color" ) ); //surfaceNode->SetProperty( "visible", parentNode->GetProperty( "visible" ) ); m_Node->SetData( m_Surface ); BaseProperty* colorProp = groupNode->GetProperty("color"); if (colorProp) m_Node->ReplaceProperty("color", colorProp->Clone()); else m_Node->SetProperty("color", ColorProperty::New(1.0, 1.0, 0.0)); bool showResult(true); GetParameter("Show result", showResult ); bool syncVisibility(false); GetParameter("Sync visibility", syncVisibility ); Image::Pointer image; GetPointerParameter("Input", image); BaseProperty* organTypeProp = image->GetProperty("organ type"); if (organTypeProp) m_Surface->SetProperty("organ type", organTypeProp); BaseProperty* visibleProp = groupNode->GetProperty("visible"); if (visibleProp && syncVisibility) m_Node->ReplaceProperty("visible", visibleProp->Clone()); else m_Node->SetProperty("visible", BoolProperty::New(showResult)); InsertBelowGroupNode(m_Node); Superclass::ThreadedUpdateSuccessful(); } } // namespace