240 lines
9.1 KiB
C++
240 lines
9.1 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 itkMutualInformationTwoImageToOneImageMetric_hxx
|
|
#define itkMutualInformationTwoImageToOneImageMetric_hxx
|
|
|
|
#include "itkImageRegionConstIteratorWithIndex.h"
|
|
#include "itkMattesMutualInformationImageToImageMetric.h"
|
|
#include "itkMutualInformationTwoImageToOneImageMetric.h"
|
|
#include "itkNearestNeighborInterpolateImageFunction.h"
|
|
#include "itkNormalizeImageFilter.h"
|
|
#include "itkTranslationTransform.h"
|
|
|
|
#include "itkProgressReporter.h"
|
|
|
|
#include <limits>
|
|
|
|
namespace itk {
|
|
|
|
template <typename TFixedImage, typename TMovingImage>
|
|
MutualInformationTwoImageToOneImageMetric<TFixedImage,
|
|
TMovingImage>::MutualInformationTwoImageToOneImageMetric()
|
|
{
|
|
m_SubtractMean = false;
|
|
m_NumberOfHistogramBins = 50;
|
|
}
|
|
|
|
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");
|
|
}
|
|
|
|
if (!this->m_Transform1) {
|
|
itkExceptionMacro(<< "m_Transform1 has not been assigned");
|
|
}
|
|
|
|
if (!this->m_Transform2) {
|
|
itkExceptionMacro(<< "m_Transform2 has not been assigned");
|
|
}
|
|
|
|
|
|
// std::cout<< "----- itkMutualInformationTwoImage GetValue -----" <<std::endl;
|
|
// std::cout<< "Mutual GetValue T1 :: center : " <<
|
|
// this->m_Transform1.GetPointer()->GetCenter() <<std::endl;
|
|
// std::cout<< "Mutual GetValue T1 :: Pars : " <<
|
|
// this->m_Transform1.GetPointer()->GetParameters() <<std::endl;
|
|
// std::cout<< "--------------------------------" <<std::endl;
|
|
|
|
this->m_Interpolator1->SetTransform(this->m_Transform1);
|
|
this->m_Interpolator1->Initialize();
|
|
this->m_Filter1->Update();
|
|
|
|
using InternalMetricType = itk::MattesMutualInformationImageToImageMetric<TFixedImage, TMovingImage>;
|
|
/*
|
|
* NearestNeighborType in our images results in not resolving for the cranio-caudal direction.
|
|
* in LPS-Z we are very much bound to the CT resolution, therefore if we take the Nearest Neighbor
|
|
* we simply do not explore that direction... and will always have an uncertainty of 1 slice thickness.
|
|
* */
|
|
//using NearestNeighborType = itk::NearestNeighborInterpolateImageFunction<TMovingImage, double>;
|
|
/*
|
|
* LinearInterpolateImageFunction performs better for us, performance degradation not noticeable.
|
|
* */
|
|
using LinearInterpType = itk::LinearInterpolateImageFunction<TMovingImage, double>;
|
|
|
|
//We need to set transform and parameters.
|
|
auto dummyTransform = TransformType::New();
|
|
auto dummyParameters = dummyTransform->GetParameters();
|
|
auto dummyInterpolator1 = LinearInterpType::New();
|
|
|
|
auto metric1 = InternalMetricType::New();
|
|
metric1->SetTransform(dummyTransform);
|
|
metric1->SetTransformParameters(dummyParameters);
|
|
metric1->SetInterpolator(dummyInterpolator1);
|
|
|
|
//auto fixedImageRegion1 = fixedImage1->GetBufferedRegion();
|
|
//metric1->SetFixedImageRegion(fixedImageRegion1);
|
|
/* We get the fixed image region from the parent class... */
|
|
metric1->SetFixedImageRegion(Superclass::GetFixedImageRegion1());
|
|
//std::cout<< "-----> Mutual :: fixedImageRegion1: "<< metric1->GetFixedImageRegion() << std::endl;
|
|
|
|
auto movingImageRegion = this->m_Filter1->GetOutput()->GetBufferedRegion();
|
|
//unsigned int numberOfPixels = movingImageRegion.GetNumberOfPixels();
|
|
|
|
//auto numberOfSamples = static_cast<unsigned int>(numberOfPixels * 0.50); //100%
|
|
// since we have a ROI, then we should not set allPixels to TRUE.
|
|
//metric1->UseAllPixelsOn();
|
|
//std::cout << "m_NumberOfHistogramBins " << m_NumberOfHistogramBins << std::endl;
|
|
metric1->SetNumberOfHistogramBins(m_NumberOfHistogramBins);
|
|
|
|
// InternalImageType::IndexType pIndex;
|
|
// pIndex[0]=200;
|
|
// pIndex[1]=300;
|
|
// InternalImageType::PixelType pPixel = fixedImage1->GetPixel(pIndex);
|
|
// std::cout<<"MATTES PIXEL: "<<pPixel <<std::endl;
|
|
|
|
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 = LinearInterpType::New();
|
|
auto metric2 = InternalMetricType::New();
|
|
metric2->SetTransform(dummyTransform);
|
|
metric2->SetTransformParameters(dummyParameters);
|
|
metric2->SetInterpolator(dummyInterpolator2);
|
|
|
|
// auto fixedImageRegion2 = fixedImage2->GetBufferedRegion();
|
|
// metric2->SetFixedImageRegion(fixedImageRegion2);
|
|
metric2->SetFixedImageRegion(Superclass::GetFixedImageRegion2());
|
|
//std::cout<< "-----> Mutual :: fixedImageRegion2: "<< metric2->GetFixedImageRegion() << std::endl;
|
|
|
|
movingImageRegion = this->m_Filter2->GetOutput()->GetBufferedRegion();
|
|
//numberOfPixels = movingImageRegion.GetNumberOfPixels();
|
|
|
|
//numberOfSamples = static_cast<unsigned int>(numberOfPixels * 0.50); //100%
|
|
//metric2->SetNumberOfSpatialSamples(numberOfSamples);
|
|
//metric2->UseAllPixelsOn();
|
|
metric2->SetNumberOfHistogramBins(m_NumberOfHistogramBins);
|
|
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
|