diff --git a/Modules/SurfaceInterpolation/mitkComputeContourSetNormalsFilter.cpp b/Modules/SurfaceInterpolation/mitkComputeContourSetNormalsFilter.cpp index 6f62e34e1d..fc51c0d25f 100644 --- a/Modules/SurfaceInterpolation/mitkComputeContourSetNormalsFilter.cpp +++ b/Modules/SurfaceInterpolation/mitkComputeContourSetNormalsFilter.cpp @@ -1,343 +1,342 @@ /*=================================================================== 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 "mitkComputeContourSetNormalsFilter.h" #include "mitkImagePixelReadAccessor.h" #include "mitkIOUtil.h" mitk::ComputeContourSetNormalsFilter::ComputeContourSetNormalsFilter() : m_SegmentationBinaryImage(NULL) , 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(); this->CreateOutputsForAllInputs(numberOfInputs); //Iterating over each input for(unsigned int i = 0; i < numberOfInputs; i++) { //Getting the inputs polydata and polygons Surface* currentSurface = const_cast( this->GetInput(i) ); vtkPolyData* polyData = currentSurface->GetVtkPolyData(); vtkSmartPointer existingPolys = polyData->GetPolys(); vtkSmartPointer existingPoints = polyData->GetPoints(); existingPolys->InitTraversal(); vtkIdType* cell (NULL); 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]; - mitk::IOUtil::SaveImage( m_SegmentationBinaryImage, "C:/tmp/interpol/SegmentationBinaryImage.nrrd" ); 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; //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; std::cout << "World coords: " << worldCoord[0] << " , " << worldCoord[1] << " , " << worldCoord[2] << std::endl; std::cout << "...for point: " << p1[0] << " , " << p1[1] << " , " << p1[2] << std::endl; std::cout << "...and normal: " << finalNormal[0] << " , " << finalNormal[1] << " , " << finalNormal[2] << std::endl; std::cout << "MaxSpacing: " << m_MaxSpacing << std::endl; double val = 0.0; mitk::ImagePixelReadAccessor readAccess(m_SegmentationBinaryImage); itk::Index<3> idx; m_SegmentationBinaryImage->GetGeometry()->WorldToIndex(worldCoord, idx); try { val = readAccess.GetPixelByIndexSafe(idx); } catch (mitk::Exception e) { // If value is outside the image's region ignore it } 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; vtkIdType id = cell[0]; normals->SetTuple(id,vertexNormal); id = cell[cellSize-1]; normals->SetTuple(id,vertexNormal); std::cout << "Negative normal counter: " << m_NegativeNormalCounter << std::endl; std::cout << "Positive normal counter: " << m_PositiveNormalCounter << std::endl; - int n = 5; + 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++) { Surface* currentSurface = const_cast( 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 (NULL); 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/mitkSurfaceInterpolationController.h b/Modules/SurfaceInterpolation/mitkSurfaceInterpolationController.h index 06275270f5..0f9960a6fc 100644 --- a/Modules/SurfaceInterpolation/mitkSurfaceInterpolationController.h +++ b/Modules/SurfaceInterpolation/mitkSurfaceInterpolationController.h @@ -1,266 +1,266 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef mitkSurfaceInterpolationController_h_Included #define mitkSurfaceInterpolationController_h_Included #include "mitkCommon.h" #include #include "mitkRestorePlanePositionOperation.h" #include "mitkSurface.h" #include "mitkInteractionConst.h" #include "mitkColorProperty.h" #include "mitkProperties.h" #include "mitkCreateDistanceImageFromSurfaceFilter.h" #include "mitkReduceContourSetFilter.h" #include "mitkComputeContourSetNormalsFilter.h" #include "mitkDataNode.h" #include "mitkDataStorage.h" #include "mitkWeakPointer.h" #include "vtkPolygon.h" #include "vtkPoints.h" #include "vtkCellArray.h" #include "vtkPolyData.h" #include "vtkSmartPointer.h" #include "vtkAppendPolyData.h" #include "vtkMarchingCubes.h" #include "vtkImageData.h" #include "mitkVtkRepresentationProperty.h" #include "vtkProperty.h" #include "mitkImageTimeSelector.h" #include "mitkProgressBar.h" namespace mitk { class MitkSurfaceInterpolation_EXPORT SurfaceInterpolationController : public itk::Object { public: mitkClassMacro(SurfaceInterpolationController, itk::Object); itkFactorylessNewMacro(Self); itkCloneMacro(Self); struct ContourPositionInformation { Surface::Pointer contour; Vector3D contourNormal; Point3D contourPoint; }; typedef std::vector ContourPositionInformationList; typedef std::vector ContourPositionInformationVec2D; //typedef std::map ContourListMap; typedef std::map ContourListMap; static SurfaceInterpolationController* GetInstance(); - void SetCurrentTimeStep( int ts ) + void SetCurrentTimeStep( unsigned int ts ) { if ( m_CurrentTimeStep != ts ) { m_CurrentTimeStep = ts; if ( m_SelectedSegmentation ) { m_TimeSelector->SetInput( m_SelectedSegmentation ); m_TimeSelector->SetTimeNr( m_CurrentTimeStep ); m_TimeSelector->SetChannelNr( 0 ); m_TimeSelector->Update(); this->ReinitializeInterpolation(); } } }; unsigned int GetCurrentTimeStep() { return m_CurrentTimeStep; }; /** * @brief Adds a new extracted contour to the list * @param newContour the contour to be added. If a contour at that position * already exists the related contour will be updated */ void AddNewContour (Surface::Pointer newContour); /** * @brief Removes the contour for a given plane for the current selected segmenation * @param contourInfo the contour which should be removed * @return true if a contour was found and removed, false if no contour was found */ bool RemoveContour (ContourPositionInformation contourInfo ); /** * @brief Adds new extracted contours to the list. If one or more contours at a given position * already exist they will be updated respectively * @param newContours the list of the contours */ void AddNewContours (std::vector newContours); /** * @brief Returns the contour for a given plane for the current selected segmenation * @param ontourInfo the contour which should be returned * @return the contour as an mitk::Surface. If no contour is available at the give position NULL is returned */ const mitk::Surface* GetContour (ContourPositionInformation contourInfo ); /** * @brief Returns the number of available contours for the current selected segmentation * @return the number of contours */ unsigned int GetNumberOfContours(); /** * Interpolates the 3D surface from the given extracted contours */ void Interpolate (); mitk::Surface::Pointer GetInterpolationResult(); /** * Sets the minimum spacing of the current selected segmentation * This is needed since the contour points we reduced before they are used to interpolate the surface */ void SetMinSpacing(double minSpacing); /** * Sets the minimum spacing of the current selected segmentation * This is needed since the contour points we reduced before they are used to interpolate the surface */ void SetMaxSpacing(double maxSpacing); /** * Sets the volume i.e. the number of pixels that the distance image should have * By evaluation we found out that 50.000 pixel delivers a good result */ void SetDistanceImageVolume(unsigned int distImageVolume); /** * @brief Get the current selected segmentation for which the interpolation is performed * @return the current segmentation image */ mitk::Image::Pointer GetCurrentSegmentation(); Surface* GetContoursAsSurface(); void SetDataStorage(DataStorage::Pointer ds); /** * Sets the current list of contourpoints which is used for the surface interpolation * @param segmentation The current selected segmentation * \deprecatedSince{2014_03} */ DEPRECATED (void SetCurrentSegmentationInterpolationList(mitk::Image::Pointer segmentation)); /** * Sets the current list of contourpoints which is used for the surface interpolation * @param segmentation The current selected segmentation */ void SetCurrentInterpolationSession(mitk::Image::Pointer currentSegmentationImage); /** * Removes the segmentation and all its contours from the list * @param segmentation The segmentation to be removed * \deprecatedSince{2014_03} */ DEPRECATED (void RemoveSegmentationFromContourList(mitk::Image* segmentation)); /** * @brief Remove interpolation session * @param segmentationImage the session to be removed */ void RemoveInterpolationSession(mitk::Image::Pointer segmentationImage); /** * Replaces the current interpolation session with a new one. All contours form the old * session will be applied to the new session. This only works if the two images have the * geometry * @param oldSession the session which should be replaced * @param newSession the new session which replaces the old one * @return true it the the replacement was successful, false if not (e.g. the image's geometry differs) */ bool ReplaceInterpolationSession(mitk::Image::Pointer oldSession, mitk::Image::Pointer newSession); /** * @brief Removes all sessions */ void RemoveAllInterpolationSessions(); /** * @brief Reinitializes the interpolation using the provided contour data * @param contours a mitk::Surface which contains the contours as polys in the vtkPolyData */ void ReinitializeInterpolation(mitk::Surface::Pointer contours); mitk::Image* GetImage(); /** * Estimates the memory which is needed to build up the equationsystem for the interpolation. * \returns The percentage of the real memory which will be used by the interpolation */ double EstimatePortionOfNeededMemory(); unsigned int GetNumberOfInterpolationSessions(); protected: SurfaceInterpolationController(); ~SurfaceInterpolationController(); template void GetImageBase(itk::Image* input, itk::ImageBase<3>::Pointer& result); private: void OnSegmentationDeleted(const itk::Object *caller, const itk::EventObject &event); void ReinitializeInterpolation(); void AddToInterpolationPipeline(ContourPositionInformation contourInfo ); ReduceContourSetFilter::Pointer m_ReduceFilter; ComputeContourSetNormalsFilter::Pointer m_NormalsFilter; CreateDistanceImageFromSurfaceFilter::Pointer m_InterpolateSurfaceFilter; Surface::Pointer m_Contours; vtkSmartPointer m_PolyData; mitk::DataStorage::Pointer m_DataStorage; ContourListMap m_ListOfInterpolationSessions; mitk::Surface::Pointer m_InterpolationResult; unsigned int m_CurrentNumberOfReducedContours; mitk::Image* m_SelectedSegmentation; std::map m_SegmentationObserverTags; - int m_CurrentTimeStep; + unsigned int m_CurrentTimeStep; mitk::ImageTimeSelector::Pointer m_TimeSelector; }; } #endif