PERF: Optimize compare image

Summary

Speed-up the performances on compare-image.

On a test case, a rough x3 speed-up has been observed.

Rationale

Comparing the images produced against their baseline is not as efficient as it could be at the end of automated tests.

Implementation Details

The improvements have been done on two axes

Hidden costs of vector images

The DifferenceImageFilter is often used on VectorImages, in particular when doing compare-image.

Creating new VariableLengthVector pixels at each iteration, or casting them will automatically end up in construction and destruction of VLV pixels, and thus allocation and liberation of heap memory. This is not efficient.

The first refactoring aims at minimizing the number of allocation by factorizing out everything that can be (pixel size, typical max value, typical zero value). Resetting value is done by assignment, which has the good property of not inducing any allocation on VLV variables.

Non-localized reductions

The reduce computation done in DifferenceImageFilter was done on some m_ThreadAccumulatedCounter[threadId] member. This has two drawbacks:

• first the compiler cannot know whether we are the only one function to work on that value, it cannot isolate it in a register nor anything similar
• moreover false sharing happens and prevent further optimizations.

As a consequence the counters are now local variables in the ThreadedGenerateData() function, and they are commited into the shared member variables at the end of the execution.

Note: this anti-pattern appears in other places like for instance StreamingCompareImageFilter. This performance related refactoring should be extended to these other cases.

Copyright

The copyright owner is CNES and has signed the ORFEO ToolBox Contributor License Agreement.

---

Check before merging:

• All discussions are resolved
• At least 2 votes from core developers, no vote.
• The feature branch is (reasonably) up-to-date with the base branch
• Dashboard is green
• Copyright owner has signed the ORFEO ToolBox Contributor License Agreement

Modules/IO/TestKernel/include/otbDifferenceImageFilter.hxx

@@ -99,11 +99,11 @@ DifferenceImageFilter<TInputImage, TOutputImage>
 {
   Superclass::GenerateOutputInformation();
   if (this->GetInput(0)->GetNumberOfComponentsPerPixel() != this->GetInput(1)->GetNumberOfComponentsPerPixel())
-    {
+  {
     itkExceptionMacro(<< "Image 1 has " << this->GetInput(
                         0)->GetNumberOfComponentsPerPixel() << " bands, whereas image 2 has " << this->GetInput(
                         1)->GetNumberOfComponentsPerPixel());
-    }
+  }
   this->GetOutput()->SetNumberOfComponentsPerPixel(this->GetInput(0)->GetNumberOfComponentsPerPixel());
 }
 //----------------------------------------------------------------------------
@@ -132,9 +132,9 @@ DifferenceImageFilter<TInputImage, TOutputImage>
 
   // Initialize the temporaries
   for (int i = 0; i < numberOfThreads; ++i)
-    {
+  {
     m_ThreadDifferenceSum.push_back(m_TotalDifference);
-    }
+  }
   m_ThreadNumberOfPixels.Fill(0);
 }
 
@@ -163,14 +163,7 @@ DifferenceImageFilter<TInputImage, TOutputImage>
 
   // Create a radius of pixels.
   RadiusType radius;
-  if (m_ToleranceRadius > 0)
-    {
-    radius.Fill(m_ToleranceRadius);
-    }
-  else
-    {
-    radius.Fill(0);
-    }
+  radius.Fill(std::max(0, m_ToleranceRadius));
 
   // Find the data-set boundary faces.
   FacesCalculator boundaryCalculator;
@@ -179,43 +172,60 @@ DifferenceImageFilter<TInputImage, TOutputImage>
   // Support progress methods/callbacks.
   itk::ProgressReporter progress(this, threadId, threadRegion.GetNumberOfPixels());
 
+  AccumulateType threadDifferenceSum;
+  itk::NumericTraits<AccumulateType>::SetLength(
+    threadDifferenceSum,
+    this->GetInput(0)->GetNumberOfComponentsPerPixel()); // @post: threadDifferenceSum=={0...}
+  unsigned long threadNumberOfPixels = 0;
+
   // Process the internal face and each of the boundary faces.
   for (FaceListIterator face = faceList.begin(); face != faceList.end(); ++face)
-    {
+  {
     SmartIterator test(radius, testImage, *face); // Iterate over test image.
     InputIterator valid(validImage, *face);       // Iterate over valid image.
     OutputIterator out(outputPtr, *face);         // Iterate over output image.
     test.OverrideBoundaryCondition(&nbc);
 
+    // Extract a typical pixel in order to know the exact size of the
+    // pixel, the typical 0, and the typical max
+    valid.GoToBegin();
+    InputPixelType t = valid.Get();
+
+    const auto pixel_size = itk::NumericTraits<InputPixelType>::GetLength(t);
+    const auto out_max    = itk::NumericTraits<OutputPixelType>::max(t);
+    const auto zero       = itk::NumericTraits<OutputPixelType>::ZeroValue(t);
+
+    OutputPixelType minimumDifference = out_max;
+
     for (valid.GoToBegin(), test.GoToBegin(), out.GoToBegin();
          !valid.IsAtEnd();
          ++valid, ++test, ++out)
-      {
+    {
       // Get the current valid pixel.
-      InputPixelType t = valid.Get();
+      t = valid.Get();
 
       // Find the closest-valued pixel in the neighborhood of the test
       // image.
-      OutputPixelType minimumDifference = itk::NumericTraits<OutputPixelType>::max(t);
+      minimumDifference = out_max; // reset by assignment, avoid allocation in VLV case
       unsigned int    neighborhoodSize = test.Size();
       for (unsigned int i = 0; i < neighborhoodSize; ++i)
-        {
-        // Use the RealType for the difference to make sure we get the
-        // sign.
-        RealType difference = static_cast<RealType>(t) - static_cast<RealType>(test.GetPixel(i));
+      {
 
-        for (unsigned int j = 0; j < itk::NumericTraits<RealType>::GetLength(difference) ; ++j)
-          {
+        for (unsigned int j = 0; j < pixel_size ; ++j)
+        {
+          // Use the RealType for the difference to make sure we get
+          // the sign.
+          // Work on component independently in order to avoid
+          // allocation of VLV pixels
           ScalarRealType d = static_cast<ScalarRealType>(
-            itk::DefaultConvertPixelTraits<RealType>::GetNthComponent(j,difference));
-          if (d < 0)
-            {
-            d *= -1;
-            }
+            itk::DefaultConvertPixelTraits<InputPixelType>::GetNthComponent(j,t)
+            - itk::DefaultConvertPixelTraits<InputPixelType>::GetNthComponent(j,test.GetPixel(i))
+            );
+          d = std::abs(d);
           ScalarRealType m = static_cast<ScalarRealType>(
             itk::DefaultConvertPixelTraits<OutputPixelType>::GetNthComponent(j,minimumDifference));
           if (d < m)
-            {
+          {
             itk::DefaultConvertPixelTraits<OutputPixelType>::SetNthComponent(
               j,
               minimumDifference,
@@ -224,55 +234,61 @@ DifferenceImageFilter<TInputImage, TOutputImage>
 //             std::cout << std::setprecision(16) << t << std::endl;
 //             std::cout << std::setprecision(16) << test.GetPixel(i) << std::endl;
 //             std::cout << "----------------------" << std::endl;
-            }
           }
         }
+      }
 
       //for complex and vector type. FIXME: module might be better
 //        ScalarRealType tMax=std::abs(t[0]);
       ScalarRealType tMax = 0.01; //Avoiding the 0 case for neighborhood computing
       // NB: still more restrictive than before for small values.
-      for (unsigned int j = 0; j < itk::NumericTraits<InputPixelType>::GetLength(t); ++j)
-        {
+      for (unsigned int j = 0; j < pixel_size ; ++j)
+      {
         ScalarRealType tc = static_cast<ScalarRealType>(
           itk::DefaultConvertPixelTraits<InputPixelType>::GetNthComponent(j,t));
-        if (std::abs(tc) > tMax) tMax = std::abs(tc);
-        }
+        tMax = std::max(std::abs(tc), tMax);
+      }
 
       // Check if difference is above threshold
       // the threshold is interpreted as relative to the value
       bool isDifferent = false;
 
-      for (unsigned int j = 0; j < itk::NumericTraits<OutputPixelType>::GetLength(minimumDifference); ++j)
-        {
+      for (unsigned int j = 0; j < pixel_size ; ++j)
+      {
         ScalarRealType m = static_cast<ScalarRealType>(
           itk::DefaultConvertPixelTraits<OutputPixelType>::GetNthComponent(j,minimumDifference));
         if (m > m_DifferenceThreshold * tMax)
-          {
+        {
 //           std::cout << std::setprecision(16) << minimumDifference[j] << std::endl;
           isDifferent = true;
-          }
+          break;
         }
+      }
 
       if (isDifferent)
-        {
+      {
         // Store the minimum difference value in the output image.
         out.Set(minimumDifference);
 
-        // Update difference image statistics.
-        m_ThreadDifferenceSum[threadId] += minimumDifference;
-        m_ThreadNumberOfPixels[threadId]++;
-        }
+        // Update local difference image statistics.
+        threadDifferenceSum += minimumDifference;
+        threadNumberOfPixels++;
+      }
       else
-        {
+      {
         // Difference is below threshold.
-        out.Set(itk::NumericTraits<OutputPixelType>::ZeroValue(minimumDifference));
-        }
+        out.Set(zero);
+      }
 
       // Update progress.
       progress.CompletedPixel();
-      }
     }
+  }
+
+  // Update global difference image statistics.
+  m_ThreadDifferenceSum[threadId]  = threadDifferenceSum;
+  m_ThreadNumberOfPixels[threadId] = threadNumberOfPixels;
+
 }
 //----------------------------------------------------------------------------
 template <class TInputImage, class TOutputImage>
@@ -283,10 +299,10 @@ DifferenceImageFilter<TInputImage, TOutputImage>
   // Set statistics about difference image.
   int numberOfThreads = this->GetNumberOfThreads();
   for (int i = 0; i < numberOfThreads; ++i)
-    {
+  {
     m_TotalDifference += m_ThreadDifferenceSum[i];
     m_NumberOfPixelsWithDifferences += m_ThreadNumberOfPixels[i];
-    }
+  }
 
   // Get the total number of pixels processed in the region.
   // This is different from the m_TotalNumberOfPixels which
@@ -296,7 +312,6 @@ DifferenceImageFilter<TInputImage, TOutputImage>
   unsigned int          numberOfPixels = region.GetNumberOfPixels();
 
   // Calculate the mean difference.
-
   m_MeanDifference = m_TotalDifference / numberOfPixels;
 }