1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkGaussianKernel.cxx
5
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 All rights reserved.
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
13
14 =========================================================================*/
15 #include "vtkGaussianKernel.h"
16 #include "vtkAbstractPointLocator.h"
17 #include "vtkObjectFactory.h"
18 #include "vtkIdList.h"
19 #include "vtkDoubleArray.h"
20 #include "vtkDataSet.h"
21 #include "vtkPointData.h"
22 #include "vtkMath.h"
23 #include "vtkMathUtilities.h"
24
25 vtkStandardNewMacro(vtkGaussianKernel);
26
27 //----------------------------------------------------------------------------
vtkGaussianKernel()28 vtkGaussianKernel::vtkGaussianKernel()
29 {
30 this->Sharpness = 2.0;
31 this->F2 = this->Sharpness / this->Radius;
32 }
33
34
35 //----------------------------------------------------------------------------
36 vtkGaussianKernel::~vtkGaussianKernel() = default;
37
38 //----------------------------------------------------------------------------
39 void vtkGaussianKernel::
Initialize(vtkAbstractPointLocator * loc,vtkDataSet * ds,vtkPointData * pd)40 Initialize(vtkAbstractPointLocator *loc, vtkDataSet *ds, vtkPointData *pd)
41 {
42 this->Superclass::Initialize(loc, ds, pd);
43
44 this->F2 = this->Sharpness / this->Radius;
45 this->F2 = this->F2 * this->F2;
46 }
47
48 //----------------------------------------------------------------------------
49 vtkIdType vtkGaussianKernel::
ComputeWeights(double x[3],vtkIdList * pIds,vtkDoubleArray * prob,vtkDoubleArray * weights)50 ComputeWeights(double x[3], vtkIdList *pIds, vtkDoubleArray *prob,
51 vtkDoubleArray *weights)
52 {
53 vtkIdType numPts = pIds->GetNumberOfIds();
54 double d2, y[3], sum = 0.0;
55 weights->SetNumberOfTuples(numPts);
56 double *p = (prob ? prob->GetPointer(0) : nullptr);
57 double *w = weights->GetPointer(0);
58 double f2=this->F2;
59
60 for (vtkIdType i=0; i<numPts; ++i)
61 {
62 vtkIdType id = pIds->GetId(i);
63 this->DataSet->GetPoint(id,y);
64 d2 = vtkMath::Distance2BetweenPoints(x,y);
65
66 if ( vtkMathUtilities::FuzzyCompare(d2, 0.0, std::numeric_limits<double>::epsilon()*256.0 )) //precise hit on existing point
67 {
68 pIds->SetNumberOfIds(1);
69 pIds->SetId(0,id);
70 weights->SetNumberOfTuples(1);
71 weights->SetValue(0,1.0);
72 return 1;
73 }
74 else
75 {
76 w[i] = (p ? p[i]*exp(-f2 * d2) : exp(-f2 * d2));
77 sum += w[i];
78 }
79 }//over all points
80
81 // Normalize
82 if ( this->NormalizeWeights && sum != 0.0 )
83 {
84 for (vtkIdType i=0; i<numPts; ++i)
85 {
86 w[i] /= sum;
87 }
88 }
89
90 return numPts;
91 }
92
93 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)94 void vtkGaussianKernel::PrintSelf(ostream& os, vtkIndent indent)
95 {
96 this->Superclass::PrintSelf(os,indent);
97
98 os << indent << "Sharpness: " << this->GetSharpness() << endl;
99 }
100