1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkSelectEnclosedPoints.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 "vtkSelectEnclosedPoints.h"
16
17 #include "vtkDataSet.h"
18 #include "vtkIdList.h"
19 #include "vtkIdTypeArray.h"
20 #include "vtkInformation.h"
21 #include "vtkInformationVector.h"
22 #include "vtkObjectFactory.h"
23 #include "vtkPointData.h"
24 #include "vtkCellData.h"
25 #include "vtkPolyData.h"
26 #include "vtkUnsignedCharArray.h"
27 #include "vtkExecutive.h"
28 #include "vtkFeatureEdges.h"
29 #include "vtkStaticCellLocator.h"
30 #include "vtkGenericCell.h"
31 #include "vtkMath.h"
32 #include "vtkGarbageCollector.h"
33 #include "vtkSMPTools.h"
34 #include "vtkSMPThreadLocal.h"
35 #include "vtkSMPThreadLocalObject.h"
36 #include "vtkRandomPool.h"
37
38 vtkStandardNewMacro(vtkSelectEnclosedPoints);
39
40 //----------------------------------------------------------------------------
41 // Classes support threading. Each point can be processed separately, so the
42 // in/out containment check is threaded.
43 namespace {
44
45 //----------------------------------------------------------------------------
46 // The threaded core of the algorithm. Thread on point type.
47 struct SelectInOutCheck
48 {
49 vtkIdType NumPts;
50 vtkDataSet *DataSet;
51 vtkPolyData *Surface;
52 double Bounds[6];
53 double Length;
54 double Tolerance;
55 vtkStaticCellLocator *Locator;
56 unsigned char *Hits;
57 vtkSelectEnclosedPoints *Selector;
58 vtkTypeBool InsideOut;
59 vtkRandomPool *Sequence;
60 vtkSMPThreadLocal<vtkIntersectionCounter> Counter;
61
62 // Don't want to allocate working arrays on every thread invocation. Thread local
63 // storage eliminates lots of new/delete.
64 vtkSMPThreadLocalObject<vtkIdList> CellIds;
65 vtkSMPThreadLocalObject<vtkGenericCell> Cell;
66
SelectInOutCheck__anonacd90eac0111::SelectInOutCheck67 SelectInOutCheck(vtkIdType numPts, vtkDataSet *ds, vtkPolyData *surface, double bds[6], double tol,
68 vtkStaticCellLocator *loc, unsigned char *hits, vtkSelectEnclosedPoints *sel,
69 vtkTypeBool io) : NumPts(numPts), DataSet(ds), Surface(surface), Tolerance(tol),
70 Locator(loc), Hits(hits), Selector(sel), InsideOut(io)
71 {
72 this->Bounds[0] = bds[0];
73 this->Bounds[1] = bds[1];
74 this->Bounds[2] = bds[2];
75 this->Bounds[3] = bds[3];
76 this->Bounds[4] = bds[4];
77 this->Bounds[5] = bds[5];
78 this->Length = sqrt( (bds[1]-bds[0])*(bds[1]-bds[0]) + (bds[3]-bds[2])*(bds[3]-bds[2]) +
79 (bds[5]-bds[4])*(bds[5]-bds[4]) );
80
81 // Precompute a sufficiently large enough random sequence
82 this->Sequence = vtkRandomPool::New();
83 this->Sequence->SetSize((numPts > 1500 ? numPts : 1500));
84 this->Sequence->GeneratePool();
85 }
86
~SelectInOutCheck__anonacd90eac0111::SelectInOutCheck87 ~SelectInOutCheck()
88 {
89 this->Sequence->Delete();
90 }
91
Initialize__anonacd90eac0111::SelectInOutCheck92 void Initialize()
93 {
94 vtkIdList*& cellIds = this->CellIds.Local();
95 cellIds->Allocate(512);
96 vtkIntersectionCounter& counter = this->Counter.Local();
97 counter.SetTolerance(this->Tolerance);
98 }
99
operator ()__anonacd90eac0111::SelectInOutCheck100 void operator() (vtkIdType ptId, vtkIdType endPtId)
101 {
102 double x[3];
103 unsigned char *hits = this->Hits + ptId;
104 vtkGenericCell*& cell = this->Cell.Local();
105 vtkIdList*& cellIds = this->CellIds.Local();
106 vtkIntersectionCounter& counter = this->Counter.Local();
107
108 for ( ; ptId < endPtId; ++ptId )
109 {
110 this->DataSet->GetPoint(ptId, x);
111
112 if ( this->Selector->IsInsideSurface(x, this->Surface, this->Bounds, this->Length,
113 this->Tolerance, this->Locator, cellIds, cell,
114 counter, this->Sequence, ptId) )
115 {
116 *hits++ = (this->InsideOut ? 0 : 1);
117 }
118 else
119 {
120 *hits++ = (this->InsideOut ? 1 : 0);
121 }
122 }
123 }
124
Reduce__anonacd90eac0111::SelectInOutCheck125 void Reduce()
126 {
127 }
128
Execute__anonacd90eac0111::SelectInOutCheck129 static void Execute(vtkIdType numPts, vtkDataSet *ds, vtkPolyData *surface,
130 double bds[6], double tol, vtkStaticCellLocator *loc,
131 unsigned char *hits, vtkSelectEnclosedPoints *sel)
132 {
133 SelectInOutCheck inOut(numPts, ds, surface, bds, tol, loc, hits, sel,
134 sel->GetInsideOut());
135 vtkSMPTools::For(0, numPts, inOut);
136 }
137 }; //SelectInOutCheck
138
139 } //anonymous namespace
140
141
142 //----------------------------------------------------------------------------
143 // Construct object.
vtkSelectEnclosedPoints()144 vtkSelectEnclosedPoints::vtkSelectEnclosedPoints()
145 {
146 this->SetNumberOfInputPorts(2);
147
148 this->CheckSurface = false;
149 this->InsideOut = 0;
150 this->Tolerance = 0.0001;
151
152 this->InsideOutsideArray = nullptr;
153
154 // These are needed to support backward compatibility
155 this->CellLocator = vtkStaticCellLocator::New();
156 this->CellIds = vtkIdList::New();
157 this->Cell = vtkGenericCell::New();
158 }
159
160 //----------------------------------------------------------------------------
~vtkSelectEnclosedPoints()161 vtkSelectEnclosedPoints::~vtkSelectEnclosedPoints()
162 {
163 if ( this->InsideOutsideArray )
164 {
165 this->InsideOutsideArray->Delete();
166 }
167
168 if ( this->CellLocator )
169 {
170 vtkAbstractCellLocator *loc = this->CellLocator;
171 this->CellLocator = nullptr;
172 loc->Delete();
173 }
174
175 this->CellIds->Delete();
176 this->Cell->Delete();
177 }
178
179 //----------------------------------------------------------------------------
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)180 int vtkSelectEnclosedPoints::RequestData(
181 vtkInformation *vtkNotUsed(request),
182 vtkInformationVector **inputVector,
183 vtkInformationVector *outputVector)
184 {
185 // get the info objects
186 vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
187 vtkInformation *in2Info = inputVector[1]->GetInformationObject(0);
188 vtkInformation *outInfo = outputVector->GetInformationObject(0);
189
190 // get the two inputs and output
191 vtkDataSet *input = vtkDataSet::SafeDownCast(
192 inInfo->Get(vtkDataObject::DATA_OBJECT()));
193 vtkPolyData *surface = vtkPolyData::SafeDownCast(
194 in2Info->Get(vtkDataObject::DATA_OBJECT()));
195 vtkDataSet *output = vtkDataSet::SafeDownCast(
196 outInfo->Get(vtkDataObject::DATA_OBJECT()));
197
198 vtkDebugMacro("Selecting enclosed points");
199
200 // If requested, check that the surface is closed
201 if ( this->CheckSurface && ! this->IsSurfaceClosed(surface) )
202 {
203 return 0;
204 }
205
206 // Initiailize search structures
207 this->Initialize(surface);
208
209 // Create array to mark inside/outside
210 if ( this->InsideOutsideArray )
211 {
212 this->InsideOutsideArray->Delete();
213 }
214 this->InsideOutsideArray = vtkUnsignedCharArray::New();
215 vtkUnsignedCharArray *hits = this->InsideOutsideArray;
216 hits->SetName("SelectedPointsArray");
217
218 // Loop over all input points determining inside/outside
219 vtkIdType numPts = input->GetNumberOfPoints();
220 hits->SetNumberOfValues(numPts);
221 unsigned char *hitsPtr = static_cast<unsigned char *>(hits->GetVoidPointer(0));
222
223 // Process the points in parallel
224 SelectInOutCheck::Execute(numPts, input, surface, this->Bounds, this->Tolerance,
225 this->CellLocator, hitsPtr, this);
226
227 // Copy all the input geometry and data to the output.
228 output->CopyStructure(input);
229 output->GetPointData()->PassData(input->GetPointData());
230 output->GetCellData()->PassData(input->GetCellData());
231
232 // Add the new scalars array to the output.
233 hits->SetName("SelectedPoints");
234 output->GetPointData()->SetScalars(hits);
235
236 // release memory
237 this->Complete();
238
239 return 1;
240 }
241
242 //----------------------------------------------------------------------------
243 int vtkSelectEnclosedPoints::
IsSurfaceClosed(vtkPolyData * surface)244 IsSurfaceClosed(vtkPolyData *surface)
245 {
246 vtkPolyData *checker = vtkPolyData::New();
247 checker->CopyStructure(surface);
248
249 vtkFeatureEdges *features = vtkFeatureEdges::New();
250 features->SetInputData(checker);
251 features->BoundaryEdgesOn();
252 features->NonManifoldEdgesOn();
253 features->ManifoldEdgesOff();
254 features->FeatureEdgesOff();
255 features->Update();
256
257 vtkIdType numCells = features->GetOutput()->GetNumberOfCells();
258 features->Delete();
259 checker->Delete();
260
261 if ( numCells > 0 )
262 {
263 return 0;
264 }
265 else
266 {
267 return 1;
268 }
269 }
270
271 //----------------------------------------------------------------------------
Initialize(vtkPolyData * surface)272 void vtkSelectEnclosedPoints::Initialize(vtkPolyData *surface)
273 {
274 if ( ! this->CellLocator )
275 {
276 this->CellLocator = vtkStaticCellLocator::New();
277 }
278
279 this->Surface = surface;
280 surface->GetBounds(this->Bounds);
281 this->Length = surface->GetLength();
282
283 // Set up structures for acceleration ray casting
284 this->CellLocator->SetDataSet(surface);
285 this->CellLocator->BuildLocator();
286 }
287
288 //----------------------------------------------------------------------------
IsInside(vtkIdType inputPtId)289 int vtkSelectEnclosedPoints::IsInside(vtkIdType inputPtId)
290 {
291 if ( !this->InsideOutsideArray ||
292 this->InsideOutsideArray->GetValue(inputPtId) == 0 )
293 {
294 return 0;
295 }
296 else
297 {
298 return 1;
299 }
300 }
301
302 //----------------------------------------------------------------------------
IsInsideSurface(double x,double y,double z)303 int vtkSelectEnclosedPoints::IsInsideSurface(double x, double y, double z)
304 {
305 double xyz[3];
306 xyz[0] = x;
307 xyz[1] = y;
308 xyz[2] = z;
309 return this->IsInsideSurface(xyz);
310 }
311
312 //----------------------------------------------------------------------------
313 // This is done to preserve backward compatibility. However it is not thread
314 // safe due to the use of the data member CellIds and Cell.
IsInsideSurface(double x[3])315 int vtkSelectEnclosedPoints::IsInsideSurface(double x[3])
316 {
317 vtkIntersectionCounter counter(this->Tolerance, this->Length);
318
319 return this->IsInsideSurface(x, this->Surface, this->Bounds, this->Length,
320 this->Tolerance, this->CellLocator, this->CellIds,
321 this->Cell, counter);
322 }
323
324 //----------------------------------------------------------------------------
325 // General method uses ray casting to determine in/out. Since this is a
326 // numerically delicate operation, we use a crude "statistical" method (based
327 // on voting) to provide a better answer. Plus there is a process to merge
328 // nearly conincident points along the intersection rays.
329 //
330 // This is a static method so it can be used by other filters; hence the
331 // many parameters used.
332 //
333 // Provision for reproducible threaded random number generation is made by
334 // supporting the precomputation of a random sequence (see vtkRandomPool).
335 //
336 #define VTK_MAX_ITER 10 //Maximum iterations for ray-firing
337 #define VTK_VOTE_THRESHOLD 2 //Vote margin for test
338
339 int vtkSelectEnclosedPoints::
IsInsideSurface(double x[3],vtkPolyData * surface,double bds[6],double length,double tolerance,vtkAbstractCellLocator * locator,vtkIdList * cellIds,vtkGenericCell * genCell,vtkIntersectionCounter & counter,vtkRandomPool * seq,vtkIdType seqIdx)340 IsInsideSurface(double x[3], vtkPolyData *surface, double bds[6],
341 double length, double tolerance,
342 vtkAbstractCellLocator *locator, vtkIdList *cellIds,
343 vtkGenericCell *genCell, vtkIntersectionCounter &counter,
344 vtkRandomPool* seq, vtkIdType seqIdx)
345 {
346 // do a quick inside bounds check against the surface bounds
347 if ( x[0] < bds[0] || x[0] > bds[1] ||
348 x[1] < bds[2] || x[1] > bds[3] ||
349 x[2] < bds[4] || x[2] > bds[5])
350 {
351 return 0;
352 }
353
354 // Shortly we are going to start firing rays. It's important that the rays
355 // are long enough to go from the test point all the way through the
356 // enclosing surface. So compute a vector from the test point to the center
357 // of the surface, and then add in the length (diagonal of bounding box) of
358 // the surface.
359 double offset[3], totalLength;
360 offset[0] = x[0] - ((bds[0]+bds[1]) / 2.0);
361 offset[1] = x[1] - ((bds[2]+bds[3]) / 2.0);
362 offset[2] = x[2] - ((bds[4]+bds[5]) / 2.0);
363 totalLength = length + vtkMath::Norm(offset);
364
365 // Perform in/out by shooting random rays. Multiple rays are fired
366 // to improve accuracy of the result.
367 //
368 // The variable iterNumber counts the number of rays fired and is
369 // limited by the defined variable VTK_MAX_ITER.
370 //
371 // The variable deltaVotes keeps track of the number of votes for
372 // "in" versus "out" of the surface. When deltaVotes > 0, more votes
373 // have counted for "in" than "out". When deltaVotes < 0, more votes
374 // have counted for "out" than "in". When the delta_vote exceeds or
375 // equals the defined variable VTK_VOTE_THRESHOLD, then the
376 // appropriate "in" or "out" status is returned.
377 //
378 double rayMag, ray[3], xray[3], t, pcoords[3], xint[3];
379 int i, numInts, iterNumber, deltaVotes, subId;
380 vtkIdType idx, numCells;
381 double tol = tolerance * length;
382
383 for (deltaVotes = 0, iterNumber = 1;
384 (iterNumber < VTK_MAX_ITER) && (abs(deltaVotes) < VTK_VOTE_THRESHOLD);
385 iterNumber++)
386 {
387 // Define a random ray to fire.
388 rayMag = 0.0;
389 while ( rayMag == 0.0 )
390 {
391 if ( seq == nullptr ) //in serial mode
392 {
393 ray[0] = vtkMath::Random(-1.0,1.0);
394 ray[1] = vtkMath::Random(-1.0,1.0);
395 ray[2] = vtkMath::Random(-1.0,1.0);
396 }
397 else //threading, have to scale sequence -1<=x<=1
398 {
399 ray[0] = 2.0 *(0.5 - seq->GetValue(seqIdx++));
400 ray[1] = 2.0 *(0.5 - seq->GetValue(seqIdx++));
401 ray[2] = 2.0 *(0.5 - seq->GetValue(seqIdx++));
402 }
403 rayMag = vtkMath::Norm(ray);
404 }
405
406 // The ray must be appropriately sized wrt the bounding box. (It has to
407 // go all the way through the bounding box. Remember though that an
408 // "inside bounds" check was done previously so diagonal length should
409 // be long enough.)
410 for (i=0; i<3; i++)
411 {
412 xray[i] = x[i] + 2.0*totalLength*(ray[i]/rayMag);
413 }
414
415 // Retrieve the candidate cells from the locator to limit the
416 // intersections to be attempted.
417 locator->FindCellsAlongLine(x,xray,tol,cellIds);
418 numCells = cellIds->GetNumberOfIds();
419
420 counter.Reset();
421 for ( idx=0; idx < numCells; idx++ )
422 {
423 surface->GetCell(cellIds->GetId(idx), genCell);
424 if ( genCell->IntersectWithLine(x, xray, tol, t, xint, pcoords, subId) )
425 {
426 counter.AddIntersection(t);
427 }
428 } //for all candidate cells along this ray
429
430 numInts = counter.CountIntersections();
431
432 if ( (numInts % 2) == 0) //if outside
433 {
434 --deltaVotes;
435 }
436 else //if inside
437 {
438 ++deltaVotes;
439 }
440 } //try another ray
441
442 // If the number of votes is positive, the point is inside
443 //
444 return ( deltaVotes < 0 ? 0 : 1 );
445 }
446
447 #undef VTK_MAX_ITER
448 #undef VTK_VOTE_THRESHOLD
449
450
451 //----------------------------------------------------------------------------
452 // Specify the second enclosing surface input via a connection
453 void vtkSelectEnclosedPoints::
SetSurfaceConnection(vtkAlgorithmOutput * algOutput)454 SetSurfaceConnection(vtkAlgorithmOutput* algOutput)
455 {
456 this->SetInputConnection(1, algOutput);
457 }
458
459 //----------------------------------------------------------------------------
460 // Specify the second enclosing surface input data
SetSurfaceData(vtkPolyData * pd)461 void vtkSelectEnclosedPoints::SetSurfaceData(vtkPolyData *pd)
462 {
463 this->SetInputData(1, pd);
464 }
465
466 //----------------------------------------------------------------------------
467 // Return the enclosing surface
GetSurface()468 vtkPolyData *vtkSelectEnclosedPoints::GetSurface()
469 {
470 return vtkPolyData::SafeDownCast(
471 this->GetExecutive()->GetInputData(1, 0));
472 }
473
474 //----------------------------------------------------------------------------
GetSurface(vtkInformationVector * sourceInfo)475 vtkPolyData* vtkSelectEnclosedPoints::GetSurface(vtkInformationVector *sourceInfo)
476 {
477 vtkInformation *info = sourceInfo->GetInformationObject(1);
478 if (!info)
479 {
480 return nullptr;
481 }
482 return vtkPolyData::SafeDownCast(info->Get(vtkDataObject::DATA_OBJECT()));
483 }
484
485 //----------------------------------------------------------------------------
FillInputPortInformation(int port,vtkInformation * info)486 int vtkSelectEnclosedPoints::FillInputPortInformation(int port, vtkInformation *info)
487 {
488 if (port == 0)
489 {
490 info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
491 }
492 else if (port == 1)
493 {
494 info->Set(vtkAlgorithm::INPUT_IS_REPEATABLE(), 0);
495 info->Set(vtkAlgorithm::INPUT_IS_OPTIONAL(), 0);
496 info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkPolyData");
497 }
498
499 return 1;
500 }
501
502 //----------------------------------------------------------------------------
Complete()503 void vtkSelectEnclosedPoints::Complete()
504 {
505 this->CellLocator->FreeSearchStructure();
506 }
507
508 //----------------------------------------------------------------------------
ReportReferences(vtkGarbageCollector * collector)509 void vtkSelectEnclosedPoints::ReportReferences(vtkGarbageCollector* collector)
510 {
511 this->Superclass::ReportReferences(collector);
512 // These filters share our input and are therefore involved in a
513 // reference loop.
514 vtkGarbageCollectorReport(collector, this->CellLocator, "CellLocator");
515 }
516
517 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)518 void vtkSelectEnclosedPoints::PrintSelf(ostream& os, vtkIndent indent)
519 {
520 this->Superclass::PrintSelf(os,indent);
521
522 os << indent << "Check Surface: "
523 << (this->CheckSurface ? "On\n" : "Off\n");
524
525 os << indent << "Inside Out: "
526 << (this->InsideOut ? "On\n" : "Off\n");
527
528 os << indent << "Tolerance: " << this->Tolerance << "\n";
529 }
530