fattori_g/reg23Topograms
0
0
Fork 0
Code Issues 1 Pull Requests Actions Packages Projects Releases Wiki Activity

R23 optimizes the isocentric IEC transform, internally the two internal transform are computed at every iteration.

Browse Source 
As a result, the feedback to regDialog is now already correct, no need for flips or sign change.
CalculateExternalTransform becomes obsolete... and has been commented out.
Methods for gui to get latest paremeters have been rewritten.
This commit is contained in:
Proton local user
2023-05-12 23:36:24 +02:00
parent 6ad0344791
commit d4c800dfcd
11 changed files with 586 additions and 474 deletions
Download Patch File Download Diff File Expand all files Collapse all files

23
reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp
View File

@ -517,6 +517,29 @@ RTGeometryMetaInformation
}
InternalTransformMetaInformation::InternalTransformMetaInformation(){
m_pCalProjCenter.Fill(0.);
m_pRTIsocenter.Fill(0.);
m_IECtoLPSDirs.SetIdentity();
}
void
InternalTransformMetaInformation
::PrintSelf(std::ostream& os, itk::Indent indent) const
{
Superclass::PrintSelf(os, indent);
}
InternalTransformMetaInformation
::~InternalTransformMetaInformation ()
{
}
R23MetaInformation::
R23MetaInformation (){

59
reg23Topograms/itkDTRrecon/DRTMetaInformation.h
View File

@ -514,6 +514,65 @@ private:
};
class InternalTransformMetaInformation :
public itk::Object{
public:
/** standard typedefs **/
typedef InternalTransformMetaInformation Self;
typedef itk::Object Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::Point<double, 3> PointType;
typedef itk::Matrix<double,3,3> TransformMatrixType;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro(InternalTransformMetaInformation, itk::Object);
/** object information streaming **/
void PrintSelf(std::ostream& os, itk::Indent indent) const;
itkSetMacro(pCalProjCenter,PointType);
itkSetMacro(pRTIsocenter,PointType);
itkSetMacro(IECtoLPSDirs,TransformMatrixType);
itkGetMacro(pCalProjCenter,PointType);
itkGetMacro(pRTIsocenter,PointType);
itkGetMacro(IECtoLPSDirs,TransformMatrixType);
protected:
PointType
m_pCalProjCenter,
m_pRTIsocenter;
TransformMatrixType
m_IECtoLPSDirs;
/** Default Constructor **/
InternalTransformMetaInformation ();
/** Default Destructor **/
virtual ~InternalTransformMetaInformation ();
private:
/** purposely not implemented **/
InternalTransformMetaInformation (const Self&);
/** purposely not implemented **/
void operator=(const Self&);
};
class R23MetaInformation :
public itk::Object{

10
reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h
View File

@ -177,6 +177,12 @@ public:
itkSetObjectMacro(TransformMetaInfo, R23MetaInformation);
itkGetConstObjectMacro(TransformMetaInfo, R23MetaInformation);
itkSetObjectMacro(internalMeta1, InternalTransformMetaInformation);
itkGetConstObjectMacro(internalMeta1, InternalTransformMetaInformation);
itkSetObjectMacro(internalMeta2, InternalTransformMetaInformation);
itkGetConstObjectMacro(internalMeta2, InternalTransformMetaInformation);
/** Set/Get the output filters. */
itkSetObjectMacro(Filter1, ChangeInformationFilterType);
itkGetConstObjectMacro(Filter1, ChangeInformationFilterType);
@ -259,6 +265,10 @@ private:
InterpolatorPointer m_Interpolator1;
InterpolatorPointer m_Interpolator2;
InternalTransformMetaInformation::Pointer
m_internalMeta1,
m_internalMeta2;
//TransformMetaInformation::Pointer
// m_TransformMetaInfo;

27
reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx
View File

@ -40,6 +40,8 @@ TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::TwoProjectionIm
m_Transform2 = nullptr; // has to be provided by the user.
m_IsocIECTransform = nullptr;
m_TransformMetaInfo = nullptr; // has to be provided by the user.
m_internalMeta1 = nullptr;
m_internalMeta2 = nullptr;
m_Interpolator1 = nullptr; // has to be provided by the user.
m_Interpolator2 = nullptr; // has to be provided by the user.
m_Metric = nullptr; // has to be provided by the user.
@ -144,6 +146,14 @@ void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::Initialize
itkExceptionMacro(<< "TransformMetaInfo is not present");
}
if (!m_internalMeta1) {
itkExceptionMacro(<< "internalMeta1 is not present");
}
if (!m_internalMeta2) {
itkExceptionMacro(<< "internalMeta1 is not present");
}
// Setup the metric
m_Metric->SetMovingImage(m_MovingImage);
m_Metric->SetFixedImage1(m_FixedImage1);
@ -152,6 +162,9 @@ void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::Initialize
m_Metric->SetTransform2(m_Transform2);
m_Metric->SetIsocTransf(m_IsocIECTransform);
m_Metric->SetTransformMetaInfo(m_TransformMetaInfo);
m_Metric->SetinternalMeta1(m_internalMeta1);
m_Metric->SetinternalMeta2(m_internalMeta2);
m_Metric->SetFilter1(m_Filter1);
m_Metric->SetFilter2(m_Filter2);
@ -206,7 +219,7 @@ void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::OnInternal
// if (m_CancelRequest)
// filter->SetAbortGenerateData(true); // may be handled by filter
std::cout << "OnInternalFilterProgressReceptor :: filter (TURE) " << std::endl;
std::cout << "OnInternalFilterProgressReceptor :: filter (TRUE) " << std::endl;
}
}
@ -217,22 +230,10 @@ template <typename TFixedImage, typename TMovingImage>
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::StartRegistration()
{
// typedef itk::MemberCommand<TwoProjectionImageRegistrationMethod> ITKCommandType;
// typedef ITKCommandType::Pointer MemberPointer ;
// ITKCommandType::Pointer cmd = ITKCommandType::New();
// itk::Command cmd= itk::Command::Pointer(this);
// cmd->SetCallbackFunction(this,
// (ITKCommandType::Pointer)&TwoProjectionImageRegistrationMethod::OnInternalFilterProgressReceptor);
itk::ProgressReporter progress(
this,
1, 1000,100);
// this->SetAbortGenerateData(true);
ParametersType empty(1);
empty.Fill(0.0);
try {

13
reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.h
View File

@ -40,6 +40,8 @@
#include <linux/limits.h>
#include <unistd.h>
#include <itkImageProcessorHelpers.h>
namespace itk {
/** \class gTwoImageToOneImageMetric
@ -194,6 +196,13 @@ public:
/** Get a pointer to the DRTGeometryMetaInfo. */
itkGetConstObjectMacro(TransformMetaInfo, R23MetaInformation);
itkSetObjectMacro(internalMeta1, InternalTransformMetaInformation);
itkGetConstObjectMacro(internalMeta1, InternalTransformMetaInformation);
itkSetObjectMacro(internalMeta2, InternalTransformMetaInformation);
itkGetConstObjectMacro(internalMeta2, InternalTransformMetaInformation);
/** Set the region over which the metric will be computed */
itkSetMacro(FixedImageRegion1, FixedImageRegionType);
@ -283,6 +292,10 @@ protected:
R23MetaInformation::Pointer
m_TransformMetaInfo;
InternalTransformMetaInformation::Pointer
m_internalMeta1,
m_internalMeta2;
ChangeInformationFilterPointer
m_Filter1,
m_Filter2;

104
reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx
View File

@ -32,6 +32,8 @@ gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::gTwoImageToOneImageMetric(
m_Transform1 = nullptr; // has to be provided by the user.
m_Transform2 = nullptr; // has to be provided by the user.
m_TransformMetaInfo = nullptr; // has to be provided by the user.
m_internalMeta1 = nullptr;
m_internalMeta2 = nullptr;
m_Interpolator1 = nullptr; // has to be provided by the user.
m_Interpolator2 = nullptr; // has to be provided by the user.
m_GradientImage = nullptr; // will receive the output of the filter;
@ -55,6 +57,13 @@ bool gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::SetTransformParameter
itkExceptionMacro(<< "Transform 2 has not been assigned");
}
if (!m_internalMeta1){
itkExceptionMacro(<< "m_internalMeta1 has not been assigned");
}
if (!m_internalMeta2){
itkExceptionMacro(<< "m_internalMeta2 has not been assigned");
}
//auto transformParameters = m_Transform1->GetParameters();
auto transformParameters = m_IsocTransf->GetParameters();
double TranslationAlongX = transformParameters[3];
@ -64,6 +73,8 @@ bool gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::SetTransformParameter
double RotationAlongY = transformParameters[1];
double RotationAlongZ = transformParameters[2];
//std::cout << "m_IsocTransf PARS: "<< transformParameters <<std::endl;
switch (m_TransformMetaInfo->GetDegreeOfFreedom()) {
case X_ONLY:
TranslationAlongX = parameters[0];
@ -102,33 +113,73 @@ bool gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::SetTransformParameter
std::cout << " Translation X = " << TranslationAlongX << " mm" << std::endl;
std::cout << " Translation Y = " << TranslationAlongY << " mm" << std::endl;
std::cout << " Translation Z = " << TranslationAlongZ << " mm" << std::endl;
std::cout << " Rotation Along X = " << RotationAlongX / dtr << " deg" << std::endl;
std::cout << " Rotation Along Y = " << RotationAlongY / dtr << " deg" << std::endl;
std::cout << " Rotation Along Z = " << RotationAlongZ / dtr << " deg" << std::endl;
std::cout << " Rotation Along X = " << RotationAlongX /*/ dtr*/ << " deg" << std::endl;
std::cout << " Rotation Along Y = " << RotationAlongY /*/ dtr*/ << " deg" << std::endl;
std::cout << " Rotation Along Z = " << RotationAlongZ /*/ dtr*/ << " deg" << std::endl;
std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
//HERE WE HAVE TO CALCULATE INTERNAL TRANSFORMS!!!
// GIOVANNI
transformParameters[0] = RotationAlongX;
transformParameters[1] = RotationAlongY;
transformParameters[2] = RotationAlongZ;
transformParameters[3] = TranslationAlongX;
transformParameters[4] = TranslationAlongY;
transformParameters[5] = TranslationAlongZ;
ImageType3D::PointType
pT;
pT[0] = TranslationAlongX;
pT[1] = TranslationAlongY;
pT[2] = TranslationAlongZ;
bool transformValid = (std::abs(TranslationAlongX) < m_MaxTranslation) &&
(std::abs(TranslationAlongY) < m_MaxTranslation) &&
(std::abs(TranslationAlongZ) < m_MaxTranslation);
ImageType3D::PointType
pR;
pR[0] = RotationAlongX;
pR[1] = RotationAlongY;
pR[2] = RotationAlongZ;
if (transformValid) {
m_Transform1->SetParameters(transformParameters);
m_Transform2->SetParameters(transformParameters);
} else {
std::cout << " Transform invalid, out of range!" << std::endl;
std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
}
TransformType::Pointer CurrTransform;
CurrTransform = CalculateInternalTransformV3(
pT,
pR,
m_internalMeta1->GetpCalProjCenter(),
m_internalMeta1->GetpRTIsocenter(),
m_internalMeta1->GetIECtoLPSDirs());
m_Transform1->SetIdentity();
m_Transform1->SetParameters(
CurrTransform->GetParameters());
m_Transform1->SetCenter(
m_internalMeta1->GetpCalProjCenter());
// std::cout << "m_Transform1 PARS: "<< m_Transform1->GetParameters()<<std::endl;
CurrTransform = CalculateInternalTransformV3(
pT,
pR,
m_internalMeta2->GetpCalProjCenter(),
m_internalMeta2->GetpRTIsocenter(),
m_internalMeta2->GetIECtoLPSDirs());
m_Transform2->SetIdentity();
m_Transform2->SetParameters(
CurrTransform->GetParameters());
m_Transform2->SetCenter(
m_internalMeta2->GetpCalProjCenter());
// std::cout << "m_Transform2 PARS: "<< m_Transform2->GetParameters()<<std::endl;
bool transformValid = true;
// transformParameters[0] = RotationAlongX;
// transformParameters[1] = RotationAlongY;
// transformParameters[2] = RotationAlongZ;
// transformParameters[3] = TranslationAlongX;
// transformParameters[4] = TranslationAlongY;
// transformParameters[5] = TranslationAlongZ;
// bool transformValid = (std::abs(TranslationAlongX) < m_MaxTranslation) &&
// (std::abs(TranslationAlongY) < m_MaxTranslation) &&
// (std::abs(TranslationAlongZ) < m_MaxTranslation);
// if (transformValid) {
// m_Transform1->SetParameters(transformParameters);
// m_Transform2->SetParameters(transformParameters);
// } else {
// std::cout << " Transform invalid, out of range!" << std::endl;
// std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
// }
return transformValid;
}
@ -136,18 +187,25 @@ template <typename TFixedImage, typename TMovingImage>
void gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::Initialize()
{
if (!m_Transform1) {
if (!m_IsocTransf) {
itkExceptionMacro(<< "Transform is not present");
}
if (!m_Transform1) {
itkExceptionMacro(<< "Transform is not present");
itkExceptionMacro(<< "Transform1 is not present");
}
if (!m_Transform2) {
itkExceptionMacro(<< "Transform is not present");
itkExceptionMacro(<< "Transform2 is not present");
}
if(!m_internalMeta1) {
itkExceptionMacro(<< "m_internalMeta1 is not present");
}
if(!m_internalMeta2) {
itkExceptionMacro(<< "m_internalMeta2 is not present");
}
if (!m_Interpolator1) {
itkExceptionMacro(<< "Interpolator1 is not present");
}
@ -244,7 +302,7 @@ unsigned int gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetNumberOfPa
break;
}
return m_Transform1->GetNumberOfParameters();
return m_IsocTransf->GetNumberOfParameters();
}
template <typename TFixedImage, typename TMovingImage>

613
reg23Topograms/itkDTRrecon/itkImageProcessor.cpp
View File

@ -68,6 +68,10 @@ itkImageProcessor::itkImageProcessor()
transform2 = TransformType::New();
transform2->SetIdentity();
m_InternalTransf1 = InternalTransformMetaInformation::New();
m_InternalTransf2 = InternalTransformMetaInformation::New();
interpolator1 = InterpolatorType::New();
interpolator2 = InterpolatorType::New();
@ -1281,174 +1285,33 @@ double itkImageProcessor::GetLocalizerDisplayWindowWidth(int iImg)
}
//itkImageProcessor::TransformType::Pointer
//itkImageProcessor::MapTransformToNewOrigin(
// ImageType3D::PointType m_COR, // Center of rotation for the proj geometry. this is my new origin.
// ImageType3D::PointType m_Translations,
// ImageType3D::PointType m_Rotations
// ){
// TransformType::Pointer InputTransform = TransformType::New();
// InputTransform->SetComputeZYX(true);
// InputTransform->SetIdentity();
// TransformType::OutputVectorType translation;
// translation[0] = m_Translations[0];
// translation[1] = m_Translations[1];
// translation[2] = m_Translations[2];
// InputTransform->SetTranslation(translation);
// const double dtr = (atan(1.0) * 4.0) / 180.0;
// InputTransform->SetRotation(
// dtr * m_Rotations[0],
// dtr * m_Rotations[1],
// dtr * m_Rotations[2]);
// ImageType3D::PointType m_TransformOrigin;
// m_TransformOrigin.Fill(0.);
// InputTransform->SetCenter(
// m_TransformOrigin );
// ImageType3D::PointType NewOriginTranslations =
// InputTransform->TransformPoint(m_COR);
// ImageType3D::PointType DeltaNewOrigin =
// NewOriginTranslations - m_COR;
// TransformType::Pointer m_OutputTransform =
// TransformType::New();
// m_OutputTransform ->SetComputeZYX(true);
// m_OutputTransform ->SetIdentity();
// translation[0] = DeltaNewOrigin[0];
// translation[1] = DeltaNewOrigin[1];
// translation[2] = DeltaNewOrigin[2];
// m_OutputTransform->SetTranslation(translation);
// m_OutputTransform->SetRotation(
// dtr * m_Rotations[0],
// dtr * m_Rotations[1],
// dtr * m_Rotations[2]);
// m_OutputTransform->SetCenter(m_COR);
// InputTransform = NULL;
// return m_OutputTransform;
//}
// This External User Transform thing need to be checked out
void itkImageProcessor::CalculateExternalUserTransform(TransformType::Pointer transform, DRTImageMetaInformation::Pointer imageMetaInfo)
{
//crude hack to get from internal transform in LPS to external transform in IEC.
//result is stored in m_TransformMetaInfo
//TODO: support for projectionTransformCenter and userTransformCenter which are not the same
//currentlz we support only pFakeIsoLPS situations.
InternalImageType::DirectionType LPStoIEC_Directions;
LPStoIEC_Directions = m_CTMetaInfo->GetLPS2IECDirections();
auto parameters = transform->GetParameters();
ImageType3D::PointType translationUser;
translationUser[0] = parameters[3];
translationUser[1] = parameters[4];
translationUser[2] = parameters[5];
ImageType3D::PointType rotationUser;
rotationUser[0] = parameters[0];
rotationUser[1] = parameters[1];
rotationUser[2] = parameters[2];
m_TransformMetaInfo->SetUserTranslations(LPStoIEC_Directions * translationUser);
m_TransformMetaInfo->SetUserRotations(LPStoIEC_Directions * rotationUser);
}
//itkImageProcessor::TransformType::Pointer
//itkImageProcessor::CalculateInternalTransformV2(
// ImageType3D::PointType m_TranslationOffset, //IEC
// ImageType3D::PointType m_RotationOffset, //IEC
// ImageType3D::PointType m_TranslationUser, //IEC
// ImageType3D::PointType m_RotationUser, //IEC
// ImageType3D::PointType m_CalibratedProjectionCenter, //LPS
// ImageType3D::PointType m_RTIsocenter, //LPS
// InternalImageType::DirectionType m_IECtoLPSDirections
// )
//void itkImageProcessor::CalculateExternalUserTransform(TransformType::Pointer transform, DRTImageMetaInformation::Pointer imageMetaInfo)
//{
// //crude hack to get from internal transform in LPS to external transform in IEC.
// //result is stored in m_TransformMetaInfo
// //Convert all inputs into LPS
// //TODO: support for projectionTransformCenter and userTransformCenter which are not the same
// //currentlz we support only pFakeIsoLPS situations.
// ImageType3D::PointType m_TOffsetLPS =
// m_IECtoLPSDirections * m_TranslationOffset;
// InternalImageType::DirectionType LPStoIEC_Directions;
// LPStoIEC_Directions = m_CTMetaInfo->GetLPS2IECDirections();
// ImageType3D::PointType m_ROffsetLPS =
// m_IECtoLPSDirections * m_RotationOffset;
// auto parameters = transform->GetParameters();
// ImageType3D::PointType m_TUserLPS =
// m_IECtoLPSDirections * m_TranslationUser;
// ImageType3D::PointType m_RUserLPS =
// m_IECtoLPSDirections * m_RotationUser;
// TransformType::OutputVectorType translation;
// translation[0] = m_TOffsetLPS[0] + m_TUserLPS[0];
// translation[1] = m_TOffsetLPS[1] + m_TUserLPS[1];
// translation[2] = m_TOffsetLPS[2] + m_TUserLPS[2];
// TransformType::OutputVectorType rotations;
// rotations[0] = m_ROffsetLPS[0] + m_RUserLPS[0];
// rotations[1] = m_ROffsetLPS[1] + m_RUserLPS[1];
// rotations[2] = m_ROffsetLPS[2] + m_RUserLPS[2];
// // Map offset to the projection center
// TransformType::Pointer m_outputTransform =
// MapTransformToNewOrigin (
// m_CalibratedProjectionCenter - m_RTIsocenter,
// translation,
// rotations
// );
// m_outputTransform->SetCenter(m_CalibratedProjectionCenter);
// return m_outputTransform;
//}
//itkImageProcessor::TransformType::Pointer
//itkImageProcessor::CalculateInternalTransformV3(
// ImageType3D::PointType m_Translation, //IEC
// ImageType3D::PointType m_Rotation, //IEC
// ImageType3D::PointType m_CalibratedProjectionCenter, //LPS
// ImageType3D::PointType m_RTIsocenter, //LPS
// InternalImageType::DirectionType m_IECtoLPSDirections
// )
//{
// //Convert all inputs into LPS
// ImageType3D::PointType m_TLPS =
// m_IECtoLPSDirections * m_Translation;
// ImageType3D::PointType m_RLPS =
// m_IECtoLPSDirections * m_Rotation;
// // Map offset to the projection center
// TransformType::Pointer m_outputTransform =
// MapTransformToNewOrigin (
// m_CalibratedProjectionCenter - m_RTIsocenter,
// m_TLPS,
// m_RLPS
// );
// m_outputTransform->SetCenter(m_CalibratedProjectionCenter);
// return m_outputTransform;
// ImageType3D::PointType translationUser;
// translationUser[0] = parameters[3];
// translationUser[1] = parameters[4];
// translationUser[2] = parameters[5];
// ImageType3D::PointType rotationUser;
// rotationUser[0] = parameters[0];
// rotationUser[1] = parameters[1];
// rotationUser[2] = parameters[2];
// m_TransformMetaInfo->SetUserTranslations(LPStoIEC_Directions * translationUser);
// m_TransformMetaInfo->SetUserRotations(LPStoIEC_Directions * rotationUser);
//}
@ -1521,103 +1384,19 @@ void itkImageProcessor::InitializeRegistration(
// registration->SetTransformMetaInfo(m_TransformMetaInfo);
registration->SetTransformMetaInfo(m_r23MetaInfo);
registration->SetinternalMeta1(m_InternalTransf1);
registration->SetinternalMeta2(m_InternalTransf2);
registration->SetFilter1(filter1);
registration->SetFilter2(filter2);
//CalculateInternalTransform(transform1, m_DRTImage1MetaInfo);
/********* BEGIN CALCULATION OF INTERNAL TRANSFORM FOR PROJ1 *********/
// ImageType3D::PointType ZeroPoint;
// ZeroPoint.Fill(0.);
// InternalImageType::DirectionType IECtoLPS_Directions;
// IECtoLPS_Directions =
// m_CTMetaInfo->GetLPS2IECDirections().GetTranspose();
// TransformType::Pointer CurrTransform;
// if(m_RTMetaInfo == NULL)
// {
// CurrTransform =
// CalculateInternalTransformV3(
// m_TransformMetaInfo->GetT(),
// m_TransformMetaInfo->GetR(),
// m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(),
// m_DRTImage1MetaInfo->GetProjectionOriginLPS(),
// IECtoLPS_Directions);
// } else {
// CurrTransform =
// CalculateInternalTransformV3(
// m_TransformMetaInfo->GetT(),
// m_TransformMetaInfo->GetR(),
// m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(),
// m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS(),
// IECtoLPS_Directions
// );
// }
// transform1->SetComputeZYX(true);
// transform1->SetIdentity();
// transform1->SetTranslation(
// CurrTransform->GetTranslation());
// transform1->SetRotation(
// CurrTransform->GetAngleX(),
// CurrTransform->GetAngleY(),
// CurrTransform->GetAngleZ()
// );
// transform1->SetCenter(
// m_DRTImage1MetaInfo->GetProjectionOriginLPSZero() );
// /********* END OF CALCULATE INTERNAL TRANSFORM FOR PROJ1 *********/
// //CalculateInternalTransform(transform2, m_DRTImage2MetaInfo);
// /********* BEGIN CALCULATION OF INTERNAL TRANSFORM FOR PROJ2 *********/
// if(m_RTMetaInfo == NULL)
// {
// CurrTransform =
// CalculateInternalTransformV3(
// m_TransformMetaInfo->GetT(),
// m_TransformMetaInfo->GetR(),
// m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(),
// m_DRTImage2MetaInfo->GetProjectionOriginLPS(),
// IECtoLPS_Directions);
// } else {
// CurrTransform =
// CalculateInternalTransformV3(
// m_TransformMetaInfo->GetT(),
// m_TransformMetaInfo->GetR(),
// m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(),
// m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS(),
// IECtoLPS_Directions
// );
// }
// transform2->SetComputeZYX(true);
// transform2->SetIdentity();
// transform2->SetTranslation(
// CurrTransform->GetTranslation());
// transform2->SetRotation(
// CurrTransform->GetAngleX(),
// CurrTransform->GetAngleY(),
// CurrTransform->GetAngleZ()
// );
// transform2->SetCenter(
// m_DRTImage2MetaInfo->GetProjectionOriginLPSZero() );
/********* END OF CALCULATE INTERNAL TRANSFORM FOR PROJ1 *********/
IsocTransf = TransformType::New();
IsocTransf->SetComputeZYX(true);
IsocTransf->SetIdentity();
IsocTransf->SetRotation(
m_TransformMetaInfo->GetR()[0],
m_TransformMetaInfo->GetR()[0],
m_TransformMetaInfo->GetR()[1],
m_TransformMetaInfo->GetR()[2]
);
@ -1690,8 +1469,8 @@ void itkImageProcessor::InitializeRegistration(
const int numberOfSteps = 25; //In each direction. Total number of steps is ((2*numberOfSteps+1))^3. For 25 -> 132651.
const double stepLength = 0.1;
auto parameters = transform1->GetParameters();
transform1->SetParameters(parameters);
//auto parameters = transform1->GetParameters();
//transform1->SetParameters(parameters);
ExhaustiveOptimizerType::StepsType steps(3);
steps[0] = numberOfSteps;
@ -1717,7 +1496,7 @@ int itkImageProcessor::StartRegistration(std::string extraInfo)
// TODO: Check if the registartion pipeline has been initialized
using ParametersType = RegistrationType::ParametersType;
auto startParameters = transform1->GetParameters();
//auto startParameters = transform1->GetParameters();
time_t t = time(0); // get time now
struct tm* now = localtime(&t);
@ -1784,54 +1563,106 @@ int itkImageProcessor::StartRegistration(std::string extraInfo)
fs.close();
m_OptmizerValue = optimizer->GetValue();
auto finalParameters = optimizer->GetCurrentPosition();
auto finalParameters = transform1->GetParameters();
std::cout<<"-->startParameters "<<startParameters<<std::endl;
std::cout<<"-->finalParameters1 "<<transform1->GetParameters()<<std::endl;
std::cout<<"-->finalParameters2 "<<transform2->GetParameters()<<std::endl;
//optimizer->GetCurrentPosition();
finalParameters = finalParameters - startParameters;
std::cout<<"startParameters "<<startParameters<<std::endl;
std::cout<<"finalParameters "<<finalParameters<<std::endl;
std::cout<<"finalParametersDelta "<<finalParameters<<std::endl;
std::cout<<"transform1 "<<transform1->GetCenter()<<std::endl;
std::cout<<"transform1 "<<transform1->GetTranslation()<<std::endl;
std::cout<<"transform2 "<<transform2->GetCenter()<<std::endl;
std::cout<<"transform2 "<<transform2->GetTranslation()<<std::endl;
std::cout<<"REG FINAL PARS: "<<finalParameters<<std::endl;
TransformType::Pointer offsetTransform = TransformType::New();
offsetTransform->SetParameters(finalParameters);
// auto finalParameters = transform1->GetParameters();
// std::cout<<"-->startParameters "<<startParameters<<std::endl;
// std::cout<<"-->finalParameters1 "<<transform1->GetParameters()<<std::endl;
// std::cout<<"-->finalParameters2 "<<transform2->GetParameters()<<std::endl;
CalculateExternalUserTransform(offsetTransform, m_DRTImage1MetaInfo);
// //optimizer->GetCurrentPosition();
// finalParameters = finalParameters - startParameters;
const double translationAlongX = finalParameters[3];
const double translationAlongY = finalParameters[4];
const double translationAlongZ = finalParameters[5];
const double rotationAlongX = finalParameters[0];
const double rotationAlongY = finalParameters[1];
const double rotationAlongZ = finalParameters[2];
// std::cout<<"startParameters "<<startParameters<<std::endl;
// std::cout<<"finalParameters "<<finalParameters<<std::endl;
// std::cout<<"finalParametersDelta "<<finalParameters<<std::endl;
// std::cout<<"transform1 "<<transform1->GetCenter()<<std::endl;
// std::cout<<"transform1 "<<transform1->GetTranslation()<<std::endl;
// std::cout<<"transform2 "<<transform2->GetCenter()<<std::endl;
// std::cout<<"transform2 "<<transform2->GetTranslation()<<std::endl;
const int numberOfIterations = optimizer->GetCurrentIteration();
std::cout << "Result = " << std::endl;
std::cout << " Rotation Along X = " << rotationAlongX / dtr << " deg" << std::endl;
std::cout << " Rotation Along Y = " << rotationAlongY / dtr << " deg" << std::endl;
std::cout << " Rotation Along Z = " << rotationAlongZ / dtr << " deg" << std::endl;
std::cout << " Translation X = " << translationAlongX << " mm" << std::endl;
std::cout << " Translation Y = " << translationAlongY << " mm" << std::endl;
std::cout << " Translation Z = " << translationAlongZ << " mm" << std::endl;
std::cout << " Number Of Iterations = " << numberOfIterations << std::endl;
std::cout << " Metric value = " << m_OptmizerValue << std::endl;
// TransformType::Pointer offsetTransform = TransformType::New();
// offsetTransform->SetParameters(finalParameters);
// CalculateExternalUserTransform(offsetTransform, m_DRTImage1MetaInfo);
// const double translationAlongX = finalParameters[3];
// const double translationAlongY = finalParameters[4];
// const double translationAlongZ = finalParameters[5];
// const double rotationAlongX = finalParameters[0];
// const double rotationAlongY = finalParameters[1];
// const double rotationAlongZ = finalParameters[2];
// const int numberOfIterations = optimizer->GetCurrentIteration();
// std::cout << "Result = " << std::endl;
// std::cout << " Rotation Along X = " << rotationAlongX / dtr << " deg" << std::endl;
// std::cout << " Rotation Along Y = " << rotationAlongY / dtr << " deg" << std::endl;
// std::cout << " Rotation Along Z = " << rotationAlongZ / dtr << " deg" << std::endl;
// std::cout << " Translation X = " << translationAlongX << " mm" << std::endl;
// std::cout << " Translation Y = " << translationAlongY << " mm" << std::endl;
// std::cout << " Translation Z = " << translationAlongZ << " mm" << std::endl;
// std::cout << " Number Of Iterations = " << numberOfIterations << std::endl;
// std::cout << " Metric value = " << m_OptmizerValue << std::endl;
std::cout << "#" << __COMPACT_PRETTY_FUNCTION__ << std::endl;
return 0;
}
Optimizer::ParametersType
itkImageProcessor::GetFinalR23Parameters(){
Optimizer::ParametersType pPars(6);
pPars.fill(0.);
if(optimizer == nullptr){
return pPars;
}
switch (m_r23MetaInfo->GetDegreeOfFreedom()) {
case THREE_DOF:
pPars[3] = optimizer->GetCurrentPosition()[0]
- m_TransformMetaInfo->GetHiddenTranslations()[0];
pPars[4] = optimizer->GetCurrentPosition()[1]
- m_TransformMetaInfo->GetHiddenTranslations()[1];
pPars[5] = optimizer->GetCurrentPosition()[2]
- m_TransformMetaInfo->GetHiddenTranslations()[2];
break;
case SIX_DOF:
pPars[0] = optimizer->GetCurrentPosition()[0]
- m_TransformMetaInfo->GetHiddenRotations()[0];
pPars[1] = optimizer->GetCurrentPosition()[1]
- m_TransformMetaInfo->GetHiddenRotations()[1];
pPars[2] = optimizer->GetCurrentPosition()[2]
- m_TransformMetaInfo->GetHiddenRotations()[2];
pPars[3] = optimizer->GetCurrentPosition()[3]
- m_TransformMetaInfo->GetHiddenTranslations()[0];
pPars[4] = optimizer->GetCurrentPosition()[4]
- m_TransformMetaInfo->GetHiddenTranslations()[1];
pPars[5] = optimizer->GetCurrentPosition()[5]
- m_TransformMetaInfo->GetHiddenTranslations()[2];
break;
default:
pPars.fill(0.);
break;
}
return
pPars;
}
void itkImageProcessor::InitializeProjector()
{
@ -1851,37 +1682,19 @@ void itkImageProcessor::InitializeProjector()
ZeroPoint.Fill(0.);
InternalImageType::DirectionType IECtoLPS_Directions;
// InternalImageType::DirectionType IECtoLPS_Directions;
IECtoLPS_Directions =
m_CTMetaInfo->GetLPS2IECDirections().GetTranspose();
// IECtoLPS_Directions =
// m_CTMetaInfo->GetLPS2IECDirections().GetTranspose();
TransformType::Pointer CurrTransform;
if(m_RTMetaInfo == NULL)
{
CurrTransform = CalculateInternalTransformV3(
m_TransformMetaInfo->GetT(),
m_TransformMetaInfo->GetR(),
m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(),
m_DRTImage1MetaInfo->GetProjectionOriginLPS(),
IECtoLPS_Directions);
} else {
CurrTransform =
CalculateInternalTransformV3(
m_TransformMetaInfo->GetT(),
m_TransformMetaInfo->GetR(),
m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(),
m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS(),
IECtoLPS_Directions
);
}
CurrTransform= CalculateInternalTransformV3(
m_TransformMetaInfo->GetT(),
m_TransformMetaInfo->GetR(),
m_InternalTransf1->GetpCalProjCenter(),
m_InternalTransf1->GetpRTIsocenter(),
m_InternalTransf1->GetIECtoLPSDirs());
transform1->SetComputeZYX(true);
@ -1912,30 +1725,12 @@ void itkImageProcessor::InitializeProjector()
interpolator1->SetTransform(transform1);
interpolator1->Initialize();
if(m_RTMetaInfo == NULL)
{
CurrTransform =
CalculateInternalTransformV3(
m_TransformMetaInfo->GetT(),
m_TransformMetaInfo->GetR(),
m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(),
m_DRTImage2MetaInfo->GetProjectionOriginLPS(),
IECtoLPS_Directions);
} else {
CurrTransform = CalculateInternalTransformV3(
m_TransformMetaInfo->GetT(),
m_TransformMetaInfo->GetR(),
m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(),
m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS(),
IECtoLPS_Directions
);
}
CurrTransform= CalculateInternalTransformV3(
m_TransformMetaInfo->GetT(),
m_TransformMetaInfo->GetR(),
m_InternalTransf2->GetpCalProjCenter(),
m_InternalTransf2->GetpRTIsocenter(),
m_InternalTransf2->GetIECtoLPSDirs());
transform2->SetComputeZYX(true);
@ -2258,6 +2053,29 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){
Std_DRT2LPS)
);
if(m_RTMetaInfo == NULL)
{
m_InternalTransf1->SetpCalProjCenter(
m_DRTImage1MetaInfo->GetProjectionOriginLPSZero());
m_InternalTransf1->SetpRTIsocenter(
m_DRTImage1MetaInfo->GetProjectionOriginLPS());
InternalImageType::DirectionType IECtoLPS_Directions;
IECtoLPS_Directions =
m_CTMetaInfo->GetLPS2IECDirections().GetTranspose();
m_InternalTransf1->SetIECtoLPSDirs(IECtoLPS_Directions);
} else {
m_InternalTransf1->SetpCalProjCenter(
m_DRTImage1MetaInfo->GetProjectionOriginLPSZero());
m_InternalTransf1->SetpRTIsocenter(
m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS());
InternalImageType::DirectionType IECtoLPS_Directions;
IECtoLPS_Directions =
m_CTMetaInfo->GetLPS2IECDirections().GetTranspose();
m_InternalTransf1->SetIECtoLPSDirs(IECtoLPS_Directions);
}
// SECOND
@ -2347,6 +2165,29 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){
Std_DRT2LPS)
);
if(m_RTMetaInfo == NULL)
{
m_InternalTransf2->SetpCalProjCenter(
m_DRTImage2MetaInfo->GetProjectionOriginLPSZero());
m_InternalTransf2->SetpRTIsocenter(
m_DRTImage2MetaInfo->GetProjectionOriginLPS());
InternalImageType::DirectionType IECtoLPS_Directions;
IECtoLPS_Directions =
m_CTMetaInfo->GetLPS2IECDirections().GetTranspose();
m_InternalTransf2->SetIECtoLPSDirs(IECtoLPS_Directions);
} else {
m_InternalTransf2->SetpCalProjCenter(
m_DRTImage2MetaInfo->GetProjectionOriginLPSZero());
m_InternalTransf2->SetpRTIsocenter(
m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS());
InternalImageType::DirectionType IECtoLPS_Directions;
IECtoLPS_Directions =
m_CTMetaInfo->GetLPS2IECDirections().GetTranspose();
m_InternalTransf2->SetIECtoLPSDirs(IECtoLPS_Directions);
}
}
@ -2398,35 +2239,13 @@ void itkImageProcessor::GetProjectionImages(){
transform1->SetComputeZYX(true);
transform1->SetIdentity();
InternalImageType::DirectionType IECtoLPS_Directions;
IECtoLPS_Directions =
m_CTMetaInfo->GetLPS2IECDirections().GetTranspose();
TransformType::Pointer CurrTransform;
if(m_RTMetaInfo == NULL)
{
CurrTransform =
CalculateInternalTransformV3(
m_TransformMetaInfo->GetT(),
m_TransformMetaInfo->GetR(),
m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(),
m_DRTImage1MetaInfo->GetProjectionOriginLPS(),
IECtoLPS_Directions);
} else {
CurrTransform =
CalculateInternalTransformV3(
m_TransformMetaInfo->GetT(),
m_TransformMetaInfo->GetR(),
m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(),
m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS(),
IECtoLPS_Directions
);
}
CurrTransform= CalculateInternalTransformV3(
m_TransformMetaInfo->GetT(),
m_TransformMetaInfo->GetR(),
m_InternalTransf1->GetpCalProjCenter(),
m_InternalTransf1->GetpRTIsocenter(),
m_InternalTransf1->GetIECtoLPSDirs());
transform1->SetComputeZYX(true);
@ -2460,32 +2279,12 @@ void itkImageProcessor::GetProjectionImages(){
transform2->SetComputeZYX(true);
transform2->SetIdentity();
if(m_RTMetaInfo == NULL)
{
CurrTransform =
CalculateInternalTransformV3(
m_TransformMetaInfo->GetT(),
m_TransformMetaInfo->GetR(),
m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(),
m_DRTImage2MetaInfo->GetProjectionOriginLPS(),
IECtoLPS_Directions);
} else {
CurrTransform = CalculateInternalTransformV3(
m_TransformMetaInfo->GetT(),
m_TransformMetaInfo->GetR(),
m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(),
m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS(),
IECtoLPS_Directions
);
}
CurrTransform = CalculateInternalTransformV3(
m_TransformMetaInfo->GetT(),
m_TransformMetaInfo->GetR(),
m_InternalTransf2->GetpCalProjCenter(),
m_InternalTransf2->GetpRTIsocenter(),
m_InternalTransf2->GetIECtoLPSDirs());
transform2->SetComputeZYX(true);
@ -3059,26 +2858,26 @@ double* itkImageProcessor::GetTransformParameters(){
return dTransfParam;
}
// THESE ARE CALLED AT THE END OF REG
void itkImageProcessor::GetFinalTranslations(double& dX, double& dY, double& dZ)
{
ImageType3D::PointType userTranslations = m_TransformMetaInfo->GetUserTranslations();
//// THESE ARE CALLED AT THE END OF REG
//void itkImageProcessor::GetFinalTranslations(double& dX, double& dY, double& dZ)
//{
// ImageType3D::PointType userTranslations = m_TransformMetaInfo->GetUserTranslations();
dX = userTranslations[0];
dY = userTranslations[1];
dZ = userTranslations[2];
}
// dX = userTranslations[0];
// dY = userTranslations[1];
// dZ = userTranslations[2];
//}
// THESE ARE CALLED AT THE END OF REG
void itkImageProcessor::GetFinalRotations(double& dRX, double& dRY, double& dRZ)
{
ImageType3D::PointType userRotations = m_TransformMetaInfo->GetUserRotations();
//// THESE ARE CALLED AT THE END OF REG
//void itkImageProcessor::GetFinalRotations(double& dRX, double& dRY, double& dRZ)
//{
// ImageType3D::PointType userRotations = m_TransformMetaInfo->GetUserRotations();
dRX = userRotations[0];
dRY = userRotations[1];
dRZ = userRotations[2];
// dRX = userRotations[0];
// dRY = userRotations[1];
// dRZ = userRotations[2];
}
//}
double itkImageProcessor::GetOptimizerValue()
{

43
reg23Topograms/itkDTRrecon/itkImageProcessor.h
View File

@ -139,13 +139,17 @@ public:
void SetInitialTranslations(double, double, double);
void SetInitialRotations(double, double, double);
/** Get transform parameters for 3D Volume */
void GetFinalTranslations(double&, double&, double&);
void GetFinalRotations(double&, double&, double&);
/** Set user transform parameters for 3D Volume in LPS*/
void SetUserTranslationsLPS(double, double, double);
void SetUserRotationsLPS(double, double, double);
Optimizer::ParametersType
GetFinalR23Parameters();
/** Get transform parameters for 3D Volume */
// void GetFinalTranslations(double&, double&, double&);
// void GetFinalRotations(double&, double&, double&);
// /** Set user transform parameters for 3D Volume in LPS*/
// void SetUserTranslationsLPS(double, double, double);
// void SetUserRotationsLPS(double, double, double);
/** Initialize the registration pipeline*/
void InitializeRegistration(double, double, eDegreeOfFreedomType);
@ -238,6 +242,8 @@ private:
using InterpolatorType = itk::gSiddonJacobsRayCastInterpolateImageFunction<InternalImageType, double>;
using ResampleFilterType = itk::ResampleImageFilter<InternalImageType, InternalImageType>;
//using gTransformType = itk::gEuler3DTransform<double>;
/** Optimizer which tries to find the minimn (Powell Optimizer)*/
using AmoebaOptimizerType = itk::AmoebaOptimizer;
/** Optimizer which samples the whol space */
@ -266,23 +272,10 @@ private:
/** ITK to VTK filter */
using ITKtoVTKFilterType = itk::ImageToVTKImageFilter<OutputImageType>;
void
CalculateExternalUserTransform(
TransformType::Pointer transform,
DRTImageMetaInformation::Pointer imageMetaInfo);
/* Calculate the transform used in siddon.
* The isocentric transform is mapped to the calibrated center of projection */
// TransformType::Pointer
// CalculateInternalTransformV2(
// ImageType3D::PointType m_TranslationOffset,
// ImageType3D::PointType m_RotationOffset,
// ImageType3D::PointType m_TranslationUser,
// ImageType3D::PointType m_RotationUser,
// ImageType3D::PointType m_CalibratedProjectionCenter,
// ImageType3D::PointType m_RTIsocenter,
// InternalImageType::DirectionType m_IECtoLPSDirections
// );
// void
// CalculateExternalUserTransform(
// TransformType::Pointer transform,
// DRTImageMetaInformation::Pointer imageMetaInfo);
TransformType::Pointer
@ -419,6 +412,10 @@ private:
R23MetaInformation::Pointer
m_r23MetaInfo;
InternalTransformMetaInformation::Pointer
m_InternalTransf1,
m_InternalTransf2;
double m_OptmizerValue;
int m_MaxNumberOfIterations;

81
reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp
View File

@ -1,7 +1,78 @@
#include "itkImageProcessorHelpers.h"
itk::TransformType::Pointer
itk::MapTransformToNewOrigin(
namespace itk {
//template <typename TParametersValueType>
//gEuler3DTransform <TParametersValueType>::gEuler3DTransform()
//{
// //m_CalibPrjCenter = nullptr; // has to be provided by the user.
// //m_RTIso = nullptr; // has to be provided by the user.
// // m_IECtoLPSDirs =; // has to be provided by the user.
//}
//template <typename TParametersValueType>
//gEuler3DTransform<TParametersValueType>::gEuler3DTransform
//(const MatrixType & matrix, const OutputPointType & offset)
//{
//}
//template <typename TParametersValueType>
//gEuler3DTransform<TParametersValueType>::gEuler3DTransform
//(unsigned int parametersDimension)
// : Superclass(parametersDimension)
//{
//}
//template <typename TParametersValueType>
//void gEuler3DTransform <TParametersValueType>::SetParameters(const ParametersType & parameters)
//{
//// if(m_CalibPrjCenter != nullptr &&
//// m_RTIso != nullptr //&&
//// //m_IECtoLPSDirs != nullptr
//// ){
//// itk::Point<double,3> pT;
//// pT[0] = parameters[3];
//// pT[1] = parameters[4];
//// pT[2] = parameters[5];
//// itk::Point<double,3> pR;
//// pR[0] = parameters[0];
//// pR[1] = parameters[1];
//// pR[2] = parameters[2];
//// TransformType::Pointer m_OutputTransform =
//// CalculateInternalTransformV3(
//// pT,
//// pR,
//// m_CalibPrjCenter,
//// m_RTIso,
//// m_IECtoLPSDirs);
//// Superclass::SetParameters(m_OutputTransform->GetParameters());
//// }else{
// Superclass::SetParameters(parameters);
//// }
//}
//template <typename TParametersValueType>
//void gEuler3DTransform <TParametersValueType>::PrintSelf(std::ostream& os, Indent indent) const
//{
// Superclass::PrintSelf(os, indent);
// //os << indent << "ComputeGradient: " << static_cast<typename NumericTraits<bool>::PrintType>(m_ComputeGradient)
// // << std::endl;
//}
TransformType::Pointer
MapTransformToNewOrigin(
ImageType3D::PointType m_COR, // Center of rotation for the proj geometry. this is my new origin.
ImageType3D::PointType m_Translations,
ImageType3D::PointType m_Rotations
@ -60,8 +131,8 @@ itk::MapTransformToNewOrigin(
itk::TransformType::Pointer
itk::CalculateInternalTransformV3(
TransformType::Pointer
CalculateInternalTransformV3(
ImageType3D::PointType m_Translation, //IEC
ImageType3D::PointType m_Rotation, //IEC
ImageType3D::PointType m_CalibratedProjectionCenter, //LPS
@ -91,3 +162,5 @@ itk::CalculateInternalTransformV3(
return m_outputTransform;
}
}

73
reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h
View File

@ -4,10 +4,80 @@
#include "itkEuler3DTransform.h"
#include "itkImage.h"
#include "itkEuler3DTransform.h"
#include "itkCommand.h"
#include "itkObject.h"
#include "itkObjectFactory.h"
namespace itk
{
//template <typename TParametersValueType = double>
//class ITK_TEMPLATE_EXPORT gEuler3DTransform :
// public Euler3DTransform <TParametersValueType>
//{
//public:
// ITK_DISALLOW_COPY_AND_ASSIGN(gEuler3DTransform);
// /** Standard class type aliases. */
// using Self = gEuler3DTransform;
// using Superclass = Euler3DTransform<TParametersValueType>;
// using Pointer = SmartPointer<Self>;
// using ConstPointer = SmartPointer<const Self>;
// /** New macro for creation of through a Smart Pointer. */
// itkNewMacro(Self);
// /** Run-time type information (and related methods). */
// itkTypeMacro(gEuler3DTransform, Euler3DTransform);
// // static constexpr unsigned int ParametersDimension = 6;
// using MatrixType = typename Superclass::MatrixType;
// using OutputPointType = typename Superclass::OutputPointType;
// using ParametersType = typename Superclass::ParametersType;
// using PointType = typename itk::Point<double, 3>;
// /** Set/Get the transformation from a container of parameters
// * This is typically used by optimizers. There are 6 parameters. The first
// * three represent the angles to rotate around the coordinate axis, and the
// * last three represents the offset. */
// void
// SetParameters(const ParametersType & parameters) override;
// itkSetMacro (CalibPrjCenter, PointType);
// itkSetMacro (RTIso, PointType);
// itkSetMacro(IECtoLPSDirs,MatrixType);
//protected:
// gEuler3DTransform(const MatrixType & matrix, const OutputPointType & offset);
// gEuler3DTransform(unsigned int parametersDimension);
// gEuler3DTransform();
// ~gEuler3DTransform() override = default;
// PointType
// m_CalibPrjCenter,
// m_RTIso;
// MatrixType
// m_IECtoLPSDirs;
// void
// PrintSelf(std::ostream & os, Indent indent) const override;
//};
using TransformType = itk::Euler3DTransform<double>;
constexpr static unsigned int Dimension = 3;
using PixelType3D = short;
using InternalPixelType = float;
@ -15,7 +85,6 @@ using InternalPixelType = float;
using ImageType3D = itk::Image<PixelType3D, Dimension>;
using InternalImageType = itk::Image<InternalPixelType, Dimension>;
using TransformType = itk::Euler3DTransform<double>;
TransformType::Pointer
MapTransformToNewOrigin(
@ -33,6 +102,8 @@ TransformType::Pointer
ImageType3D::PointType m_RTIsocenter,
InternalImageType::DirectionType m_IECtoLPSDirections
);
}
#endif

14
reg23Topograms/itkDTRrecon/itkQtIterationUpdate.h
View File

@ -118,11 +118,17 @@ public:
std::cout.precision(oldprecision);
std::cout << "Position: " << optimizer->GetCurrentPosition() << std::endl;
// objIterUpdate->onIteration(
// optimizer->GetCurrentIteration()+1,
// optimizer->GetCurrentPosition()[0],
// optimizer->GetCurrentPosition()[2],
// -optimizer->GetCurrentPosition()[1]
// );
objIterUpdate->onIteration(
optimizer->GetCurrentIteration()+1,
optimizer->GetCurrentPosition()[0],
optimizer->GetCurrentPosition()[2],
-optimizer->GetCurrentPosition()[1]
optimizer->GetCurrentPosition()[1],
optimizer->GetCurrentPosition()[2]
);
}
return;
@ -172,7 +178,9 @@ public:
os->precision(std::numeric_limits<double>::digits10 + 2);
*os << optimizer->GetCurrentValue();
os->precision(oldprecision);
*os << "\t" << position[0] << "\t" << position[2] << "\t" << -position[1] << std::endl;
//*os << "\t" << position[0] << "\t" << position[2] << "\t" << -position[1] << std::endl;
*os << "\t" << position[0] << "\t" << position[1] << "\t" << position[2] << std::endl;
}
return;
}