Index: Core/Code/Controllers/mitkSliceNavigationController.cpp =================================================================== --- Core/Code/Controllers/mitkSliceNavigationController.cpp (revision 24914) +++ Core/Code/Controllers/mitkSliceNavigationController.cpp (working copy) @@ -664,13 +664,13 @@ std::stringstream stream; - // get the position and gray value from the image and build up status bar text - mitk::Point3D p; + // get the position and gray value from the image and build up status bar text if(image3D.IsNotNull()) { + Index3D p; image3D->GetGeometry()->WorldToIndex(worldposition, p); stream<<"Position: <"< mm"; - stream<<"; Index: <"< "; + stream<<"; Index: <"< "; stream<<"; Time: " << baseRenderer->GetTime() << " ms; Pixelvalue: "<GetPixelValueByIndex(p, baseRenderer->GetTimeStep())<<" "; } else Index: Core/Code/DataManagement/mitkGeometry3D.cpp =================================================================== --- Core/Code/DataManagement/mitkGeometry3D.cpp (revision 24914) +++ Core/Code/DataManagement/mitkGeometry3D.cpp (working copy) @@ -142,15 +142,7 @@ void mitk::Geometry3D::IndexToWorld(const mitk::Point3D &pt_units, mitk::Point3D &pt_mm) const { - mitk::Point3D temp = pt_units; - if(m_ImageGeometry) - { - for (int i = 0; i < 3; i++) - { - temp[i] -= 0.5; - } - } - pt_mm = m_IndexToWorldTransform->TransformPoint(temp); + pt_mm = m_IndexToWorldTransform->TransformPoint(pt_units); } void mitk::Geometry3D::WorldToIndex(const mitk::Point3D &atPt3d_mm, const mitk::Vector3D &vec_mm, mitk::Vector3D &vec_units) const @@ -365,15 +357,6 @@ out[i] += inverse[i][j]*temp[j]; } } - - // Image index coordinates are edge-based - if(m_ImageGeometry) - { - for (i = 0; i < 3; i++) - { - out[i] += 0.5; - } - } } void mitk::Geometry3D::BackTransform(const mitk::Point3D &/*at*/, const mitk::Vector3D &in, mitk::Vector3D& out) const Index: Core/Code/DataManagement/mitkGeometry3D.h =================================================================== --- Core/Code/DataManagement/mitkGeometry3D.h (revision 24914) +++ Core/Code/DataManagement/mitkGeometry3D.h (working copy) @@ -58,16 +58,20 @@ //## GetBoundingBox() //## \li a transform to convert intrinsic coordinates into a //## world-coordinate system with coordinates in millimeters -//## and milliseconds (floating point values), to be accessed -//## by GetIndexToWorldTransform() +//## and milliseconds (all are floating point values), to +//## be accessed by GetIndexToWorldTransform() //## \li a life span, i.e. a bounding box in time in ms (with //## start and end time), to be accessed by GetTimeBounds(). -//## The default is The default is minus infinity to plus infinity. +//## The default is minus infinity to plus infinity. //## //## Geometry3D and its sub-classes allow converting between //## intrinsic coordinates (called index or unit coordinates) -//## and word-coordinates (called world or mm coordinates), +//## and world-coordinates (called world or mm coordinates), //## e.g. WorldToIndex. +//## In case you need integer index coordinates, provide an +//## mitk::Index3D (or itk::Index) as target variable to +//## WorldToIndex, otherwise you will get a continuous index +//## (floating point values). //## //## An important sub-class is SlicedGeometry3D, which descibes //## data objects consisting of slices, e.g., objects of type Image. @@ -248,7 +252,7 @@ //##Documentation //## @brief Translate the origin by a vector //## - virtual void Translate(const Vector3D & vector); + virtual void Translate(const Vector3D& vector); //##Documentation //## @brief Set the transform to identity @@ -291,32 +295,36 @@ } //##Documentation - //## @brief Convert world coordinates (in mm) of a \em point to index coordinates (in units) - void WorldToIndex(const mitk::Point3D &pt_mm, mitk::Point3D &pt_units) const; + //## @brief Convert world coordinates (in mm) of a \em point to (continuous!) index coordinates (in units) + //## \warning If you need integer index coordinates (e.g., for accessing a pixel in an image), + //## use WorldToIndex(const mitk::Point3D& pt_mm, itk::Index &index). + void WorldToIndex(const mitk::Point3D& pt_mm, mitk::Point3D& pt_units) const; //##Documentation //## @brief Convert index coordinates (in units) of a \em point to world coordinates (in mm) - void IndexToWorld(const mitk::Point3D &pt_units, mitk::Point3D &pt_mm) const; + void IndexToWorld(const mitk::Point3D& pt_units, mitk::Point3D& pt_mm) const; //##Documentation //## @brief Convert world coordinates (in mm) of a \em vector - //## \a vec_mm (at the point \a atPt3d_mm) to index coordinates (in units) - void WorldToIndex(const mitk::Point3D &atPt3d_mm, const mitk::Vector3D &vec_mm, mitk::Vector3D &vec_units) const; + //## \a vec_mm (at the point \a atPt3d_mm) to (continuous!) index coordinates (in units) + //## \warning If you need integer index coordinates (e.g., for accessing a pixel in an image), + //## use WorldToIndex(const mitk::Point3D& pt_mm, itk::Index &index). + void WorldToIndex(const mitk::Point3D& atPt3d_mm, const mitk::Vector3D& vec_mm, mitk::Vector3D& vec_units) const; //##Documentation //## @brief Convert index coordinates (in units) of a \em vector //## \a vec_units (at the point \a atPt3d_units) to world coordinates (in mm) - void IndexToWorld(const mitk::Point3D &atPt3d_units, const mitk::Vector3D &vec_units, mitk::Vector3D &vec_mm) const; + void IndexToWorld(const mitk::Point3D& atPt3d_units, const mitk::Vector3D& vec_units, mitk::Vector3D& vec_mm) const; //##Documentation //## @brief Convert world coordinates (in mm) of a \em point to index coordinates (in units). //## This method rounds to integer indices! template - void WorldToIndex(const mitk::Point3D &pt_mm, itk::Index &index) const + void WorldToIndex(const mitk::Point3D& pt_mm, itk::Index &index) const { typedef itk::Index IndexType; mitk::Point3D pt_units; - WorldToIndex(pt_mm, pt_units); + this->WorldToIndex(pt_mm, pt_units); int i, dim=index.GetIndexDimension(); if(dim>3) { @@ -324,17 +332,18 @@ dim=3; } for(i=0;i( pt_units[i] ); } } //##Documentation - //## @brief Convert world coordinates (in mm) of a \em point to + //## @brief Deprecated for use with ITK version 3.10 or newer. + //## Convert world coordinates (in mm) of a \em point to //## ITK physical coordinates (in mm, but without a possible rotation) //## //## This method is useful if you have want to access an mitk::Image - //## via an itk::Image. ITK does not support rotated (tilted) images, - //## i.e., ITK images are always parallel to the coordinate axes. + //## via an itk::Image. ITK v3.8 and older did not support rotated (tilted) + //## images, i.e., ITK images are always parallel to the coordinate axes. //## When accessing a (possibly rotated) mitk::Image via an itk::Image //## the rotational part of the transformation in the Geometry3D is //## simply discarded; in other word: only the origin and spacing is @@ -344,32 +353,41 @@ //## can be used with the ITK image as a ITK physical coordinate //## (excluding the rotation). template - void WorldToItkPhysicalPoint(const mitk::Point3D &pt_mm, + void WorldToItkPhysicalPoint(const mitk::Point3D& pt_mm, itk::Point& itkPhysicalPoint) const { - mitk::Point3D index; - WorldToIndex(pt_mm, index); - for (unsigned int i = 0 ; i < 3 ; i++) - { - itkPhysicalPoint[i] = static_cast( this->m_Spacing[i] * index[i] + this->m_Origin[i] ); - } + #if ((ITK_VERSION_MAJOR > 3) || (ITK_VERSION_MAJOR == 3 && ITK_VERSION_MINOR > 8)) + mitk::vtk2itk(pt_mm, itkPhysicalPoint); + #else + mitk::Point3D index; + WorldToIndex(pt_mm, index); + for (unsigned int i = 0 ; i < 3 ; i++) + { + itkPhysicalPoint[i] = static_cast( this->m_Spacing[i] * index[i] + this->m_Origin[i] ); + } + #endif } //##Documentation - //## @brief Convert ITK physical coordinates of a \em point (in mm, + //## @brief Deprecated for use with ITK version 3.10 or newer. + //## Convert ITK physical coordinates of a \em point (in mm, //## but without a rotation) into MITK world coordinates (in mm) //## - //## For more information, see WorldToItkPhysicalPoint + //## For more information, see WorldToItkPhysicalPoint. template void ItkPhysicalPointToWorld(const itk::Point& itkPhysicalPoint, - mitk::Point3D &pt_mm) const + mitk::Point3D& pt_mm) const { - mitk::Point3D index; - for (unsigned int i = 0 ; i < 3 ; i++) - { - index[i] = static_cast( (itkPhysicalPoint[i]- this->m_Origin[i]) / this->m_Spacing[i] ); - } - IndexToWorld(index, pt_mm); + #if ((ITK_VERSION_MAJOR > 3) || (ITK_VERSION_MAJOR == 3 && ITK_VERSION_MINOR > 8)) + mitk::vtk2itk(itkPhysicalPoint, pt_mm); + #else + mitk::Point3D index; + for (unsigned int i = 0 ; i < 3 ; i++) + { + index[i] = static_cast( (itkPhysicalPoint[i]- this->m_Origin[i]) / this->m_Spacing[i] ); + } + IndexToWorld(index, pt_mm); + #endif } //##Documentation @@ -571,8 +589,8 @@ virtual void PrintSelf(std::ostream& os, itk::Indent indent) const; - virtual void BackTransform(const mitk::Point3D &in, mitk::Point3D& out) const; - virtual void BackTransform(const mitk::Point3D &at, const mitk::Vector3D &in, mitk::Vector3D& out) const; + virtual void BackTransform(const mitk::Point3D& in, mitk::Point3D& out) const; + virtual void BackTransform(const mitk::Point3D& at, const mitk::Vector3D& in, mitk::Vector3D& out) const; //##Documentation //## @brief Set the parametric bounds Index: Core/Code/DataManagement/mitkImage.cpp =================================================================== --- Core/Code/DataManagement/mitkImage.cpp (revision 24914) +++ Core/Code/DataManagement/mitkImage.cpp (working copy) @@ -99,7 +99,30 @@ return m_CompleteData->GetData(); } -double mitk::Image::GetPixelValueByIndex(const mitk::Point3D &position, unsigned int timestep) +template +void AccessPixel(mitkIpPicDescriptor* pic, const mitk::Index3D& p, double& value, int timestep) +{ + if ( (p[0]>=0 && p[1] >=0 && p[2]>=0 && timestep>=0) && (unsigned int)p[0] < pic->n[0] && (unsigned int)p[1] < pic->n[1] && (unsigned int)p[2] < pic->n[2] && (unsigned int)timestep < pic->n[3] ) + { + if(pic->bpe!=24) + { + value = (double) (((T*) pic->data)[ p[0] + p[1]*pic->n[0] + p[2]*pic->n[0]*pic->n[1] + timestep*pic->n[0]*pic->n[1]*pic->n[2] ]); + } + else + { + double returnvalue = (((T*) pic->data)[p[0]*3 + 0 + p[1]*pic->n[0]*3 + p[2]*pic->n[0]*pic->n[1]*3 + timestep*pic->n[0]*pic->n[1]*pic->n[2]*3 ]); + returnvalue += (((T*) pic->data)[p[0]*3 + 1 + p[1]*pic->n[0]*3 + p[2]*pic->n[0]*pic->n[1]*3 + timestep*pic->n[0]*pic->n[1]*pic->n[2]*3]); + returnvalue += (((T*) pic->data)[p[0]*3 + 2 + p[1]*pic->n[0]*3 + p[2]*pic->n[0]*pic->n[1]*3 + timestep*pic->n[0]*pic->n[1]*pic->n[2]*3]); + value = returnvalue; + } + } + else + { + value = 0; + } +}; + +double mitk::Image::GetPixelValueByIndex(const mitk::Index3D &position, unsigned int timestep) { mitkIpPicDescriptor* pic = this->GetPic(); double value = 0; @@ -111,7 +134,7 @@ return value; } -double mitk::Image::GetPixelValueByWorldCoordinate(const mitk::Point3D &position, unsigned int timestep) +double mitk::Image::GetPixelValueByWorldCoordinate(const mitk::Point3D& position, unsigned int timestep) { mitkIpPicDescriptor* pic = this->GetPic(); double value = 0; @@ -120,13 +143,10 @@ timestep = this->GetTimeSteps(); } - itk::Index<3> itkIndex; + Index3D itkIndex; this->GetGeometry()->WorldToIndex(position,itkIndex); - Point3D mitkPointIndex; - mitkPointIndex[0]=itkIndex[0]; - mitkPointIndex[1]=itkIndex[1]; - mitkPointIndex[2]=itkIndex[2]; - mitkIpPicTypeMultiplex3(AccessPixel, pic, mitkPointIndex, value, timestep); + mitkIpPicTypeMultiplex3(AccessPixel, pic, itkIndex, value, timestep); + return value; } Index: Core/Code/DataManagement/mitkImage.h =================================================================== --- Core/Code/DataManagement/mitkImage.h (revision 24914) +++ Core/Code/DataManagement/mitkImage.h (working copy) @@ -54,9 +54,9 @@ //## For importing ITK images use of mitk::ITKImageImport is recommended, see //## \ref Adaptor. //## -//## Converting coordinates from the ITK physical coordinate system -//## (which does not support rotated images) to the MITK world coordinate -//## system should be performed via the Geometry3D of the Image, see +//## For ITK v3.8 and older: Converting coordinates from the ITK physical +//## coordinate system (which does not support rotated images) to the MITK world +//## coordinate system should be performed via the Geometry3D of the Image, see //## Geometry3D::WorldToItkPhysicalPoint. //## @ingroup Data class MITK_CORE_EXPORT Image : public SlicedData @@ -113,12 +113,12 @@ //## @brief Get the pixel value at one specific position. //## //## The pixel type is always being converted to double. - double GetPixelValueByIndex(const mitk::Point3D &position, unsigned int timestep = 0); + double GetPixelValueByIndex(const mitk::Index3D& position, unsigned int timestep = 0); //## @brief Get the pixel value at one specific world position. //## //## The pixel type is always being converted to double. - double GetPixelValueByWorldCoordinate(const mitk::Point3D &position, unsigned int timestep = 0); + double GetPixelValueByWorldCoordinate(const mitk::Point3D& position, unsigned int timestep = 0); //##Documentation //## @brief Get a volume at a specific time @a t of channel @a n as a vtkImageData. @@ -540,8 +540,6 @@ friend void _ComputeExtremaInItkImage(ItkImageType* itkImage, mitk::Image * mitkImage, int t); protected: - template - void AccessPixel(mitkIpPicDescriptor* pic, mitk::Point3D p, double& value, int timestep = 0); int GetSliceIndex(int s = 0, int t = 0, int n = 0) const; @@ -600,34 +598,6 @@ }; -//## @brief Get the pixel value at one specific position. -//## -//## The pixel type is always being converted to double. -template -void Image::AccessPixel(mitkIpPicDescriptor* pic, mitk::Point3D p, double& value, int timestep) -{ - itk::Point pi; - mitk::itk2vtk(p, pi); - if ( (pi[0]>=0 && pi[1] >=0 && pi[2]>=0 && timestep>=0) && (unsigned int)pi[0] < pic->n[0] && (unsigned int)pi[1] < pic->n[1] && (unsigned int)pi[2] < pic->n[2] && (unsigned int)timestep < pic->n[3] ) - { - if(pic->bpe!=24) - { - value = (double) (((T*) pic->data)[ pi[0] + pi[1]*pic->n[0] + pi[2]*pic->n[0]*pic->n[1] + timestep*pic->n[0]*pic->n[1]*pic->n[2] ]); - } - else - { - double returnvalue = (((T*) pic->data)[pi[0]*3 + 0 + pi[1]*pic->n[0]*3 + pi[2]*pic->n[0]*pic->n[1]*3 + timestep*pic->n[0]*pic->n[1]*pic->n[2]*3 ]); - returnvalue += (((T*) pic->data)[pi[0]*3 + 1 + pi[1]*pic->n[0]*3 + pi[2]*pic->n[0]*pic->n[1]*3 + timestep*pic->n[0]*pic->n[1]*pic->n[2]*3]); - returnvalue += (((T*) pic->data)[pi[0]*3 + 2 + pi[1]*pic->n[0]*3 + pi[2]*pic->n[0]*pic->n[1]*3 + timestep*pic->n[0]*pic->n[1]*pic->n[2]*3]); - value = returnvalue; - } - } - else - { - value = 0; - } -}; - //##Documentation //## @brief Cast an itk::Image (with a specific type) to an mitk::Image. //## Index: Core/Code/DataManagement/mitkPlaneGeometry.cpp =================================================================== --- Core/Code/DataManagement/mitkPlaneGeometry.cpp (revision 24914) +++ Core/Code/DataManagement/mitkPlaneGeometry.cpp (working copy) @@ -95,12 +95,6 @@ { pt_mm[0]=m_ScaleFactorMMPerUnitX*pt_units[0]; pt_mm[1]=m_ScaleFactorMMPerUnitY*pt_units[1]; - - if (m_ImageGeometry) - { - pt_mm[0]-=0.5; - pt_mm[1]-=0.5; - } } @@ -108,12 +102,6 @@ { pt_units[0]=pt_mm[0]*(1.0/m_ScaleFactorMMPerUnitX); pt_units[1]=pt_mm[1]*(1.0/m_ScaleFactorMMPerUnitY); - - if (m_ImageGeometry) - { - pt_units[0]+=0.5; - pt_units[1]+=0.5; - } } Index: Core/Code/DataManagement/mitkVector.h =================================================================== --- Core/Code/DataManagement/mitkVector.h (revision 24914) +++ Core/Code/DataManagement/mitkVector.h (working copy) @@ -22,6 +22,7 @@ #include #include #include +#include #include #include #include @@ -48,6 +49,8 @@ typedef itk::Point Point4I; typedef itk::Vector Vector2D; typedef itk::Vector Vector3D; +typedef itk::Index<3> Index3D; +typedef itk::ContinuousIndex ContinuousIndex3D; typedef vnl_quaternion Quaternion; //##Documentation Index: Core/Code/Testing/mitkGeometry3DTest.cpp =================================================================== --- Core/Code/Testing/mitkGeometry3DTest.cpp (revision 24914) +++ Core/Code/Testing/mitkGeometry3DTest.cpp (working copy) @@ -22,9 +22,9 @@ #include "mitkRotationOperation.h" #include "mitkInteractionConst.h" +#include #include "mitkTestingMacros.h" - #include bool testGetAxisVectorVariants(mitk::Geometry3D* geometry) @@ -65,8 +65,128 @@ return true; } -int mitkGeometry3DTest(int /*argc*/, char* /*argv*/[]) +// a part of the test requires axis-parallel coordinates +int testIndexAndWorldConsistency(mitk::Geometry3D* geometry3d) { + MITK_TEST_OUTPUT( << "Testing consistency of index and world coordinate systems: "); + mitk::Point3D origin = geometry3d->GetOrigin(); + mitk::Point3D dummy; + + MITK_TEST_OUTPUT( << " Testing index->world->index conversion consistency"); + geometry3d->WorldToIndex(origin, dummy); + geometry3d->IndexToWorld(dummy, dummy); + MITK_TEST_CONDITION_REQUIRED(dummy == origin, ""); + + MITK_TEST_OUTPUT( << " Testing WorldToIndex(origin, mitk::Point3D)==(0,0,0)"); + mitk::Point3D globalOrigin; + mitk::FillVector3D(globalOrigin, 0,0,0); + + mitk::Point3D originContinuousIndex; + geometry3d->WorldToIndex(origin, originContinuousIndex); + MITK_TEST_CONDITION_REQUIRED(originContinuousIndex == globalOrigin, ""); + + MITK_TEST_OUTPUT( << " Testing WorldToIndex(origin, itk::Index)==(0,0,0)"); + itk::Index<3> itkindex; + geometry3d->WorldToIndex(origin, itkindex); + itk::Index<3> globalOriginIndex; + mitk::vtk2itk(globalOrigin, globalOriginIndex); + MITK_TEST_CONDITION_REQUIRED(itkindex == globalOriginIndex, ""); + + MITK_TEST_OUTPUT( << " Testing WorldToIndex(origin-0.5*spacing, itk::Index)==(0,0,0)"); + mitk::Vector3D halfSpacingStep = geometry3d->GetSpacing()*0.5; + mitk::Matrix3D rotation; + mitk::Point3D originOffCenter = origin-halfSpacingStep; + geometry3d->WorldToIndex(originOffCenter, itkindex); + MITK_TEST_CONDITION_REQUIRED(itkindex == globalOriginIndex, ""); + + MITK_TEST_OUTPUT( << " Testing WorldToIndex(origin+0.5*spacing-eps, itk::Index)==(0,0,0)"); + originOffCenter = origin+halfSpacingStep; + originOffCenter -= 0.0001; + geometry3d->WorldToIndex( originOffCenter, itkindex); + MITK_TEST_CONDITION_REQUIRED(itkindex == globalOriginIndex, ""); + + MITK_TEST_OUTPUT( << " Testing WorldToIndex(origin+0.5*spacing, itk::Index)==(1,1,1)"); + originOffCenter = origin+halfSpacingStep; + itk::Index<3> global111; + mitk::FillVector3D(global111, 1,1,1); + geometry3d->WorldToIndex( originOffCenter, itkindex); + MITK_TEST_CONDITION_REQUIRED(itkindex == global111, ""); + + MITK_TEST_OUTPUT( << " Testing WorldToIndex(GetCenter())==BoundingBox.GetCenter: "); + mitk::Point3D center = geometry3d->GetCenter(); + mitk::Point3D centerContIndex; + geometry3d->WorldToIndex(center, centerContIndex); + mitk::BoundingBox::ConstPointer boundingBox = geometry3d->GetBoundingBox(); + mitk::BoundingBox::PointType centerBounds = boundingBox->GetCenter(); + MITK_TEST_CONDITION_REQUIRED(mitk::Equal(centerContIndex,centerBounds), ""); + + MITK_TEST_OUTPUT( << " Testing GetCenter()==IndexToWorld(BoundingBox.GetCenter): "); + center = geometry3d->GetCenter(); + mitk::Point3D centerBoundsInWorldCoords; + geometry3d->IndexToWorld(centerBounds, centerBoundsInWorldCoords); + MITK_TEST_CONDITION_REQUIRED(mitk::Equal(center,centerBoundsInWorldCoords), ""); + + return EXIT_SUCCESS; +} + +#include + +int testItkImageIsCenterBased() +{ + MITK_TEST_OUTPUT(<< "Testing whether itk::Image coordinates are center-based."); + typedef itk::Image ItkIntImage3D; + ItkIntImage3D::Pointer itkintimage = ItkIntImage3D::New(); + ItkIntImage3D::SizeType size; + size.Fill(10); + mitk::Point3D origin; + mitk::FillVector3D(origin, 2,3,7); + itkintimage->Initialize(); + itkintimage->SetRegions(size); + itkintimage->SetOrigin(origin); + std::cout<<"[PASSED]"< originContinuousIndex; + itkintimage->TransformPhysicalPointToContinuousIndex(origin, originContinuousIndex); + MITK_TEST_CONDITION_REQUIRED(originContinuousIndex == globalOrigin, ""); + + MITK_TEST_OUTPUT( << " Testing itk::Image::TransformPhysicalPointToIndex(origin)==(0,0,0)"); + itk::Index<3> itkindex; + itkintimage->TransformPhysicalPointToIndex(origin, itkindex); + itk::Index<3> globalOriginIndex; + mitk::vtk2itk(globalOrigin, globalOriginIndex); + MITK_TEST_CONDITION_REQUIRED(itkindex == globalOriginIndex, ""); + + MITK_TEST_OUTPUT( << " Testing itk::Image::TransformPhysicalPointToIndex(origin-0.5*spacing)==(0,0,0)"); + mitk::Vector3D halfSpacingStep = itkintimage->GetSpacing()*0.5; + mitk::Matrix3D rotation; + mitk::Point3D originOffCenter = origin-halfSpacingStep; + itkintimage->TransformPhysicalPointToIndex(originOffCenter, itkindex); + MITK_TEST_CONDITION_REQUIRED(itkindex == globalOriginIndex, ""); + + MITK_TEST_OUTPUT( << " Testing itk::Image::TransformPhysicalPointToIndex(origin+0.5*spacing-eps, itk::Index)==(0,0,0)"); + originOffCenter = origin+halfSpacingStep; + originOffCenter -= 0.0001; + itkintimage->TransformPhysicalPointToIndex( originOffCenter, itkindex); + MITK_TEST_CONDITION_REQUIRED(itkindex == globalOriginIndex, ""); + + MITK_TEST_OUTPUT( << " Testing itk::Image::TransformPhysicalPointToIndex(origin+0.5*spacing, itk::Index)==(1,1,1)"); + originOffCenter = origin+halfSpacingStep; + itk::Index<3> global111; + mitk::FillVector3D(global111, 1,1,1); + itkintimage->TransformPhysicalPointToIndex( originOffCenter, itkindex); + MITK_TEST_CONDITION_REQUIRED(itkindex == global111, ""); + + MITK_TEST_OUTPUT( << "=> Yes, itk::Image coordinates are center-based."); + + return EXIT_SUCCESS; +} + +int testGeometry3D(bool imageGeometry) +{ float bounds[ ] = {-10.0, 17.0, -12.0, 188.0, 13.0, 211.0}; mitk::Geometry3D::Pointer geometry3d = mitk::Geometry3D::New(); @@ -75,6 +195,10 @@ geometry3d->Initialize(); std::cout<<"[PASSED]"<SetImageGeometry(imageGeometry); + std::cout<<"[PASSED]"<SetFloatBounds(bounds); std::cout<<"[PASSED]"<SetSpacing(newspacing); + const mitk::Vector3D& spacing3 = geometry3d->GetSpacing(); + if( mitk::Equal(spacing3, newspacing) == false ) + { + std::cout<<"[FAILED]"<ExecuteOperation(op); - // Todo: Find a meaningful way to test rotation success + MITK_TEST_OUTPUT( << " Testing mitk::GetRotation() and success of rotation"); + mitk::Matrix3D rotation; + mitk::GetRotation(geometry3d, rotation); + mitk::Vector3D voxelStep=rotation*newspacing; + mitk::Vector3D voxelStepIndex; + geometry3d->WorldToIndex(geometry3d->GetOrigin(), voxelStep, voxelStepIndex); + mitk::Vector3D expectedVoxelStepIndex; + expectedVoxelStepIndex.Fill(1); + MITK_TEST_CONDITION_REQUIRED(mitk::Equal(voxelStepIndex,expectedVoxelStepIndex), ""); delete op; std::cout<<"[PASSED]"<GetImageGeometry() == imageGeometry, ""); - mitk::Point3D point = geometry3d->GetOrigin(); - mitk::Point3D dummy; + return EXIT_SUCCESS; +} - geometry3d->WorldToIndex(point, dummy); - geometry3d->IndexToWorld(dummy, dummy); - MITK_TEST_CONDITION_REQUIRED(dummy == point, ""); +int mitkGeometry3DTest(int /*argc*/, char* /*argv*/[]) +{ + MITK_TEST_BEGIN(mitkGeometry3DTest); - std::cout<<"[TEST DONE]"< +#include +#include + +#include + #include #include #include @@ -342,17 +346,24 @@ } std::cout<<"[PASSED]"<Initialize(*imgMem->GetPixelType().GetTypeId(), *imgMem->GetGeometry(), 2 /* #channels */, 0 /*tDim*/ ); vecImg->SetImportChannel(imgMem->GetData(), 0, mitk::Image::CopyMemory ); vecImg->SetImportChannel(imgMem->GetData(), 1, mitk::Image::CopyMemory ); - if( !vecImg->IsValidSlice(0,0,1)) - { - std::cout<<"[FAILED]"<IsValidSlice(0,0,1) , ""); + + MITK_TEST_OUTPUT(<< " Testing whether CopyMemory worked"); + MITK_TEST_CONDITION_REQUIRED(imgMem->GetData() != vecImg->GetData(), ""); + + MITK_TEST_OUTPUT(<< " Testing destruction after SetImportChannel"); + vecImg = NULL; + std::cout<<"[PASSED]"<SetFileName("C:/MITK-QT4/mitk_source3m/mitk/Core/Code/Testing/Data/brain.mhd"); + nodeReader->SetFileName("V:/windows/source/mitk/mitk/Core/Code/Testing/Data/brain.mhd"); nodeReader->Update(); node = nodeReader->GetOutput(); } @@ -423,28 +434,49 @@ } mitk::Image::Pointer image = dynamic_cast(node->GetData()); - - mitk::Point3D point; - + // test by index coordinates - mitk::FillVector3D(point, 55, 39, 50); - double val = image->GetPixelValueByIndex(point); - - // there is always a small rounding error - if(abs(112.225 - val) > 1.0){ - std::cout << "[FAILED]"; - return EXIT_FAILURE; - } + mitk::Index3D index; + mitk::FillVector3D(index, 55, 39, 50); + MITK_TEST_OUTPUT(<< "Testing mitk::Image::GetPixelValueByIndex"); + double val = image->GetPixelValueByIndex(index); + MITK_TEST_CONDITION_REQUIRED( mitk::Equal(val,112.22475433349609), ""); //test by world coordinates + MITK_TEST_OUTPUT(<< "Testing mitk::Image::GetPixelValueByWorldCoordinate"); + mitk::Point3D point; mitk::FillVector3D(point, -5.93752, 18.7199, 6.74218); val = image->GetPixelValueByWorldCoordinate(point); + MITK_TEST_CONDITION_REQUIRED( mitk::Equal(val,94.456184387207031), ""); - if(abs(94.4562 - val) > 1.0){ - std::cout << "[FAILED]"; - return EXIT_FAILURE; - } + MITK_TEST_OUTPUT(<< "Convert to index and access value by mitk::Image::GetPixelValueByIndex again"); + mitk::Index3D index2; + image->GetGeometry()->WorldToIndex(point, index2); + float val2 = image->GetPixelValueByWorldCoordinate(point); + MITK_TEST_CONDITION_REQUIRED( mitk::Equal(val,94.456184387207031), ""); + //access via itk + MITK_TEST_OUTPUT(<< "Test conversion to itk::Image"); + typedef itk::Image ItkFloatImage3D; + ItkFloatImage3D::Pointer itkimage; + mitk::CastToItkImage(image, itkimage); + std::cout<<"[PASSED]"<itk-physical->world consistency"); + mitk::Point3D itkPhysicalPoint; + image->GetGeometry()->WorldToItkPhysicalPoint(point, itkPhysicalPoint); + + mitk::Point3D backTransformedPoint; + image->GetGeometry()->ItkPhysicalPointToWorld(itkPhysicalPoint, backTransformedPoint); + MITK_TEST_CONDITION_REQUIRED( mitk::Equal(point,backTransformedPoint), ""); + + MITK_TEST_OUTPUT(<< "Compare value of pixel returned by mitk in comparison to itk"); + itk::Index<3> idx; + itkimage->TransformPhysicalPointToIndex(itkPhysicalPoint, idx); + float valByItk = itkimage->GetPixel(idx); + + MITK_TEST_CONDITION_REQUIRED( mitk::Equal(valByItk,94.456184387207031), ""); + MITK_TEST_END(); return EXIT_SUCCESS; Index: Documentation/Doxygen/Modules/ModuleGeometry.dox =================================================================== --- Documentation/Doxygen/Modules/ModuleGeometry.dox (revision 24914) +++ Documentation/Doxygen/Modules/ModuleGeometry.dox (working copy) @@ -39,7 +39,7 @@ Image has a TimeSlicedGeometry, which contains one or more SlicedGeometry3D instances (one for each time step), all of which contain one or more instances of (sub-classes of) Geometry2D (usually PlaneGeometry2D). -Converting coordinates from the ITK physical coordinate system (which does not support rotated images) to the MITK world coordinate system should be performed via the Geometry3D of the Image, see Geometry3D::WorldToItkPhysicalPoint. +\deprecated For ITK rev. 3.8 and earlier: Converting coordinates from the ITK physical coordinate system (which did not support rotated images for ITK v3.8 and earlier) to the MITK world coordinate system should be performed via the Geometry3D of the Image, see Geometry3D::WorldToItkPhysicalPoint. As a reminder: Geometry instances referring to images need a slightly different definition of corners, see Geometry3D::SetImageGeometry. This is usualy automatically called by Image. */ Index: Modules/MitkExt/DataManagement/mitkBoundingObject.cpp =================================================================== --- Modules/MitkExt/DataManagement/mitkBoundingObject.cpp (revision 24914) +++ Modules/MitkExt/DataManagement/mitkBoundingObject.cpp (working copy) @@ -43,7 +43,17 @@ { GetGeometry()->SetIdentity(); GetGeometry()->Compose(aGeometry3D->GetIndexToWorldTransform()); - GetGeometry()->SetOrigin(aGeometry3D->GetCenter()); + + if(aGeometry3D->GetImageGeometry()) + { + mitk::Vector3D halfSpacingStep = aGeometry3D->GetSpacing()*0.5; + GetGeometry()->SetOrigin(aGeometry3D->GetCenter()-halfSpacingStep); + } + else + { + GetGeometry()->SetOrigin(aGeometry3D->GetCenter()); + } + mitk::Vector3D size; for(unsigned int i=0; i < 3; ++i) size[i] = aGeometry3D->GetExtentInMM(i)/2.0; Index: Modules/QmitkExt/QmitkVtkLineProfileWidget.cpp =================================================================== --- Modules/QmitkExt/QmitkVtkLineProfileWidget.cpp (revision 24914) +++ Modules/QmitkExt/QmitkVtkLineProfileWidget.cpp (working copy) @@ -144,13 +144,8 @@ { const PathType::OutputType &continuousIndex = m_DerivedPath->Evaluate( t ); - mitk::Point3D indexPoint; - indexPoint[0] = continuousIndex[0]; - indexPoint[1] = continuousIndex[1]; - indexPoint[2] = continuousIndex[2]; - mitk::Point3D worldPoint; - imageGeometry->IndexToWorld( indexPoint, worldPoint ); + imageGeometry->IndexToWorld( continuousIndex, worldPoint ); if ( i == 0 ) { @@ -158,6 +153,9 @@ } distance += currentWorldPoint.EuclideanDistanceTo( worldPoint ); + + mitk::Index3D indexPoint; + imageGeometry->WorldToIndex( worldPoint, worldPoint ); double intensity = m_Image->GetPixelValueByIndex( indexPoint ); MITK_INFO << t << "/" << distance << ": " << indexPoint << " (" << intensity << ")";