Support multiple metrics and regularizations, add tests and backward compatibility ("legacy hack")

[N-Dekker]

• Cherry-picked pull-request ENH: Add support for multiple metrics and regularization #1437 from @Leirbag-gabrieL (followed by clang-format)
• Originally based on pull request ENH: Add support for multiple metrics and regularization #482 by @alvarez-pa
• Addresssing issue Multiple metrics and regularization #481 by @alvarez-pa
• Included backward compatibility ("legacy hack" support), following comment ENH: Add support for multiple metrics and regularization #1437 (comment) by @stefanklein
• Added unit tests

[N-Dekker]

@mstaring @stefanklein This unit test is intended to check that the original "legacy hack" still works. Can you please check that this is indeed the way the "legacy hack" worked? I mean, was the trick to add a dummy fixed image, a dummy moving image, an extra "FixedImagePyramid", an extra "MovingImagePyramid", an extra "ImageSampler", and an extra "Interpolator"?

[N-Dekker]

@stefanklein @mstaring Now I'm starting to doubt if a transform based metric like "TransformBendingEnergyPenalty" can work without having its own images. TransformBendingEnergyPenaltyTerm is derived from TransformPenaltyTerm<TFixedImage, TScalarType>, which is derived from AdvancedImageToImageMetric<TFixedImage, TFixedImage>. AdvancedImageToImageMetric has a pointer to a fixed image, and a pointer to a moving image.

The "GetValue" member functions of TransformBendingEnergyPenaltyTerm (for example TransformBendingEnergyPenaltyTerm::ThreadedGetValueAndDerivative) iterate over image samples.

---

Looks like the logic is implemented in elastix::MultiMetricMultiResolutionRegistration::SetComponents(), here:

elastix/Components/Registrations/MultiMetricMultiResolutionRegistration/elxMultiMetricMultiResolutionRegistration.hxx

Lines 311 to 348 in 2aafd7b

	/** Samplers are not always needed: */
	for (unsigned int i = 0; i < nrOfMetrics; ++i)
	{
	if (this->GetElastix()->GetElxMetricBase(i)->GetAdvancedMetricUseImageSampler())
	{
	/** Try the i-th sampler for the i-th metric. */
	if (this->GetElastix()->GetElxImageSamplerBase(i))
	{
	this->GetElastix()->GetElxMetricBase(i)->SetAdvancedMetricImageSampler(
	this->GetElastix()->GetElxImageSamplerBase(i)->GetAsITKBaseType());
	}
	else
	{
	/** When a different fixed image pyramid is used for each metric,
	* using one sampler for all metrics makes no sense.
	*/
	if (this->GetElastix()->GetElxFixedImagePyramidBase(i))
	{
	log::error(std::ostringstream{} << "ERROR: An ImageSamper for metric " << i << " must be provided!");
	itkExceptionMacro("Not enough ImageSamplers provided!\nProvide an ImageSampler for metric "
	<< i << ", like:\n (ImageSampler \"Random\" ... \"Random\")");
	}
	/** Try the zeroth image sampler for each metric. */
	if (this->GetElastix()->GetElxImageSamplerBase(0))
	{
	this->GetElastix()->GetElxMetricBase(i)->SetAdvancedMetricImageSampler(
	this->GetElastix()->GetElxImageSamplerBase(0)->GetAsITKBaseType());
	}
	else
	{
	log::error("ERROR: No ImageSampler has been specified.");
	itkExceptionMacro("One of the metrics requires an ImageSampler, but it is not available!");
	}
	}
	} // if sampler required by metric
	} // for loop over metrics

So when the i-th metric requires to use an image sampler (including "TransformBendingEnergyPenalty") and there is no "i-th sampler", the metric may use "the zeroth image sampler".

So then, I guess this pull request should be OK, right?

[N-Dekker]

elastix supports the following 23 metrics. I guess we should check which metrics do or don't use the input images. Because it appears that transform based metrics (like "TransformBendingEnergyPenalty") may still use the images. I think CheckOnInitialize should become less strict than the current (main/release) version, but I doubt if that should be based on whether or not metrics are transform based.

"AdvancedKappaStatistic"    
"AdvancedMattesMutualInformation"    
"AdvancedMeanSquares"    
"AdvancedNormalizedCorrelation"    
"TransformBendingEnergyPenalty"    
"CorrespondingPointsEuclideanDistanceMetric"    
"DisplacementMagnitudePenalty"    
"DistancePreservingRigidityPenalty"    
"GradientDifference"    
"Impact"    
"KNNGraphAlphaMutualInformation"    
"MissingStructurePenalty"    
"NormalizedGradientCorrelation"    
"NormalizedMutualInformation"    
"PatternIntensity"    
"PCAMetric"    
"PCAMetric2"    
"PolydataDummyPenalty"    
"TransformRigidityPenalty"    
"StatisticalShapePenalty"    
"SumOfPairwiseCorrelationCoefficientsMetric"    
"SumSquaredTissueVolumeDifference"    
"VarianceOverLastDimensionMetric"

Four metric classes are derived from TransformPenaltyTerm (which is itself derived from AdvancedImageToImageMetric):

TransformBendingEnergyPenaltyTerm
DisplacementMagnitudePenaltyTerm
DistancePreservingRigidityPenaltyTerm
TransformRigidityPenaltyTerm

I think the only thing that matters is whether elastix::MetricBase::GetAdvancedMetricUseImageSampler() returns true or false. Metric classes at the following 13 locations do SetUseImageSampler(true):

Common\CostFunctions\itkParzenWindowHistogramImageToImageMetric
Components\Metrics\AdvancedKappaStatistic\itkAdvancedKappaStatisticImageToImageMetric
Components\Metrics\AdvancedMeanSquares\itkAdvancedMeanSquaresImageToImageMetric
Components\Metrics\AdvancedNormalizedCorrelation\itkAdvancedNormalizedCorrelationImageToImageMetric
Components\Metrics\BendingEnergyPenalty\itkTransformBendingEnergyPenaltyTerm
Components\Metrics\DisplacementMagnitudePenalty\itkDisplacementMagnitudePenaltyTerm
Components\Metrics\Impact\itkImpactImageToImageMetric
Components\Metrics\KNNGraphAlphaMutualInformation\itkKNNGraphAlphaMutualInformationImageToImageMetric
Components\Metrics\PCAMetric2\itkPCAMetric2
Components\Metrics\PCAMetric\itkPCAMetric
Components\Metrics\SumOfPairwiseCorrelationsMetric\itkSumOfPairwiseCorrelationCoefficientsMetric
Components\Metrics\SumSquaredTissueVolumeDifferenceMetric\itkSumSquaredTissueVolumeDifferenceImageToImageMetric
Components\Metrics\VarianceOverLastDimension\itkVarianceOverLastDimensionImageMetric