Exposes vtkHighestDensityRegionStatistics::ComputeHDR() to the public API

Filters/Statistics/vtkHighestDensityRegionsStatistics.cxx

@@ -196,8 +196,17 @@ void vtkHighestDensityRegionsStatistics::Derive(vtkMultiBlockDataSet*)
 // ----------------------------------------------------------------------
 double vtkHighestDensityRegionsStatistics::ComputeHDR(vtkDataArray *inObs,
                                                       vtkDataArray *outDensity)
+{
+  return ComputeHDR(inObs, inObs, outDensity);
+}
+
+// ----------------------------------------------------------------------
+double vtkHighestDensityRegionsStatistics
+::ComputeHDR(vtkDataArray *inObs, vtkDataArray *inPointsOfInterest,
+             vtkDataArray *outDensity)
 {
   vtkIdType nbObservations = inObs->GetNumberOfTuples();
+  vtkIdType nbPoints = inPointsOfInterest->GetNumberOfTuples();
 
   if (nbObservations == 0)
     {
@@ -208,29 +217,31 @@ double vtkHighestDensityRegionsStatistics::ComputeHDR(vtkDataArray *inObs,
 
   double denom = 1.0 / static_cast<double>(nbObservations);
 
-  // Let's compute the HDR for each points of the observations
-  for (vtkIdType i = 0; i < nbObservations; i++)
+  // Let's compute the HDR for each points of interest
+  for (vtkIdType i = 0; i < nbPoints; i++)
     {
     double currentXi[2];
     double currentXj[2];
     double hdr = 0.0;
 
     // We are working in a bivariate model.
-    inObs->GetTuple(i, currentXi);
+    inPointsOfInterest->GetTuple(i, currentXi);
     // Sum all gaussian kernel
     for (vtkIdType j = 0; j < nbObservations; j++)
       {
+
+      inObs->GetTuple(j, currentXj);
+
+      const double deltaX = currentXi[0] - currentXj[0];
+      const double deltaY = currentXi[1] - currentXj[1];
       // Avoid case where point is compared to itself
-      if (i == j)
+      if (deltaX == 0. && deltaY == 0.)
         {
         continue;
         }
-      inObs->GetTuple(j, currentXj);
-
       hdr += this->ComputeSmoothGaussianKernel(
         inObs->GetNumberOfComponents(),
-        currentXi[0] - currentXj[0],
-        currentXi[1] - currentXj[1]);
+        deltaX, deltaY);
       }
     double d = denom * hdr;
     outDensity->SetTuple1(i, d);

Filters/Statistics/vtkHighestDensityRegionsStatistics.h

@@ -69,6 +69,22 @@ public:
   vtkGetVectorMacro(SmoothHC2, double, 2);
   vtkSetVectorMacro(SmoothHC2, double, 2);
 
+  // Description:
+  // Fill outDensity with density vector that is computed from
+  // inObservations values. This method uses a Gaussian kernel.
+  // For n observations and with X an observation point:
+  // f(X) = (1 / n) * Sum(KH(X -Xi)) for (i = 1 to n).
+  // Look ComputeSmoothGaussianKernel for KH kernel definition.
+  double ComputeHDR(vtkDataArray *inObservations, vtkDataArray *outDensity);
+
+  // Description:
+  // Fill outDensity with density vector defined by inPOI and computed from
+  // the inObs values. This method uses a Gaussian kernel.
+  // For n observations and with X an observation point:
+  // f(X) = (1 / n) * Sum(KH(X -Xi)) for (i = 1 to n).
+  // Look ComputeSmoothGaussianKernel for KH kernel definition.
+  double ComputeHDR(vtkDataArray *inObs, vtkDataArray* inPOI,
+                    vtkDataArray *outDensity);
 protected:
   vtkHighestDensityRegionsStatistics();
   ~vtkHighestDensityRegionsStatistics();
@@ -104,14 +120,6 @@ protected:
                                    AssessFunctor*&) { return; }
 //ETX
 
-  // Description:
-  // Fill outDensity with density vector that is computed from
-  // inObservations values. This method uses a Gaussian kernel.
-  // For n observations and with X an observation point:
-  // f(X) = (1 / n) * Sum(KH(X -Xi)) for (i = 1 to n).
-  // Look ComputeSmoothGaussianKernel for KH kernel definition.
-  double ComputeHDR(vtkDataArray *inObservations, vtkDataArray *outDensity);
-
   // Description:
   // Store the smooth matrix parameter H. Specify a smooth direction
   // for the Gaussian kernel.