382 lines
13 KiB
C++
382 lines
13 KiB
C++
/*=========================================================================
|
|
*
|
|
* Copyright NumFOCUS
|
|
*
|
|
* Licensed under the Apache License, Version 2.0 (the "License");
|
|
* you may not use this file except in compliance with the License.
|
|
* You may obtain a copy of the License at
|
|
*
|
|
* http://www.apache.org/licenses/LICENSE-2.0.txt
|
|
*
|
|
* Unless required by applicable law or agreed to in writing, software
|
|
* distributed under the License is distributed on an "AS IS" BASIS,
|
|
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
|
* See the License for the specific language governing permissions and
|
|
* limitations under the License.
|
|
*
|
|
*=========================================================================*/
|
|
#ifndef itkTwoProjectionImageRegistrationMethod_hxx
|
|
#define itkTwoProjectionImageRegistrationMethod_hxx
|
|
|
|
#include "itkCommand.h"
|
|
|
|
#include "itkTwoProjectionImageRegistrationMethod.h"
|
|
|
|
namespace itk {
|
|
|
|
// TwoProjectionImageRegistrationMethod from Jian Wu.
|
|
// This part is not really needed and functionality should be moved to
|
|
// to the itkImageProcessor.cpp. This is just abstraction with an additional step.
|
|
|
|
template <typename TFixedImage, typename TMovingImage>
|
|
TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::TwoProjectionImageRegistrationMethod()
|
|
{
|
|
this->SetNumberOfRequiredOutputs(1); // for the Transform
|
|
|
|
m_FixedImage1 = nullptr; // has to be provided by the user.
|
|
m_FixedImage2 = nullptr; // has to be provided by the user.
|
|
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.
|
|
m_Optimizer = nullptr; // has to be provided by the user.
|
|
|
|
m_InitialOptimizerParameters = ParametersType(1);
|
|
m_LastOptimizerParameters = ParametersType(1);
|
|
|
|
m_InitialOptimizerParameters.Fill(0.0f);
|
|
m_LastOptimizerParameters.Fill(0.0f);
|
|
|
|
m_FixedImageRegionDefined1 = false;
|
|
m_FixedImageRegionDefined2 = false;
|
|
|
|
TransformOutputPointer transformDecorator = static_cast<TransformOutputType*>(this->MakeOutput(0).GetPointer());
|
|
|
|
this->ProcessObject::SetNthOutput(0, transformDecorator.GetPointer());
|
|
}
|
|
|
|
/*
|
|
* Set the region of the fixed image 1 to be considered for registration
|
|
*/
|
|
template <typename TFixedImage, typename TMovingImage>
|
|
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::SetFixedImageRegion1(
|
|
const FixedImageRegionType& region1)
|
|
{
|
|
m_FixedImageRegion1 = region1;
|
|
m_FixedImageRegionDefined1 = true;
|
|
}
|
|
|
|
/*
|
|
* Set the region of the fixed image 2 to be considered for registration
|
|
*/
|
|
template <typename TFixedImage, typename TMovingImage>
|
|
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::SetFixedImageRegion2(
|
|
const FixedImageRegionType& region2)
|
|
{
|
|
m_FixedImageRegion2 = region2;
|
|
m_FixedImageRegionDefined2 = true;
|
|
}
|
|
|
|
|
|
/*
|
|
* Initialize by setting the interconnects between components.
|
|
*/
|
|
template <typename TFixedImage, typename TMovingImage>
|
|
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::Initialize()
|
|
{
|
|
|
|
if (!m_FixedImage1) {
|
|
itkExceptionMacro(<< "FixedImage1 is not present");
|
|
}
|
|
|
|
if (!m_FixedImage2) {
|
|
itkExceptionMacro(<< "FixedImage2 is not present");
|
|
}
|
|
|
|
if (!m_MovingImage) {
|
|
itkExceptionMacro(<< "MovingImage is not present");
|
|
}
|
|
|
|
if (!m_Metric) {
|
|
itkExceptionMacro(<< "Metric is not present");
|
|
}
|
|
|
|
if (!m_Optimizer) {
|
|
itkExceptionMacro(<< "Optimizer is not present");
|
|
}
|
|
|
|
if (!m_Transform1) {
|
|
itkExceptionMacro(<< "Transform1 is not present");
|
|
}
|
|
|
|
if (!m_Transform2) {
|
|
itkExceptionMacro(<< "Transform2 is not present");
|
|
}
|
|
|
|
//
|
|
// Connect the transform to the Decorator.
|
|
auto* transformOutput = static_cast<TransformOutputType*>(this->ProcessObject::GetOutput(0));
|
|
|
|
//transformOutput->Set(m_Transform1.GetPointer());
|
|
transformOutput->Set(m_IsocIECTransform.GetPointer());
|
|
|
|
if (!m_Interpolator1) {
|
|
itkExceptionMacro(<< "Interpolator1 is not present");
|
|
}
|
|
|
|
if (!m_Interpolator2) {
|
|
itkExceptionMacro(<< "Interpolator2 is not present");
|
|
}
|
|
|
|
if (!m_Filter1) {
|
|
itkExceptionMacro(<< "Filter1 is not present");
|
|
}
|
|
|
|
if (!m_Filter1) {
|
|
itkExceptionMacro(<< "Filter2 is not present");
|
|
}
|
|
|
|
if (!m_TransformMetaInfo) {
|
|
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);
|
|
|
|
m_Metric->SetInterpolator1(m_Interpolator1);
|
|
m_Metric->SetInterpolator2(m_Interpolator2);
|
|
|
|
if (m_FixedImageRegionDefined1) {
|
|
m_Metric->SetFixedImageRegion1(m_FixedImageRegion1);
|
|
std::cout << "m_FixedImageRegionDefined1 is defined " << std::endl;
|
|
} else {
|
|
m_Metric->SetFixedImageRegion1(m_FixedImage1->GetBufferedRegion());
|
|
}
|
|
|
|
if (m_FixedImageRegionDefined2) {
|
|
m_Metric->SetFixedImageRegion2(m_FixedImageRegion2);
|
|
std::cout << "m_FixedImageRegionDefined2 is defined " << std::endl;
|
|
} else {
|
|
m_Metric->SetFixedImageRegion2(m_FixedImage2->GetBufferedRegion());
|
|
}
|
|
|
|
m_Metric->Initialize();
|
|
|
|
// Setup the optimizer
|
|
m_Optimizer->SetScales(m_Metric->GetWeightings());
|
|
m_Optimizer->SetCostFunction(m_Metric);
|
|
|
|
m_InitialOptimizerParameters =
|
|
m_Metric->GetParameters();
|
|
|
|
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);
|
|
|
|
// 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 (m_CancelRequest)
|
|
//// filter->SetAbortGenerateData(true); // may be handled by filter
|
|
|
|
// std::cout << "OnInternalFilterProgressReceptor :: filter (TRUE) " << std::endl;
|
|
// }
|
|
//}
|
|
|
|
/*
|
|
* Starts the Registration Process
|
|
*/
|
|
template <typename TFixedImage, typename TMovingImage>
|
|
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::StartRegistration()
|
|
{
|
|
|
|
itk::ProgressReporter progress(
|
|
this,
|
|
1, 1000,100);
|
|
|
|
ParametersType empty(1);
|
|
empty.Fill(0.0);
|
|
try {
|
|
// initialize the interconnects between components
|
|
this->Initialize();
|
|
if (GetDebug()) {
|
|
m_Metric->Print(std::cout);
|
|
m_Optimizer->Print(std::cout);
|
|
this->Print(std::cout);
|
|
}
|
|
} catch (ExceptionObject& err) {
|
|
m_LastOptimizerParameters = empty;
|
|
|
|
// pass exception to caller
|
|
throw err;
|
|
}
|
|
|
|
this->StartOptimization();
|
|
}
|
|
|
|
/*
|
|
* Starts the Optimization process
|
|
*/
|
|
template <typename TFixedImage, typename TMovingImage>
|
|
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::StartOptimization()
|
|
{
|
|
try {
|
|
// do the optimization
|
|
m_Optimizer->StartOptimization();
|
|
|
|
} catch (ExceptionObject& err) {
|
|
// An error has occurred in the optimization.
|
|
// Update the parameters
|
|
m_LastOptimizerParameters = m_Optimizer->GetCurrentPosition();
|
|
|
|
// Pass exception to caller
|
|
throw err;
|
|
}
|
|
|
|
// get the results
|
|
m_LastOptimizerParameters = m_Optimizer->GetCurrentPosition();
|
|
|
|
//Update transform1 and transform2 with the last transform parameters from the optimizer.
|
|
m_Metric->SetTransformParameters(m_LastOptimizerParameters);
|
|
}
|
|
|
|
template <typename TFixedImage, typename TMovingImage>
|
|
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::PrintSelf(std::ostream& os, Indent indent) const
|
|
{
|
|
Superclass::PrintSelf(os, indent);
|
|
os << indent << "Metric: " << m_Metric.GetPointer() << std::endl;
|
|
os << indent << "Optimizer: " << m_Optimizer.GetPointer() << std::endl;
|
|
os << indent << "Transform1: " << m_Transform1.GetPointer() << std::endl;
|
|
os << indent << "Transform2: " << m_Transform2.GetPointer() << std::endl;
|
|
os << indent << "Interpolator 1: " << m_Interpolator1.GetPointer() << std::endl;
|
|
os << indent << "Interpolator 2: " << m_Interpolator2.GetPointer() << std::endl;
|
|
os << indent << "Fixed Image 1: " << m_FixedImage1.GetPointer() << std::endl;
|
|
os << indent << "Fixed Image 2: " << m_FixedImage2.GetPointer() << std::endl;
|
|
os << indent << "Moving Image: " << m_MovingImage.GetPointer() << std::endl;
|
|
os << indent << "Transform MetaInfo: " << m_TransformMetaInfo.GetPointer() << std::endl;
|
|
os << indent << "Fixed Image 1 Region Defined: " << m_FixedImageRegionDefined1 << std::endl;
|
|
os << indent << "Fixed Image 2 Region Defined: " << m_FixedImageRegionDefined2 << std::endl;
|
|
os << indent << "Fixed Image 1 Region: " << m_FixedImageRegion1 << std::endl;
|
|
os << indent << "Fixed Image 2 Region: " << m_FixedImageRegion2 << std::endl;
|
|
os << indent << "Initial Optimizer Parameters: " << m_InitialOptimizerParameters << std::endl;
|
|
os << indent << "Last Optimizer Parameters: " << m_LastOptimizerParameters << std::endl;
|
|
}
|
|
|
|
template <typename TFixedImage, typename TMovingImage>
|
|
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::GenerateData()
|
|
{
|
|
|
|
this->StartRegistration();
|
|
}
|
|
|
|
template <typename TFixedImage, typename TMovingImage>
|
|
const typename TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::TransformOutputType*
|
|
TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::GetOutput() const
|
|
{
|
|
return static_cast<const TransformOutputType*>(this->ProcessObject::GetOutput(0));
|
|
}
|
|
|
|
template <typename TFixedImage, typename TMovingImage>
|
|
DataObject::Pointer
|
|
TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::MakeOutput(DataObjectPointerArraySizeType output)
|
|
{
|
|
switch (output) {
|
|
case 0:
|
|
return static_cast<DataObject*>(TransformOutputType::New().GetPointer());
|
|
break;
|
|
default:
|
|
itkExceptionMacro("MakeOutput request for an output number larger than the expected number of outputs");
|
|
return nullptr;
|
|
}
|
|
}
|
|
|
|
template <typename TFixedImage, typename TMovingImage>
|
|
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::SetFixedImage1(const FixedImageType* fixedImage1)
|
|
{
|
|
itkDebugMacro("setting Fixed Image 1 to " << fixedImage1);
|
|
|
|
if (this->m_FixedImage1.GetPointer() != fixedImage1) {
|
|
this->m_FixedImage1 = fixedImage1;
|
|
|
|
// Process object is not const-correct so the const_cast is required here
|
|
this->ProcessObject::SetNthInput(0, const_cast<FixedImageType*>(fixedImage1));
|
|
|
|
this->Modified();
|
|
}
|
|
}
|
|
|
|
template <typename TFixedImage, typename TMovingImage>
|
|
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::SetFixedImage2(const FixedImageType* fixedImage2)
|
|
{
|
|
itkDebugMacro("setting Fixed Image 2 to " << fixedImage2);
|
|
|
|
if (this->m_FixedImage2.GetPointer() != fixedImage2) {
|
|
this->m_FixedImage2 = fixedImage2;
|
|
|
|
// Process object is not const-correct so the const_cast is required here
|
|
this->ProcessObject::SetNthInput(0, const_cast<FixedImageType*>(fixedImage2));
|
|
|
|
this->Modified();
|
|
}
|
|
}
|
|
|
|
template <typename TFixedImage, typename TMovingImage>
|
|
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::SetMovingImage(const MovingImageType* movingImage)
|
|
{
|
|
itkDebugMacro("setting Moving Image to " << movingImage);
|
|
|
|
if (this->m_MovingImage.GetPointer() != movingImage) {
|
|
this->m_MovingImage = movingImage;
|
|
|
|
// Process object is not const-correct so the const_cast is required here
|
|
this->ProcessObject::SetNthInput(1, const_cast<MovingImageType*>(movingImage));
|
|
|
|
this->Modified();
|
|
}
|
|
}
|
|
|
|
} // end namespace itk
|
|
|
|
#endif
|