diff --git a/Modules/Segmentation/Algorithms/mitkCreateDistanceImageFromSurfaceFilter.cpp b/Modules/Segmentation/Algorithms/mitkCreateDistanceImageFromSurfaceFilter.cpp index 6c836a7d43..18eed7ea10 100644 --- a/Modules/Segmentation/Algorithms/mitkCreateDistanceImageFromSurfaceFilter.cpp +++ b/Modules/Segmentation/Algorithms/mitkCreateDistanceImageFromSurfaceFilter.cpp @@ -1,517 +1,592 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) -Copyright (c) German Cancer Research Center, +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 +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 "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); + 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(); - //Determine 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]; + // Determine the bounds of the input points in index- and world-coordinates + DistanceImageType::PointType minPointInWorldCoordinates, maxPointInWorldCoordinates; + DistanceImageType::IndexType minPointInIndexCoordinates, maxPointInIndexCoordinates; - CenterList::iterator centerIter = m_Centers.begin(); + DetermineBounds( minPointInWorldCoordinates, maxPointInWorldCoordinates, + minPointInIndexCoordinates, maxPointInIndexCoordinates ); - for ( ++centerIter; centerIter!=m_Centers.end(); centerIter++) - { - if (xmin > (*centerIter)[0]) - { - xmin = (*centerIter)[0]; - } - if (ymin > (*centerIter)[1]) - { - ymin = (*centerIter)[1]; - } - if (zmin > (*centerIter)[2]) - { - zmin = (*centerIter)[2]; - } - if (xmax < (*centerIter)[0]) - { - xmax = (*centerIter)[0]; - } - if (ymax < (*centerIter)[1]) - { - ymax = (*centerIter)[1]; - } - if (zmax < (*centerIter)[2]) - { - zmax = (*centerIter)[2]; - } - } + // 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; - extentMM[0] = xmax - xmin + 5; - extentMM[1] = ymax - ymin + 5; - extentMM[2] = zmax - zmin + 5; - - //Shifting the distance image's offset to achieve an exact distance calculation - xmin = xmin - 2; - ymin = ymin - 2; - zmin = zmin - 2; + extentMM[0] = ( std::max(abs(maxPointInIndexCoordinates[0] - minPointInIndexCoordinates[0]), (long) 1) ) * m_ReferenceImage->GetSpacing()[0]; + extentMM[1] = ( std::max(abs(maxPointInIndexCoordinates[1] - minPointInIndexCoordinates[1]), (long) 1) ) * m_ReferenceImage->GetSpacing()[1]; + extentMM[2] = ( std::max(abs(maxPointInIndexCoordinates[2] - minPointInIndexCoordinates[2]), (long) 1) ) * m_ReferenceImage->GetSpacing()[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) + * Now create an empty distance image. The create image will always have the same sizeOfRegion, 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; + // 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; - DistanceImageType::SizeType size; + // 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 + 4; + sizeOfRegion[1] = numberOfYPixel + 4; + sizeOfRegion[2] = numberOfZPixel + 4; - //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; + // The region starts at index 0,0,0 + DistanceImageType::IndexType initialOriginAsIndex; + initialOriginAsIndex.Fill(0); - DistanceImageType::IndexType start; - start[0] = 0; - start[1] = 0; - start[2] = 0; + DistanceImageType::PointType originAsWorld = minPointInWorldCoordinates; DistanceImageType::RegionType lpRegion; - - lpRegion.SetSize(size); - lpRegion.SetIndex(start); - + 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 + distanceImg->SetOrigin( originAsWorld ); + distanceImg->SetDirection( m_ReferenceImage->GetDirection() ); 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 + // Now we move the origin of the distanceImage 2 index-Coordinates + // in all directions + DistanceImageType::IndexType originAsIndex; + distanceImg->TransformPhysicalPointToIndex( originAsWorld, originAsIndex ); + originAsIndex[0] -= 2; + originAsIndex[1] -= 2; + originAsIndex[2] -= 2; + distanceImg->TransformIndexToPhysicalPoint( originAsIndex, originAsWorld ); + distanceImg->SetOrigin( originAsWorld ); - This is done until the narrowband_point_list is empty. + /* + * 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. */ 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; - currentIndex[0] = ( currentPoint[0]-xmin ) / m_DistanceImageSpacing; - currentIndex[1] = ( currentPoint[1]-ymin ) / m_DistanceImageSpacing; - currentIndex[2] = ( currentPoint[2]-zmin ) / m_DistanceImageSpacing; + distanceImg->TransformPhysicalPointToIndex( currentPointAsPoint, currentIndex ); 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(); + unsigned int* relativeNb = &relativeNbIdx[0]; for (int i = 0; i < 6; i++) { - nIt.GetPixel(relativeNbIdx[i], isInBounds); - if( isInBounds && nIt.GetPixel(relativeNbIdx[i]) == 10) + nIt.GetPixel(*relativeNb, isInBounds); + if( isInBounds && nIt.GetPixel(*relativeNb) == 10) { - currentIndex = nIt.GetIndex(relativeNbIdx[i]); + currentIndex = nIt.GetIndex(*relativeNb); + + // Transform the currently checked point from index-coordinates to + // world-coordinates + distanceImg->TransformIndexToPhysicalPoint( currentIndex, currentPointAsPoint ); - currentPoint[0] = currentIndex[0]*m_DistanceImageSpacing + xmin; - currentPoint[1] = currentIndex[1]*m_DistanceImageSpacing + ymin; - currentPoint[2] = currentIndex[2]*m_DistanceImageSpacing + zmin; + // 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 ( abs(distance) <= m_DistanceImageSpacing*2 ) { - nIt.SetPixel(relativeNbIdx[i], distance); + nIt.SetPixel(*relativeNb, distance); narrowbandPoints.push(currentIndex); } } + relativeNb++; } } 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) }; + unsigned int _size[3] = { (unsigned int)(sizeOfRegion[0] - 1), (unsigned int)(sizeOfRegion[1] - 1), (unsigned int)(sizeOfRegion[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] + 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] + 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; - + // Cast the created distance-Image from itk::Image to the mitk::Image + // that is our output. CastToMitkImage(distanceImg, resultImage); - resultImage->GetGeometry()->SetOrigin(origin); } double mitk::CreateDistanceImageFromSurfaceFilter::CalculateDistanceValue(PointType p) { double distanceValue (0); PointType p1; PointType p2; double norm; CenterList::iterator centerIter; InterpolationWeights::iterator weightsIter; unsigned int i = 0; for ( centerIter=m_Centers.begin(), weightsIter=m_Weights.begin(); centerIter!=m_Centers.end(), weightsIter!=m_Weights.end(); centerIter++, weightsIter++ ) { p1 = *centerIter; p2 = p-p1; norm = p2.two_norm(); distanceValue = distanceValue + norm* (*weightsIter); } 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; } + +void mitk::CreateDistanceImageFromSurfaceFilter::SetReferenceImage( itk::Image::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 + DistanceImageType::IndexType tmpIndex; + m_ReferenceImage->TransformPhysicalPointToIndex( tmpPoint, tmpIndex ); + + // initialize the variables with this first point + int xmin = tmpIndex[0]; + int ymin = tmpIndex[1]; + int zmin = tmpIndex[2]; + int xmax = tmpIndex[0]; + int ymax = tmpIndex[1]; + int zmax = tmpIndex[2]; + + // iterate over the rest of the points + CenterList::iterator 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->TransformPhysicalPointToIndex( 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/Segmentation/Algorithms/mitkCreateDistanceImageFromSurfaceFilter.h b/Modules/Segmentation/Algorithms/mitkCreateDistanceImageFromSurfaceFilter.h index d6980e4074..510f1d0f74 100644 --- a/Modules/Segmentation/Algorithms/mitkCreateDistanceImageFromSurfaceFilter.h +++ b/Modules/Segmentation/Algorithms/mitkCreateDistanceImageFromSurfaceFilter.h @@ -1,161 +1,189 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) -Copyright (c) German Cancer Research Center, +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 +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 mitkCreateDistanceImageFromSurfaceFilter_h_Included #define mitkCreateDistanceImageFromSurfaceFilter_h_Included #include "SegmentationExports.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: + 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 Segmentation_EXPORT CreateDistanceImageFromSurfaceFilter : public ImageSource { public: typedef vnl_vector_fixed PointType; + typedef itk::Image DistanceImageType; + //typedef DistanceImageType::PointType PointType; + typedef DistanceImageType::IndexType IndexType; + 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); + void SetReferenceImage( itk::Image::Pointer referenceImage ); + + protected: CreateDistanceImageFromSurfaceFilter(); virtual ~CreateDistanceImageFromSurfaceFilter(); virtual void GenerateData(); virtual void GenerateOutputInformation(); private: void CreateSolutionMatrixAndFunctionValues(); double CalculateDistanceValue(PointType p); void CreateDistanceImage (); + /** + * \brief This method fills the given variables with the minimum and + * maximum coordinates that contain all input-points in index- and + * world-coordinates. + * + * This method iterates over all input-points and transforms them from + * world-coordinates to index-coordinates using the transform of the + * reference-Image. + * Next, the minimal and maximal index-coordinates are determined that + * span an area that contains all given input-points. + * These index-coordinates are then transformed back to world-coordinates. + * + * These minimal and maximal points are then set to the given variables. + */ + void DetermineBounds( DistanceImageType::PointType &minPointInWorldCoordinates, + DistanceImageType::PointType &maxPointInWorldCoordinates, + DistanceImageType::IndexType &minPointInIndexCoordinates, + DistanceImageType::IndexType &maxPointInIndexCoordinates ); + //Datastructures for the interpolation CenterList m_Centers; NormalList m_Normals; FunctionValues m_FunctionValues; InterpolationWeights m_Weights; SolutionMatrix m_SolutionMatrix; double m_DistanceImageSpacing; + itk::Image::Pointer m_ReferenceImage; + unsigned int m_DistanceImageVolume; bool m_UseProgressBar; unsigned int m_ProgressStepSize; }; }//namespace #endif