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 }
156
~vtkTemporalStatistics()157 vtkTemporalStatistics::~vtkTemporalStatistics()
158 {
159 }
160
PrintSelf(ostream & os,vtkIndent indent)161 void vtkTemporalStatistics::PrintSelf(ostream &os, vtkIndent indent)
162 {
163 this->Superclass::PrintSelf(os, indent);
164
165 os << indent << "ComputeAverage: " << this->ComputeAverage << endl;
166 os << indent << "ComputeMinimum: " << this->ComputeMinimum << endl;
167 os << indent << "ComputeMaximum: " << this->ComputeMaximum << endl;
168 os << indent << "ComputeStandardDeviation: " <<
169 this->ComputeStandardDeviation << endl;
170 }
171
172 //-----------------------------------------------------------------------------
FillInputPortInformation(int vtkNotUsed (port),vtkInformation * info)173 int vtkTemporalStatistics::FillInputPortInformation(int vtkNotUsed(port),
174 vtkInformation *info)
175 {
176 info->Remove(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE());
177 info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
178 info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkGraph");
179 info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkCompositeDataSet");
180 return 1;
181 }
182
183 //-----------------------------------------------------------------------------
RequestInformation(vtkInformation * vtkNotUsed (request),vtkInformationVector ** vtkNotUsed (inputVector),vtkInformationVector * outputVector)184 int vtkTemporalStatistics::RequestInformation(
185 vtkInformation *vtkNotUsed(request),
186 vtkInformationVector **vtkNotUsed(inputVector),
187 vtkInformationVector *outputVector)
188 {
189 vtkInformation *outInfo = outputVector->GetInformationObject(0);
190
191 // The output data of this filter has no time associated with it. It is the
192 // result of computations that happen over all time.
193 outInfo->Remove(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
194 outInfo->Remove(vtkStreamingDemandDrivenPipeline::TIME_RANGE());
195
196 return 1;
197 }
198
199 //-----------------------------------------------------------------------------
RequestDataObject(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)200 int vtkTemporalStatistics::RequestDataObject(
201 vtkInformation *vtkNotUsed(request),
202 vtkInformationVector **inputVector,
203 vtkInformationVector *outputVector)
204 {
205 vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
206 vtkInformation *outInfo = outputVector->GetInformationObject(0);
207
208 vtkDataObject *input = vtkDataObject::GetData(inInfo);
209 vtkDataObject *output = vtkDataObject::GetData(outInfo);
210
211 if (!input)
212 {
213 return 0;
214 }
215
216 vtkSmartPointer<vtkDataObject> newOutput;
217
218 if (!output || !output->IsA(input->GetClassName()))
219 {
220 newOutput.TakeReference(input->NewInstance());
221 }
222
223
224 if (newOutput)
225 {
226 outInfo->Set(vtkDataObject::DATA_OBJECT(), newOutput);
227 }
228
229 return 1;
230 }
231
232 //-----------------------------------------------------------------------------
RequestUpdateExtent(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * vtkNotUsed (outputVector))233 int vtkTemporalStatistics::RequestUpdateExtent(
234 vtkInformation *vtkNotUsed(request),
235 vtkInformationVector **inputVector,
236 vtkInformationVector *vtkNotUsed(outputVector))
237 {
238 vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
239
240 // The RequestData method will tell the pipeline executive to iterate the
241 // upstream pipeline to get each time step in order. The executive in turn
242 // will call this method to get the extent request for each iteration (in this
243 // case the time step).
244 double *inTimes = inInfo->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
245 if (inTimes)
246 {
247 inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP(), inTimes[this->CurrentTimeIndex]);
248 }
249
250 return 1;
251 }
252
253 //-----------------------------------------------------------------------------
RequestData(vtkInformation * request,vtkInformationVector ** inputVector,vtkInformationVector * outputVector)254 int vtkTemporalStatistics::RequestData(vtkInformation *request,
255 vtkInformationVector **inputVector,
256 vtkInformationVector *outputVector)
257 {
258 vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
259 vtkInformation *outInfo = outputVector->GetInformationObject(0);
260
261 vtkDataObject *input = vtkDataObject::GetData(inInfo);
262 vtkDataObject *output = vtkDataObject::GetData(outInfo);
263
264 if (this->CurrentTimeIndex == 0)
265 {
266 // First execution, initialize arrays.
267 this->InitializeStatistics(input, output);
268 }
269 else
270 {
271 // Subsequent execution, accumulate new data.
272 this->AccumulateStatistics(input, output);
273 }
274
275 this->CurrentTimeIndex++;
276
277 if ( this->CurrentTimeIndex
278 < inInfo->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS()))
279 {
280 // There is still more to do.
281 request->Set(vtkStreamingDemandDrivenPipeline::CONTINUE_EXECUTING(), 1);
282 }
283 else
284 {
285 // We are done. Finish up.
286 this->PostExecute(input, output);
287 request->Remove(vtkStreamingDemandDrivenPipeline::CONTINUE_EXECUTING());
288 this->CurrentTimeIndex = 0;
289 }
290
291 return 1;
292 }
293
294 //-----------------------------------------------------------------------------
InitializeStatistics(vtkDataObject * input,vtkDataObject * output)295 void vtkTemporalStatistics::InitializeStatistics(vtkDataObject *input,
296 vtkDataObject *output)
297 {
298 if (input->IsA("vtkDataSet"))
299 {
300 this->InitializeStatistics(vtkDataSet::SafeDownCast(input),
301 vtkDataSet::SafeDownCast(output));
302 return;
303 }
304
305 if (input->IsA("vtkGraph"))
306 {
307 this->InitializeStatistics(vtkGraph::SafeDownCast(input),
308 vtkGraph::SafeDownCast(output));
309 return;
310 }
311
312 if (input->IsA("vtkCompositeDataSet"))
313 {
314 this->InitializeStatistics(vtkCompositeDataSet::SafeDownCast(input),
315 vtkCompositeDataSet::SafeDownCast(output));
316 return;
317 }
318
319 vtkWarningMacro(<< "Unsupported input type: " << input->GetClassName());
320 }
321
322 //-----------------------------------------------------------------------------
InitializeStatistics(vtkDataSet * input,vtkDataSet * output)323 void vtkTemporalStatistics::InitializeStatistics(vtkDataSet *input,
324 vtkDataSet *output)
325 {
326 output->CopyStructure(input);
327 this->InitializeArrays(input->GetFieldData(), output->GetFieldData());
328 this->InitializeArrays(input->GetPointData(), output->GetPointData());
329 this->InitializeArrays(input->GetCellData(), output->GetCellData());
330 }
331
332 //-----------------------------------------------------------------------------
InitializeStatistics(vtkGraph * input,vtkGraph * output)333 void vtkTemporalStatistics::InitializeStatistics(vtkGraph *input,
334 vtkGraph *output)
335 {
336 output->CopyStructure(input);
337 this->InitializeArrays(input->GetFieldData(), output->GetFieldData());
338 this->InitializeArrays(input->GetVertexData(), output->GetVertexData());
339 this->InitializeArrays(input->GetEdgeData(), output->GetEdgeData());
340 }
341
342 //-----------------------------------------------------------------------------
InitializeStatistics(vtkCompositeDataSet * input,vtkCompositeDataSet * output)343 void vtkTemporalStatistics::InitializeStatistics(vtkCompositeDataSet *input,
344 vtkCompositeDataSet *output)
345 {
346 output->CopyStructure(input);
347
348 vtkSmartPointer<vtkCompositeDataIterator> inputItr;
349 inputItr.TakeReference(input->NewIterator());
350
351 for (inputItr->InitTraversal(); !inputItr->IsDoneWithTraversal();
352 inputItr->GoToNextItem())
353 {
354 vtkDataObject *inputObj = inputItr->GetCurrentDataObject();
355
356 vtkSmartPointer<vtkDataObject> outputObj;
357 outputObj.TakeReference(inputObj->NewInstance());
358
359 output->SetDataSet(inputItr, outputObj);
360
361 this->InitializeStatistics(inputObj, outputObj);
362 }
363 }
364
365 //-----------------------------------------------------------------------------
InitializeArrays(vtkFieldData * inFd,vtkFieldData * outFd)366 void vtkTemporalStatistics::InitializeArrays(vtkFieldData *inFd,
367 vtkFieldData *outFd)
368 {
369 // Because we need to do mathematical operations, we require all arrays we
370 // process to be numeric data (i.e. a vtkDataArray). We also handle global
371 // ids and petigree ids special (we just pass them). Ideally would just let
372 // vtkFieldData or vtkDataSetAttributes handle this for us, but no such method
373 // that fits our needs here. Thus, we pass data a bit differently then other
374 // filters. If I miss something important, it should be added here.
375
376 outFd->Initialize();
377
378 vtkDataSetAttributes *inDsa = vtkDataSetAttributes::SafeDownCast(inFd);
379 vtkDataSetAttributes *outDsa = vtkDataSetAttributes::SafeDownCast(outFd);
380 if (inDsa)
381 {
382 vtkDataArray *globalIds = inDsa->GetGlobalIds();
383 vtkAbstractArray *pedigreeIds = inDsa->GetPedigreeIds();
384 if (globalIds) outDsa->SetGlobalIds(globalIds);
385 if (pedigreeIds) outDsa->SetPedigreeIds(pedigreeIds);
386 }
387
388 int numArrays = inFd->GetNumberOfArrays();
389 for (int i = 0; i < numArrays; i++)
390 {
391 vtkDataArray *array = inFd->GetArray(i);
392 if (!array) continue; // Array not numeric.
393 if (outFd->HasArray(array->GetName())) continue; // Must be Ids.
394
395 this->InitializeArray(array, outFd);
396 }
397 }
398
399 //-----------------------------------------------------------------------------
InitializeArray(vtkDataArray * array,vtkFieldData * outFd)400 void vtkTemporalStatistics::InitializeArray(vtkDataArray *array,
401 vtkFieldData *outFd)
402 {
403 if (this->ComputeAverage || this->ComputeStandardDeviation)
404 {
405 vtkSmartPointer<vtkDataArray> newArray;
406 newArray.TakeReference(vtkDataArray::SafeDownCast(
407 vtkAbstractArray::CreateArray(array->GetDataType())));
408 newArray->DeepCopy(array);
409 newArray->SetName(vtkTemporalStatisticsMangleName(array->GetName(),
410 AVERAGE_SUFFIX));
411 if (outFd->HasArray(newArray->GetName()))
412 {
413 vtkWarningMacro(<< "Input has two arrays named " << array->GetName()
414 << ". Output statistics will probably be wrong.");
415 return;
416 }
417 outFd->AddArray(newArray);
418 }
419
420 if (this->ComputeMinimum)
421 {
422 vtkSmartPointer<vtkDataArray> newArray;
423 newArray.TakeReference(vtkDataArray::SafeDownCast(
424 vtkAbstractArray::CreateArray(array->GetDataType())));
425 newArray->DeepCopy(array);
426 newArray->SetName(vtkTemporalStatisticsMangleName(array->GetName(),
427 MINIMUM_SUFFIX));
428 outFd->AddArray(newArray);
429 }
430
431 if (this->ComputeMaximum)
432 {
433 vtkSmartPointer<vtkDataArray> newArray;
434 newArray.TakeReference(vtkDataArray::SafeDownCast(
435 vtkAbstractArray::CreateArray(array->GetDataType())));
436 newArray->DeepCopy(array);
437 newArray->SetName(vtkTemporalStatisticsMangleName(array->GetName(),
438 MAXIMUM_SUFFIX));
439 outFd->AddArray(newArray);
440 }
441
442 if (this->ComputeStandardDeviation)
443 {
444 vtkSmartPointer<vtkDataArray> newArray;
445 newArray.TakeReference(vtkDataArray::SafeDownCast(
446 vtkAbstractArray::CreateArray(array->GetDataType())));
447 newArray->SetName(vtkTemporalStatisticsMangleName(array->GetName(),
448 STANDARD_DEVIATION_SUFFIX));
449
450 newArray->SetNumberOfComponents(array->GetNumberOfComponents());
451 newArray->CopyComponentNames( array );
452
453 newArray->SetNumberOfTuples(array->GetNumberOfTuples());
454 switch (array->GetDataType())
455 {
456 vtkTemplateMacro(vtkTemporalStatisticsInitializeStdDev(
457 static_cast<VTK_TT*>(newArray->GetVoidPointer(0)),
458 array->GetNumberOfComponents()
459 *array->GetNumberOfTuples()));
460 }
461 outFd->AddArray(newArray);
462 }
463 }
464
465 //-----------------------------------------------------------------------------
AccumulateStatistics(vtkDataObject * input,vtkDataObject * output)466 void vtkTemporalStatistics::AccumulateStatistics(vtkDataObject *input,
467 vtkDataObject *output)
468 {
469 if (input->IsA("vtkDataSet"))
470 {
471 this->AccumulateStatistics(vtkDataSet::SafeDownCast(input),
472 vtkDataSet::SafeDownCast(output));
473 return;
474 }
475
476 if (input->IsA("vtkGraph"))
477 {
478 this->AccumulateStatistics(vtkGraph::SafeDownCast(input),
479 vtkGraph::SafeDownCast(output));
480 return;
481 }
482
483 if (input->IsA("vtkCompositeDataSet"))
484 {
485 this->AccumulateStatistics(vtkCompositeDataSet::SafeDownCast(input),
486 vtkCompositeDataSet::SafeDownCast(output));
487 }
488 }
489
490 //-----------------------------------------------------------------------------
AccumulateStatistics(vtkDataSet * input,vtkDataSet * output)491 void vtkTemporalStatistics::AccumulateStatistics(vtkDataSet *input,
492 vtkDataSet *output)
493 {
494 this->AccumulateArrays(input->GetFieldData(), output->GetFieldData());
495 this->AccumulateArrays(input->GetPointData(), output->GetPointData());
496 this->AccumulateArrays(input->GetCellData(), output->GetCellData());
497 }
498
499 //-----------------------------------------------------------------------------
AccumulateStatistics(vtkGraph * input,vtkGraph * output)500 void vtkTemporalStatistics::AccumulateStatistics(vtkGraph *input,
501 vtkGraph *output)
502 {
503 this->AccumulateArrays(input->GetFieldData(), output->GetFieldData());
504 this->AccumulateArrays(input->GetVertexData(), output->GetVertexData());
505 this->AccumulateArrays(input->GetEdgeData(), output->GetEdgeData());
506 }
507
508 //-----------------------------------------------------------------------------
AccumulateStatistics(vtkCompositeDataSet * input,vtkCompositeDataSet * output)509 void vtkTemporalStatistics::AccumulateStatistics(vtkCompositeDataSet *input,
510 vtkCompositeDataSet *output)
511 {
512 vtkSmartPointer<vtkCompositeDataIterator> inputItr;
513 inputItr.TakeReference(input->NewIterator());
514
515 for (inputItr->InitTraversal(); !inputItr->IsDoneWithTraversal();
516 inputItr->GoToNextItem())
517 {
518 vtkDataObject *inputObj = inputItr->GetCurrentDataObject();
519 vtkDataObject *outputObj = output->GetDataSet(inputItr);
520
521 this->AccumulateStatistics(inputObj, outputObj);
522 }
523 }
524
525 //-----------------------------------------------------------------------------
AccumulateArrays(vtkFieldData * inFd,vtkFieldData * outFd)526 void vtkTemporalStatistics::AccumulateArrays(vtkFieldData *inFd,
527 vtkFieldData *outFd)
528 {
529 int numArrays = inFd->GetNumberOfArrays();
530 for (int i = 0; i < numArrays; i++)
531 {
532 vtkDataArray *inArray = inFd->GetArray(i);
533 vtkDataArray *outArray;
534 if (!inArray) continue;
535
536 outArray = this->GetArray(outFd, inArray, AVERAGE_SUFFIX);
537 if (outArray)
538 {
539
540
541 vtkDataArray* stdevOutArray =
542 this->GetArray(outFd, inArray, STANDARD_DEVIATION_SUFFIX);
543 if (stdevOutArray)
544 {
545 switch (inArray->GetDataType())
546 {
547 // standard deviation must be called before average since the one-pass
548 // algorithm uses the average up to the previous time step
549 vtkTemplateMacro(vtkTemporalStatisticsAccumulateStdDev(
550 static_cast<const VTK_TT*>(inArray->GetVoidPointer(0)),
551 static_cast<VTK_TT*>(stdevOutArray->GetVoidPointer(0)),
552 static_cast<const VTK_TT*>(outArray->GetVoidPointer(0)),
553 inArray->GetNumberOfComponents()*inArray->GetNumberOfTuples(),
554 this->CurrentTimeIndex));
555 }
556 // Alert change in data.
557 stdevOutArray->DataChanged();
558 }
559
560
561
562
563 switch (inArray->GetDataType())
564 {
565 vtkTemplateMacro(vtkTemporalStatisticsAccumulateAverage(
566 static_cast<const VTK_TT*>(inArray->GetVoidPointer(0)),
567 static_cast<VTK_TT*>(outArray->GetVoidPointer(0)),
568 inArray->GetNumberOfComponents()
569 *inArray->GetNumberOfTuples()));
570 }
571 // Alert change in data.
572 outArray->DataChanged();
573 }
574
575 outArray = this->GetArray(outFd, inArray, MINIMUM_SUFFIX);
576 if (outArray)
577 {
578 switch (inArray->GetDataType())
579 {
580 vtkTemplateMacro(vtkTemporalStatisticsAccumulateMinimum(
581 static_cast<const VTK_TT*>(inArray->GetVoidPointer(0)),
582 static_cast<VTK_TT*>(outArray->GetVoidPointer(0)),
583 inArray->GetNumberOfComponents()
584 *inArray->GetNumberOfTuples()));
585 }
586 // Alert change in data.
587 outArray->DataChanged();
588 }
589
590 outArray = this->GetArray(outFd, inArray, MAXIMUM_SUFFIX);
591 if (outArray)
592 {
593 switch (inArray->GetDataType())
594 {
595 vtkTemplateMacro(vtkTemporalStatisticsAccumulateMaximum(
596 static_cast<const VTK_TT*>(inArray->GetVoidPointer(0)),
597 static_cast<VTK_TT*>(outArray->GetVoidPointer(0)),
598 inArray->GetNumberOfComponents()
599 *inArray->GetNumberOfTuples()));
600 }
601 // Alert change in data.
602 outArray->DataChanged();
603 }
604 }
605 }
606
607 //-----------------------------------------------------------------------------
PostExecute(vtkDataObject * input,vtkDataObject * output)608 void vtkTemporalStatistics::PostExecute(vtkDataObject *input,
609 vtkDataObject *output)
610 {
611 if (input->IsA("vtkDataSet"))
612 {
613 this->PostExecute(vtkDataSet::SafeDownCast(input),
614 vtkDataSet::SafeDownCast(output));
615 return;
616 }
617
618 if (input->IsA("vtkGraph"))
619 {
620 this->PostExecute(vtkGraph::SafeDownCast(input),
621 vtkGraph::SafeDownCast(output));
622 return;
623 }
624
625 if (input->IsA("vtkCompositeDataSet"))
626 {
627 this->PostExecute(vtkCompositeDataSet::SafeDownCast(input),
628 vtkCompositeDataSet::SafeDownCast(output));
629 }
630 }
631
632 //-----------------------------------------------------------------------------
PostExecute(vtkDataSet * input,vtkDataSet * output)633 void vtkTemporalStatistics::PostExecute(vtkDataSet *input, vtkDataSet *output)
634 {
635 this->FinishArrays(input->GetFieldData(), output->GetFieldData());
636 this->FinishArrays(input->GetPointData(), output->GetPointData());
637 this->FinishArrays(input->GetCellData(), output->GetCellData());
638 }
639
640 //-----------------------------------------------------------------------------
PostExecute(vtkGraph * input,vtkGraph * output)641 void vtkTemporalStatistics::PostExecute(vtkGraph *input, vtkGraph *output)
642 {
643 this->FinishArrays(input->GetFieldData(), output->GetFieldData());
644 this->FinishArrays(input->GetVertexData(), output->GetVertexData());
645 this->FinishArrays(input->GetEdgeData(), output->GetEdgeData());
646 }
647
648 //-----------------------------------------------------------------------------
PostExecute(vtkCompositeDataSet * input,vtkCompositeDataSet * output)649 void vtkTemporalStatistics::PostExecute(vtkCompositeDataSet *input,
650 vtkCompositeDataSet *output)
651 {
652 vtkSmartPointer<vtkCompositeDataIterator> inputItr;
653 inputItr.TakeReference(input->NewIterator());
654
655 for (inputItr->InitTraversal(); !inputItr->IsDoneWithTraversal();
656 inputItr->GoToNextItem())
657 {
658 vtkDataObject *inputObj = inputItr->GetCurrentDataObject();
659 vtkDataObject *outputObj = output->GetDataSet(inputItr);
660
661 this->PostExecute(inputObj, outputObj);
662 }
663 }
664
665 //-----------------------------------------------------------------------------
FinishArrays(vtkFieldData * inFd,vtkFieldData * outFd)666 void vtkTemporalStatistics::FinishArrays(vtkFieldData *inFd,
667 vtkFieldData *outFd)
668 {
669 int numArrays = inFd->GetNumberOfArrays();
670 for (int i = 0; i < numArrays; i++)
671 {
672 vtkDataArray *inArray = inFd->GetArray(i);
673 vtkDataArray *outArray;
674 if (!inArray) continue;
675
676 outArray = this->GetArray(outFd, inArray, AVERAGE_SUFFIX);
677 if (outArray)
678 {
679 switch (inArray->GetDataType())
680 {
681 vtkTemplateMacro(vtkTemporalStatisticsFinishAverage(
682 static_cast<VTK_TT*>(outArray->GetVoidPointer(0)),
683 inArray->GetNumberOfComponents()
684 *inArray->GetNumberOfTuples(),
685 this->CurrentTimeIndex));
686 }
687 }
688 vtkDataArray *avgArray = outArray;
689
690 // No post processing on minimum.
691 // No post processing on maximum.
692
693 outArray = this->GetArray(outFd, inArray, STANDARD_DEVIATION_SUFFIX);
694 if (outArray)
695 {
696 if (!avgArray)
697 {
698 vtkWarningMacro(<< "Average not computed for " << inArray->GetName()
699 << ", standard deviation skipped.");
700 outFd->RemoveArray(outArray->GetName());
701 }
702 else
703 {
704 switch (inArray->GetDataType())
705 {
706 vtkTemplateMacro(vtkTemporalStatisticsFinishStdDev(
707 static_cast<VTK_TT*>(outArray->GetVoidPointer(0)),
708 inArray->GetNumberOfComponents()
709 *inArray->GetNumberOfTuples(),
710 this->CurrentTimeIndex));
711 }
712 if (!this->ComputeAverage)
713 {
714 outFd->RemoveArray(avgArray->GetName());
715 }
716 }
717 }
718 }
719 }
720
721 //-----------------------------------------------------------------------------
GetArray(vtkFieldData * fieldData,vtkDataArray * inArray,const char * nameSuffix)722 vtkDataArray *vtkTemporalStatistics::GetArray(vtkFieldData *fieldData,
723 vtkDataArray *inArray,
724 const char *nameSuffix)
725 {
726 vtkStdString outArrayName
727 = vtkTemporalStatisticsMangleName(inArray->GetName(), nameSuffix);
728 vtkDataArray *outArray = fieldData->GetArray(outArrayName.c_str());
729 if (!outArray) return NULL;
730
731 if ( (inArray->GetNumberOfComponents() != outArray->GetNumberOfComponents())
732 || (inArray->GetNumberOfTuples() != outArray->GetNumberOfTuples()) )
733 {
734 vtkWarningMacro(<< "Size of array " << outArray->GetName()
735 << " has changed. Does the source change the topology "
736 << " over time?");
737 fieldData->RemoveArray(outArray->GetName());
738 return NULL;
739 }
740
741 return outArray;
742 }
743