Comparing ITK images efficiently.

Beginner Questions
1

Hello everyone,

I was wondering what would be the fastest way in comparing 2 ITK images irrespective of how large the images are. For instance if I have 2 images both 1000x1000 and I would like to make sure they are exactly the same, would it be naive to do such thing as depicted below where I get the region for both images and manually check each pixel?

// 'InputPixelType' -> is just a template for the primitive type ie: int32_t, double, float et al.
// both readers are just
// typename itk::ImageSeriesReader<itk::Image<InputPixelType, kDimension>>::Pointer
// kDimension = 3 (defined in my code)

typename itk::Image<InputPixelType, kDimension>::Pointer image = reader->GetOutput();
typename itk::Image<InputPixelType, kDimension>::Pointer image2 = reader2->GetOutput();

// then I basically check both their spacings, origins and sizes (code omitted here). 
// and Finally I do something like this

typedef itk::ImageRegionConstIterator<itk::Image<InputPixelType, kDimension>> ConstIter;
							
ConstIter iter1(image, image->GetRequestedRegion());
ConstIter iter2(image2, image2->GetRequestedRegion());
							
for (iter1.GoToBegin(), iter2.GoToBegin(); !iter1.IsAtEnd(); ++iter1, ++iter2)
{
       // IsAlmostEqual -> custom function to compare using a fudge factor
       if(!IsAlmostEqual(iter1.Value(), iter2.Value(), DBL_EPSILON))
	   {
           // boolean success to indicate failure.
		  success = false;
		  break;
	   }
}

If so could someone point me in the right direction as to the best way to approach this issue. I have had a look at ‘InsightSoftwareGuide-Book1-5.0.0’ but couldn’t figure out the best approach. Thanks

2

You should also check their directions. Other than that, you code looks fine. Optionally parallelize your code.

You could also use itk::Testing::ComparisonImageFilter.

3

Hello Dženan,

Thanks for the expeditious response. Is it safe to assume that ‘itk::Testing::ComparisonImageFilter’ is just another approach to my code above and it wouldn’t necessarily be any faster? It seems to me parallelization will produce faster comparisons. Oh and thanks for highlighting ‘checking the directions’ of both images (I’m actually doing so in my code but forgot to put it as part of my comment above) Thanks.

4

ComparisonImageFilter is multi-threaded. You should set tolerance above zero for float images.

5

Great, thanks a bunch. I’ll try your suggestions out and get back to you.

6

Hello Dženan

So I was able to try the above suggested ideas and ‘itk::Testing::ComparisonImageFilter’ basically was slightly slower than my initial code (unless I’m doing something wrong)

// 'InputPixelType' -> is just a template for the primitive type ie: int32_t, double, float et al.
// both readers are just
// typename itk::ImageSeriesReader<itk::Image<InputPixelType, kDimension>>::Pointer
// kDimension = 3 (defined in my code)

using ImageType = itk::Image<InputPixelType, kDimension>;
using DiffType = itk::Testing::ComparisonImageFilter<ImageType, ImageType>;
							
typename DiffType::Pointer difference = DiffType::New();
difference->SetValidInput(image);
difference->SetTestInput(image2);
difference->SetDifferenceThreshold(DBL_EPSILON);
difference->SetToleranceRadius(0.1); // set above 0.1 for float 
difference->UpdateLargestPossibleRegion();
							
bool differenceFailed = false;
int numberOfPixelsTolerance = 5;
							
const double averageIntensityDifference = difference->GetTotalDifference();
							
const unsigned long numberOfPixelsWithDifferences = difference->GetNumberOfPixelsWithDifferences();
							
if (averageIntensityDifference > 0.0)
{
	if (static_cast<int>(numberOfPixelsWithDifferences) > numberOfPixelsTolerance)
	{
		differenceFailed = true;
	}
	else
	{
		differenceFailed = false;
	}
}
else
{
	differenceFailed = false;
}
							
if (differenceFailed)
{
	success = false;
}

whereas when I tried parallelizing it produced slightly faster than my initial code but not that impressive, I’m not sure if there are some settings that needs to be executed to reap the benefits of the MultiThreader. Parallelize code can be seen below.

// 'InputPixelType' -> is just a template for the primitive type ie: int32_t, double, float et al.
// both readers are just
// typename itk::ImageSeriesReader<itk::Image<InputPixelType, kDimension>>::Pointer
// kDimension = 3 (defined in my code)


itk::MultiThreaderBase::Pointer	multiThread = itk::MultiThreaderBase::New();
							
// Typically this is 16 by default.
//itk::ThreadIdType numGlobalThread = multiThread->GetGlobalDefaultNumberOfThreads();
itk::ThreadIdType maxGlobalThread = multiThread->GetGlobalMaximumNumberOfThreads();
							
// trying 'maxGlobalThread'
itk::ThreadIdType globalThreadsToUse = maxGlobalThread ; // numGlobalThread * 4;
globalThreadsToUse = (maxGlobalThread < globalThreadsToUse) ? maxGlobalThread - 2 : globalThreadsToUse;
							
std::cout << "\n\n\t # of Threads: " << globalThreadsToUse << "\n" << endl;
							
multiThread->SetGlobalDefaultNumberOfThreads(globalThreadsToUse);
							
using ImageType = itk::Image<InputPixelType, kDimension>;
							
auto compare = [&success, image, image2](const typename ImageType::RegionType& region) -> void
{
	// Check voxels to see if they match.
	typedef itk::ImageRegionConstIterator<ImageType> ConstIter;
	ConstIter iter1(image, region);
	ConstIter iter2(image2, region);
								
	for (; !iter1.IsAtEnd(); ++iter1, ++iter2)
	{
		if(!IsAlmostEqual(iter1.Value(), iter2.Value(), DBL_EPSILON))
		{
	        std::cout << "\n[Parallelizing]: In comparison (" << iter1.Value() << " is != " << iter2.Value() << ")" << std::endl;
			success = false;
			break;
		}
	}
								
	if (!success){
		return;
	}
};
							
multiThread->ParallelizeImageRegion<kDimension>(image->GetBufferedRegion(), compare, nullptr);
7

I guess the main conclusion here is that this problem is memory speed limited (usually called memory bound), not CPU bound, so it does not benefit a lot from parallelization.

ComparisonImageFilter does more things than your original code.

This is what might be making it slower.

8

Thanks for the feedback, I did try just 0.0 as the tolerance default and it made no difference. I guess I’m not quite understanding what the tolerance does.

9

ToleranceRadius allows to compare pixel’s value to neighboring pixels, within the radius. I guess this needs to be >=1 to have an effect.

ComparisonImageFilter computes mean of the difference, so that might be the reason for running slower.

10

Thanks a lot, I’ll keep that in mind and test with a tolerance larger than 1.0. Hopefully I’ll see better results.