diff --git a/Common/CostFunctions/itkParzenWindowHistogramImageToImageMetric.h b/Common/CostFunctions/itkParzenWindowHistogramImageToImageMetric.h
index fefa871e6..7bad0cc3b 100644
--- a/Common/CostFunctions/itkParzenWindowHistogramImageToImageMetric.h
+++ b/Common/CostFunctions/itkParzenWindowHistogramImageToImageMetric.h
@@ -20,6 +20,7 @@
 
 #include "itkAdvancedImageToImageMetric.h"
 #include "itkKernelFunctionBase2.h"
+#include <mutex>
 #include <vector>
 
 
@@ -504,6 +505,8 @@ class ITK_TEMPLATE_EXPORT ParzenWindowHistogramImageToImageMetric
   bool          m_UseExplicitPDFDerivatives{ true };
   bool          m_UseFiniteDifferenceDerivative{ false };
   double        m_FiniteDifferencePerturbation{ 1.0 };
+
+  std::mutex m_Mutex{};
 };
 
 } // end namespace itk
diff --git a/Common/CostFunctions/itkParzenWindowHistogramImageToImageMetric.hxx b/Common/CostFunctions/itkParzenWindowHistogramImageToImageMetric.hxx
index 572d4b7ea..b8f9d8b93 100644
--- a/Common/CostFunctions/itkParzenWindowHistogramImageToImageMetric.hxx
+++ b/Common/CostFunctions/itkParzenWindowHistogramImageToImageMetric.hxx
@@ -195,7 +195,7 @@ ParzenWindowHistogramImageToImageMetric<TFixedImage, TMovingImage>::InitializeHi
    */
   m_JointPDF = JointPDFType::New();
   m_JointPDF->SetRegions(JointPDFSizeType{ m_NumberOfMovingHistogramBins, m_NumberOfFixedHistogramBins });
-  m_JointPDF->Allocate();
+  m_JointPDF->AllocateInitialized();
 
   if (this->GetUseDerivative())
   {
@@ -1046,6 +1046,21 @@ ParzenWindowHistogramImageToImageMetric<TFixedImage, TMovingImage>::ThreadedComp
   m_ParzenWindowHistogramGetValueAndDerivativePerThreadVariables[threadId].st_NumberOfPixelsCounted =
     numberOfPixelsCounted;
 
+  {
+    const std::lock_guard<std::mutex> lockGuard(m_Mutex);
+
+    /** Accumulate joint histogram. */
+    const ImageBufferRange imageBufferRange(*m_JointPDF);
+
+    auto pdfValueIteratorOfThread = ImageBufferRange(*jointPDF).cbegin();
+
+    for (PDFValueType & pdfValue : imageBufferRange)
+    {
+      pdfValue += *pdfValueIteratorOfThread;
+      ++pdfValueIteratorOfThread;
+    }
+  }
+
 } // end ThreadedComputePDFs()
 
 
@@ -1075,21 +1090,6 @@ ParzenWindowHistogramImageToImageMetric<TFixedImage, TMovingImage>::AfterThreade
   /** Compute alpha. */
   m_Alpha = 1.0 / static_cast<double>(Superclass::m_NumberOfPixelsCounted);
 
-  /** Accumulate joint histogram. */
-  const ImageBufferRange imageBufferRange(*m_JointPDF);
-
-  std::fill(imageBufferRange.begin(), imageBufferRange.end(), PDFValueType{});
-
-  for (const auto & perThreadStruct : m_ParzenWindowHistogramGetValueAndDerivativePerThreadVariables)
-  {
-    auto pdfValueIteratorPerThread = ImageBufferRange(*perThreadStruct.st_JointPDF).cbegin();
-
-    for (PDFValueType & pdfValue : imageBufferRange)
-    {
-      pdfValue += *pdfValueIteratorPerThread;
-      ++pdfValueIteratorPerThread;
-    }
-  }
 } // end AfterThreadedComputePDFs()