Adding new files for Automatic Registration
This commit is contained in:
Proton local user
32
reg23Topograms/itkDTRrecon/autoreg/CMakeLists.txt
Normal file
32
reg23Topograms/itkDTRrecon/autoreg/CMakeLists.txt
Normal file
@ -0,0 +1,32 @@
|
||||
SET(LIB_NAME autoreg)
|
||||
SET(CMAKE_CXX_FLAGS "-std=c++11 -fPIC")
|
||||
|
||||
|
||||
find_package(ITK REQUIRED)
|
||||
include(${ITK_USE_FILE})
|
||||
INCLUDE_DIRECTORIES(${ITK_INCLUDE_DIRS})
|
||||
|
||||
|
||||
SET(HDR
|
||||
itkDRTHelpers.h
|
||||
itkgTwoImageToOneImageMetric.h
|
||||
itkgTwoImageToOneImageMetric.hxx
|
||||
itkMutualInformationTwoImageToOneImageMetric.h
|
||||
itkMutualInformationTwoImageToOneImageMetric.hxx
|
||||
itkNormalizedCorrelationTwoImageToOneImageMetric.h
|
||||
itkNormalizedCorrelationTwoImageToOneImageMetric.hxx
|
||||
itkTwoProjectionImageRegistrationMethod.h
|
||||
itkTwoProjectionImageRegistrationMethod.hxx
|
||||
)
|
||||
|
||||
ADD_LIBRARY(${LIB_NAME} ${HDR})
|
||||
|
||||
target_link_libraries(
|
||||
${LIB_NAME}
|
||||
${ITK_LIBRARIES}
|
||||
)
|
||||
|
||||
target_include_directories (${LIB_NAME}
|
||||
PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
|
||||
|
||||
|
||||
93
reg23Topograms/itkDTRrecon/autoreg/itkDRTHelpers.h
Normal file
93
reg23Topograms/itkDTRrecon/autoreg/itkDRTHelpers.h
Normal file
@ -0,0 +1,93 @@
|
||||
#ifndef ITKDRTHELPERS_H
|
||||
#define ITKDRTHELPERS_H
|
||||
|
||||
#include "itkMath.h"
|
||||
#include <cstddef>
|
||||
#include <cstring>
|
||||
#include <iterator>
|
||||
|
||||
//Function to get method and class name for logging purposes.
|
||||
template <size_t FL, size_t PFL>
|
||||
const char* computeMethodName(const char (&function)[FL], const char (&prettyFunction)[PFL])
|
||||
{
|
||||
using reverse_ptr = std::reverse_iterator<const char*>;
|
||||
thread_local static char result[PFL];
|
||||
const char* locFuncName = std::search(prettyFunction, prettyFunction + PFL - 1, function, function + FL - 1);
|
||||
const char* locClassName = std::find(reverse_ptr(locFuncName), reverse_ptr(prettyFunction), ' ').base();
|
||||
const char* endFuncName = std::find(locFuncName, prettyFunction + PFL - 1, '(');
|
||||
result[0] = '\0';
|
||||
std::strncat(result, locClassName, endFuncName - locClassName + 1);
|
||||
std::strcat(result, ")");
|
||||
return result;
|
||||
}
|
||||
#define __COMPACT_PRETTY_FUNCTION__ computeMethodName(__FUNCTION__, __PRETTY_FUNCTION__)
|
||||
|
||||
namespace itk {
|
||||
/* reference string required for comparison with tag values */
|
||||
static const char* ImageOrientationStrings[] = {
|
||||
"NotDefined",
|
||||
"HFS",
|
||||
"FFS",
|
||||
"HFP",
|
||||
"FFP",
|
||||
};
|
||||
|
||||
// constant for converting degrees into radians
|
||||
const float dtr = itk::Math::pi_over_180;
|
||||
|
||||
static const bool verbose = true;
|
||||
|
||||
/* this is in the end a IEC to HFS...
|
||||
* but we keep this for ourselves...
|
||||
*/
|
||||
static double Standard_DRT2LPS[9] = {
|
||||
1, 0, 0,
|
||||
0, 0, -1,
|
||||
0, 1, 0
|
||||
};
|
||||
|
||||
static double PAElementsIEC[9] = {
|
||||
1, 0, 0,
|
||||
0, -1, 0,
|
||||
0, 0, -1
|
||||
};
|
||||
|
||||
static double LATElementsIEC[9] = {
|
||||
0, 0, -1,
|
||||
0, -1, 0,
|
||||
-1, 0, 0
|
||||
};
|
||||
|
||||
static double HFS2IEC[9] = {
|
||||
1, 0, 0,
|
||||
0, 0, 1,
|
||||
0, -1, 0
|
||||
};
|
||||
|
||||
static double FFS2IEC[9] = {
|
||||
-1, 0, 0,
|
||||
0, 0, -1,
|
||||
0, -1, 0
|
||||
};
|
||||
|
||||
static double HFP2IEC[9] = {
|
||||
-1, 0, 0,
|
||||
0, 0, 1,
|
||||
0, 1, 0
|
||||
};
|
||||
|
||||
static double FFP2IEC[9] = {
|
||||
1, 0, 0,
|
||||
0, 0, -1,
|
||||
0, 1, 0
|
||||
};
|
||||
|
||||
static double PAT2IEC[9] = {
|
||||
1, 0, 0,
|
||||
0, 0, 1,
|
||||
0, -1, 0
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif // ITKDRTHELPERS_H
|
||||
@ -0,0 +1,117 @@
|
||||
/*=========================================================================
|
||||
*
|
||||
* 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 itkMutualInformationTwoImageToOneImageMetric_h
|
||||
#define itkMutualInformationTwoImageToOneImageMetric_h
|
||||
|
||||
#include "itkCovariantVector.h"
|
||||
#include "itkMacro.h"
|
||||
#include "itkPoint.h"
|
||||
#include "itkgTwoImageToOneImageMetric.h"
|
||||
|
||||
namespace itk {
|
||||
/** \class NormalizedCorrelationTwoImageToOneImageMetric
|
||||
* \brief Computes similarity between two fixed images and one moving image
|
||||
*
|
||||
* This metric computes the correlation between pixels in the two fixed images
|
||||
* and pixels in the moving image. The spatial correspondance between
|
||||
* two fixed images and the moving image is established through a Transform. Pixel values are
|
||||
* taken from the fixed images, their positions are mapped to the moving
|
||||
* image and result in general in non-grid position on it. Values at these
|
||||
* non-grid position of the moving image are interpolated using user-selected
|
||||
* Interpolators. The correlation is normalized by the autocorrelations of both
|
||||
* the fixed and moving images.
|
||||
*
|
||||
* \ingroup RegistrationMetrics
|
||||
* \ingroup TwoProjectionRegistration
|
||||
*/
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
class MutualInformationTwoImageToOneImageMetric : public gTwoImageToOneImageMetric<TFixedImage, TMovingImage> {
|
||||
public:
|
||||
ITK_DISALLOW_COPY_AND_ASSIGN(MutualInformationTwoImageToOneImageMetric);
|
||||
|
||||
/** Standard class type alias. */
|
||||
using Self = MutualInformationTwoImageToOneImageMetric;
|
||||
using Superclass = gTwoImageToOneImageMetric<TFixedImage, TMovingImage>;
|
||||
|
||||
using Pointer = SmartPointer<Self>;
|
||||
using ConstPointer = SmartPointer<const Self>;
|
||||
|
||||
/** Method for creation through the object factory. */
|
||||
itkNewMacro(Self);
|
||||
|
||||
/** Run-time type information (and related methods). */
|
||||
itkTypeMacro(MutualInformationTwoImageToOneImageMetric, Object);
|
||||
|
||||
/** Types transferred from the base class */
|
||||
using RealType = typename Superclass::RealType;
|
||||
using TransformType = typename Superclass::TransformType;
|
||||
using TransformPointer = typename Superclass::TransformPointer;
|
||||
using TransformParametersType = typename Superclass::TransformParametersType;
|
||||
using TransformJacobianType = typename Superclass::TransformJacobianType;
|
||||
using GradientPixelType = typename Superclass::GradientPixelType;
|
||||
|
||||
using MeasureType = typename Superclass::MeasureType;
|
||||
using DerivativeType = typename Superclass::DerivativeType;
|
||||
using FixedImageType = typename Superclass::FixedImageType;
|
||||
using MovingImageType = typename Superclass::MovingImageType;
|
||||
using FixedImageConstPointer = typename Superclass::FixedImageConstPointer;
|
||||
using MovingImageConstPointer = typename Superclass::MovingImageConstPointer;
|
||||
|
||||
/** Get the derivatives of the match measure. */
|
||||
void
|
||||
GetDerivative(const TransformParametersType& parameters, DerivativeType& Derivative) const override;
|
||||
|
||||
/** Get the value for single valued optimizers. */
|
||||
MeasureType
|
||||
GetValue(const TransformParametersType& parameters) const override;
|
||||
|
||||
/** Get the value using the current transforms. */
|
||||
MeasureType
|
||||
GetValue() const;
|
||||
|
||||
/** Get value and derivatives for multiple valued optimizers. */
|
||||
void
|
||||
GetValueAndDerivative(const TransformParametersType& parameters,
|
||||
MeasureType& Value,
|
||||
DerivativeType& Derivative) const override;
|
||||
|
||||
/** Set/Get SubtractMean boolean. If true, the sample mean is subtracted
|
||||
* from the sample values in the cross-correlation formula and
|
||||
* typically results in narrower valleys in the cost fucntion.
|
||||
* Default value is false. */
|
||||
itkSetMacro(SubtractMean, bool);
|
||||
itkGetConstReferenceMacro(SubtractMean, bool);
|
||||
itkBooleanMacro(SubtractMean);
|
||||
|
||||
protected:
|
||||
MutualInformationTwoImageToOneImageMetric();
|
||||
~MutualInformationTwoImageToOneImageMetric() override = default;
|
||||
void
|
||||
PrintSelf(std::ostream& os, Indent indent) const override;
|
||||
|
||||
private:
|
||||
bool m_SubtractMean;
|
||||
};
|
||||
|
||||
} // end namespace itk
|
||||
|
||||
#ifndef ITK_MANUAL_INSTANTIATION
|
||||
#include "itkMutualInformationTwoImageToOneImageMetric.hxx"
|
||||
#endif
|
||||
|
||||
#endif
|
||||
@ -0,0 +1,201 @@
|
||||
/*=========================================================================
|
||||
*
|
||||
* 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 itkMutualInformationTwoImageToOneImageMetric_hxx
|
||||
#define itkMutualInformationTwoImageToOneImageMetric_hxx
|
||||
|
||||
#include "itkImageRegionConstIteratorWithIndex.h"
|
||||
#include "itkMattesMutualInformationImageToImageMetric.h"
|
||||
#include "itkMutualInformationTwoImageToOneImageMetric.h"
|
||||
#include "itkNearestNeighborInterpolateImageFunction.h"
|
||||
#include "itkNormalizeImageFilter.h"
|
||||
#include "itkTranslationTransform.h"
|
||||
|
||||
#include <limits>
|
||||
|
||||
namespace itk {
|
||||
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
MutualInformationTwoImageToOneImageMetric<TFixedImage,
|
||||
TMovingImage>::MutualInformationTwoImageToOneImageMetric()
|
||||
{
|
||||
m_SubtractMean = false;
|
||||
}
|
||||
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
typename MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::MeasureType
|
||||
MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValue(
|
||||
const TransformParametersType& parameters) const
|
||||
{
|
||||
FixedImageConstPointer fixedImage1 = this->m_FixedImage1;
|
||||
|
||||
if (!fixedImage1) {
|
||||
itkExceptionMacro(<< "Fixed image1 has not been assigned");
|
||||
}
|
||||
|
||||
FixedImageConstPointer fixedImage2 = this->m_FixedImage2;
|
||||
|
||||
if (!fixedImage2) {
|
||||
itkExceptionMacro(<< "Fixed image2 has not been assigned");
|
||||
}
|
||||
|
||||
if (!this->m_Filter1) {
|
||||
itkExceptionMacro(<< "Fixed filter1 has not been assigned");
|
||||
}
|
||||
|
||||
if (!this->m_Filter2) {
|
||||
itkExceptionMacro(<< "Fixed filter1 has not been assigned");
|
||||
}
|
||||
|
||||
this->SetTransformParameters(parameters);
|
||||
|
||||
return GetValue();
|
||||
}
|
||||
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
typename MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::MeasureType
|
||||
MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValue() const
|
||||
{
|
||||
|
||||
FixedImageConstPointer fixedImage1 = this->m_FixedImage1;
|
||||
|
||||
if (!fixedImage1) {
|
||||
itkExceptionMacro(<< "Fixed image1 has not been assigned");
|
||||
}
|
||||
|
||||
FixedImageConstPointer fixedImage2 = this->m_FixedImage2;
|
||||
|
||||
if (!fixedImage2) {
|
||||
itkExceptionMacro(<< "Fixed image2 has not been assigned");
|
||||
}
|
||||
|
||||
if (!this->m_Filter1) {
|
||||
itkExceptionMacro(<< "Fixed filter1 has not been assigned");
|
||||
}
|
||||
|
||||
if (!this->m_Filter2) {
|
||||
itkExceptionMacro(<< "Fixed filter1 has not been assigned");
|
||||
}
|
||||
|
||||
this->m_Interpolator1->SetTransform(this->m_Transform1);
|
||||
this->m_Interpolator1->Initialize();
|
||||
this->m_Filter1->Update();
|
||||
|
||||
using InternalMetricType = itk::MattesMutualInformationImageToImageMetric<TFixedImage, TMovingImage>;
|
||||
|
||||
using NearestNeighborType = itk::NearestNeighborInterpolateImageFunction<TMovingImage, double>;
|
||||
|
||||
//We need to set transform and parameters.
|
||||
auto dummyTransform = TransformType::New();
|
||||
auto dummyParameters = dummyTransform->GetParameters();
|
||||
auto dummyInterpolator1 = NearestNeighborType::New();
|
||||
|
||||
auto metric1 = InternalMetricType::New();
|
||||
metric1->SetTransform(dummyTransform);
|
||||
metric1->SetTransformParameters(dummyParameters);
|
||||
metric1->SetInterpolator(dummyInterpolator1);
|
||||
|
||||
auto fixedImageRegion1 = fixedImage1->GetBufferedRegion();
|
||||
metric1->SetFixedImageRegion(fixedImageRegion1);
|
||||
|
||||
auto movingImageRegion = this->m_Filter1->GetOutput()->GetBufferedRegion();
|
||||
unsigned int numberOfPixels = movingImageRegion.GetNumberOfPixels();
|
||||
|
||||
//auto numberOfSamples = static_cast<unsigned int>(numberOfPixels * 0.50); //100%
|
||||
metric1->UseAllPixelsOn();
|
||||
metric1->SetNumberOfHistogramBins(30);
|
||||
|
||||
metric1->SetFixedImage(fixedImage1);
|
||||
metric1->SetMovingImage(this->m_Filter1->GetOutput());
|
||||
|
||||
metric1->Initialize();
|
||||
|
||||
auto measure1 = metric1->GetValue(dummyParameters);
|
||||
|
||||
this->m_Interpolator2->SetTransform(this->m_Transform2);
|
||||
this->m_Interpolator2->Initialize();
|
||||
this->m_Filter2->Update();
|
||||
|
||||
auto dummyInterpolator2 = NearestNeighborType::New();
|
||||
auto metric2 = InternalMetricType::New();
|
||||
metric2->SetTransform(dummyTransform);
|
||||
metric2->SetTransformParameters(dummyParameters);
|
||||
metric2->SetInterpolator(dummyInterpolator2);
|
||||
|
||||
auto fixedImageRegion2 = fixedImage1->GetBufferedRegion();
|
||||
metric2->SetFixedImageRegion(fixedImageRegion2);
|
||||
|
||||
movingImageRegion = this->m_Filter2->GetOutput()->GetBufferedRegion();
|
||||
numberOfPixels = movingImageRegion.GetNumberOfPixels();
|
||||
|
||||
//numberOfSamples = static_cast<unsigned int>(numberOfPixels * 0.50); //100%
|
||||
//%metric2->SetNumberOfSpatialSamples(numberOfSamples);
|
||||
metric2->UseAllPixelsOn();
|
||||
metric2->SetNumberOfHistogramBins(30);
|
||||
|
||||
metric2->SetFixedImage(fixedImage2);
|
||||
metric2->SetMovingImage(this->m_Filter2->GetOutput());
|
||||
|
||||
metric2->Initialize();
|
||||
|
||||
auto measure2 = metric2->GetValue(dummyParameters);
|
||||
|
||||
//std::cout << "movingValueMin: " << movingValueMin << " movingValueMax: " << movingValueMax << std::endl;
|
||||
//std::cout << "fixedValueMin: " << fixedValueMin << " fixedValueMax: " << fixedValueMax << std::endl;
|
||||
|
||||
// Calculate the measure value between fixed image 2 and the moving image
|
||||
|
||||
auto measure = (measure1 + measure2) / 2.0;
|
||||
|
||||
auto oldprecision = std::cout.precision();
|
||||
std::cout.precision(std::numeric_limits<double>::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;
|
||||
}
|
||||
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
void MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetDerivative(
|
||||
const TransformParametersType& itkNotUsed(parameters),
|
||||
DerivativeType& itkNotUsed(derivative)) const
|
||||
{
|
||||
// under construction
|
||||
}
|
||||
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
void MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValueAndDerivative(
|
||||
const TransformParametersType& itkNotUsed(parameters),
|
||||
MeasureType& itkNotUsed(value),
|
||||
DerivativeType& itkNotUsed(derivative)) const
|
||||
{
|
||||
// under construction
|
||||
}
|
||||
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
void MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::PrintSelf(std::ostream& os,
|
||||
Indent indent) const
|
||||
{
|
||||
Superclass::PrintSelf(os, indent);
|
||||
os << indent << "SubtractMean: " << m_SubtractMean << std::endl;
|
||||
}
|
||||
|
||||
} // end namespace itk
|
||||
|
||||
#endif
|
||||
@ -0,0 +1,127 @@
|
||||
/*=========================================================================
|
||||
*
|
||||
* 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 itkNormalizedCorrelationTwoImageToOneImageMetric_h
|
||||
#define itkNormalizedCorrelationTwoImageToOneImageMetric_h
|
||||
|
||||
#include "itkCovariantVector.h"
|
||||
#include "itkMacro.h"
|
||||
#include "itkNormalizedCorrelationImageToImageMetric.hxx"
|
||||
#include "itkPoint.h"
|
||||
#include "itkgTwoImageToOneImageMetric.h"
|
||||
|
||||
namespace itk {
|
||||
/** \class NormalizedCorrelationTwoImageToOneImageMetric
|
||||
* \brief Computes similarity between two fixed images and one moving image
|
||||
*
|
||||
* This metric computes the correlation between pixels in the two fixed images
|
||||
* and pixels in the moving image. The spatial correspondance between
|
||||
* two fixed images and the moving image is established through a Transform. Pixel values are
|
||||
* taken from the fixed images, their positions are mapped to the moving
|
||||
* image and result in general in non-grid position on it. Values at these
|
||||
* non-grid position of the moving image are interpolated using user-selected
|
||||
* Interpolators. The correlation is normalized by the autocorrelations of both
|
||||
* the fixed and moving images.
|
||||
*
|
||||
* \ingroup RegistrationMetrics
|
||||
* \ingroup TwoProjectionRegistration
|
||||
*/
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
class NormalizedCorrelationTwoImageToOneImageMetric : public gTwoImageToOneImageMetric<TFixedImage, TMovingImage> {
|
||||
public:
|
||||
ITK_DISALLOW_COPY_AND_ASSIGN(NormalizedCorrelationTwoImageToOneImageMetric);
|
||||
|
||||
/** Standard class type alias. */
|
||||
using Self = NormalizedCorrelationTwoImageToOneImageMetric;
|
||||
using Superclass = gTwoImageToOneImageMetric<TFixedImage, TMovingImage>;
|
||||
|
||||
using Pointer = SmartPointer<Self>;
|
||||
using ConstPointer = SmartPointer<const Self>;
|
||||
|
||||
/** Method for creation through the object factory. */
|
||||
itkNewMacro(Self);
|
||||
|
||||
/** Run-time type information (and related methods). */
|
||||
itkTypeMacro(NormalizedCorrelationTwoImageToOneImageMetric, Object);
|
||||
|
||||
/** Types transferred from the base class */
|
||||
using RealType = typename Superclass::RealType;
|
||||
using TransformType = typename Superclass::TransformType;
|
||||
using TransformPointer = typename Superclass::TransformPointer;
|
||||
using TransformParametersType = typename Superclass::TransformParametersType;
|
||||
using TransformJacobianType = typename Superclass::TransformJacobianType;
|
||||
using GradientPixelType = typename Superclass::GradientPixelType;
|
||||
|
||||
using MeasureType = typename Superclass::MeasureType;
|
||||
using DerivativeType = typename Superclass::DerivativeType;
|
||||
|
||||
using FixedImageType = typename Superclass::FixedImageType;
|
||||
using FixedImageConstPointer = typename Superclass::FixedImageConstPointer;
|
||||
using FixedImageMaskPointer = typename Superclass::FixedImageMaskPointer;
|
||||
using FixedImageRegionType = typename Superclass::FixedImageRegionType;
|
||||
|
||||
using MovingImageType = typename Superclass::MovingImageType;
|
||||
using MovingImageConstPointer = typename Superclass::MovingImageConstPointer;
|
||||
using MovingImageMaskPointer = typename Superclass::MovingImageMaskPointer;
|
||||
|
||||
using ChangeInformationFilterPointer = typename Superclass::ChangeInformationFilterPointer;
|
||||
|
||||
/** Get the derivatives of the match measure. */
|
||||
void
|
||||
GetDerivative(const TransformParametersType& parameters, DerivativeType& Derivative) const override;
|
||||
|
||||
/** Get the value for single valued optimizers, after applying paramters to the transforms. */
|
||||
MeasureType
|
||||
GetValue(const TransformParametersType& parameters) const override;
|
||||
|
||||
/** Get the value using the current transforms. */
|
||||
MeasureType
|
||||
GetValue() const;
|
||||
|
||||
/** Get value and derivatives for multiple valued optimizers. */
|
||||
void
|
||||
GetValueAndDerivative(const TransformParametersType& parameters,
|
||||
MeasureType& Value,
|
||||
DerivativeType& Derivative) const override;
|
||||
|
||||
/** Set/Get SubtractMean boolean. If true, the sample mean is subtracted
|
||||
* from the sample values in the cross-correlation formula and
|
||||
* typically results in narrower valleys in the cost fucntion.
|
||||
* Default value is false. */
|
||||
itkSetMacro(SubtractMean, bool);
|
||||
itkGetConstReferenceMacro(SubtractMean, bool);
|
||||
itkBooleanMacro(SubtractMean);
|
||||
|
||||
protected:
|
||||
NormalizedCorrelationTwoImageToOneImageMetric();
|
||||
~NormalizedCorrelationTwoImageToOneImageMetric() override = default;
|
||||
void
|
||||
PrintSelf(std::ostream& os, Indent indent) const override;
|
||||
|
||||
private:
|
||||
bool m_SubtractMean;
|
||||
|
||||
MeasureType CalculateMeasure(int imageId) const;
|
||||
};
|
||||
|
||||
} // end namespace itk
|
||||
|
||||
#ifndef ITK_MANUAL_INSTANTIATION
|
||||
#include "itkNormalizedCorrelationTwoImageToOneImageMetric.hxx"
|
||||
#endif
|
||||
|
||||
#endif
|
||||
@ -0,0 +1,321 @@
|
||||
/*=========================================================================
|
||||
*
|
||||
* 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 itkNormalizedCorrelationTwoImageToOneImageMetric_hxx
|
||||
#define itkNormalizedCorrelationTwoImageToOneImageMetric_hxx
|
||||
|
||||
#include "itkImageRegionConstIteratorWithIndex.h"
|
||||
#include "itkNormalizedCorrelationTwoImageToOneImageMetric.h"
|
||||
|
||||
#
|
||||
|
||||
namespace itk {
|
||||
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage,
|
||||
TMovingImage>::NormalizedCorrelationTwoImageToOneImageMetric()
|
||||
{
|
||||
m_SubtractMean = false;
|
||||
}
|
||||
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
typename NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::MeasureType
|
||||
NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::CalculateMeasure(int imageId) const
|
||||
{
|
||||
|
||||
using FixedIteratorType = itk::ImageRegionConstIteratorWithIndex<FixedImageType>;
|
||||
using MovingIteratorType = itk::ImageRegionConstIteratorWithIndex<MovingImageType>;
|
||||
using AccumulateType = typename NumericTraits<MeasureType>::AccumulateType;
|
||||
|
||||
FixedImageConstPointer fixedImage;
|
||||
ChangeInformationFilterPointer filter;
|
||||
FixedImageMaskPointer fixedImageMask;
|
||||
FixedImageRegionType fixedImageRegion;
|
||||
MeasureType measure;
|
||||
|
||||
typename FixedImageType::IndexType indexFixed;
|
||||
typename MovingImageType::IndexType indexMoving;
|
||||
typename Superclass::InputPointType inputPoint;
|
||||
|
||||
if (imageId == 1) {
|
||||
this->m_Interpolator1->SetTransform(this->m_Transform1);
|
||||
this->m_Interpolator1->Initialize();
|
||||
|
||||
fixedImage = this->m_FixedImage1;
|
||||
fixedImageMask = this->m_FixedImageMask1;
|
||||
fixedImageRegion = this->GetFixedImageRegion1();
|
||||
filter = this->m_Filter1;
|
||||
} else if (imageId == 2) {
|
||||
this->m_Interpolator2->SetTransform(this->m_Transform2);
|
||||
this->m_Interpolator2->Initialize();
|
||||
|
||||
fixedImage = this->m_FixedImage2;
|
||||
fixedImageMask = this->m_FixedImageMask2;
|
||||
fixedImageRegion = this->GetFixedImageRegion2();
|
||||
filter = this->m_Filter2;
|
||||
} else {
|
||||
itkExceptionMacro(<< "Non supported image id.");
|
||||
}
|
||||
|
||||
FixedIteratorType ti(fixedImage, fixedImageRegion);
|
||||
|
||||
this->m_NumberOfPixelsCounted = 0;
|
||||
int numberOfPixelsVisited = 0;
|
||||
|
||||
//filter->Update();
|
||||
filter->Update();
|
||||
MovingImageConstPointer imageDRTIn = filter->GetOutput();
|
||||
//this->WriteImage(imageDRTIn, "", "");
|
||||
|
||||
AccumulateType sff = NumericTraits<AccumulateType>::ZeroValue();
|
||||
AccumulateType smm = NumericTraits<AccumulateType>::ZeroValue();
|
||||
AccumulateType sfm = NumericTraits<AccumulateType>::ZeroValue();
|
||||
AccumulateType sf = NumericTraits<AccumulateType>::ZeroValue();
|
||||
AccumulateType sm = NumericTraits<AccumulateType>::ZeroValue();
|
||||
|
||||
RealType movingValueMin = NumericTraits<RealType>::max();
|
||||
RealType fixedValueMin = NumericTraits<RealType>::max();
|
||||
RealType movingValueMax = NumericTraits<RealType>::NonpositiveMin();
|
||||
RealType fixedValueMax = NumericTraits<RealType>::NonpositiveMin();
|
||||
|
||||
RealType maxX = NumericTraits<RealType>::NonpositiveMin();
|
||||
RealType maxY = NumericTraits<RealType>::NonpositiveMin();
|
||||
RealType maxZ = NumericTraits<RealType>::NonpositiveMin();
|
||||
RealType minX = NumericTraits<RealType>::max();
|
||||
RealType minY = NumericTraits<RealType>::max();
|
||||
RealType minZ = NumericTraits<RealType>::max();
|
||||
|
||||
AccumulateType m_meanCurr_f = NumericTraits<AccumulateType>::ZeroValue();
|
||||
AccumulateType m_meanCurr_m = NumericTraits<AccumulateType>::ZeroValue();
|
||||
AccumulateType m_meanPrev_f = NumericTraits<AccumulateType>::ZeroValue();
|
||||
AccumulateType m_meanPrev_m = NumericTraits<AccumulateType>::ZeroValue();
|
||||
AccumulateType m_sPrev_f = NumericTraits<AccumulateType>::ZeroValue();
|
||||
AccumulateType m_sPrev_m = NumericTraits<AccumulateType>::ZeroValue();
|
||||
AccumulateType m_sCurr_f = NumericTraits<AccumulateType>::ZeroValue();
|
||||
AccumulateType m_sCurr_m = NumericTraits<AccumulateType>::ZeroValue();
|
||||
AccumulateType m_varianceCurr_f = NumericTraits<AccumulateType>::ZeroValue();
|
||||
AccumulateType m_varianceCurr_m = NumericTraits<AccumulateType>::ZeroValue();
|
||||
|
||||
while (!ti.IsAtEnd()) {
|
||||
|
||||
indexFixed = ti.GetIndex();
|
||||
++numberOfPixelsVisited;
|
||||
|
||||
fixedImage->TransformIndexToPhysicalPoint(indexFixed, inputPoint);
|
||||
|
||||
if (fixedImageMask && !fixedImageMask->IsInsideInWorldSpace(inputPoint)) {
|
||||
++ti;
|
||||
continue;
|
||||
}
|
||||
|
||||
if (this->m_MovingImageMask && !this->m_MovingImageMask->IsInsideInWorldSpace(inputPoint)) {
|
||||
++ti;
|
||||
continue;
|
||||
}
|
||||
|
||||
const RealType fixedValue = ti.Get();
|
||||
|
||||
if (fixedValue == -1024) {
|
||||
++ti;
|
||||
continue;
|
||||
}
|
||||
|
||||
if (imageDRTIn->TransformPhysicalPointToIndex(inputPoint, indexMoving)) {
|
||||
|
||||
double x = inputPoint[0];
|
||||
double y = inputPoint[1];
|
||||
double z = inputPoint[2];
|
||||
maxX = maxX > x ? maxX : x;
|
||||
minX = minX < x ? minX : x;
|
||||
maxY = maxY > y ? maxY : y;
|
||||
minY = minY < y ? minY : y;
|
||||
maxZ = maxZ > z ? maxZ : z;
|
||||
minZ = minZ < z ? minZ : z;
|
||||
|
||||
const RealType movingValue = imageDRTIn->GetPixel(indexMoving);
|
||||
|
||||
if (NumericTraits<short>::NonpositiveMin() == movingValue) {
|
||||
++ti;
|
||||
continue;
|
||||
}
|
||||
|
||||
movingValueMin = movingValue < movingValueMin ? movingValue : movingValueMin;
|
||||
movingValueMax = movingValue > movingValueMax ? movingValue : movingValueMax;
|
||||
fixedValueMin = fixedValue < fixedValueMin ? fixedValue : fixedValueMin;
|
||||
fixedValueMax = fixedValue > fixedValueMax ? fixedValue : fixedValueMax;
|
||||
|
||||
sff += fixedValue * fixedValue;
|
||||
smm += movingValue * movingValue;
|
||||
sfm += fixedValue * movingValue;
|
||||
if (this->m_SubtractMean) {
|
||||
sf += fixedValue;
|
||||
sm += movingValue;
|
||||
}
|
||||
this->m_NumberOfPixelsCounted++;
|
||||
|
||||
if (this->m_NumberOfPixelsCounted == 1) {
|
||||
m_meanCurr_f = fixedValue;
|
||||
m_meanCurr_m = movingValue;
|
||||
} else {
|
||||
m_meanPrev_f = m_meanCurr_f;
|
||||
m_meanPrev_m = m_meanCurr_m;
|
||||
m_sPrev_f = m_sCurr_f;
|
||||
m_sPrev_m = m_sCurr_m;
|
||||
|
||||
m_meanCurr_f = m_meanPrev_f + (fixedValue - m_meanPrev_f) / this->m_NumberOfPixelsCounted;
|
||||
m_sCurr_f = m_sPrev_f + (fixedValue - m_meanPrev_f) * (fixedValue - m_meanCurr_f);
|
||||
m_varianceCurr_f = m_sCurr_f / (this->m_NumberOfPixelsCounted - 1);
|
||||
|
||||
m_meanCurr_m = m_meanPrev_m + (movingValue - m_meanPrev_m) / this->m_NumberOfPixelsCounted;
|
||||
m_sCurr_m = m_sPrev_m + (movingValue - m_meanPrev_m) * (movingValue - m_meanCurr_m);
|
||||
m_varianceCurr_m = m_sCurr_m / (this->m_NumberOfPixelsCounted - 1);
|
||||
}
|
||||
}
|
||||
|
||||
++ti;
|
||||
}
|
||||
|
||||
// sfm = sfm - m_meanCurr_f * sm - m_meanCurr_m * sf + this->m_NumberOfPixelsCounted * m_meanCurr_m * m_meanCurr_f;
|
||||
// const RealType denom = -1 * sqrt(m_varianceCurr_f * m_varianceCurr_m) * this->m_NumberOfPixelsCounted;
|
||||
|
||||
// if (this->m_NumberOfPixelsCounted > 0 && denom != 0.0) {
|
||||
// measure = sfm / denom;
|
||||
// } else {
|
||||
// measure = NumericTraits<MeasureType>::ZeroValue();
|
||||
// }
|
||||
|
||||
if (this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0) {
|
||||
sff -= (sf * sf / this->m_NumberOfPixelsCounted);
|
||||
smm -= (sm * sm / this->m_NumberOfPixelsCounted);
|
||||
sfm -= (sf * sm / this->m_NumberOfPixelsCounted);
|
||||
}
|
||||
|
||||
RealType denom = -1.0 * sqrt(sff * smm);
|
||||
|
||||
if (this->m_NumberOfPixelsCounted > 0 && denom != 0.0) {
|
||||
measure = sfm / denom;
|
||||
} else {
|
||||
measure = NumericTraits<MeasureType>::Zero;
|
||||
}
|
||||
std::cout << "movingValueMin: " << movingValueMin << " movingValueMax: " << movingValueMax << std::endl;
|
||||
std::cout << "fixedValueMin: " << fixedValueMin << " fixedValueMax: " << fixedValueMax << std::endl;
|
||||
|
||||
std::cout << "max coordinates: " << maxX << " " << maxX << " " << maxZ << std::endl;
|
||||
std::cout << "min coordinates: " << minX << " " << minY << " " << minZ << std::endl;
|
||||
|
||||
return measure;
|
||||
}
|
||||
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
typename NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::MeasureType
|
||||
NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValue(
|
||||
const TransformParametersType& parameters) const
|
||||
{
|
||||
if (!this->m_FixedImage1) {
|
||||
itkExceptionMacro(<< "Fixed image1 has not been assigned");
|
||||
}
|
||||
|
||||
if (!this->m_FixedImage2) {
|
||||
itkExceptionMacro(<< "Fixed image2 has not been assigned");
|
||||
}
|
||||
|
||||
if (!this->m_Filter1) {
|
||||
itkExceptionMacro(<< "Fixed filter1 has not been assigned");
|
||||
}
|
||||
|
||||
if (!this->m_Filter2) {
|
||||
itkExceptionMacro(<< "Fixed filter1 has not been assigned");
|
||||
}
|
||||
|
||||
if (!this->SetTransformParameters(parameters)) {
|
||||
std::cout << "Measure: 0.0, transform Invalid!" << std::endl;
|
||||
std::cout << "=======================================================" << std::endl;
|
||||
|
||||
return 0.0;
|
||||
}
|
||||
return GetValue();
|
||||
}
|
||||
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
typename NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::MeasureType
|
||||
NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValue() const
|
||||
{
|
||||
|
||||
if (!this->m_FixedImage1) {
|
||||
itkExceptionMacro(<< "Fixed image1 has not been assigned");
|
||||
}
|
||||
|
||||
if (!this->m_FixedImage2) {
|
||||
itkExceptionMacro(<< "Fixed image2 has not been assigned");
|
||||
}
|
||||
|
||||
if (!this->m_Filter1) {
|
||||
itkExceptionMacro(<< "Fixed filter1 has not been assigned");
|
||||
}
|
||||
|
||||
if (!this->m_Filter2) {
|
||||
itkExceptionMacro(<< "Fixed filter1 has not been assigned");
|
||||
}
|
||||
|
||||
// Calculate the measure value between fixed image 1 and the moving image
|
||||
auto oldprecision = std::cout.precision();
|
||||
|
||||
MeasureType measure1 = CalculateMeasure(1);
|
||||
std::cout.precision(std::numeric_limits<double>::digits10 + 2);
|
||||
std::cout << "Measure1: " << measure1 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl;
|
||||
std::cout.precision(oldprecision);
|
||||
|
||||
MeasureType measure2 = CalculateMeasure(2);
|
||||
std::cout.precision(std::numeric_limits<double>::digits10 + 2);
|
||||
std::cout << "Measure2: " << measure2 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl;
|
||||
|
||||
auto measure = (measure1 + measure2) / 2.0;
|
||||
std::cout << "Measure: " << measure << std::endl;
|
||||
std::cout << "=======================================================" << std::endl;
|
||||
std::cout.precision(oldprecision);
|
||||
|
||||
return measure;
|
||||
}
|
||||
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
void NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetDerivative(
|
||||
const TransformParametersType& itkNotUsed(parameters),
|
||||
DerivativeType& itkNotUsed(derivative)) const
|
||||
{
|
||||
// under construction
|
||||
}
|
||||
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
void NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValueAndDerivative(
|
||||
const TransformParametersType& itkNotUsed(parameters),
|
||||
MeasureType& itkNotUsed(value),
|
||||
DerivativeType& itkNotUsed(derivative)) const
|
||||
{
|
||||
// under construction
|
||||
}
|
||||
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
void NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::PrintSelf(std::ostream& os,
|
||||
Indent indent) const
|
||||
{
|
||||
Superclass::PrintSelf(os, indent);
|
||||
os << indent << "SubtractMean: " << m_SubtractMean << std::endl;
|
||||
}
|
||||
|
||||
} // end namespace itk
|
||||
|
||||
#endif
|
||||
@ -0,0 +1,269 @@
|
||||
/*=========================================================================
|
||||
*
|
||||
* 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_h
|
||||
#define itkTwoProjectionImageRegistrationMethod_h
|
||||
|
||||
#include "itkChangeInformationImageFilter.h"
|
||||
#include "itkDataObjectDecorator.h"
|
||||
#include "itkImage.h"
|
||||
#include "itkMacro.h"
|
||||
#include "itkProcessObject.h"
|
||||
#include "itkSingleValuedNonLinearOptimizer.h"
|
||||
#include "itkgTwoImageToOneImageMetric.h"
|
||||
#include <iostream>
|
||||
|
||||
namespace itk {
|
||||
|
||||
/** \class TwoProjectionImageRegistrationMethod
|
||||
* \brief Base class for Image Registration Methods
|
||||
*
|
||||
* This Class define the generic interface for a registration method.
|
||||
*
|
||||
* This class is templated over the type of the two image to be
|
||||
* registered. A generic Transform is used by this class. That allows
|
||||
* to select at run time the particular type of transformation that
|
||||
* is to be applied for registering the images.
|
||||
*
|
||||
* This method use a generic Metric in order to compare the two images.
|
||||
* the final goal of the registration method is to find the set of
|
||||
* parameters of the Transformation that optimizes the metric.
|
||||
*
|
||||
* The registration method also support a generic optimizer that can
|
||||
* be selected at run-time. The only restriction for the optimizer is
|
||||
* that it should be able to operate in single-valued cost functions
|
||||
* given that the metrics used to compare images provide a single
|
||||
* value as output.
|
||||
*
|
||||
* The terms : Fixed image and Moving image are used in this class
|
||||
* to indicate what image is being mapped by the transform.
|
||||
*
|
||||
* This class uses the coordinate system of the Fixed image as a reference
|
||||
* and searchs for a Transform that will map points from the space of the
|
||||
* Fixed image to the space of the Moving image.
|
||||
*
|
||||
* For doing so, a Metric will be continously applied to compare the Fixed
|
||||
* image with the Transformed Moving image. This process also requires to
|
||||
* interpolate values from the Moving image.
|
||||
*
|
||||
* \ingroup RegistrationFilters
|
||||
* \ingroup TwoProjectionRegistration
|
||||
*/
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
class TwoProjectionImageRegistrationMethod : public ProcessObject {
|
||||
public:
|
||||
ITK_DISALLOW_COPY_AND_ASSIGN(TwoProjectionImageRegistrationMethod);
|
||||
|
||||
/** Standard class type alias. */
|
||||
using Self = TwoProjectionImageRegistrationMethod;
|
||||
using Superclass = ProcessObject;
|
||||
using Pointer = SmartPointer<Self>;
|
||||
using ConstPointer = SmartPointer<const Self>;
|
||||
|
||||
/** Method for creation through the object factory. */
|
||||
itkNewMacro(Self);
|
||||
|
||||
/** Run-time type information (and related methods). */
|
||||
itkTypeMacro(TwoProjectionImageRegistrationMethod, ProcessObject);
|
||||
|
||||
/** Type of the Fixed image. */
|
||||
using FixedImageType = TFixedImage;
|
||||
using FixedImageConstPointer = typename FixedImageType::ConstPointer;
|
||||
|
||||
/** Type of the Moving image. */
|
||||
using MovingImageType = TMovingImage;
|
||||
using MovingImageConstPointer = typename MovingImageType::ConstPointer;
|
||||
using ChangeInformationFilterType = ChangeInformationImageFilter<MovingImageType>;
|
||||
using ChangeInformationFilterPointer = typename ChangeInformationFilterType::Pointer;
|
||||
|
||||
/** Type of the metric. */
|
||||
using MetricType = gTwoImageToOneImageMetric<FixedImageType, MovingImageType>;
|
||||
using MetricPointer = typename MetricType::Pointer;
|
||||
using FixedImageRegionType = typename MetricType::FixedImageRegionType;
|
||||
|
||||
/** Type of the Transform . */
|
||||
using TransformType = typename MetricType::TransformType;
|
||||
using TransformPointer = typename TransformType::Pointer;
|
||||
|
||||
/** Type for the output: Using Decorator pattern for enabling
|
||||
* the Transform to be passed in the data pipeline */
|
||||
using TransformOutputType = DataObjectDecorator<TransformType>;
|
||||
using TransformOutputPointer = typename TransformOutputType::Pointer;
|
||||
using TransformOutputConstPointer = typename TransformOutputType::ConstPointer;
|
||||
|
||||
/** Type of the Interpolator. */
|
||||
using InterpolatorType = typename MetricType::InterpolatorType;
|
||||
using InterpolatorPointer = typename InterpolatorType::Pointer;
|
||||
|
||||
/** Type of the optimizer. */
|
||||
using OptimizerType = SingleValuedNonLinearOptimizer;
|
||||
|
||||
/** Type of the Transformation parameters This is the same type used to
|
||||
* represent the search space of the optimization algorithm */
|
||||
using ParametersType = typename MetricType::TransformParametersType;
|
||||
|
||||
/** Smart Pointer type to a DataObject. */
|
||||
using DataObjectPointer = typename DataObject::Pointer;
|
||||
|
||||
/** Method that initiates the registration. This will Initialize and ensure
|
||||
* that all inputs the registration needs are in place, via a call to
|
||||
* Initialize() will then start the optimization process via a call to
|
||||
* StartOptimization() */
|
||||
void
|
||||
StartRegistration();
|
||||
|
||||
/** Method that initiates the optimization process. */
|
||||
void
|
||||
StartOptimization();
|
||||
|
||||
/** Set/Get the Fixed images. */
|
||||
void
|
||||
SetFixedImage1(const FixedImageType* fixedImage1);
|
||||
void
|
||||
SetFixedImage2(const FixedImageType* fixedImage2);
|
||||
itkGetConstObjectMacro(FixedImage1, FixedImageType);
|
||||
itkGetConstObjectMacro(FixedImage2, FixedImageType);
|
||||
|
||||
/** Set/Get the Moving image. */
|
||||
void
|
||||
SetMovingImage(const MovingImageType* movingImage);
|
||||
itkGetConstObjectMacro(MovingImage, MovingImageType);
|
||||
|
||||
/** Set/Get the Optimizer. */
|
||||
itkSetObjectMacro(Optimizer, OptimizerType);
|
||||
itkGetConstObjectMacro(Optimizer, OptimizerType);
|
||||
|
||||
/** Set/Get the Metric. */
|
||||
itkSetObjectMacro(Metric, MetricType);
|
||||
itkGetConstObjectMacro(Metric, MetricType);
|
||||
|
||||
/** Set/Get the Transfrom1. */
|
||||
itkSetObjectMacro(Transform1, TransformType);
|
||||
itkGetConstObjectMacro(Transform1, TransformType);
|
||||
|
||||
/** Set/Get the Transfrom2. */
|
||||
itkSetObjectMacro(Transform2, TransformType);
|
||||
itkGetConstObjectMacro(Transform2, TransformType);
|
||||
|
||||
/** Set/Get the Interpolators. */
|
||||
itkSetObjectMacro(Interpolator1, InterpolatorType);
|
||||
itkSetObjectMacro(Interpolator2, InterpolatorType);
|
||||
itkGetConstObjectMacro(Interpolator1, InterpolatorType);
|
||||
itkGetConstObjectMacro(Interpolator2, InterpolatorType);
|
||||
|
||||
/** Set/Get the Meta informations. */
|
||||
itkSetObjectMacro(TransformMetaInfo, TransformMetaInformation);
|
||||
itkGetConstObjectMacro(TransformMetaInfo, TransformMetaInformation);
|
||||
|
||||
/** Set/Get the output filters. */
|
||||
itkSetObjectMacro(Filter1, ChangeInformationFilterType);
|
||||
itkGetConstObjectMacro(Filter1, ChangeInformationFilterType);
|
||||
itkSetObjectMacro(Filter2, ChangeInformationFilterType);
|
||||
itkGetConstObjectMacro(Filter2, ChangeInformationFilterType);
|
||||
|
||||
/** Set the region of the fixed image to be considered as region of
|
||||
interest during the registration. This region will be passed to
|
||||
the ImageMetric in order to restrict the metric computation to
|
||||
consider only this region.
|
||||
\warning The same region can also be set directly into the metric.
|
||||
please avoid to set the region in both places since this can lead
|
||||
to inconsistent configurations. */
|
||||
void
|
||||
SetFixedImageRegion1(const FixedImageRegionType& region1);
|
||||
void
|
||||
SetFixedImageRegion2(const FixedImageRegionType& region2);
|
||||
/** Get the region of the fixed image to be considered as region of
|
||||
interest during the registration. This region will be passed to
|
||||
the ImageMetric in order to restrict the metric computation to
|
||||
consider only this region. */
|
||||
itkGetConstReferenceMacro(FixedImageRegion1, FixedImageRegionType);
|
||||
itkGetConstReferenceMacro(FixedImageRegion2, FixedImageRegionType);
|
||||
/** True if a region has been defined for the fixed image to which
|
||||
the ImageMetric will limit its computation */
|
||||
itkGetMacro(FixedImageRegionDefined1, bool);
|
||||
itkGetMacro(FixedImageRegionDefined2, bool);
|
||||
/** Turn on/off the use of a fixed image region to which
|
||||
the ImageMetric will limit its computation.
|
||||
\warning The region must have been previously defined using the
|
||||
SetFixedImageRegion member function */
|
||||
itkSetMacro(FixedImageRegionDefined1, bool);
|
||||
itkSetMacro(FixedImageRegionDefined2, bool);
|
||||
|
||||
/** Initialize by setting the interconnects between the components. */
|
||||
virtual void
|
||||
Initialize();
|
||||
|
||||
/** Returns the transform resulting from the registration process */
|
||||
const TransformOutputType*
|
||||
GetOutput() const;
|
||||
|
||||
/** Make a DataObject of the correct type to be used as the specified
|
||||
* output. */
|
||||
using Superclass::MakeOutput;
|
||||
DataObjectPointer
|
||||
MakeOutput(DataObjectPointerArraySizeType idx) override;
|
||||
|
||||
protected:
|
||||
TwoProjectionImageRegistrationMethod();
|
||||
~TwoProjectionImageRegistrationMethod() override = default;
|
||||
void
|
||||
PrintSelf(std::ostream& os, Indent indent) const override;
|
||||
|
||||
/** Method invoked by the pipeline in order to trigger the computation of
|
||||
* the registration. */
|
||||
void
|
||||
GenerateData() override;
|
||||
|
||||
/** Provides derived classes with the ability to set this private var */
|
||||
itkSetMacro(LastOptimizerParameters, ParametersType);
|
||||
|
||||
private:
|
||||
MetricPointer m_Metric;
|
||||
OptimizerType::Pointer m_Optimizer;
|
||||
|
||||
MovingImageConstPointer m_MovingImage;
|
||||
FixedImageConstPointer m_FixedImage1;
|
||||
FixedImageConstPointer m_FixedImage2;
|
||||
|
||||
TransformPointer m_Transform1;
|
||||
TransformPointer m_Transform2;
|
||||
InterpolatorPointer m_Interpolator1;
|
||||
InterpolatorPointer m_Interpolator2;
|
||||
|
||||
TransformMetaInformation::Pointer
|
||||
m_TransformMetaInfo;
|
||||
|
||||
ChangeInformationFilterPointer
|
||||
m_Filter1,
|
||||
m_Filter2;
|
||||
|
||||
ParametersType m_InitialOptimizerParameters;
|
||||
ParametersType m_LastOptimizerParameters;
|
||||
|
||||
bool m_FixedImageRegionDefined1;
|
||||
bool m_FixedImageRegionDefined2;
|
||||
FixedImageRegionType m_FixedImageRegion1;
|
||||
FixedImageRegionType m_FixedImageRegion2;
|
||||
};
|
||||
|
||||
} // end namespace itk
|
||||
|
||||
#ifndef ITK_MANUAL_INSTANTIATION
|
||||
#include "itkTwoProjectionImageRegistrationMethod.hxx"
|
||||
#endif
|
||||
|
||||
#endif
|
||||
@ -0,0 +1,327 @@
|
||||
/*=========================================================================
|
||||
*
|
||||
* 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 "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_TransformMetaInfo = nullptr; // has to be provided by the user.
|
||||
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());
|
||||
|
||||
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");
|
||||
}
|
||||
|
||||
// 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->SetTransformMetaInfo(m_TransformMetaInfo);
|
||||
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);
|
||||
} else {
|
||||
m_Metric->SetFixedImageRegion1(m_FixedImage1->GetBufferedRegion());
|
||||
}
|
||||
|
||||
if (m_FixedImageRegionDefined2) {
|
||||
m_Metric->SetFixedImageRegion2(m_FixedImageRegion2);
|
||||
} 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);
|
||||
}
|
||||
|
||||
/*
|
||||
* Starts the Registration Process
|
||||
*/
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::StartRegistration()
|
||||
{
|
||||
|
||||
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
|
||||
@ -0,0 +1,300 @@
|
||||
/*=========================================================================
|
||||
*
|
||||
* 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 itkgTwoImageToOneImageMetric_h
|
||||
#define itkgTwoImageToOneImageMetric_h
|
||||
|
||||
#include "DRTMetaInformation.h"
|
||||
#include "itkChangeInformationImageFilter.h"
|
||||
#include "itkDRTHelpers.h"
|
||||
#include "itkGradientRecursiveGaussianImageFilter.h"
|
||||
#include "itkImageBase.h"
|
||||
#include "itkInterpolateImageFunction.h"
|
||||
#include "itkOptimizer.h"
|
||||
#include "itkResampleImageFilter.h"
|
||||
#include "itkSingleValuedCostFunction.h"
|
||||
#include "itkSpatialObject.h"
|
||||
#include "itkTransform.h"
|
||||
#include "itkgSiddonJacobsRayCastInterpolateImageFunction.h"
|
||||
|
||||
#include "itkImageFileWriter.h"
|
||||
#include "itkRescaleIntensityImageFilter.h"
|
||||
|
||||
#include <cstring>
|
||||
#include <iostream>
|
||||
#include <libgen.h>
|
||||
#include <linux/limits.h>
|
||||
#include <unistd.h>
|
||||
|
||||
namespace itk {
|
||||
|
||||
/** \class gTwoImageToOneImageMetric
|
||||
* \brief Computes similarity between two fixed images and one fixed image.
|
||||
*
|
||||
* This Class is templated over the type of the two input images.
|
||||
* It expects a Transform and two Interpolators to be plugged in.
|
||||
* This particular class is the base class for a hierarchy of
|
||||
* similarity metrics.
|
||||
*
|
||||
* This class computes a value that measures the similarity
|
||||
* between two Fixed image and the transformed Moving images.
|
||||
* The Interpolators are used to compute intensity values on
|
||||
* non-grid positions resulting from mapping points through
|
||||
* the Transform.
|
||||
*
|
||||
*
|
||||
* \ingroup RegistrationMetrics
|
||||
* \ingroup TwoProjectionRegistration
|
||||
*
|
||||
*/
|
||||
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
class gTwoImageToOneImageMetric : public SingleValuedCostFunction {
|
||||
public:
|
||||
ITK_DISALLOW_COPY_AND_ASSIGN(gTwoImageToOneImageMetric);
|
||||
|
||||
/** Standard class type alias. */
|
||||
using Self = gTwoImageToOneImageMetric;
|
||||
using Superclass = SingleValuedCostFunction;
|
||||
using Pointer = SmartPointer<Self>;
|
||||
using ConstPointer = SmartPointer<const Self>;
|
||||
|
||||
/** Type used for representing point components */
|
||||
using CoordinateRepresentationType = Superclass::ParametersValueType;
|
||||
|
||||
/** Run-time type information (and related methods). */
|
||||
itkTypeMacro(gTwoImageToOneImageMetric, SingleValuedCostFunction);
|
||||
|
||||
/** Type of the moving Image. */
|
||||
using MovingImageType = TMovingImage;
|
||||
using MovingImagePixelType = typename TMovingImage::PixelType;
|
||||
using MovingImageConstPointer = typename MovingImageType::ConstPointer;
|
||||
|
||||
/** Type of the fixed Image. */
|
||||
using FixedImageType = TFixedImage;
|
||||
using FixedImageConstPointer = typename FixedImageType::ConstPointer;
|
||||
using FixedImageRegionType = typename FixedImageType::RegionType;
|
||||
|
||||
/** Constants for the image dimensions */
|
||||
static constexpr unsigned int MovingImageDimension = TMovingImage::ImageDimension;
|
||||
static constexpr unsigned int FixedImageDimension = TFixedImage::ImageDimension;
|
||||
|
||||
/** Type of the Transform Base class */
|
||||
using TransformType = itk::Euler3DTransform<double>;
|
||||
|
||||
using TransformPointer = typename TransformType::Pointer;
|
||||
using InputPointType = typename TransformType::InputPointType;
|
||||
using OutputPointType = typename TransformType::OutputPointType;
|
||||
using TransformParametersType = typename TransformType::ParametersType;
|
||||
using TransformJacobianType = typename TransformType::JacobianType;
|
||||
|
||||
/** Type of the Interpolator Base class */
|
||||
using InterpolatorType = gSiddonJacobsRayCastInterpolateImageFunction<MovingImageType, CoordinateRepresentationType>;
|
||||
|
||||
/** Gaussian filter to compute the gradient of the Moving Image */
|
||||
using RealType = typename NumericTraits<MovingImagePixelType>::RealType;
|
||||
using GradientPixelType = CovariantVector<RealType, itkGetStaticConstMacro(MovingImageDimension)>;
|
||||
using GradientImageType = Image<GradientPixelType, itkGetStaticConstMacro(MovingImageDimension)>;
|
||||
using GradientImagePointer = SmartPointer<GradientImageType>;
|
||||
using GradientImageFilterType = GradientRecursiveGaussianImageFilter<MovingImageType, GradientImageType>;
|
||||
using GradientImageFilterPointer = typename GradientImageFilterType::Pointer;
|
||||
using InterpolatorPointer = typename InterpolatorType::Pointer;
|
||||
|
||||
using ResampleImageFilterType = ResampleImageFilter<MovingImageType, MovingImageType>;
|
||||
using ResampleImageFilterPointer = typename ResampleImageFilterType::Pointer;
|
||||
using ChangeInformationFilterType = ChangeInformationImageFilter<MovingImageType>;
|
||||
using ChangeInformationFilterPointer = typename ChangeInformationFilterType::Pointer;
|
||||
|
||||
/** Type for the mask of the fixed image. Only pixels that are "inside"
|
||||
this mask will be considered for the computation of the metric */
|
||||
typedef SpatialObject<itkGetStaticConstMacro(FixedImageDimension)> FixedImageMaskType;
|
||||
using FixedImageMaskPointer = typename FixedImageMaskType::Pointer;
|
||||
|
||||
/** Type for the mask of the moving image. Only pixels that are "inside"
|
||||
this mask will be considered for the computation of the metric */
|
||||
typedef SpatialObject<itkGetStaticConstMacro(MovingImageDimension)> MovingImageMaskType;
|
||||
using MovingImageMaskPointer = typename MovingImageMaskType::Pointer;
|
||||
|
||||
/** Type of the measure. */
|
||||
using MeasureType = Superclass::MeasureType;
|
||||
|
||||
/** Type of the derivative. */
|
||||
using DerivativeType = Superclass::DerivativeType;
|
||||
|
||||
/** Type of the Transformation parameters This is the same type used to
|
||||
* represent the search space of the optimization algorithm */
|
||||
using ParametersType = Superclass::ParametersType;
|
||||
|
||||
/** Connect the Fixed Image. */
|
||||
itkSetConstObjectMacro(FixedImage1, FixedImageType);
|
||||
|
||||
/** Connect the Fixed Image. */
|
||||
itkSetConstObjectMacro(FixedImage2, FixedImageType);
|
||||
|
||||
/** Get the Fixed Image. */
|
||||
itkGetConstObjectMacro(FixedImage1, FixedImageType);
|
||||
|
||||
/** Get the Fixed Image. */
|
||||
itkGetConstObjectMacro(FixedImage2, FixedImageType);
|
||||
|
||||
/** Connect the Moving Image. */
|
||||
itkSetConstObjectMacro(MovingImage, MovingImageType);
|
||||
|
||||
/** Get the Moving Image. */
|
||||
itkGetConstObjectMacro(MovingImage, MovingImageType);
|
||||
|
||||
/** Connect the Transform1. */
|
||||
itkSetObjectMacro(Transform1, TransformType);
|
||||
|
||||
/** Get a pointer to the Transform1. */
|
||||
itkGetConstObjectMacro(Transform1, TransformType);
|
||||
|
||||
/** Connect the Transform2. */
|
||||
itkSetObjectMacro(Transform2, TransformType);
|
||||
|
||||
/** Get a pointer to the Transform2. */
|
||||
itkGetConstObjectMacro(Transform2, TransformType);
|
||||
|
||||
/** Connect the Interpolator. */
|
||||
itkSetObjectMacro(Interpolator1, InterpolatorType);
|
||||
|
||||
/** Get a pointer to the Interpolator. */
|
||||
itkGetConstObjectMacro(Interpolator1, InterpolatorType);
|
||||
|
||||
/** Connect the Interpolator. */
|
||||
itkSetObjectMacro(Interpolator2, InterpolatorType);
|
||||
|
||||
/** Get a pointer to the Interpolator. */
|
||||
itkGetConstObjectMacro(Interpolator2, InterpolatorType);
|
||||
|
||||
/** Get the number of pixels considered in the computation. */
|
||||
itkGetConstReferenceMacro(NumberOfPixelsCounted, unsigned long);
|
||||
|
||||
/** Connect the DRTGeometryMetaInfo. */
|
||||
itkSetObjectMacro(TransformMetaInfo, TransformMetaInformation);
|
||||
|
||||
/** Get a pointer to the DRTGeometryMetaInfo. */
|
||||
itkGetConstObjectMacro(TransformMetaInfo, TransformMetaInformation);
|
||||
|
||||
/** Set the region over which the metric will be computed */
|
||||
itkSetMacro(FixedImageRegion1, FixedImageRegionType);
|
||||
|
||||
/** Set the region over which the metric will be computed */
|
||||
itkSetMacro(FixedImageRegion2, FixedImageRegionType);
|
||||
|
||||
/** Get the region over which the metric will be computed */
|
||||
itkGetConstReferenceMacro(FixedImageRegion1, FixedImageRegionType);
|
||||
|
||||
/** Get the region over which the metric will be computed */
|
||||
itkGetConstReferenceMacro(FixedImageRegion2, FixedImageRegionType);
|
||||
|
||||
/** Set/Get the output filters. */
|
||||
itkSetObjectMacro(Filter1, ChangeInformationFilterType);
|
||||
itkGetConstObjectMacro(Filter1, ChangeInformationFilterType);
|
||||
itkSetObjectMacro(Filter2, ChangeInformationFilterType);
|
||||
itkGetConstObjectMacro(Filter2, ChangeInformationFilterType);
|
||||
|
||||
/** Set/Get the moving image mask. */
|
||||
itkSetObjectMacro(MovingImageMask, MovingImageMaskType);
|
||||
itkGetConstObjectMacro(MovingImageMask, MovingImageMaskType);
|
||||
|
||||
/** Set/Get the fixed image mask. */
|
||||
itkSetObjectMacro(FixedImageMask1, FixedImageMaskType);
|
||||
itkSetObjectMacro(FixedImageMask2, FixedImageMaskType);
|
||||
itkGetConstObjectMacro(FixedImageMask1, FixedImageMaskType);
|
||||
itkGetConstObjectMacro(FixedImageMask2, FixedImageMaskType);
|
||||
|
||||
/** Set/Get gradient computation. */
|
||||
itkSetMacro(ComputeGradient, bool);
|
||||
itkGetConstReferenceMacro(ComputeGradient, bool);
|
||||
itkBooleanMacro(ComputeGradient);
|
||||
|
||||
itkSetMacro(MaxTranslation, double);
|
||||
itkGetConstReferenceMacro(MaxTranslation, double);
|
||||
|
||||
/** Get Gradient Image. */
|
||||
itkGetConstObjectMacro(GradientImage, GradientImageType);
|
||||
|
||||
/** Set the parameters defining the Transform. */
|
||||
bool
|
||||
SetTransformParameters(const ParametersType& parameters) const;
|
||||
|
||||
/** Return the number of parameters required by the Transform */
|
||||
|
||||
unsigned int GetNumberOfParameters() const;
|
||||
itk::Optimizer::ScalesType GetWeightings() const;
|
||||
ParametersType GetParameters() const;
|
||||
|
||||
/** Initialize the Metric by making sure that all the components
|
||||
* are present and plugged together correctly */
|
||||
virtual void
|
||||
Initialize();
|
||||
|
||||
protected:
|
||||
gTwoImageToOneImageMetric();
|
||||
~gTwoImageToOneImageMetric() override = default;
|
||||
void
|
||||
PrintSelf(std::ostream& os, Indent indent) const override;
|
||||
|
||||
constexpr static unsigned int Dimension = 3;
|
||||
|
||||
mutable unsigned long m_NumberOfPixelsCounted;
|
||||
mutable double m_bestMeassure;
|
||||
|
||||
FixedImageConstPointer m_FixedImage1;
|
||||
FixedImageConstPointer m_FixedImage2;
|
||||
MovingImageConstPointer m_MovingImage;
|
||||
|
||||
mutable TransformPointer m_Transform1;
|
||||
mutable TransformPointer m_Transform2;
|
||||
InterpolatorPointer m_Interpolator1;
|
||||
InterpolatorPointer m_Interpolator2;
|
||||
|
||||
bool m_ComputeGradient;
|
||||
GradientImagePointer m_GradientImage;
|
||||
|
||||
using OutputPixelType = unsigned char;
|
||||
using OutputImageType = itk::Image<OutputPixelType, Dimension>;
|
||||
|
||||
mutable FixedImageMaskPointer m_FixedImageMask1;
|
||||
mutable FixedImageMaskPointer m_FixedImageMask2;
|
||||
mutable MovingImageMaskPointer m_MovingImageMask;
|
||||
|
||||
TransformMetaInformation::Pointer
|
||||
m_TransformMetaInfo;
|
||||
|
||||
ChangeInformationFilterPointer
|
||||
m_Filter1,
|
||||
m_Filter2;
|
||||
|
||||
double m_MaxTranslation;
|
||||
|
||||
/* Writes the image to the disk. If path is empty it will be saved to /img/out/*/
|
||||
void WriteImage(MovingImageConstPointer img, std::string fileName, std::string path) const;
|
||||
|
||||
private:
|
||||
FixedImageRegionType m_FixedImageRegion1;
|
||||
FixedImageRegionType m_FixedImageRegion2;
|
||||
};
|
||||
|
||||
} // end namespace itk
|
||||
|
||||
#ifndef ITK_MANUAL_INSTANTIATION
|
||||
#include "itkgTwoImageToOneImageMetric.hxx"
|
||||
#endif
|
||||
|
||||
#endif
|
||||
@ -0,0 +1,395 @@
|
||||
/*=========================================================================
|
||||
*
|
||||
* 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 itkgTwoImageToOneImageMetric_hxx
|
||||
#define itkgTwoImageToOneImageMetric_hxx
|
||||
|
||||
#include "itkgTwoImageToOneImageMetric.h"
|
||||
#include <limits>
|
||||
|
||||
namespace itk {
|
||||
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::gTwoImageToOneImageMetric()
|
||||
{
|
||||
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_TransformMetaInfo = nullptr; // has to be provided by the user.
|
||||
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;
|
||||
m_ComputeGradient = true; // metric computes gradient by default
|
||||
m_NumberOfPixelsCounted = 0; // initialize to zero
|
||||
m_bestMeassure = std::numeric_limits<short>::max();
|
||||
m_MaxTranslation = std::numeric_limits<short>::max();
|
||||
}
|
||||
|
||||
/*
|
||||
* Set the parameters that define a unique transform
|
||||
*/
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
bool gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::SetTransformParameters(const ParametersType& parameters) const
|
||||
{
|
||||
|
||||
if (!m_Transform1) {
|
||||
itkExceptionMacro(<< "Transform 1 has not been assigned");
|
||||
}
|
||||
if (!m_Transform2) {
|
||||
itkExceptionMacro(<< "Transform 2 has not been assigned");
|
||||
}
|
||||
|
||||
auto transformParameters = m_Transform1->GetParameters();
|
||||
double TranslationAlongX = transformParameters[3];
|
||||
double TranslationAlongY = transformParameters[4];
|
||||
double TranslationAlongZ = transformParameters[5];
|
||||
double RotationAlongX = transformParameters[0];
|
||||
double RotationAlongY = transformParameters[1];
|
||||
double RotationAlongZ = transformParameters[2];
|
||||
|
||||
switch (m_TransformMetaInfo->GetDegreeOfFreedom()) {
|
||||
case X_ONLY:
|
||||
TranslationAlongX = parameters[0];
|
||||
break;
|
||||
case Y_ONLY:
|
||||
TranslationAlongY = parameters[0];
|
||||
break;
|
||||
case Z_ONLY:
|
||||
TranslationAlongZ = parameters[0];
|
||||
break;
|
||||
case THREE_DOF:
|
||||
TranslationAlongX = parameters[0];
|
||||
TranslationAlongY = parameters[1];
|
||||
TranslationAlongZ = parameters[2];
|
||||
break;
|
||||
case ROTATION_ONLY:
|
||||
RotationAlongX = parameters[0];
|
||||
RotationAlongY = parameters[1];
|
||||
RotationAlongZ = parameters[2];
|
||||
break;
|
||||
case SIX_DOF:
|
||||
TranslationAlongX = parameters[0];
|
||||
TranslationAlongY = parameters[1];
|
||||
TranslationAlongZ = parameters[2];
|
||||
RotationAlongX = parameters[3];
|
||||
RotationAlongY = parameters[4];
|
||||
RotationAlongZ = parameters[5];
|
||||
break;
|
||||
|
||||
default:
|
||||
itkExceptionMacro(<< "Unsupported degree of freedeom");
|
||||
break;
|
||||
}
|
||||
|
||||
std::cout << "New Transform Parameters = " << 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 << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++" << 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);
|
||||
|
||||
if (transformValid) {
|
||||
m_Transform1->SetParameters(transformParameters);
|
||||
m_Transform2->SetParameters(transformParameters);
|
||||
} else {
|
||||
std::cout << " Transform invalid, out of range!" << std::endl;
|
||||
std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
|
||||
}
|
||||
return transformValid;
|
||||
}
|
||||
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
void gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::Initialize()
|
||||
{
|
||||
|
||||
if (!m_Transform1) {
|
||||
itkExceptionMacro(<< "Transform is not present");
|
||||
}
|
||||
|
||||
if (!m_Transform1) {
|
||||
itkExceptionMacro(<< "Transform is not present");
|
||||
}
|
||||
|
||||
if (!m_Transform2) {
|
||||
itkExceptionMacro(<< "Transform is not present");
|
||||
}
|
||||
|
||||
if (!m_Interpolator1) {
|
||||
itkExceptionMacro(<< "Interpolator1 is not present");
|
||||
}
|
||||
if (!m_Interpolator2) {
|
||||
itkExceptionMacro(<< "Interpolator2 is not present");
|
||||
}
|
||||
|
||||
if (!m_FixedImage1) {
|
||||
itkExceptionMacro(<< "FixedImage1 is not present");
|
||||
}
|
||||
|
||||
if (!m_FixedImage2) {
|
||||
itkExceptionMacro(<< "FixedImage2 is not present");
|
||||
}
|
||||
|
||||
if (m_FixedImageRegion1.GetNumberOfPixels() == 0) {
|
||||
itkExceptionMacro(<< "FixedImageRegion1 is empty");
|
||||
}
|
||||
|
||||
if (m_FixedImageRegion2.GetNumberOfPixels() == 0) {
|
||||
itkExceptionMacro(<< "FixedImageRegion2 is empty");
|
||||
}
|
||||
|
||||
// If the image is provided by a source, update the source.
|
||||
if (m_MovingImage->GetSource()) {
|
||||
m_MovingImage->GetSource()->Update();
|
||||
}
|
||||
|
||||
// If the image is provided by a source, update the source.
|
||||
if (m_FixedImage1->GetSource()) {
|
||||
m_FixedImage1->GetSource()->Update();
|
||||
}
|
||||
|
||||
if (m_FixedImage2->GetSource()) {
|
||||
m_FixedImage2->GetSource()->Update();
|
||||
}
|
||||
|
||||
// Make sure the FixedImageRegion is within the FixedImage buffered region
|
||||
if (!m_FixedImageRegion1.Crop(m_FixedImage1->GetBufferedRegion())) {
|
||||
itkExceptionMacro(<< "FixedImageRegion1 does not overlap the fixed image buffered region");
|
||||
}
|
||||
|
||||
if (!m_FixedImageRegion2.Crop(m_FixedImage2->GetBufferedRegion())) {
|
||||
itkExceptionMacro(<< "FixedImageRegion2 does not overlap the fixed image buffered region");
|
||||
}
|
||||
|
||||
this->SetTransformParameters(this->GetParameters());
|
||||
|
||||
if (m_ComputeGradient) {
|
||||
|
||||
GradientImageFilterPointer gradientFilter = GradientImageFilterType::New();
|
||||
|
||||
gradientFilter->SetInput(m_MovingImage);
|
||||
|
||||
const typename MovingImageType::SpacingType& spacing = m_MovingImage->GetSpacing();
|
||||
double maximumSpacing = 0.0;
|
||||
for (unsigned int i = 0; i < MovingImageDimension; i++) {
|
||||
if (spacing[i] > maximumSpacing) {
|
||||
maximumSpacing = spacing[i];
|
||||
}
|
||||
}
|
||||
gradientFilter->SetSigma(maximumSpacing);
|
||||
gradientFilter->SetNormalizeAcrossScale(true);
|
||||
|
||||
gradientFilter->Update();
|
||||
|
||||
m_GradientImage = gradientFilter->GetOutput();
|
||||
}
|
||||
|
||||
// If there are any observers on the metric, call them to give the
|
||||
// user code a chance to set parameters on the metric
|
||||
this->InvokeEvent(InitializeEvent());
|
||||
}
|
||||
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
unsigned int gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetNumberOfParameters() const
|
||||
{
|
||||
if (!m_TransformMetaInfo) {
|
||||
itkExceptionMacro(<< "TransformMetaInfo is missing");
|
||||
}
|
||||
|
||||
switch (m_TransformMetaInfo->GetDegreeOfFreedom()) {
|
||||
case X_ONLY:
|
||||
case Y_ONLY:
|
||||
case Z_ONLY:
|
||||
return 1;
|
||||
case THREE_DOF:
|
||||
case ROTATION_ONLY:
|
||||
return 3;
|
||||
case SIX_DOF:
|
||||
return 6;
|
||||
default:
|
||||
itkExceptionMacro(<< "Unsupported degree of freedeom");
|
||||
break;
|
||||
}
|
||||
|
||||
return m_Transform1->GetNumberOfParameters();
|
||||
}
|
||||
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
itk::Optimizer::ScalesType gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetWeightings() const
|
||||
{
|
||||
itk::Optimizer::ScalesType weightings(this->GetNumberOfParameters());
|
||||
|
||||
switch (m_TransformMetaInfo->GetDegreeOfFreedom()) {
|
||||
case X_ONLY:
|
||||
case Y_ONLY:
|
||||
case Z_ONLY:
|
||||
weightings[0] = 1.0;
|
||||
break;
|
||||
;
|
||||
case THREE_DOF:
|
||||
weightings[0] = 1.;
|
||||
weightings[1] = 1.;
|
||||
weightings[2] = 1.;
|
||||
break;
|
||||
case ROTATION_ONLY:
|
||||
weightings[0] = 1. / dtr;
|
||||
weightings[1] = 1. / dtr;
|
||||
weightings[2] = 1. / dtr;
|
||||
break;
|
||||
case SIX_DOF:
|
||||
weightings[0] = 1.;
|
||||
weightings[1] = 1.;
|
||||
weightings[2] = 1.;
|
||||
weightings[3] = 1. / dtr;
|
||||
weightings[4] = 1. / dtr;
|
||||
weightings[5] = 1. / dtr;
|
||||
break;
|
||||
|
||||
default:
|
||||
itkExceptionMacro(<< "Unsupported degree of freedeom");
|
||||
break;
|
||||
}
|
||||
return weightings;
|
||||
}
|
||||
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
typename gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::ParametersType gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetParameters() const
|
||||
{
|
||||
|
||||
ParametersType parametersTransform = m_Transform1->GetParameters(); //angleX, angleY, angleZ, transX, transY, transZ
|
||||
ParametersType parameters(this->GetNumberOfParameters());
|
||||
switch (m_TransformMetaInfo->GetDegreeOfFreedom()) {
|
||||
case X_ONLY:
|
||||
parameters[0] = parametersTransform[3];
|
||||
break;
|
||||
case Y_ONLY:
|
||||
parameters[0] = parametersTransform[4];
|
||||
break;
|
||||
case Z_ONLY:
|
||||
parameters[0] = parametersTransform[5];
|
||||
break;
|
||||
case THREE_DOF:
|
||||
parameters[0] = parametersTransform[3];
|
||||
parameters[1] = parametersTransform[4];
|
||||
parameters[2] = parametersTransform[5];
|
||||
break;
|
||||
case ROTATION_ONLY:
|
||||
parameters[0] = parametersTransform[0];
|
||||
parameters[1] = parametersTransform[1];
|
||||
parameters[2] = parametersTransform[2];
|
||||
break;
|
||||
case SIX_DOF:
|
||||
parameters[0] = parametersTransform[3];
|
||||
parameters[1] = parametersTransform[4];
|
||||
parameters[2] = parametersTransform[5];
|
||||
parameters[3] = parametersTransform[0];
|
||||
parameters[4] = parametersTransform[1];
|
||||
parameters[5] = parametersTransform[2];
|
||||
break;
|
||||
|
||||
default:
|
||||
itkExceptionMacro(<< "Unsupported degree of freedeom");
|
||||
break;
|
||||
}
|
||||
return parameters;
|
||||
}
|
||||
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
void gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::PrintSelf(std::ostream& os, Indent indent) const
|
||||
{
|
||||
Superclass::PrintSelf(os, indent);
|
||||
os << indent << "ComputeGradient: " << static_cast<typename NumericTraits<bool>::PrintType>(m_ComputeGradient)
|
||||
<< std::endl;
|
||||
os << indent << "Moving Image: " << m_MovingImage.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 << "Gradient Image: " << m_GradientImage.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 << "Filter 1: " << m_Filter1.GetPointer() << std::endl;
|
||||
os << indent << "Filter 2: " << m_Filter2.GetPointer() << std::endl;
|
||||
os << indent << "Transform MetaInfoe: " << m_TransformMetaInfo.GetPointer() << std::endl;
|
||||
os << indent << "FixedImageRegion 1: " << m_FixedImageRegion1 << std::endl;
|
||||
os << indent << "FixedImageRegion 2: " << m_FixedImageRegion2 << std::endl;
|
||||
os << indent << "Moving Image Mask: " << m_MovingImageMask.GetPointer() << std::endl;
|
||||
os << indent << "Fixed Image Mask 1: " << m_FixedImageMask1.GetPointer() << std::endl;
|
||||
os << indent << "Fixed Image Mask 2: " << m_FixedImageMask2.GetPointer() << std::endl;
|
||||
os << indent << "Number of Pixels Counted: " << m_NumberOfPixelsCounted << std::endl;
|
||||
}
|
||||
|
||||
template <typename TFixedImage, typename TMovingImage>
|
||||
void gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::WriteImage(MovingImageConstPointer img, std::string fileName, std::string path) const
|
||||
{
|
||||
|
||||
using RescaleFilterType = itk::RescaleIntensityImageFilter<MovingImageType, OutputImageType>;
|
||||
using WriterType = itk::ImageFileWriter<OutputImageType>;
|
||||
|
||||
if (img) {
|
||||
|
||||
// Rescale the intensity of the projection images to 0-255 for output.
|
||||
typename RescaleFilterType::Pointer rescaler1 = RescaleFilterType::New();
|
||||
rescaler1->SetOutputMinimum(0);
|
||||
rescaler1->SetOutputMaximum(255);
|
||||
rescaler1->SetInput(img);
|
||||
|
||||
WriterType::Pointer writer1 = WriterType::New();
|
||||
|
||||
if (fileName.empty()) {
|
||||
|
||||
char result[PATH_MAX];
|
||||
ssize_t count = readlink("/proc/self/exe", result, PATH_MAX);
|
||||
if (count != -1) {
|
||||
path = dirname(result);
|
||||
path += "/img/out/";
|
||||
}
|
||||
}
|
||||
|
||||
if (fileName.empty()) {
|
||||
|
||||
fileName = path;
|
||||
fileName += "2D1.tif";
|
||||
}
|
||||
writer1->SetFileName(fileName);
|
||||
writer1->SetInput(rescaler1->GetOutput());
|
||||
|
||||
try {
|
||||
std::cout << "Writing image: " << fileName << std::endl;
|
||||
writer1->Update();
|
||||
} catch (itk::ExceptionObject& err) {
|
||||
std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
|
||||
std::cerr << err << std::endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
} // end namespace itk
|
||||
|
||||
#endif
|
||||
Reference in New Issue
Block a user