Skip to content

Commit

Permalink
ENH: Calculate parallel geometry parameters from projection matrix
Browse files Browse the repository at this point in the history
  • Loading branch information
Mikhail Polkovnikov committed Jul 2, 2022
1 parent bcc1a8b commit 8df460e
Show file tree
Hide file tree
Showing 4 changed files with 173 additions and 6 deletions.
21 changes: 19 additions & 2 deletions include/rtkThreeDCircularProjectionGeometry.h
Original file line number Diff line number Diff line change
Expand Up @@ -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);

Expand Down Expand Up @@ -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<double> &
GetCollimationUInf() const
Expand Down
56 changes: 52 additions & 4 deletions src/rtkThreeDCircularProjectionGeometry.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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];
Expand Down Expand Up @@ -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<double>::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<double>;
EulerType::Pointer euler = EulerType::New();
euler->SetComputeZYX(false); // ZXY order

// Extract angle using parent method without orthogonality check
euler->itk::MatrixOffsetTransformBase<double>::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)
Expand Down Expand Up @@ -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;
Expand All @@ -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;

Expand Down
1 change: 1 addition & 0 deletions test/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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}
Expand Down
101 changes: 101 additions & 0 deletions test/rtkparallelgeometryfrommatrixtest.cxx
Original file line number Diff line number Diff line change
@@ -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;
}

0 comments on commit 8df460e

Please sign in to comment.