diff --git a/Core/Code/Algorithms/mitkClippedSurfaceBoundsCalculator.cpp b/Core/Code/Algorithms/mitkClippedSurfaceBoundsCalculator.cpp index 1b8f763fab..ae8c58e11c 100644 --- a/Core/Code/Algorithms/mitkClippedSurfaceBoundsCalculator.cpp +++ b/Core/Code/Algorithms/mitkClippedSurfaceBoundsCalculator.cpp @@ -1,211 +1,276 @@ #include "mitkClippedSurfaceBoundsCalculator.h" #include "mitkLine.h" -#define ROUND_P(x) ((int)((x)+0.5)) +#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) + 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->SetInput(geometry, image); } +mitk::ClippedSurfaceBoundsCalculator::ClippedSurfaceBoundsCalculator( + const mitk::Geometry3D* geometry, + 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->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 + } +} + +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 } } 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() { - // SEE HEADER DOCUMENTATION for explanation - - typedef std::vector< std::pair > EdgesVector; - 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))); + } +} + +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 - m_PlaneGeometry->IntersectionPoint(line, intersectionWorldPoint); + geometry->IntersectionPoint(line, intersectionWorldPoint); double t = -1.0; - bool isIntersectionPointOnLine; - isIntersectionPointOnLine = m_PlaneGeometry->IntersectionPointParam(line, t); + + 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(0.0 <= t && t <= 1.0 && isIntersectionPointOnLine) + 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 } } } } } diff --git a/Core/Code/Algorithms/mitkClippedSurfaceBoundsCalculator.h b/Core/Code/Algorithms/mitkClippedSurfaceBoundsCalculator.h index 89201b1520..b5eca2af94 100644 --- a/Core/Code/Algorithms/mitkClippedSurfaceBoundsCalculator.h +++ b/Core/Code/Algorithms/mitkClippedSurfaceBoundsCalculator.h @@ -1,102 +1,107 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Module: $RCSfile: mitkPropertyManager.cpp,v $ Language: C++ Date: $Date: 2005/06/28 12:37:25 $ Version: $Revision: 1.12 $ 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 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: ClippedSurfaceBoundsCalculator(const mitk::PlaneGeometry* geometry = NULL, mitk::Image::Pointer image = NULL); + ClippedSurfaceBoundsCalculator(const mitk::Geometry3D* geometry, + mitk::Image::Pointer image); virtual ~ClippedSurfaceBoundsCalculator(); void SetInput(const mitk::PlaneGeometry* geometry, mitk::Image* image); + void SetInput(const mitk::Geometry3D *geometry, 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); mitk::PlaneGeometry::ConstPointer m_PlaneGeometry; + mitk::Geometry3D::ConstPointer m_Geometry3D; mitk::Image::Pointer m_Image; std::vector< OutputType > m_MinMaxOutput; }; } //namespace mitk #endif diff --git a/Core/Code/Testing/mitkClippedSurfaceBoundsCalculatorTest.cpp b/Core/Code/Testing/mitkClippedSurfaceBoundsCalculatorTest.cpp index b363685972..b9a0d7b22b 100644 --- a/Core/Code/Testing/mitkClippedSurfaceBoundsCalculatorTest.cpp +++ b/Core/Code/Testing/mitkClippedSurfaceBoundsCalculatorTest.cpp @@ -1,388 +1,428 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date: 2008-02-25 17:27:17 +0100 (Mo, 25 Feb 2008) $ Version: $Revision: 7837 $ 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 "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); //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); //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); //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); + + 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"); +} + + 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); + /** ToDo: * test also rotated 3D geometry! */ MITK_TEST_END(); }