diff --git a/reg23Topograms/itkDTRrecon/CMakeLists.txt b/reg23Topograms/itkDTRrecon/CMakeLists.txt index fd530ce..e97344c 100644 --- a/reg23Topograms/itkDTRrecon/CMakeLists.txt +++ b/reg23Topograms/itkDTRrecon/CMakeLists.txt @@ -16,6 +16,7 @@ INCLUDE_DIRECTORIES(${GDCM_INCLUDE_DIRS}) SET(SRCS itkImageProcessor.cpp + itkImageProcessorHelpers.cpp vtkContourTopogramProjectionFilter.cxx DRTMetaInformation.cpp @@ -23,6 +24,7 @@ SET(SRCS SET(HDR itkImageProcessor.h + itkImageProcessorHelpers.h itkQtIterationUpdate.h itkgSiddonJacobsRayCastInterpolateImageFunction.h itkgSiddonJacobsRayCastInterpolateImageFunction.hxx diff --git a/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp b/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp index 55ada69..ff6ae1c 100644 --- a/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp +++ b/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp @@ -474,6 +474,8 @@ DRTProjectionGeometryImageMetaInformation(){ this->m_ProjectionCenter.Fill(0.); + this->m_UseRotationsForHiddenTransform = true; + } void @@ -517,25 +519,98 @@ 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 (){ + + m_DegreeOfFreedom = tDegreeOfFreedomEnum::UNKNOWN; + +} + + +void +R23MetaInformation +::PrintSelf(std::ostream& os, itk::Indent indent) const +{ + Superclass::PrintSelf(os, indent); +} + + +R23MetaInformation +::~R23MetaInformation () +{ + +} TransformMetaInformation :: TransformMetaInformation (){ - m_Translations.Fill(0.); + m_HiddenTranslations.Fill(0.); - m_Rotations.Fill(0.); + m_HiddenRotations.Fill(0.); m_UserRotations.Fill(0.); m_UserTranslations.Fill(0.); - m_ReferenceTransform.Fill(0.); + //m_ReferenceTransform.Fill(0.); - m_CurrentTransform.Fill(0.); + //m_CurrentTransform.Fill(0.); - m_DegreeOfFreedom = tDegreeOfFreedomEnum::UNKNOWN; + // m_DegreeOfFreedom = tDegreeOfFreedomEnum::UNKNOWN; + +} + +TransformMetaInformation::PointType +TransformMetaInformation::GetT(){ + + PointType + pT; + + pT[0] = m_UserTranslations[0] + m_HiddenTranslations[0]; + pT[1] = m_UserTranslations[1] + m_HiddenTranslations[1]; + pT[2] = m_UserTranslations[2] + m_HiddenTranslations[2]; + + return + pT; + +} + +TransformMetaInformation::PointType +TransformMetaInformation::TransformMetaInformation::GetR(){ + + PointType + pR; + + pR[0] = m_UserRotations[0] + m_HiddenRotations[0]; + pR[1] = m_UserRotations[1] + m_HiddenRotations[1]; + pR[2] = m_UserRotations[2] + m_HiddenRotations[2]; + + return + pR; } diff --git a/reg23Topograms/itkDTRrecon/DRTMetaInformation.h b/reg23Topograms/itkDTRrecon/DRTMetaInformation.h index 1f77f24..0e6e9cd 100644 --- a/reg23Topograms/itkDTRrecon/DRTMetaInformation.h +++ b/reg23Topograms/itkDTRrecon/DRTMetaInformation.h @@ -408,6 +408,10 @@ public: itkSetMacro(ProjectionCenter, PointType); itkGetMacro(ProjectionCenter, PointType); + itkSetMacro(UseRotationsForHiddenTransform, bool); + itkGetMacro(UseRotationsForHiddenTransform, bool); + + protected: double @@ -447,6 +451,9 @@ protected: m_ProjectionCenterOffset2; + bool + m_UseRotationsForHiddenTransform; + /** Default Constructor **/ DRTProjectionGeometryImageMetaInformation (); /** Default Destructor **/ @@ -514,6 +521,107 @@ private: }; + + + +class InternalTransformMetaInformation : + public itk::Object{ + +public: + /** standard typedefs **/ + typedef InternalTransformMetaInformation Self; + typedef itk::Object Superclass; + typedef itk::SmartPointer Pointer; + + typedef itk::Point PointType; + typedef itk::Matrix 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{ + +public: + /** standard typedefs **/ + typedef R23MetaInformation Self; + typedef itk::Object Superclass; + typedef itk::SmartPointer Pointer; + + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(R23MetaInformation, itk::Object); + + /** object information streaming **/ + void PrintSelf(std::ostream& os, itk::Indent indent) const; + + itkSetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum); + itkGetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum); + + +protected: + + tDegreeOfFreedomEnum + m_DegreeOfFreedom; + + /** Default Constructor **/ + R23MetaInformation (); + /** Default Destructor **/ + virtual ~R23MetaInformation (); + +private: + /** purposely not implemented **/ + R23MetaInformation (const Self&); + + /** purposely not implemented **/ + void operator=(const Self&); + +}; + class TransformMetaInformation : public itk::Object{ @@ -535,11 +643,11 @@ public: void PrintSelf(std::ostream& os, itk::Indent indent) const; - itkSetMacro(Translations,PointType); - itkGetMacro(Translations,PointType); + itkSetMacro(HiddenTranslations,PointType); + itkGetMacro(HiddenTranslations,PointType); - itkSetMacro(Rotations,PointType); - itkGetMacro(Rotations,PointType); + itkSetMacro(HiddenRotations,PointType); + itkGetMacro(HiddenRotations,PointType); itkSetMacro(UserTranslations,PointType); itkGetMacro(UserTranslations,PointType); @@ -547,29 +655,17 @@ public: itkSetMacro(UserRotations,PointType); itkGetMacro(UserRotations,PointType); - itkSetMacro(ReferenceTransform,TransformMatrixType); - itkGetMacro(ReferenceTransform,TransformMatrixType); - - itkSetMacro(CurrentTransform,TransformMatrixType); - itkGetMacro(CurrentTransform,TransformMatrixType); - - itkSetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum); - itkGetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum); + PointType GetT(); + PointType GetR(); protected: PointType - m_Translations, - m_Rotations, + m_HiddenTranslations, + m_HiddenRotations, m_UserTranslations, m_UserRotations; - TransformMatrixType - m_ReferenceTransform, - m_CurrentTransform; - - tDegreeOfFreedomEnum - m_DegreeOfFreedom; /** Default Constructor **/ TransformMetaInformation (); diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx b/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx index 089c7f9..2e18525 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx +++ b/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx @@ -93,33 +93,61 @@ MutualInformationTwoImageToOneImageMetric::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 -----" <m_Transform1.GetPointer()->GetCenter() <m_Transform1.GetPointer()->GetParameters() <m_Interpolator1->SetTransform(this->m_Transform1); this->m_Interpolator1->Initialize(); this->m_Filter1->Update(); using InternalMetricType = itk::MattesMutualInformationImageToImageMetric; - - using NearestNeighborType = itk::NearestNeighborInterpolateImageFunction; + /* + * 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; + /* + * LinearInterpolateImageFunction performs better for us, performance degradation not noticeable. + * */ + using LinearInterpType = itk::LinearInterpolateImageFunction; //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(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::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(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); @@ -163,13 +191,13 @@ MutualInformationTwoImageToOneImageMetric::GetValue() auto measure = (measure1 + measure2) / 2.0; - auto oldprecision = std::cout.precision(); - std::cout.precision(std::numeric_limits::digits10 + 2); - std::cout << "Measure1: " << measure1 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl; - std::cout << "Measure2: " << measure2 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl; - std::cout << "Measure: " << measure << std::endl; - std::cout << "=======================================================" << std::endl; - std::cout.precision(oldprecision); +// auto oldprecision = std::cout.precision(); +// std::cout.precision(std::numeric_limits::digits10 + 2); +// std::cout << "Measure1: " << measure1 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl; +// std::cout << "Measure2: " << measure2 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl; +// std::cout << "Measure: " << measure << std::endl; +// std::cout << "=======================================================" << std::endl; +// std::cout.precision(oldprecision); return measure; } diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h index 3dd6241..9ddf934 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h +++ b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h @@ -163,6 +163,10 @@ public: itkSetObjectMacro(Transform2, TransformType); itkGetConstObjectMacro(Transform2, TransformType); + /** Set/Get the IsocIECTransform. */ + itkSetObjectMacro(IsocIECTransform, TransformType); + itkGetConstObjectMacro(IsocIECTransform, TransformType); + /** Set/Get the Interpolators. */ itkSetObjectMacro(Interpolator1, InterpolatorType); itkSetObjectMacro(Interpolator2, InterpolatorType); @@ -170,8 +174,14 @@ public: itkGetConstObjectMacro(Interpolator2, InterpolatorType); /** Set/Get the Meta informations. */ - itkSetObjectMacro(TransformMetaInfo, TransformMetaInformation); - itkGetConstObjectMacro(TransformMetaInfo, TransformMetaInformation); + 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); @@ -236,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: @@ -248,12 +258,18 @@ private: FixedImageConstPointer m_FixedImage1; FixedImageConstPointer m_FixedImage2; + TransformPointer m_IsocIECTransform; + TransformPointer m_Transform1; TransformPointer m_Transform2; InterpolatorPointer m_Interpolator1; InterpolatorPointer m_Interpolator2; - TransformMetaInformation::Pointer + InternalTransformMetaInformation::Pointer + m_internalMeta1, + m_internalMeta2; + + R23MetaInformation::Pointer m_TransformMetaInfo; ChangeInformationFilterPointer diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx index b6f4b12..2c9b6e0 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx +++ b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx @@ -38,7 +38,10 @@ TwoProjectionImageRegistrationMethod::TwoProjectionIm m_MovingImage = nullptr; // has to be provided by the user. m_Transform1 = nullptr; // has to be provided by the user. 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. @@ -80,6 +83,7 @@ void TwoProjectionImageRegistrationMethod::SetFixedIm m_FixedImageRegionDefined2 = true; } + /* * Initialize by setting the interconnects between components. */ @@ -119,7 +123,8 @@ void TwoProjectionImageRegistrationMethod::Initialize // Connect the transform to the Decorator. auto* transformOutput = static_cast(this->ProcessObject::GetOutput(0)); - transformOutput->Set(m_Transform1.GetPointer()); + //transformOutput->Set(m_Transform1.GetPointer()); + transformOutput->Set(m_IsocIECTransform.GetPointer()); if (!m_Interpolator1) { itkExceptionMacro(<< "Interpolator1 is not present"); @@ -141,13 +146,26 @@ void TwoProjectionImageRegistrationMethod::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); 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); + m_Metric->SetinternalMeta2(m_internalMeta2); + m_Metric->SetFilter1(m_Filter1); m_Metric->SetFilter2(m_Filter2); @@ -174,36 +192,37 @@ void TwoProjectionImageRegistrationMethod::Initialize m_Optimizer->SetScales(m_Metric->GetWeightings()); m_Optimizer->SetCostFunction(m_Metric); - m_InitialOptimizerParameters = m_Metric->GetParameters(); + m_InitialOptimizerParameters = + m_Metric->GetParameters(); m_Optimizer->SetInitialPosition(m_InitialOptimizerParameters); } -template -void TwoProjectionImageRegistrationMethod::OnInternalFilterProgressReceptor( - itk::Object *caller, const itk::EventObject &event) -{ - itk::ProcessObject *filter = dynamic_cast(caller); +//template +//void TwoProjectionImageRegistrationMethod::OnInternalFilterProgressReceptor( +// itk::Object *caller, const itk::EventObject &event) +//{ +// itk::ProcessObject *filter = dynamic_cast(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 << "Porca Madonna " << std::endl; - } -} +// std::cout << "OnInternalFilterProgressReceptor :: filter (TRUE) " << std::endl; +// } +//} /* * Starts the Registration Process @@ -212,22 +231,10 @@ template void TwoProjectionImageRegistrationMethod::StartRegistration() { -// typedef itk::MemberCommand 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 { diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.h b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.h index 9cb93cb..2c08287 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.h +++ b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.h @@ -40,6 +40,8 @@ #include #include +#include + namespace itk { /** \class gTwoImageToOneImageMetric @@ -169,6 +171,10 @@ public: /** Get a pointer to the Transform2. */ itkGetConstObjectMacro(Transform2, TransformType); + itkSetObjectMacro(IsocTransf, TransformType); + itkGetObjectMacro(IsocTransf, TransformType); + + /** Connect the Interpolator. */ itkSetObjectMacro(Interpolator1, InterpolatorType); @@ -185,10 +191,17 @@ public: itkGetConstReferenceMacro(NumberOfPixelsCounted, unsigned long); /** Connect the DRTGeometryMetaInfo. */ - itkSetObjectMacro(TransformMetaInfo, TransformMetaInformation); + itkSetObjectMacro(TransformMetaInfo, R23MetaInformation); /** Get a pointer to the DRTGeometryMetaInfo. */ - itkGetConstObjectMacro(TransformMetaInfo, TransformMetaInformation); + 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); @@ -261,6 +274,8 @@ protected: mutable TransformPointer m_Transform1; mutable TransformPointer m_Transform2; + mutable TransformPointer m_IsocTransf; + InterpolatorPointer m_Interpolator1; InterpolatorPointer m_Interpolator2; @@ -274,9 +289,13 @@ protected: mutable FixedImageMaskPointer m_FixedImageMask2; mutable MovingImageMaskPointer m_MovingImageMask; - TransformMetaInformation::Pointer + R23MetaInformation::Pointer m_TransformMetaInfo; + InternalTransformMetaInformation::Pointer + m_internalMeta1, + m_internalMeta2; + ChangeInformationFilterPointer m_Filter1, m_Filter2; diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx index 90a9548..770ba54 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx +++ b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx @@ -32,6 +32,8 @@ gTwoImageToOneImageMetric::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; @@ -54,8 +56,18 @@ bool gTwoImageToOneImageMetric::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"); + } + if (!m_IsocTransf){ + itkExceptionMacro(<< "m_IsocTransf has not been assigned"); + } - auto transformParameters = m_Transform1->GetParameters(); + //auto transformParameters = m_Transform1->GetParameters(); + auto transformParameters = m_IsocTransf->GetParameters(); double TranslationAlongX = transformParameters[3]; double TranslationAlongY = transformParameters[4]; double TranslationAlongZ = transformParameters[5]; @@ -63,6 +75,8 @@ bool gTwoImageToOneImageMetric::SetTransformParameter double RotationAlongY = transformParameters[1]; double RotationAlongZ = transformParameters[2]; + //std::cout << "m_IsocTransf PARS: "<< transformParameters <GetDegreeOfFreedom()) { case X_ONLY: TranslationAlongX = parameters[0]; @@ -97,31 +111,82 @@ bool gTwoImageToOneImageMetric::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; -// 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; - 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); + // individual R and T components of the IsoC transform + ImageType3D::PointType + pT; + pT[0] = TranslationAlongX; + pT[1] = TranslationAlongY; + pT[2] = TranslationAlongZ; + + ImageType3D::PointType + pR; + pR[0] = RotationAlongX; + pR[1] = RotationAlongY; + pR[2] = RotationAlongZ; + + // 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); if (transformValid) { - m_Transform1->SetParameters(transformParameters); - m_Transform2->SetParameters(transformParameters); + /* + 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<<"----- itkgTwoImageToOne SetTransformParameters -----"<GetCenter()<GetParameters()<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; } @@ -129,16 +194,23 @@ template void gTwoImageToOneImageMetric::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) { @@ -237,7 +309,7 @@ unsigned int gTwoImageToOneImageMetric::GetNumberOfPa break; } - return m_Transform1->GetNumberOfParameters(); + return m_IsocTransf->GetNumberOfParameters(); } template @@ -263,6 +335,7 @@ itk::Optimizer::ScalesType gTwoImageToOneImageMetric: weightings[2] = 1. / dtr; break; case SIX_DOF: + //G: ITHINKTHOSEARE FLIPPED!! weightings[0] = 1.; weightings[1] = 1.; weightings[2] = 1.; @@ -279,10 +352,12 @@ itk::Optimizer::ScalesType gTwoImageToOneImageMetric: } template -typename gTwoImageToOneImageMetric::ParametersType gTwoImageToOneImageMetric::GetParameters() const +typename gTwoImageToOneImageMetric:: +ParametersType gTwoImageToOneImageMetric::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()) { case X_ONLY: diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp index c1c6b8c..ecaf8ed 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp @@ -36,31 +36,9 @@ gfattori 08.11.2021 namespace itk { - -//FUNCTION DECLARATION NECESSARY FOR COPY/PASTE FROM vtkGDCMImageReader -// const char *gGetStringValueFromTag(const gdcm::Tag& t, const gdcm::DataSet& ds); - -const char *gGetStringValueFromTag(const gdcm::Tag& t, const gdcm::DataSet& ds) -{ - static std::string buffer; - buffer = ""; // cleanup previous call - if( ds.FindDataElement( t ) ) - { - const gdcm::DataElement& de = ds.GetDataElement( t ); - const gdcm::ByteValue *bv = de.GetByteValue(); - if( bv ) // Can be Type 2 - { - buffer = std::string( bv->GetPointer(), bv->GetLength() ); - // Will be padded with at least one \0 - } - } - // Since return is a const char* the very first \0 will be considered - return buffer.c_str(); -}; - - itkImageProcessor::itkImageProcessor() { + iNWUnits = 1; // std::cout<<"itkImageProcessor::itkImageProcessor contructor"<SetIdentity(); + m_InternalTransf1 = InternalTransformMetaInformation::New(); + m_InternalTransf2 = InternalTransformMetaInformation::New(); + + interpolator1 = InterpolatorType::New(); interpolator2 = InterpolatorType::New(); @@ -118,6 +100,9 @@ itkImageProcessor::itkImageProcessor() m_TransformMetaInfo = NULL; m_TransformMetaInfo = TransformMetaInformation::New(); + m_r23MetaInfo = NULL; + m_r23MetaInfo = R23MetaInformation::New(); + /* Initialise the projection geoemtry with defaults */ m_DRTGeometryMetaInfo = NULL; m_DRTGeometryMetaInfo @@ -133,7 +118,7 @@ itkImageProcessor::itkImageProcessor() Point3D[1]=0.; Point3D[2]=175.; m_DRTGeometryMetaInfo->SetProjectionCenter(Point3D); -// m_DRTGeometryMetaInfo->SetProjectionCenter(Point3D); + // m_DRTGeometryMetaInfo->SetProjectionCenter(Point3D); ImageType3D::SizeType ImageSize; ImageSize[0]=512; @@ -156,43 +141,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(); @@ -232,6 +180,37 @@ itkImageProcessor::~itkImageProcessor() } +void itkImageProcessor::SetNumberOfWorkingUnits(int iN){ + if(iN>0 && iN<50){ + iNWUnits = iN; + } +} + + +const CTVolumeImageMetaInformation::Pointer +itkImageProcessor::GetCTMetaInfo( + ){ + + return m_CTMetaInfo; + +} + +const RTGeometryMetaInformation::Pointer +itkImageProcessor::GetRTMetaInfo( + ){ + + return m_RTMetaInfo; + +} + +const DRTProjectionGeometryImageMetaInformation::Pointer +itkImageProcessor::GetDRTGeoMetaInfo(){ + + return m_DRTGeometryMetaInfo; + +} + + void itkImageProcessor::PrintSelf(std::ostream& os, Indent indent) const { Superclass::PrintSelf(os, indent); @@ -284,15 +263,27 @@ double itkImageProcessor::GetPanelOffset2(){ void itkImageProcessor::SetDegreeOfFreedom(tDegreeOfFreedomEnum dof) { - m_TransformMetaInfo->SetDegreeOfFreedom(dof); + m_r23MetaInfo->SetDegreeOfFreedom(dof); } tDegreeOfFreedomEnum itkImageProcessor::GetDegreeOfFreedom() { - return m_TransformMetaInfo->GetDegreeOfFreedom(); + return + m_r23MetaInfo->GetDegreeOfFreedom(); } +void itkImageProcessor::SetUseRotationsForImportOffset(bool bVal){ + + if(m_DRTGeometryMetaInfo == NULL){ + return; + } + m_DRTGeometryMetaInfo->SetUseRotationsForHiddenTransform(bVal); + + return; +} + + void itkImageProcessor::SetCustom_ImportTransform(double dTx, double dTy, double dTz, @@ -428,6 +419,7 @@ int itkImageProcessor::unload3DVolumeAndMeta(){ m_TransformMetaInfo = NULL; m_TransformMetaInfo = TransformMetaInformation::New(); + return 1; } @@ -471,7 +463,7 @@ int itkImageProcessor::load3DSerieFromFiles( std::vector 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(); @@ -516,11 +508,11 @@ int itkImageProcessor::load3DSerieFromFiles( std::vector v_fnames){ if (m_VolumeSourceDupli == NULL) std::cout << "NEVER HERE m_VolumeSourceDupli" << std::endl; -// if (m_VolumeSourceDupli) { -// std::cout << "NEVER HERE m_VolumeSourceDupli" << std::endl; -// m_VolumeSourceDupli = NULL; -// m_VolumeSourceDupli = DuplicatorType::New(); -// } + // if (m_VolumeSourceDupli) { + // std::cout << "NEVER HERE m_VolumeSourceDupli" << std::endl; + // m_VolumeSourceDupli = NULL; + // m_VolumeSourceDupli = DuplicatorType::New(); + // } m_VolumeSourceDupli->SetInputImage(caster3D->GetOutput()); m_VolumeSourceDupli->Update(); @@ -531,142 +523,7 @@ int itkImageProcessor::load3DSerieFromFiles( std::vector v_fnames){ } catch (itk::ExceptionObject & ex) { -// std::cout << ex << std::endl; - return EXIT_FAILURE; - } - - InternalImageType::Pointer m_InputImage = - m_VolumeSourceDupli->GetOutput(); - return - this->fill3DVolumeMeta(m_InputImage, m_PatOrientation); - -} - -int itkImageProcessor::load3DSerieFromFolder(const char * pcDirName) -{ - - using NamesGeneratorType = itk::GDCMSeriesFileNames; - NamesGeneratorType::Pointer nameGenerator = NamesGeneratorType::New(); - - nameGenerator->SetUseSeriesDetails(true); - nameGenerator->AddSeriesRestriction("0008|0021"); - nameGenerator->SetGlobalWarningDisplay(false); - nameGenerator->SetDirectory(pcDirName); - - - tPatOrientation m_PatOrientation - = tPatOrientation::NotDefined; - - try - { - using SeriesIdContainer = std::vector; - const SeriesIdContainer & seriesUID = nameGenerator->GetSeriesUIDs(); - auto seriesItr = seriesUID.begin(); - auto seriesEnd = seriesUID.end(); - - if (seriesItr != seriesEnd) - { - std::cout << "The directory: "; - std::cout << pcDirName << std::endl; - std::cout << "Contains the following DICOM Series: "; - std::cout << std::endl; - } - else - { - std::cout << "No DICOMs in: " << pcDirName << std::endl; - return EXIT_FAILURE; - } - - while (seriesItr != seriesEnd) - { - std::cout << seriesItr->c_str() << std::endl; - ++seriesItr; - } - - /* we expect only one serie at this point... */ - if(seriesUID.size() != 1){ - std::cout<<"More than one serie in the selected folder. returning..."<c_str(); - seriesItr++; - - std::cout << "\nReading: "; - std::cout << seriesIdentifier << std::endl; - using FileNamesContainer = std::vector; - FileNamesContainer fileNames = nameGenerator->GetFileNames(seriesIdentifier); - - using ImageIOType = itk::GDCMImageIO; - ImageIOType::Pointer dicomIO = ImageIOType::New(); - ImageSeriesReaderType3D::Pointer - imageSeriesReader3D = ImageSeriesReaderType3D::New(); - - imageSeriesReader3D->SetImageIO(dicomIO); - imageSeriesReader3D->SetFileNames(fileNames); - imageSeriesReader3D->SetNumberOfWorkUnits(12); - imageSeriesReader3D->ForceOrthogonalDirectionOff(); // properly read CTs with gantry tilt - - imageSeriesReader3D->Update(); - - if(fileNames.size()>0){ - - gdcm::Reader R; - R.SetFileName(fileNames.at(0).c_str()); - if(!R.Read()) - { - cerr<<"ERROR: cannot read file: "<SetInput(imageSeriesReader3D->GetOutput()); - caster3D->Update(); - - if (m_VolumeSourceDupli) { - m_VolumeSourceDupli = NULL; - m_VolumeSourceDupli = DuplicatorType::New(); - } - - m_VolumeSourceDupli->SetInputImage(caster3D->GetOutput()); - m_VolumeSourceDupli->Update(); - - caster3D = NULL; - - } - - } - catch (itk::ExceptionObject & ex) - { - std::cout << ex << std::endl; + // std::cout << ex << std::endl; return EXIT_FAILURE; } @@ -748,22 +605,22 @@ int itkImageProcessor::fill3DVolumeMeta( m_DRTGeometryMetaInfo->GetSCD2Offset() ); -// std::cout<< "///////////////// 3D VOLUME BEG ///////////////" <GetOriginLPS() <GetCOV() <GetOriginLPS() <GetCOV() <GetImportOffset() <GetLPS2IECDirections() <GetSize() <GetSpacing() <GetLPS2IECDirections() <GetSize() <GetSpacing() <GetProjectionAngleLPS() <GetSCD() <GetProjectionAngleLPS() <GetSCD() <GetProjectionAngleLPS() <GetSCD() <GetProjectionAngleLPS() <GetSCD() <UpdateProjectionGeometryMeta(); @@ -904,7 +761,7 @@ void itkImageProcessor::SetProjectionAngleOffsetIEC(double dOff1, double dOff2) m_DRTGeometryMetaInfo->SetProjectionAngle1OffsetIEC(dOff1); m_DRTGeometryMetaInfo->SetProjectionAngle2OffsetIEC(dOff2); - std::cout << "SetProjectionAngleOffsetIEC " << dOff1 << " " << dOff2 << std::endl; + std::cout << "SetProjectionAngleOffsetIEC " << dOff1 << " " << dOff2 << std::endl; } @@ -979,7 +836,7 @@ int itkImageProcessor::load2D(const char * pcFName){ strcpy( sTmpString,gGetStringValueFromTag( gdcm::Tag(0x0008,0x0008), ds)); std::string sLocalizerString = "LOCALIZER"; if (std::string(sTmpString).find(sLocalizerString) != std::string::npos) { -// std::cout << "Loading Localizer: "<GetOutput()->GetDirection(); -// std::cout<<"LocalizerImagDirectionDCM " <Update(); -// std::cout<< " ////////////// 2D Topo BEG " <GetOutput()->GetDirection()<GetOutput()->GetOrigin()<GetOutput()->GetSpacing()<GetOutput()->GetLargestPossibleRegion().GetSize()<GetOutput()->GetDirection()<GetOutput()->GetOrigin()<GetOutput()->GetSpacing()<GetOutput()->GetLargestPossibleRegion().GetSize()<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; - - // m_OutputTransform.Print(std::cout); - 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 - ) -{ - - //Convert all inputs into LPS - - ImageType3D::PointType m_TOffsetLPS = - m_IECtoLPSDirections * m_TranslationOffset; - - ImageType3D::PointType m_ROffsetLPS = - m_IECtoLPSDirections * m_RotationOffset; - - 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; - -} void itkImageProcessor::InitializeRegistration( @@ -1436,7 +1157,9 @@ void itkImageProcessor::InitializeRegistration( metric->SetSubtractMean(true); metric->SetMaxTranslation(maxTranslation); - m_TransformMetaInfo->SetDegreeOfFreedom(dof); + + // m_TransformMetaInfo->SetDegreeOfFreedom(dof); + m_r23MetaInfo->SetDegreeOfFreedom(dof); mimetric->ComputeGradientOff(); @@ -1461,6 +1184,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; @@ -1476,116 +1201,38 @@ void itkImageProcessor::InitializeRegistration( registration->SetFixedImageRegion1(region1); registration->SetFixedImageRegion2(region2); + mimetric->SetFixedImageRegion1(roiAutoReg1); + mimetric->SetFixedImageRegion2(roiAutoReg2); } - registration->SetTransformMetaInfo(m_TransformMetaInfo); + // 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) - { - //AAAAAA TODOOOOO CERCA - ImageType3D::PointType pFakeIsoLPS = - m_DRTImage1MetaInfo->GetProjectionOriginLPS(); + IsocTransf = TransformType::New(); + IsocTransf->SetComputeZYX(true); + IsocTransf->SetIdentity(); - CurrTransform = - CalculateInternalTransformV2( - ZeroPoint, - ZeroPoint, - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(), - m_DRTImage1MetaInfo->GetProjectionOriginLPS(), - IECtoLPS_Directions - ); + IsocTransf->SetRotation( + m_TransformMetaInfo->GetR()[0], + m_TransformMetaInfo->GetR()[1], + m_TransformMetaInfo->GetR()[2] + ); - } else { + TransformType::OutputVectorType TranslV; + TranslV[0] = m_TransformMetaInfo->GetT()[0]; + TranslV[1] = m_TransformMetaInfo->GetT()[1]; + TranslV[2] = m_TransformMetaInfo->GetT()[2]; - CurrTransform = CalculateInternalTransformV2( - ZeroPoint , - m_DRTGeometryMetaInfo->GetIECS2IECScannerR(), - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - 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) - { - ImageType3D::PointType pFakeIsoLPS = - m_DRTImage2MetaInfo->GetProjectionOriginLPS(); - - CurrTransform = - CalculateInternalTransformV2( - ZeroPoint, - ZeroPoint, - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(), - m_DRTImage2MetaInfo->GetProjectionOriginLPS(), - IECtoLPS_Directions - ); - - } else { - - - CurrTransform = CalculateInternalTransformV2( - ZeroPoint, - m_DRTGeometryMetaInfo->GetIECS2IECScannerR(), - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - 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->SetTranslation(TranslV); + registration->SetIsocIECTransform(IsocTransf); if (verbose) { registration->DebugOn(); @@ -1646,8 +1293,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; @@ -1673,7 +1320,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); @@ -1706,10 +1353,10 @@ int itkImageProcessor::StartRegistration(std::string extraInfo) } try { -// timer.Start("Registration"); + // timer.Start("Registration"); // Start the registration. registration->StartRegistration(); -// timer.Stop("Registration"); + // timer.Stop("Registration"); } catch (itk::ExceptionObject& err) { registration->ResetPipeline(); @@ -1740,41 +1387,111 @@ int itkImageProcessor::StartRegistration(std::string extraInfo) fs.close(); + + m_OptmizerValue = optimizer->GetValue(); + auto finalParameters = optimizer->GetCurrentPosition(); - auto finalParameters = transform1->GetParameters(); - optimizer->GetCurrentPosition(); - finalParameters = finalParameters - startParameters; + std::cout<<"REG FINAL PARS: "<SetParameters(finalParameters); - CalculateExternalUserTransform(offsetTransform, m_DRTImage1MetaInfo); + // auto finalParameters = transform1->GetParameters(); + // std::cout<<"-->startParameters "<finalParameters1 "<GetParameters()<finalParameters2 "<GetParameters()<GetCurrentPosition(); + // finalParameters = finalParameters - startParameters; - const int numberOfIterations = optimizer->GetCurrentIteration(); + // std::cout<<"startParameters "<GetCenter()<GetTranslation()<GetCenter()<GetTranslation()<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() { @@ -1784,6 +1501,7 @@ void itkImageProcessor::InitializeProjector() m_TransformMetaInfo == NULL ){ return; //TODO + return; } const double dtr = (atan(1.0) * 4.0) / 180.0; @@ -1793,41 +1511,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 = CalculateInternalTransformV2( - ZeroPoint , - ZeroPoint , - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(), - m_DRTImage1MetaInfo->GetProjectionOriginLPS(), - IECtoLPS_Directions - ); - - } else { - - - CurrTransform = CalculateInternalTransformV2( - ZeroPoint , - m_DRTGeometryMetaInfo->GetIECS2IECScannerR(), - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - 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); @@ -1858,33 +1554,12 @@ void itkImageProcessor::InitializeProjector() interpolator1->SetTransform(transform1); interpolator1->Initialize(); - if(m_RTMetaInfo == NULL) - { - - CurrTransform = - CalculateInternalTransformV2( - ZeroPoint, - ZeroPoint, - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(), - m_DRTImage2MetaInfo->GetProjectionOriginLPS(), - IECtoLPS_Directions - ); - - - } else { - - CurrTransform = CalculateInternalTransformV2( - ZeroPoint, - m_DRTGeometryMetaInfo->GetIECS2IECScannerR(), - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - 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); @@ -1921,7 +1596,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 @@ -1940,7 +1615,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 @@ -1970,7 +1645,7 @@ void itkImageProcessor::InitializeProjector() filter1->ChangeOriginOn(); filter1->UpdateOutputInformation(); -// std::cout<< "itkImageProcessor::InitializeProjector() " <Update(); @@ -2018,7 +1693,13 @@ int itkImageProcessor::unloadRTPlanAndMeta(){ std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << std::endl; m_RTMetaInfo = NULL; - return 1; + ImageType3D::PointType pZero; + pZero.Fill(0.); + + m_CTMetaInfo->SetImportOffset(pZero); + m_TransformMetaInfo->SetHiddenRotations(pZero); + + return 0; } @@ -2057,16 +1738,27 @@ void itkImageProcessor::loadRTPlanInfo( m_RTMetaInfo->SetIsocenterIECS(Point); + + ImageType3D::PointType + pZero (3); + pZero.Fill(0.); + m_CTMetaInfo->SetImportOffset( CalcImportVolumeOffset( m_RTMetaInfo->GetIsocenterIECS(), m_CTMetaInfo->GetLPS2IECDirections(), m_RTMetaInfo->GetIsocenterLPS(), m_DRTGeometryMetaInfo->GetIECS2IECScannerT(), - m_DRTGeometryMetaInfo->GetIECS2IECScannerR() + (m_DRTGeometryMetaInfo->GetUseRotationsForHiddenTransform() ? + m_DRTGeometryMetaInfo->GetIECS2IECScannerR() : pZero) ) ); + m_TransformMetaInfo->SetHiddenRotations( + (m_DRTGeometryMetaInfo->GetUseRotationsForHiddenTransform() ? + m_DRTGeometryMetaInfo->GetIECS2IECScannerR() : pZero) + ); + this->UpdateProjectionGeometryMeta(); this->InitializeProjector(); @@ -2118,11 +1810,11 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ IsocenterOffsetLPS = m_CTMetaInfo->ConvertIECPointToLPSPoint( m_DRTGeometryMetaInfo->GetProjectionCenterOffset1()); -// std::cout<< "///////////////// PGEOM META BEG ///////////////" <GetProjectionCenter() <GetProjectionCenter() <SetProjectionAngleLPS( this->CalcProjectionAngleLPS( @@ -2185,7 +1877,6 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ this->CalcProjectionAngleLPS( m_CTMetaInfo->GetPatientOrientation(), m_DRTGeometryMetaInfo->GetProjectionAngle1IEC() ), - //m_DRTGeometryMetaInfo->GetProjectionAngle1IEC(),// m_DRTImage1MetaInfo->GetProjectionAngleLPS(), Std_DRT2LPS ) @@ -2196,10 +1887,32 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ this->CalcProjectionAngleLPS( m_CTMetaInfo->GetPatientOrientation(), m_DRTGeometryMetaInfo->GetProjectionAngle1IEC()), - //m_DRTGeometryMetaInfo->GetProjectionAngle1IEC(),//m_DRTImage1MetaInfo->GetProjectionAngleLPS(), 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 @@ -2264,7 +1977,6 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ this->CalcProjectionAngleLPS( m_CTMetaInfo->GetPatientOrientation(), m_DRTGeometryMetaInfo->GetProjectionAngle2IEC()) , - //m_DRTGeometryMetaInfo->GetProjectionAngle2IEC(),//m_DRTImage2MetaInfo->GetProjectionAngleLPS(), Std_DRT2LPS )); @@ -2277,7 +1989,6 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ this->CalcProjectionAngleLPS( m_CTMetaInfo->GetPatientOrientation(), m_DRTGeometryMetaInfo->GetProjectionAngle2IEC()) , - //m_DRTGeometryMetaInfo->GetProjectionAngle2IEC(),//m_DRTImage2MetaInfo->GetProjectionAngleLPS(), Std_DRT2LPS ) @@ -2288,10 +1999,32 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ this->CalcProjectionAngleLPS( m_CTMetaInfo->GetPatientOrientation(), m_DRTGeometryMetaInfo->GetProjectionAngle2IEC()) , - //m_DRTGeometryMetaInfo->GetProjectionAngle2IEC(),//m_DRTImage2MetaInfo->GetProjectionAngleLPS(), 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); + + } + } @@ -2324,16 +2057,17 @@ itkImageProcessor::CalcDRTImageOrigin( void itkImageProcessor::GetProjectionImages(){ -// std::cout<< "itkImageProcessor::GetProjectionImages" <SetComputeZYX(true); transform1->SetIdentity(); - InternalImageType::DirectionType IECtoLPS_Directions; - - IECtoLPS_Directions = - m_CTMetaInfo->GetLPS2IECDirections().GetTranspose(); TransformType::Pointer CurrTransform; - - - if(m_RTMetaInfo == NULL) - { - CurrTransform = - CalculateInternalTransformV2( - ZeroPoint, - ZeroPoint, - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(), - m_DRTImage1MetaInfo->GetProjectionOriginLPS(), - IECtoLPS_Directions - ); - - } else { - - CurrTransform = CalculateInternalTransformV2( - ZeroPoint , - m_DRTGeometryMetaInfo->GetIECS2IECScannerR(), - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - 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); @@ -2390,6 +2098,16 @@ void itkImageProcessor::GetProjectionImages(){ transform1 ->SetCenter( m_DRTImage1MetaInfo->GetProjectionOriginLPSZero()); +// std::cout<< "CurrTransform Rotations: "<< +// CurrTransform->GetAngleX() <<" "<< +// CurrTransform->GetAngleY() <<" "<< +// CurrTransform->GetAngleZ() <<" "<< std::endl; + +// std::cout<< "CurrTransform Translations: "<< +// CurrTransform->GetTranslation()<GetCenter()<Print(std::cout); // The parameters of interpolator1, such as ProjectionAngle and FocalPointToIsocenterDistance @@ -2401,37 +2119,12 @@ void itkImageProcessor::GetProjectionImages(){ transform2->SetComputeZYX(true); transform2->SetIdentity(); - - - if(m_RTMetaInfo == NULL) - { - - CurrTransform = - CalculateInternalTransformV2( - ZeroPoint, - ZeroPoint, - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(), - m_DRTImage2MetaInfo->GetProjectionOriginLPS(), - IECtoLPS_Directions - ); - - } else { - - CurrTransform = CalculateInternalTransformV2( - ZeroPoint, - m_DRTGeometryMetaInfo->GetIECS2IECScannerR(), - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - 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); @@ -2616,22 +2309,22 @@ vtkImageData* itkImageProcessor::GetLocalizer1VTK() ImageType3D::PointType LastVoxelPosition = origin3D + (imagDirection * Dest); -// std::cout<<" -------- Localizer 1 ITK --------"<GetOutput()->GetBounds(); -// std::cout<< "-------- Localizer 1 VTK --------" <SetInput( imageDRT1In ); rescaler1->Update(); -// using ImageCalculatorFilterType = itk::MinimumMaximumImageCalculator; -// auto imageCalculatorFilter = ImageCalculatorFilterType::New(); -// imageCalculatorFilter->SetImage(imageDRT1In); -// imageCalculatorFilter->Compute(); + // using ImageCalculatorFilterType = itk::MinimumMaximumImageCalculator; + // auto imageCalculatorFilter = ImageCalculatorFilterType::New(); + // imageCalculatorFilter->SetImage(imageDRT1In); + // imageCalculatorFilter->Compute(); -// std::cout<< "itkImageProcessor::imageDRT1In() " << -// imageCalculatorFilter <; -// auto imageCalculatorFilter2 = ImageCalculatorFilterType2::New(); -// imageCalculatorFilter2->SetImage(rescaler1->GetOutput()); -// imageCalculatorFilter2->Compute(); -// std::cout<< "itkImageProcessor::imageDRT2In() " << -// imageCalculatorFilter2 <; + // auto imageCalculatorFilter2 = ImageCalculatorFilterType2::New(); + // imageCalculatorFilter2->SetImage(rescaler1->GetOutput()); + // imageCalculatorFilter2->Compute(); + // std::cout<< "itkImageProcessor::imageDRT2In() " << + // imageCalculatorFilter2 <SetInput(rescaler1->GetOutput()); toVTK2D1->Update(); @@ -2890,12 +2583,12 @@ vtkImageData* itkImageProcessor::GetProjection2VTK() - //double* dBounds = toVTK2D2->GetOutput()->GetBounds(); -// std::cout<< "-------- Proj 2 --------" <GetOutput()->GetBounds(); + // std::cout<< "-------- Proj 2 --------" <GetOutput(); @@ -2929,7 +2622,7 @@ void itkImageProcessor::WriteProjectionImages() try { -// std::cout << "Writing image 1 " << std::endl; + // std::cout << "Writing image 1 " << std::endl; writer1->Update(); } catch (itk::ExceptionObject & err) @@ -2943,7 +2636,7 @@ void itkImageProcessor::WriteProjectionImages() try { -// std::cout << "Writing image 2" << std::endl; + // std::cout << "Writing image 2" << std::endl; writer2->Update(); } catch (itk::ExceptionObject & err) @@ -2965,13 +2658,12 @@ void itkImageProcessor::SetInitialTranslations(double dX, double dY, double dZ) //todo } - ImageType3D::PointType Translations; Translations[0]= dX; Translations[1]= dY; Translations[2]= dZ; - m_TransformMetaInfo->SetTranslations(Translations); + m_TransformMetaInfo->SetUserTranslations(Translations); } void itkImageProcessor::SetInitialRotations(double dX, double dY, double dZ) @@ -2981,49 +2673,96 @@ void itkImageProcessor::SetInitialRotations(double dX, double dY, double dZ) Rotations[1]= dY; Rotations[2]= dZ; - m_TransformMetaInfo->SetRotations(Rotations); + m_TransformMetaInfo->SetUserRotations(Rotations); } -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->GetTranslations(); - Rotations = m_TransformMetaInfo->GetRotations(); + 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; } -// Doubts about this user thing -void itkImageProcessor::GetFinalTranslations(double& dX, double& dY, double& dZ) -{ - ImageType3D::PointType translations = m_TransformMetaInfo->GetTranslations(); - ImageType3D::PointType userTranslations = m_TransformMetaInfo->GetUserTranslations(); - dX = translations[0] + userTranslations[0]; - dY = translations[1] + userTranslations[1]; - dZ = translations[2] + userTranslations[2]; -} -void itkImageProcessor::GetFinalRotations(double& dRX, double& dRY, double& dRZ) -{ - ImageType3D::PointType rotations = m_TransformMetaInfo->GetRotations(); - ImageType3D::PointType userRotations = m_TransformMetaInfo->GetUserRotations(); +TransformType::Pointer +itkImageProcessor::GetCompleteIsocentricTransform +(/*ImageType3D::PointType TransformCenter*/){ - dRX = rotations[0] + userRotations[0]; - dRY = rotations[1] + userRotations[1]; - dRZ = rotations[2] + userRotations[2]; + TransformType::Pointer + regTransform = TransformType::New(); + regTransform->SetComputeZYX(true); + regTransform->SetIdentity(); + + if(m_TransformMetaInfo == nullptr || + m_CTMetaInfo == nullptr){ + return nullptr; + } + +// std::cout<<" --- GetCompleteIsocentricTransform --- "<GetR()[0] <<" "<< +// m_TransformMetaInfo->GetR()[1] <<" "<< +// m_TransformMetaInfo->GetR()[2] <<" "<GetT()[0] <<" "<< +// m_TransformMetaInfo->GetT()[1] <<" "<< +// m_TransformMetaInfo->GetT()[2] <<" "<GetImportOffset()[0]<<" "<< +// m_CTMetaInfo->GetImportOffset()[1]<<" "<< +// m_CTMetaInfo->GetImportOffset()[2]<<" "<GetLPS2IECDirections() * m_RTMetaInfo->GetIsocenterLPS(); + + mISORTIEC[0] *= -1; + mISORTIEC[1] *= -1; + mISORTIEC[2] *= -1; + + TransformType::Pointer mTransf= + MapTransformToNewOrigin( + mISORTIEC, + m_TransformMetaInfo->GetT(), + m_TransformMetaInfo->GetR() + ); + + itk::TransformMetaInformation::PointType + pTM = m_CTMetaInfo->GetLPS2IECDirections() * m_CTMetaInfo->GetImportOffset(); + + TransformType::OutputVectorType pTransl; + + pTransl[0] = mTransf->GetTranslation()[0] - pTM[0]; + pTransl[1] = mTransf->GetTranslation()[1] - pTM[1]; + pTransl[2] = mTransf->GetTranslation()[2] - pTM[2]; + + mTransf->SetTranslation(pTransl); + itk::PointpZero; + pZero.Fill(0.); + mTransf->SetCenter(pZero); + + return mTransf; } @@ -3049,10 +2788,6 @@ double itkImageProcessor::GetCurrentMetricValue() } } -void itkImageProcessor::SetROI(double, double, double, double) -{ - //TODO: Not implemanted yet. ROI are hard coded and can be enabled using SetFullROI. -} void itkImageProcessor::SetOptimizer(std::string optimizer) { @@ -3111,6 +2846,7 @@ void itkImageProcessor::SetRegionFixed1( std::cout << "itkImageProcessor " << std::endl; std::cout << roiAutoReg1 << std::endl; + } void itkImageProcessor::SetRegionFixed2( diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.h b/reg23Topograms/itkDTRrecon/itkImageProcessor.h index 553212e..df9f45b 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.h +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.h @@ -51,6 +51,8 @@ gfattori 08.11.2021 #include "itkQtIterationUpdate.h" +#include "itkImageProcessorHelpers.h" + namespace itk { @@ -76,15 +78,16 @@ public: CommandIterationUpdate::Pointer optimizerObserver; + /** Set number of logic CPU to be made available to interpolators*/ + void SetNumberOfWorkingUnits(int iN); /** Input data load methods*/ - int load3DSerieFromFolder(const char* ); int load3DSerieFromFiles( std::vector ); int unload3DVolumeAndMeta(); int load2D(const char *); int unload2DAndMeta(int); - int query2DimageReconstructionDiameter(const char*); + //int query2DimageReconstructionDiameter(const char*); void loadRTPlanInfo(double, double, double, double, double ,double); int unloadRTPlanAndMeta(); @@ -118,7 +121,6 @@ public: tDegreeOfFreedomEnum GetDegreeOfFreedom(); void SetCustom_ProjectionCenterOffsetFixedAxes_IEC(double,double,double,double,double,double); - void SetCustom_ProjectionCenterFixedAxes_IEC(double,double,double); /** Intensity threshold used for image projection, @@ -128,8 +130,16 @@ public: /** Custom settings of the projection images */ void SetCustom_2Dres(double,double,double,double); void SetCustom_2Dsize(int,int,int,int); + + /** Custom transform to map IEC of imaging device into IEC-Support + * Tx Ty Tz [mm] Rx Ry Rz [deg] + */ void SetCustom_ImportTransform(double, double, double, double, double, double); + + /** Should rotations be used when computing import offset and hidden transform? */ + void SetUseRotationsForImportOffset(bool bVal); + void SetCustom_UpdateMetaInfo(); @@ -137,21 +147,42 @@ 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); + /** Returns user components only of the autoReg + * */ + Optimizer::ParametersType + GetFinalR23Parameters(); /** Initialize the registration pipeline*/ void InitializeRegistration(double, double, eDegreeOfFreedomType); /** Start the registration process*/ int StartRegistration(std::string extraInfo); - /** Get transform parameters for 3D Volume */ - double* GetTransformParameters(); + ///** Get transform parameters for 3D Volume */ + /** 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. Import offset and hidden not included. + * This can be use to update UI transform display + * ParametersType :: Rx Ry Rz Tx Ty Tz + */ + Optimizer::ParametersType + GetUserIsocentricTransform(); + + /** Return useful meta info for external use + * Consumer should check for nullptr... + */ + const CTVolumeImageMetaInformation::Pointer + GetCTMetaInfo(); + const RTGeometryMetaInformation::Pointer + GetRTMetaInfo(); + const DRTProjectionGeometryImageMetaInformation::Pointer + GetDRTGeoMetaInfo(); /** Initialize projection geometry */ void InitializeProjector(); @@ -184,24 +215,26 @@ public: void WriteProjectionImages(); void Write2DImages(); + /** Auto Reg23 methods */ + /** Get the current cost function value from the optimizer*/ double GetOptimizerValue(); /** Get the cost function value for the current transform*/ double GetCurrentMetricValue(); - void SetROI(double, double, double, double); - + /** Set Optimizer type */ void SetOptimizer(std::string); + /** Set Metric type */ void SetMetric(std::string); + /** Use full image as ROI */ void SetFullROI(bool); + /** Maximum number of iterations for the optimizer */ void SetMaxNumberOfIterations(int); + /** Define the region to be used as ROI on fixed images */ void SetRegionFixed1(int,int,int,int); void SetRegionFixed2(int,int,int,int); - /** Optimizer which tries to find the minimn (Powell Optimizer)*/ - using OptimizerType = itk::PowellOptimizer; - protected: /** Various pixel types */ using InternalPixelType = float; @@ -232,10 +265,11 @@ 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; using InterpolatorType = itk::gSiddonJacobsRayCastInterpolateImageFunction; using ResampleFilterType = itk::ResampleImageFilter; + /** Optimizer which tries to find the minimun (Powell Optimizer)*/ + using OptimizerType = itk::PowellOptimizer; /** Optimizer which tries to find the minimn (Powell Optimizer)*/ using AmoebaOptimizerType = itk::AmoebaOptimizer; /** Optimizer which samples the whol space */ @@ -264,41 +298,10 @@ private: /** ITK to VTK filter */ using ITKtoVTKFilterType = itk::ImageToVTKImageFilter; - void - CalculateExternalUserTransform( - TransformType::Pointer transform, - DRTImageMetaInformation::Pointer imageMetaInfo); - -// TransformType::Pointer -// CalculateInternalTransform( -// ImageType3D::PointType m_TranslationOffset, -// ImageType3D::PointType m_RotationOffset, -// ImageType3D::PointType m_TranslationUser, -// ImageType3D::PointType m_RotationUser, -// ImageType3D::PointType m_ProjectionTransformCenter, -// ImageType3D::PointType m_UserTransformCenter, -// ImageType3D::PointType m_OffsetTransformCenter, -// InternalImageType::DirectionType m_IECtoLPSDirections -// ); - - /* 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 - ); - - - TransformType::Pointer transform1, - transform2; + transform2, + IsocTransf; InterpolatorType::Pointer interpolator1, @@ -360,20 +363,6 @@ private: ); -// ImageType3D::PointType -// CalcDRTImageOffset( -// ImageType3D::PointType m_DRTOffset, -// double dAngle, -// InternalImageType::DirectionType stdDRT2LPS -// ); - - TransformType::Pointer - MapTransformToNewOrigin( - ImageType3D::PointType m_COR, - ImageType3D::PointType m_Translations, - ImageType3D::PointType m_Rotations - ); - double CalcProjectionAngleLPS( tPatOrientation pOrient, @@ -407,9 +396,6 @@ private: * m_Projection1VTKTransform, * m_Projection2VTKTransform; - /*Transformation Parameters */ - double dTransfParam[6]; - /** * Many meta containers */ @@ -433,6 +419,12 @@ private: TransformMetaInformation::Pointer m_TransformMetaInfo; + R23MetaInformation::Pointer + m_r23MetaInfo; + + InternalTransformMetaInformation::Pointer + m_InternalTransf1, + m_InternalTransf2; double m_OptmizerValue; int m_MaxNumberOfIterations; @@ -448,6 +440,9 @@ private: bool m_UseMutualInformation; bool m_UseDumptoFile; + int + iNWUnits; + }; diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp new file mode 100644 index 0000000..bb93324 --- /dev/null +++ b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp @@ -0,0 +1,165 @@ +#include "itkImageProcessorHelpers.h" +#include +#include + +namespace itk { + + +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 + ){ + + 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; +} + + + +TransformType::Pointer +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; + +} + +const char* queryStationName(const char* pcFName){ + + static std::string buffer; + + /* Check if we can open the file */ + gdcm::Reader R; + R.SetFileName(pcFName); + if(!R.Read()) + {std::cerr<<"ERROR: cannot read file: "<GetPointer(), bv->GetLength() ); + // Will be padded with at least one \0 + } + } + // Since return is a const char* the very first \0 will be considered + return buffer.c_str(); +} + +} diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h new file mode 100644 index 0000000..ea16394 --- /dev/null +++ b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h @@ -0,0 +1,55 @@ + +#ifndef ITKIMAGEPROCESSORHELPERS_H +#define ITKIMAGEPROCESSORHELPERS_H + +#include "itkEuler3DTransform.h" +#include "itkImage.h" +#include "itkEuler3DTransform.h" + +#include "itkCommand.h" +#include "itkObject.h" +#include "itkObjectFactory.h" + +#include +#include + +namespace itk +{ + + +using TransformType = itk::Euler3DTransform; +constexpr static unsigned int Dimension = 3; +using PixelType3D = short; +using InternalPixelType = float; + +using ImageType3D = itk::Image; +using InternalImageType = itk::Image; + + +TransformType::Pointer +MapTransformToNewOrigin( + ImageType3D::PointType m_COR, + ImageType3D::PointType m_Translations, + ImageType3D::PointType m_Rotations + ); + + +TransformType::Pointer + CalculateInternalTransformV3( + ImageType3D::PointType m_Translation, + ImageType3D::PointType m_Rotation, + ImageType3D::PointType m_CalibratedProjectionCenter, + ImageType3D::PointType m_RTIsocenter, + InternalImageType::DirectionType m_IECtoLPSDirections + ); + +const char *gGetStringValueFromTag(const gdcm::Tag& t, const gdcm::DataSet& ds); + +int query2DimageReconstructionDiameter(const char*); + +const char* queryStationName(const char*); + +} + + +#endif diff --git a/reg23Topograms/itkDTRrecon/itkQtIterationUpdate.h b/reg23Topograms/itkDTRrecon/itkQtIterationUpdate.h index 91e982a..89b231d 100644 --- a/reg23Topograms/itkDTRrecon/itkQtIterationUpdate.h +++ b/reg23Topograms/itkDTRrecon/itkQtIterationUpdate.h @@ -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::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; }