// -*- c++ -*- /*========================================================================= Program: Visualization Toolkit Module: vtkTemporalStatistics.cxx Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen All rights reserved. See Copyright.txt or http://www.kitware.com/Copyright.htm for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notice for more information. =========================================================================*/ /* * Copyright 2008 Sandia Corporation. * Under the terms of Contract DE-AC04-94AL85000, there is a non-exclusive * license for use of this work by or on behalf of the * U.S. Government. Redistribution and use in source and binary forms, with * or without modification, are permitted provided that this Notice and any * statement of authorship are reproduced on all copies. */ #include "vtkTemporalStatistics.h" #include "vtkCellData.h" #include "vtkCompositeDataIterator.h" #include "vtkCompositeDataSet.h" #include "vtkDataArray.h" #include "vtkDataSet.h" #include "vtkGraph.h" #include "vtkInformation.h" #include "vtkInformationVector.h" #include "vtkMultiBlockDataSet.h" #include "vtkObjectFactory.h" #include "vtkPointData.h" #include "vtkStdString.h" #include "vtkStreamingDemandDrivenPipeline.h" #include "vtkSmartPointer.h" #include //============================================================================= vtkStandardNewMacro(vtkTemporalStatistics); //============================================================================= const char * const AVERAGE_SUFFIX = "average"; const char * const MINIMUM_SUFFIX = "minimum"; const char * const MAXIMUM_SUFFIX = "maximum"; const char * const STANDARD_DEVIATION_SUFFIX = "stddev"; inline vtkStdString vtkTemporalStatisticsMangleName(const char *originalName, const char *suffix) { if (!originalName) return suffix; return vtkStdString(originalName) + "_" + vtkStdString(suffix); } //----------------------------------------------------------------------------- // The interm stddev array keeps a sum of squares. template inline void vtkTemporalStatisticsInitializeStdDev(T *outArray, vtkIdType arraySize) { for (vtkIdType i = 0; i < arraySize; i++) { outArray[i] = 0; } } //----------------------------------------------------------------------------- template inline void vtkTemporalStatisticsAccumulateAverage(const T *inArray, T *outArray, vtkIdType arraySize) { for (vtkIdType i = 0; i < arraySize; i++) { outArray[i] += inArray[i]; } } template inline void vtkTemporalStatisticsAccumulateMinimum(const T *inArray, T *outArray, vtkIdType arraySize) { for (vtkIdType i = 0; i < arraySize; i++) { if (outArray[i] > inArray[i]) outArray[i] = inArray[i]; } } template inline void vtkTemporalStatisticsAccumulateMaximum(const T *inArray, T *outArray, vtkIdType arraySize) { for (vtkIdType i = 0; i < arraySize; i++) { if (outArray[i] < inArray[i]) outArray[i] = inArray[i]; } } // standard deviation one-pass algorithm from // http://www.cs.berkeley.edu/~mhoemmen/cs194/Tutorials/variance.pdf // this is numerically stable! template inline void vtkTemporalStatisticsAccumulateStdDev( const T *inArray, T *outArray, const T *previousAverage, vtkIdType arraySize, int pass) { for (vtkIdType i = 0; i < arraySize; i++) { double temp = inArray[i]-previousAverage[i]/static_cast(pass); outArray[i] = outArray[i] + static_cast( pass*temp*temp/static_cast(pass+1) ); } } //----------------------------------------------------------------------------- template inline void vtkTemporalStatisticsFinishAverage(T *outArray, vtkIdType arraySize, int sumSize) { for (vtkIdType i = 0; i < arraySize; i++) { outArray[i] /= sumSize; } } template inline void vtkTemporalStatisticsFinishStdDev(T *outArray, vtkIdType arraySize, int sumSize) { for (vtkIdType i = 0; i < arraySize; i++) { outArray[i] = static_cast(sqrt(static_cast(outArray[i])/sumSize)); } } //============================================================================= vtkTemporalStatistics::vtkTemporalStatistics() { this->ComputeAverage = 1; this->ComputeMinimum = 1; this->ComputeMaximum = 1; this->ComputeStandardDeviation = 1; this->CurrentTimeIndex = 0; } vtkTemporalStatistics::~vtkTemporalStatistics() { } void vtkTemporalStatistics::PrintSelf(ostream &os, vtkIndent indent) { this->Superclass::PrintSelf(os, indent); os << indent << "ComputeAverage: " << this->ComputeAverage << endl; os << indent << "ComputeMinimum: " << this->ComputeMinimum << endl; os << indent << "ComputeMaximum: " << this->ComputeMaximum << endl; os << indent << "ComputeStandardDeviation: " << this->ComputeStandardDeviation << endl; } //----------------------------------------------------------------------------- int vtkTemporalStatistics::FillInputPortInformation(int vtkNotUsed(port), vtkInformation *info) { info->Remove(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE()); info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet"); info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkGraph"); info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkCompositeDataSet"); return 1; } //----------------------------------------------------------------------------- int vtkTemporalStatistics::RequestInformation( vtkInformation *vtkNotUsed(request), vtkInformationVector **vtkNotUsed(inputVector), vtkInformationVector *outputVector) { vtkInformation *outInfo = outputVector->GetInformationObject(0); // The output data of this filter has no time assoicated with it. It is the // result of computations that happen over all time. outInfo->Remove(vtkStreamingDemandDrivenPipeline::TIME_STEPS()); outInfo->Remove(vtkStreamingDemandDrivenPipeline::TIME_RANGE()); return 1; } //----------------------------------------------------------------------------- int vtkTemporalStatistics::RequestDataObject( vtkInformation *vtkNotUsed(request), vtkInformationVector **inputVector, vtkInformationVector *outputVector) { vtkInformation *inInfo = inputVector[0]->GetInformationObject(0); vtkInformation *outInfo = outputVector->GetInformationObject(0); if (!inInfo) return 0; vtkDataObject *input = vtkDataObject::GetData(inInfo); vtkDataObject *output = vtkDataObject::GetData(outInfo); if (!input) return 0; vtkSmartPointer newOutput; if (input->IsA("vtkTemporalDataSet")) { // Special case: there are multiple problems with outputing a temporal data // set. The first is that this class removes all temporal-ness and // outputing a temporal data set makes no sense. The second is that the // pipeline has a nasty habit of automatically changing the output to // vtkDataObject. The easiest thing to do is just copy the output into a // multi block data set instead. if (!output || !output->IsA("vtkCompositeDataSet")) { newOutput = vtkSmartPointer::New(); } } else { if (!output || !output->IsA(input->GetClassName())) { newOutput.TakeReference(input->NewInstance()); } } if (newOutput) { newOutput->SetPipelineInformation(outInfo); } return 1; } //----------------------------------------------------------------------------- int vtkTemporalStatistics::RequestUpdateExtent( vtkInformation *vtkNotUsed(request), vtkInformationVector **inputVector, vtkInformationVector *vtkNotUsed(outputVector)) { vtkInformation *inInfo = inputVector[0]->GetInformationObject(0); // The RequestData method will tell the pipeline executive to iterate the // upstream pipeline to get each time step in order. The executive in turn // will call this method to get the extent request for each iteration (in this // case the time step). double *inTimes = inInfo->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS()); if (inTimes) { inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEPS(), &inTimes[this->CurrentTimeIndex], 1); } return 1; } //----------------------------------------------------------------------------- int vtkTemporalStatistics::RequestData(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector) { vtkInformation *inInfo = inputVector[0]->GetInformationObject(0); vtkInformation *outInfo = outputVector->GetInformationObject(0); vtkDataObject *input = vtkDataObject::GetData(inInfo); vtkDataObject *output = vtkDataObject::GetData(outInfo); if (this->CurrentTimeIndex == 0) { // First execution, initialize arrays. this->InitializeStatistics(input, output); } else { // Subsequent execution, accumulate new data. this->AccumulateStatistics(input, output); } this->CurrentTimeIndex++; if ( this->CurrentTimeIndex < inInfo->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS())) { // There is still more to do. request->Set(vtkStreamingDemandDrivenPipeline::CONTINUE_EXECUTING(), 1); } else { // We are done. Finish up. this->PostExecute(input, output); request->Remove(vtkStreamingDemandDrivenPipeline::CONTINUE_EXECUTING()); this->CurrentTimeIndex = 0; } return 1; } //----------------------------------------------------------------------------- void vtkTemporalStatistics::InitializeStatistics(vtkDataObject *input, vtkDataObject *output) { if (input->IsA("vtkDataSet")) { this->InitializeStatistics(vtkDataSet::SafeDownCast(input), vtkDataSet::SafeDownCast(output)); return; } if (input->IsA("vtkGraph")) { this->InitializeStatistics(vtkGraph::SafeDownCast(input), vtkGraph::SafeDownCast(output)); return; } if (input->IsA("vtkCompositeDataSet")) { this->InitializeStatistics(vtkCompositeDataSet::SafeDownCast(input), vtkCompositeDataSet::SafeDownCast(output)); return; } vtkWarningMacro(<< "Unsupported input type: " << input->GetClassName()); } //----------------------------------------------------------------------------- void vtkTemporalStatistics::InitializeStatistics(vtkDataSet *input, vtkDataSet *output) { output->CopyStructure(input); this->InitializeArrays(input->GetFieldData(), output->GetFieldData()); this->InitializeArrays(input->GetPointData(), output->GetPointData()); this->InitializeArrays(input->GetCellData(), output->GetCellData()); } //----------------------------------------------------------------------------- void vtkTemporalStatistics::InitializeStatistics(vtkGraph *input, vtkGraph *output) { output->CopyStructure(input); this->InitializeArrays(input->GetFieldData(), output->GetFieldData()); this->InitializeArrays(input->GetVertexData(), output->GetVertexData()); this->InitializeArrays(input->GetEdgeData(), output->GetEdgeData()); } //----------------------------------------------------------------------------- void vtkTemporalStatistics::InitializeStatistics(vtkCompositeDataSet *input, vtkCompositeDataSet *output) { output->CopyStructure(input); vtkSmartPointer inputItr; inputItr.TakeReference(input->NewIterator()); for (inputItr->InitTraversal(); !inputItr->IsDoneWithTraversal(); inputItr->GoToNextItem()) { vtkDataObject *inputObj = inputItr->GetCurrentDataObject(); vtkSmartPointer outputObj; outputObj.TakeReference(inputObj->NewInstance()); output->SetDataSet(inputItr, outputObj); this->InitializeStatistics(inputObj, outputObj); } } //----------------------------------------------------------------------------- void vtkTemporalStatistics::InitializeArrays(vtkFieldData *inFd, vtkFieldData *outFd) { // Because we need to do mathematical operations, we require all arrays we // process to be numeric data (i.e. a vtkDataArray). We also handle global // ids and petigree ids special (we just pass them). Ideally would just let // vtkFieldData or vtkDataSetAttributes handle this for us, but no such method // that fits our needs here. Thus, we pass data a bit differently then other // filters. If I miss something important, it should be added here. outFd->Initialize(); vtkDataSetAttributes *inDsa = vtkDataSetAttributes::SafeDownCast(inFd); vtkDataSetAttributes *outDsa = vtkDataSetAttributes::SafeDownCast(outFd); if (inDsa) { vtkDataArray *globalIds = inDsa->GetGlobalIds(); vtkAbstractArray *pedigreeIds = inDsa->GetPedigreeIds(); if (globalIds) outDsa->SetGlobalIds(globalIds); if (pedigreeIds) outDsa->SetPedigreeIds(pedigreeIds); } int numArrays = inFd->GetNumberOfArrays(); for (int i = 0; i < numArrays; i++) { vtkDataArray *array = inFd->GetArray(i); if (!array) continue; // Array not numeric. if (outFd->HasArray(array->GetName())) continue; // Must be Ids. this->InitializeArray(array, outFd); } } //----------------------------------------------------------------------------- void vtkTemporalStatistics::InitializeArray(vtkDataArray *array, vtkFieldData *outFd) { if (this->ComputeAverage || this->ComputeStandardDeviation) { vtkSmartPointer newArray; newArray.TakeReference(vtkDataArray::SafeDownCast( vtkAbstractArray::CreateArray(array->GetDataType()))); newArray->DeepCopy(array); newArray->SetName(vtkTemporalStatisticsMangleName(array->GetName(), AVERAGE_SUFFIX)); if (outFd->HasArray(newArray->GetName())) { vtkWarningMacro(<< "Input has two arrays named " << array->GetName() << ". Output statistics will probably be wrong."); return; } outFd->AddArray(newArray); } if (this->ComputeMinimum) { vtkSmartPointer newArray; newArray.TakeReference(vtkDataArray::SafeDownCast( vtkAbstractArray::CreateArray(array->GetDataType()))); newArray->DeepCopy(array); newArray->SetName(vtkTemporalStatisticsMangleName(array->GetName(), MINIMUM_SUFFIX)); outFd->AddArray(newArray); } if (this->ComputeMaximum) { vtkSmartPointer newArray; newArray.TakeReference(vtkDataArray::SafeDownCast( vtkAbstractArray::CreateArray(array->GetDataType()))); newArray->DeepCopy(array); newArray->SetName(vtkTemporalStatisticsMangleName(array->GetName(), MAXIMUM_SUFFIX)); outFd->AddArray(newArray); } if (this->ComputeStandardDeviation) { vtkSmartPointer newArray; newArray.TakeReference(vtkDataArray::SafeDownCast( vtkAbstractArray::CreateArray(array->GetDataType()))); newArray->SetName(vtkTemporalStatisticsMangleName(array->GetName(), STANDARD_DEVIATION_SUFFIX)); newArray->SetNumberOfComponents(array->GetNumberOfComponents()); newArray->CopyComponentNames( array ); newArray->SetNumberOfTuples(array->GetNumberOfTuples()); switch (array->GetDataType()) { vtkTemplateMacro(vtkTemporalStatisticsInitializeStdDev( static_cast(newArray->GetVoidPointer(0)), array->GetNumberOfComponents() *array->GetNumberOfTuples())); } outFd->AddArray(newArray); } } //----------------------------------------------------------------------------- void vtkTemporalStatistics::AccumulateStatistics(vtkDataObject *input, vtkDataObject *output) { if (input->IsA("vtkDataSet")) { this->AccumulateStatistics(vtkDataSet::SafeDownCast(input), vtkDataSet::SafeDownCast(output)); return; } if (input->IsA("vtkGraph")) { this->AccumulateStatistics(vtkGraph::SafeDownCast(input), vtkGraph::SafeDownCast(output)); return; } if (input->IsA("vtkCompositeDataSet")) { this->AccumulateStatistics(vtkCompositeDataSet::SafeDownCast(input), vtkCompositeDataSet::SafeDownCast(output)); } } //----------------------------------------------------------------------------- void vtkTemporalStatistics::AccumulateStatistics(vtkDataSet *input, vtkDataSet *output) { this->AccumulateArrays(input->GetFieldData(), output->GetFieldData()); this->AccumulateArrays(input->GetPointData(), output->GetPointData()); this->AccumulateArrays(input->GetCellData(), output->GetCellData()); } //----------------------------------------------------------------------------- void vtkTemporalStatistics::AccumulateStatistics(vtkGraph *input, vtkGraph *output) { this->AccumulateArrays(input->GetFieldData(), output->GetFieldData()); this->AccumulateArrays(input->GetVertexData(), output->GetVertexData()); this->AccumulateArrays(input->GetEdgeData(), output->GetEdgeData()); } //----------------------------------------------------------------------------- void vtkTemporalStatistics::AccumulateStatistics(vtkCompositeDataSet *input, vtkCompositeDataSet *output) { vtkSmartPointer inputItr; inputItr.TakeReference(input->NewIterator()); for (inputItr->InitTraversal(); !inputItr->IsDoneWithTraversal(); inputItr->GoToNextItem()) { vtkDataObject *inputObj = inputItr->GetCurrentDataObject(); vtkDataObject *outputObj = output->GetDataSet(inputItr); this->AccumulateStatistics(inputObj, outputObj); } } //----------------------------------------------------------------------------- void vtkTemporalStatistics::AccumulateArrays(vtkFieldData *inFd, vtkFieldData *outFd) { int numArrays = inFd->GetNumberOfArrays(); for (int i = 0; i < numArrays; i++) { vtkDataArray *inArray = inFd->GetArray(i); vtkDataArray *outArray; if (!inArray) continue; outArray = this->GetArray(outFd, inArray, AVERAGE_SUFFIX); if (outArray) { vtkDataArray* stdevOutArray = this->GetArray(outFd, inArray, STANDARD_DEVIATION_SUFFIX); if (stdevOutArray) { switch (inArray->GetDataType()) { // standard deviation must be called before average since the one-pass // algorithm uses the average up to the previous time step vtkTemplateMacro(vtkTemporalStatisticsAccumulateStdDev( static_cast(inArray->GetVoidPointer(0)), static_cast(stdevOutArray->GetVoidPointer(0)), static_cast(outArray->GetVoidPointer(0)), inArray->GetNumberOfComponents()*inArray->GetNumberOfTuples(), this->CurrentTimeIndex)); } // Alert change in data. stdevOutArray->DataChanged(); } switch (inArray->GetDataType()) { vtkTemplateMacro(vtkTemporalStatisticsAccumulateAverage( static_cast(inArray->GetVoidPointer(0)), static_cast(outArray->GetVoidPointer(0)), inArray->GetNumberOfComponents() *inArray->GetNumberOfTuples())); } // Alert change in data. outArray->DataChanged(); } outArray = this->GetArray(outFd, inArray, MINIMUM_SUFFIX); if (outArray) { switch (inArray->GetDataType()) { vtkTemplateMacro(vtkTemporalStatisticsAccumulateMinimum( static_cast(inArray->GetVoidPointer(0)), static_cast(outArray->GetVoidPointer(0)), inArray->GetNumberOfComponents() *inArray->GetNumberOfTuples())); } // Alert change in data. outArray->DataChanged(); } outArray = this->GetArray(outFd, inArray, MAXIMUM_SUFFIX); if (outArray) { switch (inArray->GetDataType()) { vtkTemplateMacro(vtkTemporalStatisticsAccumulateMaximum( static_cast(inArray->GetVoidPointer(0)), static_cast(outArray->GetVoidPointer(0)), inArray->GetNumberOfComponents() *inArray->GetNumberOfTuples())); } // Alert change in data. outArray->DataChanged(); } } } //----------------------------------------------------------------------------- void vtkTemporalStatistics::PostExecute(vtkDataObject *input, vtkDataObject *output) { if (input->IsA("vtkDataSet")) { this->PostExecute(vtkDataSet::SafeDownCast(input), vtkDataSet::SafeDownCast(output)); return; } if (input->IsA("vtkGraph")) { this->PostExecute(vtkGraph::SafeDownCast(input), vtkGraph::SafeDownCast(output)); return; } if (input->IsA("vtkCompositeDataSet")) { this->PostExecute(vtkCompositeDataSet::SafeDownCast(input), vtkCompositeDataSet::SafeDownCast(output)); } } //----------------------------------------------------------------------------- void vtkTemporalStatistics::PostExecute(vtkDataSet *input, vtkDataSet *output) { this->FinishArrays(input->GetFieldData(), output->GetFieldData()); this->FinishArrays(input->GetPointData(), output->GetPointData()); this->FinishArrays(input->GetCellData(), output->GetCellData()); } //----------------------------------------------------------------------------- void vtkTemporalStatistics::PostExecute(vtkGraph *input, vtkGraph *output) { this->FinishArrays(input->GetFieldData(), output->GetFieldData()); this->FinishArrays(input->GetVertexData(), output->GetVertexData()); this->FinishArrays(input->GetEdgeData(), output->GetEdgeData()); } //----------------------------------------------------------------------------- void vtkTemporalStatistics::PostExecute(vtkCompositeDataSet *input, vtkCompositeDataSet *output) { vtkSmartPointer inputItr; inputItr.TakeReference(input->NewIterator()); for (inputItr->InitTraversal(); !inputItr->IsDoneWithTraversal(); inputItr->GoToNextItem()) { vtkDataObject *inputObj = inputItr->GetCurrentDataObject(); vtkDataObject *outputObj = output->GetDataSet(inputItr); this->PostExecute(inputObj, outputObj); } } //----------------------------------------------------------------------------- void vtkTemporalStatistics::FinishArrays(vtkFieldData *inFd, vtkFieldData *outFd) { int numArrays = inFd->GetNumberOfArrays(); for (int i = 0; i < numArrays; i++) { vtkDataArray *inArray = inFd->GetArray(i); vtkDataArray *outArray; if (!inArray) continue; outArray = this->GetArray(outFd, inArray, AVERAGE_SUFFIX); if (outArray) { switch (inArray->GetDataType()) { vtkTemplateMacro(vtkTemporalStatisticsFinishAverage( static_cast(outArray->GetVoidPointer(0)), inArray->GetNumberOfComponents() *inArray->GetNumberOfTuples(), this->CurrentTimeIndex)); } } vtkDataArray *avgArray = outArray; // No post processing on minimum. // No post processing on maximum. outArray = this->GetArray(outFd, inArray, STANDARD_DEVIATION_SUFFIX); if (outArray) { if (!avgArray) { vtkWarningMacro(<< "Average not computed for " << inArray->GetName() << ", standard deviation skipped."); outFd->RemoveArray(outArray->GetName()); } else { switch (inArray->GetDataType()) { vtkTemplateMacro(vtkTemporalStatisticsFinishStdDev( static_cast(outArray->GetVoidPointer(0)), inArray->GetNumberOfComponents() *inArray->GetNumberOfTuples(), this->CurrentTimeIndex)); } if (!this->ComputeAverage) { outFd->RemoveArray(avgArray->GetName()); } } } } } //----------------------------------------------------------------------------- vtkDataArray *vtkTemporalStatistics::GetArray(vtkFieldData *fieldData, vtkDataArray *inArray, const char *nameSuffix) { vtkStdString outArrayName = vtkTemporalStatisticsMangleName(inArray->GetName(), nameSuffix); vtkDataArray *outArray = fieldData->GetArray(outArrayName.c_str()); if (!outArray) return NULL; if ( (inArray->GetNumberOfComponents() != outArray->GetNumberOfComponents()) || (inArray->GetNumberOfTuples() != outArray->GetNumberOfTuples()) ) { vtkWarningMacro(<< "Size of array " << outArray->GetName() << " has changed. Does the source change the topology " << " over time?"); fieldData->RemoveArray(outArray->GetName()); return NULL; } return outArray; }