1 // -*- c++ -*-
2 /*=========================================================================
3
4 Program: Visualization Toolkit
5 Module: vtkTemporalStatistics.cxx
6
7 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
8 All rights reserved.
9 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
10
11 This software is distributed WITHOUT ANY WARRANTY; without even
12 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
13 PURPOSE. See the above copyright notice for more information.
14
15 =========================================================================*/
16
17 /*
18 * Copyright 2008 Sandia Corporation.
19 * Under the terms of Contract DE-AC04-94AL85000, there is a non-exclusive
20 * license for use of this work by or on behalf of the
21 * U.S. Government. Redistribution and use in source and binary forms, with
22 * or without modification, are permitted provided that this Notice and any
23 * statement of authorship are reproduced on all copies.
24 */
25
26 #include "vtkTemporalStatistics.h"
27
28 #include "vtkCellData.h"
29 #include "vtkCompositeDataIterator.h"
30 #include "vtkCompositeDataSet.h"
31 #include "vtkDataArray.h"
32 #include "vtkDataSet.h"
33 #include "vtkGraph.h"
34 #include "vtkInformation.h"
35 #include "vtkInformationVector.h"
36 #include "vtkMultiBlockDataSet.h"
37 #include "vtkObjectFactory.h"
38 #include "vtkPointData.h"
39 #include "vtkStdString.h"
40 #include "vtkStreamingDemandDrivenPipeline.h"
41
42 #include "vtkSmartPointer.h"
43
44 #include <algorithm>
45
46 //=============================================================================
47 vtkStandardNewMacro(vtkTemporalStatistics);
48
49 //=============================================================================
50 const char * const AVERAGE_SUFFIX = "average";
51 const char * const MINIMUM_SUFFIX = "minimum";
52 const char * const MAXIMUM_SUFFIX = "maximum";
53 const char * const STANDARD_DEVIATION_SUFFIX = "stddev";
54
vtkTemporalStatisticsMangleName(const char * originalName,const char * suffix)55 inline vtkStdString vtkTemporalStatisticsMangleName(const char *originalName,
56 const char *suffix)
57 {
58 if (!originalName) return suffix;
59 return vtkStdString(originalName) + "_" + vtkStdString(suffix);
60 }
61
62 //-----------------------------------------------------------------------------
63 // The interim stddev array keeps a sum of squares.
64 template<class T>
vtkTemporalStatisticsInitializeStdDev(T * outArray,vtkIdType arraySize)65 inline void vtkTemporalStatisticsInitializeStdDev(T *outArray,
66 vtkIdType arraySize)
67 {
68 for (vtkIdType i = 0; i < arraySize; i++)
69 {
70 outArray[i] = 0;
71 }
72 }
73
74 //-----------------------------------------------------------------------------
75 template<class T>
vtkTemporalStatisticsAccumulateAverage(const T * inArray,T * outArray,vtkIdType arraySize)76 inline void vtkTemporalStatisticsAccumulateAverage(const T *inArray,
77 T *outArray,
78 vtkIdType arraySize)
79 {
80 for (vtkIdType i = 0; i < arraySize; i++)
81 {
82 outArray[i] += inArray[i];
83 }
84 }
85
86 template<class T>
vtkTemporalStatisticsAccumulateMinimum(const T * inArray,T * outArray,vtkIdType arraySize)87 inline void vtkTemporalStatisticsAccumulateMinimum(const T *inArray,
88 T *outArray,
89 vtkIdType arraySize)
90 {
91 for (vtkIdType i = 0; i < arraySize; i++)
92 {
93 if (outArray[i] > inArray[i]) outArray[i] = inArray[i];
94 }
95 }
96
97 template<class T>
vtkTemporalStatisticsAccumulateMaximum(const T * inArray,T * outArray,vtkIdType arraySize)98 inline void vtkTemporalStatisticsAccumulateMaximum(const T *inArray,
99 T *outArray,
100 vtkIdType arraySize)
101 {
102 for (vtkIdType i = 0; i < arraySize; i++)
103 {
104 if (outArray[i] < inArray[i]) outArray[i] = inArray[i];
105 }
106 }
107
108 // standard deviation one-pass algorithm from
109 // http://www.cs.berkeley.edu/~mhoemmen/cs194/Tutorials/variance.pdf
110 // this is numerically stable!
111 template<class T>
vtkTemporalStatisticsAccumulateStdDev(const T * inArray,T * outArray,const T * previousAverage,vtkIdType arraySize,int pass)112 inline void vtkTemporalStatisticsAccumulateStdDev(
113 const T *inArray, T *outArray, const T *previousAverage,
114 vtkIdType arraySize, int pass)
115 {
116 for (vtkIdType i = 0; i < arraySize; i++)
117 {
118 double temp = inArray[i]-previousAverage[i]/static_cast<double>(pass);
119 outArray[i] = outArray[i] + static_cast<T>(
120 pass*temp*temp/static_cast<double>(pass+1) );
121 }
122 }
123
124 //-----------------------------------------------------------------------------
125 template<class T>
vtkTemporalStatisticsFinishAverage(T * outArray,vtkIdType arraySize,int sumSize)126 inline void vtkTemporalStatisticsFinishAverage(T *outArray, vtkIdType arraySize,
127 int sumSize)
128 {
129 for (vtkIdType i = 0; i < arraySize; i++)
130 {
131 outArray[i] /= sumSize;
132 }
133 }
134
135 template<class T>
vtkTemporalStatisticsFinishStdDev(T * outArray,vtkIdType arraySize,int sumSize)136 inline void vtkTemporalStatisticsFinishStdDev(T *outArray,
137 vtkIdType arraySize, int sumSize)
138 {
139 for (vtkIdType i = 0; i < arraySize; i++)
140 {
141 outArray[i] =
142 static_cast<T>(sqrt(static_cast<double>(outArray[i])/sumSize));
143 }
144 }
145
146 //=============================================================================
vtkTemporalStatistics()147 vtkTemporalStatistics::vtkTemporalStatistics()
148 {
149 this->ComputeAverage = 1;
150 this->ComputeMinimum = 1;
151 this->ComputeMaximum = 1;
152 this->ComputeStandardDeviation = 1;
153
154 this->CurrentTimeIndex = 0;
155 this->GeneratedChangingTopologyWarning = false;
156 }
157
158 vtkTemporalStatistics::~vtkTemporalStatistics() = default;
159
PrintSelf(ostream & os,vtkIndent indent)160 void vtkTemporalStatistics::PrintSelf(ostream &os, vtkIndent indent)
161 {
162 this->Superclass::PrintSelf(os, indent);
163
164 os << indent << "ComputeAverage: " << this->ComputeAverage << endl;
165 os << indent << "ComputeMinimum: " << this->ComputeMinimum << endl;
166 os << indent << "ComputeMaximum: " << this->ComputeMaximum << endl;
167 os << indent << "ComputeStandardDeviation: " <<
168 this->ComputeStandardDeviation << endl;
169 }
170
171 //-----------------------------------------------------------------------------
FillInputPortInformation(int vtkNotUsed (port),vtkInformation * info)172 int vtkTemporalStatistics::FillInputPortInformation(int vtkNotUsed(port),
173 vtkInformation *info)
174 {
175 info->Remove(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE());
176 info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
177 info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkGraph");
178 info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkCompositeDataSet");
179 return 1;
180 }
181
182 //-----------------------------------------------------------------------------
RequestInformation(vtkInformation * vtkNotUsed (request),vtkInformationVector ** vtkNotUsed (inputVector),vtkInformationVector * outputVector)183 int vtkTemporalStatistics::RequestInformation(
184 vtkInformation *vtkNotUsed(request),
185 vtkInformationVector **vtkNotUsed(inputVector),
186 vtkInformationVector *outputVector)
187 {
188 vtkInformation *outInfo = outputVector->GetInformationObject(0);
189
190 // The output data of this filter has no time associated with it. It is the
191 // result of computations that happen over all time.
192 outInfo->Remove(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
193 outInfo->Remove(vtkStreamingDemandDrivenPipeline::TIME_RANGE());
194
195 return 1;
196 }
197
198 //-----------------------------------------------------------------------------
RequestDataObject(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)199 int vtkTemporalStatistics::RequestDataObject(
200 vtkInformation *vtkNotUsed(request),
201 vtkInformationVector **inputVector,
202 vtkInformationVector *outputVector)
203 {
204 vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
205 vtkInformation *outInfo = outputVector->GetInformationObject(0);
206
207 vtkDataObject *input = vtkDataObject::GetData(inInfo);
208 vtkDataObject *output = vtkDataObject::GetData(outInfo);
209
210 if (!input)
211 {
212 return 0;
213 }
214
215 vtkSmartPointer<vtkDataObject> newOutput;
216
217 if (!output || !output->IsA(input->GetClassName()))
218 {
219 newOutput.TakeReference(input->NewInstance());
220 }
221
222
223 if (newOutput)
224 {
225 outInfo->Set(vtkDataObject::DATA_OBJECT(), newOutput);
226 }
227
228 return 1;
229 }
230
231 //-----------------------------------------------------------------------------
RequestUpdateExtent(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * vtkNotUsed (outputVector))232 int vtkTemporalStatistics::RequestUpdateExtent(
233 vtkInformation *vtkNotUsed(request),
234 vtkInformationVector **inputVector,
235 vtkInformationVector *vtkNotUsed(outputVector))
236 {
237 vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
238
239 // The RequestData method will tell the pipeline executive to iterate the
240 // upstream pipeline to get each time step in order. The executive in turn
241 // will call this method to get the extent request for each iteration (in this
242 // case the time step).
243 double *inTimes = inInfo->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
244 if (inTimes)
245 {
246 inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP(), inTimes[this->CurrentTimeIndex]);
247 }
248
249 return 1;
250 }
251
252 //-----------------------------------------------------------------------------
RequestData(vtkInformation * request,vtkInformationVector ** inputVector,vtkInformationVector * outputVector)253 int vtkTemporalStatistics::RequestData(vtkInformation *request,
254 vtkInformationVector **inputVector,
255 vtkInformationVector *outputVector)
256 {
257 vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
258 vtkInformation *outInfo = outputVector->GetInformationObject(0);
259
260 vtkDataObject *input = vtkDataObject::GetData(inInfo);
261 vtkDataObject *output = vtkDataObject::GetData(outInfo);
262
263 if (this->CurrentTimeIndex == 0)
264 {
265 // First execution, initialize arrays.
266 this->InitializeStatistics(input, output);
267 }
268 else
269 {
270 // Subsequent execution, accumulate new data.
271 this->AccumulateStatistics(input, output);
272 }
273
274 this->CurrentTimeIndex++;
275
276 if ( this->CurrentTimeIndex
277 < inInfo->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS()))
278 {
279 // There is still more to do.
280 request->Set(vtkStreamingDemandDrivenPipeline::CONTINUE_EXECUTING(), 1);
281 }
282 else
283 {
284 // We are done. Finish up.
285 this->PostExecute(input, output);
286 request->Remove(vtkStreamingDemandDrivenPipeline::CONTINUE_EXECUTING());
287 this->CurrentTimeIndex = 0;
288 }
289
290 return 1;
291 }
292
293 //-----------------------------------------------------------------------------
InitializeStatistics(vtkDataObject * input,vtkDataObject * output)294 void vtkTemporalStatistics::InitializeStatistics(vtkDataObject *input,
295 vtkDataObject *output)
296 {
297 if (input->IsA("vtkDataSet"))
298 {
299 this->InitializeStatistics(vtkDataSet::SafeDownCast(input),
300 vtkDataSet::SafeDownCast(output));
301 return;
302 }
303
304 if (input->IsA("vtkGraph"))
305 {
306 this->InitializeStatistics(vtkGraph::SafeDownCast(input),
307 vtkGraph::SafeDownCast(output));
308 return;
309 }
310
311 if (input->IsA("vtkCompositeDataSet"))
312 {
313 this->InitializeStatistics(vtkCompositeDataSet::SafeDownCast(input),
314 vtkCompositeDataSet::SafeDownCast(output));
315 return;
316 }
317
318 vtkWarningMacro(<< "Unsupported input type: " << input->GetClassName());
319 }
320
321 //-----------------------------------------------------------------------------
InitializeStatistics(vtkDataSet * input,vtkDataSet * output)322 void vtkTemporalStatistics::InitializeStatistics(vtkDataSet *input,
323 vtkDataSet *output)
324 {
325 output->CopyStructure(input);
326 this->InitializeArrays(input->GetFieldData(), output->GetFieldData());
327 this->InitializeArrays(input->GetPointData(), output->GetPointData());
328 this->InitializeArrays(input->GetCellData(), output->GetCellData());
329 }
330
331 //-----------------------------------------------------------------------------
InitializeStatistics(vtkGraph * input,vtkGraph * output)332 void vtkTemporalStatistics::InitializeStatistics(vtkGraph *input,
333 vtkGraph *output)
334 {
335 output->CopyStructure(input);
336 this->InitializeArrays(input->GetFieldData(), output->GetFieldData());
337 this->InitializeArrays(input->GetVertexData(), output->GetVertexData());
338 this->InitializeArrays(input->GetEdgeData(), output->GetEdgeData());
339 }
340
341 //-----------------------------------------------------------------------------
InitializeStatistics(vtkCompositeDataSet * input,vtkCompositeDataSet * output)342 void vtkTemporalStatistics::InitializeStatistics(vtkCompositeDataSet *input,
343 vtkCompositeDataSet *output)
344 {
345 output->CopyStructure(input);
346
347 vtkSmartPointer<vtkCompositeDataIterator> inputItr;
348 inputItr.TakeReference(input->NewIterator());
349
350 for (inputItr->InitTraversal(); !inputItr->IsDoneWithTraversal();
351 inputItr->GoToNextItem())
352 {
353 vtkDataObject *inputObj = inputItr->GetCurrentDataObject();
354
355 vtkSmartPointer<vtkDataObject> outputObj;
356 outputObj.TakeReference(inputObj->NewInstance());
357
358 this->InitializeStatistics(inputObj, outputObj);
359 output->SetDataSet(inputItr, outputObj);
360 }
361 }
362
363 //-----------------------------------------------------------------------------
InitializeArrays(vtkFieldData * inFd,vtkFieldData * outFd)364 void vtkTemporalStatistics::InitializeArrays(vtkFieldData *inFd,
365 vtkFieldData *outFd)
366 {
367 // Because we need to do mathematical operations, we require all arrays we
368 // process to be numeric data (i.e. a vtkDataArray). We also handle global
369 // ids and petigree ids special (we just pass them). Ideally would just let
370 // vtkFieldData or vtkDataSetAttributes handle this for us, but no such method
371 // that fits our needs here. Thus, we pass data a bit differently then other
372 // filters. If I miss something important, it should be added here.
373
374 outFd->Initialize();
375
376 vtkDataSetAttributes *inDsa = vtkDataSetAttributes::SafeDownCast(inFd);
377 vtkDataSetAttributes *outDsa = vtkDataSetAttributes::SafeDownCast(outFd);
378 if (inDsa)
379 {
380 vtkDataArray *globalIds = inDsa->GetGlobalIds();
381 vtkAbstractArray *pedigreeIds = inDsa->GetPedigreeIds();
382 if (globalIds) outDsa->SetGlobalIds(globalIds);
383 if (pedigreeIds) outDsa->SetPedigreeIds(pedigreeIds);
384 }
385
386 int numArrays = inFd->GetNumberOfArrays();
387 for (int i = 0; i < numArrays; i++)
388 {
389 vtkDataArray *array = inFd->GetArray(i);
390 if (!array) continue; // Array not numeric.
391 if (outFd->HasArray(array->GetName())) continue; // Must be Ids.
392
393 this->InitializeArray(array, outFd);
394 }
395 }
396
397 //-----------------------------------------------------------------------------
InitializeArray(vtkDataArray * array,vtkFieldData * outFd)398 void vtkTemporalStatistics::InitializeArray(vtkDataArray *array,
399 vtkFieldData *outFd)
400 {
401 if (this->ComputeAverage || this->ComputeStandardDeviation)
402 {
403 vtkSmartPointer<vtkDataArray> newArray;
404 newArray.TakeReference(vtkArrayDownCast<vtkDataArray>(
405 vtkAbstractArray::CreateArray(array->GetDataType())));
406 newArray->DeepCopy(array);
407 newArray->SetName(vtkTemporalStatisticsMangleName(array->GetName(),
408 AVERAGE_SUFFIX));
409 if (outFd->HasArray(newArray->GetName()))
410 {
411 vtkWarningMacro(<< "Input has two arrays named " << array->GetName()
412 << ". Output statistics will probably be wrong.");
413 return;
414 }
415 outFd->AddArray(newArray);
416 }
417
418 if (this->ComputeMinimum)
419 {
420 vtkSmartPointer<vtkDataArray> newArray;
421 newArray.TakeReference(vtkArrayDownCast<vtkDataArray>(
422 vtkAbstractArray::CreateArray(array->GetDataType())));
423 newArray->DeepCopy(array);
424 newArray->SetName(vtkTemporalStatisticsMangleName(array->GetName(),
425 MINIMUM_SUFFIX));
426 outFd->AddArray(newArray);
427 }
428
429 if (this->ComputeMaximum)
430 {
431 vtkSmartPointer<vtkDataArray> newArray;
432 newArray.TakeReference(vtkArrayDownCast<vtkDataArray>(
433 vtkAbstractArray::CreateArray(array->GetDataType())));
434 newArray->DeepCopy(array);
435 newArray->SetName(vtkTemporalStatisticsMangleName(array->GetName(),
436 MAXIMUM_SUFFIX));
437 outFd->AddArray(newArray);
438 }
439
440 if (this->ComputeStandardDeviation)
441 {
442 vtkSmartPointer<vtkDataArray> newArray;
443 newArray.TakeReference(vtkArrayDownCast<vtkDataArray>(
444 vtkAbstractArray::CreateArray(array->GetDataType())));
445 newArray->SetName(vtkTemporalStatisticsMangleName(array->GetName(),
446 STANDARD_DEVIATION_SUFFIX));
447
448 newArray->SetNumberOfComponents(array->GetNumberOfComponents());
449 newArray->CopyComponentNames( array );
450
451 newArray->SetNumberOfTuples(array->GetNumberOfTuples());
452 switch (array->GetDataType())
453 {
454 vtkTemplateMacro(vtkTemporalStatisticsInitializeStdDev(
455 static_cast<VTK_TT*>(newArray->GetVoidPointer(0)),
456 array->GetNumberOfComponents()
457 *array->GetNumberOfTuples()));
458 }
459 outFd->AddArray(newArray);
460 }
461 }
462
463 //-----------------------------------------------------------------------------
AccumulateStatistics(vtkDataObject * input,vtkDataObject * output)464 void vtkTemporalStatistics::AccumulateStatistics(vtkDataObject *input,
465 vtkDataObject *output)
466 {
467 if (input->IsA("vtkDataSet"))
468 {
469 this->AccumulateStatistics(vtkDataSet::SafeDownCast(input),
470 vtkDataSet::SafeDownCast(output));
471 return;
472 }
473
474 if (input->IsA("vtkGraph"))
475 {
476 this->AccumulateStatistics(vtkGraph::SafeDownCast(input),
477 vtkGraph::SafeDownCast(output));
478 return;
479 }
480
481 if (input->IsA("vtkCompositeDataSet"))
482 {
483 this->AccumulateStatistics(vtkCompositeDataSet::SafeDownCast(input),
484 vtkCompositeDataSet::SafeDownCast(output));
485 }
486 }
487
488 //-----------------------------------------------------------------------------
AccumulateStatistics(vtkDataSet * input,vtkDataSet * output)489 void vtkTemporalStatistics::AccumulateStatistics(vtkDataSet *input,
490 vtkDataSet *output)
491 {
492 this->AccumulateArrays(input->GetFieldData(), output->GetFieldData());
493 this->AccumulateArrays(input->GetPointData(), output->GetPointData());
494 this->AccumulateArrays(input->GetCellData(), output->GetCellData());
495 }
496
497 //-----------------------------------------------------------------------------
AccumulateStatistics(vtkGraph * input,vtkGraph * output)498 void vtkTemporalStatistics::AccumulateStatistics(vtkGraph *input,
499 vtkGraph *output)
500 {
501 this->AccumulateArrays(input->GetFieldData(), output->GetFieldData());
502 this->AccumulateArrays(input->GetVertexData(), output->GetVertexData());
503 this->AccumulateArrays(input->GetEdgeData(), output->GetEdgeData());
504 }
505
506 //-----------------------------------------------------------------------------
AccumulateStatistics(vtkCompositeDataSet * input,vtkCompositeDataSet * output)507 void vtkTemporalStatistics::AccumulateStatistics(vtkCompositeDataSet *input,
508 vtkCompositeDataSet *output)
509 {
510 vtkSmartPointer<vtkCompositeDataIterator> inputItr;
511 inputItr.TakeReference(input->NewIterator());
512
513 for (inputItr->InitTraversal(); !inputItr->IsDoneWithTraversal();
514 inputItr->GoToNextItem())
515 {
516 vtkDataObject *inputObj = inputItr->GetCurrentDataObject();
517 vtkDataObject *outputObj = output->GetDataSet(inputItr);
518
519 this->AccumulateStatistics(inputObj, outputObj);
520 }
521 }
522
523 //-----------------------------------------------------------------------------
AccumulateArrays(vtkFieldData * inFd,vtkFieldData * outFd)524 void vtkTemporalStatistics::AccumulateArrays(vtkFieldData *inFd,
525 vtkFieldData *outFd)
526 {
527 int numArrays = inFd->GetNumberOfArrays();
528 for (int i = 0; i < numArrays; i++)
529 {
530 vtkDataArray *inArray = inFd->GetArray(i);
531 vtkDataArray *outArray;
532 if (!inArray) continue;
533
534 outArray = this->GetArray(outFd, inArray, AVERAGE_SUFFIX);
535 if (outArray)
536 {
537
538
539 vtkDataArray* stdevOutArray =
540 this->GetArray(outFd, inArray, STANDARD_DEVIATION_SUFFIX);
541 if (stdevOutArray)
542 {
543 switch (inArray->GetDataType())
544 {
545 // standard deviation must be called before average since the one-pass
546 // algorithm uses the average up to the previous time step
547 vtkTemplateMacro(vtkTemporalStatisticsAccumulateStdDev(
548 static_cast<const VTK_TT*>(inArray->GetVoidPointer(0)),
549 static_cast<VTK_TT*>(stdevOutArray->GetVoidPointer(0)),
550 static_cast<const VTK_TT*>(outArray->GetVoidPointer(0)),
551 inArray->GetNumberOfComponents()*inArray->GetNumberOfTuples(),
552 this->CurrentTimeIndex));
553 }
554 // Alert change in data.
555 stdevOutArray->DataChanged();
556 }
557
558
559
560
561 switch (inArray->GetDataType())
562 {
563 vtkTemplateMacro(vtkTemporalStatisticsAccumulateAverage(
564 static_cast<const VTK_TT*>(inArray->GetVoidPointer(0)),
565 static_cast<VTK_TT*>(outArray->GetVoidPointer(0)),
566 inArray->GetNumberOfComponents()
567 *inArray->GetNumberOfTuples()));
568 }
569 // Alert change in data.
570 outArray->DataChanged();
571 }
572
573 outArray = this->GetArray(outFd, inArray, MINIMUM_SUFFIX);
574 if (outArray)
575 {
576 switch (inArray->GetDataType())
577 {
578 vtkTemplateMacro(vtkTemporalStatisticsAccumulateMinimum(
579 static_cast<const VTK_TT*>(inArray->GetVoidPointer(0)),
580 static_cast<VTK_TT*>(outArray->GetVoidPointer(0)),
581 inArray->GetNumberOfComponents()
582 *inArray->GetNumberOfTuples()));
583 }
584 // Alert change in data.
585 outArray->DataChanged();
586 }
587
588 outArray = this->GetArray(outFd, inArray, MAXIMUM_SUFFIX);
589 if (outArray)
590 {
591 switch (inArray->GetDataType())
592 {
593 vtkTemplateMacro(vtkTemporalStatisticsAccumulateMaximum(
594 static_cast<const VTK_TT*>(inArray->GetVoidPointer(0)),
595 static_cast<VTK_TT*>(outArray->GetVoidPointer(0)),
596 inArray->GetNumberOfComponents()
597 *inArray->GetNumberOfTuples()));
598 }
599 // Alert change in data.
600 outArray->DataChanged();
601 }
602 }
603 }
604
605 //-----------------------------------------------------------------------------
PostExecute(vtkDataObject * input,vtkDataObject * output)606 void vtkTemporalStatistics::PostExecute(vtkDataObject *input,
607 vtkDataObject *output)
608 {
609 if (input->IsA("vtkDataSet"))
610 {
611 this->PostExecute(vtkDataSet::SafeDownCast(input),
612 vtkDataSet::SafeDownCast(output));
613 return;
614 }
615
616 if (input->IsA("vtkGraph"))
617 {
618 this->PostExecute(vtkGraph::SafeDownCast(input),
619 vtkGraph::SafeDownCast(output));
620 return;
621 }
622
623 if (input->IsA("vtkCompositeDataSet"))
624 {
625 this->PostExecute(vtkCompositeDataSet::SafeDownCast(input),
626 vtkCompositeDataSet::SafeDownCast(output));
627 }
628 }
629
630 //-----------------------------------------------------------------------------
PostExecute(vtkDataSet * input,vtkDataSet * output)631 void vtkTemporalStatistics::PostExecute(vtkDataSet *input, vtkDataSet *output)
632 {
633 this->FinishArrays(input->GetFieldData(), output->GetFieldData());
634 this->FinishArrays(input->GetPointData(), output->GetPointData());
635 this->FinishArrays(input->GetCellData(), output->GetCellData());
636 }
637
638 //-----------------------------------------------------------------------------
PostExecute(vtkGraph * input,vtkGraph * output)639 void vtkTemporalStatistics::PostExecute(vtkGraph *input, vtkGraph *output)
640 {
641 this->FinishArrays(input->GetFieldData(), output->GetFieldData());
642 this->FinishArrays(input->GetVertexData(), output->GetVertexData());
643 this->FinishArrays(input->GetEdgeData(), output->GetEdgeData());
644 }
645
646 //-----------------------------------------------------------------------------
PostExecute(vtkCompositeDataSet * input,vtkCompositeDataSet * output)647 void vtkTemporalStatistics::PostExecute(vtkCompositeDataSet *input,
648 vtkCompositeDataSet *output)
649 {
650 vtkSmartPointer<vtkCompositeDataIterator> inputItr;
651 inputItr.TakeReference(input->NewIterator());
652
653 for (inputItr->InitTraversal(); !inputItr->IsDoneWithTraversal();
654 inputItr->GoToNextItem())
655 {
656 vtkDataObject *inputObj = inputItr->GetCurrentDataObject();
657 vtkDataObject *outputObj = output->GetDataSet(inputItr);
658
659 this->PostExecute(inputObj, outputObj);
660 }
661 }
662
663 //-----------------------------------------------------------------------------
FinishArrays(vtkFieldData * inFd,vtkFieldData * outFd)664 void vtkTemporalStatistics::FinishArrays(vtkFieldData *inFd,
665 vtkFieldData *outFd)
666 {
667 int numArrays = inFd->GetNumberOfArrays();
668 for (int i = 0; i < numArrays; i++)
669 {
670 vtkDataArray *inArray = inFd->GetArray(i);
671 vtkDataArray *outArray;
672 if (!inArray) continue;
673
674 outArray = this->GetArray(outFd, inArray, AVERAGE_SUFFIX);
675 if (outArray)
676 {
677 switch (inArray->GetDataType())
678 {
679 vtkTemplateMacro(vtkTemporalStatisticsFinishAverage(
680 static_cast<VTK_TT*>(outArray->GetVoidPointer(0)),
681 inArray->GetNumberOfComponents()
682 *inArray->GetNumberOfTuples(),
683 this->CurrentTimeIndex));
684 }
685 }
686 vtkDataArray *avgArray = outArray;
687
688 // No post processing on minimum.
689 // No post processing on maximum.
690
691 outArray = this->GetArray(outFd, inArray, STANDARD_DEVIATION_SUFFIX);
692 if (outArray)
693 {
694 if (!avgArray)
695 {
696 vtkWarningMacro(<< "Average not computed for " << inArray->GetName()
697 << ", standard deviation skipped.");
698 outFd->RemoveArray(outArray->GetName());
699 }
700 else
701 {
702 switch (inArray->GetDataType())
703 {
704 vtkTemplateMacro(vtkTemporalStatisticsFinishStdDev(
705 static_cast<VTK_TT*>(outArray->GetVoidPointer(0)),
706 inArray->GetNumberOfComponents()
707 *inArray->GetNumberOfTuples(),
708 this->CurrentTimeIndex));
709 }
710 if (!this->ComputeAverage)
711 {
712 outFd->RemoveArray(avgArray->GetName());
713 }
714 }
715 }
716 }
717 }
718
719 //-----------------------------------------------------------------------------
GetArray(vtkFieldData * fieldData,vtkDataArray * inArray,const char * nameSuffix)720 vtkDataArray *vtkTemporalStatistics::GetArray(vtkFieldData *fieldData,
721 vtkDataArray *inArray,
722 const char *nameSuffix)
723 {
724 vtkStdString outArrayName
725 = vtkTemporalStatisticsMangleName(inArray->GetName(), nameSuffix);
726 vtkDataArray *outArray = fieldData->GetArray(outArrayName.c_str());
727 if (!outArray) return nullptr;
728
729 if ( (inArray->GetNumberOfComponents() != outArray->GetNumberOfComponents())
730 || (inArray->GetNumberOfTuples() != outArray->GetNumberOfTuples()) )
731 {
732 if(!this->GeneratedChangingTopologyWarning)
733 {
734 std::string fieldType = vtkCellData::SafeDownCast(fieldData) == nullptr ?
735 "points" : "cells";
736 vtkWarningMacro("The number of " << fieldType << " has changed between time "
737 << "steps. No arrays of this type will be output since this "
738 << "filter can not handle grids that change over time.");
739 this->GeneratedChangingTopologyWarning = true;
740 }
741 fieldData->RemoveArray(outArray->GetName());
742 return nullptr;
743 }
744
745 return outArray;
746 }
747