ENH: antsAtroposN4.sh - exposed more subcommand parameters as arguments

[br-cpvc]

when trying to optimize antsAtroposN4.sh execution speed contra output precision it would be great to expose these internal variables as external cli parameters

[cookpa]

Thanks for this.

To do:

• Verify script is compatible with existing antsCorticalThickness.sh call (same options get used)
• Verify non-default options get passed through correctly (eg, Atropos -g 1 works and doesn't have to be -g [ 1 ])

[gdevenyi]

when trying to optimize antsAtroposN4.sh execution speed contra output precision

Do you have any docs/experiments on why these Atropos options speed things up/what they do? Atropos is a bit of black box as it has a ton of things not mentioned in the paper which introduces it.

[br-cpvc]

@gdevenyi That is a very good question :) I would also very much appreciate more information/documentation on what the different Atropos parameters actually affect. Unfortunately I do not have the enough insights or knowledge to fully explain it, but I can provide what insight I have gained during my work.

We have been testing out different combinations of parameters and have currently found that the following parameters yield the best trade off between execution speed and output quality for our input images: Atropos parameters "-e 1 -g 0" (for the N4BiasFieldCorrection we use the parameters "-s 4 -c [35x35x35x35,0] -b [200]").

In regards to why the "-g 0" Atropos paramter accelerate the execution speed: This sets the m_UseAsynchronousUpdating=false in the antsAtroposSegmentationImageFilter.hxx filter used by Atropos, which skips the computations inside the if statement in:

ANTs/ImageSegmentation/antsAtroposSegmentationImageFilter.hxx

Lines 1153 to 1206 in 77b8f2a

	if (this->m_UseAsynchronousUpdating)
	{
	if (this->m_MaximumICMCode == 0)
	{
	this->ComputeICMCodeImage();
	}
	typename NeighborhoodIterator<ClassifiedImageType>::RadiusType radius;
	for (unsigned int d = 0; d < ImageDimension; d++)
	{
	radius[d] = this->m_MRFRadius[d];
	}
	NeighborhoodIterator<ClassifiedImageType> ItO(radius, this->GetOutput(), this->GetOutput()->GetRequestedRegion());
	maxPosteriorSum = 0.0;
	RealType oldMaxPosteriorSum = -1.0;
	unsigned int numberOfIterations = 0;
	while (maxPosteriorSum > oldMaxPosteriorSum && numberOfIterations++ < this->m_MaximumNumberOfICMIterations)
	{
	itkDebugMacro("ICM iteration: " << numberOfIterations);
	oldMaxPosteriorSum = maxPosteriorSum;
	maxPosteriorSum = 0.0;
	// Iterate randomly through the ICM codes so we don't visit the same codes
	// in the same order every iteration. We use the Fisher-Yates shuffle to
	// create a random shuffling between 1 and m_MaximumICMCode.
	Array<LabelType> icmCodeSet(this->m_MaximumICMCode);
	for (unsigned int n = 1; n <= this->m_MaximumICMCode; n++)
	{
	icmCodeSet[n - 1] = static_cast<LabelType>(n);
	}
	for (int i = icmCodeSet.Size() - 1; i > 0; i--)
	{
	unsigned int j = this->m_Randomizer->GetIntegerVariate(i);
	LabelType tmp = icmCodeSet[i];
	icmCodeSet[i] = icmCodeSet[j];
	icmCodeSet[j] = tmp;
	}
	for (unsigned int n = 0; n < icmCodeSet.Size(); n++)
	{
	for (ItO.GoToBegin(); !ItO.IsAtEnd(); ++ItO)
	{
	if (this->m_ICMCodeImage->GetPixel(ItO.GetIndex()) == icmCodeSet[n])
	{
	maxPosteriorSum += this->PerformLocalLabelingUpdate(ItO);
	}
	}
	}
	itkDebugMacro("ICM posterior probability sum: " << maxPosteriorSum);
	}
	}

In relation to the "-e 1" parameter. As I understand the "Atropos --help" text:

ANTs/Examples/Atropos.cxx

Lines 1536 to 1544 in 77b8f2a

	std::string description = std::string("The propagation of each prior label can be controlled ") +
	std::string("by the lambda and boundary probability parameters. The ") +
	std::string("latter parameter is the probability (in the range ") +
	std::string("[0,1]) of the label on the boundary which increases linearly ") +
	std::string("to a maximum value of 1.0 in the interior of the labeled ") +
	std::string("region. The former parameter dictates the exponential ") +
	std::string("decay of probability propagation outside the labeled ") +
	std::string("region from the boundary probability, i.e. ") +
	std::string("boundaryProbability*exp( -lambda * distance ).");

and from the code in:

ANTs/ImageSegmentation/antsAtroposSegmentationImageFilter.h

Lines 611 to 617 in 77b8f2a

	/**
	* Set the label propagation type. Euclidean distance uses the Maurer distance
	* transform to calculate the distance transform image. Otherwise the fast
	* marching filter is used to produce the geodesic distance. The former option
	* is faster but for non-Euclidean shapes (such as the cortex), it might be
	* more accurate to use the latter option. Default = false.
	*/

the best results are found when "-e 1", which is in line what we have observed.

The exact difference in computation between the two can be seen here:

ANTs/ImageSegmentation/antsAtroposSegmentationImageFilter.hxx

Lines 2433 to 2513 in 77b8f2a

	if (this->m_UseEuclideanDistanceForPriorLabels)
	{
	typedef SignedMaurerDistanceMapImageFilter<RealImageType, RealImageType> DistancerType;
	typename DistancerType::Pointer distancer = DistancerType::New();
	distancer->SetInput(thresholder->GetOutput());
	distancer->SetSquaredDistance(false);
	distancer->SetUseImageSpacing(true);
	distancer->SetInsideIsPositive(false);
	distancer->Update();
	distanceImage = distancer->GetOutput();
	}
	else
	{
	typedef BinaryContourImageFilter<RealImageType, RealImageType> ContourFilterType;
	typename ContourFilterType::Pointer contour = ContourFilterType::New();
	contour->SetInput(thresholder->GetOutput());
	contour->FullyConnectedOff();
	contour->SetBackgroundValue(0);
	contour->SetForegroundValue(1);
	contour->Update();
	typedef FastMarchingImageFilter<RealImageType, RealImageType> FastMarchingFilterType;
	typename FastMarchingFilterType::Pointer fastMarching = FastMarchingFilterType::New();
	typedef CastImageFilter<MaskImageType, RealImageType> CasterType;
	typename CasterType::Pointer caster = CasterType::New();
	if (this->GetMaskImage())
	{
	caster->SetInput(const_cast<MaskImageType *>(this->GetMaskImage()));
	caster->Update();
	fastMarching->SetInput(caster->GetOutput());
	}
	else
	{
	fastMarching->SetSpeedConstant(1.0);
	fastMarching->SetOverrideOutputInformation(true);
	fastMarching->SetOutputOrigin(this->GetOutput()->GetOrigin());
	fastMarching->SetOutputSpacing(this->GetOutput()->GetSpacing());
	fastMarching->SetOutputRegion(this->GetOutput()->GetRequestedRegion());
	fastMarching->SetOutputDirection(this->GetOutput()->GetDirection());
	}
	typedef typename FastMarchingFilterType::NodeContainer NodeContainer;
	typedef typename FastMarchingFilterType::NodeType NodeType;
	typename NodeContainer::Pointer trialPoints = NodeContainer::New();
	trialPoints->Initialize();
	unsigned long trialCount = 0;
	ImageRegionIteratorWithIndex<RealImageType> ItC(contour->GetOutput(),
	contour->GetOutput()->GetRequestedRegion());
	for (ItC.GoToBegin(); !ItC.IsAtEnd(); ++ItC)
	{
	if (Math::FloatAlmostEqual(ItC.Get(), contour->GetForegroundValue()))
	{
	NodeType node;
	node.SetValue(0.0);
	node.SetIndex(ItC.GetIndex());
	trialPoints->InsertElement(trialCount++, node);
	}
	}
	fastMarching->SetTrialPoints(trialPoints);
	fastMarching->SetStoppingValue(NumericTraits<RealType>::max());
	// fastMarching->SetTopologyCheck( FastMarchingFilterType::None );
	fastMarching->Update();
	ImageRegionIterator<RealImageType> ItT(thresholder->GetOutput(),
	thresholder->GetOutput()->GetRequestedRegion());
	ImageRegionIterator<RealImageType> ItF(fastMarching->GetOutput(),
	fastMarching->GetOutput()->GetRequestedRegion());
	for (ItT.GoToBegin(), ItF.GoToBegin(); !ItT.IsAtEnd(); ++ItT, ++ItF)
	{
	if (Math::FloatAlmostEqual(ItT.Get(), NumericTraits<float>::OneValue()))
	{
	ItF.Set(-ItF.Get());
	}
	}
	distanceImage = fastMarching->GetOutput();
	}

[gdevenyi]

Thanks @br-cpvc as I suspected reading the code is still the only way to determine some of the effects here. The comment about the distance calculation is particularly helpful. Cheers.

[br-cpvc]

@cookpa I initially thought that the "To do list" that you inserted was for your self or someone on your team, but as no updates have been posted I have started to wonder if it was ment for me? Are you waiting for me to do something or what is the current status?

[cookpa]

It's for anyone who has time to do it, you're welcome to take it on if you'd like.

Since you have exposed some N4 parameters, it might be nice to add the initial mesh resolution / spacing (currently hard coded as N4_BSPLINE_PARAMS). That would make the script more applicable to smaller brains.

BTW, if you're benchmarking performance, there have been recent changes in ITK that might affect multi-threading performance in N4 (#1017).

[br-cpvc]

@cookpa thank you for the comment, and thank you for making me aware of the recent patch (#1017) affect multi-threading performance in N4.

I have looked into what it takes to align the input parameters between the two scripts.

With the current getopts argument parser it seems impossible to add and/or align the input parameters between antsAtroposN4.sh and antsCorticalThickness.sh - this is because of two things:

• getopts only supports single-letter options parsing
• antsCorticalThickness.sh have already used up all single-letter option letters (except -h)

I think the best will be to make it possible to use long options, but as getopts does not support this, another approach to parsing the arguments must be used.

I first looked at the getopt parameter parser, which supports long options. But of what I could find, the getopt package is not recommended, see: http://mywiki.wooledge.org/ComplexOptionParsing the only other option I could find was to implement the parsing directly/manually inside the script, which I have tried based on the example: https://mywiki.wooledge.org/BashFAQ/035#Manual_loop

I have redone my implementation based on the current head of the ANTs main branch, and included a new cli option for setting the N4_BSPLINE_PARAMS as you suggested. The reimplementation is currently in master...br-cpvc:ANTs:upstream_20221110_0945

What do you think of the approach? Is this in the right direction or should it be done another way?

[gdevenyi]

The solution for complex bash parsing is https://argbash.readthedocs.io/en/latest/ and its fantastic.

[cookpa]

Yeah, I'm not sure it's feasible to make every option used within antsCorticalThickness.sh available, without completely rewriting the option parsing. I really just want to verify that changes to antsAtroposN4.sh don't result in changes how that script behaves when called from within antsCorticalThickness.sh

[br-cpvc]

@cookpa Ah, okay, then I completely misunderstood what you were after :)

I have now taken another go at it.

I have pushed another patch to the current PR branch which adds the N4_BSPLINE_PARAMS parameters as you previously suggested. After that, I made a test trying to show that any script that uses antsAtroposN4.sh should expect the same results as before the changes.

The test shows that antsAtroposN4.sh produces the same results (it does not alter the output). By executing it two times, before and after the change to the Atropos parameters and their default values (ATROPOS_SEGMENTATION_ICM=1 and ATROPOS_SEGMENTATION_USE_EUCLIDEAN_DISTANCE=0) which were added to the exe_segmentation statement in antsAtroposN4.sh master...br-cpvc:ANTs:upstream_20220906_1020#diff-b35ac52557562254e75ede77f4c181598f938e248bcfcd8d7c4836c7a7d1048dL502-L503

The test:

• execute antsAtroposN4.sh with -u 0 using git checkout before the changes (e043d8a)
• execute antsAtroposN4.sh with -u 0 using git checkout after the changes (fcb7aa6)
• compare the output using md5sum

Execution of the test

wget https://openneuro.org/crn/datasets/ds004288/snapshots/1.0.1/files/sub-0303:anat:sub-0303_T1w.nii.gz

git checkout e043d8a690be769cff7656f3747bb803cca7d24a
mkdir -p processing_output/baseline
./Scripts/antsAtroposN4.sh -d 3 -u 0 -a sub-0303:anat:sub-0303_T1w.nii.gz -x sub-0303:anat:sub-0303_T1w.nii.gz -o processing_output/baseline/


git checkout fcb7aa60c9dd036f884e3d6f28e2c0f6fef488d0
mkdir -p processing_output/upstream_20220906_1020
./Scripts/antsAtroposN4.sh -d 3 -u 0 -a sub-0303:anat:sub-0303_T1w.nii.gz -x sub-0303:anat:sub-0303_T1w.nii.gz -o processing_output/upstream_20220906_1020/

md5sum processing_output/*/*
4404b487734935d2bb5a455e5f772e7b  processing_output/baseline/Segmentation.nii.gz
b2f3be2cf952cb124edc94096d64991e  processing_output/baseline/Segmentation0N4.nii.gz
dc4780156a987295ee98ba66afb30b59  processing_output/baseline/SegmentationConvergence.txt
67772085703eba7367614918f4548792  processing_output/baseline/SegmentationPosteriors1.nii.gz
1385edf67ea6c11d05d061bfea1c4476  processing_output/baseline/SegmentationPosteriors2.nii.gz
472f5cc135ed93b0fae8a615966086d9  processing_output/baseline/SegmentationPosteriors3.nii.gz

4404b487734935d2bb5a455e5f772e7b  processing_output/upstream_20220906_1020/Segmentation.nii.gz
b2f3be2cf952cb124edc94096d64991e  processing_output/upstream_20220906_1020/Segmentation0N4.nii.gz
dc4780156a987295ee98ba66afb30b59  processing_output/upstream_20220906_1020/SegmentationConvergence.txt
67772085703eba7367614918f4548792  processing_output/upstream_20220906_1020/SegmentationPosteriors1.nii.gz
1385edf67ea6c11d05d061bfea1c4476  processing_output/upstream_20220906_1020/SegmentationPosteriors2.nii.gz
472f5cc135ed93b0fae8a615966086d9  processing_output/upstream_20220906_1020/SegmentationPosteriors3.nii.gz

Both outputs are exacly the same.

I also tried changing the ICM parameter to see that it actually change the results:

./Scripts/antsAtroposN4.sh -d 3 -u 0 -a sub-0303:anat:sub-0303_T1w.nii.gz -x sub-0303:anat:sub-0303_T1w.nii.gz -o processing_output/upstream_20220906_1020_disabled/ -i 0

md5sum processing_output/upstream_20220906_1020_disabled/*
215573439fa491dd3b49bd9617d0494f  processing_output/upstream_20220906_1020_disabled/Segmentation.nii.gz
b2f3be2cf952cb124edc94096d64991e  processing_output/upstream_20220906_1020_disabled/Segmentation0N4.nii.gz
8b8cd3587688ea93b81b30f45151b7b8  processing_output/upstream_20220906_1020_disabled/SegmentationConvergence.txt
d9e21f674b42a3d4a3014da272fbe70d  processing_output/upstream_20220906_1020_disabled/SegmentationPosteriors1.nii.gz
be2190411d00843e6d78c54e7e5f7144  processing_output/upstream_20220906_1020_disabled/SegmentationPosteriors2.nii.gz
503d74d7485987bc634ff394d421e9d5  processing_output/upstream_20220906_1020_disabled/SegmentationPosteriors3.nii.gz

What do you think, is this test deep enough to cover what you were asking?

[cookpa]

Yes that looks good to me, thanks. I did a similar test and got identical results with the new script. So I think this is suitable to merge.