diff --git a/Modules/MitkExt/Algorithms/mitkComputeContourSetNormalsFilter.cpp b/Modules/MitkExt/Algorithms/mitkComputeContourSetNormalsFilter.cpp index cdf50a0cb1..0ced965fb5 100644 --- a/Modules/MitkExt/Algorithms/mitkComputeContourSetNormalsFilter.cpp +++ b/Modules/MitkExt/Algorithms/mitkComputeContourSetNormalsFilter.cpp @@ -1,287 +1,303 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Module: $RCSfile$ Language: C++ Date: $Date: $ Version: $Revision: $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #include "mitkComputeContourSetNormalsFilter.h" mitk::ComputeContourSetNormalsFilter::ComputeContourSetNormalsFilter() { m_MaxSpacing = 5; + this->m_UseProgressBar = false; + this->m_ProgressStepSize = 1; } mitk::ComputeContourSetNormalsFilter::~ComputeContourSetNormalsFilter() { } void mitk::ComputeContourSetNormalsFilter::GenerateData() { unsigned int numberOfInputs = this->GetNumberOfInputs(); 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; //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 (unsigned int 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 (unsigned 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 (j == 0 && 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 = m_SegmentationBinaryImage->GetPixelValueByWorldCoordinate(worldCoord); if (val == 1.0) { ++m_PositiveNormalCounter; } else { ++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); int normalDirection(-1); if(m_NegativeNormalCounter < m_PositiveNormalCounter) { normalDirection = 1; } for(unsigned int n = 0; n < normals->GetNumberOfTuples(); n++) { double normal[3]; normals->GetTuple(n,normal); normal[0] = normalDirection*normal[0]; normal[1] = normalDirection*normal[1]; normal[2] = normalDirection*normal[2]; } }//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 = vtkPolyData::New(); vtkSmartPointer newLines = vtkCellArray::New(); vtkSmartPointer newPoints = vtkPoints::New(); unsigned int idCounter (0); //Debug end for (unsigned int i = 0; i < this->GetNumberOfOutputs(); 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 ( unsigned int j = 0; j < cellSize; j++ ) { double currentNormal[3]; currentCellNormals->GetTuple(cell[j], currentNormal); vtkSmartPointer line = vtkLine::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->GetNumberOfInputs(); i++) { this->PopBackInput(); } this->SetNumberOfInputs(0); this->SetNumberOfOutputs(0); } + +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/MitkExt/Algorithms/mitkComputeContourSetNormalsFilter.h b/Modules/MitkExt/Algorithms/mitkComputeContourSetNormalsFilter.h index 460a75faf5..b3467bd15f 100644 --- a/Modules/MitkExt/Algorithms/mitkComputeContourSetNormalsFilter.h +++ b/Modules/MitkExt/Algorithms/mitkComputeContourSetNormalsFilter.h @@ -1,87 +1,106 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Module: $RCSfile$ Language: C++ Date: $Date: $ Version: $Revision: $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #ifndef mitkComputeContourSetNormalsFilter_h_Included #define mitkComputeContourSetNormalsFilter_h_Included #include "MitkExtExports.h" #include "mitkSurfaceToSurfaceFilter.h" +#include "mitkProgressBar.h" +#include "mitkSurface.h" + #include "vtkCellArray.h" #include "vtkPolyData.h" #include "vtkSmartPointer.h" -#include "mitkSurface.h" #include "vtkDoubleArray.h" #include "vtkMath.h" #include "vtkCellData.h" #include "vtkLine.h" #include "mitkImage.h" namespace mitk { /** \brief Filter to compute the normales for contours based on vtkPolygons This filter takes a number of extracted contours and computes the normals for each contour edge point. The normals can be accessed by calling: filter->GetOutput(i)->GetVtkPolyData()->GetCellData()->GetNormals(); See also the method GetNormalsAsSurface() Note: If a segmentation binary image is provided this filter assures that the computed normals do not point into the segmentation image $Author: fetzer$ */ class MitkExt_EXPORT ComputeContourSetNormalsFilter : public SurfaceToSurfaceFilter { public: mitkClassMacro(ComputeContourSetNormalsFilter,SurfaceToSurfaceFilter); itkNewMacro(Self); itkSetMacro(SegmentationBinaryImage, mitk::Image::Pointer); /* \brief Returns the computed normals as a surface */ mitk::Surface::Pointer GetNormalsAsSurface(); //Resets the filter, i.e. removes all inputs and outputs void Reset(); void SetMaxSpacing(double); + /** + \brief Set whether the mitkProgressBar should be used + + \a Parameter true for using the progress bar, false otherwise + */ + void SetUseProgressBar(bool); + + /** + \brief Set the stepsize which the progress bar should proceed + + \a Parameter The stepsize for progressing + */ + void SetProgressStepSize(unsigned int stepSize); + protected: ComputeContourSetNormalsFilter(); virtual ~ComputeContourSetNormalsFilter(); virtual void GenerateData(); virtual void GenerateOutputInformation(); private: //The segmentation out of which the contours were extracted. Can be used to determine the direction of the normals mitk::Image::Pointer m_SegmentationBinaryImage; double m_MaxSpacing; unsigned int m_NegativeNormalCounter; unsigned int m_PositiveNormalCounter; + bool m_UseProgressBar; + unsigned int m_ProgressStepSize; + };//class }//namespace #endif diff --git a/Modules/MitkExt/Algorithms/mitkCreateDistanceImageFromSurfaceFilter.cpp b/Modules/MitkExt/Algorithms/mitkCreateDistanceImageFromSurfaceFilter.cpp index 2f9a508af4..3bbf2851df 100644 --- a/Modules/MitkExt/Algorithms/mitkCreateDistanceImageFromSurfaceFilter.cpp +++ b/Modules/MitkExt/Algorithms/mitkCreateDistanceImageFromSurfaceFilter.cpp @@ -1,488 +1,508 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Module: $RCSfile$ Language: C++ Date: $Date: $ Version: $Revision: $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #include "mitkCreateDistanceImageFromSurfaceFilter.h" mitk::CreateDistanceImageFromSurfaceFilter::CreateDistanceImageFromSurfaceFilter() { m_DistanceImageVolume = 50000; + this->m_UseProgressBar = false; + this->m_ProgressStepSize = 5; } mitk::CreateDistanceImageFromSurfaceFilter::~CreateDistanceImageFromSurfaceFilter() { } void mitk::CreateDistanceImageFromSurfaceFilter::GenerateData() { //First of all we have to build the equation-system from the existing contour-edge-points this->CreateSolutionMatrixAndFunctionValues(); //Then we solve the equation-system via QR - decomposition. The interpolation weights are obtained in that way vnl_qr solver (m_SolutionMatrix); m_Weights = solver.solve(m_FunctionValues); + //Setting progressbar + if (this->m_UseProgressBar) + mitk::ProgressBar::GetInstance()->Progress(2); + //The last step is to create the distance map with the interpolated distance function this->CreateDistanceImage(); m_Centers.clear(); m_FunctionValues.clear(); m_Normals.clear(); m_Weights.clear(); m_SolutionMatrix.clear(); + + //Setting progressbar + if (this->m_UseProgressBar) + mitk::ProgressBar::GetInstance()->Progress(3); } void mitk::CreateDistanceImageFromSurfaceFilter::CreateSolutionMatrixAndFunctionValues() { unsigned int numberOfInputs = this->GetNumberOfInputs(); 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 Surface* currentSurface; vtkSmartPointer polyData; vtkSmartPointer currentCellNormals; vtkSmartPointer existingPolys; vtkSmartPointer existingPoints; double p[3]; PointType currentPoint; PointType normal; for (unsigned int i = 0; i < numberOfInputs; i++) { currentSurface = const_cast( 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 (NULL); vtkIdType cellSize (0); for( existingPolys->InitTraversal(); existingPolys->GetNextCell(cellSize, cell);) { for ( unsigned int 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 //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.set_size(numberOfCenters*3); m_FunctionValues.fill(0); //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]; currentPoint[1] = currentPoint[1] - normal[1]; currentPoint[2] = currentPoint[2] - normal[2]; m_Centers.push_back(currentPoint); m_FunctionValues.put(numberOfCenters+i, -1); } //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]; currentPoint[1] = currentPoint[1] + normal[1]; currentPoint[2] = currentPoint[2] + normal[2]; m_Centers.push_back(currentPoint); m_FunctionValues.put(numberOfCenters*2+i, 1); } //Now we have created all centers and all function values. Next step is to create the solution matrix numberOfCenters = m_Centers.size(); m_SolutionMatrix.set_size(numberOfCenters, numberOfCenters); m_Weights.set_size(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::CreateDistanceImage() { typedef itk::Image DistanceImageType; typedef itk::ImageRegionIteratorWithIndex ImageIterator; typedef itk::NeighborhoodIterator NeighborhoodImageIterator; DistanceImageType::Pointer distanceImg = DistanceImageType::New(); //Determin the bounding box of the delineated contours double xmin = m_Centers.at(0)[0]; double ymin = m_Centers.at(0)[1]; double zmin = m_Centers.at(0)[2]; double xmax = m_Centers.at(0)[0]; double ymax = m_Centers.at(0)[1]; double zmax = m_Centers.at(0)[2]; for (unsigned int i = 1; i < m_Centers.size(); i++) { if (xmin > m_Centers.at(i)[0]) { xmin = m_Centers.at(i)[0]; } if (ymin > m_Centers.at(i)[1]) { ymin = m_Centers.at(i)[1]; } if (zmin > m_Centers.at(i)[2]) { zmin = m_Centers.at(i)[2]; } if (xmax < m_Centers.at(i)[0]) { xmax = m_Centers.at(i)[0]; } if (ymax < m_Centers.at(i)[1]) { ymax = m_Centers.at(i)[1]; } if (zmax < m_Centers.at(i)[2]) { zmax = m_Centers.at(i)[2]; } } Vector3D extentMM; extentMM[0] = xmax - xmin + 2; extentMM[1] = ymax - ymin + 2; extentMM[2] = zmax - zmin + 2; //Shifting the distance image's offest to achieve an exact distance calculation xmin = xmin - 2; ymin = ymin - 2; zmin = zmin - 2; /* Now create an empty distance image. The create image will always have the same size, 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 = ( 500000 / extentX*extentY*extentZ )^(1/3) */ double basis = (extentMM[0]*extentMM[1]*extentMM[2]) / m_DistanceImageVolume; double exponent = 1.0/3.0; double distImgSpacing = pow(basis, exponent); int tempSpacing = (distImgSpacing+0.05)*10; m_DistanceImageSpacing = (double)tempSpacing/10.0; unsigned int numberOfXPixel = extentMM[0] / m_DistanceImageSpacing; unsigned int numberOfYPixel = extentMM[1] / m_DistanceImageSpacing; unsigned int numberOfZPixel = extentMM[2] / m_DistanceImageSpacing; DistanceImageType::SizeType size; //Increase the distance image's size a little bit to achieve an exact distance calculation size[0] = numberOfXPixel + 5; size[1] = numberOfYPixel + 5; size[2] = numberOfZPixel + 5; DistanceImageType::IndexType start; start[0] = 0; start[1] = 0; start[2] = 0; DistanceImageType::RegionType lpRegion; lpRegion.SetSize(size); lpRegion.SetIndex(start); distanceImg->SetRegions( lpRegion ); distanceImg->SetSpacing( m_DistanceImageSpacing ); distanceImg->Allocate(); //First of all the image is initialized with the value 10 for each pixel distanceImg->FillBuffer(10); /* Now we must caculate 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 start with 1. again This is done until the narrowband_point_list is empty. */ std::queue narrowbandPoints; PointType currentPoint = m_Centers.at(0); double distance = this->CalculateDistanceValue(currentPoint); DistanceImageType::IndexType currentIndex; currentIndex[0] = ( currentPoint[0]-xmin ) / m_DistanceImageSpacing; currentIndex[1] = ( currentPoint[1]-ymin ) / m_DistanceImageSpacing; currentIndex[2] = ( currentPoint[2]-zmin ) / m_DistanceImageSpacing; narrowbandPoints.push(currentIndex); distanceImg->SetPixel(currentIndex, distance); NeighborhoodImageIterator::RadiusType radius; radius.Fill(1); NeighborhoodImageIterator nIt(radius, distanceImg, distanceImg->GetLargestPossibleRegion()); unsigned int relativeNbIdx[] = {4, 10, 12, 14, 16, 22}; bool isInBounds = false; while ( !narrowbandPoints.empty() ) { nIt.SetLocation(narrowbandPoints.front()); narrowbandPoints.pop(); for (int i = 0; i < 6; i++) { nIt.GetPixel(relativeNbIdx[i], isInBounds); if( isInBounds && nIt.GetPixel(relativeNbIdx[i]) == 10) { currentIndex = nIt.GetIndex(relativeNbIdx[i]); currentPoint[0] = currentIndex[0]*m_DistanceImageSpacing + xmin; currentPoint[1] = currentIndex[1]*m_DistanceImageSpacing + ymin; currentPoint[2] = currentIndex[2]*m_DistanceImageSpacing + zmin; distance = this->CalculateDistanceValue(currentPoint); if ( abs(distance) <= m_DistanceImageSpacing*2 ) { nIt.SetPixel(relativeNbIdx[i], distance); narrowbandPoints.push(currentIndex); } } } } ImageIterator imgRegionIterator (distanceImg, distanceImg->GetLargestPossibleRegion()); imgRegionIterator.GoToBegin(); double prevPixelVal = 1; unsigned int _size[3] = { (unsigned int)(size[0] - 1), (unsigned int)(size[1] - 1), (unsigned int)(size[2] - 1) }; //Set every pixel inside the surface to -10 except the edge point (so that the received surface is closed) while (!imgRegionIterator.IsAtEnd()) { if ( imgRegionIterator.Get() == 10 && prevPixelVal < 0 ) { while (imgRegionIterator.Get() == 10) { 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(10); prevPixelVal = 10; ++imgRegionIterator; break; } else { imgRegionIterator.Set(-10); ++imgRegionIterator; prevPixelVal = -10; } } } 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(10); prevPixelVal = 10; ++imgRegionIterator; } else { prevPixelVal = imgRegionIterator.Get(); ++imgRegionIterator; } } Image::Pointer resultImage = this->GetOutput(); Point3D origin; origin[0] = xmin; origin[1] = ymin; origin[2] = zmin; CastToMitkImage(distanceImg, resultImage); resultImage->GetGeometry()->SetOrigin(origin); resultImage->SetOrigin(origin); } double mitk::CreateDistanceImageFromSurfaceFilter::CalculateDistanceValue(PointType p) { double distanceValue (0); PointType p1; PointType p2; double norm; for (unsigned int i = 0; i < m_Centers.size(); i++) { p1 = m_Centers.at(i); p2 = p-p1; norm = p2.two_norm(); distanceValue = distanceValue + norm*m_Weights.get(i); } return distanceValue; } void mitk::CreateDistanceImageFromSurfaceFilter::GenerateOutputInformation() { } void mitk::CreateDistanceImageFromSurfaceFilter::PrintEquationSystem() { std::ofstream esfile; esfile.open("C:/Users/fetzer/Desktop/equationSystem/es.txt"); esfile<<"Nummber of rows: "<SetInput( 0, const_cast( 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->GetNumberOfInputs() < 1) return NULL; return static_cast(this->ProcessObject::GetInput(0)); } const mitk::Surface* mitk::CreateDistanceImageFromSurfaceFilter::GetInput( unsigned int idx) { if (this->GetNumberOfInputs() < 1) return NULL; return static_cast(this->ProcessObject::GetInput(idx)); } void mitk::CreateDistanceImageFromSurfaceFilter::RemoveInputs(mitk::Surface* input) { this->RemoveInput(input); } void mitk::CreateDistanceImageFromSurfaceFilter::Reset() { for (unsigned int i = 0; i < this->GetNumberOfInputs(); i++) { this->PopBackInput(); } this->SetNumberOfInputs(0); this->SetNumberOfOutputs(1); } + +void mitk::CreateDistanceImageFromSurfaceFilter::SetUseProgressBar(bool status) +{ + this->m_UseProgressBar = status; +} + +void mitk::CreateDistanceImageFromSurfaceFilter::SetProgressStepSize(unsigned int stepSize) +{ + this->m_ProgressStepSize = stepSize; +} diff --git a/Modules/MitkExt/Algorithms/mitkCreateDistanceImageFromSurfaceFilter.h b/Modules/MitkExt/Algorithms/mitkCreateDistanceImageFromSurfaceFilter.h index 56c9b789b9..3113cd11b9 100644 --- a/Modules/MitkExt/Algorithms/mitkCreateDistanceImageFromSurfaceFilter.h +++ b/Modules/MitkExt/Algorithms/mitkCreateDistanceImageFromSurfaceFilter.h @@ -1,142 +1,159 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Module: $RCSfile$ Language: C++ Date: $Date: $ Version: $Revision: $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #ifndef mitkCreateDistanceImageFromSurfaceFilter_h_Included #define mitkCreateDistanceImageFromSurfaceFilter_h_Included #include "MitkExtExports.h" #include "mitkImageSource.h" #include "mitkSurface.h" +#include "mitkProgressBar.h" #include "vtkSmartPointer.h" #include "vtkDoubleArray.h" #include "vtkCellArray.h" #include "vtkCellData.h" #include "vtkPolyData.h" #include "vnl/vnl_matrix.h" #include "vnl/vnl_vector.h" #include "vnl/vnl_vector_fixed.h" #include "vnl/algo/vnl_qr.h" #include "itkImage.h" #include "itkImageRegionIteratorWithIndex.h" #include "itkNeighborhoodIterator.h" #include namespace mitk { /** \brief This filter interpolates the 3D surface for a segmented area. The basis for the interpolation are the edge-points of contours that are drawn into an image. The interpolation itself is performed via Radial Basis Function Interpolation. ATTENTION: This filter needs beside the edge points of the delineated contours additionally the normals for each edge point. \sa mitkSurfaceInterpolationController Based on the contour edge points and their normal this filter calculates a distance function with the following properties: - Putting a point into the distance function that lies inside the considered surface gives a negativ scalar value - Putting a point into the distance function that lies outside the considered surface gives a positive scalar value - Putting a point into the distance function that lies exactly on the considered surface gives the value zero With this interpolated distance function a distance image will be created. The desired surface can then be extract e.g. with the marching cubes algorithm. (Within the distance image the surface goes exactly where the pixelvalues are zero) Note that the obtained distance image has always an isotropig spacing. The size (in this case volume) of the image can be adjusted by calling SetDistanceImageVolume(unsigned int volume) which specifies the number ob pixels enclosed by the image. \ingroup Process $Author: fetzer$ */ class MitkExt_EXPORT CreateDistanceImageFromSurfaceFilter : public ImageSource { public: typedef vnl_vector_fixed PointType; typedef std::vector< PointType > NormalList; typedef std::vector< PointType > CenterList; typedef vnl_matrix SolutionMatrix; typedef vnl_vector FunctionValues; typedef vnl_vector InterpolationWeights; typedef std::vector SurfaceList; mitkClassMacro(CreateDistanceImageFromSurfaceFilter,ImageSource); itkNewMacro(Self); //Methods copied from mitkSurfaceToSurfaceFilter virtual void SetInput( const mitk::Surface* surface ); virtual void SetInput( unsigned int idx, const mitk::Surface* surface ); virtual const mitk::Surface* GetInput(); virtual const mitk::Surface* GetInput( unsigned int idx ); virtual void RemoveInputs(mitk::Surface* input); /* \brief Set the size of the output distance image. The size is specified by the image's volume (i.e. in this case how many pixels are enclosed by the image) If non is set, the volume will be 500000 pixels. */ itkSetMacro(DistanceImageVolume, unsigned int); void PrintEquationSystem(); //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); + + /** + \brief Set the stepsize which the progress bar should proceed + + \a Parameter The stepsize for progressing + */ + void SetProgressStepSize(unsigned int stepSize); + protected: CreateDistanceImageFromSurfaceFilter(); virtual ~CreateDistanceImageFromSurfaceFilter(); virtual void GenerateData(); virtual void GenerateOutputInformation(); private: void CreateSolutionMatrixAndFunctionValues(); double CalculateDistanceValue(PointType p); void CreateDistanceImage (); //Datastructures for the interpolation CenterList m_Centers; NormalList m_Normals; FunctionValues m_FunctionValues; InterpolationWeights m_Weights; SolutionMatrix m_SolutionMatrix; double m_DistanceImageSpacing; unsigned int m_DistanceImageVolume; + bool m_UseProgressBar; + unsigned int m_ProgressStepSize; }; }//namespace #endif diff --git a/Modules/MitkExt/Algorithms/mitkImageToContourFilter.cpp b/Modules/MitkExt/Algorithms/mitkImageToContourFilter.cpp index cf73d38924..72a235cb8a 100644 --- a/Modules/MitkExt/Algorithms/mitkImageToContourFilter.cpp +++ b/Modules/MitkExt/Algorithms/mitkImageToContourFilter.cpp @@ -1,120 +1,136 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Module: $RCSfile$ Language: C++ Date: $Date: $ Version: $Revision: $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #include "mitkImageToContourFilter.h" #include "mitkImageCast.h" #include "vtkSmartPointer.h" #include "vtkMatrix4x4.h" #include "vtkLinearTransform.h" #include "mitkVtkRepresentationProperty.h" #include "vtkProperty.h" #include "mitkImageAccessByItk.h" mitk::ImageToContourFilter::ImageToContourFilter() { + this->m_UseProgressBar = false; + this->m_ProgressStepSize = 1; } mitk::ImageToContourFilter::~ImageToContourFilter() { } void mitk::ImageToContourFilter::GenerateData() { mitk::Image::ConstPointer sliceImage = ImageToSurfaceFilter::GetInput(); if ( !sliceImage ) { MITK_ERROR << "mitk::ImageToContourFilter: No input available. Please set the input!" << std::endl; itkExceptionMacro("mitk::ImageToContourFilter: No input available. Please set the input!"); return; } if ( sliceImage->GetDimension() > 2 || sliceImage->GetDimension() < 2) { MITK_ERROR << "mitk::ImageToImageFilter::GenerateData() works only with 2D images. Please assure that your input image is 2D!" << std::endl; itkExceptionMacro("mitk::ImageToImageFilter::GenerateData() works only with 2D images. Please assure that your input image is 2D!"); return; } m_SliceGeometry = sliceImage->GetGeometry(); AccessFixedDimensionByItk(sliceImage, Itk2DContourExtraction, 2); + + //Setting progressbar + if (this->m_UseProgressBar) + mitk::ProgressBar::GetInstance()->Progress(this->m_ProgressStepSize); } template void mitk::ImageToContourFilter::Itk2DContourExtraction (itk::Image* sliceImage) { typedef itk::Image ImageType; typedef itk::ContourExtractor2DImageFilter ContourExtractor; typename ContourExtractor::Pointer contourExtractor = ContourExtractor::New(); contourExtractor->SetInput(sliceImage); contourExtractor->SetContourValue(0.5); contourExtractor->Update(); unsigned int foundPaths = contourExtractor->GetNumberOfOutputs(); vtkSmartPointer contourSurface = vtkSmartPointer::New(); vtkSmartPointer points = vtkSmartPointer::New(); vtkSmartPointer polygons = vtkSmartPointer::New(); unsigned int pointId (0); for (unsigned int i = 0; i < foundPaths; i++) { const ContourPath* currentPath = contourExtractor->GetOutput(i)->GetVertexList(); vtkSmartPointer polygon = vtkSmartPointer::New(); polygon->GetPointIds()->SetNumberOfIds(currentPath->Size()); Point3D currentPoint; Point3D currentWorldPoint; for (unsigned int j = 0; j < currentPath->Size(); j++) { currentPoint[0] = currentPath->ElementAt(j)[0] + 0.5; currentPoint[1] = currentPath->ElementAt(j)[1] + 0.5; currentPoint[2] = 0; m_SliceGeometry->IndexToWorld(currentPoint, currentWorldPoint); points->InsertPoint(pointId, currentWorldPoint[0],currentWorldPoint[1],currentWorldPoint[2]); polygon->GetPointIds()->SetId(j,pointId); pointId++; }//for2 polygons->InsertNextCell(polygon); }//for1 contourSurface->SetPoints( points ); contourSurface->SetPolys( polygons ); contourSurface->BuildLinks(); Surface::Pointer finalSurface = this->GetOutput(); finalSurface->SetVtkPolyData( contourSurface ); } void mitk::ImageToContourFilter::GenerateOutputInformation() { Superclass::GenerateOutputInformation(); } +void mitk::ImageToContourFilter::SetUseProgressBar(bool status) +{ + this->m_UseProgressBar = status; +} + +void mitk::ImageToContourFilter::SetProgressStepSize(unsigned int stepSize) +{ + this->m_ProgressStepSize = stepSize; +} + diff --git a/Modules/MitkExt/Algorithms/mitkImageToContourFilter.h b/Modules/MitkExt/Algorithms/mitkImageToContourFilter.h index 0a1aecd843..e7022b7470 100644 --- a/Modules/MitkExt/Algorithms/mitkImageToContourFilter.h +++ b/Modules/MitkExt/Algorithms/mitkImageToContourFilter.h @@ -1,74 +1,92 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Module: $RCSfile$ Language: C++ Date: $Date: $ Version: $Revision: $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #ifndef mitkImageToContourFilter_h_Included #define mitkImageToContourFilter_h_Included //#include "MitkSBExports.h" #include "MitkExtExports.h" #include "itkImage.h" #include "mitkImage.h" #include "itkContourExtractor2DImageFilter.h" #include "mitkImageToSurfaceFilter.h" #include "mitkSurface.h" #include "vtkPolyData.h" #include "vtkPolygon.h" #include "vtkCellArray.h" +#include "mitkProgressBar.h" + namespace mitk { /** \brief A filter that can extract contours out of a 2D binary image This class takes an 2D mitk::Image as input and extracts all contours which are drawn it. The contour extraction is done by using the itk::ContourExtractor2DImageFilter. The output is a mitk::Surface. $Author: fetzer$ */ class MitkExt_EXPORT ImageToContourFilter : public ImageToSurfaceFilter { public: mitkClassMacro(ImageToContourFilter,ImageToSurfaceFilter); itkNewMacro(Self); /** \brief Set macro for the geometry of the slice. If it is not set explicitly the geometry will be taken from the slice \a Parameter The slice`s geometry */ itkSetMacro(SliceGeometry, Geometry3D*); //typedef unsigned int VDimension; typedef itk::PolyLineParametricPath<2> PolyLineParametricPath2D; typedef PolyLineParametricPath2D::VertexListType ContourPath; + + /** + \brief Set whether the mitkProgressBar should be used + + \a Parameter true for using the progress bar, false otherwise + */ + void SetUseProgressBar(bool); + + /** + \brief Set the stepsize which the progress bar should proceed + + \a Parameter The stepsize for progressing + */ + void SetProgressStepSize(unsigned int stepSize); protected: ImageToContourFilter(); virtual ~ImageToContourFilter(); virtual void GenerateData(); virtual void GenerateOutputInformation(); private: const Geometry3D* m_SliceGeometry; + bool m_UseProgressBar; + unsigned int m_ProgressStepSize; template void Itk2DContourExtraction (itk::Image* sliceImage); };//class }//namespace #endif diff --git a/Modules/MitkExt/Algorithms/mitkReduceContourSetFilter.cpp b/Modules/MitkExt/Algorithms/mitkReduceContourSetFilter.cpp index a8e9a3f883..3b4ee410e2 100644 --- a/Modules/MitkExt/Algorithms/mitkReduceContourSetFilter.cpp +++ b/Modules/MitkExt/Algorithms/mitkReduceContourSetFilter.cpp @@ -1,473 +1,489 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Module: $RCSfile$ Language: C++ Date: $Date: $ Version: $Revision: $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #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; } mitk::ReduceContourSetFilter::~ReduceContourSetFilter() { } void mitk::ReduceContourSetFilter::GenerateData() { unsigned int numberOfInputs = this->GetNumberOfInputs(); unsigned int numberOfOutputs (0); this->CreateOutputsForAllInputs(numberOfInputs); vtkSmartPointer newPolyData; vtkSmartPointer newPolygons; vtkSmartPointer newPoints; //For the purpose of evaluation // unsigned int numberOfPointsBefore (0); // unsigned int numberOfPointsAfter (0); for(unsigned int i = 0; i < numberOfInputs; i++) { mitk::Surface* currentSurface = const_cast( this->GetInput(i) ); vtkSmartPointer polyData = currentSurface->GetVtkPolyData(); newPolyData = vtkPolyData::New(); newPolygons = vtkCellArray::New(); newPoints = vtkPoints::New(); vtkSmartPointer existingPolys = polyData->GetPolys(); vtkSmartPointer existingPoints = polyData->GetPoints(); existingPolys->InitTraversal(); vtkIdType* cell (NULL); 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 = vtkPolygon::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; // numberOfPointsAfter += newPolygon->GetPointIds()->GetNumberOfIds(); } if (newPolygons->GetNumberOfCells() != 0) { newPolyData->SetPolys(newPolygons); newPolyData->SetPoints(newPoints); newPolyData->BuildLinks(); Surface::Pointer surface = this->GetOutput(numberOfOutputs); surface->SetVtkPolyData(newPolyData); numberOfOutputs++; } } // MITK_INFO<<"Points before: "<SetNumberOfOutputs(numberOfOutputs); + + //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) { 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 (unsigned int 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) { //If the cell is too small to obtain a reduced polygon with the given stepsize return if (cellSize <= 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, 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->GetNumberOfInputs(); i++) { //Don't check for intersection with the polygon itself if (i == currentInputIndex) continue; //Get the next polydata to check for intersection vtkSmartPointer poly = const_cast( this->GetInput(i) )->GetVtkPolyData(); vtkSmartPointer polygonArray = poly->GetPolys(); polygonArray->InitTraversal(); vtkIdType anotherInputPolygonSize (0); vtkIdType* anotherInputPolygonIDs(NULL); /* 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 (unsigned int 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 (unsigned int 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->GetNumberOfInputs(); i++) { this->PopBackInput(); } this->SetNumberOfInputs(0); this->SetNumberOfOutputs(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/MitkExt/Algorithms/mitkReduceContourSetFilter.h b/Modules/MitkExt/Algorithms/mitkReduceContourSetFilter.h index d245886653..02b2c5171f 100644 --- a/Modules/MitkExt/Algorithms/mitkReduceContourSetFilter.h +++ b/Modules/MitkExt/Algorithms/mitkReduceContourSetFilter.h @@ -1,104 +1,122 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Module: $RCSfile$ Language: C++ Date: $Date: $ Version: $Revision: $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #ifndef mitkReduceContourSetFilter_h_Included #define mitkReduceContourSetFilter_h_Included #include "MitkExtExports.h" #include "mitkSurfaceToSurfaceFilter.h" #include "mitkSurface.h" +#include "mitkProgressBar.h" #include "vtkSmartPointer.h" #include "vtkPolyData.h" #include "vtkPolygon.h" #include "vtkPoints.h" #include "vtkCellArray.h" #include "vtkMath.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 MitkExt_EXPORT ReduceContourSetFilter : public SurfaceToSurfaceFilter { public: enum Reduction_Type { NTH_POINT, DOUGLAS_PEUCKER }; struct LineSegment { unsigned int StartIndex; unsigned int EndIndex; }; mitkClassMacro(ReduceContourSetFilter,SurfaceToSurfaceFilter); itkNewMacro(Self); itkSetMacro(MinSpacing, double); itkSetMacro(MaxSpacing, double); itkSetMacro(ReductionType, Reduction_Type); itkSetMacro(StepSize, unsigned int); itkSetMacro(Tolerance, double); //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); + + /** + \brief Set the stepsize which the progress bar should proceed + + \a Parameter The stepsize for progressing + */ + void SetProgressStepSize(unsigned int stepSize); + protected: ReduceContourSetFilter(); virtual ~ReduceContourSetFilter(); virtual void GenerateData(); virtual void GenerateOutputInformation(); private: void ReduceNumberOfPointsByNthPoint (vtkIdType cellSize, vtkIdType* cell, vtkPoints* points, vtkPolygon* reducedPolygon, vtkPoints* reducedPoints); void ReduceNumberOfPointsByDouglasPeucker (vtkIdType cellSize, vtkIdType* cell, vtkPoints* points, vtkPolygon* reducedPolygon, vtkPoints* reducedPoints); bool CheckForIntersection (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; };//class }//namespace #endif diff --git a/Modules/MitkExt/Controllers/mitkSurfaceInterpolationController.cpp b/Modules/MitkExt/Controllers/mitkSurfaceInterpolationController.cpp index f4a2a55a34..86517908fb 100644 --- a/Modules/MitkExt/Controllers/mitkSurfaceInterpolationController.cpp +++ b/Modules/MitkExt/Controllers/mitkSurfaceInterpolationController.cpp @@ -1,211 +1,221 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date$ Version: $Revision$ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #include "mitkSurfaceInterpolationController.h" mitk::SurfaceInterpolationController::SurfaceInterpolationController() : m_CurrentContourListID (0) { m_ReduceFilter = ReduceContourSetFilter::New(); m_NormalsFilter = ComputeContourSetNormalsFilter::New(); m_InterpolateSurfaceFilter = CreateDistanceImageFromSurfaceFilter::New(); + m_ReduceFilter->SetUseProgressBar(true); + m_NormalsFilter->SetUseProgressBar(true); + m_InterpolateSurfaceFilter->SetUseProgressBar(true); + m_Contours = Surface::New(); m_PolyData = vtkSmartPointer::New(); m_PolyData->SetPoints(vtkPoints::New()); m_InterpolationResult = 0; } mitk::SurfaceInterpolationController::~SurfaceInterpolationController() { for (unsigned int i = 0; i < m_ListOfContourLists.size(); ++i) { for (unsigned int j = 0; j < m_ListOfContourLists.at(i).size(); ++j) { delete(m_ListOfContourLists.at(i).at(j).position); } } } mitk::SurfaceInterpolationController* mitk::SurfaceInterpolationController::GetInstance() { static mitk::SurfaceInterpolationController* m_Instance; if ( m_Instance == 0) { m_Instance = new SurfaceInterpolationController(); } return m_Instance; } void mitk::SurfaceInterpolationController::AddNewContour (mitk::Surface::Pointer newContour ,RestorePlanePositionOperation* op) { AffineTransform3D::Pointer transform = AffineTransform3D::New(); transform = op->GetTransform(); mitk::Vector3D direction = op->GetDirectionVector(); int pos (-1); for (unsigned int i = 0; i < m_ListOfContourLists.at(m_CurrentContourListID).size(); i++) { itk::Matrix diffM = transform->GetMatrix()-m_ListOfContourLists.at(m_CurrentContourListID).at(i).position->GetTransform()->GetMatrix(); bool isSameMatrix(true); for (unsigned int j = 0; j < 3; j++) { if (fabs(diffM[j][0]) > 0.0001 && fabs(diffM[j][1]) > 0.0001 && fabs(diffM[j][2]) > 0.0001) { isSameMatrix = false; break; } } itk::Vector diffV = m_ListOfContourLists.at(m_CurrentContourListID).at(i).position->GetTransform()->GetOffset()-transform->GetOffset(); if ( isSameMatrix && m_ListOfContourLists.at(m_CurrentContourListID).at(i).position->GetPos() == op->GetPos() && (fabs(diffV[0]) < 0.0001 && fabs(diffV[1]) < 0.0001 && fabs(diffV[2]) < 0.0001) ) { pos = i; break; } } if (pos == -1) { //MITK_INFO<<"New Contour"; mitk::RestorePlanePositionOperation* newOp = new mitk::RestorePlanePositionOperation (OpRESTOREPLANEPOSITION, op->GetWidth(), op->GetHeight(), op->GetSpacing(), op->GetPos(), direction, transform); ContourPositionPair newData; newData.contour = newContour; newData.position = newOp; m_ReduceFilter->SetInput(m_ListOfContourLists.at(m_CurrentContourListID).size(), newContour); m_NormalsFilter->SetInput(m_ListOfContourLists.at(m_CurrentContourListID).size(), m_ReduceFilter->GetOutput(m_ListOfContourLists.at(m_CurrentContourListID).size())); m_InterpolateSurfaceFilter->SetInput(m_ListOfContourLists.at(m_CurrentContourListID).size(), m_NormalsFilter->GetOutput(m_ListOfContourLists.at(m_CurrentContourListID).size())); m_ListOfContourLists.at(m_CurrentContourListID).push_back(newData); } else { //MITK_INFO<<"Modified Contour"; m_ListOfContourLists.at(m_CurrentContourListID).at(pos).contour = newContour; m_ReduceFilter->SetInput(pos, newContour); m_NormalsFilter->SetInput(pos, m_ReduceFilter->GetOutput(pos)); m_InterpolateSurfaceFilter->SetInput(pos, m_NormalsFilter->GetOutput(pos)); } this->Modified(); } void mitk::SurfaceInterpolationController::Interpolate() { if (m_ListOfContourLists.at(m_CurrentContourListID).size() < 2) return; + //Setting up progress bar + mitk::ProgressBar::GetInstance()->AddStepsToDo(8); + m_InterpolateSurfaceFilter->Update(); Image::Pointer distanceImage = m_InterpolateSurfaceFilter->GetOutput(); vtkSmartPointer mcFilter = vtkMarchingCubes::New(); mcFilter->SetInput(distanceImage->GetVtkImageData()); mcFilter->SetValue(0,0); mcFilter->Update(); m_InterpolationResult = 0; m_InterpolationResult = mitk::Surface::New(); m_InterpolationResult->SetVtkPolyData(mcFilter->GetOutput()); m_InterpolationResult->GetGeometry()->SetOrigin(distanceImage->GetGeometry()->GetOrigin()); vtkSmartPointer polyDataAppender = vtkSmartPointer::New(); for (unsigned int i = 0; i < m_ReduceFilter->GetNumberOfOutputs(); i++) { polyDataAppender->AddInput(m_ReduceFilter->GetOutput(i)->GetVtkPolyData()); } polyDataAppender->Update(); m_Contours->SetVtkPolyData(polyDataAppender->GetOutput()); + //Last progress step + mitk::ProgressBar::GetInstance()->Progress(); + m_InterpolationResult->DisconnectPipeline(); } mitk::Surface::Pointer mitk::SurfaceInterpolationController::GetInterpolationResult() { return m_InterpolationResult; } mitk::Surface* mitk::SurfaceInterpolationController::GetContoursAsSurface() { return m_Contours; } void mitk::SurfaceInterpolationController::SetDataStorage(DataStorage &ds) { m_DataStorage = &ds; } unsigned int mitk::SurfaceInterpolationController::CreateNewContourList() { unsigned int newID = m_ListOfContourLists.size(); ContourPositionPairList newList; m_ListOfContourLists.push_back(newList); this->SetCurrentListID(newID); m_InterpolationResult = 0; return m_CurrentContourListID; } void mitk::SurfaceInterpolationController::SetCurrentListID ( unsigned int ID ) { if (ID == m_CurrentContourListID ) return; m_CurrentContourListID = ID; m_ReduceFilter->Reset(); m_NormalsFilter->Reset(); m_InterpolateSurfaceFilter->Reset(); for (unsigned int i = 0; i < m_ListOfContourLists.at(m_CurrentContourListID).size(); i++) { m_ReduceFilter->SetInput(i, m_ListOfContourLists.at(m_CurrentContourListID).at(i).contour); m_NormalsFilter->SetInput(i,m_ReduceFilter->GetOutput(i)); m_InterpolateSurfaceFilter->SetInput(i,m_NormalsFilter->GetOutput(i)); } } void mitk::SurfaceInterpolationController::SetMinSpacing(double minSpacing) { m_ReduceFilter->SetMinSpacing(minSpacing); } void mitk::SurfaceInterpolationController::SetMaxSpacing(double maxSpacing) { m_ReduceFilter->SetMaxSpacing(maxSpacing); m_NormalsFilter->SetMaxSpacing(maxSpacing); } void mitk::SurfaceInterpolationController::SetDistanceImageVolume(unsigned int distImgVolume) { m_InterpolateSurfaceFilter->SetDistanceImageVolume(distImgVolume); } void mitk::SurfaceInterpolationController::SetWorkingImage(Image* workingImage) { m_NormalsFilter->SetSegmentationBinaryImage(workingImage); } mitk::Image* mitk::SurfaceInterpolationController::GetImage() { return m_InterpolateSurfaceFilter->GetOutput(); } diff --git a/Modules/MitkExt/Controllers/mitkSurfaceInterpolationController.h b/Modules/MitkExt/Controllers/mitkSurfaceInterpolationController.h index afebd1b5d3..583eca813f 100644 --- a/Modules/MitkExt/Controllers/mitkSurfaceInterpolationController.h +++ b/Modules/MitkExt/Controllers/mitkSurfaceInterpolationController.h @@ -1,131 +1,132 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date$ Version: $Revision$ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #ifndef mitkSurfaceInterpolationController_h_Included #define mitkSurfaceInterpolationController_h_Included #include "mitkCommon.h" #include "MitkExtExports.h" #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 "mitkProgressBar.h" namespace mitk { class MitkExt_EXPORT SurfaceInterpolationController : public itk::Object { public: mitkClassMacro(SurfaceInterpolationController, itk::Object); itkNewMacro(Self); static SurfaceInterpolationController* GetInstance(); /* * Adds a new extracted contour to the list */ void AddNewContour(Surface::Pointer newContour, RestorePlanePositionOperation *op); /* * Interpolates the 3D surface from the given extracted contours */ void Interpolate (); mitk::Surface::Pointer GetInterpolationResult(); void SetMinSpacing(double minSpacing); void SetMaxSpacing(double maxSpacing); void SetDistanceImageVolume(unsigned int distImageVolume); void SetWorkingImage(Image* workingImage); Surface* GetContoursAsSurface(); void SetDataStorage(DataStorage &ds); unsigned int CreateNewContourList(); void SetCurrentListID (unsigned int ID); mitk::Image* GetImage(); protected: SurfaceInterpolationController(); ~SurfaceInterpolationController(); private: struct ContourPositionPair { Surface::Pointer contour; RestorePlanePositionOperation* position; }; typedef std::vector ContourPositionPairList; ContourPositionPairList::iterator m_Iterator; ReduceContourSetFilter::Pointer m_ReduceFilter; ComputeContourSetNormalsFilter::Pointer m_NormalsFilter; CreateDistanceImageFromSurfaceFilter::Pointer m_InterpolateSurfaceFilter; double m_MinSpacing; double m_MaxSpacing; const Image* m_WorkingImage; Surface::Pointer m_Contours; vtkSmartPointer m_PolyData; unsigned int m_DistImageVolume; mitk::WeakPointer m_DataStorage; std::vector m_ListOfContourLists; unsigned int m_CurrentContourListID; mitk::Surface::Pointer m_InterpolationResult; }; } #endif