mitkImageToSurfaceFilter.cpp

Go to the documentation of this file.

/*============================================================================

The Medical Imaging Interaction Toolkit (MITK)

Copyright (c) German Cancer Research Center (DKFZ)
All rights reserved.

Use of this source code is governed by a 3-clause BSD license that can be
found in the LICENSE file.

============================================================================*/

#include "mitkException.h"
#include <mitkImageToSurfaceFilter.h>
#include <vtkDecimatePro.h>
#include <vtkImageChangeInformation.h>
#include <vtkImageData.h>
#include <vtkLinearTransform.h>
#include <vtkMath.h>
#include <vtkMatrix4x4.h>
#include <vtkQuadricDecimation.h>

#include <vtkCleanPolyData.h>
#include <vtkPolyDataNormals.h>
#include <vtkSmartPointer.h>

#include "mitkProgressBar.h"

mitk::ImageToSurfaceFilter::ImageToSurfaceFilter()
  : m_Smooth(false),
    m_Decimate(NoDecimation),
    m_Threshold(1.0),
    m_TargetReduction(0.95f),
    m_SmoothIteration(50),
    m_SmoothRelaxation(0.1)
{
}

mitk::ImageToSurfaceFilter::~ImageToSurfaceFilter()
{
}

void mitk::ImageToSurfaceFilter::CreateSurface(int time,
                                               vtkImageData *vtkimage,
                                               mitk::Surface *surface,
                                               const ScalarType threshold)
{
  vtkImageChangeInformation *indexCoordinatesImageFilter = vtkImageChangeInformation::New();
  indexCoordinatesImageFilter->SetInputData(vtkimage);
  indexCoordinatesImageFilter->SetOutputOrigin(0.0, 0.0, 0.0);

  // MarchingCube -->create Surface
  vtkSmartPointer<vtkMarchingCubes> skinExtractor = vtkSmartPointer<vtkMarchingCubes>::New();
  skinExtractor->ComputeScalarsOff();
  skinExtractor->SetInputConnection(indexCoordinatesImageFilter->GetOutputPort()); // RC++
  indexCoordinatesImageFilter->Delete();
  skinExtractor->SetValue(0, threshold);

  vtkPolyData *polydata;
  skinExtractor->Update();
  polydata = skinExtractor->GetOutput();
  polydata->Register(nullptr); // RC++

  if (m_Smooth)
  {
    vtkSmoothPolyDataFilter *smoother = vtkSmoothPolyDataFilter::New();
    // read poly1 (poly1 can be the original polygon, or the decimated polygon)
    smoother->SetInputConnection(skinExtractor->GetOutputPort()); // RC++
    smoother->SetNumberOfIterations(m_SmoothIteration);
    smoother->SetRelaxationFactor(m_SmoothRelaxation);
    smoother->SetFeatureAngle(60);
    smoother->FeatureEdgeSmoothingOff();
    smoother->BoundarySmoothingOff();
    smoother->SetConvergence(0);
    smoother->Update();

    polydata->Delete(); // RC--
    polydata = smoother->GetOutput();
    polydata->Register(nullptr); // RC++
    smoother->Delete();
  }
  ProgressBar::GetInstance()->Progress();

  // decimate = to reduce number of polygons
  if (m_Decimate == DecimatePro)
  {
    vtkDecimatePro *decimate = vtkDecimatePro::New();
    decimate->SplittingOff();
    decimate->SetErrorIsAbsolute(5);
    decimate->SetFeatureAngle(30);
    decimate->PreserveTopologyOn();
    decimate->BoundaryVertexDeletionOff();
    decimate->SetDegree(10); // std-value is 25!

    decimate->SetInputData(polydata); // RC++
    decimate->SetTargetReduction(m_TargetReduction);
    decimate->SetMaximumError(0.002);
    decimate->Update();

    polydata->Delete(); // RC--
    polydata = decimate->GetOutput();
    polydata->Register(nullptr); // RC++
    decimate->Delete();
  }
  else if (m_Decimate == QuadricDecimation)
  {
    vtkQuadricDecimation *decimate = vtkQuadricDecimation::New();
    decimate->SetTargetReduction(m_TargetReduction);

    decimate->SetInputData(polydata);
    decimate->Update();
    polydata->Delete();
    polydata = decimate->GetOutput();
    polydata->Register(nullptr);
    decimate->Delete();
  }

  ProgressBar::GetInstance()->Progress();

  if (polydata->GetNumberOfPoints() > 0)
  {
    mitk::Vector3D spacing = GetInput()->GetGeometry(time)->GetSpacing();

    vtkPoints *points = polydata->GetPoints();
    vtkMatrix4x4 *vtkmatrix = vtkMatrix4x4::New();
    GetInput()->GetGeometry(time)->GetVtkTransform()->GetMatrix(vtkmatrix);
    double(*matrix)[4] = vtkmatrix->Element;

    unsigned int i, j;
    for (i = 0; i < 3; ++i)
      for (j = 0; j < 3; ++j)
        matrix[i][j] /= spacing[j];

    unsigned int n = points->GetNumberOfPoints();
    double point[3];

    for (i = 0; i < n; i++)
    {
      points->GetPoint(i, point);
      mitkVtkLinearTransformPoint(matrix, point, point);
      points->SetPoint(i, point);
    }
    vtkmatrix->Delete();
  }
  ProgressBar::GetInstance()->Progress();

  // determine point_data normals for the poly data points.
  vtkSmartPointer<vtkPolyDataNormals> normalsGenerator = vtkSmartPointer<vtkPolyDataNormals>::New();
  normalsGenerator->SetInputData(polydata);
  normalsGenerator->FlipNormalsOn();

  vtkSmartPointer<vtkCleanPolyData> cleanPolyDataFilter = vtkSmartPointer<vtkCleanPolyData>::New();
  cleanPolyDataFilter->SetInputConnection(normalsGenerator->GetOutputPort());
  cleanPolyDataFilter->PieceInvariantOff();
  cleanPolyDataFilter->ConvertLinesToPointsOff();
  cleanPolyDataFilter->ConvertPolysToLinesOff();
  cleanPolyDataFilter->ConvertStripsToPolysOff();
  cleanPolyDataFilter->PointMergingOn();
  cleanPolyDataFilter->Update();

  surface->SetVtkPolyData(cleanPolyDataFilter->GetOutput(), time);
  polydata->UnRegister(nullptr);
}

void mitk::ImageToSurfaceFilter::GenerateData()
{
  mitk::Surface *surface = this->GetOutput();
  auto *image = (mitk::Image *)GetInput();
  if (image == nullptr || !image->IsInitialized())
    mitkThrow() << "No input image set, please set an valid input image!";

  mitk::Image::RegionType outputRegion = image->GetRequestedRegion();

  int tstart = outputRegion.GetIndex(3);
  int tmax = tstart + outputRegion.GetSize(3); // GetSize()==1 - will aber 0 haben, wenn nicht zeitaufgeloest

  if ((tmax - tstart) > 0)
  {
    ProgressBar::GetInstance()->AddStepsToDo(4 * (tmax - tstart));
  }

  int t;
  for (t = tstart; t < tmax; ++t)
  {
    vtkImageData *vtkimagedata = image->GetVtkImageData(t);
    CreateSurface(t, vtkimagedata, surface, m_Threshold);
    ProgressBar::GetInstance()->Progress();
  }
}

void mitk::ImageToSurfaceFilter::SetSmoothIteration(int smoothIteration)
{
  m_SmoothIteration = smoothIteration;
}

void mitk::ImageToSurfaceFilter::SetSmoothRelaxation(float smoothRelaxation)
{
  m_SmoothRelaxation = smoothRelaxation;
}

void mitk::ImageToSurfaceFilter::SetInput(const mitk::Image *image)
{
  // Process object is not const-correct so the const_cast is required here
  this->ProcessObject::SetNthInput(0, const_cast<mitk::Image *>(image));
}

const mitk::Image *mitk::ImageToSurfaceFilter::GetInput(void)
{
  if (this->GetNumberOfInputs() < 1)
  {
    return nullptr;
  }

  return static_cast<const mitk::Image *>(this->ProcessObject::GetInput(0));
}

void mitk::ImageToSurfaceFilter::GenerateOutputInformation()
{
  mitk::Image::ConstPointer inputImage = (mitk::Image *)this->GetInput();

  // mitk::Image *inputImage = (mitk::Image*)this->GetImage();
  mitk::Surface::Pointer output = this->GetOutput();

  itkDebugMacro(<< "GenerateOutputInformation()");

  if (inputImage.IsNull())
    return;

  // Set Data
}