Skip to content

Commit

Permalink
Merge pull request #525 from SimonRit/ImagingRing
Browse files Browse the repository at this point in the history 
 Ora geometry enhancement
  • Loading branch information
@SimonRit
SimonRit authored Dec 16, 2022
2 parents 432417f + 40a2fc6 commit cfcf604
Show file tree
Hide file tree
Showing 13 changed files with 187 additions and 24 deletions.
1 change: 1 addition & 0 deletions applications/rtkorageometry/rtkorageometry.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@ main(int argc, char * argv[])
rtk::OraGeometryReader::Pointer oraReader = rtk::OraGeometryReader::New();
oraReader->SetProjectionsFileNames(rtk::GetProjectionsFileNamesFromGgo(args_info));
oraReader->SetCollimationMargin(margin);
oraReader->SetOptiTrackObjectID(args_info.optitrack_arg);
TRY_AND_EXIT_ON_ITK_EXCEPTION(oraReader->UpdateOutputData())

// Write
Expand Down
1 change: 1 addition & 0 deletions applications/rtkorageometry/rtkorageometry.ggo
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@ option "verbose" v "Verbose execution" flag off
option "config" - "Config file" string no
option "output" o "Output file name" string yes
option "margin" m "Collimation margin (uinf, usup, vinf, vsup)" double no multiple default="0."
option "optitrack" - "OptiTrack object ID (unused by default)" int no default="-1"

section "Projections"
option "path" p "Path containing projections" string yes
Expand Down
6 changes: 6 additions & 0 deletions include/rtkOraGeometryReader.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,6 +89,11 @@ class RTK_EXPORT OraGeometryReader : public itk::LightProcessObject
itkGetMacro(CollimationMargin, MarginVectorType);
itkSetMacro(CollimationMargin, MarginVectorType);

/** Selected OptiTrack object ID. Default is -1, i.e., do not use the
* OptiTrack object ID. * */
itkGetMacro(OptiTrackObjectID, int);
itkSetMacro(OptiTrackObjectID, int);


protected:
OraGeometryReader()
Expand All @@ -104,6 +109,7 @@ class RTK_EXPORT OraGeometryReader : public itk::LightProcessObject
GeometryType::Pointer m_Geometry;
FileNamesContainer m_ProjectionsFileNames;
MarginVectorType m_CollimationMargin;
int m_OptiTrackObjectID{ -1 };
};

} // namespace rtk
Expand Down
3 changes: 3 additions & 0 deletions include/rtkOraXMLFileReader.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,6 +76,9 @@ class RTK_EXPORT OraXMLFileReader : public itk::XMLReader<itk::MetaDataDictionar
EncapsulateMatrix3x3(const char * metaName, const char * name);
void
EncapsulateDouble(const char * metaName, const char * name);
template <typename TValue>
void
EncapsulateVector(const char * metaName, const char * name);
void
EncapsulateString(const char * metaName, const char * name);

Expand Down
107 changes: 105 additions & 2 deletions src/rtkOraGeometryReader.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@
#include <itkImageIOBase.h>
#include <itkImageIOFactory.h>
#include <itkVersorRigid3DTransform.h>
#include <itkQuaternionRigidTransform.h>

namespace rtk
{
Expand All @@ -37,6 +38,8 @@ OraGeometryReader::GenerateData()
{
m_Geometry = GeometryType::New();
RegisterIOFactories();
itk::QuaternionRigidTransform<double>::Pointer firstQuaternionsX{ nullptr };
itk::Vector<double, 3> firstTranslation;
for (const std::string & projectionsFileName : m_ProjectionsFileNames)
{
itk::ImageIOBase::Pointer reader;
Expand All @@ -54,6 +57,8 @@ OraGeometryReader::GenerateData()
using MetaDataVectorType = itk::MetaDataObject<VectorType>;
using MetaDataMatrixType = itk::MetaDataObject<Matrix3x3Type>;
using MetaDataDoubleType = itk::MetaDataObject<double>;
using MetaDataVectorDoubleType = itk::MetaDataObject<std::vector<double>>;
using MetaDataVectorIntType = itk::MetaDataObject<std::vector<int>>;

// Source position
MetaDataVectorType * spMeta = dynamic_cast<MetaDataVectorType *>(dic["SourcePosition"].GetPointer());
Expand Down Expand Up @@ -103,7 +108,7 @@ OraGeometryReader::GenerateData()

// Ring tilt (only available in some versions)
MetaDataDoubleType * tiltLeftMeta = dynamic_cast<MetaDataDoubleType *>(dic["tiltleft_deg"].GetPointer());
if (tiltLeftMeta != nullptr)
if (tiltLeftMeta != nullptr && m_OptiTrackObjectID < 0)
{
double tiltLeft = tiltLeftMeta->GetMetaDataObjectValue();
MetaDataDoubleType * tiltRightMeta = dynamic_cast<MetaDataDoubleType *>(dic["tiltright_deg"].GetPointer());
Expand All @@ -127,8 +132,106 @@ OraGeometryReader::GenerateData()
v = tiltTransform->TransformVector(v);
}

// Ring yaw (only available in some versions)
MetaDataDoubleType * yawMeta = dynamic_cast<MetaDataDoubleType *>(dic["room_cs_yaw_deg"].GetPointer());
if (yawMeta != nullptr && m_OptiTrackObjectID < 0)
{
double yaw = yawMeta->GetMetaDataObjectValue();
auto tiltTransform = itk::VersorRigid3DTransform<double>::New();
const double deg2rad = std::atan(1.0) / 45.0;
tiltTransform->SetRotation(itk::MakeVector(0., 0., 1.), yaw * deg2rad);

// Set center of rotation
MetaDataDoubleType * yvecMeta =
dynamic_cast<MetaDataDoubleType *>(dic["ydistancebaseunitcs2imagingcs_cm"].GetPointer());
double yvec = yvecMeta->GetMetaDataObjectValue();
tiltTransform->SetCenter(itk::MakePoint(0., -10. * yvec, 0.));

sp = tiltTransform->TransformPoint(sp);
dp = tiltTransform->TransformPoint(dp);
u = tiltTransform->TransformVector(u);
v = tiltTransform->TransformVector(v);
}

// OptiTrack objects (objects tracked with infrared cameras)
if (m_OptiTrackObjectID >= 0)
{
// Find ID index of the OptiTrack object
MetaDataVectorIntType * idsMeta = dynamic_cast<MetaDataVectorIntType *>(dic["optitrack_object_ids"].GetPointer());
if (idsMeta == nullptr)
itkExceptionMacro("Could not find optitrack_object_ids in " << projectionsFileName);
const std::vector<int> ids = idsMeta->GetMetaDataObjectValue();
auto idIt = std::find(ids.begin(), ids.end(), m_OptiTrackObjectID);
int idIdx = idIt - ids.begin();

// Translation
MetaDataVectorDoubleType * posMeta =
dynamic_cast<MetaDataVectorDoubleType *>(dic["optitrack_positions"].GetPointer());
if (posMeta == nullptr)
itkExceptionMacro("Could not find optitrack_positions in " << projectionsFileName);
const std::vector<double> p = posMeta->GetMetaDataObjectValue();
if (p.size() < 3 * (idIdx + 1))
itkExceptionMacro("Not enough values in optitrack_positions of " << projectionsFileName);
itk::Vector<double, 3> translation = 10. * itk::MakeVector(p[idIdx * 3], p[idIdx * 3 + 1], p[idIdx * 3 + 2]);

// Rotation
MetaDataVectorDoubleType * rotMeta =
dynamic_cast<MetaDataVectorDoubleType *>(dic["optitrack_rotations"].GetPointer());
if (rotMeta == nullptr)
itkExceptionMacro("Could not find optitrack_rotations in " << projectionsFileName);
const std::vector<double> optitrackRotations = rotMeta->GetMetaDataObjectValue();
if (optitrackRotations.size() < 4 * (idIdx + 1))
itkExceptionMacro("Not enough values in optitrack_rotations of " << projectionsFileName);
auto quaternionsX = itk::QuaternionRigidTransform<double>::New();
itk::QuaternionRigidTransform<double>::ParametersType quaternionsXParam(7);
quaternionsXParam[3] = optitrackRotations[idIdx * 4];
quaternionsXParam[0] = optitrackRotations[idIdx * 4 + 1];
quaternionsXParam[1] = optitrackRotations[idIdx * 4 + 2];
quaternionsXParam[2] = optitrackRotations[idIdx * 4 + 3];
quaternionsXParam[4] = 0.;
quaternionsXParam[5] = 0.;
quaternionsXParam[6] = 0.;
quaternionsX->SetParameters(quaternionsXParam);

// Set center of rotation
MetaDataDoubleType * yvecMeta =
dynamic_cast<MetaDataDoubleType *>(dic["ydistancebaseunitcs2imagingcs_cm"].GetPointer());
double yvec = yvecMeta->GetMetaDataObjectValue();
MetaDataDoubleType * zvecMeta =
dynamic_cast<MetaDataDoubleType *>(dic["zdistancebaseunitcs2imagingcs_cm"].GetPointer());
double zvec = zvecMeta->GetMetaDataObjectValue();
if (firstQuaternionsX.GetPointer() == nullptr)
{
firstQuaternionsX = quaternionsX;
firstTranslation = translation;
}
else
{
itk::MatrixOffsetTransformBase<double, 3, 3>::InverseTransformBasePointer invQuaternionsX =
quaternionsX->GetInverseTransform();

sp = sp - translation;
dp = dp - translation;
sp = invQuaternionsX->TransformPoint(sp);
dp = invQuaternionsX->TransformPoint(dp);
u = invQuaternionsX->TransformVector(u);
v = invQuaternionsX->TransformVector(v);

sp = firstQuaternionsX->TransformPoint(sp);
dp = firstQuaternionsX->TransformPoint(dp);
u = firstQuaternionsX->TransformVector(u);
v = firstQuaternionsX->TransformVector(v);
sp = sp + firstTranslation;
dp = dp + firstTranslation;
}
}

// Got it, add to geometry
m_Geometry->AddProjection(sp, dp, u, v);
if (!m_Geometry->AddProjection(sp, dp, u, v))
{
itkWarningMacro("Could not add " << projectionsFileName << " with sp=" << sp << ", dp=" << dp << ", u=" << u
<< " and v=" << v);
}

// Now add the collimation
// longitudinalposition_cm
Expand Down
42 changes: 33 additions & 9 deletions src/rtkOraXMLFileReader.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,13 +23,13 @@

namespace rtk
{
OraXMLFileReader ::OraXMLFileReader()
OraXMLFileReader::OraXMLFileReader()
{
m_OutputObject = &m_Dictionary;
}

int
OraXMLFileReader ::CanReadFile(const char * name)
OraXMLFileReader::CanReadFile(const char * name)
{
if (!itksys::SystemTools::FileExists(name) || itksys::SystemTools::FileIsDirectory(name) ||
itksys::SystemTools::FileLength(name) == 0)
Expand All @@ -38,13 +38,13 @@ OraXMLFileReader ::CanReadFile(const char * name)
}

void
OraXMLFileReader ::StartElement(const char * itkNotUsed(name), const char ** itkNotUsed(atts))
OraXMLFileReader::StartElement(const char * itkNotUsed(name), const char ** itkNotUsed(atts))
{
m_CurCharacterData = "";
}

void
OraXMLFileReader ::EndElement(const char * name)
OraXMLFileReader::EndElement(const char * name)
{
EncapsulatePoint("SourcePosition", name);
EncapsulatePoint("Origin", name);
Expand All @@ -60,19 +60,23 @@ OraXMLFileReader ::EndElement(const char * name)
EncapsulateDouble("xrayy2_cm", name);
EncapsulateDouble("tiltleft_deg", name);
EncapsulateDouble("tiltright_deg", name);
EncapsulateDouble("room_cs_yaw_deg", name);
EncapsulateDouble("ydistancebaseunitcs2imagingcs_cm", name);
EncapsulateDouble("zdistancebaseunitcs2imagingcs_cm", name);
EncapsulateVector<int>("optitrack_object_ids", name);
EncapsulateVector<double>("optitrack_positions", name);
EncapsulateVector<double>("optitrack_rotations", name);
}

void
OraXMLFileReader ::CharacterDataHandler(const char * inData, int inLength)
OraXMLFileReader::CharacterDataHandler(const char * inData, int inLength)
{
for (int i = 0; i < inLength; i++)
m_CurCharacterData = m_CurCharacterData + inData[i];
}

void
OraXMLFileReader ::EncapsulatePoint(const char * metaName, const char * name)
OraXMLFileReader::EncapsulatePoint(const char * metaName, const char * name)
{
if (itksys::SystemTools::Strucmp(name, metaName) == 0)
{
Expand All @@ -89,7 +93,7 @@ OraXMLFileReader ::EncapsulatePoint(const char * metaName, const char * name)
}

void
OraXMLFileReader ::EncapsulateMatrix3x3(const char * metaName, const char * name)
OraXMLFileReader::EncapsulateMatrix3x3(const char * metaName, const char * name)
{
if (itksys::SystemTools::Strucmp(name, metaName) == 0)
{
Expand All @@ -109,7 +113,7 @@ OraXMLFileReader ::EncapsulateMatrix3x3(const char * metaName, const char * name
}

void
OraXMLFileReader ::EncapsulateDouble(const char * metaName, const char * name)
OraXMLFileReader::EncapsulateDouble(const char * metaName, const char * name)
{
if (itksys::SystemTools::Strucmp(name, metaName) == 0)
{
Expand All @@ -118,8 +122,28 @@ OraXMLFileReader ::EncapsulateDouble(const char * metaName, const char * name)
}
}

template <typename TValue>
void
OraXMLFileReader ::EncapsulateString(const char * metaName, const char * name)
OraXMLFileReader::EncapsulateVector(const char * metaName, const char * name)
{
if (itksys::SystemTools::Strucmp(name, metaName) == 0)
{
std::vector<TValue> v;
std::istringstream iss(m_CurCharacterData);
TValue d;
iss >> d;
while (!iss.fail())
{
v.push_back(d);
iss.ignore(1);
iss >> d;
}
itk::EncapsulateMetaData<std::vector<TValue>>(m_Dictionary, metaName, v);
}
}

void
OraXMLFileReader::EncapsulateString(const char * metaName, const char * name)
{
if (itksys::SystemTools::Strucmp(name, metaName) == 0)
{
Expand Down
1 change: 1 addition & 0 deletions test/Baseline/Ora/geometry_optitrack.xml.sha512
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
5bab8df7a936505c56a94881b266791e65d4c4fde6a3ea96043efdb7f8a2a37ecb482d1b7158e843be2d0ee16878df4a3c6c97e7f61369b29334d86422b4f06c
1 change: 0 additions & 1 deletion test/Baseline/Ora/geometry_tilt.xml.sha512
View file

This file was deleted.

1 change: 1 addition & 0 deletions test/Baseline/Ora/geometry_yawtilt.xml.sha512
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
9f56034b740d0765a507ee3a0d8dc7467cb569f35f80d737b11077d5ffd0f426b3e8d079e73658e19d278890a89ab746a6b85fdb7c86a035fbfce863947e2bbd
5 changes: 3 additions & 2 deletions test/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -290,8 +290,9 @@ rtk_add_test(rtkParallelGeometryFromMatrixTest rtkparallelgeometryfrommatrixtest
rtk_add_test(rtkOraTest rtkoratest.cxx
DATA{Input/Ora/0_afterLog.ora.xml,0_afterLog.mhd,0_afterLog.raw}
DATA{Baseline/Ora/geometry.xml}
DATA{Input/Ora/084183_20211217170607335.ora.xml,084183_20211217170607335.mhd}
DATA{Baseline/Ora/geometry_tilt.xml}
DATA{Input/Ora/2006137_20220918183246810.ora.xml,2006137_20220918183246810.mhd}
DATA{Baseline/Ora/geometry_yawtilt.xml}
DATA{Baseline/Ora/geometry_optitrack.xml}
DATA{Baseline/Ora/attenuation.mha})

rtk_add_test(rtkBioscanTest rtkbioscantest.cxx
Expand Down
1 change: 1 addition & 0 deletions test/Input/Ora/2006137_20220918183246810.mhd.sha512
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
7dbe56d87d6026dfa9c350fe2f1f6e9558697044e5acf90c6dcb2a31eb106bb3b24a6a5ae39a32aa9ebe4532af4ec8e5b26cc7f3ad49c56e0eba289a2d39d6c8
1 change: 1 addition & 0 deletions test/Input/Ora/2006137_20220918183246810.ora.xml.sha512
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
20f97eb60e40586308bf576e4873e9d534aedbf7148151d2babb30a1556db61c76ba5c0c5c0fa29bb7c8bd93ad94761ddfe7cee7565006ed2fae15d6add2f4d3
41 changes: 31 additions & 10 deletions test/rtkoratest.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,9 @@ main(int argc, char * argv[])
{
std::cerr << "Usage: " << std::endl;
std::cerr << argv[0] << "oraGeometry.xml refGeometry.xml"
<< " oraGeometry_tilt.xml refGeometry_tilt.xml"
<< " oraGeometry_yawtilt.xml refGeometry_yawtilt.xml"
<< " oraGeometry_yaw.xml refGeometry_yaw.xml"
<< " refGeometry_optitrack.xml"
<< " reference.mha" << std::endl;
return EXIT_FAILURE;
}
Expand All @@ -48,23 +50,42 @@ main(int argc, char * argv[])
// Check geometries
CheckGeometries(geoTargReader->GetGeometry(), geoRefReader->GetOutputObject());

std::cout << "Testing geometry with tilt..." << std::endl;
std::cout << "Testing geometry with tilt and yaw..." << std::endl;

// Ora geometry
std::vector<std::string> filenames_tilt;
filenames_tilt.emplace_back(argv[3]);
rtk::OraGeometryReader::Pointer geoTargReader_tilt;
geoTargReader_tilt = rtk::OraGeometryReader::New();
geoTargReader_tilt->SetProjectionsFileNames(filenames_tilt);
TRY_AND_EXIT_ON_ITK_EXCEPTION(geoTargReader_tilt->UpdateOutputData());
std::vector<std::string> filenames_yawtilt;
filenames_yawtilt.emplace_back(argv[3]);
rtk::OraGeometryReader::Pointer geoTargReader_yawtilt;
geoTargReader_yawtilt = rtk::OraGeometryReader::New();
geoTargReader_yawtilt->SetProjectionsFileNames(filenames_yawtilt);
TRY_AND_EXIT_ON_ITK_EXCEPTION(geoTargReader_yawtilt->UpdateOutputData());

// Reference geometry
geoRefReader = rtk::ThreeDCircularProjectionGeometryXMLFileReader::New();
geoRefReader->SetFilename(argv[4]);
TRY_AND_EXIT_ON_ITK_EXCEPTION(geoRefReader->GenerateOutputInformation())

// Check geometries
CheckGeometries(geoTargReader_tilt->GetGeometry(), geoRefReader->GetOutputObject());
CheckGeometries(geoTargReader_yawtilt->GetGeometry(), geoRefReader->GetOutputObject());

std::cout << "Testing geometry with optitrack..." << std::endl;

// Ora geometry
std::vector<std::string> filenames_opti;
filenames_opti.emplace_back(argv[3]);
rtk::OraGeometryReader::Pointer geoTargReader_opti;
geoTargReader_opti = rtk::OraGeometryReader::New();
geoTargReader_opti->SetProjectionsFileNames(filenames_opti);
geoTargReader_opti->SetOptiTrackObjectID(2);
TRY_AND_EXIT_ON_ITK_EXCEPTION(geoTargReader_opti->UpdateOutputData());

// Reference geometry
geoRefReader = rtk::ThreeDCircularProjectionGeometryXMLFileReader::New();
geoRefReader->SetFilename(argv[5]);
TRY_AND_EXIT_ON_ITK_EXCEPTION(geoRefReader->GenerateOutputInformation())

// Check geometries
CheckGeometries(geoTargReader_opti->GetGeometry(), geoRefReader->GetOutputObject());

// ******* COMPARING projections *******
std::cout << "Testing attenuation conversion..." << std::endl;
Expand Down Expand Up @@ -98,7 +119,7 @@ main(int argc, char * argv[])
// Reference projections reader
ReaderType::Pointer readerRef = ReaderType::New();
filenames.clear();
filenames.emplace_back(argv[5]);
filenames.emplace_back(argv[6]);
readerRef->SetFileNames(filenames);
TRY_AND_EXIT_ON_ITK_EXCEPTION(readerRef->Update());

Expand Down

0 comments on commit cfcf604

Please sign in to comment.