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

a number of autoREG and geometry fixes

Browse Source 
- automatic registration actually makes use of ROI. Before were ignored.
- cleanup of user and hidden transform and pass to UI after registration.
- change of internal image interpolator used in autoREG. Now autoREG performs well.
- general number of working units variable.
- GetFinalRotations, GetFinalTranslations and GetTransformParameters become obsolete.
- GetFinalR23Parameters and GetCompleteIsocentricTransform and GetUserIsocentricTransform are provided.
- only GetCompleteIsocentricTransform is missing implementation.
This commit is contained in:
Proton local user
2023-05-15 17:02:18 +02:00
parent d4c800dfcd
commit 9b616b3b55
8 changed files with 250 additions and 318 deletions
Download Patch File Download Diff File Expand all files Collapse all files

58
reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx
View File

@ -93,33 +93,61 @@ MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValue()
itkExceptionMacro(<< "Fixed filter1 has not been assigned");
}
if (!this->m_Transform1) {
itkExceptionMacro(<< "m_Transform1 has not been assigned");
}
if (!this->m_Transform2) {
itkExceptionMacro(<< "m_Transform2 has not been assigned");
}
std::cout<< "----- itkMutualInformationTwoImage GetValue -----" <<std::endl;
std::cout<< "Mutual GetValue T1 :: center : " <<
this->m_Transform1.GetPointer()->GetCenter() <<std::endl;
std::cout<< "Mutual GetValue T1 :: Pars : " <<
this->m_Transform1.GetPointer()->GetParameters() <<std::endl;
std::cout<< "--------------------------------" <<std::endl;
this->m_Interpolator1->SetTransform(this->m_Transform1);
this->m_Interpolator1->Initialize();
this->m_Filter1->Update();
using InternalMetricType = itk::MattesMutualInformationImageToImageMetric<TFixedImage, TMovingImage>;
using NearestNeighborType = itk::NearestNeighborInterpolateImageFunction<TMovingImage, double>;
/*
* NearestNeighborType in our images results in not resolving for the cranio-caudal direction.
* in LPS-Z we are very much bound to the CT resolution, therefore if we take the Nearest Neighbor
* we simply do not explore that direction... and will always have an uncertainty of 1 slice thickness.
* */
//using NearestNeighborType = itk::NearestNeighborInterpolateImageFunction<TMovingImage, double>;
/*
* LinearInterpolateImageFunction performs better for us, performance degradation not noticeable.
* */
using LinearInterpType = itk::LinearInterpolateImageFunction<TMovingImage, double>;
//We need to set transform and parameters.
auto dummyTransform = TransformType::New();
auto dummyParameters = dummyTransform->GetParameters();
auto dummyInterpolator1 = NearestNeighborType::New();
auto dummyInterpolator1 = LinearInterpType::New();
auto metric1 = InternalMetricType::New();
metric1->SetTransform(dummyTransform);
metric1->SetTransformParameters(dummyParameters);
metric1->SetInterpolator(dummyInterpolator1);
auto fixedImageRegion1 = fixedImage1->GetBufferedRegion();
metric1->SetFixedImageRegion(fixedImageRegion1);
//auto fixedImageRegion1 = fixedImage1->GetBufferedRegion();
//metric1->SetFixedImageRegion(fixedImageRegion1);
/* We get the fixed image region from the parent class... */
metric1->SetFixedImageRegion(Superclass::GetFixedImageRegion1());
std::cout<< "-----> Mutual :: fixedImageRegion1: "<< metric1->GetFixedImageRegion() << std::endl;
auto movingImageRegion = this->m_Filter1->GetOutput()->GetBufferedRegion();
unsigned int numberOfPixels = movingImageRegion.GetNumberOfPixels();
//auto numberOfSamples = static_cast<unsigned int>(numberOfPixels * 0.50); //100%
metric1->UseAllPixelsOn();
metric1->SetNumberOfHistogramBins(30);
// since we have a ROI, then we should not set allPixels to TRUE.
//metric1->UseAllPixelsOn();
metric1->SetNumberOfHistogramBins(50);
metric1->SetFixedImage(fixedImage1);
metric1->SetMovingImage(this->m_Filter1->GetOutput());
@ -132,26 +160,26 @@ MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValue()
this->m_Interpolator2->Initialize();
this->m_Filter2->Update();
auto dummyInterpolator2 = NearestNeighborType::New();
auto dummyInterpolator2 = LinearInterpType::New();
auto metric2 = InternalMetricType::New();
metric2->SetTransform(dummyTransform);
metric2->SetTransformParameters(dummyParameters);
metric2->SetInterpolator(dummyInterpolator2);
auto fixedImageRegion2 = fixedImage1->GetBufferedRegion();
metric2->SetFixedImageRegion(fixedImageRegion2);
// auto fixedImageRegion2 = fixedImage2->GetBufferedRegion();
// metric2->SetFixedImageRegion(fixedImageRegion2);
metric2->SetFixedImageRegion(Superclass::GetFixedImageRegion2());
std::cout<< "-----> Mutual :: fixedImageRegion2: "<< metric2->GetFixedImageRegion() << std::endl;
movingImageRegion = this->m_Filter2->GetOutput()->GetBufferedRegion();
numberOfPixels = movingImageRegion.GetNumberOfPixels();
//numberOfSamples = static_cast<unsigned int>(numberOfPixels * 0.50); //100%
//%metric2->SetNumberOfSpatialSamples(numberOfSamples);
metric2->UseAllPixelsOn();
metric2->SetNumberOfHistogramBins(30);
//metric2->SetNumberOfSpatialSamples(numberOfSamples);
//metric2->UseAllPixelsOn();
metric2->SetNumberOfHistogramBins(50);
metric2->SetFixedImage(fixedImage2);
metric2->SetMovingImage(this->m_Filter2->GetOutput());
metric2->Initialize();
auto measure2 = metric2->GetValue(dummyParameters);

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

@ -246,8 +246,8 @@ protected:
itkSetMacro(LastOptimizerParameters, ParametersType);
/** Entry-point for internal ITK filter observer. **/
void OnInternalFilterProgressReceptor(itk::Object *caller,
const itk::EventObject &event);
// void OnInternalFilterProgressReceptor(itk::Object *caller,
// const itk::EventObject &event);
private:
@ -269,9 +269,6 @@ private:
m_internalMeta1,
m_internalMeta2;
//TransformMetaInformation::Pointer
// m_TransformMetaInfo;
R23MetaInformation::Pointer
m_TransformMetaInfo;

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

@ -160,6 +160,7 @@ void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::Initialize
m_Metric->SetFixedImage2(m_FixedImage2);
m_Metric->SetTransform1(m_Transform1);
m_Metric->SetTransform2(m_Transform2);
m_Metric->SetIsocTransf(m_IsocIECTransform);
m_Metric->SetTransformMetaInfo(m_TransformMetaInfo);
m_Metric->SetinternalMeta1(m_internalMeta1);
@ -197,31 +198,31 @@ void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::Initialize
m_Optimizer->SetInitialPosition(m_InitialOptimizerParameters);
}
template <typename TFixedImage, typename TMovingImage>
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::OnInternalFilterProgressReceptor(
itk::Object *caller, const itk::EventObject &event)
{
itk::ProcessObject *filter = dynamic_cast<itk::ProcessObject *>(caller);
//template <typename TFixedImage, typename TMovingImage>
//void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::OnInternalFilterProgressReceptor(
// itk::Object *caller, const itk::EventObject &event)
//{
// itk::ProcessObject *filter = dynamic_cast<itk::ProcessObject *>(caller);
if(!itk::ProgressEvent().CheckEvent(&event) || !filter)
return;
// if(!itk::ProgressEvent().CheckEvent(&event) || !filter)
// return;
if (filter)
{
// double p = m_CurrentProgressStart;
// // filter->GetProgress() in [0;1]
// p += m_CurrentProgressSpan * filter->GetProgress();
// // NOTE: filter is buggy - not in [0;1] if multi-threading is activated!
// if (p > (m_CurrentProgressStart + m_CurrentProgressSpan))
// p = m_CurrentProgressStart + m_CurrentProgressSpan;
// TaskProgressInfo(m_CurrentProgressDirection, p);
// if (filter)
// {
//// double p = m_CurrentProgressStart;
//// // filter->GetProgress() in [0;1]
//// p += m_CurrentProgressSpan * filter->GetProgress();
//// // NOTE: filter is buggy - not in [0;1] if multi-threading is activated!
//// if (p > (m_CurrentProgressStart + m_CurrentProgressSpan))
//// p = m_CurrentProgressStart + m_CurrentProgressSpan;
//// TaskProgressInfo(m_CurrentProgressDirection, p);
// if (m_CancelRequest)
// filter->SetAbortGenerateData(true); // may be handled by filter
//// if (m_CancelRequest)
//// filter->SetAbortGenerateData(true); // may be handled by filter
std::cout << "OnInternalFilterProgressReceptor :: filter (TRUE) " << std::endl;
}
}
// std::cout << "OnInternalFilterProgressReceptor :: filter (TRUE) " << std::endl;
// }
//}
/*
* Starts the Registration Process

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

@ -56,13 +56,15 @@ bool gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::SetTransformParameter
if (!m_Transform2) {
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");
}
itkExceptionMacro(<< "m_internalMeta2 has not been assigned");
}
if (!m_IsocTransf){
itkExceptionMacro(<< "m_IsocTransf has not been assigned");
}
//auto transformParameters = m_Transform1->GetParameters();
auto transformParameters = m_IsocTransf->GetParameters();
@ -109,7 +111,7 @@ bool gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::SetTransformParameter
break;
}
std::cout << "New Transform Parameters = " << std::endl;
std::cout << "New Transform Parameters ISOC IEC = " << 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;
@ -119,67 +121,72 @@ bool gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::SetTransformParameter
std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
//HERE WE HAVE TO CALCULATE INTERNAL TRANSFORMS!!!
// GIOVANNI
// individual R and T components of the IsoC transform
ImageType3D::PointType
pT;
pT;
pT[0] = TranslationAlongX;
pT[1] = TranslationAlongY;
pT[2] = TranslationAlongZ;
ImageType3D::PointType
pR;
pR;
pR[0] = RotationAlongX;
pR[1] = RotationAlongY;
pR[2] = RotationAlongZ;
TransformType::Pointer CurrTransform;
CurrTransform = CalculateInternalTransformV3(
pT,
pR,
m_internalMeta1->GetpCalProjCenter(),
m_internalMeta1->GetpRTIsocenter(),
m_internalMeta1->GetIECtoLPSDirs());
// check if we are within tolerance
bool transformValid =
(std::abs(pT[0]) < m_MaxTranslation) &&
(std::abs(pT[1]) < m_MaxTranslation) &&
(std::abs(pT[2]) < m_MaxTranslation);
m_Transform1->SetIdentity();
m_Transform1->SetParameters(
CurrTransform->GetParameters());
m_Transform1->SetCenter(
if (transformValid) {
/*
Here we calculate the transform for eack projection corresponding to the isocentric
one we are optimizing.
This calculation takes into account calibrated center of projection for each view.
*/
// Calculate internal transform1
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;
std::cout<<"----- itkgTwoImageToOne SetTransformParameters -----"<<std::endl;
std::cout << "m_Transform1 :: Center: "<< m_Transform1->GetCenter()<<std::endl;
std::cout << "m_Transform1 :: Pars: "<< m_Transform1->GetParameters()<<std::endl;
std::cout<< "--------------------------------" <<std::endl;
// 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;
// }
// Calculate internal transform2
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());
} else {
std::cout << " Transform invalid, out of range!" << std::endl;
std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
}
return transformValid;
}
@ -203,8 +210,8 @@ void gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::Initialize()
itkExceptionMacro(<< "m_internalMeta1 is not present");
}
if(!m_internalMeta2) {
itkExceptionMacro(<< "m_internalMeta2 is not present");
}
itkExceptionMacro(<< "m_internalMeta2 is not present");
}
if (!m_Interpolator1) {
itkExceptionMacro(<< "Interpolator1 is not present");
@ -328,6 +335,7 @@ itk::Optimizer::ScalesType gTwoImageToOneImageMetric<TFixedImage, TMovingImage>:
weightings[2] = 1. / dtr;
break;
case SIX_DOF:
//G: ITHINKTHOSEARE FLIPPED!!
weightings[0] = 1.;
weightings[1] = 1.;
weightings[2] = 1.;
@ -348,7 +356,7 @@ typename gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::
ParametersType gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetParameters() const
{
// ParametersType parametersTransform = m_Transform1->GetParameters(); //angleX, angleY, angleZ, transX, transY, transZ
// ParametersType parametersTransform = m_Transform1->GetParameters(); //angleX, angleY, angleZ, transX, transY, transZ
ParametersType parametersTransform = m_IsocTransf->GetParameters(); //angleX, angleY, angleZ, transX, transY, transZ
ParametersType parameters(this->GetNumberOfParameters());
switch (m_TransformMetaInfo->GetDegreeOfFreedom()) {

181
reg23Topograms/itkDTRrecon/itkImageProcessor.cpp
View File

@ -61,6 +61,7 @@ const char *gGetStringValueFromTag(const gdcm::Tag& t, const gdcm::DataSet& ds)
itkImageProcessor::itkImageProcessor()
{
iNWUnits = 1;
// std::cout<<"itkImageProcessor::itkImageProcessor contructor"<<std::endl;
transform1 = TransformType::New();
@ -163,43 +164,6 @@ itkImageProcessor::itkImageProcessor()
m_DRTGeometryMetaInfo->SetProjectionCenterOffset1(Point3D);
m_DRTGeometryMetaInfo->SetProjectionCenterOffset2(Point3D);
// // TEST MAPPING
// {
// ImageType3D::PointType m_COR;
// ImageType3D::PointType m_TransformOrigin; // Origin of the user defined trasnform. likely isocenter.
// ImageType3D::PointType m_Translations;
// ImageType3D::PointType m_Rotations;
// TransformType::Pointer m_OutputTransform ;
// m_COR[0] = 4.70067;
// m_COR[1] = -775.183;
// m_COR[2] = 202.321 ;
// m_TransformOrigin.Fill(0.);
// m_Translations[0] = 85.33068531;
// m_Translations[1] = 61.87811729;
// m_Translations[2] = -1.13774466 ;
// m_Rotations[0] = -0.09535216;
// m_Rotations[1] = 0.15113885;
// m_Rotations[2] = -5.19765723;
// m_OutputTransform=
// this->CalculateIsocentricTransform(
// m_COR,
// m_TransformOrigin,
// m_Translations,
// m_Rotations
// );
// m_OutputTransform->Print(std::cout);
// };
// setup the Automatic Registration
metric = MetricType::New();
mimetric = MIMetricType::New();
@ -239,6 +203,13 @@ itkImageProcessor::~itkImageProcessor()
}
void itkImageProcessor::SetNumberOfWorkingUnits(int iN){
if(iN>0){
iNWUnits = iN;
}
}
void itkImageProcessor::PrintSelf(std::ostream& os, Indent indent) const
{
Superclass::PrintSelf(os, indent);
@ -482,7 +453,7 @@ int itkImageProcessor::load3DSerieFromFiles( std::vector<std::string> v_fnames){
imageSeriesReader3D->SetImageIO(dicomIO);
imageSeriesReader3D->SetFileNames(v_sortedFnames);
imageSeriesReader3D->SetNumberOfWorkUnits(8);
imageSeriesReader3D->SetNumberOfWorkUnits(iNWUnits);
imageSeriesReader3D->ForceOrthogonalDirectionOff(); // properly read CTs with gantry tilt
//imageSeriesReader3D->SetSpacingWarningRelThreshold();
@ -621,7 +592,7 @@ int itkImageProcessor::load3DSerieFromFolder(const char * pcDirName)
imageSeriesReader3D->SetImageIO(dicomIO);
imageSeriesReader3D->SetFileNames(fileNames);
imageSeriesReader3D->SetNumberOfWorkUnits(12);
imageSeriesReader3D->SetNumberOfWorkUnits(iNWUnits);
imageSeriesReader3D->ForceOrthogonalDirectionOff(); // properly read CTs with gantry tilt
imageSeriesReader3D->Update();
@ -1287,33 +1258,6 @@ double itkImageProcessor::GetLocalizerDisplayWindowWidth(int iImg)
// 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);
//}
void itkImageProcessor::InitializeRegistration(
double stepLength,
@ -1338,6 +1282,7 @@ void itkImageProcessor::InitializeRegistration(
metric->SetSubtractMean(true);
metric->SetMaxTranslation(maxTranslation);
// m_TransformMetaInfo->SetDegreeOfFreedom(dof);
m_r23MetaInfo->SetDegreeOfFreedom(dof);
@ -1364,6 +1309,8 @@ void itkImageProcessor::InitializeRegistration(
registration->SetFixedImageRegion1(roiAutoReg1);
registration->SetFixedImageRegion2(roiAutoReg2);
mimetric->SetFixedImageRegion1(roiAutoReg1);
mimetric->SetFixedImageRegion2(roiAutoReg2);
std::cout << "using ROIs for Auto Reg " << std::endl;
@ -1379,6 +1326,8 @@ void itkImageProcessor::InitializeRegistration(
registration->SetFixedImageRegion1(region1);
registration->SetFixedImageRegion2(region2);
mimetric->SetFixedImageRegion1(roiAutoReg1);
mimetric->SetFixedImageRegion2(roiAutoReg2);
}
@ -1620,6 +1569,11 @@ int itkImageProcessor::StartRegistration(std::string extraInfo)
Optimizer::ParametersType
itkImageProcessor::GetFinalR23Parameters(){
if (!m_TransformMetaInfo) {
itkExceptionMacro(<< "m_TransformMetaInfo not present");
}
Optimizer::ParametersType pPars(6);
pPars.fill(0.);
@ -1767,7 +1721,7 @@ void itkImageProcessor::InitializeProjector()
resampleFilter1->SetInput(
image3DIn);
resampleFilter1->SetDefaultPixelValue(0);
resampleFilter1->SetNumberOfWorkUnits(8);
resampleFilter1->SetNumberOfWorkUnits(iNWUnits);
// The parameters of interpolator1, such as ProjectionAngle and FocalPointToIsocenterDistance
// have been set before registration. Here we only need to replace the initial
@ -1786,7 +1740,7 @@ void itkImageProcessor::InitializeProjector()
resampleFilter2->SetInput(
image3DIn);
resampleFilter2->SetDefaultPixelValue(0);
resampleFilter2->SetNumberOfWorkUnits(8);
resampleFilter2->SetNumberOfWorkUnits(iNWUnits);
// The parameters of interpolator2, such as ProjectionAngle and FocalPointToIsocenterDistance
@ -2268,6 +2222,9 @@ void itkImageProcessor::GetProjectionImages(){
std::cout<< "CurrTransform Translations: "<<
CurrTransform->GetTranslation()<<std::endl;
std::cout<< "CurrTransform center: "<<
CurrTransform->GetCenter()<<std::endl;
//transform1 ->Print(std::cout);
// The parameters of interpolator1, such as ProjectionAngle and FocalPointToIsocenterDistance
@ -2838,26 +2795,87 @@ void itkImageProcessor::SetInitialRotations(double dX, double dY, double dZ)
}
// THIS IS READ FOR EXPORT. TO FIX.
double* itkImageProcessor::GetTransformParameters(){
Optimizer::ParametersType
itkImageProcessor::GetUserIsocentricTransform(){
ImageType3D::PointType Translations;
ImageType3D::PointType Rotations;
Optimizer::ParametersType
m_Pars(6);
m_Pars.fill(0.);
Translations = m_TransformMetaInfo->GetUserTranslations();
Rotations = m_TransformMetaInfo->GetUserRotations();
if(!m_TransformMetaInfo){
return m_Pars;
}
dTransfParam[0] = Translations[0];
dTransfParam[1] = Translations[1];
dTransfParam[2] = Translations[2];
m_Pars[0] = m_TransformMetaInfo->GetUserRotations()[0];
m_Pars[1] = m_TransformMetaInfo->GetUserRotations()[1];
m_Pars[2] = m_TransformMetaInfo->GetUserRotations()[2];
dTransfParam[3] = Rotations[0];
dTransfParam[4] = Rotations[1];
dTransfParam[5] = Rotations[2];
m_Pars[3] = m_TransformMetaInfo->GetUserTranslations()[0];
m_Pars[4] = m_TransformMetaInfo->GetUserTranslations()[1];
m_Pars[5] = m_TransformMetaInfo->GetUserTranslations()[2];
return
m_Pars;
return dTransfParam;
}
TransformType::Pointer
itkImageProcessor::GetCompleteIsocentricTransform
(ImageType3D::PointType TransformCenter){
TransformType::Pointer
regTransform = TransformType::New();
regTransform->SetComputeZYX(true);
regTransform->SetIdentity();
if(m_TransformMetaInfo == nullptr ||
m_CTMetaInfo == nullptr){
return regTransform;
}
std::cout<<" --- GetCompleteIsocentricTransform --- "<<std::endl;
std::cout<< "R: "<<
m_TransformMetaInfo->GetR()[0] <<" "<<
m_TransformMetaInfo->GetR()[1] <<" "<<
m_TransformMetaInfo->GetR()[2] <<" "<<std::endl;
std::cout<< "T: "<<
m_TransformMetaInfo->GetT()[0] <<" "<<
m_TransformMetaInfo->GetT()[1] <<" "<<
m_TransformMetaInfo->GetT()[2] <<" "<<std::endl;
std::cout<< "Import Offset: "<<
m_CTMetaInfo->GetImportOffset()[0]<<" "<<
m_CTMetaInfo->GetImportOffset()[1]<<" "<<
m_CTMetaInfo->GetImportOffset()[2]<<" "<<std::endl;
std::cout<<" --- --- --- "<<std::endl;
}
// THIS IS READ FOR EXPORT. TO FIX.
//double* itkImageProcessor::GetTransformParameters(){
// ImageType3D::PointType Translations;
// ImageType3D::PointType Rotations;
// Translations = m_TransformMetaInfo->GetUserTranslations();
// Rotations = m_TransformMetaInfo->GetUserRotations();
// dTransfParam[0] = Translations[0];
// dTransfParam[1] = Translations[1];
// dTransfParam[2] = Translations[2];
// dTransfParam[3] = Rotations[0];
// dTransfParam[4] = Rotations[1];
// dTransfParam[5] = Rotations[2];
// return dTransfParam;
//}
//// THESE ARE CALLED AT THE END OF REG
//void itkImageProcessor::GetFinalTranslations(double& dX, double& dY, double& dZ)
//{
@ -2963,6 +2981,7 @@ void itkImageProcessor::SetRegionFixed1(
std::cout << "itkImageProcessor " << std::endl;
std::cout << roiAutoReg1 << std::endl;
}
void itkImageProcessor::SetRegionFixed2(

34
reg23Topograms/itkDTRrecon/itkImageProcessor.h
View File

@ -78,6 +78,7 @@ public:
CommandIterationUpdate::Pointer
optimizerObserver;
void SetNumberOfWorkingUnits(int iN);
/** Input data load methods*/
int load3DSerieFromFolder(const char* );
@ -156,8 +157,21 @@ public:
/** Start the registration process*/
int StartRegistration(std::string extraInfo);
/** Get transform parameters for 3D Volume */
double* GetTransformParameters();
///** Get transform parameters for 3D Volume */
//double* GetTransformParameters();
/** Get the complete transform, likely for export to SRO.
* This includes user and hidden contributions.
* Transform center has to be specified, Zero would result
* in Isocentric transform with center at the RT Iso.
*/
TransformType::Pointer
GetCompleteIsocentricTransform(ImageType3D::PointType TransformCenter);
/** Get only the User component of the transform
* as parameters' list */
Optimizer::ParametersType
GetUserIsocentricTransform();
/** Initialize projection geometry */
void InitializeProjector();
@ -238,7 +252,7 @@ private:
/** The following lines define each of the components used in the
registration: The transform, optimizer, metric, interpolator and
the registration method itself. */
using TransformType = itk::Euler3DTransform<double>;
//using TransformType = itk::Euler3DTransform<double>;
using InterpolatorType = itk::gSiddonJacobsRayCastInterpolateImageFunction<InternalImageType, double>;
using ResampleFilterType = itk::ResampleImageFilter<InternalImageType, InternalImageType>;
@ -343,13 +357,6 @@ private:
);
// TransformType::Pointer
// MapTransformToNewOrigin(
// ImageType3D::PointType m_COR,
// ImageType3D::PointType m_Translations,
// ImageType3D::PointType m_Rotations
// );
double
CalcProjectionAngleLPS(
tPatOrientation pOrient,
@ -383,8 +390,8 @@ private:
* m_Projection1VTKTransform,
* m_Projection2VTKTransform;
/*Transformation Parameters */
double dTransfParam[6];
///*Transformation Parameters */
// double dTransfParam[6];
/**
* Many meta containers
@ -430,6 +437,9 @@ private:
bool m_UseMutualInformation;
bool m_UseDumptoFile;
int
iNWUnits;
};

68
reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp
View File

@ -2,74 +2,6 @@
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(

63
reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h
View File

@ -14,69 +14,6 @@ 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;