diff --git a/Modules/Core/include/mitkPlaneGeometryDataToSurfaceFilter.h b/Modules/Core/include/mitkPlaneGeometryDataToSurfaceFilter.h index 35a6e6b779..6f66b14c63 100644 --- a/Modules/Core/include/mitkPlaneGeometryDataToSurfaceFilter.h +++ b/Modules/Core/include/mitkPlaneGeometryDataToSurfaceFilter.h @@ -1,223 +1,220 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef mitkPlaneGeometryDataToSurfaceFilter_h #define mitkPlaneGeometryDataToSurfaceFilter_h #include "mitkGeometry3D.h" #include "mitkSurfaceSource.h" #include "vtkSystemIncludes.h" class vtkPlaneSource; class vtkTransformPolyDataFilter; class vtkCubeSource; class vtkTransform; class vtkPlane; class vtkCutter; -class vtkStripper; class vtkPolyData; class vtkPPolyDataNormals; -class vtkTriangleFilter; class vtkTextureMapToPlane; class vtkBox; class vtkClipPolyData; +class vtkContourTriangulator; namespace mitk { class PlaneGeometryData; class PlaneGeometryDataToSurfaceFilter; /** \deprecatedSince{2014_10} This class is deprecated. Please use PlaneGeometryDataToSurfaceFilter instead. */ DEPRECATED(typedef PlaneGeometryDataToSurfaceFilter Geometry2DDataToSurfaceFilter); /** \brief Superclass of all classes having a PlaneGeometryData as input and * generating Images as output * * Currently implemented for PlaneGeometry and AbstractTransformGeometry. * Currently, this class does not really have subclasses, but does the job * for itself. It checks which kind of PlaneGeometry is stored in the * PlaneGeometryData and - if it knows how - it generates the respective * Surface. Of course, this has the disadvantage that for any new type of * PlaneGeometry this class (PlaneGeometryDataToSurfaceFilter) has to be * changed/extended. The idea is to move the type specific generation code in * subclasses, and internally (within this class) use a factory to create an * instance of the required subclass and delegate the surface generation to * it. * * \sa mitk::DeformablePlane * \todo make extension easier * \ingroup Process */ class MITKCORE_EXPORT PlaneGeometryDataToSurfaceFilter : public SurfaceSource { public: mitkClassMacro(PlaneGeometryDataToSurfaceFilter, SurfaceSource); itkFactorylessNewMacro(Self); itkCloneMacro(Self); void GenerateOutputInformation() override; void GenerateData() override; const PlaneGeometryData *GetInput(void); const PlaneGeometryData *GetInput(unsigned int idx); virtual void SetInput(const PlaneGeometryData *image); using itk::ProcessObject::SetInput; virtual void SetInput(unsigned int index, const PlaneGeometryData *image); /** \brief If \a true (default), use Geometry3D::GetParametricBounds() to define the resolution in parameter space, * otherwise use m_XResolution and m_YResolution */ itkGetMacro(UseGeometryParametricBounds, bool); /** \brief If \a true (default), use Geometry3D::GetParametricBounds() to define the resolution in parameter space, * otherwise use m_XResolution and m_YResolution */ itkSetMacro(UseGeometryParametricBounds, bool); /** \brief Get x-resolution in parameter space * * The m_PlaneSource will create this many sub-rectangles * in x-direction (see vtkPlaneSource::SetXResolution) * \note Only used, when GetUseGeometryParametricBounds() is \a false, otherwise the * the x-bounds of Geometry3D::GetParametricBounds() are used. * \sa m_XResolution */ itkGetMacro(XResolution, int); /** \brief Set x-resolution in parameter space * * The m_PlaneSource will create this many sub-rectangles * in x-direction (see vtkPlaneSource::SetXResolution) * \note Only used, when GetUseGeometryParametricBounds() is \a false, otherwise the * the x-bounds of Geometry3D::GetParametricBounds() are used. * \sa m_XResolution */ itkSetMacro(XResolution, int); /** \brief Get y-resolution in parameter space * * The m_PlaneSource will create this many sub-rectangles * in y-direction (see vtkPlaneSource::SetYResolution) * \note Only used, when GetUseGeometryParametricBounds() is \a false, otherwise the * the y-bounds of Geometry3D::GetParametricBounds() are used. * \sa m_YResolution */ itkGetMacro(YResolution, int); /** \brief Set y-resolution in parameter space * * The m_PlaneSource will create this many sub-rectangles * in y-direction (see vtkPlaneSource::SetYResolution) * \note Only used, when GetUseGeometryParametricBounds() is \a false, otherwise the * the y-bounds of Geometry3D::GetParametricBounds() are used. * \sa m_YResolution */ itkSetMacro(YResolution, int); /** \brief Get whether the Surface is at the origin and placed using the Geometry * * Default is \a false, i.e., the transform of the Geometry is the identity, thus * the points within the Surface are at their final position. Otherwise * (m_PlaceByGeometry==\a true), the first cornerpoint of the created Surface is * at the origin and the actual position is determined by the transform of the * Geometry. * \sa m_PlaceByGeometry */ itkGetConstMacro(PlaceByGeometry, bool); /** \brief Set whether the Surface is at the origin and placed using the Geometry * * Default is \a false, i.e., the transform of the Geometry is the identity, thus * the points within the Surface are at their final position. Otherwise * (m_PlaceByGeometry==\a true), the first cornerpoint of the created Surface is * at the origin and the actual position is determined by the transform of the * Geometry. * \sa m_PlaceByGeometry */ itkSetMacro(PlaceByGeometry, bool); itkBooleanMacro(PlaceByGeometry); itkGetConstMacro(UseBoundingBox, bool); itkSetMacro(UseBoundingBox, bool); itkBooleanMacro(UseBoundingBox); void SetBoundingBox(const BoundingBox *boundingBox); const BoundingBox *GetBoundingBox() const; protected: PlaneGeometryDataToSurfaceFilter(); ~PlaneGeometryDataToSurfaceFilter() override; /** \brief Source to create the vtk-representation of the parameter space rectangle of the PlaneGeometry */ vtkPlaneSource *m_PlaneSource; /** \brief Filter to create the vtk-representation of the PlaneGeometry, which is a * transformation of the m_PlaneSource */ vtkTransformPolyDataFilter *m_VtkTransformPlaneFilter; /** \brief If \a true, use Geometry3D::GetParametricBounds() to define the resolution in parameter space, * otherwise use m_XResolution and m_YResolution */ bool m_UseGeometryParametricBounds; /** \brief X-resolution in parameter space * * The m_PlaneSource will create this many sub-rectangles * in x-direction (see vtkPlaneSource::SetXResolution) * \note Only used, when GetUseGeometryParametricBounds() is \a false, otherwise the * the x-bounds of Geometry3D::GetParametricBounds() are used. * \sa m_XResolution */ int m_XResolution; /** \brief Y-resolution in parameter space * * The m_PlaneSource will create this many sub-rectangles * in y-direction (see vtkPlaneSource::SetYResolution) * \note Only used, when GetUseGeometryParametricBounds() is \a false, otherwise the * the y-bounds of Geometry3D::GetParametricBounds() are used. */ int m_YResolution; /** \brief Define whether the Surface is at the origin and placed using the Geometry * * Default is \a false, i.e., the transform of the Geometry is the identity, thus * the points within the Surface are at their final position. Otherwise * (m_PlaceByGeometry==\a true), the first cornerpoint of the created Surface is * at the origin and the actual position is determined by the transform of the * Geometry. */ bool m_PlaceByGeometry; bool m_UseBoundingBox; BoundingBox::ConstPointer m_BoundingBox; vtkCubeSource *m_CubeSource; vtkTransform *m_Transform; vtkTransformPolyDataFilter *m_PolyDataTransformer; vtkPlane *m_Plane; vtkCutter *m_PlaneCutter; - vtkStripper *m_PlaneStripper; - vtkPolyData *m_PlanePolyData; + vtkContourTriangulator *m_ContourTriangulator; vtkPPolyDataNormals *m_NormalsUpdater; - vtkTriangleFilter *m_PlaneTriangler; vtkTextureMapToPlane *m_TextureMapToPlane; vtkBox *m_Box; vtkClipPolyData *m_PlaneClipper; }; } // namespace mitk #endif diff --git a/Modules/Core/src/Algorithms/mitkPlaneGeometryDataToSurfaceFilter.cpp b/Modules/Core/src/Algorithms/mitkPlaneGeometryDataToSurfaceFilter.cpp index 988afe14aa..a7510e655b 100644 --- a/Modules/Core/src/Algorithms/mitkPlaneGeometryDataToSurfaceFilter.cpp +++ b/Modules/Core/src/Algorithms/mitkPlaneGeometryDataToSurfaceFilter.cpp @@ -1,391 +1,382 @@ /*============================================================================ 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 "mitkPlaneGeometryDataToSurfaceFilter.h" #include "mitkAbstractTransformGeometry.h" #include "mitkGeometry3D.h" #include "mitkPlaneGeometry.h" #include "mitkPlaneGeometryData.h" #include "mitkSurface.h" #include <vtkPlaneSource.h> #include <vtkPolyData.h> #include <vtkTransformPolyDataFilter.h> #include <vtkBox.h> #include <vtkClipPolyData.h> #include <vtkCubeSource.h> #include <vtkCutter.h> #include <vtkGeneralTransform.h> #include <vtkPPolyDataNormals.h> #include <vtkPlane.h> -#include <vtkStripper.h> #include <vtkTextureMapToPlane.h> #include <vtkTransform.h> #include <vtkTransformPolyDataFilter.h> -#include <vtkTriangleFilter.h> +#include <vtkContourTriangulator.h> mitk::PlaneGeometryDataToSurfaceFilter::PlaneGeometryDataToSurfaceFilter() : m_UseGeometryParametricBounds(true), m_XResolution(10), m_YResolution(10), m_PlaceByGeometry(false), m_UseBoundingBox(false) { m_PlaneSource = vtkPlaneSource::New(); m_Transform = vtkTransform::New(); m_CubeSource = vtkCubeSource::New(); m_PolyDataTransformer = vtkTransformPolyDataFilter::New(); m_Plane = vtkPlane::New(); m_PlaneCutter = vtkCutter::New(); - m_PlaneStripper = vtkStripper::New(); - m_PlanePolyData = vtkPolyData::New(); + m_ContourTriangulator = vtkContourTriangulator::New(); m_NormalsUpdater = vtkPPolyDataNormals::New(); - m_PlaneTriangler = vtkTriangleFilter::New(); m_TextureMapToPlane = vtkTextureMapToPlane::New(); m_Box = vtkBox::New(); m_PlaneClipper = vtkClipPolyData::New(); m_VtkTransformPlaneFilter = vtkTransformPolyDataFilter::New(); m_VtkTransformPlaneFilter->SetInputConnection(m_PlaneSource->GetOutputPort()); } mitk::PlaneGeometryDataToSurfaceFilter::~PlaneGeometryDataToSurfaceFilter() { m_PlaneSource->Delete(); m_Transform->Delete(); m_CubeSource->Delete(); m_PolyDataTransformer->Delete(); m_Plane->Delete(); m_PlaneCutter->Delete(); - m_PlaneStripper->Delete(); - m_PlanePolyData->Delete(); + m_ContourTriangulator->Delete(); m_NormalsUpdater->Delete(); - m_PlaneTriangler->Delete(); m_TextureMapToPlane->Delete(); m_Box->Delete(); m_PlaneClipper->Delete(); m_VtkTransformPlaneFilter->Delete(); } void mitk::PlaneGeometryDataToSurfaceFilter::GenerateOutputInformation() { mitk::PlaneGeometryData::ConstPointer input = this->GetInput(); mitk::Surface::Pointer output = this->GetOutput(); if (input.IsNull() || (input->GetPlaneGeometry() == nullptr) || (input->GetPlaneGeometry()->IsValid() == false) || (m_UseBoundingBox && (m_BoundingBox.IsNull() || (m_BoundingBox->GetDiagonalLength2() < mitk::eps)))) { return; } Point3D origin; Point3D right, bottom; vtkPolyData *planeSurface = nullptr; // Does the PlaneGeometryData contain an AbstractTransformGeometry? if (auto *abstractGeometry = dynamic_cast<AbstractTransformGeometry *>(input->GetPlaneGeometry())) { // In the case of an AbstractTransformGeometry (which holds a possibly // non-rigid transform), we proceed slightly differently: since the // plane can be arbitrarily deformed, we need to transform it by the // abstract transform before clipping it. The setup for this is partially // done in the constructor. origin = abstractGeometry->GetPlane()->GetOrigin(); right = origin + abstractGeometry->GetPlane()->GetAxisVector(0); bottom = origin + abstractGeometry->GetPlane()->GetAxisVector(1); // Define the plane m_PlaneSource->SetOrigin(origin[0], origin[1], origin[2]); m_PlaneSource->SetPoint1(right[0], right[1], right[2]); m_PlaneSource->SetPoint2(bottom[0], bottom[1], bottom[2]); // Set the plane's resolution (unlike for non-deformable planes, the plane // grid needs to have a certain resolution so that the deformation has the // desired effect). if (m_UseGeometryParametricBounds) { m_PlaneSource->SetXResolution((int)abstractGeometry->GetParametricExtent(0)); m_PlaneSource->SetYResolution((int)abstractGeometry->GetParametricExtent(1)); } else { m_PlaneSource->SetXResolution(m_XResolution); m_PlaneSource->SetYResolution(m_YResolution); } if (m_PlaceByGeometry) { // Let the output use the input geometry to appropriately transform the // coordinate system. mitk::Geometry3D::TransformType *affineTransform = abstractGeometry->GetIndexToWorldTransform(); TimeGeometry *timeGeometry = output->GetTimeGeometry(); BaseGeometry *g3d = timeGeometry->GetGeometryForTimeStep(0); g3d->SetIndexToWorldTransform(affineTransform); vtkGeneralTransform *composedResliceTransform = vtkGeneralTransform::New(); composedResliceTransform->Identity(); composedResliceTransform->Concatenate(abstractGeometry->GetVtkTransform()->GetLinearInverse()); composedResliceTransform->Concatenate(abstractGeometry->GetVtkAbstractTransform()); // Use the non-rigid transform for transforming the plane. m_VtkTransformPlaneFilter->SetTransform(composedResliceTransform); } else { // Use the non-rigid transform for transforming the plane. m_VtkTransformPlaneFilter->SetTransform(abstractGeometry->GetVtkAbstractTransform()); } if (m_UseBoundingBox) { mitk::BoundingBox::PointType boundingBoxMin = m_BoundingBox->GetMinimum(); mitk::BoundingBox::PointType boundingBoxMax = m_BoundingBox->GetMaximum(); // mitk::BoundingBox::PointType boundingBoxCenter = m_BoundingBox->GetCenter(); m_Box->SetXMin(boundingBoxMin[0], boundingBoxMin[1], boundingBoxMin[2]); m_Box->SetXMax(boundingBoxMax[0], boundingBoxMax[1], boundingBoxMax[2]); } else { // Plane will not be clipped m_Box->SetXMin(-10000.0, -10000.0, -10000.0); m_Box->SetXMax(10000.0, 10000.0, 10000.0); } m_Transform->Identity(); m_Transform->Concatenate(input->GetPlaneGeometry()->GetVtkTransform()); m_Transform->PreMultiply(); m_Box->SetTransform(m_Transform); m_PlaneClipper->SetInputConnection(m_VtkTransformPlaneFilter->GetOutputPort()); m_PlaneClipper->SetClipFunction(m_Box); m_PlaneClipper->GenerateClippedOutputOff(); // important to NOT generate normals data for clipped part m_PlaneClipper->InsideOutOn(); m_PlaneClipper->SetValue(0.0); m_PlaneClipper->Update(); planeSurface = m_PlaneClipper->GetOutput(); } // Does the PlaneGeometryData contain a PlaneGeometry? else if (dynamic_cast<PlaneGeometry *>(input->GetPlaneGeometry()) != nullptr) { auto *planeGeometry = dynamic_cast<PlaneGeometry *>(input->GetPlaneGeometry()); if (m_PlaceByGeometry) { // Let the output use the input geometry to appropriately transform the // coordinate system. mitk::Geometry3D::TransformType *affineTransform = planeGeometry->GetIndexToWorldTransform(); TimeGeometry *timeGeometry = output->GetTimeGeometry(); BaseGeometry *geometrie3d = timeGeometry->GetGeometryForTimeStep(0); geometrie3d->SetIndexToWorldTransform(affineTransform); } if (!m_UseBoundingBox) { // We do not have a bounding box, so no clipping is required. if (m_PlaceByGeometry) { // Derive coordinate axes and origin from input geometry extent origin.Fill(0.0); FillVector3D(right, planeGeometry->GetExtent(0), 0.0, 0.0); FillVector3D(bottom, 0.0, planeGeometry->GetExtent(1), 0.0); } else { // Take the coordinate axes and origin directly from the input geometry. origin = planeGeometry->GetOrigin(); right = planeGeometry->GetCornerPoint(false, true); bottom = planeGeometry->GetCornerPoint(true, false); } // Since the plane is planar, there is no need to subdivide the grid // (cf. AbstractTransformGeometry case) m_PlaneSource->SetXResolution(1); m_PlaneSource->SetYResolution(1); m_PlaneSource->SetOrigin(origin[0], origin[1], origin[2]); m_PlaneSource->SetPoint1(right[0], right[1], right[2]); m_PlaneSource->SetPoint2(bottom[0], bottom[1], bottom[2]); m_PlaneSource->Update(); planeSurface = m_PlaneSource->GetOutput(); } else { // Set up a cube with the extent and origin of the bounding box. This // cube will be clipped by a plane later on. The intersection of the // cube and the plane will be the surface we are interested in. Note // that the bounding box needs to be explicitly specified by the user // of this class, since it is not necessarily clear from the data // available herein which bounding box to use. In most cases, this // would be the bounding box of the input geometry's reference // geometry, but this is not an inevitable requirement. mitk::BoundingBox::PointType boundingBoxMin = m_BoundingBox->GetMinimum(); mitk::BoundingBox::PointType boundingBoxMax = m_BoundingBox->GetMaximum(); mitk::BoundingBox::PointType boundingBoxCenter = m_BoundingBox->GetCenter(); m_CubeSource->SetXLength(boundingBoxMax[0] - boundingBoxMin[0]); m_CubeSource->SetYLength(boundingBoxMax[1] - boundingBoxMin[1]); m_CubeSource->SetZLength(boundingBoxMax[2] - boundingBoxMin[2]); m_CubeSource->SetCenter(boundingBoxCenter[0], boundingBoxCenter[1], boundingBoxCenter[2]); // Now we have to transform the cube, so that it will cut our plane // appropriately. (As can be seen below, the plane corresponds to the // z-plane in the coordinate system and is *not* transformed.) Therefore, // we get the inverse of the plane geometry's transform and concatenate // it with the transform of the reference geometry, if available. m_Transform->Identity(); m_Transform->Concatenate(planeGeometry->GetVtkTransform()->GetLinearInverse()); const BaseGeometry *referenceGeometry = planeGeometry->GetReferenceGeometry(); if (referenceGeometry) { m_Transform->Concatenate(referenceGeometry->GetVtkTransform()); } // Transform the cube accordingly (s.a.) m_PolyDataTransformer->SetInputConnection(m_CubeSource->GetOutputPort()); m_PolyDataTransformer->SetTransform(m_Transform); // Initialize the plane to clip the cube with, as lying on the z-plane m_Plane->SetOrigin(0.0, 0.0, 0.0); m_Plane->SetNormal(0.0, 0.0, 1.0); // Cut the plane with the cube. m_PlaneCutter->SetInputConnection(m_PolyDataTransformer->GetOutputPort()); + m_PlaneCutter->GenerateTrianglesOff(); m_PlaneCutter->SetCutFunction(m_Plane); // The output of the cutter must be converted into appropriate poly data. - m_PlaneStripper->SetInputConnection(m_PlaneCutter->GetOutputPort()); - m_PlaneStripper->Update(); + m_ContourTriangulator->SetInputConnection(m_PlaneCutter->GetOutputPort()); + m_ContourTriangulator->Update(); - if (m_PlaneStripper->GetOutput()->GetNumberOfPoints() < 3) + if (m_ContourTriangulator->GetTriangulationError()) { return; } - m_PlanePolyData->SetPoints(m_PlaneStripper->GetOutput()->GetPoints()); - m_PlanePolyData->SetPolys(m_PlaneStripper->GetOutput()->GetLines()); - - m_PlaneTriangler->SetInputData(m_PlanePolyData); - // Get bounds of the resulting surface and use it to generate the texture // mapping information - m_PlaneTriangler->Update(); - m_PlaneTriangler->GetOutput()->ComputeBounds(); - double *surfaceBounds = m_PlaneTriangler->GetOutput()->GetBounds(); + m_ContourTriangulator->Update(); + m_ContourTriangulator->GetOutput()->ComputeBounds(); + double *surfaceBounds = m_ContourTriangulator->GetOutput()->GetBounds(); origin[0] = surfaceBounds[0]; origin[1] = surfaceBounds[2]; origin[2] = surfaceBounds[4]; right[0] = surfaceBounds[1]; right[1] = surfaceBounds[2]; right[2] = surfaceBounds[4]; bottom[0] = surfaceBounds[0]; bottom[1] = surfaceBounds[3]; bottom[2] = surfaceBounds[4]; // Now we tell the data how it shall be textured afterwards; // description see above. - m_TextureMapToPlane->SetInputConnection(m_PlaneTriangler->GetOutputPort()); + m_TextureMapToPlane->SetInputConnection(m_ContourTriangulator->GetOutputPort()); m_TextureMapToPlane->AutomaticPlaneGenerationOn(); m_TextureMapToPlane->SetOrigin(origin[0], origin[1], origin[2]); m_TextureMapToPlane->SetPoint1(right[0], right[1], right[2]); m_TextureMapToPlane->SetPoint2(bottom[0], bottom[1], bottom[2]); // Need to call update so that output data and bounds are immediately // available m_TextureMapToPlane->Update(); // Return the output of this generation process planeSurface = dynamic_cast<vtkPolyData *>(m_TextureMapToPlane->GetOutput()); } } m_NormalsUpdater->SetInputData(planeSurface); m_NormalsUpdater->AutoOrientNormalsOn(); // that's the trick! Brings consistency between // normals direction and front/back faces direction (see bug 1440) m_NormalsUpdater->ComputePointNormalsOn(); m_NormalsUpdater->Update(); output->SetVtkPolyData(m_NormalsUpdater->GetOutput()); output->CalculateBoundingBox(); } void mitk::PlaneGeometryDataToSurfaceFilter::GenerateData() { mitk::Surface::Pointer output = this->GetOutput(); if (output.IsNull()) return; if (output->GetVtkPolyData() == nullptr) return; // output->GetVtkPolyData()->Update(); //VTK6_TODO vtk pipeline } const mitk::PlaneGeometryData *mitk::PlaneGeometryDataToSurfaceFilter::GetInput() { if (this->GetNumberOfInputs() < 1) { return nullptr; } return static_cast<const mitk::PlaneGeometryData *>(this->ProcessObject::GetInput(0)); } const mitk::PlaneGeometryData *mitk::PlaneGeometryDataToSurfaceFilter::GetInput(unsigned int idx) { return static_cast<const mitk::PlaneGeometryData *>(this->ProcessObject::GetInput(idx)); } void mitk::PlaneGeometryDataToSurfaceFilter::SetInput(const mitk::PlaneGeometryData *input) { // Process object is not const-correct so the const_cast is required here this->ProcessObject::SetNthInput(0, const_cast<mitk::PlaneGeometryData *>(input)); } void mitk::PlaneGeometryDataToSurfaceFilter::SetInput(unsigned int index, const mitk::PlaneGeometryData *input) { if (index + 1 > this->GetNumberOfInputs()) { this->SetNumberOfRequiredInputs(index + 1); } // Process object is not const-correct so the const_cast is required here this->ProcessObject::SetNthInput(index, const_cast<mitk::PlaneGeometryData *>(input)); } void mitk::PlaneGeometryDataToSurfaceFilter::SetBoundingBox(const mitk::BoundingBox *boundingBox) { m_BoundingBox = boundingBox; this->UseBoundingBoxOn(); } const mitk::BoundingBox *mitk::PlaneGeometryDataToSurfaceFilter::GetBoundingBox() const { return m_BoundingBox.GetPointer(); }