diff --git a/Core/Code/Algorithms/mitkClippedSurfaceBoundsCalculator.cpp b/Core/Code/Algorithms/mitkClippedSurfaceBoundsCalculator.cpp index c0b0b1d0e0..8cf5763fd9 100644 --- a/Core/Code/Algorithms/mitkClippedSurfaceBoundsCalculator.cpp +++ b/Core/Code/Algorithms/mitkClippedSurfaceBoundsCalculator.cpp @@ -1,291 +1,354 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkClippedSurfaceBoundsCalculator.h" #include "mitkLine.h" #define ROUND_P(x) ((x)>=0?(int)((x)+0.5):(int)((x)-0.5)) mitk::ClippedSurfaceBoundsCalculator::ClippedSurfaceBoundsCalculator( const mitk::PlaneGeometry* geometry, - mitk::Image::Pointer image): - m_PlaneGeometry(NULL), - m_Geometry3D(NULL), - m_Image(NULL) + mitk::Image::Pointer image) +: m_PlaneGeometry(NULL) +, m_Geometry3D(NULL) +, m_Image(NULL) { - // initialize with meaningless slice indices - m_MinMaxOutput.clear(); - for(int i = 0; i < 3; i++) - { - m_MinMaxOutput.push_back( - OutputType( std::numeric_limits::max() , - std::numeric_limits::min() )); - } - + this->InitializeOutput(); this->SetInput(geometry, image); } mitk::ClippedSurfaceBoundsCalculator::ClippedSurfaceBoundsCalculator( const mitk::Geometry3D* geometry, - mitk::Image::Pointer image): - m_PlaneGeometry(NULL), - m_Geometry3D(NULL), - m_Image(NULL) + mitk::Image::Pointer image) +: m_PlaneGeometry(NULL) +, m_Geometry3D(NULL) +, m_Image(NULL) +{ + this->InitializeOutput(); + + this->SetInput(geometry, image); +} + +mitk::ClippedSurfaceBoundsCalculator::ClippedSurfaceBoundsCalculator( const PointListType pointlist, mitk::Image::Pointer image ) +: m_PlaneGeometry(NULL) +, m_Geometry3D(NULL) +, m_Image(image) +{ + this->InitializeOutput(); + + m_ObjectPointsInWorldCoordinates = pointlist; +} + +void mitk::ClippedSurfaceBoundsCalculator::InitializeOutput() { // initialize with meaningless slice indices m_MinMaxOutput.clear(); for(int i = 0; i < 3; i++) { m_MinMaxOutput.push_back( - OutputType( std::numeric_limits::max() , - std::numeric_limits::min() )); + OutputType( std::numeric_limits::max() , + std::numeric_limits::min() )); } - - - this->SetInput(geometry, image); } mitk::ClippedSurfaceBoundsCalculator::~ClippedSurfaceBoundsCalculator() { } void mitk::ClippedSurfaceBoundsCalculator::SetInput( const mitk::PlaneGeometry* geometry, mitk::Image* image) { if(geometry && image) { this->m_PlaneGeometry = geometry; this->m_Image = image; this->m_Geometry3D = NULL; //Not possible to set both + m_ObjectPointsInWorldCoordinates.clear(); } } void mitk::ClippedSurfaceBoundsCalculator::SetInput( const mitk::Geometry3D* geometry, mitk::Image* image) { if(geometry && image) { this->m_Geometry3D = geometry; this->m_Image = image; this->m_PlaneGeometry = NULL; //Not possible to set both + m_ObjectPointsInWorldCoordinates.clear(); + } +} + +void mitk::ClippedSurfaceBoundsCalculator::SetInput( const std::vector pointlist, mitk::Image *image ) +{ + if ( !pointlist.empty() && image ) + { + m_Geometry3D = NULL; + m_PlaneGeometry = NULL; + m_Image = image; + m_ObjectPointsInWorldCoordinates = pointlist; } } mitk::ClippedSurfaceBoundsCalculator::OutputType mitk::ClippedSurfaceBoundsCalculator::GetMinMaxSpatialDirectionX() { return this->m_MinMaxOutput[0]; } mitk::ClippedSurfaceBoundsCalculator::OutputType mitk::ClippedSurfaceBoundsCalculator::GetMinMaxSpatialDirectionY() { return this->m_MinMaxOutput[1]; } mitk::ClippedSurfaceBoundsCalculator::OutputType mitk::ClippedSurfaceBoundsCalculator::GetMinMaxSpatialDirectionZ() { return this->m_MinMaxOutput[2]; } void mitk::ClippedSurfaceBoundsCalculator::Update() { this->m_MinMaxOutput.clear(); for(int i = 0; i < 3; i++) { this->m_MinMaxOutput.push_back(OutputType( std::numeric_limits::max() , std::numeric_limits::min() )); } if(m_PlaneGeometry.IsNotNull()) { this->CalculateIntersectionPoints(m_PlaneGeometry); } else if(m_Geometry3D.IsNotNull()) { // go through all slices of the image, ... const mitk::SlicedGeometry3D* slicedGeometry3D = dynamic_cast( m_Geometry3D.GetPointer() ); int allSlices = slicedGeometry3D->GetSlices(); this->CalculateIntersectionPoints(dynamic_cast(slicedGeometry3D->GetGeometry2D(0))); this->CalculateIntersectionPoints(dynamic_cast(slicedGeometry3D->GetGeometry2D(allSlices-1))); } + else if( !m_ObjectPointsInWorldCoordinates.empty() ) + { + this->CalculateIntersectionPoints( m_ObjectPointsInWorldCoordinates ); + } } void mitk::ClippedSurfaceBoundsCalculator::CalculateIntersectionPoints(const mitk::PlaneGeometry* geometry) { // SEE HEADER DOCUMENTATION for explanation typedef std::vector< std::pair > EdgesVector; Point3D origin; Vector3D xDirection, yDirection, zDirection; const Vector3D spacing = m_Image->GetGeometry()->GetSpacing(); origin = m_Image->GetGeometry()->GetOrigin(); //Left, bottom, front //Get axis vector for the spatial directions xDirection = m_Image->GetGeometry()->GetAxisVector(1); yDirection = m_Image->GetGeometry()->GetAxisVector(0); zDirection = m_Image->GetGeometry()->GetAxisVector(2); /* * For the calculation of the intersection points we need as corner points the center-based image coordinates. * With the method GetCornerPoint() of the class Geometry3D we only get the corner-based coordinates. * Therefore we need to calculate the center-based corner points here. For that we add/substract the corner- * based coordinates with the spacing of the geometry3D. */ for( int i = 0; i < 3; i++ ) { if(xDirection[i] < 0) { xDirection[i] += spacing[i]; } else if( xDirection[i] > 0 ) { xDirection[i] -= spacing[i]; } if(yDirection[i] < 0) { yDirection[i] += spacing[i]; } else if( yDirection[i] > 0 ) { yDirection[i] -= spacing[i]; } if(zDirection[i] < 0) { zDirection[i] += spacing[i]; } else if( zDirection[i] > 0 ) { zDirection[i] -= spacing[i]; } } Point3D leftBottomFront, leftTopFront, leftBottomBack, leftTopBack; Point3D rightBottomFront, rightTopFront, rightBottomBack, rightTopBack; leftBottomFront = origin; leftTopFront = origin + yDirection; leftBottomBack = origin + zDirection; leftTopBack = origin + yDirection + zDirection; rightBottomFront = origin + xDirection; rightTopFront = origin + xDirection + yDirection; rightBottomBack = origin + xDirection + zDirection; rightTopBack = origin + xDirection + yDirection + zDirection; EdgesVector edgesOf3DBox; edgesOf3DBox.push_back(std::make_pair(leftBottomBack, // x = left=xfront, y=bottom=yfront, z=front=zfront leftTopFront)); // left, top, front edgesOf3DBox.push_back(std::make_pair(leftBottomFront, // left, bottom, front leftBottomBack)); // left, bottom, back edgesOf3DBox.push_back(std::make_pair(leftBottomFront, // left, bottom, front rightBottomFront)); // right, bottom, front edgesOf3DBox.push_back(std::make_pair(leftTopFront, // left, top, front rightTopFront)); // right, top, front edgesOf3DBox.push_back(std::make_pair(leftTopFront, // left, top, front leftTopBack)); // left, top, back edgesOf3DBox.push_back(std::make_pair(rightTopFront, // right, top, front rightTopBack)); // right, top, back edgesOf3DBox.push_back(std::make_pair(rightTopFront, // right, top, front rightBottomFront)); // right, bottom, front edgesOf3DBox.push_back(std::make_pair(rightBottomFront, // right, bottom, front rightBottomBack)); // right, bottom, back edgesOf3DBox.push_back(std::make_pair(rightBottomBack, // right, bottom, back leftBottomBack)); // left, bottom, back edgesOf3DBox.push_back(std::make_pair(rightBottomBack, // right, bottom, back rightTopBack)); // right, top, back edgesOf3DBox.push_back(std::make_pair(rightTopBack, // right, top, back leftTopBack)); // left, top, back edgesOf3DBox.push_back(std::make_pair(leftTopBack, // left, top, back leftBottomBack)); // left, bottom, back for (EdgesVector::iterator iterator = edgesOf3DBox.begin(); iterator != edgesOf3DBox.end();iterator++) { Point3D startPoint = (*iterator).first; // start point of the line Point3D endPoint = (*iterator).second; // end point of the line Vector3D lineDirection = endPoint - startPoint; mitk::Line3D line(startPoint, lineDirection); Point3D intersectionWorldPoint; intersectionWorldPoint.Fill(std::numeric_limits::min()); // Get intersection point of line and plane geometry geometry->IntersectionPoint(line, intersectionWorldPoint); double t = -1.0; bool doesLineIntersectWithPlane(false); if(line.GetDirection().GetNorm() < mitk::eps && geometry->Distance(line.GetPoint1()) < mitk::sqrteps) { t = 1.0; doesLineIntersectWithPlane = true; intersectionWorldPoint = line.GetPoint1(); } else { geometry->IntersectionPoint(line, intersectionWorldPoint); doesLineIntersectWithPlane = geometry->IntersectionPointParam(line, t); } mitk::Point3D intersectionIndexPoint; //Get index point m_Image->GetGeometry()->WorldToIndex(intersectionWorldPoint, intersectionIndexPoint); if ( doesLineIntersectWithPlane && -mitk::sqrteps <= t && t <= 1.0 + mitk::sqrteps ) { for(int dim = 0; dim < 3; dim++) { // minimum //If new point value is lower than old if( this->m_MinMaxOutput[dim].first > ROUND_P(intersectionIndexPoint[dim]) ) { this->m_MinMaxOutput[dim].first = ROUND_P(intersectionIndexPoint[dim]); //set new value } // maximum //If new point value is higher than old if( this->m_MinMaxOutput[dim].second < ROUND_P(intersectionIndexPoint[dim]) ) { this->m_MinMaxOutput[dim].second = ROUND_P(intersectionIndexPoint[dim]); //set new value } } + + this->EnforceImageBounds(); } } } +void mitk::ClippedSurfaceBoundsCalculator::CalculateIntersectionPoints( PointListType pointList ) +{ + PointListType::iterator pointIterator; + + mitk::SlicedGeometry3D::Pointer imageGeometry = m_Image->GetSlicedGeometry(); + + + for ( pointIterator = pointList.begin(); pointIterator != pointList.end(); pointIterator++ ) + { + mitk::Point3D pntInIndexCoordinates; + imageGeometry->WorldToIndex( (*pointIterator), pntInIndexCoordinates ); + + m_MinMaxOutput[0].first = pntInIndexCoordinates[0] < m_MinMaxOutput[0].first ? ROUND_P(pntInIndexCoordinates[0]) : m_MinMaxOutput[0].first; + m_MinMaxOutput[0].second = pntInIndexCoordinates[0] > m_MinMaxOutput[0].second ? ROUND_P(pntInIndexCoordinates[0]) : m_MinMaxOutput[0].second; + + m_MinMaxOutput[1].first = pntInIndexCoordinates[1] < m_MinMaxOutput[1].first ? ROUND_P(pntInIndexCoordinates[1]) : m_MinMaxOutput[1].first; + m_MinMaxOutput[1].second = pntInIndexCoordinates[1] > m_MinMaxOutput[1].second ? ROUND_P(pntInIndexCoordinates[1]) : m_MinMaxOutput[1].second; + + m_MinMaxOutput[2].first = pntInIndexCoordinates[2] < m_MinMaxOutput[2].first ? ROUND_P(pntInIndexCoordinates[2]) : m_MinMaxOutput[2].first; + m_MinMaxOutput[2].second = pntInIndexCoordinates[2] > m_MinMaxOutput[2].second ? ROUND_P(pntInIndexCoordinates[2]) : m_MinMaxOutput[2].second; + } + + this->EnforceImageBounds(); + +} + +void mitk::ClippedSurfaceBoundsCalculator::EnforceImageBounds() +{ + m_MinMaxOutput[0].first = std::max( m_MinMaxOutput[0].first, 0 ); + m_MinMaxOutput[1].first = std::max( m_MinMaxOutput[1].first, 0 ); + m_MinMaxOutput[2].first = std::max( m_MinMaxOutput[2].first, 0 ); + + m_MinMaxOutput[0].second = std::min( m_MinMaxOutput[0].second, (int) m_Image->GetDimension(0)-1 ); + m_MinMaxOutput[1].second = std::min( m_MinMaxOutput[1].second, (int) m_Image->GetDimension(1)-1 ); + m_MinMaxOutput[2].second = std::min( m_MinMaxOutput[2].second, (int) m_Image->GetDimension(2)-1 ); +} + + diff --git a/Core/Code/Algorithms/mitkClippedSurfaceBoundsCalculator.h b/Core/Code/Algorithms/mitkClippedSurfaceBoundsCalculator.h index 1341ad26ef..690c10e554 100644 --- a/Core/Code/Algorithms/mitkClippedSurfaceBoundsCalculator.h +++ b/Core/Code/Algorithms/mitkClippedSurfaceBoundsCalculator.h @@ -1,105 +1,122 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef ClippedSurfaceBoundsCalculator_h_included #define ClippedSurfaceBoundsCalculator_h_included #include "mitkImage.h" #include "mitkPlaneGeometry.h" #include /** * \brief Find image slices visible on a given plane. * * The class name is not helpful in finding this class. Good suggestions welcome. * * Given a PlaneGeometry (e.g. the 2D plane of a render window), this class * calculates which slices of an mitk::Image are visible on this plane. * Calculation is done for X, Y, and Z direction, the result is available in * form of a pair (minimum,maximum) slice index. * * Such calculations are useful if you want to display information about the * currently visible slice (overlays, statistics, ...) and you don't want to * depend on any prior information about hat the renderwindow is currently showing. * * \warning The interface attempts to look like an ITK filter but it is far from being one. */ namespace mitk { class MITK_CORE_EXPORT ClippedSurfaceBoundsCalculator { public: + typedef std::vector PointListType; + ClippedSurfaceBoundsCalculator(const mitk::PlaneGeometry* geometry = NULL, mitk::Image::Pointer image = NULL); ClippedSurfaceBoundsCalculator(const mitk::Geometry3D* geometry, mitk::Image::Pointer image); + ClippedSurfaceBoundsCalculator(const PointListType pointlist, + mitk::Image::Pointer image); + + void InitializeOutput(); + virtual ~ClippedSurfaceBoundsCalculator(); void SetInput(const mitk::PlaneGeometry* geometry, mitk::Image* image); void SetInput(const mitk::Geometry3D *geometry, mitk::Image *image); + void SetInput(const PointListType pointlist, mitk::Image *image); /** \brief Request calculation. How cut/visible slice indices are determined: 1. construct a bounding box of the image. This is the box that connect the outer voxel centers(!). 2. check the edges of this box. 3. intersect each edge with the plane geometry - if the intersection point is within the image box, we update the visible/cut slice indices for all dimensions. - else we ignore the cut */ void Update(); /** \brief Minimum (first) and maximum (second) slice index. */ typedef std::pair OutputType; /** \brief What X coordinates (slice indices) are cut/visible in given plane. */ OutputType GetMinMaxSpatialDirectionX(); /** \brief What Y coordinates (slice indices) are cut/visible in given plane. */ OutputType GetMinMaxSpatialDirectionY(); /** \brief What Z coordinates (slice indices) are cut/visible in given plane. */ OutputType GetMinMaxSpatialDirectionZ(); protected: void CalculateIntersectionPoints(const mitk::PlaneGeometry* geometry); + void CalculateIntersectionPoints( PointListType pointList ); + + /** + * \brief Clips the resulting index-coordinates to make sure they do + * not exceed the imagebounds. + */ + void EnforceImageBounds(); + mitk::PlaneGeometry::ConstPointer m_PlaneGeometry; mitk::Geometry3D::ConstPointer m_Geometry3D; mitk::Image::Pointer m_Image; + std::vector m_ObjectPointsInWorldCoordinates; std::vector< OutputType > m_MinMaxOutput; }; } //namespace mitk #endif diff --git a/Core/Code/Testing/mitkClippedSurfaceBoundsCalculatorTest.cpp b/Core/Code/Testing/mitkClippedSurfaceBoundsCalculatorTest.cpp index 5c8b7456c1..13d97a5ff8 100644 --- a/Core/Code/Testing/mitkClippedSurfaceBoundsCalculatorTest.cpp +++ b/Core/Code/Testing/mitkClippedSurfaceBoundsCalculatorTest.cpp @@ -1,427 +1,484 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTestingMacros.h" #include #include "mitkClippedSurfaceBoundsCalculator.h" #include "mitkGeometry3D.h" #include "mitkGeometry2D.h" #include "mitkVector.h" static void CheckPlanesInsideBoundingBoxOnlyOnOneSlice(mitk::Geometry3D::Pointer geometry3D) { //Check planes which are inside the bounding box mitk::ClippedSurfaceBoundsCalculator* calculator = new mitk::ClippedSurfaceBoundsCalculator(); mitk::Image::Pointer image = mitk::Image::New(); - image->SetGeometry(geometry3D); + image->Initialize( mitk::MakePixelType(), *(geometry3D.GetPointer()) ); //Check planes which are only on one slice: //Slice 0 mitk::Point3D origin; origin[0] = 511; origin[1] = 0; origin[2] = 0; mitk::Vector3D normal; mitk::FillVector3D(normal, 0, 0, 1); mitk::PlaneGeometry::Pointer planeOnSliceZero = mitk::PlaneGeometry::New(); planeOnSliceZero->InitializePlane(origin, normal); calculator->SetInput( planeOnSliceZero , image); calculator->Update(); mitk::ClippedSurfaceBoundsCalculator::OutputType minMax = calculator->GetMinMaxSpatialDirectionZ(); MITK_TEST_CONDITION(minMax.first == minMax.second, "Check if plane is only on one slice"); MITK_TEST_CONDITION(minMax.first == 0 && minMax.second == 0, "Check if plane is on slice 0"); //Slice 3 origin[2] = 3; mitk::PlaneGeometry::Pointer planeOnSliceThree = mitk::PlaneGeometry::New(); planeOnSliceThree->InitializePlane(origin, normal); planeOnSliceThree->SetImageGeometry(false); calculator->SetInput( planeOnSliceThree , image); calculator->Update(); minMax = calculator->GetMinMaxSpatialDirectionZ(); MITK_TEST_CONDITION(minMax.first == minMax.second, "Check if plane is only on one slice"); MITK_TEST_CONDITION(minMax.first == 3 && minMax.second == 3, "Check if plane is on slice 3"); //Slice 17 origin[2] = 17; mitk::PlaneGeometry::Pointer planeOnSliceSeventeen = mitk::PlaneGeometry::New(); planeOnSliceSeventeen->InitializePlane(origin, normal); calculator->SetInput( planeOnSliceSeventeen , image); calculator->Update(); minMax = calculator->GetMinMaxSpatialDirectionZ(); MITK_TEST_CONDITION(minMax.first == minMax.second, "Check if plane is only on one slice"); MITK_TEST_CONDITION(minMax.first == 17 && minMax.second == 17, "Check if plane is on slice 17"); //Slice 20 origin[2] = 19; mitk::PlaneGeometry::Pointer planeOnSliceTwenty = mitk::PlaneGeometry::New(); planeOnSliceTwenty->InitializePlane(origin, normal); calculator->SetInput( planeOnSliceTwenty , image); calculator->Update(); minMax = calculator->GetMinMaxSpatialDirectionZ(); MITK_TEST_CONDITION(minMax.first == minMax.second, "Check if plane is only on one slice"); MITK_TEST_CONDITION(minMax.first == 19 && minMax.second == 19, "Check if plane is on slice 19"); delete calculator; } static void CheckPlanesInsideBoundingBox(mitk::Geometry3D::Pointer geometry3D) { //Check planes which are inside the bounding box mitk::ClippedSurfaceBoundsCalculator* calculator = new mitk::ClippedSurfaceBoundsCalculator(); mitk::Image::Pointer image = mitk::Image::New(); - image->SetGeometry(geometry3D); + image->Initialize( mitk::MakePixelType(), *(geometry3D.GetPointer()) ); //Check planes which are only on one slice: //Slice 0 mitk::Point3D origin; origin[0] = 511; // Set to 511.9 so that the intersection point is inside the bounding box origin[1] = 0; origin[2] = 0; mitk::Vector3D normal; mitk::FillVector3D(normal, 1, 0, 0); mitk::PlaneGeometry::Pointer planeSagittalOne = mitk::PlaneGeometry::New(); planeSagittalOne->InitializePlane(origin, normal); calculator->SetInput( planeSagittalOne , image); calculator->Update(); mitk::ClippedSurfaceBoundsCalculator::OutputType minMax = calculator->GetMinMaxSpatialDirectionZ(); MITK_TEST_CONDITION(minMax.first == 0 && minMax.second == 19, "Check if plane is from slice 0 to slice 19"); //Slice 3 origin[0] = 256; MITK_INFO << "Case1 origin: " << origin; mitk::PlaneGeometry::Pointer planeSagittalTwo = mitk::PlaneGeometry::New(); planeSagittalTwo->InitializePlane(origin, normal); MITK_INFO << "PlaneNormal: " << planeSagittalTwo->GetNormal(); MITK_INFO << "PlaneOrigin: " << planeSagittalTwo->GetOrigin(); calculator->SetInput( planeSagittalTwo , image); calculator->Update(); minMax = calculator->GetMinMaxSpatialDirectionZ(); MITK_INFO << "min: " << minMax.first << " max: " << minMax.second; MITK_TEST_CONDITION(minMax.first == 0 && minMax.second == 19, "Check if plane is from slice 0 to slice 19"); //Slice 17 origin[0] = 0; // Set to 0.1 so that the intersection point is inside the bounding box mitk::PlaneGeometry::Pointer planeOnSliceSeventeen = mitk::PlaneGeometry::New(); planeOnSliceSeventeen->InitializePlane(origin, normal); calculator->SetInput( planeOnSliceSeventeen , image); calculator->Update(); minMax = calculator->GetMinMaxSpatialDirectionZ(); MITK_TEST_CONDITION(minMax.first == 0 && minMax.second == 19, "Check if plane is from slice 0 to slice 19"); //Crooked planes: origin[0] = 0; origin[1] = 507; origin[2] = 0; normal[0] = 1; normal[1] = -1; normal[2] = 1; mitk::PlaneGeometry::Pointer planeCrookedOne = mitk::PlaneGeometry::New(); planeCrookedOne->InitializePlane(origin, normal); calculator->SetInput( planeCrookedOne , image); calculator->Update(); minMax = calculator->GetMinMaxSpatialDirectionZ(); MITK_INFO << "min: " << minMax.first << " max: " << minMax.second; MITK_TEST_CONDITION(minMax.first == 0 && minMax.second == 4, "Check if plane is from slice 0 to slice 4 with inclined plane"); origin[0] = 512; origin[1] = 0; origin[2] = 16; mitk::PlaneGeometry::Pointer planeCrookedTwo = mitk::PlaneGeometry::New(); planeCrookedTwo->InitializePlane(origin, normal); calculator->SetInput( planeCrookedTwo , image); calculator->Update(); minMax = calculator->GetMinMaxSpatialDirectionZ(); MITK_INFO << "min: " << minMax.first << " max: " << minMax.second; MITK_TEST_CONDITION(minMax.first == 17 && minMax.second == 19, "Check if plane is from slice 17 to slice 19 with inclined plane"); origin[0] = 511; origin[1] = 0; origin[2] = 0; normal[1] = 0; normal[2] = 0.04; mitk::PlaneGeometry::Pointer planeCrookedThree = mitk::PlaneGeometry::New(); planeCrookedThree->InitializePlane(origin, normal); calculator->SetInput( planeCrookedThree , image); calculator->Update(); minMax = calculator->GetMinMaxSpatialDirectionZ(); MITK_INFO << "min: " << minMax.first << " max: " << minMax.second; MITK_TEST_CONDITION(minMax.first == 0 && minMax.second == 19, "Check if plane is from slice 0 to slice 19 with inclined plane"); delete calculator; } static void CheckPlanesOutsideOfBoundingBox(mitk::Geometry3D::Pointer geometry3D) { //Check planes which are outside of the bounding box mitk::ClippedSurfaceBoundsCalculator* calculator = new mitk::ClippedSurfaceBoundsCalculator(); mitk::Image::Pointer image = mitk::Image::New(); - image->SetGeometry(geometry3D); + image->Initialize( mitk::MakePixelType(), *(geometry3D.GetPointer()) ); //In front of the bounding box mitk::Point3D origin; origin[0] = 511; origin[1] = 0; origin[2] = -5; mitk::Vector3D normal; mitk::FillVector3D(normal, 0, 0, 1); mitk::PlaneGeometry::Pointer planeInFront = mitk::PlaneGeometry::New(); planeInFront->InitializePlane(origin, normal); calculator->SetInput( planeInFront , image); calculator->Update(); mitk::ClippedSurfaceBoundsCalculator::OutputType minMax = calculator->GetMinMaxSpatialDirectionZ(); MITK_TEST_CONDITION(minMax.first == std::numeric_limits::max(), "Check if min value hasn't been set"); MITK_TEST_CONDITION(minMax.second == std::numeric_limits::min(), "Check if max value hasn't been set"); //Behind the bounding box origin[2] = 515; mitk::PlaneGeometry::Pointer planeBehind = mitk::PlaneGeometry::New(); planeBehind->InitializePlane(origin, normal); calculator->SetInput( planeBehind , image); calculator->Update(); minMax = calculator->GetMinMaxSpatialDirectionZ(); MITK_TEST_CONDITION(minMax.first == std::numeric_limits::max(), "Check if min value hasn't been set"); MITK_TEST_CONDITION(minMax.second == std::numeric_limits::min(), "Check if max value hasn't been set"); //Above origin[1] = 515; mitk::FillVector3D(normal, 0, 1, 0); mitk::PlaneGeometry::Pointer planeAbove = mitk::PlaneGeometry::New(); planeAbove->InitializePlane(origin, normal); calculator->SetInput( planeAbove , image); calculator->Update(); minMax = calculator->GetMinMaxSpatialDirectionZ(); MITK_TEST_CONDITION(minMax.first == std::numeric_limits::max(), "Check if min value hasn't been set"); MITK_TEST_CONDITION(minMax.second == std::numeric_limits::min(), "Check if max value hasn't been set"); //Below origin[1] = -5; mitk::PlaneGeometry::Pointer planeBelow = mitk::PlaneGeometry::New(); planeBelow->InitializePlane(origin, normal); calculator->SetInput( planeBelow , image); calculator->Update(); minMax = calculator->GetMinMaxSpatialDirectionZ(); MITK_TEST_CONDITION(minMax.first == std::numeric_limits::max(), "Check if min value hasn't been set"); MITK_TEST_CONDITION(minMax.second == std::numeric_limits::min(), "Check if max value hasn't been set"); //Left side origin[0] = -5; mitk::FillVector3D(normal, 1, 0, 0); mitk::PlaneGeometry::Pointer planeLeftSide = mitk::PlaneGeometry::New(); planeLeftSide->InitializePlane(origin, normal); calculator->SetInput( planeLeftSide , image); calculator->Update(); minMax = calculator->GetMinMaxSpatialDirectionZ(); MITK_TEST_CONDITION(minMax.first == std::numeric_limits::max(), "Check if min value hasn't been set"); MITK_TEST_CONDITION(minMax.second == std::numeric_limits::min(), "Check if max value hasn't been set"); //Right side origin[1] = 515; mitk::PlaneGeometry::Pointer planeRightSide = mitk::PlaneGeometry::New(); planeRightSide->InitializePlane(origin, normal); calculator->SetInput( planeRightSide , image); calculator->Update(); minMax = calculator->GetMinMaxSpatialDirectionZ(); MITK_TEST_CONDITION(minMax.first == std::numeric_limits::max(), "Check if min value hasn't been set"); MITK_TEST_CONDITION(minMax.second == std::numeric_limits::min(), "Check if max value hasn't been set"); delete calculator; } static void CheckIntersectionPointsOfTwoGeometry3D(mitk::Geometry3D::Pointer firstGeometry3D, mitk::Geometry3D::Pointer secondGeometry3D) { mitk::ClippedSurfaceBoundsCalculator* calculator = new mitk::ClippedSurfaceBoundsCalculator(); mitk::Image::Pointer firstImage = mitk::Image::New(); - firstImage->SetGeometry(firstGeometry3D); + firstImage->Initialize( mitk::MakePixelType(), *(firstGeometry3D.GetPointer()) ); calculator->SetInput( secondGeometry3D, firstImage); calculator->Update(); mitk::ClippedSurfaceBoundsCalculator::OutputType minMax = calculator->GetMinMaxSpatialDirectionZ(); minMax = calculator->GetMinMaxSpatialDirectionZ(); MITK_INFO << "min: " << minMax.first << " max: " << minMax.second; MITK_TEST_CONDITION(minMax.first == 0 && minMax.second == 19, "Check if plane is from slice 0 to slice 19"); } +static void CheckIntersectionWithPointCloud( mitk::Geometry3D::Pointer geometry3D ) +{ + //Check planes which are inside the bounding box + + mitk::Image::Pointer image = mitk::Image::New(); + image->Initialize( mitk::MakePixelType(), *(geometry3D.GetPointer()) ); + + { + mitk::Point3D pnt1, pnt2; + pnt1[0] = 3; + pnt1[1] = 5; + pnt1[2] = 3; + pnt2[0] = 8; + pnt2[1] = 3; + pnt2[2] = 8; + + mitk::ClippedSurfaceBoundsCalculator::PointListType pointlist; + pointlist.push_back( pnt1 ); + pointlist.push_back( pnt2 ); + + mitk::ClippedSurfaceBoundsCalculator calculator; + calculator.SetInput( pointlist, image ); + calculator.Update(); + + mitk::ClippedSurfaceBoundsCalculator::OutputType minMaxZ = calculator.GetMinMaxSpatialDirectionZ(); + MITK_TEST_CONDITION(minMaxZ.first == 3 && minMaxZ.second == 8, "Check if points span from slice 3 to slice 8 in axial"); + + mitk::ClippedSurfaceBoundsCalculator::OutputType minMaxX = calculator.GetMinMaxSpatialDirectionX(); + MITK_TEST_CONDITION(minMaxX.first == 3 && minMaxX.second == 5, "Check if points span from slice 3 to slice 5 in sagittal"); + } + + { + mitk::Point3D pnt1, pnt2; + pnt1.Fill( -3 ); + pnt2.Fill( 600 ); + + mitk::ClippedSurfaceBoundsCalculator::PointListType pointlist; + pointlist.push_back( pnt1 ); + pointlist.push_back( pnt2 ); + + mitk::ClippedSurfaceBoundsCalculator calculator; + calculator.SetInput( pointlist, image ); + calculator.Update(); + + mitk::ClippedSurfaceBoundsCalculator::OutputType minMaxZ = calculator.GetMinMaxSpatialDirectionZ(); + MITK_TEST_CONDITION(minMaxZ.first == 0 && minMaxZ.second == 19, "Check if points are correctly clipped to slice 0 and slice 19 in axial"); + + mitk::ClippedSurfaceBoundsCalculator::OutputType minMaxX = calculator.GetMinMaxSpatialDirectionX(); + MITK_TEST_CONDITION(minMaxX.first == 0 && minMaxX.second == 511, "Check if points are correctly clipped to slice 0 and slice 511 in sagittal"); + } + + + +} + + int mitkClippedSurfaceBoundsCalculatorTest(int, char* []) { // always start with this! MITK_TEST_BEGIN("ClippedSurfaceBoundsCalculator"); /** The class mitkClippedSurfaceBoundsCalculator calculates the intersection points of a PlaneGeometry and a Geometry3D. * This unittest checks if the correct min and max values for the three spatial directions (x, y, z) * are calculated. To test this we define artifical PlaneGeometries and Geometry3Ds and test different * scenarios: * * 1. planes which are inside the bounding box of a 3D geometry but only on one slice * 2. planes which are outside of the bounding box * 3. planes which are inside the bounding box but over more than one slice * * Note: Currently rotated geometries are not tested! */ /********************* Define Geometry3D ***********************/ //Define origin: mitk::Point3D origin; origin[0] = 511; origin[1] = 0; origin[2] = 0; //Define normal: mitk::Vector3D normal; mitk::FillVector3D(normal, 0, 0, 1); //Initialize PlaneGeometry: mitk::PlaneGeometry::Pointer planeGeometry = mitk::PlaneGeometry::New(); planeGeometry->InitializePlane(origin, normal); //Set Bounds: mitk::BoundingBox::BoundsArrayType bounds = planeGeometry->GetBounds(); bounds[0] = 0; bounds[1] = 512; bounds[2] = 0; bounds[3] = 512; bounds[4] = 0; bounds[5] = 1; planeGeometry->SetBounds(bounds); //Initialize SlicedGeometry3D: mitk::SlicedGeometry3D::Pointer slicedGeometry3D = mitk::SlicedGeometry3D::New(); slicedGeometry3D->InitializeEvenlySpaced(dynamic_cast(planeGeometry.GetPointer()), 20); mitk::Geometry3D::Pointer geometry3D = dynamic_cast< mitk::Geometry3D* > ( slicedGeometry3D.GetPointer() ); geometry3D->SetImageGeometry(true); //Define origin for second Geometry3D; mitk::Point3D origin2; origin2[0] = 511; origin2[1] = 60; origin2[2] = 0; //Define normal: mitk::Vector3D normal2; mitk::FillVector3D(normal2, 0, 1, 0); //Initialize PlaneGeometry: mitk::PlaneGeometry::Pointer planeGeometry2 = mitk::PlaneGeometry::New(); planeGeometry2->InitializePlane(origin2, normal2); //Initialize SlicedGeometry3D: mitk::SlicedGeometry3D::Pointer secondSlicedGeometry3D = mitk::SlicedGeometry3D::New(); secondSlicedGeometry3D->InitializeEvenlySpaced(dynamic_cast(planeGeometry2.GetPointer()), 20); mitk::Geometry3D::Pointer secondGeometry3D = dynamic_cast< mitk::Geometry3D* > ( secondSlicedGeometry3D.GetPointer() ); secondGeometry3D->SetImageGeometry(true); /***************************************************************/ CheckPlanesInsideBoundingBoxOnlyOnOneSlice(geometry3D); CheckPlanesOutsideOfBoundingBox(geometry3D); CheckPlanesInsideBoundingBox(geometry3D); CheckIntersectionPointsOfTwoGeometry3D(geometry3D, secondGeometry3D); + CheckIntersectionWithPointCloud( geometry3D ); /** ToDo: * test also rotated 3D geometry! */ MITK_TEST_END(); }