diff --git a/Modules/Core/test/mitkGeometry3DTest.cpp b/Modules/Core/test/mitkGeometry3DTest.cpp
index c4395e1a28..a699e40249 100644
--- a/Modules/Core/test/mitkGeometry3DTest.cpp
+++ b/Modules/Core/test/mitkGeometry3DTest.cpp
@@ -1,622 +1,620 @@
/*============================================================================
The Medical Imaging Interaction Toolkit (MITK)
Copyright (c) German Cancer Research Center (DKFZ)
All rights reserved.
Use of this source code is governed by a 3-clause BSD license that can be
found in the LICENSE file.
============================================================================*/
#include "mitkGeometry3D.h"
#include <vnl/vnl_quaternion.h>
#include <vnl/vnl_quaternion.hxx>
#include "mitkInteractionConst.h"
#include "mitkRotationOperation.h"
#include <mitkImageCast.h>
#include <mitkMatrixConvert.h>
#include "mitkTestingMacros.h"
#include <fstream>
#include <mitkNumericTypes.h>
bool testGetAxisVectorVariants(mitk::Geometry3D *geometry)
{
int direction;
for (direction = 0; direction < 3; ++direction)
{
mitk::Vector3D frontToBack(0.);
switch (direction)
{
case 0:
frontToBack = geometry->GetCornerPoint(false, false, false) - geometry->GetCornerPoint(true, false, false);
break; // 7-3
case 1:
frontToBack = geometry->GetCornerPoint(false, false, false) - geometry->GetCornerPoint(false, true, false);
break; // 7-5
case 2:
frontToBack = geometry->GetCornerPoint(false, false, false) - geometry->GetCornerPoint(false, false, true);
break; // 7-2
}
std::cout << "Testing GetAxisVector(int) vs GetAxisVector(bool, bool, bool): ";
if (mitk::Equal(geometry->GetAxisVector(direction), frontToBack) == false)
{
std::cout << "[FAILED]" << std::endl;
return false;
}
std::cout << "[PASSED]" << std::endl;
}
return true;
}
bool testGetAxisVectorExtent(mitk::Geometry3D *geometry)
{
int direction;
for (direction = 0; direction < 3; ++direction)
{
if (mitk::Equal(geometry->GetAxisVector(direction).GetNorm(), geometry->GetExtentInMM(direction)) == false)
{
std::cout << "[FAILED]" << std::endl;
return false;
}
std::cout << "[PASSED]" << std::endl;
}
return true;
}
// 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 -= mitk::Vector(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();
// itk assumes corner based geometry. If our geometry is center based (imageGoe == true), everything needs to be
// shifted
if (geometry3d->GetImageGeometry())
{
centerBounds[0] -= 0.5;
centerBounds[1] -= 0.5;
centerBounds[2] -= 0.5;
}
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;
}
int testIndexAndWorldConsistencyForVectors(mitk::Geometry3D *geometry3d)
{
MITK_TEST_OUTPUT(<< "Testing consistency of index and world coordinate systems for vectors: ");
mitk::Vector3D xAxisMM = geometry3d->GetAxisVector(0);
mitk::Vector3D xAxisContinuousIndex;
mitk::Vector3D xAxisContinuousIndexDeprecated;
mitk::Point3D p, pIndex, origin;
origin = geometry3d->GetOrigin();
p[0] = xAxisMM[0];
p[1] = xAxisMM[1];
p[2] = xAxisMM[2];
geometry3d->WorldToIndex(p, pIndex);
geometry3d->WorldToIndex(xAxisMM, xAxisContinuousIndexDeprecated);
geometry3d->WorldToIndex(xAxisMM, xAxisContinuousIndex);
MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[0] == pIndex[0], "");
MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[1] == pIndex[1], "");
MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[2] == pIndex[2], "");
MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[0] == xAxisContinuousIndexDeprecated[0], "");
MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[1] == xAxisContinuousIndexDeprecated[1], "");
MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[2] == xAxisContinuousIndexDeprecated[2], "");
MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndexDeprecated[0] == pIndex[0], "");
MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndexDeprecated[1] == pIndex[1], "");
MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndexDeprecated[2] == pIndex[2], "");
geometry3d->IndexToWorld(xAxisContinuousIndex, xAxisContinuousIndex);
geometry3d->IndexToWorld(xAxisContinuousIndexDeprecated, xAxisContinuousIndexDeprecated);
geometry3d->IndexToWorld(pIndex, p);
MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex == xAxisMM, "");
MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[0] == p[0], "");
MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[1] == p[1], "");
MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[2] == p[2], "");
MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndexDeprecated == xAxisMM, "");
MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndexDeprecated[0] == p[0], "");
MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndexDeprecated[1] == p[1], "");
MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndexDeprecated[2] == p[2], "");
return EXIT_SUCCESS;
}
int testIndexAndWorldConsistencyForIndex(mitk::Geometry3D *geometry3d)
{
MITK_TEST_OUTPUT(<< "Testing consistency of index and world coordinate systems: ");
// creating testing data
itk::Index<4> itkIndex4, itkIndex4b;
itk::Index<3> itkIndex3, itkIndex3b;
itk::Index<2> itkIndex2, itkIndex2b;
itk::Index<3> mitkIndex, mitkIndexb;
itkIndex4[0] = itkIndex4[1] = itkIndex4[2] = itkIndex4[3] = 4;
itkIndex3[0] = itkIndex3[1] = itkIndex3[2] = 6;
itkIndex2[0] = itkIndex2[1] = 2;
mitkIndex[0] = mitkIndex[1] = mitkIndex[2] = 13;
// check for consistency
mitk::Point3D point;
geometry3d->IndexToWorld(itkIndex2, point);
geometry3d->WorldToIndex(point, itkIndex2b);
MITK_TEST_CONDITION_REQUIRED(((itkIndex2b[0] == itkIndex2[0]) && (itkIndex2b[1] == itkIndex2[1])),
"Testing itk::index<2> for IndexToWorld/WorldToIndex consistency");
geometry3d->IndexToWorld(itkIndex3, point);
geometry3d->WorldToIndex(point, itkIndex3b);
MITK_TEST_CONDITION_REQUIRED(
((itkIndex3b[0] == itkIndex3[0]) && (itkIndex3b[1] == itkIndex3[1]) && (itkIndex3b[2] == itkIndex3[2])),
"Testing itk::index<3> for IndexToWorld/WorldToIndex consistency");
geometry3d->IndexToWorld(itkIndex4, point);
geometry3d->WorldToIndex(point, itkIndex4b);
MITK_TEST_CONDITION_REQUIRED(((itkIndex4b[0] == itkIndex4[0]) && (itkIndex4b[1] == itkIndex4[1]) &&
(itkIndex4b[2] == itkIndex4[2]) && (itkIndex4b[3] == 0)),
"Testing itk::index<3> for IndexToWorld/WorldToIndex consistency");
geometry3d->IndexToWorld(mitkIndex, point);
geometry3d->WorldToIndex(point, mitkIndexb);
MITK_TEST_CONDITION_REQUIRED(
((mitkIndexb[0] == mitkIndex[0]) && (mitkIndexb[1] == mitkIndex[1]) && (mitkIndexb[2] == mitkIndex[2])),
"Testing mitk::Index for IndexToWorld/WorldToIndex consistency");
return EXIT_SUCCESS;
}
#include <itkImage.h>
int testItkImageIsCenterBased()
{
MITK_TEST_OUTPUT(<< "Testing whether itk::Image coordinates are center-based.");
typedef itk::Image<int, 3> 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]" << std::endl;
MITK_TEST_OUTPUT(<< " Testing itk::Image::TransformPhysicalPointToContinuousIndex(origin)==(0,0,0)");
mitk::Point3D globalOrigin;
mitk::FillVector3D(globalOrigin, 0, 0, 0);
- itk::ContinuousIndex<mitk::ScalarType, 3> originContinuousIndex;
- itkintimage->TransformPhysicalPointToContinuousIndex(origin, originContinuousIndex);
+ auto originContinuousIndex = itkintimage->TransformPhysicalPointToContinuousIndex<mitk::ScalarType, mitk::ScalarType>(origin);
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> itkindex = itkintimage->TransformPhysicalPointToIndex(origin);
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);
+ itkindex = itkintimage->TransformPhysicalPointToIndex(originOffCenter);
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 -= mitk::Vector(0.0001);
- itkintimage->TransformPhysicalPointToIndex(originOffCenter, itkindex);
+ itkindex = itkintimage->TransformPhysicalPointToIndex(originOffCenter);
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);
+ itkindex = itkintimage->TransformPhysicalPointToIndex(originOffCenter);
MITK_TEST_CONDITION_REQUIRED(itkindex == global111, "");
MITK_TEST_OUTPUT(<< "=> Yes, itk::Image coordinates are center-based.");
return EXIT_SUCCESS;
}
int testGeometry3D(bool imageGeometry)
{
// Build up a new image Geometry
mitk::Geometry3D::Pointer geometry3d = mitk::Geometry3D::New();
float bounds[] = {-10.0, 17.0, -12.0, 188.0, 13.0, 211.0};
MITK_TEST_OUTPUT(<< "Initializing");
geometry3d->Initialize();
MITK_TEST_OUTPUT(<< "Setting ImageGeometry to " << imageGeometry);
geometry3d->SetImageGeometry(imageGeometry);
MITK_TEST_OUTPUT(<< "Setting bounds by SetFloatBounds(): " << bounds);
geometry3d->SetFloatBounds(bounds);
MITK_TEST_OUTPUT(<< "Testing AxisVectors");
if (testGetAxisVectorVariants(geometry3d) == false)
return EXIT_FAILURE;
if (testGetAxisVectorExtent(geometry3d) == false)
return EXIT_FAILURE;
MITK_TEST_OUTPUT(<< "Creating an AffineTransform3D transform");
mitk::AffineTransform3D::MatrixType matrix;
matrix.SetIdentity();
matrix(1, 1) = 2;
mitk::AffineTransform3D::Pointer transform;
transform = mitk::AffineTransform3D::New();
transform->SetMatrix(matrix);
MITK_TEST_OUTPUT(<< "Testing a SetIndexToWorldTransform");
geometry3d->SetIndexToWorldTransform(transform);
MITK_TEST_OUTPUT(<< "Testing correctness of value returned by GetSpacing");
const mitk::Vector3D &spacing1 = geometry3d->GetSpacing();
mitk::Vector3D expectedSpacing;
expectedSpacing.Fill(1.0);
expectedSpacing[1] = 2;
if (mitk::Equal(spacing1, expectedSpacing) == false)
{
MITK_TEST_OUTPUT(<< " [FAILED]");
return EXIT_FAILURE;
}
MITK_TEST_OUTPUT(<< "Testing a Compose(transform)");
geometry3d->Compose(transform);
MITK_TEST_OUTPUT(<< "Testing correctness of value returned by GetSpacing");
const mitk::Vector3D &spacing2 = geometry3d->GetSpacing();
expectedSpacing[1] = 4;
if (mitk::Equal(spacing2, expectedSpacing) == false)
{
MITK_TEST_OUTPUT(<< " [FAILED]");
return EXIT_FAILURE;
}
MITK_TEST_OUTPUT(<< "Testing correctness of SetSpacing");
mitk::Vector3D newspacing;
mitk::FillVector3D(newspacing, 1.5, 2.5, 3.5);
geometry3d->SetSpacing(newspacing);
const mitk::Vector3D &spacing3 = geometry3d->GetSpacing();
if (mitk::Equal(spacing3, newspacing) == false)
{
MITK_TEST_OUTPUT(<< " [FAILED]");
return EXIT_FAILURE;
}
// Separate Test function for Index and World consistency
testIndexAndWorldConsistency(geometry3d);
testIndexAndWorldConsistencyForVectors(geometry3d);
testIndexAndWorldConsistencyForIndex(geometry3d);
MITK_TEST_OUTPUT(<< "Testing a rotation of the geometry");
double angle = 35.0;
mitk::Vector3D rotationVector;
mitk::FillVector3D(rotationVector, 1, 0, 0);
mitk::Point3D center = geometry3d->GetCenter();
auto op = new mitk::RotationOperation(mitk::OpROTATE, center, rotationVector, angle);
geometry3d->ExecuteOperation(op);
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(voxelStep, voxelStepIndex);
mitk::Vector3D expectedVoxelStepIndex;
expectedVoxelStepIndex.Fill(1);
MITK_TEST_CONDITION_REQUIRED(mitk::Equal(voxelStepIndex, expectedVoxelStepIndex), "");
delete op;
std::cout << "[PASSED]" << std::endl;
MITK_TEST_OUTPUT(<< "Testing that ImageGeometry is still " << imageGeometry);
MITK_TEST_CONDITION_REQUIRED(geometry3d->GetImageGeometry() == imageGeometry, "");
// Test if the translate function moves the origin correctly.
mitk::Point3D oldOrigin = geometry3d->GetOrigin();
// use some random values for translation
mitk::Vector3D translationVector;
translationVector.SetElement(0, 17.5f);
translationVector.SetElement(1, -32.3f);
translationVector.SetElement(2, 4.0f);
// compute ground truth
mitk::Point3D tmpResult = geometry3d->GetOrigin() + translationVector;
geometry3d->Translate(translationVector);
MITK_TEST_CONDITION(mitk::Equal(geometry3d->GetOrigin(), tmpResult), "Testing if origin was translated.");
translationVector *= -1; // vice versa
geometry3d->Translate(translationVector);
MITK_TEST_CONDITION(mitk::Equal(geometry3d->GetOrigin(), oldOrigin),
"Testing if the translation could be done vice versa.");
return EXIT_SUCCESS;
}
int testGeometryAfterCasting()
{
// Epsilon. Allowed difference for rotationvalue
float eps = 0.0001;
// Cast ITK and MITK images and see if geometry stays
typedef itk::Image<char, 2> Image2DType;
typedef itk::Image<char, 3> Image3DType;
// Create 3D ITK Image from Scratch, cast to 3D MITK image, compare Geometries
Image3DType::Pointer image3DItk = Image3DType::New();
Image3DType::RegionType myRegion;
Image3DType::SizeType mySize;
Image3DType::IndexType myIndex;
Image3DType::SpacingType mySpacing;
Image3DType::DirectionType myDirection, rotMatrixX, rotMatrixY, rotMatrixZ;
mySpacing[0] = 31;
mySpacing[1] = 0.1;
mySpacing[2] = 2.9;
myIndex[0] = -15;
myIndex[1] = 15;
myIndex[2] = 12;
mySize[0] = 10;
mySize[1] = 2;
mySize[2] = 5;
myRegion.SetSize(mySize);
myRegion.SetIndex(myIndex);
image3DItk->SetSpacing(mySpacing);
image3DItk->SetRegions(myRegion);
image3DItk->Allocate();
image3DItk->FillBuffer(0);
myDirection.SetIdentity();
rotMatrixX.SetIdentity();
rotMatrixY.SetIdentity();
rotMatrixZ.SetIdentity();
mitk::Image::Pointer mitkImage;
// direction [row] [column]
MITK_TEST_OUTPUT(<< "Casting a rotated 3D ITK Image to a MITK Image and check if Geometry is still same");
for (double rotX = 0; rotX < itk::Math::pi * 2; rotX += itk::Math::pi * 0.4)
{
// Set Rotation X
rotMatrixX[1][1] = cos(rotX);
rotMatrixX[1][2] = -sin(rotX);
rotMatrixX[2][1] = sin(rotX);
rotMatrixX[2][2] = cos(rotX);
for (double rotY = 0; rotY < itk::Math::pi * 2; rotY += itk::Math::pi * 0.3)
{
// Set Rotation Y
rotMatrixY[0][0] = cos(rotY);
rotMatrixY[0][2] = sin(rotY);
rotMatrixY[2][0] = -sin(rotY);
rotMatrixY[2][2] = cos(rotY);
for (double rotZ = 0; rotZ < itk::Math::pi * 2; rotZ += itk::Math::pi * 0.2)
{
// Set Rotation Z
rotMatrixZ[0][0] = cos(rotZ);
rotMatrixZ[0][1] = -sin(rotZ);
rotMatrixZ[1][0] = sin(rotZ);
rotMatrixZ[1][1] = cos(rotZ);
// Multiply matrizes
myDirection = myDirection * rotMatrixX * rotMatrixY * rotMatrixZ;
image3DItk->SetDirection(myDirection);
mitk::CastToMitkImage(image3DItk, mitkImage);
const mitk::AffineTransform3D::MatrixType &matrix =
mitkImage->GetGeometry()->GetIndexToWorldTransform()->GetMatrix();
for (int row = 0; row < 3; row++)
{
for (int col = 0; col < 3; col++)
{
double mitkValue = matrix[row][col] / mitkImage->GetGeometry()->GetSpacing()[col];
double itkValue = myDirection[row][col];
double diff = mitkValue - itkValue;
// if you decrease this value, you can see that there might be QUITE high inaccuracy!!!
if (diff > eps) // need to check, how exact it SHOULD be .. since it is NOT EXACT!
{
std::cout << "Had a difference of : " << diff;
std::cout << "Error: Casting altered Geometry!";
std::cout << "ITK Matrix:\n" << myDirection;
std::cout << "Mitk Matrix (With Spacing):\n" << matrix;
std::cout << "Mitk Spacing: " << mitkImage->GetGeometry()->GetSpacing();
MITK_TEST_CONDITION_REQUIRED(false == true, "");
return false;
}
}
}
}
}
}
// Create 2D ITK Image from Scratch, cast to 2D MITK image, compare Geometries
Image2DType::Pointer image2DItk = Image2DType::New();
Image2DType::RegionType myRegion2D;
Image2DType::SizeType mySize2D;
Image2DType::IndexType myIndex2D;
Image2DType::SpacingType mySpacing2D;
Image2DType::DirectionType myDirection2D, rotMatrix;
mySpacing2D[0] = 31;
mySpacing2D[1] = 0.1;
myIndex2D[0] = -15;
myIndex2D[1] = 15;
mySize2D[0] = 10;
mySize2D[1] = 2;
myRegion2D.SetSize(mySize2D);
myRegion2D.SetIndex(myIndex2D);
image2DItk->SetSpacing(mySpacing2D);
image2DItk->SetRegions(myRegion2D);
image2DItk->Allocate();
image2DItk->FillBuffer(0);
myDirection2D.SetIdentity();
rotMatrix.SetIdentity();
// direction [row] [column]
MITK_TEST_OUTPUT(<< "Casting a rotated 2D ITK Image to a MITK Image and check if Geometry is still same");
for (double rotTheta = 0; rotTheta < itk::Math::pi * 2; rotTheta += itk::Math::pi * 0.2)
{
// Set Rotation
rotMatrix[0][0] = cos(rotTheta);
rotMatrix[0][1] = -sin(rotTheta);
rotMatrix[1][0] = sin(rotTheta);
rotMatrix[1][1] = cos(rotTheta);
// Multiply matrizes
myDirection2D = myDirection2D * rotMatrix;
image2DItk->SetDirection(myDirection2D);
mitk::CastToMitkImage(image2DItk, mitkImage);
const mitk::AffineTransform3D::MatrixType &matrix =
mitkImage->GetGeometry()->GetIndexToWorldTransform()->GetMatrix();
// Compare MITK and ITK matrix
for (int row = 0; row < 3; row++)
{
for (int col = 0; col < 3; col++)
{
double mitkValue = matrix[row][col] / mitkImage->GetGeometry()->GetSpacing()[col];
if ((row == 2) && (col == row))
{
if (mitkValue != 1)
{
MITK_TEST_OUTPUT(<< "After casting a 2D ITK to 3D MITK images, MITK matrix values for 0|2, 1|2, 2|0, 2|1 "
"MUST be 0 and value for 2|2 must be 1");
return false;
}
}
else if ((row == 2) || (col == 2))
{
if (mitkValue != 0)
{
MITK_TEST_OUTPUT(<< "After casting a 2D ITK to 3D MITK images, MITK matrix values for 0|2, 1|2, 2|0, 2|1 "
"MUST be 0 and value for 2|2 must be 1");
return false;
}
}
else
{
double itkValue = myDirection2D[row][col];
double diff = mitkValue - itkValue;
// if you decrease this value, you can see that there might be QUITE high inaccuracy!!!
if (diff > eps) // need to check, how exact it SHOULD be .. since it is NOT EXACT!
{
std::cout << "Had a difference of : " << diff;
std::cout << "Error: Casting altered Geometry!";
std::cout << "ITK Matrix:\n" << myDirection2D;
std::cout << "Mitk Matrix (With Spacing):\n" << matrix;
std::cout << "Mitk Spacing: " << mitkImage->GetGeometry()->GetSpacing();
MITK_TEST_CONDITION_REQUIRED(false == true, "");
return false;
}
}
}
}
}
// THIS WAS TESTED:
// 2D ITK -> 2D MITK,
// 3D ITK -> 3D MITK,
// Still need to test: 2D MITK Image with ADDITIONAL INFORMATION IN MATRIX -> 2D ITK
// 1. Possibility: 3x3 MITK matrix can be converted without loss into 2x2 ITK matrix
// 2. Possibility: 3x3 MITK matrix can only be converted with loss into 2x2 ITK matrix
// .. before implementing this, we wait for further development in geometry classes (e.g. Geoemtry3D::SetRotation(..))
return EXIT_SUCCESS;
}
int mitkGeometry3DTest(int /*argc*/, char * /*argv*/ [])
{
MITK_TEST_BEGIN(mitkGeometry3DTest);
int result;
MITK_TEST_CONDITION_REQUIRED((result = testItkImageIsCenterBased()) == EXIT_SUCCESS, "");
MITK_TEST_OUTPUT(<< "Running main part of test with ImageGeometry = false");
MITK_TEST_CONDITION_REQUIRED((result = testGeometry3D(false)) == EXIT_SUCCESS, "");
MITK_TEST_OUTPUT(<< "Running main part of test with ImageGeometry = true");
MITK_TEST_CONDITION_REQUIRED((result = testGeometry3D(true)) == EXIT_SUCCESS, "");
MITK_TEST_OUTPUT(<< "Running test to see if Casting MITK to ITK and the other way around destroys geometry");
MITK_TEST_CONDITION_REQUIRED((result = testGeometryAfterCasting()) == EXIT_SUCCESS, "");
MITK_TEST_END();
return EXIT_SUCCESS;
}
diff --git a/Modules/ImageStatistics/mitkMaskUtilities.tpp b/Modules/ImageStatistics/mitkMaskUtilities.tpp
index ab740cb9f7..645499f9f7 100644
--- a/Modules/ImageStatistics/mitkMaskUtilities.tpp
+++ b/Modules/ImageStatistics/mitkMaskUtilities.tpp
@@ -1,194 +1,191 @@
/*============================================================================
The Medical Imaging Interaction Toolkit (MITK)
Copyright (c) German Cancer Research Center (DKFZ)
All rights reserved.
Use of this source code is governed by a 3-clause BSD license that can be
found in the LICENSE file.
============================================================================*/
#ifndef MITKMASKUTIL_TPP
#define MITKMASKUTIL_TPP
#include <mitkMaskUtilities.h>
#include <mitkImageAccessByItk.h>
#include <itkExtractImageFilter.h>
#include <itkChangeInformationImageFilter.h>
#include <mitkITKImageImport.h>
namespace mitk
{
template <class TPixel, unsigned int VImageDimension>
void MaskUtilities<TPixel, VImageDimension>::SetImage(const ImageType* image)
{
if (image != m_Image)
{
m_Image = image;
}
}
template <class TPixel, unsigned int VImageDimension>
void MaskUtilities<TPixel, VImageDimension>::SetMask(const MaskType* mask)
{
if (mask != m_Mask)
{
m_Mask = mask;
}
}
template <class TPixel, unsigned int VImageDimension>
bool MaskUtilities<TPixel, VImageDimension>::CheckMaskSanity()
{
if (m_Mask==nullptr || m_Image==nullptr)
{
MITK_ERROR << "Set an image and a mask first";
}
typedef itk::Image< TPixel, VImageDimension > ImageType;
typedef typename ImageType::PointType PointType;
typedef typename ImageType::DirectionType DirectionType;
bool maskSanity = true;
if (m_Mask==nullptr)
{
MITK_ERROR << "Something went wrong when casting the mitk mask image to an itk mask image. Do the mask and the input image have the same dimension?";
// note to self: We could try to convert say a 2d mask to a 3d mask if the image is 3d. (mask and image dimension have to match.)
}
// check direction
DirectionType imageDirection = m_Image->GetDirection();
DirectionType maskDirection = m_Mask->GetDirection();
for(unsigned int i = 0; i < imageDirection.ColumnDimensions; ++i )
{
for(unsigned int j = 0; j < imageDirection.ColumnDimensions; ++j )
{
double differenceDirection = imageDirection[i][j] - maskDirection[i][j];
if (fabs(differenceDirection) > MASK_SUITABILITY_TOLERANCE_DIRECTION)
{
maskSanity = false;
MITK_INFO << "Mask needs to have same direction as image! (Image direction: " << imageDirection << "; Mask direction: " << maskDirection << ")";
}
}
}
// check spacing
PointType imageSpacing(m_Image->GetSpacing().GetDataPointer());
PointType maskSpacing(m_Mask->GetSpacing().GetDataPointer());
for (unsigned int i = 0; i < VImageDimension; i++)
{
if ( fabs( maskSpacing[i] - imageSpacing[i] ) > MASK_SUITABILITY_TOLERANCE_COORDINATE )
{
maskSanity = false;
MITK_INFO << "Spacing of mask and image is not equal. Mask: " << maskSpacing << " image: " << imageSpacing;
}
}
// check alignment
// Make sure that the voxels of mask and image are correctly "aligned", i.e., voxel boundaries are the same in both images
PointType imageOrigin = m_Image->GetOrigin();
PointType maskOrigin = m_Mask->GetOrigin();
- typedef itk::ContinuousIndex<double, VImageDimension> ContinousIndexType;
- ContinousIndexType maskOriginContinousIndex, imageOriginContinousIndex;
-
- m_Image->TransformPhysicalPointToContinuousIndex(maskOrigin, maskOriginContinousIndex);
- m_Image->TransformPhysicalPointToContinuousIndex(imageOrigin, imageOriginContinousIndex);
+ auto maskOriginContinousIndex = m_Image->TransformPhysicalPointToContinuousIndex<PointType::ValueType, double>(maskOrigin);
+ auto imageOriginContinousIndex = m_Image->TransformPhysicalPointToContinuousIndex<PointType::ValueType, double>(imageOrigin);
for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i )
{
double misalignment = maskOriginContinousIndex[i] - floor( maskOriginContinousIndex[i] + 0.5 );
// misalignment must be a multiple (int) of spacing in that direction
if ( fmod(misalignment,imageSpacing[i]) > MASK_SUITABILITY_TOLERANCE_COORDINATE)
{
maskSanity = false;
MITK_INFO << "Pixels/voxels of mask and image are not sufficiently aligned! (Misalignment: " << fmod(misalignment,imageSpacing[i]) << ")";
}
}
// mask must be completely inside image region
// Make sure that mask region is contained within image region
if ( m_Mask!=nullptr &&
!m_Image->GetLargestPossibleRegion().IsInside( m_Mask->GetLargestPossibleRegion() ) )
{
maskSanity = false;
MITK_INFO << "Mask region needs to be inside of image region! (Image region: "
<< m_Image->GetLargestPossibleRegion() << "; Mask region: " << m_Mask->GetLargestPossibleRegion() << ")";
}
return maskSanity;
}
template <class TPixel, unsigned int VImageDimension>
typename MaskUtilities<TPixel, VImageDimension >::ImageType::ConstPointer MaskUtilities<TPixel, VImageDimension>::ExtractMaskImageRegion()
{
if (m_Mask==nullptr || m_Image==nullptr)
{
MITK_ERROR << "Set an image and a mask first";
}
bool maskSanity = CheckMaskSanity();
if (!maskSanity)
{
MITK_ERROR << "Mask and image are not compatible";
}
typedef itk::ExtractImageFilter< ImageType, ImageType > ExtractImageFilterType;
typename ImageType::SizeType imageSize = m_Image->GetBufferedRegion().GetSize();
typename ImageType::SizeType maskSize = m_Mask->GetBufferedRegion().GetSize();
typename itk::Image<TPixel, VImageDimension>::ConstPointer resultImg;
bool maskSmallerImage = false;
for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i )
{
if ( maskSize[i] < imageSize[i] )
{
maskSmallerImage = true;
}
}
if ( maskSmallerImage )
{
typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New();
typename MaskType::PointType maskOrigin = m_Mask->GetOrigin();
typename ImageType::PointType imageOrigin = m_Image->GetOrigin();
typename MaskType::SpacingType maskSpacing = m_Mask->GetSpacing();
typename ImageType::RegionType extractionRegion;
typename ImageType::IndexType extractionRegionIndex;
for (unsigned int i=0; i < maskOrigin.GetPointDimension(); i++)
{
extractionRegionIndex[i] = (maskOrigin[i] - imageOrigin[i]) / maskSpacing[i];
}
extractionRegion.SetIndex(extractionRegionIndex);
extractionRegion.SetSize(m_Mask->GetLargestPossibleRegion().GetSize());
extractImageFilter->SetInput( m_Image );
extractImageFilter->SetExtractionRegion( extractionRegion );
extractImageFilter->SetCoordinateTolerance(MASK_SUITABILITY_TOLERANCE_COORDINATE);
extractImageFilter->SetDirectionTolerance(MASK_SUITABILITY_TOLERANCE_DIRECTION);
extractImageFilter->Update();
auto extractedImg = extractImageFilter->GetOutput();
extractedImg->SetOrigin(m_Mask->GetOrigin());
extractedImg->SetLargestPossibleRegion(m_Mask->GetLargestPossibleRegion());
extractedImg->SetBufferedRegion(m_Mask->GetBufferedRegion());
resultImg = extractedImg;
}
else
{
resultImg = m_Image;
}
return resultImg;
}
}
#endif
diff --git a/Modules/SurfaceInterpolation/mitkCreateDistanceImageFromSurfaceFilter.cpp b/Modules/SurfaceInterpolation/mitkCreateDistanceImageFromSurfaceFilter.cpp
index 16952fd4a7..f3d4924082 100644
--- a/Modules/SurfaceInterpolation/mitkCreateDistanceImageFromSurfaceFilter.cpp
+++ b/Modules/SurfaceInterpolation/mitkCreateDistanceImageFromSurfaceFilter.cpp
@@ -1,617 +1,614 @@
/*============================================================================
The Medical Imaging Interaction Toolkit (MITK)
Copyright (c) German Cancer Research Center (DKFZ)
All rights reserved.
Use of this source code is governed by a 3-clause BSD license that can be
found in the LICENSE file.
============================================================================*/
#include "mitkCreateDistanceImageFromSurfaceFilter.h"
#include "mitkImageCast.h"
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkDoubleArray.h"
#include "vtkPolyData.h"
#include "vtkSmartPointer.h"
#include "itkImageRegionIteratorWithIndex.h"
#include "itkNeighborhoodIterator.h"
#include <queue>
void mitk::CreateDistanceImageFromSurfaceFilter::CreateEmptyDistanceImage()
{
// Determine the bounds of the input points in index- and world-coordinates
DistanceImageType::PointType minPointInWorldCoordinates, maxPointInWorldCoordinates;
DistanceImageType::IndexType minPointInIndexCoordinates, maxPointInIndexCoordinates;
DetermineBounds(
minPointInWorldCoordinates, maxPointInWorldCoordinates, minPointInIndexCoordinates, maxPointInIndexCoordinates);
// 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;
for (unsigned int dim = 0; dim < 3; ++dim)
{
extentMM[dim] = (std::abs(maxPointInIndexCoordinates[dim] - minPointInIndexCoordinates[dim])) *
m_ReferenceImage->GetSpacing()[dim];
}
/*
* Now create an empty distance image. The created image will always have the same number of pixels, 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 = ( extentX*extentY*extentZ / 500000 )^(1/3)
*/
double basis = (extentMM[0] * extentMM[1] * extentMM[2]) / m_DistanceImageVolume;
double exponent = 1.0 / 3.0;
m_DistanceImageSpacing = pow(basis, exponent);
// 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;
// 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 + 8;
sizeOfRegion[1] = numberOfYPixel + 8;
sizeOfRegion[2] = numberOfZPixel + 8;
// The region starts at index 0,0,0
DistanceImageType::IndexType initialOriginAsIndex;
initialOriginAsIndex.Fill(0);
DistanceImageType::PointType originAsWorld = minPointInWorldCoordinates;
DistanceImageType::RegionType lpRegion;
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
m_DistanceImageITK = DistanceImageType::New();
m_DistanceImageITK->SetOrigin(originAsWorld);
m_DistanceImageITK->SetDirection(m_ReferenceImage->GetDirection());
m_DistanceImageITK->SetRegions(lpRegion);
m_DistanceImageITK->SetSpacing(itk::Vector<double, 3>(m_DistanceImageSpacing));
m_DistanceImageITK->Allocate();
// First of all the image is initialized with the value 10*m_DistanceImageSpacing for each pixel
m_DistanceImageDefaultBufferValue = 10 * m_DistanceImageSpacing;
m_DistanceImageITK->FillBuffer(m_DistanceImageDefaultBufferValue);
// Now we move the origin of the distanceImage 2 index-Coordinates
// in all directions
- DistanceImageType::IndexType originAsIndex;
- m_DistanceImageITK->TransformPhysicalPointToIndex(originAsWorld, originAsIndex);
+ auto originAsIndex = m_DistanceImageITK->TransformPhysicalPointToIndex(originAsWorld);
originAsIndex[0] -= 2;
originAsIndex[1] -= 2;
originAsIndex[2] -= 2;
m_DistanceImageITK->TransformIndexToPhysicalPoint(originAsIndex, originAsWorld);
m_DistanceImageITK->SetOrigin(originAsWorld);
}
mitk::CreateDistanceImageFromSurfaceFilter::CreateDistanceImageFromSurfaceFilter()
: m_DistanceImageSpacing(0.0), m_DistanceImageDefaultBufferValue(0.0)
{
m_DistanceImageVolume = 50000;
this->m_UseProgressBar = false;
this->m_ProgressStepSize = 5;
mitk::Image::Pointer output = mitk::Image::New();
this->SetNthOutput(0, output.GetPointer());
}
mitk::CreateDistanceImageFromSurfaceFilter::~CreateDistanceImageFromSurfaceFilter()
{
}
void mitk::CreateDistanceImageFromSurfaceFilter::GenerateData()
{
this->PreprocessContourPoints();
this->CreateEmptyDistanceImage();
// First of all we have to build the equation-system from the existing contour-edge-points
this->CreateSolutionMatrixAndFunctionValues();
if (this->m_UseProgressBar)
mitk::ProgressBar::GetInstance()->Progress(1);
m_Weights = m_SolutionMatrix.partialPivLu().solve(m_FunctionValues);
if (this->m_UseProgressBar)
mitk::ProgressBar::GetInstance()->Progress(2);
// The last step is to create the distance map with the interpolated distance function
this->FillDistanceImage();
if (this->m_UseProgressBar)
mitk::ProgressBar::GetInstance()->Progress(2);
m_Centers.clear();
m_Normals.clear();
}
void mitk::CreateDistanceImageFromSurfaceFilter::PreprocessContourPoints()
{
unsigned int numberOfInputs = this->GetNumberOfIndexedInputs();
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
vtkSmartPointer<vtkPolyData> polyData;
vtkSmartPointer<vtkDoubleArray> currentCellNormals;
vtkSmartPointer<vtkCellArray> existingPolys;
vtkSmartPointer<vtkPoints> existingPoints;
double p[3];
PointType currentPoint;
PointType normal;
for (unsigned int i = 0; i < numberOfInputs; i++)
{
auto currentSurface = 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();
const vtkIdType *cell(nullptr);
vtkIdType cellSize(0);
for (existingPolys->InitTraversal(); existingPolys->GetNextCell(cellSize, cell);)
{
for (vtkIdType 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
}
void mitk::CreateDistanceImageFromSurfaceFilter::CreateSolutionMatrixAndFunctionValues()
{
// 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.resize(numberOfCenters * 3);
m_FunctionValues.fill(0);
PointType currentPoint;
PointType normal;
// 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] * m_DistanceImageSpacing;
currentPoint[1] = currentPoint[1] - normal[1] * m_DistanceImageSpacing;
currentPoint[2] = currentPoint[2] - normal[2] * m_DistanceImageSpacing;
m_Centers.push_back(currentPoint);
m_FunctionValues[numberOfCenters + i] = -m_DistanceImageSpacing;
}
// 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] * m_DistanceImageSpacing;
currentPoint[1] = currentPoint[1] + normal[1] * m_DistanceImageSpacing;
currentPoint[2] = currentPoint[2] + normal[2] * m_DistanceImageSpacing;
m_Centers.push_back(currentPoint);
m_FunctionValues[numberOfCenters * 2 + i] = m_DistanceImageSpacing;
}
// Now we have created all centers and all function values. Next step is to create the solution matrix
numberOfCenters = m_Centers.size();
m_SolutionMatrix.resize(numberOfCenters, numberOfCenters);
m_Weights.resize(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::FillDistanceImage()
{
/*
* 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.
*/
typedef itk::ImageRegionIteratorWithIndex<DistanceImageType> ImageIterator;
typedef itk::NeighborhoodIterator<DistanceImageType> NeighborhoodImageIterator;
std::queue<DistanceImageType::IndexType> 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;
- m_DistanceImageITK->TransformPhysicalPointToIndex(currentPointAsPoint, currentIndex);
+ auto currentIndex = m_DistanceImageITK->TransformPhysicalPointToIndex(currentPointAsPoint);
assert(
m_DistanceImageITK->GetLargestPossibleRegion().IsInside(currentIndex)); // we are quite certain this should hold
narrowbandPoints.push(currentIndex);
m_DistanceImageITK->SetPixel(currentIndex, distance);
NeighborhoodImageIterator::RadiusType radius;
radius.Fill(1);
NeighborhoodImageIterator nIt(radius, m_DistanceImageITK, m_DistanceImageITK->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(*relativeNb, isInBounds);
if (isInBounds && nIt.GetPixel(*relativeNb) == m_DistanceImageDefaultBufferValue)
{
currentIndex = nIt.GetIndex(*relativeNb);
// Transform the currently checked point from index-coordinates to
// world-coordinates
m_DistanceImageITK->TransformIndexToPhysicalPoint(currentIndex, currentPointAsPoint);
// 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 (std::fabs(distance) <= m_DistanceImageSpacing * 2)
{
nIt.SetPixel(*relativeNb, distance);
narrowbandPoints.push(currentIndex);
}
}
relativeNb++;
}
}
ImageIterator imgRegionIterator(m_DistanceImageITK, m_DistanceImageITK->GetLargestPossibleRegion());
imgRegionIterator.GoToBegin();
double prevPixelVal = 1;
DistanceImageType::IndexType _size;
_size.Fill(-1);
_size += m_DistanceImageITK->GetLargestPossibleRegion().GetSize();
// Set every pixel inside the surface to -m_DistanceImageDefaultBufferValue except the edge point (so that the
// received surface is closed)
while (!imgRegionIterator.IsAtEnd())
{
if (imgRegionIterator.Get() == m_DistanceImageDefaultBufferValue && prevPixelVal < 0)
{
while (imgRegionIterator.Get() == m_DistanceImageDefaultBufferValue)
{
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(m_DistanceImageDefaultBufferValue);
prevPixelVal = m_DistanceImageDefaultBufferValue;
++imgRegionIterator;
break;
}
else
{
imgRegionIterator.Set((-1) * m_DistanceImageDefaultBufferValue);
++imgRegionIterator;
prevPixelVal = (-1) * m_DistanceImageDefaultBufferValue;
}
}
}
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(m_DistanceImageDefaultBufferValue);
prevPixelVal = m_DistanceImageDefaultBufferValue;
++imgRegionIterator;
}
else
{
prevPixelVal = imgRegionIterator.Get();
++imgRegionIterator;
}
}
Image::Pointer resultImage = this->GetOutput();
// Cast the created distance-Image from itk::Image to the mitk::Image
// that is our output.
CastToMitkImage(m_DistanceImageITK, resultImage);
}
double mitk::CreateDistanceImageFromSurfaceFilter::CalculateDistanceValue(PointType p)
{
double distanceValue(0);
PointType p1;
PointType p2;
double norm;
CenterList::iterator centerIter;
unsigned int count(0);
for (centerIter = m_Centers.begin(); centerIter != m_Centers.end(); centerIter++)
{
p1 = *centerIter;
p2 = p - p1;
norm = p2.two_norm();
distanceValue = distanceValue + (norm * m_Weights[count]);
++count;
}
return distanceValue;
}
void mitk::CreateDistanceImageFromSurfaceFilter::GenerateOutputInformation()
{
}
void mitk::CreateDistanceImageFromSurfaceFilter::PrintEquationSystem()
{
std::stringstream out;
out << "Nummber of rows: " << m_SolutionMatrix.rows() << " ****** Number of columns: " << m_SolutionMatrix.cols()
<< endl;
out << "[ ";
for (int i = 0; i < m_SolutionMatrix.rows(); i++)
{
for (int j = 0; j < m_SolutionMatrix.cols(); j++)
{
out << m_SolutionMatrix(i, j) << " ";
}
out << ";" << endl;
}
out << " ]\n\n\n";
for (unsigned int i = 0; i < m_Centers.size(); i++)
{
out << m_Centers.at(i) << ";" << endl;
}
std::cout << "Equation system: \n\n\n" << out.str();
}
void mitk::CreateDistanceImageFromSurfaceFilter::SetInput(const mitk::Surface *surface)
{
this->SetInput(0, surface);
}
void mitk::CreateDistanceImageFromSurfaceFilter::SetInput(unsigned int idx, const mitk::Surface *surface)
{
if (this->GetInput(idx) != surface)
{
this->SetNthInput(idx, const_cast<mitk::Surface *>(surface));
this->Modified();
}
}
const mitk::Surface *mitk::CreateDistanceImageFromSurfaceFilter::GetInput()
{
if (this->GetNumberOfIndexedInputs() < 1)
return nullptr;
return static_cast<const mitk::Surface *>(this->ProcessObject::GetInput(0));
}
const mitk::Surface *mitk::CreateDistanceImageFromSurfaceFilter::GetInput(unsigned int idx)
{
if (this->GetNumberOfIndexedInputs() < 1)
return nullptr;
return static_cast<const mitk::Surface *>(this->ProcessObject::GetInput(idx));
}
void mitk::CreateDistanceImageFromSurfaceFilter::RemoveInputs(mitk::Surface *input)
{
DataObjectPointerArraySizeType nb = this->GetNumberOfIndexedInputs();
for (DataObjectPointerArraySizeType i = 0; i < nb; i++)
{
if (this->GetInput(i) == input)
{
this->RemoveInput(i);
return;
}
}
}
void mitk::CreateDistanceImageFromSurfaceFilter::Reset()
{
for (unsigned int i = 0; i < this->GetNumberOfIndexedInputs(); i++)
{
this->PopBackInput();
}
this->SetNumberOfIndexedInputs(0);
this->SetNumberOfIndexedOutputs(1);
mitk::Image::Pointer output = mitk::Image::New();
this->SetNthOutput(0, output.GetPointer());
}
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::ImageBase<3>::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
- itk::ContinuousIndex<double, 3> tmpIndex;
- m_ReferenceImage->TransformPhysicalPointToContinuousIndex(tmpPoint, tmpIndex);
+ auto tmpIndex = m_ReferenceImage->TransformPhysicalPointToContinuousIndex<DistanceImageType::PointType::ValueType>(tmpPoint);
// initialize the variables with this first point
DistanceImageType::IndexValueType xmin = tmpIndex[0];
DistanceImageType::IndexValueType ymin = tmpIndex[1];
DistanceImageType::IndexValueType zmin = tmpIndex[2];
DistanceImageType::IndexValueType xmax = tmpIndex[0];
DistanceImageType::IndexValueType ymax = tmpIndex[1];
DistanceImageType::IndexValueType zmax = tmpIndex[2];
// iterate over the rest of the points
auto 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->TransformPhysicalPointToContinuousIndex(tmpPoint, tmpIndex);
+ tmpIndex = m_ReferenceImage->TransformPhysicalPointToContinuousIndex<DistanceImageType::PointType::ValueType>(tmpPoint);
// 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);
}