Sequential run of Elastix producing different result

[Chayan010]

Hi,
I was trying to run Elastix more than once in sequence with two different fixed and moving images. The result of first run is fine but for second run the output is not same, and even the transformation matrix is different
Even I have tried passing same fixed and moving image to check but second result is different.

Kindly suggest.
Below is the model code I have included in my Code:

#include
#include

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkRoundImageFilter.h"
#include "itkImage.h"
#include "itkCastImageFilter.h"
#include "itkParameterMapInterface.h"
#include "itkParameterFileParser.h"
#include "elastixlib.h"

typedef itk::ParameterFileParser ParameterFileParserType;
std::map < std::string, std::vector< std::string > > ParameterFileParser();
void test1();
void test2();

int main(int argc, char* argv[])
{
test1();
test2();

}

void test1()
{
const unsigned int dimension = 3;
typedef float InputPixelType;
typedef float OutputPixelType;
typedef int IntPixelType;
typedef itk::Image<InputPixelType, dimension> InputImageType;
typedef itk::Image<OutputPixelType, dimension> OutputImageType;
typedef itk::Image<IntPixelType, dimension> IntImageType;
typedef itk::ImageFileReader ReaderType;
typedef itk::ImageFileWriter FloatWriterType;
typedef itk::ImageFileWriter UShortWriterType;
typedef itk::RoundImageFilter<OutputImageType, OutputImageType> RoundFilterType;

std::map < std::string, std::vector< std::string > > parameters = ParameterFileParser();

ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName("D:/tests_001/registration_flair_to_mprage/in/moving.nii.gz");
reader->Update();
InputImageType::Pointer fixed_img = reader->GetOutput();
fixed_img->Print(std::cout);

reader->SetFileName("D:/tests_001/registration_flair_to_mprage/in/fixed.nii.gz");
reader->Update();
InputImageType::Pointer moving_img = reader->GetOutput();
moving_img->Print(std::cout);


elastix::ELASTIX* elastixFilter = new elastix::ELASTIX();
std::string outPutDirPath = "D:/tests_001/registration_flair_to_mprage/WMLDOut";
int error = 0;


try
{
	error = elastixFilter->RegisterImages(static_cast<typename itk::DataObject::Pointer>(fixed_img), static_cast<typename itk::DataObject::Pointer>(moving_img), parameters, outPutDirPath, true, true);
}
catch (itk::ExceptionObject& e)
{
	std::cerr << "Parameters/Properties file not suppported" << e << std::endl;

}


OutputImageType* output_image;

if (error == 0)
{
	if (elastixFilter->GetResultImage().IsNotNull())
	{
		output_image = static_cast<OutputImageType*>(elastixFilter->GetResultImage().GetPointer());

		std::string outFile = outPutDirPath + "\\result.nii.gz";

	}
}

delete elastixFilter;

}

void test2()
{
const unsigned int dimension = 3;
typedef float InputPixelType;
typedef float OutputPixelType;
typedef int IntPixelType;
typedef itk::Image<InputPixelType, dimension> InputImageType;
typedef itk::Image<OutputPixelType, dimension> OutputImageType;
typedef itk::Image<IntPixelType, dimension> IntImageType;
typedef itk::ImageFileReader ReaderType;
typedef itk::ImageFileWriter FloatWriterType;
typedef itk::ImageFileWriter UShortWriterType;
typedef itk::RoundImageFilter<OutputImageType, OutputImageType> RoundFilterType;

std::map < std::string, std::vector< std::string > > parameters = ParameterFileParser();

ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName("D:/tests_001/registration_flair_to_mprage/in/moving.nii.gz");
reader->Update();
InputImageType::Pointer fixed_img = reader->GetOutput();
fixed_img->Print(std::cout);

reader->SetFileName("D:/tests_001/registration_flair_to_mprage/in/fixed.nii.gz");
reader->Update();
InputImageType::Pointer moving_img = reader->GetOutput();
moving_img->Print(std::cout);

elastix::ELASTIX* elastixFilter = new elastix::ELASTIX();
std::string outPutDirPath = "D:/tests_001/registration_flair_to_mprage/WMLDOut2";
int error = 0;

try
{
	error = elastixFilter->RegisterImages(static_cast<typename itk::DataObject::Pointer>(fixed_img), static_cast<typename itk::DataObject::Pointer>(moving_img), parameters, outPutDirPath, true, true);
}
catch (itk::ExceptionObject& e)
{
	std::cerr << "Parameters/Properties file not suppported" << e << std::endl;

}

OutputImageType* output_image;

if (error == 0)
{
	if (elastixFilter->GetResultImage().IsNotNull())
	{
		output_image = static_cast<OutputImageType*>(elastixFilter->GetResultImage().GetPointer());

		std::string outFile = outPutDirPath + "\\result.nii.gz";

		
	}
	
}

delete elastixFilter;

}

std::map < std::string, std::vector< std::string > > ParameterFileParser()
{
std::map < std::string, std::vector< std::string > > parameters;
constexpr auto ImageDimension = 3;

const std::pair<std::string, std::string> parameterArray[] =
{
	{"FixedInternalImagePixelType", "float"},
	{"MovingInternalImagePixelType", "float"},
	{"FixedImageDimension", std::to_string(ImageDimension)},
	{"MovingImageDimension", std::to_string(ImageDimension)},
	{"UseDirectionCosines", "true"},
	{"Registration", "MultiResolutionRegistration"},
	{"Interpolator", "BSplineInterpolator"},
	{"ResampleInterpolator", "FinalBSplineInterpolator"},
	{"Resampler", "DefaultResampler"},
	{"FixedImagePyramid", "FixedSmoothingImagePyramid"},
	{"MovingImagePyramid", "MovingSmoothingImagePyramid"},
	{"Optimizer", "AdaptiveStochasticGradientDescent"},
	{"Transform", "AffineTransform"},
	{"Metric", "AdvancedMattesMutualInformation"},
	{"AutomaticScalesEstimation", "true"},
	{"AutomaticTransformInitialization", "true"},
	{"HowToCombineTransforms", "Compose"},
	{"NumberOfHistogramBins", "32"},
	{"ErodeMask", "false"},
	{"NumberOfResolutions", "4"},
	{"ImagePyramidSchedule", "8 8 8  4 4 4  2 2 2  1 1 1"},
	{"MaximumNumberOfIterations", "1200"},
	{"NumberOfSpatialSamples", "2048"},
	{"NewSamplesEveryIteration", "true"},
	{"ImageSampler", "Random"},
	{"BSplineInterpolationOrder", "1"},
	{"FinalBSplineInterpolationOrder", "3"},
	{"DefaultPixelValue", "0"},
	{"WriteResultImage", "true"},
	{"ResultImagePixelType", "int"},
	{"ResultImageFormat", "nii"},
};

for (const auto& pair : parameterArray)
{
	parameters.insert({ pair.first, {pair.second} });
}

return parameters;

}

regards,
Chayan

Below is the link to the issue:
https://groups.google.com/g/elastix-imageregistration/c/0_36mGQrPv8

[rcorredorj]

Hi again. @N-Dekker, @mstaring maybe you guys have ideas regarding this?

I also tried what @Chayan010 indicated in the previous post (also here more details: https://groups.google.com/g/elastix-imageregistration/c/0_36mGQrPv8)

I tested in my system and the issue is very weird. The code that I uploaded in the repository https://github.com/rcorredorj/elastix-cpp-test shows exactly an example that reproduces this issue:

• In the code I perform two sequential calls to the RegisterImages(...) function after doing some simple ITK based operations.
• Despite we use exactly the same images, even by deleting the first elastix instance and creating a new instance for the second call, we get different results! . It is a difference that we cannot explain because** it does not happen when we run twice the official executable on the same images**.

Here the differences in the transformation parameters:

We'll appreciate your thoughts on this issue and thank you again for your support!

Cheers,

RaC

[mstaring]

I suspect that the seed we set in the executable is not set in the library. We will put it on our list of todo's ...

[mstaring]

note however that this is not very problematic behaviour, it is simply related to the randomness in the sampler used by the stochastic optimizer. Stochastic optimizers will walk to the same minimum for different seeds, up to some epsilon. You observe the epsilon here.

[rcorredorj]

@mstaring Thanks again for your answer.
Could you provide us some indications of that seed set in the executable and not in the library?
I can understand the effects of the randomness in the computation, but it is a bit unclear the fact of getting same results while running via the executable vs calling the library despite the fact that both run the same stochastic processes.

[hankst69]

@mstaring
I currently work together with @rcorredorj and @Chayan010 on a project using Elasitix registration as a library.
I also investigated the problem that I would call "reproducability issue" and I got more or less insight in what is going on.
The following relates to Elastix4.90 and ITK4.11 and might help to fix it future Elastix version:

Your suspect is right, that the problem is caused by the use of the RandomSampler in combination with the itk random number generator (MersenneTwisterRandomVariateGenerator).

In order to get reproducable registration results in subsequent elastixLib::RegisterImages() usages we would need to:

• either have the same pseudo random number sequence with every RegisterImages() call
• or use a break condition of the registration iterations that is based on a given quality (currently we always stop on MaximumNumberOfIterations reached but I realized there is also something prepared to stop iterations at MinimumStepSize reached)

Maybe you can give advice how to use the MinimumStepSize stop condition?

Regarding the reset of the random number generator before every use of RegisterImages():
Due implementation details of MersenneTwisterRandomVariateGenerator and its usage from elastix side this is currently not possible to achieve (without hacks)

Your guess that just the initialization of the MersenneTwisterRandomVariateGenerator with specific seed is missing in elastixLib scenario is partly true. In fact that seed initialization is already executed in ealstixLib scenario (within ElastixBase::BeforeAllBase) but it is executed a bit to late because the MersenneTwisterRandomVariateGenerator singleton is already accessed during component instantiation.

But that is not the only problem. The real problem is, that elastix uses the MersenneTwisterRandomVariateGenerator singleton in 2 different ways ~ as random number generator and as seed generator ~.
Only the seed can be reseted but not the increment counter used in MersenneTwisterRandomVariateGenerator::GetNextSeed().

[mstaring]

Great analysis, thank you!

Regarding the minimum stepsize, that parameter is most relevant for the non-stochastic optimizers. For stochastic ones there will be noise on the step size and it is perhaps better to look at a moving average. But that is not currently implemented in elastix, so we always use the max number of iterations.

Regarding the seed, good to hear that it is actually called in the lib as well. Do I understand correctly that in the lib, because the object itself stays alive unlike the command line version, the seed is the same between two consecutive runs, but this increment you mention is not?

[hankst69]

Thanks for the accolade.

Sad to hear that a quality based stop condition is no option.

Yes you understood correctly.
There is that increment counter within the itk MersenneTwisterRandomVariateGenerator that is a private static. So it gets set to zero when the Elastix/Itk module gets loaded into memory and then only ever gets up with every call to GetNextSeed().

So in lib scenario with every new run of Elastix registration the Seed value within the MersenneTwisterRandomVariateGenerator Singleton is the same (due the fix in ElastixBase::BeforeAllBase) and hence that Singleton produces same random number sequence as in first run. But all the seed values fetched from that Singleton will be much higher than with the first run.

You might not be aware to use that GetNextSeed() function because you just inheritate that usage via the itk::ImageRandomConstIteratorWithIndex. Every new instance of this itk class creates a new instance of the MersenneTwisterRandomVariateGenerator and initializes that new instance with the next seed delivered from the global Singleton instance of that same random generator class. So these random generators (which are used to generate the random sample path through the images) will produce different random numbers in every new run beause the seed value will increase with eyery run.

This is how GetNextSeed() is implemented:

MersenneTwisterRandomVariateGenerator::IntegerType
MersenneTwisterRandomVariateGenerator
::GetNextSeed()
{
  IntegerType newSeed = GetInstance()->GetSeed();
  {
    MutexLockHolder<SimpleFastMutexLock> mutexHolder(m_StaticInstanceLock);
    newSeed += m_StaticDiffer++;
  }
  return newSeed;
}

As already mentioned m_StaticDiffer is defined as private static.

So there is no way to access it and reset it from outside.
Also derivation from MersenneTwisterRandomVariateGenerator and injecting this derivate via ItkObjectFactory mechanism does not help any further because the only real interface (virtual) member of that random generator class is just the "virtual double GetVariate()" and this specific function is never used. So all the access to the random generator are going directly to that specific MersenneTwisterRandomVariateGenerator implementation.

The best and cleanest solution would be to change the itk code, means to extend the MersenneTwisterRandomVariateGenerator with a method that allows to reset that increment counter.

Another but not so clean fix could be to change the Elastix implementation in AdaptiveStochasticGradientDescent and to modify/extend the random seed fix in ElastixBase::BeforeAllBase. The idea behind that possible fix would take longer to explain. You can call me back if interested in details.

[mstaring]

Then indeed I would suggest as a first step to talk to the itk people if the increment counter could be exposed by adding set and get functions to this class, or simply a ResetStaticDiffer() function. Would be an easy addition, if they agree this would be useful and does not break something. In my experience they are very welcoming to new ideas and code enhancements.

[rcorredorj]

Hi @N-Dekker,
I found that you had a long thread with the ITK people a bit related to this topic (InsightSoftwareConsortium/ITK#970, not exactly but touching part of this code).
Would you have any additional recommendation for us?
Any way to reset that m_StaticDiffer to make the consecutive calls reproducible?
Any way to trigger the reinitialization of this nested class that keeps this value?
Everywhere is mentioned that it is important to set the Seed to keep some reproducibility (also noted in ITK mailing list and ANTs docs: https://itk.org/pipermail/community/2017-March/012830.html, https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues), but no one mentions this counter that is incremented and seems to keep incrementing in consecutive calls.

Any additional suggestion is very much appreciated.
Thanks Elastix-ers !
RaC

[rcorredorj]

Dear @mstaring and Elastix-ers,

We have a sort of related question to this thread.

When we run elastix from C++, we get the transformation parameters with a high precision number, a float number as expected.
However, when we run the executable, the values of the TransformParameters entry in the TransformParameters.0.txt file seem to be bounded to have only 6 decimals. From some tests we did, we had the impression that the values are not rounded or cropped to the 6th decimal ... we were checking the elastix code, but we couldn't conclude exactly how those values get actually printed. Seems like the Spacing, Origin, Direction, CenterOfRotationPoint and others, were actually treated properly by defining a specific precision of 10 decimals, but not for TransformParameters

Do you have a better understanding of how this is managed for the TransformParameters to get from C++ the same values we get in the txt file ?

Thanks in advance!
RaC

[N-Dekker]

When we run elastix from C++, we get the transformation parameters with a high precision number, a float number as expected.
However, when we run the executable, the values of the TransformParameters entry in the TransformParameters.0.txt file seem to be bounded to have only 6 decimals.

@rcorredorj Hi Ricardo, thanks for elaborating on this issue. Did you already try the latest version from the "develop" branch. Because it should ensure that floating point transform parameters now always have the full precision of a 64-bit double. Both in memory and in the TransformParameters.0.txt file! Specifically by these commits:

1. e1f5dae ENH: Implemented ToString(double) by itk::NumberToString (issue Reconsider string representation of floating points for parameter files and parameter maps #383)
2. af382ee ENH: Sync Transform WriteToFile with CreateTransformParametersMap

HTH, Niels

[mstaring]

Then indeed I would suggest as a first step to talk to the itk people if the increment counter could be exposed by adding set and get functions to this class, or simply a ResetStaticDiffer() function. Would be an easy addition, if they agree this would be useful and does not break something. In my experience they are very welcoming to new ideas and code enhancements.

@N-Dekker just submitted a PR related to this:
InsightSoftwareConsortium/ITK#4247

Do you think this may be useful for this issue as well?

[N-Dekker]

@hankst69 @rcorredorj @Chayan010 @mstaring This week I finally proposed a function for ITK's MersenneTwisterRandomVariateGenerator to reset its m_StaticDiffer. I named it ResetNextSeed():

• ENH: Add MersenneTwisterRandomVariateGenerator::ResetNextSeed() InsightSoftwareConsortium/ITK#4247

(Please let me know if you would prefer to name it "ResetStaticDiffer()".) Sorry for being so late! The function is now merged into the master branch of ITK (InsightSoftwareConsortium/ITK@ecd7bda), and is intended to be included with the next major ITK release (v5.4). Thanks to you all for your feedback!