diff --git a/Modules/MapperExt/src/mitkUnstructuredGridMapper2D.cpp b/Modules/MapperExt/src/mitkUnstructuredGridMapper2D.cpp index a4e89fbf06..a4db90ea37 100644 --- a/Modules/MapperExt/src/mitkUnstructuredGridMapper2D.cpp +++ b/Modules/MapperExt/src/mitkUnstructuredGridMapper2D.cpp @@ -1,555 +1,555 @@ /*============================================================================ 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 "mitkUnstructuredGridMapper2D.h" #include #include "mitkAbstractTransformGeometry.h" #include "mitkBaseRenderer.h" #include "mitkColorProperty.h" #include "mitkPlaneGeometry.h" #include "mitkProperties.h" #include "mitkTransferFunction.h" #include "mitkTransferFunctionProperty.h" #include "mitkUnstructuredGrid.h" #include "mitkVtkMapper3D.h" #include "mitkVtkScalarModeProperty.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "vtkPointSetSlicer.h" void mitk::UnstructuredGridMapper2D::GenerateDataForRenderer(mitk::BaseRenderer *renderer) { BaseLocalStorage *ls = m_LSH.GetLocalStorage(renderer); bool needGenerateData = ls->IsGenerateDataRequired(renderer, this, GetDataNode()); if (needGenerateData) { ls->UpdateGenerateDataTime(); mitk::DataNode::ConstPointer node = this->GetDataNode(); if (node.IsNull()) return; if (!node->GetProperty(m_ScalarMode, "scalar mode")) { m_ScalarMode = mitk::VtkScalarModeProperty::New(0); } if (!node->GetProperty(m_ScalarVisibility, "scalar visibility")) { m_ScalarVisibility = mitk::BoolProperty::New(true); } if (!node->GetProperty(m_Outline, "outline polygons")) { m_Outline = mitk::BoolProperty::New(false); } if (!node->GetProperty(m_Color, "color")) { m_Color = mitk::ColorProperty::New(1.0f, 1.0f, 1.0f); } if (!node->GetProperty(m_LineWidth, "line width")) { m_LineWidth = mitk::IntProperty::New(1); } } mitk::BaseData::Pointer input = GetDataNode()->GetData(); assert(input); input->Update(); if (m_VtkPointSet) m_VtkPointSet->UnRegister(nullptr); m_VtkPointSet = this->GetVtkPointSet(renderer, this->GetTimestep()); assert(m_VtkPointSet); m_VtkPointSet->Register(nullptr); if (m_ScalarVisibility->GetValue()) { mitk::DataNode::ConstPointer node = this->GetDataNode(); mitk::TransferFunctionProperty::Pointer transferFuncProp; node->GetProperty(transferFuncProp, "TransferFunction", renderer); if (transferFuncProp.IsNotNull()) { mitk::TransferFunction::Pointer tf = transferFuncProp->GetValue(); if (m_ScalarsToColors) m_ScalarsToColors->UnRegister(nullptr); m_ScalarsToColors = static_cast(tf->GetColorTransferFunction()); m_ScalarsToColors->Register(nullptr); if (m_ScalarsToOpacity) m_ScalarsToOpacity->UnRegister(nullptr); m_ScalarsToOpacity = tf->GetScalarOpacityFunction(); m_ScalarsToOpacity->Register(nullptr); } else { if (m_ScalarsToColors) m_ScalarsToColors->UnRegister(nullptr); m_ScalarsToColors = this->GetVtkLUT(renderer); assert(m_ScalarsToColors); m_ScalarsToColors->Register(nullptr); float opacity; node->GetOpacity(opacity, renderer); if (m_ScalarsToOpacity) m_ScalarsToOpacity->UnRegister(nullptr); m_ScalarsToOpacity = vtkPiecewiseFunction::New(); double range[2]; m_VtkPointSet->GetScalarRange(range); m_ScalarsToOpacity->AddSegment(range[0], opacity, range[1], opacity); } } } void mitk::UnstructuredGridMapper2D::Paint(mitk::BaseRenderer *renderer) { bool visible = true; GetDataNode()->GetVisibility(visible, renderer, "visible"); if (!visible) return; vtkLinearTransform *vtktransform = GetDataNode()->GetVtkTransform(); vtkLinearTransform *inversetransform = vtktransform->GetLinearInverse(); PlaneGeometry::ConstPointer worldGeometry = renderer->GetCurrentWorldPlaneGeometry(); PlaneGeometry::ConstPointer worldPlaneGeometry = dynamic_cast(worldGeometry.GetPointer()); Point3D point; Vector3D normal; if (worldPlaneGeometry.IsNotNull()) { // set up vtkPlane according to worldGeometry point = worldPlaneGeometry->GetOrigin(); normal = worldPlaneGeometry->GetNormal(); normal.Normalize(); m_Plane->SetTransform((vtkAbstractTransform *)nullptr); } else { //@FIXME: does not work correctly. Does m_Plane->SetTransform really transforms a "plane plane" into a "curved //plane"? return; AbstractTransformGeometry::ConstPointer worldAbstractGeometry = dynamic_cast(renderer->GetCurrentWorldPlaneGeometry()); if (worldAbstractGeometry.IsNotNull()) { // set up vtkPlane according to worldGeometry point = worldAbstractGeometry->GetParametricBoundingBox()->GetMinimum(); FillVector3D(normal, 0, 0, 1); m_Plane->SetTransform(worldAbstractGeometry->GetVtkAbstractTransform()->GetInverse()); } else return; } double vp[3], vnormal[3]; vnl2vtk(point.GetVnlVector(), vp); vnl2vtk(normal.GetVnlVector(), vnormal); // normally, we would need to transform the surface and cut the transformed surface with the cutter. // This might be quite slow. Thus, the idea is, to perform an inverse transform of the plane instead. //@todo It probably does not work for scaling operations yet:scaling operations have to be // dealed with after the cut is performed by scaling the contour. inversetransform->TransformPoint(vp, vp); inversetransform->TransformNormalAtPoint(vp, vnormal, vnormal); m_Plane->SetOrigin(vp); m_Plane->SetNormal(vnormal); // set data into cutter m_Slicer->SetInputData(m_VtkPointSet); // m_Cutter->GenerateCutScalarsOff(); // m_Cutter->SetSortByToSortByCell(); // calculate the cut m_Slicer->Update(); // apply color and opacity read from the PropertyList ApplyColorAndOpacityProperties(renderer); // traverse the cut contour vtkPolyData *contour = m_Slicer->GetOutput(); vtkPoints *vpoints = contour->GetPoints(); vtkCellArray *vlines = contour->GetLines(); vtkCellArray *vpolys = contour->GetPolys(); vtkPointData *vpointdata = contour->GetPointData(); vtkDataArray *vscalars = vpointdata->GetScalars(); vtkCellData *vcelldata = contour->GetCellData(); vtkDataArray *vcellscalars = vcelldata->GetScalars(); const int numberOfLines = contour->GetNumberOfLines(); const int numberOfPolys = contour->GetNumberOfPolys(); const bool useCellData = m_ScalarMode->GetVtkScalarMode() == VTK_SCALAR_MODE_DEFAULT || m_ScalarMode->GetVtkScalarMode() == VTK_SCALAR_MODE_USE_CELL_DATA; const bool usePointData = m_ScalarMode->GetVtkScalarMode() == VTK_SCALAR_MODE_USE_POINT_DATA; Point3D p; Point2D p2d; vlines->InitTraversal(); vpolys->InitTraversal(); mitk::Color outlineColor = m_Color->GetColor(); glLineWidth((float)m_LineWidth->GetValue()); for (int i = 0; i < numberOfLines; ++i) { - vtkIdType *cell(nullptr); + const vtkIdType *cell(nullptr); vtkIdType cellSize(0); vlines->GetNextCell(cellSize, cell); float rgba[4] = {outlineColor[0], outlineColor[1], outlineColor[2], 1.0f}; if (m_ScalarVisibility->GetValue() && vcellscalars) { if (useCellData) { // color each cell according to cell data double scalar = vcellscalars->GetComponent(i, 0); double rgb[3] = {1.0f, 1.0f, 1.0f}; m_ScalarsToColors->GetColor(scalar, rgb); rgba[0] = (float)rgb[0]; rgba[1] = (float)rgb[1]; rgba[2] = (float)rgb[2]; rgba[3] = (float)m_ScalarsToOpacity->GetValue(scalar); } else if (usePointData) { double scalar = vscalars->GetComponent(i, 0); double rgb[3] = {1.0f, 1.0f, 1.0f}; m_ScalarsToColors->GetColor(scalar, rgb); rgba[0] = (float)rgb[0]; rgba[1] = (float)rgb[1]; rgba[2] = (float)rgb[2]; rgba[3] = (float)m_ScalarsToOpacity->GetValue(scalar); } } glColor4fv(rgba); glBegin(GL_LINE_LOOP); for (int j = 0; j < cellSize; ++j) { vpoints->GetPoint(cell[j], vp); // take transformation via vtktransform into account vtktransform->TransformPoint(vp, vp); vtk2itk(vp, p); // convert 3D point (in mm) to display coordinates (units ) renderer->WorldToDisplay(p, p2d); // convert display coordinates ( (0,0) is top-left ) in GL coordinates ( (0,0) is bottom-left ) // p2d[1]=toGL-p2d[1]; // add the current vertex to the line glVertex2f(p2d[0], p2d[1]); } glEnd(); } bool polyOutline = m_Outline->GetValue(); bool scalarVisibility = m_ScalarVisibility->GetValue(); // cache the transformed points // a fixed size array is way faster than 'new' // slices through 3d cells usually do not generated // polygons with more than 6 vertices const int maxPolySize = 10; auto *cachedPoints = new Point2D[maxPolySize * numberOfPolys]; glEnable(GL_BLEND); glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); // only draw polygons if there are cell scalars // or the outline property is set to true if (scalarVisibility && vcellscalars) { glPolygonMode(GL_FRONT_AND_BACK, GL_FILL); for (int i = 0; i < numberOfPolys; ++i) { - vtkIdType *cell(nullptr); + const vtkIdType *cell(nullptr); vtkIdType cellSize(0); vpolys->GetNextCell(cellSize, cell); float rgba[4] = {1.0f, 1.0f, 1.0f, 0}; if (scalarVisibility && vcellscalars) { if (useCellData) { // color each cell according to cell data double scalar = vcellscalars->GetComponent(i + numberOfLines, 0); double rgb[3] = {1.0f, 1.0f, 1.0f}; m_ScalarsToColors->GetColor(scalar, rgb); rgba[0] = (float)rgb[0]; rgba[1] = (float)rgb[1]; rgba[2] = (float)rgb[2]; rgba[3] = (float)m_ScalarsToOpacity->GetValue(scalar); } else if (usePointData) { double scalar = vscalars->GetComponent(i, 0); double rgb[3] = {1.0f, 1.0f, 1.0f}; m_ScalarsToColors->GetColor(scalar, rgb); rgba[0] = (float)rgb[0]; rgba[1] = (float)rgb[1]; rgba[2] = (float)rgb[2]; rgba[3] = (float)m_ScalarsToOpacity->GetValue(scalar); } } glColor4fv(rgba); glBegin(GL_POLYGON); for (int j = 0; j < cellSize; ++j) { vpoints->GetPoint(cell[j], vp); // take transformation via vtktransform into account vtktransform->TransformPoint(vp, vp); vtk2itk(vp, p); // convert 3D point (in mm) to display coordinates (units ) renderer->WorldToDisplay(p, p2d); // convert display coordinates ( (0,0) is top-left ) in GL coordinates ( (0,0) is bottom-left ) // p2d[1]=toGL-p2d[1]; cachedPoints[i * 10 + j][0] = p2d[0]; cachedPoints[i * 10 + j][1] = p2d[1]; // add the current vertex to the line glVertex2f(p2d[0], p2d[1]); } glEnd(); } if (polyOutline) { vpolys->InitTraversal(); glColor4f(outlineColor[0], outlineColor[1], outlineColor[2], 1.0f); glPolygonMode(GL_FRONT_AND_BACK, GL_LINE); for (int i = 0; i < numberOfPolys; ++i) { - vtkIdType *cell(nullptr); + const vtkIdType *cell(nullptr); vtkIdType cellSize(0); vpolys->GetNextCell(cellSize, cell); glBegin(GL_POLYGON); // glPolygonOffset(1.0, 1.0); for (int j = 0; j < cellSize; ++j) { // add the current vertex to the line glVertex2f(cachedPoints[i * 10 + j][0], cachedPoints[i * 10 + j][1]); } glEnd(); } } } glDisable(GL_BLEND); delete[] cachedPoints; } vtkAbstractMapper3D *mitk::UnstructuredGridMapper2D::GetVtkAbstractMapper3D(mitk::BaseRenderer *renderer) { // MITK_INFO << "GETVTKABSTRACTMAPPER3D\n"; mitk::DataNode::ConstPointer node = this->GetDataNode(); if (node.IsNull()) return nullptr; mitk::VtkMapper::Pointer mitkMapper = dynamic_cast(node->GetMapper(2)); if (mitkMapper.IsNull()) { return nullptr; } mitkMapper->Update(renderer); auto *assembly = dynamic_cast(mitkMapper->GetVtkProp(renderer)); if (assembly) { vtkProp3DCollection *collection = assembly->GetParts(); collection->InitTraversal(); vtkProp3D *prop3d = nullptr; do { prop3d = collection->GetNextProp3D(); auto *actor = dynamic_cast(prop3d); if (actor) { return dynamic_cast(actor->GetMapper()); } auto *volume = dynamic_cast(prop3d); if (volume) { return dynamic_cast(volume->GetMapper()); } } while (prop3d != collection->GetLastProp3D()); } else { auto *actor = dynamic_cast(mitkMapper->GetVtkProp(renderer)); if (actor) { return dynamic_cast(actor->GetMapper()); } auto *volume = dynamic_cast(mitkMapper->GetVtkProp(renderer)); if (volume) { return dynamic_cast(volume->GetMapper()); } } return nullptr; } vtkPointSet *mitk::UnstructuredGridMapper2D::GetVtkPointSet(mitk::BaseRenderer *renderer, int time) { // MITK_INFO << "GETVTKPOINTSET\n"; vtkAbstractMapper3D *abstractMapper = GetVtkAbstractMapper3D(renderer); if (abstractMapper == nullptr) { // try to get data from the node mitk::DataNode::ConstPointer node = this->GetDataNode(); if (node.IsNull()) return nullptr; mitk::BaseData::Pointer data = node->GetData(); mitk::UnstructuredGrid::Pointer grid = dynamic_cast(data.GetPointer()); if (!grid.IsNull()) return static_cast(grid->GetVtkUnstructuredGrid(time)); return nullptr; } else { auto *mapper = dynamic_cast(abstractMapper); if (mapper) { return dynamic_cast(mapper->GetInput()); } auto *volMapper = dynamic_cast(abstractMapper); if (volMapper) { return dynamic_cast(volMapper->GetDataSetInput()); } } return nullptr; } vtkScalarsToColors *mitk::UnstructuredGridMapper2D::GetVtkLUT(mitk::BaseRenderer *renderer) { // MITK_INFO << "GETVTKLUT\n"; auto *mapper = dynamic_cast(GetVtkAbstractMapper3D(renderer)); if (mapper) return mapper->GetLookupTable(); else { mitk::DataNode::ConstPointer node = this->GetDataNode(); if (node.IsNull()) return nullptr; mitk::VtkMapper::Pointer mitkMapper = dynamic_cast(node->GetMapper(2)); if (mitkMapper.IsNull()) { // MITK_INFO << "mitkMapper is null\n"; return nullptr; } mitkMapper->Update(renderer); auto *volume = dynamic_cast(mitkMapper->GetVtkProp(renderer)); if (volume) { // MITK_INFO << "found volume prop\n"; return static_cast(volume->GetProperty()->GetRGBTransferFunction()); } auto *assembly = dynamic_cast(mitkMapper->GetVtkProp(renderer)); if (assembly) { // MITK_INFO << "found assembly prop\n"; mitk::TransferFunctionProperty::Pointer transferFuncProp; node->GetProperty(transferFuncProp, "TransferFunction", nullptr); if (transferFuncProp.IsNotNull()) { MITK_INFO << "return colortransferfunction\n"; return static_cast(transferFuncProp->GetValue()->GetColorTransferFunction()); } } return nullptr; } } bool mitk::UnstructuredGridMapper2D::IsConvertibleToVtkPointSet(mitk::BaseRenderer *renderer) { return (GetVtkPointSet(renderer, this->GetTimestep()) != nullptr); } mitk::UnstructuredGridMapper2D::UnstructuredGridMapper2D() { m_Plane = vtkPlane::New(); m_Slicer = vtkPointSetSlicer::New(); m_Slicer->SetSlicePlane(m_Plane); m_ScalarsToColors = nullptr; m_ScalarsToOpacity = nullptr; m_VtkPointSet = nullptr; // m_LUT = vtkLookupTable::New(); // m_LUT->SetTableRange( 0, 255 ); // m_LUT->SetNumberOfColors( 255 ); // m_LUT->SetRampToLinear (); // m_LUT->Build(); } mitk::UnstructuredGridMapper2D::~UnstructuredGridMapper2D() { m_Slicer->Delete(); m_Plane->Delete(); if (m_ScalarsToOpacity != nullptr) m_ScalarsToOpacity->UnRegister(nullptr); if (m_ScalarsToColors != nullptr) m_ScalarsToColors->UnRegister(nullptr); if (m_VtkPointSet != nullptr) m_VtkPointSet->UnRegister(nullptr); } diff --git a/Modules/SurfaceInterpolation/mitkComputeContourSetNormalsFilter.cpp b/Modules/SurfaceInterpolation/mitkComputeContourSetNormalsFilter.cpp index 5793bd362f..7589742cdf 100644 --- a/Modules/SurfaceInterpolation/mitkComputeContourSetNormalsFilter.cpp +++ b/Modules/SurfaceInterpolation/mitkComputeContourSetNormalsFilter.cpp @@ -1,343 +1,343 @@ /*============================================================================ 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 "mitkComputeContourSetNormalsFilter.h" #include "mitkIOUtil.h" #include "mitkImagePixelReadAccessor.h" mitk::ComputeContourSetNormalsFilter::ComputeContourSetNormalsFilter() : m_SegmentationBinaryImage(nullptr), m_MaxSpacing(5), m_NegativeNormalCounter(0), m_PositiveNormalCounter(0), m_UseProgressBar(false), m_ProgressStepSize(1) { mitk::Surface::Pointer output = mitk::Surface::New(); this->SetNthOutput(0, output.GetPointer()); } mitk::ComputeContourSetNormalsFilter::~ComputeContourSetNormalsFilter() { } void mitk::ComputeContourSetNormalsFilter::GenerateData() { unsigned int numberOfInputs = this->GetNumberOfIndexedInputs(); // Iterating over each input for (unsigned int i = 0; i < numberOfInputs; i++) { // Getting the inputs polydata and polygons auto *currentSurface = this->GetInput(i); vtkPolyData *polyData = currentSurface->GetVtkPolyData(); vtkSmartPointer existingPolys = polyData->GetPolys(); vtkSmartPointer existingPoints = polyData->GetPoints(); existingPolys->InitTraversal(); - vtkIdType *cell(nullptr); + const vtkIdType *cell(nullptr); vtkIdType cellSize(0); // The array that contains all the vertex normals of the current polygon vtkSmartPointer normals = vtkSmartPointer::New(); normals->SetNumberOfComponents(3); normals->SetNumberOfTuples(polyData->GetNumberOfPoints()); // If the current contour is an inner contour then the direction is -1 // A contour lies inside another one if the pixel values in the direction of the normal is 1 m_NegativeNormalCounter = 0; m_PositiveNormalCounter = 0; vtkIdType offSet(0); // Iterating over each polygon for (existingPolys->InitTraversal(); existingPolys->GetNextCell(cellSize, cell);) { if (cellSize < 3) continue; // First we calculate the current polygon's normal double polygonNormal[3] = {0.0}; double p1[3]; double p2[3]; double v1[3]; double v2[3]; existingPoints->GetPoint(cell[0], p1); unsigned int index = cellSize * 0.5; existingPoints->GetPoint(cell[index], p2); v1[0] = p2[0] - p1[0]; v1[1] = p2[1] - p1[1]; v1[2] = p2[2] - p1[2]; for (vtkIdType k = 2; k < cellSize; k++) { index = cellSize * 0.25; existingPoints->GetPoint(cell[index], p1); index = cellSize * 0.75; existingPoints->GetPoint(cell[index], p2); v2[0] = p2[0] - p1[0]; v2[1] = p2[1] - p1[1]; v2[2] = p2[2] - p1[2]; vtkMath::Cross(v1, v2, polygonNormal); if (vtkMath::Norm(polygonNormal) != 0) break; } vtkMath::Normalize(polygonNormal); // Now we start computing the normal for each vertex double vertexNormalTemp[3]; existingPoints->GetPoint(cell[0], p1); existingPoints->GetPoint(cell[1], p2); v1[0] = p2[0] - p1[0]; v1[1] = p2[1] - p1[1]; v1[2] = p2[2] - p1[2]; vtkMath::Cross(v1, polygonNormal, vertexNormalTemp); vtkMath::Normalize(vertexNormalTemp); double vertexNormal[3]; for (vtkIdType j = 0; j < cellSize - 2; j++) { existingPoints->GetPoint(cell[j + 1], p1); existingPoints->GetPoint(cell[j + 2], p2); v1[0] = p2[0] - p1[0]; v1[1] = p2[1] - p1[1]; v1[2] = p2[2] - p1[2]; vtkMath::Cross(v1, polygonNormal, vertexNormal); vtkMath::Normalize(vertexNormal); double finalNormal[3]; finalNormal[0] = (vertexNormal[0] + vertexNormalTemp[0]) * 0.5; finalNormal[1] = (vertexNormal[1] + vertexNormalTemp[1]) * 0.5; finalNormal[2] = (vertexNormal[2] + vertexNormalTemp[2]) * 0.5; vtkMath::Normalize(finalNormal); // Here we determine the direction of the normal if (m_SegmentationBinaryImage) { Point3D worldCoord; worldCoord[0] = p1[0] + finalNormal[0] * m_MaxSpacing; worldCoord[1] = p1[1] + finalNormal[1] * m_MaxSpacing; worldCoord[2] = p1[2] + finalNormal[2] * m_MaxSpacing; double val = 0.0; itk::Index<3> idx; m_SegmentationBinaryImage->GetGeometry()->WorldToIndex(worldCoord, idx); try { if (m_SegmentationBinaryImage->GetImageDescriptor() ->GetChannelDescriptor() .GetPixelType() .GetComponentType() == itk::ImageIOBase::UCHAR) { mitk::ImagePixelReadAccessor readAccess(m_SegmentationBinaryImage); val = readAccess.GetPixelByIndexSafe(idx); } else if (m_SegmentationBinaryImage->GetImageDescriptor() ->GetChannelDescriptor() .GetPixelType() .GetComponentType() == itk::ImageIOBase::USHORT) { mitk::ImagePixelReadAccessor readAccess(m_SegmentationBinaryImage); val = readAccess.GetPixelByIndexSafe(idx); } } catch (const mitk::Exception &e) { // If value is outside the image's region ignore it MITK_WARN << e.what(); } if (val == 0.0) { // MITK_INFO << "val equals zero."; ++m_PositiveNormalCounter; } else { // MITK_INFO << "val does not equal zero."; ++m_NegativeNormalCounter; } } vertexNormalTemp[0] = vertexNormal[0]; vertexNormalTemp[1] = vertexNormal[1]; vertexNormalTemp[2] = vertexNormal[2]; vtkIdType id = cell[j + 1]; normals->SetTuple(id, finalNormal); } existingPoints->GetPoint(cell[0], p1); existingPoints->GetPoint(cell[1], p2); v1[0] = p2[0] - p1[0]; v1[1] = p2[1] - p1[1]; v1[2] = p2[2] - p1[2]; vtkMath::Cross(v1, polygonNormal, vertexNormal); vtkMath::Normalize(vertexNormal); vertexNormal[0] = (vertexNormal[0] + vertexNormalTemp[0]) * 0.5; vertexNormal[1] = (vertexNormal[1] + vertexNormalTemp[1]) * 0.5; vertexNormal[2] = (vertexNormal[2] + vertexNormalTemp[2]) * 0.5; vtkMath::Normalize(vertexNormal); vtkIdType id = cell[0]; normals->SetTuple(id, vertexNormal); id = cell[cellSize - 1]; normals->SetTuple(id, vertexNormal); if (m_NegativeNormalCounter > m_PositiveNormalCounter) { for (vtkIdType n = 0; n < cellSize; n++) { double normal[3]; normals->GetTuple(offSet + n, normal); normal[0] = (-1) * normal[0]; normal[1] = (-1) * normal[1]; normal[2] = (-1) * normal[2]; normals->SetTuple(offSet + n, normal); } } m_NegativeNormalCounter = 0; m_PositiveNormalCounter = 0; offSet += cellSize; } // end for all cells Surface::Pointer surface = this->GetOutput(i); surface->GetVtkPolyData()->GetCellData()->SetNormals(normals); } // end for all inputs // Setting progressbar if (this->m_UseProgressBar) mitk::ProgressBar::GetInstance()->Progress(this->m_ProgressStepSize); } mitk::Surface::Pointer mitk::ComputeContourSetNormalsFilter::GetNormalsAsSurface() { // Just for debugging: vtkSmartPointer newPolyData = vtkSmartPointer::New(); vtkSmartPointer newLines = vtkSmartPointer::New(); vtkSmartPointer newPoints = vtkSmartPointer::New(); unsigned int idCounter(0); // Debug end for (unsigned int i = 0; i < this->GetNumberOfIndexedOutputs(); i++) { auto *currentSurface = this->GetOutput(i); vtkPolyData *polyData = currentSurface->GetVtkPolyData(); vtkSmartPointer currentCellNormals = vtkDoubleArray::SafeDownCast(polyData->GetCellData()->GetNormals()); vtkSmartPointer existingPolys = polyData->GetPolys(); vtkSmartPointer existingPoints = polyData->GetPoints(); existingPolys->InitTraversal(); - vtkIdType *cell(nullptr); + const vtkIdType *cell(nullptr); vtkIdType cellSize(0); for (existingPolys->InitTraversal(); existingPolys->GetNextCell(cellSize, cell);) { for (vtkIdType j = 0; j < cellSize; j++) { double currentNormal[3]; currentCellNormals->GetTuple(cell[j], currentNormal); vtkSmartPointer line = vtkSmartPointer::New(); line->GetPointIds()->SetNumberOfIds(2); double newPoint[3]; double p0[3]; existingPoints->GetPoint(cell[j], p0); newPoint[0] = p0[0] + currentNormal[0]; newPoint[1] = p0[1] + currentNormal[1]; newPoint[2] = p0[2] + currentNormal[2]; line->GetPointIds()->SetId(0, idCounter); newPoints->InsertPoint(idCounter, p0); idCounter++; line->GetPointIds()->SetId(1, idCounter); newPoints->InsertPoint(idCounter, newPoint); idCounter++; newLines->InsertNextCell(line); } // end for all points } // end for all cells } // end for all outputs newPolyData->SetPoints(newPoints); newPolyData->SetLines(newLines); newPolyData->BuildCells(); mitk::Surface::Pointer surface = mitk::Surface::New(); surface->SetVtkPolyData(newPolyData); return surface; } void mitk::ComputeContourSetNormalsFilter::SetMaxSpacing(double maxSpacing) { m_MaxSpacing = maxSpacing; } void mitk::ComputeContourSetNormalsFilter::GenerateOutputInformation() { Superclass::GenerateOutputInformation(); } void mitk::ComputeContourSetNormalsFilter::Reset() { for (unsigned int i = 0; i < this->GetNumberOfIndexedInputs(); i++) { this->PopBackInput(); } this->SetNumberOfIndexedInputs(0); this->SetNumberOfIndexedOutputs(0); mitk::Surface::Pointer output = mitk::Surface::New(); this->SetNthOutput(0, output.GetPointer()); } void mitk::ComputeContourSetNormalsFilter::SetUseProgressBar(bool status) { this->m_UseProgressBar = status; } void mitk::ComputeContourSetNormalsFilter::SetProgressStepSize(unsigned int stepSize) { this->m_ProgressStepSize = stepSize; } diff --git a/Modules/SurfaceInterpolation/mitkCreateDistanceImageFromSurfaceFilter.cpp b/Modules/SurfaceInterpolation/mitkCreateDistanceImageFromSurfaceFilter.cpp index 76bd92b724..87de1560a1 100644 --- a/Modules/SurfaceInterpolation/mitkCreateDistanceImageFromSurfaceFilter.cpp +++ b/Modules/SurfaceInterpolation/mitkCreateDistanceImageFromSurfaceFilter.cpp @@ -1,617 +1,617 @@ /*============================================================================ 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 "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); + const 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/mitkReduceContourSetFilter.cpp b/Modules/SurfaceInterpolation/mitkReduceContourSetFilter.cpp index d8b3ccf9f3..23d47270e0 100644 --- a/Modules/SurfaceInterpolation/mitkReduceContourSetFilter.cpp +++ b/Modules/SurfaceInterpolation/mitkReduceContourSetFilter.cpp @@ -1,519 +1,519 @@ /*============================================================================ 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 "mitkReduceContourSetFilter.h" mitk::ReduceContourSetFilter::ReduceContourSetFilter() { m_MaxSegmentLenght = 0; m_StepSize = 10; m_Tolerance = -1; m_ReductionType = DOUGLAS_PEUCKER; m_MaxSpacing = -1; m_MinSpacing = -1; this->m_UseProgressBar = false; this->m_ProgressStepSize = 1; m_NumberOfPointsAfterReduction = 0; mitk::Surface::Pointer output = mitk::Surface::New(); this->SetNthOutput(0, output.GetPointer()); } mitk::ReduceContourSetFilter::~ReduceContourSetFilter() { } void mitk::ReduceContourSetFilter::SetInput(unsigned int idx, const mitk::Surface *surface) { this->SetNthInput(idx, const_cast(surface)); this->Modified(); } void mitk::ReduceContourSetFilter::SetInput(const mitk::Surface *surface) { this->SetInput(0, surface); } void mitk::ReduceContourSetFilter::GenerateData() { unsigned int numberOfInputs = this->GetNumberOfIndexedInputs(); unsigned int numberOfOutputs(0); vtkSmartPointer newPolyData; vtkSmartPointer newPolygons; vtkSmartPointer newPoints; // For the purpose of evaluation // unsigned int numberOfPointsBefore (0); m_NumberOfPointsAfterReduction = 0; for (unsigned int i = 0; i < numberOfInputs; i++) { auto *currentSurface = this->GetInput(i); vtkSmartPointer polyData = currentSurface->GetVtkPolyData(); newPolyData = vtkSmartPointer::New(); newPolygons = vtkSmartPointer::New(); newPoints = vtkSmartPointer::New(); vtkSmartPointer existingPolys = polyData->GetPolys(); vtkSmartPointer existingPoints = polyData->GetPoints(); existingPolys->InitTraversal(); - vtkIdType *cell(nullptr); + const vtkIdType *cell(nullptr); vtkIdType cellSize(0); for (existingPolys->InitTraversal(); existingPolys->GetNextCell(cellSize, cell);) { bool incorporatePolygon = this->CheckForIntersection(cell, cellSize, existingPoints, /*numberOfIntersections, intersectionPoints, */ i); if (!incorporatePolygon) continue; vtkSmartPointer newPolygon = vtkSmartPointer::New(); if (m_ReductionType == NTH_POINT) { this->ReduceNumberOfPointsByNthPoint(cellSize, cell, existingPoints, newPolygon, newPoints); if (newPolygon->GetPointIds()->GetNumberOfIds() != 0) { newPolygons->InsertNextCell(newPolygon); } } else if (m_ReductionType == DOUGLAS_PEUCKER) { this->ReduceNumberOfPointsByDouglasPeucker(cellSize, cell, existingPoints, newPolygon, newPoints); if (newPolygon->GetPointIds()->GetNumberOfIds() > 3) { newPolygons->InsertNextCell(newPolygon); } } // Again for evaluation // numberOfPointsBefore += cellSize; m_NumberOfPointsAfterReduction += newPolygon->GetPointIds()->GetNumberOfIds(); } if (newPolygons->GetNumberOfCells() != 0) { newPolyData->SetPolys(newPolygons); newPolyData->SetPoints(newPoints); newPolyData->BuildLinks(); this->SetNumberOfIndexedOutputs(numberOfOutputs + 1); mitk::Surface::Pointer surface = mitk::Surface::New(); this->SetNthOutput(numberOfOutputs, surface.GetPointer()); surface->SetVtkPolyData(newPolyData); numberOfOutputs++; } } // MITK_INFO<<"Points before: "<SetNumberOfIndexedOutputs(numberOfOutputs); if (numberOfOutputs == 0) { mitk::Surface::Pointer tmp_output = mitk::Surface::New(); tmp_output->SetVtkPolyData(vtkPolyData::New()); this->SetNthOutput(0, tmp_output.GetPointer()); } // Setting progressbar if (this->m_UseProgressBar) mitk::ProgressBar::GetInstance()->Progress(this->m_ProgressStepSize); } void mitk::ReduceContourSetFilter::ReduceNumberOfPointsByNthPoint( - vtkIdType cellSize, vtkIdType *cell, vtkPoints *points, vtkPolygon *reducedPolygon, vtkPoints *reducedPoints) + vtkIdType cellSize, const vtkIdType *cell, vtkPoints *points, vtkPolygon *reducedPolygon, vtkPoints *reducedPoints) { unsigned int newNumberOfPoints(0); unsigned int mod = cellSize % m_StepSize; if (mod == 0) { newNumberOfPoints = cellSize / m_StepSize; } else { newNumberOfPoints = ((cellSize - mod) / m_StepSize) + 1; } if (newNumberOfPoints <= 3) { return; } reducedPolygon->GetPointIds()->SetNumberOfIds(newNumberOfPoints); reducedPolygon->GetPoints()->SetNumberOfPoints(newNumberOfPoints); for (vtkIdType i = 0; i < cellSize; i++) { if (i % m_StepSize == 0) { double point[3]; points->GetPoint(cell[i], point); vtkIdType id = reducedPoints->InsertNextPoint(point); reducedPolygon->GetPointIds()->SetId(i / m_StepSize, id); } } vtkIdType id = cell[0]; double point[3]; points->GetPoint(id, point); id = reducedPoints->InsertNextPoint(point); reducedPolygon->GetPointIds()->SetId(newNumberOfPoints - 1, id); } void mitk::ReduceContourSetFilter::ReduceNumberOfPointsByDouglasPeucker( - vtkIdType cellSize, vtkIdType *cell, vtkPoints *points, vtkPolygon *reducedPolygon, vtkPoints *reducedPoints) + vtkIdType cellSize, const vtkIdType *cell, vtkPoints *points, vtkPolygon *reducedPolygon, vtkPoints *reducedPoints) { // If the cell is too small to obtain a reduced polygon with the given stepsize return if (cellSize <= static_cast(m_StepSize * 3)) return; /* What we do now is (see the Douglas Peucker Algorithm): 1. Divide the current contour in two line segments (start - middle; middle - end), put them into the stack 2. Fetch first line segment and create the following vectors: - v1 = (start;end) - v2 = (start;currentPoint) -> for each point of the current line segment! 3. Calculate the distance from the currentPoint to v1: a. Determine the length of the orthogonal projection of v2 to v1 by: l = v2 * (normalized v1) b. There a three possibilities for the distance then: d = sqrt(lenght(v2)^2 - l^2) if l > 0 and l < length(v1) d = lenght(v2-v1) if l > 0 and l > lenght(v1) d = length(v2) if l < 0 because v2 is then pointing in a different direction than v1 4. Memorize the point with the biggest distance and create two new line segments with it at the end of the iteration and put it into the stack 5. If the distance value D <= m_Tolerance, then add the start and end index and the corresponding points to the reduced ones */ // First of all set tolerance if none is specified if (m_Tolerance < 0) { if (m_MaxSpacing > 0) { m_Tolerance = m_MinSpacing; } else { m_Tolerance = 1.5; } } std::stack lineSegments; // 1. Divide in line segments LineSegment ls2; ls2.StartIndex = cell[cellSize / 2]; ls2.EndIndex = cell[cellSize - 1]; lineSegments.push(ls2); LineSegment ls1; ls1.StartIndex = cell[0]; ls1.EndIndex = cell[cellSize / 2]; lineSegments.push(ls1); LineSegment currentSegment; double v1[3]; double v2[3]; double tempV[3]; double lenghtV1; double currentMaxDistance(0); vtkIdType currentMaxDistanceIndex(0); double l; double d; vtkIdType pointId(0); // Add the start index to the reduced points. From now on just the end indices will be added pointId = reducedPoints->InsertNextPoint(points->GetPoint(cell[0])); reducedPolygon->GetPointIds()->InsertNextId(pointId); while (!lineSegments.empty()) { currentSegment = lineSegments.top(); lineSegments.pop(); // 2. Create vectors points->GetPoint(currentSegment.EndIndex, tempV); points->GetPoint(currentSegment.StartIndex, v1); v1[0] = tempV[0] - v1[0]; v1[1] = tempV[1] - v1[1]; v1[2] = tempV[2] - v1[2]; lenghtV1 = vtkMath::Norm(v1); vtkMath::Normalize(v1); int range = currentSegment.EndIndex - currentSegment.StartIndex; for (int i = 1; i < abs(range); ++i) { points->GetPoint(currentSegment.StartIndex + i, tempV); points->GetPoint(currentSegment.StartIndex, v2); v2[0] = tempV[0] - v2[0]; v2[1] = tempV[1] - v2[1]; v2[2] = tempV[2] - v2[2]; // 3. Calculate the distance l = vtkMath::Dot(v2, v1); d = vtkMath::Norm(v2); if (l > 0 && l < lenghtV1) { d = sqrt((d * d - l * l)); } else if (l > 0 && l > lenghtV1) { tempV[0] = lenghtV1 * v1[0] - v2[0]; tempV[1] = lenghtV1 * v1[1] - v2[1]; tempV[2] = lenghtV1 * v1[2] - v2[2]; d = vtkMath::Norm(tempV); } // 4. Memorize maximum distance if (d > currentMaxDistance) { currentMaxDistance = d; currentMaxDistanceIndex = currentSegment.StartIndex + i; } } // 4. & 5. if (currentMaxDistance <= m_Tolerance) { // double temp[3]; int segmentLenght = currentSegment.EndIndex - currentSegment.StartIndex; if (segmentLenght > (int)m_MaxSegmentLenght) { m_MaxSegmentLenght = (unsigned int)segmentLenght; } // MITK_INFO<<"Lenght: "< 25) { unsigned int newLenght(segmentLenght); while (newLenght > 25) { newLenght = newLenght * 0.5; } unsigned int divisions = abs(segmentLenght) / newLenght; // MITK_INFO<<"Divisions: "<InsertNextPoint(points->GetPoint(currentSegment.StartIndex + newLenght * i)); reducedPolygon->GetPointIds()->InsertNextId(pointId); } } // MITK_INFO<<"Inserting END: "<InsertNextPoint(points->GetPoint(currentSegment.EndIndex)); reducedPolygon->GetPointIds()->InsertNextId(pointId); } else { ls2.StartIndex = currentMaxDistanceIndex; ls2.EndIndex = currentSegment.EndIndex; lineSegments.push(ls2); ls1.StartIndex = currentSegment.StartIndex; ls1.EndIndex = currentMaxDistanceIndex; lineSegments.push(ls1); } currentMaxDistance = 0; } } bool mitk::ReduceContourSetFilter::CheckForIntersection( - vtkIdType *currentCell, + const vtkIdType *currentCell, vtkIdType currentCellSize, vtkPoints *currentPoints, /* vtkIdType numberOfIntersections, vtkIdType* intersectionPoints,*/ unsigned int currentInputIndex) { /* If we check the current cell for intersections then we have to consider three possibilies: 1. There is another cell among all the other input surfaces which intersects the current polygon: - That means we have to save the intersection points because these points should not be eliminated 2. There current polygon exists just because of an intersection of another polygon with the current plane defined by the current polygon - That means the current polygon should not be incorporated and all of its points should be eliminated 3. There is no intersection - That mean we can just reduce the current polygons points without considering any intersections */ for (unsigned int i = 0; i < this->GetNumberOfIndexedInputs(); i++) { // Don't check for intersection with the polygon itself if (i == currentInputIndex) continue; // Get the next polydata to check for intersection vtkSmartPointer poly = this->GetInput(i)->GetVtkPolyData(); vtkSmartPointer polygonArray = poly->GetPolys(); polygonArray->InitTraversal(); vtkIdType anotherInputPolygonSize(0); - vtkIdType *anotherInputPolygonIDs(nullptr); + const vtkIdType *anotherInputPolygonIDs(nullptr); /* The procedure is: - Create the equation of the plane, defined by the points of next input - Calculate the distance of each point of the current polygon to the plane - If the maximum distance is not bigger than 1.5 of the maximum spacing AND the minimal distance is not bigger than 0.5 of the minimum spacing then the current contour is an intersection contour */ for (polygonArray->InitTraversal(); polygonArray->GetNextCell(anotherInputPolygonSize, anotherInputPolygonIDs);) { // Choosing three plane points to calculate the plane vectors double p1[3]; double p2[3]; double p3[3]; // The plane vectors double v1[3]; double v2[3] = {0}; // The plane normal double normal[3]; // Create first Vector poly->GetPoint(anotherInputPolygonIDs[0], p1); poly->GetPoint(anotherInputPolygonIDs[1], p2); v1[0] = p2[0] - p1[0]; v1[1] = p2[1] - p1[1]; v1[2] = p2[2] - p1[2]; // Find 3rd point for 2nd vector (The angle between the two plane vectors should be bigger than 30 degrees) double maxDistance(0); double minDistance(10000); for (vtkIdType j = 2; j < anotherInputPolygonSize; j++) { poly->GetPoint(anotherInputPolygonIDs[j], p3); v2[0] = p3[0] - p1[0]; v2[1] = p3[1] - p1[1]; v2[2] = p3[2] - p1[2]; // Calculate the angle between the two vector for the current point double dotV1V2 = vtkMath::Dot(v1, v2); double absV1 = sqrt(vtkMath::Dot(v1, v1)); double absV2 = sqrt(vtkMath::Dot(v2, v2)); double cosV1V2 = dotV1V2 / (absV1 * absV2); double arccos = acos(cosV1V2); double degree = vtkMath::DegreesFromRadians(arccos); // If angle is bigger than 30 degrees break if (degree > 30) break; } // for (to find 3rd point) // Calculate normal of the plane by taking the cross product of the two vectors vtkMath::Cross(v1, v2, normal); vtkMath::Normalize(normal); // Determine position of the plane double lambda = vtkMath::Dot(normal, p1); /* Calculate the distance to the plane for each point of the current polygon If the distance is zero then save the currentPoint as intersection point */ for (vtkIdType k = 0; k < currentCellSize; k++) { double currentPoint[3]; currentPoints->GetPoint(currentCell[k], currentPoint); double tempPoint[3]; tempPoint[0] = normal[0] * currentPoint[0]; tempPoint[1] = normal[1] * currentPoint[1]; tempPoint[2] = normal[2] * currentPoint[2]; double temp = tempPoint[0] + tempPoint[1] + tempPoint[2] - lambda; double distance = fabs(temp); if (distance > maxDistance) { maxDistance = distance; } if (distance < minDistance) { minDistance = distance; } } // for (to calculate distance and intersections with currentPolygon) if (maxDistance < 1.5 * m_MaxSpacing && minDistance < 0.5 * m_MinSpacing) { return false; } // Because we are considering the plane defined by the acual input polygon only one iteration is sufficient // We do not need to consider each cell of the plane break; } // for (to traverse through all cells of actualInputPolyData) } // for (to iterate through all inputs) return true; } void mitk::ReduceContourSetFilter::GenerateOutputInformation() { Superclass::GenerateOutputInformation(); } void mitk::ReduceContourSetFilter::Reset() { for (unsigned int i = 0; i < this->GetNumberOfIndexedInputs(); i++) { this->PopBackInput(); } this->SetNumberOfIndexedInputs(0); this->SetNumberOfIndexedOutputs(0); // BUG XXXXX Fix mitk::Surface::Pointer output = mitk::Surface::New(); this->SetNthOutput(0, output.GetPointer()); m_NumberOfPointsAfterReduction = 0; } void mitk::ReduceContourSetFilter::SetUseProgressBar(bool status) { this->m_UseProgressBar = status; } void mitk::ReduceContourSetFilter::SetProgressStepSize(unsigned int stepSize) { this->m_ProgressStepSize = stepSize; } diff --git a/Modules/SurfaceInterpolation/mitkReduceContourSetFilter.h b/Modules/SurfaceInterpolation/mitkReduceContourSetFilter.h index 1f18b80580..89c447afe4 100644 --- a/Modules/SurfaceInterpolation/mitkReduceContourSetFilter.h +++ b/Modules/SurfaceInterpolation/mitkReduceContourSetFilter.h @@ -1,134 +1,134 @@ /*============================================================================ 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 mitkReduceContourSetFilter_h_Included #define mitkReduceContourSetFilter_h_Included #include "mitkProgressBar.h" #include "mitkSurface.h" #include "mitkSurfaceToSurfaceFilter.h" #include #include "vtkCellArray.h" #include "vtkMath.h" #include "vtkPoints.h" #include "vtkPolyData.h" #include "vtkPolygon.h" #include "vtkSmartPointer.h" #include namespace mitk { /** \brief A filter that reduces the number of points of contours represented by a mitk::Surface The type of the reduction can be set via SetReductionType. The two ways provided by this filter is the \li NTH_POINT Algorithm which reduces the contours according to a certain stepsize \li DOUGLAS_PEUCKER Algorithm which incorpates an error tolerance into the reduction. Stepsize and error tolerance can be set via SetStepSize and SetTolerance. Additional if more than one input contour is provided this filter tries detect contours which occur just because of an intersection. These intersection contours are eliminated. In oder to ensure a correct elimination the min and max spacing of the original image must be provided. The output is a mitk::Surface. $Author: fetzer$ */ class MITKSURFACEINTERPOLATION_EXPORT ReduceContourSetFilter : public SurfaceToSurfaceFilter { public: enum Reduction_Type { NTH_POINT, DOUGLAS_PEUCKER }; struct LineSegment { unsigned int StartIndex; unsigned int EndIndex; }; mitkClassMacro(ReduceContourSetFilter, SurfaceToSurfaceFilter); itkFactorylessNewMacro(Self); itkCloneMacro(Self); itkSetMacro(MinSpacing, double); itkSetMacro(MaxSpacing, double); itkSetMacro(ReductionType, Reduction_Type); itkSetMacro(StepSize, unsigned int); itkSetMacro(Tolerance, double); itkGetMacro(NumberOfPointsAfterReduction, unsigned int); // Resets the filter, i.e. removes all inputs and outputs void Reset(); /** \brief Set whether the mitkProgressBar should be used \a Parameter true for using the progress bar, false otherwise */ void SetUseProgressBar(bool); using itk::ProcessObject::SetInput; void SetInput(const mitk::Surface *surface) override; void SetInput(unsigned int idx, const mitk::Surface *surface) override; /** \brief Set the stepsize which the progress bar should proceed \a Parameter The stepsize for progressing */ void SetProgressStepSize(unsigned int stepSize); protected: ReduceContourSetFilter(); ~ReduceContourSetFilter() override; void GenerateData() override; void GenerateOutputInformation() override; private: void ReduceNumberOfPointsByNthPoint( - vtkIdType cellSize, vtkIdType *cell, vtkPoints *points, vtkPolygon *reducedPolygon, vtkPoints *reducedPoints); + vtkIdType cellSize, const vtkIdType *cell, vtkPoints *points, vtkPolygon *reducedPolygon, vtkPoints *reducedPoints); void ReduceNumberOfPointsByDouglasPeucker( - vtkIdType cellSize, vtkIdType *cell, vtkPoints *points, vtkPolygon *reducedPolygon, vtkPoints *reducedPoints); + vtkIdType cellSize, const vtkIdType *cell, vtkPoints *points, vtkPolygon *reducedPolygon, vtkPoints *reducedPoints); bool CheckForIntersection( - vtkIdType *currentCell, + const vtkIdType *currentCell, vtkIdType currentCellSize, vtkPoints *currentPoints, /*vtkIdType numberOfIntersections, vtkIdType* intersectionPoints,*/ unsigned int currentInputIndex); double m_MinSpacing; double m_MaxSpacing; Reduction_Type m_ReductionType; unsigned int m_StepSize; double m_Tolerance; unsigned int m_MaxSegmentLenght; bool m_UseProgressBar; unsigned int m_ProgressStepSize; unsigned int m_NumberOfPointsAfterReduction; }; // class } // namespace #endif