477 lines
17 KiB
C++
477 lines
17 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 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_internalMeta1 = nullptr;
|
|
m_internalMeta2 = nullptr;
|
|
m_Interpolator1 = nullptr; // has to be provided by the user.
|
|
m_Interpolator2 = nullptr; // has to be provided by the user.
|
|
m_GradientImage = nullptr; // will receive the output of the filter;
|
|
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");
|
|
}
|
|
if (!m_internalMeta1){
|
|
itkExceptionMacro(<< "m_internalMeta1 has not been assigned");
|
|
}
|
|
if (!m_internalMeta2){
|
|
itkExceptionMacro(<< "m_internalMeta2 has not been assigned");
|
|
}
|
|
if (!m_IsocTransf){
|
|
itkExceptionMacro(<< "m_IsocTransf has not been assigned");
|
|
}
|
|
|
|
//auto transformParameters = m_Transform1->GetParameters();
|
|
auto transformParameters = m_IsocTransf->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];
|
|
|
|
//std::cout << "m_IsocTransf PARS: "<< transformParameters <<std::endl;
|
|
|
|
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 ISOC IEC = " << std::endl;
|
|
// std::cout << " Translation X = " << TranslationAlongX << " mm" << std::endl;
|
|
// std::cout << " Translation Y = " << TranslationAlongY << " mm" << std::endl;
|
|
// std::cout << " Translation Z = " << TranslationAlongZ << " mm" << std::endl;
|
|
// std::cout << " Rotation Along X = " << RotationAlongX /*/ dtr*/ << " deg" << std::endl;
|
|
// std::cout << " Rotation Along Y = " << RotationAlongY /*/ dtr*/ << " deg" << std::endl;
|
|
// std::cout << " Rotation Along Z = " << RotationAlongZ /*/ dtr*/ << " deg" << std::endl;
|
|
// std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
|
|
|
|
|
|
// individual R and T components of the IsoC transform
|
|
ImageType3D::PointType
|
|
pT;
|
|
pT[0] = TranslationAlongX;
|
|
pT[1] = TranslationAlongY;
|
|
pT[2] = TranslationAlongZ;
|
|
|
|
ImageType3D::PointType
|
|
pR;
|
|
pR[0] = RotationAlongX;
|
|
pR[1] = RotationAlongY;
|
|
pR[2] = RotationAlongZ;
|
|
|
|
// check if we are within tolerance
|
|
bool transformValid =
|
|
(std::abs(pT[0]) < m_MaxTranslation) &&
|
|
(std::abs(pT[1]) < m_MaxTranslation) &&
|
|
(std::abs(pT[2]) < m_MaxTranslation);
|
|
|
|
if (transformValid) {
|
|
/*
|
|
Here we calculate the transform for eack projection corresponding to the isocentric
|
|
one we are optimizing.
|
|
This calculation takes into account calibrated center of projection for each view.
|
|
*/
|
|
|
|
// Calculate internal transform1
|
|
TransformType::Pointer CurrTransform;
|
|
CurrTransform = CalculateInternalTransformV3(
|
|
pT,
|
|
pR,
|
|
m_internalMeta1->GetpCalProjCenter(),
|
|
m_internalMeta1->GetpRTIsocenter(),
|
|
m_internalMeta1->GetIECtoLPSDirs());
|
|
|
|
m_Transform1->SetIdentity();
|
|
m_Transform1->SetParameters(
|
|
CurrTransform->GetParameters());
|
|
m_Transform1->SetCenter(
|
|
m_internalMeta1->GetpCalProjCenter());
|
|
|
|
// std::cout<<"----- itkgTwoImageToOne SetTransformParameters -----"<<std::endl;
|
|
// std::cout << "m_Transform1 :: Center: "<< m_Transform1->GetCenter()<<std::endl;
|
|
// std::cout << "m_Transform1 :: Pars: "<< m_Transform1->GetParameters()<<std::endl;
|
|
// std::cout<< "--------------------------------" <<std::endl;
|
|
|
|
|
|
// Calculate internal transform2
|
|
CurrTransform = CalculateInternalTransformV3(
|
|
pT,
|
|
pR,
|
|
m_internalMeta2->GetpCalProjCenter(),
|
|
m_internalMeta2->GetpRTIsocenter(),
|
|
m_internalMeta2->GetIECtoLPSDirs());
|
|
|
|
m_Transform2->SetIdentity();
|
|
m_Transform2->SetParameters(
|
|
CurrTransform->GetParameters());
|
|
m_Transform2->SetCenter(
|
|
m_internalMeta2->GetpCalProjCenter());
|
|
|
|
} else {
|
|
std::cout << " Transform invalid, out of range!" << std::endl;
|
|
std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
|
|
}
|
|
|
|
return transformValid;
|
|
}
|
|
|
|
template <typename TFixedImage, typename TMovingImage>
|
|
void gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::Initialize()
|
|
{
|
|
|
|
if (!m_IsocTransf) {
|
|
itkExceptionMacro(<< "Transform is not present");
|
|
}
|
|
|
|
if (!m_Transform1) {
|
|
itkExceptionMacro(<< "Transform1 is not present");
|
|
}
|
|
|
|
if (!m_Transform2) {
|
|
itkExceptionMacro(<< "Transform2 is not present");
|
|
}
|
|
|
|
if(!m_internalMeta1) {
|
|
itkExceptionMacro(<< "m_internalMeta1 is not present");
|
|
}
|
|
if(!m_internalMeta2) {
|
|
itkExceptionMacro(<< "m_internalMeta2 is not present");
|
|
}
|
|
|
|
if (!m_Interpolator1) {
|
|
itkExceptionMacro(<< "Interpolator1 is not present");
|
|
}
|
|
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();
|
|
}
|
|
|
|
//std::cout<<"FixedImageRegion1 " <<this->GetFixedImageRegion1() <<std::endl;
|
|
//std::cout<<"m_FixedImage1 buffered " <<this->m_FixedImage1->GetBufferedRegion() <<std::endl;
|
|
//std::cout<<"Moving buffered " <<this->m_Filter1->GetOutput()->GetBufferedRegion() <<std::endl;
|
|
//std::cout<<"GetFixedImage1 buffered " <<this->GetFixedImage1()->GetBufferedRegion() <<std::endl;
|
|
//std::cout<<"Filter1 buffered " <<this->m_Filter1->GetOutput()->GetBufferedRegion() <<std::endl;
|
|
|
|
// 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_IsocTransf->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:
|
|
//G: ITHINKTHOSEARE FLIPPED!!
|
|
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 parametersTransform = m_IsocTransf->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
|