1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkStatisticalOutlierRemoval.cxx
5
6 Copyright (c) Kitware, Inc.
7 All rights reserved.
8 See LICENSE file 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 "vtkStatisticalOutlierRemoval.h"
16
17 #include "vtkObjectFactory.h"
18 #include "vtkAbstractPointLocator.h"
19 #include "vtkStaticPointLocator.h"
20 #include "vtkPointSet.h"
21 #include "vtkPoints.h"
22 #include "vtkIdList.h"
23 #include "vtkSMPTools.h"
24 #include "vtkSMPThreadLocalObject.h"
25 #include "vtkMath.h"
26
27 vtkStandardNewMacro(vtkStatisticalOutlierRemoval);
28 vtkCxxSetObjectMacro(vtkStatisticalOutlierRemoval,Locator,vtkAbstractPointLocator);
29
30 //----------------------------------------------------------------------------
31 // Helper classes to support efficient computing, and threaded execution.
32 namespace {
33
34 //----------------------------------------------------------------------------
35 // The threaded core of the algorithm (first pass)
36 template <typename T>
37 struct ComputeMeanDistance
38 {
39 const T *Points;
40 vtkAbstractPointLocator *Locator;
41 int SampleSize;
42 float *Distance;
43 double Mean;
44
45 // Don't want to allocate working arrays on every thread invocation. Thread local
46 // storage lots of new/delete.
47 vtkSMPThreadLocalObject<vtkIdList> PIds;
48 vtkSMPThreadLocal<double> ThreadMean;
49 vtkSMPThreadLocal<vtkIdType> ThreadCount;
50
ComputeMeanDistance__anon3454176f0111::ComputeMeanDistance51 ComputeMeanDistance(T *points, vtkAbstractPointLocator *loc, int size, float *d) :
52 Points(points), Locator(loc), SampleSize(size), Distance(d), Mean(0.0)
53 {
54 }
55
56 // Just allocate a little bit of memory to get started.
Initialize__anon3454176f0111::ComputeMeanDistance57 void Initialize()
58 {
59 vtkIdList*& pIds = this->PIds.Local();
60 pIds->Allocate(128); //allocate some memory
61
62 double &threadMean = this->ThreadMean.Local();
63 threadMean = 0.0;
64
65 vtkIdType &threadCount = this->ThreadCount.Local();
66 threadCount = 0;
67 }
68
69 // Compute average distance for each point, plus accumulate summation of
70 // mean distances and count (for averaging in the Reduce() method).
operator ()__anon3454176f0111::ComputeMeanDistance71 void operator() (vtkIdType ptId, vtkIdType endPtId)
72 {
73 const T *px = this->Points + 3*ptId;
74 const T *py;
75 double x[3], y[3];
76 vtkIdList*& pIds = this->PIds.Local();
77 double &threadMean = this->ThreadMean.Local();
78 vtkIdType &threadCount = this->ThreadCount.Local();
79
80 for ( ; ptId < endPtId; ++ptId)
81 {
82 x[0] = static_cast<double>(*px++);
83 x[1] = static_cast<double>(*px++);
84 x[2] = static_cast<double>(*px++);
85
86 // The method FindClosestNPoints will include the current point, so
87 // we increase the sample size by one.
88 this->Locator->FindClosestNPoints(this->SampleSize+1, x, pIds);
89 vtkIdType numPts = pIds->GetNumberOfIds();
90
91 double sum = 0.0;
92 vtkIdType nei;
93 for (int sample=0; sample < numPts; ++sample)
94 {
95 nei = pIds->GetId(sample);
96 if ( nei != ptId ) //exclude ourselves
97 {
98 py = this->Points + 3*nei;
99 y[0] = static_cast<double>(*py++);
100 y[1] = static_cast<double>(*py++);
101 y[2] = static_cast<double>(*py);
102 sum += sqrt( vtkMath::Distance2BetweenPoints(x,y) );
103 }
104 }//sum the lengths of all samples exclusing current point
105
106 // Average the lengths; again exclude ourselves
107 if ( numPts > 0 )
108 {
109 this->Distance[ptId] = sum / static_cast<double>(numPts-1);
110 threadMean += this->Distance[ptId];
111 threadCount++;
112 }
113 else //ignore if no points are found, something bad has happened
114 {
115 this->Distance[ptId] = VTK_FLOAT_MAX; //the effect is to eliminate it
116 }
117 }
118 }
119
120 // Compute the mean by compositing all threads
Reduce__anon3454176f0111::ComputeMeanDistance121 void Reduce()
122 {
123 double mean=0.0;
124 vtkIdType count=0;
125
126 vtkSMPThreadLocal<double>::iterator mItr;
127 vtkSMPThreadLocal<double>::iterator mEnd = this->ThreadMean.end();
128 for ( mItr=this->ThreadMean.begin(); mItr != mEnd; ++mItr )
129 {
130 mean += *mItr;
131 }
132
133 vtkSMPThreadLocal<vtkIdType>::iterator cItr;
134 vtkSMPThreadLocal<vtkIdType>::iterator cEnd = this->ThreadCount.end();
135 for ( cItr=this->ThreadCount.begin(); cItr != cEnd; ++cItr )
136 {
137 count += *cItr;
138 }
139
140 count = ( count < 1 ? 1 : count);
141 this->Mean = mean / static_cast<double>(count);
142 }
143
Execute__anon3454176f0111::ComputeMeanDistance144 static void Execute(vtkStatisticalOutlierRemoval *self, vtkIdType numPts,
145 T *points, float *distances, double& mean)
146 {
147 ComputeMeanDistance compute(points, self->GetLocator(),
148 self->GetSampleSize(), distances);
149 vtkSMPTools::For(0, numPts, compute);
150 mean = compute.Mean;
151 }
152
153 }; //ComputeMeanDistance
154
155 //----------------------------------------------------------------------------
156 // Now that the mean is known, compute the standard deviation
157 struct ComputeStdDev
158 {
159 float *Distances;
160 double Mean;
161 double StdDev;
162 vtkSMPThreadLocal<double> ThreadSigma;
163 vtkSMPThreadLocal<vtkIdType> ThreadCount;
164
ComputeStdDev__anon3454176f0111::ComputeStdDev165 ComputeStdDev(float *d, double mean) : Distances(d), Mean(mean), StdDev(0.0)
166 {
167 }
168
Initialize__anon3454176f0111::ComputeStdDev169 void Initialize()
170 {
171 double &threadSigma = this->ThreadSigma.Local();
172 threadSigma = 0.0;
173
174 vtkIdType &threadCount = this->ThreadCount.Local();
175 threadCount = 0;
176 }
177
operator ()__anon3454176f0111::ComputeStdDev178 void operator() (vtkIdType ptId, vtkIdType endPtId)
179 {
180 double &threadSigma = this->ThreadSigma.Local();
181 vtkIdType &threadCount = this->ThreadCount.Local();
182 float d;
183
184 for ( ; ptId < endPtId; ++ptId)
185 {
186 d = this->Distances[ptId];
187 if ( d < VTK_FLOAT_MAX )
188 {
189 threadSigma += (this->Mean - d) * (this->Mean - d);
190 threadCount++;
191 }
192 else
193 {
194 continue; //skip bad point
195 }
196 }
197 }
198
Reduce__anon3454176f0111::ComputeStdDev199 void Reduce()
200 {
201 double sigma=0.0;
202 vtkIdType count=0;
203
204 vtkSMPThreadLocal<double>::iterator sItr;
205 vtkSMPThreadLocal<double>::iterator sEnd = this->ThreadSigma.end();
206 for ( sItr=this->ThreadSigma.begin(); sItr != sEnd; ++sItr )
207 {
208 sigma += *sItr;
209 }
210
211 vtkSMPThreadLocal<vtkIdType>::iterator cItr;
212 vtkSMPThreadLocal<vtkIdType>::iterator cEnd = this->ThreadCount.end();
213 for ( cItr=this->ThreadCount.begin(); cItr != cEnd; ++cItr )
214 {
215 count += *cItr;
216 }
217
218 this->StdDev = sqrt( sigma / static_cast<double>(count) );
219 }
220
Execute__anon3454176f0111::ComputeStdDev221 static void Execute(vtkIdType numPts, float *distances,
222 double mean, double& sigma)
223 {
224 ComputeStdDev stdDev(distances, mean);
225 vtkSMPTools::For(0, numPts, stdDev);
226 sigma = stdDev.StdDev;
227 }
228
229 }; //ComputeStdDev
230
231 //----------------------------------------------------------------------------
232 // Statistics are computed, now filter the points
233 struct RemoveOutliers
234 {
235 double Mean;
236 double Sigma;
237 float *Distances;
238 vtkIdType *PointMap;
239
RemoveOutliers__anon3454176f0111::RemoveOutliers240 RemoveOutliers(double mean, double sigma, float *distances, vtkIdType *map) :
241 Mean(mean), Sigma(sigma), Distances(distances), PointMap(map)
242 {
243 }
244
operator ()__anon3454176f0111::RemoveOutliers245 void operator() (vtkIdType ptId, vtkIdType endPtId)
246 {
247 vtkIdType *map = this->PointMap + ptId;
248 float *d = this->Distances + ptId;
249 double mean=this->Mean, sigma=this->Sigma;
250
251 for ( ; ptId < endPtId; ++ptId)
252 {
253 *map++ = ( fabs(*d++ - mean) <= sigma ? 1 : -1 );
254 }
255 }
256
Execute__anon3454176f0111::RemoveOutliers257 static void Execute(vtkIdType numPts, float *distances, double mean,
258 double sigma, vtkIdType *map)
259 {
260 RemoveOutliers remove(mean, sigma, distances, map);
261 vtkSMPTools::For(0, numPts, remove);
262 }
263
264 }; //RemoveOutliers
265
266
267 } //anonymous namespace
268
269
270 //================= Begin class proper =======================================
271 //----------------------------------------------------------------------------
vtkStatisticalOutlierRemoval()272 vtkStatisticalOutlierRemoval::vtkStatisticalOutlierRemoval()
273 {
274 this->SampleSize = 25;
275 this->StandardDeviationFactor = 1.0;
276 this->Locator = vtkStaticPointLocator::New();
277
278 this->ComputedMean = 0.0;
279 this->ComputedStandardDeviation = 0.0;
280
281 }
282
283 //----------------------------------------------------------------------------
~vtkStatisticalOutlierRemoval()284 vtkStatisticalOutlierRemoval::~vtkStatisticalOutlierRemoval()
285 {
286 this->SetLocator(nullptr);
287 }
288
289 //----------------------------------------------------------------------------
290 // Traverse all the input points and gather statistics about average distance
291 // between them, and the standard deviation of variation. Then filter points
292 // within a specified deviation from the mean.
FilterPoints(vtkPointSet * input)293 int vtkStatisticalOutlierRemoval::FilterPoints(vtkPointSet *input)
294 {
295 // Perform the point removal
296 // Start by building the locator
297 if ( !this->Locator )
298 {
299 vtkErrorMacro(<<"Point locator required\n");
300 return 0;
301 }
302 this->Locator->SetDataSet(input);
303 this->Locator->BuildLocator();
304
305 // Compute statistics across the point cloud. Start my computing
306 // mean distance to N closest neighbors.
307 vtkIdType numPts = input->GetNumberOfPoints();
308 float *dist = new float [numPts];
309 void *inPtr = input->GetPoints()->GetVoidPointer(0);
310 double mean=0.0, sigma=0.0;
311 switch (input->GetPoints()->GetDataType())
312 {
313 vtkTemplateMacro(ComputeMeanDistance<VTK_TT>::
314 Execute(this, numPts, (VTK_TT *)inPtr, dist, mean));
315 }
316
317 // At this point the mean distance for each point, and across the point
318 // cloud is known. Now compute global standard deviation.
319 ComputeStdDev::Execute(numPts, dist, mean, sigma);
320
321 // Finally filter the points based on specified deviation range.
322 RemoveOutliers::Execute(numPts, dist, mean,
323 this->StandardDeviationFactor*sigma, this->PointMap);
324
325 // Assign derived variable
326 this->ComputedMean = mean;
327 this->ComputedStandardDeviation = sigma;
328
329 // Clean up
330 delete [] dist;
331
332 return 1;
333 }
334
335 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)336 void vtkStatisticalOutlierRemoval::PrintSelf(ostream& os, vtkIndent indent)
337 {
338 this->Superclass::PrintSelf(os,indent);
339
340 os << indent << "Sample Size: " << this->SampleSize << "\n";
341 os << indent << "Standard Deviation Factor: "
342 << this->StandardDeviationFactor << "\n";
343 os << indent << "Locator: " << this->Locator << "\n";
344
345 os << indent << "Computed Mean: " << this->ComputedMean << "\n";
346 os << indent << "Computed Standard Deviation: "
347 << this->ComputedStandardDeviation << "\n";
348 }
349