From 8df460e4e8ce924e01dd6bd2a83782c9cf1fdd9d Mon Sep 17 00:00:00 2001 From: Mikhail Polkovnikov Date: Sat, 2 Jul 2022 07:18:30 +0300 Subject: [PATCH] ENH: Calculate parallel geometry parameters from projection matrix --- include/rtkThreeDCircularProjectionGeometry.h | 21 +++- src/rtkThreeDCircularProjectionGeometry.cxx | 56 +++++++++- test/CMakeLists.txt | 1 + test/rtkparallelgeometryfrommatrixtest.cxx | 101 ++++++++++++++++++ 4 files changed, 173 insertions(+), 6 deletions(-) create mode 100644 test/rtkparallelgeometryfrommatrixtest.cxx diff --git a/include/rtkThreeDCircularProjectionGeometry.h b/include/rtkThreeDCircularProjectionGeometry.h index 09ca5c2f4..642676118 100644 --- a/include/rtkThreeDCircularProjectionGeometry.h +++ b/include/rtkThreeDCircularProjectionGeometry.h @@ -115,7 +115,25 @@ class RTK_EXPORT ThreeDCircularProjectionGeometry : public ProjectionGeometry<3> /** Add projection from a projection matrix. A projection matrix is defined * up to a scaling factor. The function here Assumes that the input matrix * pMat is normalized such that pMat*(x,y,z,1)'=(u,v,1)'. - * This code assumes that the SourceToDetectorDistance is positive. */ + * This code assumes that the SourceToDetectorDistance is positive. + * + * \warning For parallel geometry the matrix has a special form + * + * | R11 R12 R13 -ProjectionOffsetX | + * Mp =| R21 R22 R23 -ProjectionOffsetY | + * | 0 0 0 1 | + * + * where + * + * | R11 R12 R13 | + * | R21 R22 R23 | - cropped rotation 3x3 matrix + * | 0 0 0 | + * + * outOfPlaneAngle must be within (-90, 90) degrees range, + * otherwise the rotation angles could be miscalculated. + * + * If projection matrix is for parallel geometry, then SID is set to 1000. mm + */ bool AddProjection(const HomogeneousProjectionMatrixType & pMat); @@ -266,7 +284,6 @@ class RTK_EXPORT ThreeDCircularProjectionGeometry : public ProjectionGeometry<3> return this->m_MagnificationMatrices[i]; } - /** Get the vector containing the collimation jaw parameters. */ const std::vector & GetCollimationUInf() const diff --git a/src/rtkThreeDCircularProjectionGeometry.cxx b/src/rtkThreeDCircularProjectionGeometry.cxx index 0d6371d15..4f408f179 100644 --- a/src/rtkThreeDCircularProjectionGeometry.cxx +++ b/src/rtkThreeDCircularProjectionGeometry.cxx @@ -165,7 +165,7 @@ rtk::ThreeDCircularProjectionGeometry::AddProjection(const PointType & sourcePo const PointType & S = sourcePosition; // source pos const PointType & R = detectorPosition; // detector pos - if (itk::Math::abs(r * c) > 1e-6) // non-orthogonal row/column vectors + if (fabs(r * c) > 1e-6) // non-orthogonal row/column vectors return false; // Euler angles (ZXY convention) from detector orientation in IEC-based WCS: @@ -222,7 +222,11 @@ rtk::ThreeDCircularProjectionGeometry::AddProjection(const PointType & sourcePo double SID = n[0] * S[0] + n[1] * S[1] + n[2] * S[2]; // SDD: distance from source to detector along detector normal double SDD = n[0] * (S[0] - R[0]) + n[1] * (S[1] - R[1]) + n[2] * (S[2] - R[2]); - + if (itk::Math::abs(SDD) < 1.0E-06) + { + itkDebugMacro(<< "SourceDetectorDistance is less than 1.0E-06. For parallel geometry SDD is set to 0.0"); + SDD = 0.0; + } // source offset: compute source's "in-plane" x/y shift off isocenter VectorType Sv; Sv[0] = S[0]; @@ -264,6 +268,50 @@ rtk::ThreeDCircularProjectionGeometry::AddProjection(const HomogeneousProjection double d = pMat(0, 0) * pMat(1, 1) * pMat(2, 2) + pMat(0, 1) * pMat(1, 2) * pMat(2, 0) + pMat(0, 2) * pMat(1, 0) * pMat(2, 1) - pMat(0, 0) * pMat(1, 2) * pMat(2, 1) - pMat(0, 1) * pMat(1, 0) * pMat(2, 2) - pMat(0, 2) * pMat(1, 1) * pMat(2, 0); + + if (itk::Math::abs(d) < std::numeric_limits::epsilon()) // parallel geometry + { + // setup rotation matrix from three rotated unit vectors + // v1, v2 - standard unit vectors, n - normal (cross product) [v1,v2] + VectorType v1, v2; + v1[0] = A(0, 0); + v1[1] = A(0, 1); + v1[2] = A(0, 2); + v2[0] = A(1, 0); + v2[1] = A(1, 1); + v2[2] = A(1, 2); + VectorType n = itk::CrossProduct(v1, v2); + A(2, 0) = n[0]; + A(2, 1) = n[1]; + A(2, 2) = n[2]; + + // Declare a 3D euler transform in order to properly extract angles + using EulerType = itk::Euler3DTransform; + EulerType::Pointer euler = EulerType::New(); + euler->SetComputeZYX(false); // ZXY order + + // Extract angle using parent method without orthogonality check + euler->itk::MatrixOffsetTransformBase::SetMatrix(A); + double oa = euler->GetAngleX(); + double ga = euler->GetAngleY(); + double ia = euler->GetAngleZ(); + + // verify that extracted ZXY angles result in the *desired* matrix: + // (at some angle constellations we may run into numerical troubles, therefore, + // verify angles and try to fix instabilities) + if (!VerifyAngles(oa, ga, ia, A)) + { + if (!FixAngles(oa, ga, ia, A)) + { + itkWarningMacro(<< "Failed to AddProjection"); + return false; + } + } + + this->AddProjectionInRadians(1000., 0., -1. * ga, -1. * p[0], -1. * p[1], -1. * oa, -1. * ia); + return true; + } + d = -1. * d / itk::Math::abs(d); // Extract intrinsic parameters u0, v0 and f (f is chosen to be positive at that point) @@ -637,7 +685,7 @@ rtk::ThreeDCircularProjectionGeometry::VerifyAngles(const double outOfP for (int i = 0; i < 3; i++) // check whether matrices match for (int j = 0; j < 3; j++) - if (itk::Math::abs(rm[i][j] - m[i][j]) > EPSILON) + if (fabs(rm[i][j] - m[i][j]) > EPSILON) return false; return true; @@ -652,7 +700,7 @@ rtk::ThreeDCircularProjectionGeometry::FixAngles(double & outOfPlan const Matrix3x3Type & rm = referenceMatrix; // shortcut const double EPSILON = 1e-6; // internal tolerance for comparison - if (itk::Math::abs(itk::Math::abs(rm[2][1]) - 1.) > EPSILON) + if (fabs(fabs(rm[2][1]) - 1.) > EPSILON) { double oa = NAN, ga = NAN, ia = NAN; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 19231e61a..3f68dff8e 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -285,6 +285,7 @@ rtk_add_test(rtkArgsInfoManagerTest rtkargsinfomanagertest.cxx) rtk_add_test(rtkGeometryCloneTest rtkgeometryclonetest.cxx) rtk_add_test(rtkGeometryFromMatrixTest rtkgeometryfrommatrixtest.cxx) +rtk_add_test(rtkParallelGeometryFromMatrixTest rtkparallelgeometryfrommatrixtest.cxx) rtk_add_test(rtkOraTest rtkoratest.cxx DATA{Input/Ora/0_afterLog.ora.xml,0_afterLog.mhd,0_afterLog.raw} diff --git a/test/rtkparallelgeometryfrommatrixtest.cxx b/test/rtkparallelgeometryfrommatrixtest.cxx new file mode 100644 index 000000000..4df118c75 --- /dev/null +++ b/test/rtkparallelgeometryfrommatrixtest.cxx @@ -0,0 +1,101 @@ +// RTK +#include "rtkThreeDCircularProjectionGeometry.h" + + +int +main(int, char **) +{ + rtk::ThreeDCircularProjectionGeometry::Pointer geometry1, geometry2; + geometry1 = rtk::ThreeDCircularProjectionGeometry::New(); + geometry2 = rtk::ThreeDCircularProjectionGeometry::New(); + + const double sid = 1000.; + const double sdd = 0.; // Parallel geometry + for (double oa = -60.; oa < 65.; oa += 60.) + { + for (double ia = -360.; ia < 360.; ia += 60.) + { + for (double ga = -360.; ga < 360.; ga += 45.) + { + for (double px = -50.; px < 55.; px += 25.) + { + for (double py = -50.; py < 55.; py += 25.) + { + geometry1->AddProjection(sid, sdd, ga, px, py, oa, ia); + geometry2->AddProjection(geometry1->GetMatrices().back()); + + double g1 = geometry1->GetGantryAngles().back(); + double g2 = geometry2->GetGantryAngles().back(); + double d = g1 - g2; + d = itk::Math::abs(geometry1->ConvertAngleBetweenMinusAndPlusPIRadians(d)); + if (d > 0.01) + { + std::cerr << "ERROR: GantryAngles 1 is " << g1 << " and GantryAngles 2 is " << g2; + return EXIT_FAILURE; + } + + double ia1 = geometry1->GetInPlaneAngles().back(); + double ia2 = geometry2->GetInPlaneAngles().back(); + d = ia1 - ia2; + d = itk::Math::abs(geometry1->ConvertAngleBetweenMinusAndPlusPIRadians(d)); + if (d > 0.01) + { + std::cerr << "ERROR: InPlaneAngles 1 is " << ia1 << " and InPlaneAngles 2 is " << ia2; + return EXIT_FAILURE; + } + + double oa1 = geometry1->GetOutOfPlaneAngles().back(); + double oa2 = geometry2->GetOutOfPlaneAngles().back(); + d = oa1 - oa2; + d = itk::Math::abs(geometry1->ConvertAngleBetweenMinusAndPlusPIRadians(d)); + if (d > 0.01) + { + std::cerr << "ERROR: OutOfPlaneAngle 1 is " << oa1 << " and OutOfPlaneAngle 2 is " << oa2; + return EXIT_FAILURE; + } + + double py1 = geometry1->GetProjectionOffsetsY().back(); + double py2 = geometry2->GetProjectionOffsetsY().back(); + d = itk::Math::abs(py1 - py2); + if (d > 0.01) + { + std::cerr << "ERROR: ProjectionOffsetsY 1 is " << py1 << " and ProjectionOffsetsY 2 is " << py2; + return EXIT_FAILURE; + } + + double px1 = geometry1->GetProjectionOffsetsX().back(); + double px2 = geometry2->GetProjectionOffsetsX().back(); + d = itk::Math::abs(px1 - px2); + if (d > 0.01) + { + std::cerr << "ERROR: ProjectionOffsetsX 1 is " << px1 << " and ProjectionOffsetsX 2 is " << px2; + return EXIT_FAILURE; + } + + double sid1 = geometry1->GetSourceToIsocenterDistances().back(); + double sid2 = geometry2->GetSourceToIsocenterDistances().back(); + d = itk::Math::abs(sid1 - sid2); + if (d > 0.01) + { + std::cerr << "ERROR: GetSourceToIsocenterDistances 1 is " << sid1 + << " and GetSourceToIsocenterDistances 2 is " << sid2; + return EXIT_FAILURE; + } + + double sdd1 = geometry1->GetSourceToDetectorDistances().back(); + double sdd2 = geometry2->GetSourceToDetectorDistances().back(); + d = itk::Math::abs(sdd1 - sdd2); + if (d > 0.01) + { + std::cerr << "ERROR: GetSourceToDetectorDistances 1 is " << sdd1 + << " and GetSourceToDetectorDistances 2 is " << sdd2; + return EXIT_FAILURE; + } + } + } + } + } + } + std::cout << "Test PASSED! " << std::endl; + return EXIT_SUCCESS; +}