vtkProbeFilter.cxx 28 KB
Newer Older
Will Schroeder's avatar
Will Schroeder committed
1
2
/*=========================================================================

Ken Martin's avatar
Ken Martin committed
3
  Program:   Visualization Toolkit
Ken Martin's avatar
Ken Martin committed
4
  Module:    vtkProbeFilter.cxx
Will Schroeder's avatar
Will Schroeder committed
5

6
  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7
8
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
Ken Martin's avatar
Ken Martin committed
9

10
11
     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12
     PURPOSE.  See the above copyright notice for more information.
Will Schroeder's avatar
Will Schroeder committed
13
14

=========================================================================*/
Ken Martin's avatar
Ken Martin committed
15
#include "vtkProbeFilter.h"
16

17
#include "vtkBoundingBox.h"
Utkarsh Ayachit's avatar
Utkarsh Ayachit committed
18
#include "vtkCell.h"
19
#include "vtkCellData.h"
Utkarsh Ayachit's avatar
Utkarsh Ayachit committed
20
#include "vtkCharArray.h"
Sujin Philip's avatar
Sujin Philip committed
21
#include "vtkGenericCell.h"
22
#include "vtkIdTypeArray.h"
23
#include "vtkImageData.h"
24
25
#include "vtkInformation.h"
#include "vtkInformationVector.h"
26
#include "vtkMath.h"
27
28
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
29
#include "vtkPointSet.h"
Sujin Philip's avatar
Sujin Philip committed
30
#include "vtkSmartPointer.h"
Sujin Philip's avatar
Sujin Philip committed
31
32
#include "vtkSMPTools.h"
#include "vtkSMPThreadLocal.h"
Sujin Philip's avatar
Sujin Philip committed
33
#include "vtkStaticCellLocator.h"
34
#include "vtkStreamingDemandDrivenPipeline.h"
35

Sujin Philip's avatar
Sujin Philip committed
36
#include <algorithm>
37
#include <vector>
38

Brad King's avatar
Brad King committed
39
vtkStandardNewMacro(vtkProbeFilter);
Sujin Philip's avatar
Sujin Philip committed
40
vtkCxxSetObjectMacro(vtkProbeFilter, CellLocatorPrototype, vtkAbstractCellLocator);
41

42
43
#define CELL_TOLERANCE_FACTOR_SQR  1e-6

44
class vtkProbeFilter::vtkVectorOfArrays :
45
  public std::vector<vtkDataArray*>
46
47
48
{
};

49
//----------------------------------------------------------------------------
Ken Martin's avatar
Ken Martin committed
50
vtkProbeFilter::vtkProbeFilter()
51
{
52
  this->CategoricalData = 0;
53
  this->SpatialMatch = 0;
54
  this->ValidPoints = vtkIdTypeArray::New();
55
  this->MaskPoints = nullptr;
56
  this->SetNumberOfInputPorts(2);
57
  this->ValidPointMaskArrayName = nullptr;
Utkarsh Ayachit's avatar
Utkarsh Ayachit committed
58
  this->SetValidPointMaskArrayName("vtkValidPointMask");
59
  this->CellArrays = new vtkVectorOfArrays();
Sujin Philip's avatar
Sujin Philip committed
60
  this->CellLocatorPrototype = nullptr;
61

62
63
  this->PointList = nullptr;
  this->CellList = nullptr;
64

Burlen Loring's avatar
Burlen Loring committed
65
66
  this->PassCellArrays = 0;
  this->PassPointArrays = 0;
67
  this->PassFieldArrays = 1;
68
69
  this->Tolerance = 1.0;
  this->ComputeTolerance = 1;
70
71
}

72
//----------------------------------------------------------------------------
73
74
vtkProbeFilter::~vtkProbeFilter()
{
75
  if (this->MaskPoints)
76
  {
77
    this->MaskPoints->Delete();
78
  }
79
  this->ValidPoints->Delete();
80

Sujin Philip's avatar
Sujin Philip committed
81
82
83
84
  this->vtkProbeFilter::SetValidPointMaskArrayName(nullptr);
  this->vtkProbeFilter::SetCellLocatorPrototype(nullptr);

  delete this->CellArrays;
85
86
  delete this->PointList;
  delete this->CellList;
87
88
}

89
90
91
92
93
//----------------------------------------------------------------------------
void vtkProbeFilter::SetSourceConnection(vtkAlgorithmOutput* algOutput)
{
  this->SetInputConnection(1, algOutput);
}
94

95
//----------------------------------------------------------------------------
96
void vtkProbeFilter::SetSourceData(vtkDataObject *input)
97
{
98
  this->SetInputData(1, input);
99
100
101
}

//----------------------------------------------------------------------------
102
vtkDataObject *vtkProbeFilter::GetSource()
103
{
104
  if (this->GetNumberOfInputConnections(1) < 1)
105
  {
106
    return nullptr;
107
  }
108

109
  return this->GetExecutive()->GetInputData(1, 0);
110
111
}

Sujin Philip's avatar
Sujin Philip committed
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
//----------------------------------------------------------------------------
vtkIdTypeArray* vtkProbeFilter::GetValidPoints()
{
  if (this->MaskPoints &&
      this->MaskPoints->GetMTime() > this->ValidPoints->GetMTime())
  {
    char* maskArray = this->MaskPoints->GetPointer(0);
    vtkIdType numPts = this->MaskPoints->GetNumberOfTuples();
    vtkIdType numValidPoints = std::count(maskArray, maskArray + numPts,
                                          static_cast<char>(1));
    this->ValidPoints->Allocate(numValidPoints);
    for (vtkIdType i = 0; i < numPts; ++i)
    {
      if (maskArray[i])
      {
        this->ValidPoints->InsertNextValue(i);
      }
    }
    this->ValidPoints->Modified();
  }

  return this->ValidPoints;
}

136
//----------------------------------------------------------------------------
137
138
139
140
int vtkProbeFilter::RequestData(
  vtkInformation *vtkNotUsed(request),
  vtkInformationVector **inputVector,
  vtkInformationVector *outputVector)
Will Schroeder's avatar
Will Schroeder committed
141
{
142
143
144
145
146
  // get the info objects
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation *sourceInfo = inputVector[1]->GetInformationObject(0);
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

147
  // get the input and output
148
149
150
151
152
153
154
  vtkDataSet *input = vtkDataSet::SafeDownCast(
    inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkDataSet *source = vtkDataSet::SafeDownCast(
    sourceInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkDataSet *output = vtkDataSet::SafeDownCast(
    outInfo->Get(vtkDataObject::DATA_OBJECT()));

155
156
157
  // First, copy the input to the output as a starting point
  output->CopyStructure(input);

158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
  if (this->CategoricalData == 1)
  {
    // If the categorical data flag is enabled, then a) there must be scalars
    // to treat as categorical data, and b) the scalars must have one component.
    if (!source->GetPointData()->GetScalars())
    {
      vtkErrorMacro(<<"No input scalars!");
      return 1;
    }
    if (source->GetPointData()->GetScalars()->GetNumberOfComponents() != 1)
    {
      vtkErrorMacro(<<"Source scalars have more than one component! Cannot categorize!");
      return 1;
    }

    // Set the scalar to interpolate via nearest neighbor. That way, we won't
    // get any false values (for example, a zone 4 cell appearing on the
    // boundary of zone 3 and zone 5).
    output->GetPointData()->SetCopyAttribute(vtkDataSetAttributes::SCALARS, 2,
                                             vtkDataSetAttributes::INTERPOLATE);
  }

180
  if (source)
181
  {
182
    this->Probe(input, source, output);
183
  }
184

185
186
187
188
189
190
191
192
  this->PassAttributeData(input, source, output);
  return 1;
}

//----------------------------------------------------------------------------
void vtkProbeFilter::PassAttributeData(
  vtkDataSet* input, vtkDataObject* vtkNotUsed(source), vtkDataSet* output)
{
Burlen Loring's avatar
Burlen Loring committed
193
194
  // copy point data arrays
  if (this->PassPointArrays)
195
  {
Burlen Loring's avatar
Burlen Loring committed
196
197
    int numPtArrays = input->GetPointData()->GetNumberOfArrays();
    for (int i=0; i<numPtArrays; ++i)
198
    {
199
200
      vtkDataArray *da = input->GetPointData()->GetArray(i);
      if (!output->GetPointData()->HasArray(da->GetName()))
201
      {
202
        output->GetPointData()->AddArray(da);
Burlen Loring's avatar
Burlen Loring committed
203
      }
204
    }
205
206
207

    // Set active attributes in the output to the active attributes in the input
    for (int i = 0; i < vtkDataSetAttributes::NUM_ATTRIBUTES; ++i)
208
    {
209
      vtkAbstractArray* da = input->GetPointData()->GetAttribute(i);
210
      if (da && da->GetName() && !output->GetPointData()->GetAttribute(i))
211
      {
212
213
        output->GetPointData()->SetAttribute(da, i);
      }
Burlen Loring's avatar
Burlen Loring committed
214
    }
215
  }
Burlen Loring's avatar
Burlen Loring committed
216
217
218

  // copy cell data arrays
  if (this->PassCellArrays)
219
  {
Burlen Loring's avatar
Burlen Loring committed
220
221
    int numCellArrays = input->GetCellData()->GetNumberOfArrays();
    for (int i=0; i<numCellArrays; ++i)
222
    {
223
224
      vtkDataArray *da = input->GetCellData()->GetArray(i);
      if (!output->GetCellData()->HasArray(da->GetName()))
225
      {
226
        output->GetCellData()->AddArray(da);
Burlen Loring's avatar
Burlen Loring committed
227
      }
228
    }
229
230
231

    // Set active attributes in the output to the active attributes in the input
    for (int i = 0; i < vtkDataSetAttributes::NUM_ATTRIBUTES; ++i)
232
    {
233
      vtkAbstractArray* da = input->GetCellData()->GetAttribute(i);
234
      if (da && da->GetName() && !output->GetCellData()->GetAttribute(i))
235
      {
236
237
        output->GetCellData()->SetAttribute(da, i);
      }
Burlen Loring's avatar
Burlen Loring committed
238
    }
239
  }
Burlen Loring's avatar
Burlen Loring committed
240

241
  if (this->PassFieldArrays)
242
  {
243
    // nothing to do, vtkDemandDrivenPipeline takes care of that.
244
  }
245
  else
246
  {
247
    output->GetFieldData()->Initialize();
248
  }
249
250
}

251
252
253
254
255
256
257
258
259
260
261
262
263
//----------------------------------------------------------------------------
void vtkProbeFilter::BuildFieldList(vtkDataSet* source)
{
  delete this->PointList;
  delete this->CellList;

  this->PointList = new vtkDataSetAttributes::FieldList(1);
  this->PointList->InitializeFieldList(source->GetPointData());

  this->CellList = new vtkDataSetAttributes::FieldList(1);
  this->CellList->InitializeFieldList(source->GetCellData());
}

264
265
266
267
//----------------------------------------------------------------------------
// * input -- dataset probed with
// * source -- dataset probed into
// * output - output.
268
void vtkProbeFilter::InitializeForProbing(vtkDataSet* input,
269
270
  vtkDataSet* output)
{
271
  if (!this->PointList || !this->CellList)
272
  {
273
274
    vtkErrorMacro("BuildFieldList() must be called before calling this method.");
    return;
275
  }
276

277
278
  vtkIdType numPts = input->GetNumberOfPoints();

279
280
281
282
  // if this is repeatedly called by the pipeline for a composite mesh,
  // you need a new array for each block
  // (that is you need to reinitialize the object)
  if (this->MaskPoints)
283
  {
284
    this->MaskPoints->Delete();
285
  }
286
287
  this->MaskPoints = vtkCharArray::New();
  this->MaskPoints->SetNumberOfComponents(1);
288
  this->MaskPoints->SetNumberOfTuples(numPts);
Sujin Philip's avatar
Sujin Philip committed
289
  this->MaskPoints->FillValue(0);
290
  this->MaskPoints->SetName(this->ValidPointMaskArrayName?
291
292
293
294
295
    this->ValidPointMaskArrayName: "vtkValidPointMask");

  // Allocate storage for output PointData
  // All input PD is passed to output as PD. Those arrays in input CD that are
  // not present in output PD will be passed as output PD.
296
  vtkPointData* outPD = output->GetPointData();
297
298
299
  outPD->InterpolateAllocate((*this->PointList), numPts, numPts);

  vtkCellData* tempCellData = vtkCellData::New();
300
301
302
303
  // We're okay with copying global ids for cells. we just don't flag them as
  // such.
  tempCellData->CopyAllOn(vtkDataSetAttributes::COPYTUPLE);
  tempCellData->CopyAllocate((*this->CellList), numPts, numPts);
304
305

  this->CellArrays->clear();
306
  int numCellArrays = tempCellData->GetNumberOfArrays();
307
  for (int cc=0; cc < numCellArrays; cc++)
308
  {
309
    vtkDataArray* inArray = tempCellData->GetArray(cc);
310
    if (inArray && inArray->GetName() && !outPD->GetArray(inArray->GetName()))
311
    {
312
313
      outPD->AddArray(inArray);
      this->CellArrays->push_back(inArray);
314
    }
315
  }
316
  tempCellData->Delete();
317

Sujin Philip's avatar
Sujin Philip committed
318
  this->InitializeOutputArrays(outPD, numPts);
319
  outPD->AddArray(this->MaskPoints);
Sujin Philip's avatar
Sujin Philip committed
320
}
321

Sujin Philip's avatar
Sujin Philip committed
322
323
324
325
326
327
328
329
330
331
332
333
//----------------------------------------------------------------------------
void vtkProbeFilter::InitializeOutputArrays(vtkPointData *outPD, vtkIdType numPts)
{
  for (int i = 0; i < outPD->GetNumberOfArrays(); ++i)
  {
    vtkDataArray* da = outPD->GetArray(i);
    if (da)
    {
      da->SetNumberOfTuples(numPts);
      da->Fill(0);
    }
  }
334
335
}

Mathieu Malaterre's avatar
Mathieu Malaterre committed
336
//----------------------------------------------------------------------------
337
338
void vtkProbeFilter::DoProbing(vtkDataSet *input, int srcIdx, vtkDataSet *source,
                               vtkDataSet *output)
339
{
340
341
342
343
344
345
346
  vtkBoundingBox sbox(source->GetBounds());
  vtkBoundingBox ibox(input->GetBounds());
  if (!sbox.Intersects(ibox))
  {
    return;
  }

347
  if (vtkImageData::SafeDownCast(input))
348
  {
349
350
    vtkImageData *inImage = vtkImageData::SafeDownCast(input);
    vtkImageData *outImage = vtkImageData::SafeDownCast(output);
351
    this->ProbePointsImageData(inImage, srcIdx, source, outImage);
352
  }
353
  else
354
  {
355
    this->ProbeEmptyPoints(input, srcIdx, source, output);
356
  }
357
358
359
360
361
362
363
364
365
}

//----------------------------------------------------------------------------
void vtkProbeFilter::Probe(vtkDataSet *input, vtkDataSet *source,
                           vtkDataSet *output)
{
  this->BuildFieldList(source);
  this->InitializeForProbing(input, output);
  this->DoProbing(input, 0, source, output);
366
367
368
}

//----------------------------------------------------------------------------
369
void vtkProbeFilter::ProbeEmptyPoints(vtkDataSet *input,
370
371
  int srcIdx,
  vtkDataSet *source, vtkDataSet *output)
372
{
Amy Squillacote's avatar
Amy Squillacote committed
373
  vtkIdType ptId, numPts;
Ken Martin's avatar
Ken Martin committed
374
  double x[3], tol2;
Ken Martin's avatar
Ken Martin committed
375
  vtkPointData *pd, *outPD;
376
  vtkCellData* cd;
Amy Squillacote's avatar
Amy Squillacote committed
377
  int subId;
Ken Martin's avatar
Ken Martin committed
378
379
  double pcoords[3], *weights;
  double fastweights[256];
Will Schroeder's avatar
Will Schroeder committed
380

Ken Martin's avatar
Ken Martin committed
381
  vtkDebugMacro(<<"Probing data");
Will Schroeder's avatar
Will Schroeder committed
382

383
  pd = source->GetPointData();
384
  cd = source->GetCellData();
Mathieu Malaterre's avatar
Mathieu Malaterre committed
385

386
387
388
  // lets use a stack allocated array if possible for performance reasons
  int mcs = source->GetMaxCellSize();
  if (mcs<=256)
389
  {
390
    weights = fastweights;
391
  }
392
  else
393
  {
Ken Martin's avatar
Ken Martin committed
394
    weights = new double[mcs];
395
  }
396

Ken Martin's avatar
Ken Martin committed
397
  numPts = input->GetNumberOfPoints();
Ken Martin's avatar
Ken Martin committed
398
  outPD = output->GetPointData();
399
400

  char* maskArray = this->MaskPoints->GetPointer(0);
401

402
403
  tol2 = this->ComputeTolerance ? VTK_DOUBLE_MAX :
         (this->Tolerance * this->Tolerance);
404

405
406
407
  // vtkPointSet based datasets do not have an implicit structure to their
  // points. A cell locator performs better here than using the dataset's
  // FindCell function.
Sujin Philip's avatar
Sujin Philip committed
408
  vtkSmartPointer<vtkAbstractCellLocator> cellLocator;
409
410
411
412
413
414
415
416
  if (vtkPointSet::SafeDownCast(source) != nullptr)
  {
    cellLocator.TakeReference(this->CellLocatorPrototype ?
                              this->CellLocatorPrototype->NewInstance() :
                              vtkStaticCellLocator::New());
    cellLocator->SetDataSet(source);
    cellLocator->Update();
  }
Sujin Philip's avatar
Sujin Philip committed
417

418
419
  // Loop over all input points, interpolating source data
  //
Sujin Philip's avatar
Sujin Philip committed
420
  vtkNew<vtkGenericCell> gcell;
421
  int abort=0;
Amy Squillacote's avatar
Amy Squillacote committed
422
  vtkIdType progressInterval=numPts/20 + 1;
423
  for (ptId=0; ptId < numPts && !abort; ptId++)
424
  {
425
    if ( !(ptId % progressInterval) )
426
    {
Francois Bertel's avatar
Francois Bertel committed
427
      this->UpdateProgress(static_cast<double>(ptId)/numPts);
428
      abort = GetAbortExecute();
429
    }
430

431
    if (maskArray[ptId] == static_cast<char>(1))
432
    {
433
      // skip points which have already been probed with success.
434
      // This is helpful for multiblock dataset probing.
435
      continue;
436
    }
437

438
    // Get the xyz coordinate of the point in the input dataset
439
    input->GetPoint(ptId, x);
440
441

    // Find the cell that contains xyz and get it
442
443
444
    vtkIdType cellId = (cellLocator.Get() != nullptr) ?
      cellLocator->FindCell(x, tol2, gcell.GetPointer(), pcoords, weights) :
      source->FindCell(x, nullptr, -1, tol2, subId, pcoords, weights);
Sujin Philip's avatar
Sujin Philip committed
445
446

    vtkCell* cell = nullptr;
447
    if (cellId >= 0)
448
    {
449
      cell = source->GetCell(cellId);
450
      if (this->ComputeTolerance)
451
      {
452
453
454
455
456
457
        // If ComputeTolerance is set, compute a tolerance proportional to the
        // cell length.
        double dist2;
        double closestPoint[3];
        cell->EvaluatePosition(x, closestPoint, subId, pcoords, dist2, weights);
        if (dist2 > (cell->GetLength2() * CELL_TOLERANCE_FACTOR_SQR))
458
        {
Sujin Philip's avatar
Sujin Philip committed
459
          continue;
460
        }
461
      }
462
    }
463

464
    if (cell)
465
    {
466
      // Interpolate the point data
467
468
469
470
471
      outPD->InterpolatePoint((*this->PointList), pd, srcIdx, ptId,
        cell->PointIds, weights);
      vtkVectorOfArrays::iterator iter;
      for (iter = this->CellArrays->begin(); iter != this->CellArrays->end();
        ++iter)
472
      {
473
474
        vtkDataArray* inArray = cd->GetArray((*iter)->GetName());
        if (inArray)
475
        {
476
          outPD->CopyTuple(inArray, *iter, cellId, ptId);
477
        }
Will Schroeder's avatar
Will Schroeder committed
478
      }
479
480
481
      maskArray[ptId] = static_cast<char>(1);
    }
  }
Utkarsh Ayachit's avatar
Utkarsh Ayachit committed
482

Sujin Philip's avatar
Sujin Philip committed
483
  this->MaskPoints->Modified();
484
  if (mcs>256)
485
  {
486
    delete [] weights;
487
  }
Will Schroeder's avatar
Will Schroeder committed
488
489
}

Sujin Philip's avatar
Sujin Philip committed
490
//---------------------------------------------------------------------------
491
492
493
494
static void GetPointIdsInRange(double rangeMin, double rangeMax, double start,
  double stepsize, int numSteps, int &minid, int &maxid)
{
  if (stepsize == 0)
495
  {
496
497
    minid = maxid = 0;
    return;
498
  }
499
500
501

  minid = vtkMath::Ceil((rangeMin - start)/stepsize);
  if (minid < 0)
502
  {
503
    minid = 0;
504
  }
505
506
507

  maxid = vtkMath::Floor((rangeMax - start)/stepsize);
  if (maxid > numSteps-1)
508
  {
509
    maxid = numSteps-1;
510
  }
511
512
}

Sujin Philip's avatar
Sujin Philip committed
513
514
515
516
517
518
519
520
void vtkProbeFilter::ProbeImagePointsInCell(vtkCell *cell, vtkIdType cellId,
                                            vtkDataSet *source, int srcBlockId,
                                            const double start[3],
                                            const double spacing[3],
                                            const int dim[3],
                                            vtkPointData *outPD,
                                            char *maskArray,
                                            double *wtsBuff)
521
{
Sujin Philip's avatar
Sujin Philip committed
522
523
  vtkPointData *pd = source->GetPointData();
  vtkCellData *cd = source->GetCellData();
524

Sujin Philip's avatar
Sujin Philip committed
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
  // get coordinates of sampling grids
  double cellBounds[6];
  cell->GetBounds(cellBounds);

  int idxBounds[6];
  GetPointIdsInRange(cellBounds[0], cellBounds[1], start[0], spacing[0],
    dim[0], idxBounds[0], idxBounds[1]);
  GetPointIdsInRange(cellBounds[2], cellBounds[3], start[1], spacing[1],
    dim[1], idxBounds[2], idxBounds[3]);
  GetPointIdsInRange(cellBounds[4], cellBounds[5], start[2], spacing[2],
    dim[2], idxBounds[4], idxBounds[5]);

  if ((idxBounds[1] - idxBounds[0]) < 0 ||
      (idxBounds[3] - idxBounds[2]) < 0 ||
      (idxBounds[5] - idxBounds[4]) < 0)
540
  {
Sujin Philip's avatar
Sujin Philip committed
541
    return;
542
  }
Sujin Philip's avatar
Sujin Philip committed
543
544
545
546
547

  double cpbuf[3];
  double dist2 = 0;
  double *closestPoint = cpbuf;
  if (cell->IsA("vtkCell3D"))
548
  {
Sujin Philip's avatar
Sujin Philip committed
549
    // we only care about closest point and its distance for 2D cells
550
    closestPoint = nullptr;
551
  }
552

Sujin Philip's avatar
Sujin Philip committed
553
554
555
556
557
558
559
560
561
562
563
564
  double userTol2 = this->Tolerance * this->Tolerance;
  for (int iz=idxBounds[4]; iz<=idxBounds[5]; iz++)
  {
    double p[3];
    p[2] = start[2] + iz*spacing[2];
    for (int iy=idxBounds[2]; iy<=idxBounds[3]; iy++)
    {
      p[1] = start[1] + iy*spacing[1];
      for (int ix=idxBounds[0]; ix<=idxBounds[1]; ix++)
      {
        // For each grid point within the cell bound, interpolate values
        p[0] = start[0] + ix*spacing[0];
565

Sujin Philip's avatar
Sujin Philip committed
566
567
568
569
        double pcoords[3];
        int subId;
        int inside = cell->EvaluatePosition(p, closestPoint, subId, pcoords,
                                            dist2, wtsBuff);
570

Sujin Philip's avatar
Sujin Philip committed
571
572
573
574
575
        // If ComputeTolerance is set, compute a tolerance proportional to the
        // cell length. Otherwise, use the user specified absolute tolerance.
        double tol2 = this->ComputeTolerance ?
                      (CELL_TOLERANCE_FACTOR_SQR * cell->GetLength2()) :
                      userTol2;
576

Sujin Philip's avatar
Sujin Philip committed
577
578
579
        if ((inside == 1) && (dist2 <= tol2))
        {
          vtkIdType ptId = ix + dim[0]*(iy + dim[1]*iz);
580

Sujin Philip's avatar
Sujin Philip committed
581
582
583
          // Interpolate the point data
          outPD->InterpolatePoint((*this->PointList), pd, srcBlockId, ptId,
                                  cell->PointIds, wtsBuff);
584

Sujin Philip's avatar
Sujin Philip committed
585
586
587
588
589
590
591
592
593
594
595
          // Assign cell data
          vtkVectorOfArrays::iterator iter;
          for (iter = this->CellArrays->begin();
              iter != this->CellArrays->end(); ++iter)
          {
            vtkDataArray* inArray = cd->GetArray((*iter)->GetName());
            if (inArray)
            {
              outPD->CopyTuple(inArray, *iter, cellId, ptId);
            }
          }
596

Sujin Philip's avatar
Sujin Philip committed
597
598
599
600
601
602
          maskArray[ptId] = static_cast<char>(1);
        }
      }
    }
  }
}
603

Sujin Philip's avatar
Sujin Philip committed
604
namespace {
605

Sujin Philip's avatar
Sujin Philip committed
606
607
608
609
610
class CellStorage
{
public:
  CellStorage()
  {
Sujin Philip's avatar
Sujin Philip committed
611
    this->Initialize();
Sujin Philip's avatar
Sujin Philip committed
612
  }
613

Sujin Philip's avatar
Sujin Philip committed
614
615
  ~CellStorage()
  {
Sujin Philip's avatar
Sujin Philip committed
616
617
618
619
620
621
622
623
624
625
    this->Clear();
  }

  // Copying does not make sense for this class but vtkSMPThreadLocal needs
  // these functions to compile. Just initialize the object.
  CellStorage(const CellStorage&)
  {
    this->Initialize();
  }

626
  CellStorage& operator=(const CellStorage&)
Sujin Philip's avatar
Sujin Philip committed
627
628
629
630
  {
    this->Clear();
    this->Initialize();
    return *this;
Sujin Philip's avatar
Sujin Philip committed
631
  }
632

Sujin Philip's avatar
Sujin Philip committed
633
634
635
636
637
  vtkCell* GetCell(vtkDataSet *dataset, vtkIdType cellId)
  {
    int celltype = dataset->GetCellType(cellId);
    vtkGenericCell* &gc = this->Cells[celltype];
    if (!gc)
638
    {
Sujin Philip's avatar
Sujin Philip committed
639
640
641
642
643
      gc = vtkGenericCell::New();
    }
    dataset->GetCell(cellId, gc);
    return gc->GetRepresentativeCell();
  }
644

Sujin Philip's avatar
Sujin Philip committed
645
private:
Sujin Philip's avatar
Sujin Philip committed
646
647
  void Initialize()
  {
648
    vtkGenericCell *null = nullptr;
Sujin Philip's avatar
Sujin Philip committed
649
650
651
652
653
654
655
656
657
658
659
660
661
662
    std::fill(this->Cells, this->Cells + VTK_NUMBER_OF_CELL_TYPES, null);
  }

  void Clear()
  {
    for (int i = 0; i < VTK_NUMBER_OF_CELL_TYPES; ++i)
    {
      if (this->Cells[i])
      {
        this->Cells[i]->Delete();
      }
    }
  }

Sujin Philip's avatar
Sujin Philip committed
663
664
  vtkGenericCell *Cells[VTK_NUMBER_OF_CELL_TYPES];
};
665

Sujin Philip's avatar
Sujin Philip committed
666
} // anonymous namespace
667

Sujin Philip's avatar
Sujin Philip committed
668
669
670
671
672
673
674
675
676
677
678
679
class vtkProbeFilter::ProbeImageDataWorklet
{
public:
  ProbeImageDataWorklet(vtkProbeFilter *probeFilter,
                        vtkDataSet *source, int srcBlockId,
                        const double start[3], const double spacing[3],
                        const int dim[3], vtkPointData *outPD, char *maskArray,
                        int maxCellSize)
    : ProbeFilter(probeFilter), Source(source), SrcBlockId(srcBlockId),
      Start(start), Spacing(spacing), Dim(dim), OutPointData(outPD),
      MaskArray(maskArray), MaxCellSize(maxCellSize)
  {
680
  }
681

Sujin Philip's avatar
Sujin Philip committed
682
  void operator()(vtkIdType cellBegin, vtkIdType cellEnd)
683
  {
Sujin Philip's avatar
Sujin Philip committed
684
685
686
687
688
689
690
    double fastweights[256];
    double *weights;
    if (this->MaxCellSize <= 256)
    {
      weights = fastweights;
    }
    else
691
    {
Sujin Philip's avatar
Sujin Philip committed
692
693
694
      std::vector<double> &dynamicweights = this->WeightsBuffer.Local();
      dynamicweights.resize(this->MaxCellSize);
      weights = &dynamicweights[0];
695
    }
Sujin Philip's avatar
Sujin Philip committed
696
697
698

    CellStorage &cs = this->Cells.Local();
    for (vtkIdType cellId = cellBegin; cellId < cellEnd; ++cellId)
699
    {
Sujin Philip's avatar
Sujin Philip committed
700
701
702
703
      vtkCell *cell = cs.GetCell(this->Source, cellId);
      this->ProbeFilter->ProbeImagePointsInCell(cell, cellId, this->Source,
        this->SrcBlockId, this->Start, this->Spacing, this->Dim,
        this->OutPointData, this->MaskArray, weights);
704
    }
705
  }
Sujin Philip's avatar
Sujin Philip committed
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750

private:
  vtkProbeFilter *ProbeFilter;
  vtkDataSet *Source;
  int SrcBlockId;
  const double *Start;
  const double *Spacing;
  const int *Dim;
  vtkPointData *OutPointData;
  char *MaskArray;
  int MaxCellSize;

  vtkSMPThreadLocal<std::vector<double> > WeightsBuffer;
  vtkSMPThreadLocal<CellStorage> Cells;
};

//----------------------------------------------------------------------------
void vtkProbeFilter::ProbePointsImageData(vtkImageData *input,
  int srcIdx, vtkDataSet *source, vtkImageData *output)
{
  vtkPointData *outPD = output->GetPointData();
  char* maskArray = this->MaskPoints->GetPointer(0);

  //----------------------------------------
  double spacing[3];
  input->GetSpacing(spacing);
  int extent[6];
  input->GetExtent(extent);
  int dim[3];
  input->GetDimensions(dim);
  double start[3];
  input->GetOrigin(start);
  start[0] += static_cast<double>(extent[0]) * spacing[0];
  start[1] += static_cast<double>(extent[2]) * spacing[1];
  start[2] += static_cast<double>(extent[4]) * spacing[2];

  vtkIdType numSrcCells = source->GetNumberOfCells();

   // dummy call required before multithreaded calls
  static_cast<void>(source->GetCellType(0));
  ProbeImageDataWorklet worklet(this, source, srcIdx, start, spacing, dim,
                                outPD, maskArray, source->GetMaxCellSize());
  vtkSMPTools::For(0, numSrcCells, worklet);

  this->MaskPoints->Modified();
751
752
}

753
//----------------------------------------------------------------------------
754
755
756
757
int vtkProbeFilter::RequestInformation(
  vtkInformation *vtkNotUsed(request),
  vtkInformationVector **inputVector,
  vtkInformationVector *outputVector)
758
{
759
760
761
762
  // get the info objects
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation *sourceInfo = inputVector[1]->GetInformationObject(0);
  vtkInformation *outInfo = outputVector->GetInformationObject(0);
Charles Law's avatar
Charles Law committed
763

764
  outInfo->CopyEntry(sourceInfo,
765
                     vtkStreamingDemandDrivenPipeline::TIME_STEPS());
766
  outInfo->CopyEntry(sourceInfo,
767
                     vtkStreamingDemandDrivenPipeline::TIME_RANGE());
768

769
770
771
  outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),
               inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT()),
               6);
772

773
774
775
776
  // A variation of the bug fix from John Biddiscombe.
  // Make sure that the scalar type and number of components
  // are propagated from the source not the input.
  if (vtkImageData::HasScalarType(sourceInfo))
777
  {
778
779
    vtkImageData::SetScalarType(vtkImageData::GetScalarType(sourceInfo),
                                outInfo);
780
  }
781
  if (vtkImageData::HasNumberOfScalarComponents(sourceInfo))
782
  {
783
784
785
    vtkImageData::SetNumberOfScalarComponents(
      vtkImageData::GetNumberOfScalarComponents(sourceInfo),
      outInfo);
786
  }
787

788
789
  return 1;
}
790
791

//----------------------------------------------------------------------------
792
793
794
795
int vtkProbeFilter::RequestUpdateExtent(
  vtkInformation *vtkNotUsed(request),
  vtkInformationVector **inputVector,
  vtkInformationVector *outputVector)
796
{
797
798
799
800
801
  // get the info objects
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation *sourceInfo = inputVector[1]->GetInformationObject(0);
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

802
  int usePiece = 0;
803

Kyle Lutz's avatar
Kyle Lutz committed
804
  // What ever happened to CopyUpdateExtent in vtkDataObject?
805
806
  // Copying both piece and extent could be bad.  Setting the piece
  // of a structured data set will affect the extent.
807
808
809
810
  vtkDataObject* output = outInfo->Get(vtkDataObject::DATA_OBJECT());
  if (output &&
      (!strcmp(output->GetClassName(), "vtkUnstructuredGrid") ||
       !strcmp(output->GetClassName(), "vtkPolyData")))
811
  {
812
    usePiece = 1;
813
  }
814

815
  inInfo->Set(vtkStreamingDemandDrivenPipeline::EXACT_EXTENT(), 1);
816

817
818
  sourceInfo->Remove(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT());
  if (sourceInfo->Has(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT()))
819
  {
820
821
    sourceInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(),
      sourceInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT()), 6);
822
  }
823

824
  if ( ! this->SpatialMatch)
825
  {
826
827
828
829
830
831
    sourceInfo->Set(
      vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER(), 0);
    sourceInfo->Set(
      vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES(), 1);
    sourceInfo->Set(
      vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(), 0);
832
  }
833
  else if (this->SpatialMatch == 1)
834
  {
835
    if (usePiece)
836
    {
837
838
839
      // Request an extra ghost level because the probe
      // gets external values with computation prescision problems.
      // I think the probe should be changed to have an epsilon ...
840
841
842
843
844
845
846
847
848
      sourceInfo->Set(
        vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER(),
        outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER()));
      sourceInfo->Set(
        vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES(),
        outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES()));
      sourceInfo->Set(
        vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(),
        outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS())+1);
849
    }
850
    else
851
    {
852
853
854
      sourceInfo->Set(
        vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(),
        outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT()), 6);
855
    }
856
  }
857

858
  if (usePiece)
859
  {
860
861
862
863
864
865
866
867
868
    inInfo->Set(
      vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER(),
      outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER()));
    inInfo->Set(
      vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES(),
      outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES()));
    inInfo->Set(
      vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(),
      outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS()));
Kitware Robot's avatar