HistogramAnalysisAdaptor.cxx 6.04 KB
Newer Older
1 2
#include "HistogramAnalysisAdaptor.h"

3 4 5 6 7 8 9 10 11 12
#include <vtkAOSDataArrayTemplate.h>
#include <vtkArrayDispatch.h>
#include <vtkCompositeDataIterator.h>
#include <vtkCompositeDataSet.h>
#include <vtkDataArray.h>
#include <vtkDataObject.h>
#include <vtkFieldData.h>
#include <vtkInsituDataAdaptor.h>
#include <vtkObjectFactory.h>
#include <vtkSmartPointer.h>
13

14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
#include <algorithm>
#include <vector>
namespace
{

// Private worker for histogram method. Computes the local histogram on
// array (passed to operator()). To be used with vtkArrayDispatch.
//
// Inputs:
// range: Global range of data
// bins: Number of histogram bins
// array: Local data.
//
// Outputs:
// Histogram: The histogram of the local data.
struct HistogramWorker
{
  const double *Range;
  int Bins;
  std::vector<unsigned int> Histogram;
  HistogramWorker(const double *range, int bins) : Range(range), Bins(bins) {}

  template <typename ArrayT>
  void operator()(ArrayT *array)
  {
    assert(array);
    assert(array->GetNumberOfComponents() == 1);

    typedef typename ArrayT::ValueType ValueType;
    ValueType width = static_cast<ValueType>((this->Range[1] - this->Range[0]) / this->Bins);
    ValueType min = static_cast<ValueType>(this->Range[0]);
    vtkIdType numTuples = array->GetNumberOfTuples();

    // + 1 to store val == max. These will be moved to this last bin before
    // returning (Avoids having to branch in the loop below);
    this->Histogram.resize(this->Bins + 1, 0);
    for (vtkIdType tIdx = 0; tIdx < numTuples; ++tIdx)
      {
      int bin = static_cast<int>((array->GetComponent(tIdx, 0) - min) / width);
      ++this->Histogram[bin];
      }

    // Merge the last two bins (the last is only when val == max)
    this->Histogram[this->Bins-1] += this->Histogram[this->Bins];
    this->Histogram.resize(this->Bins);
  }
};

class vtkHistogram
{
  double Range[2];
  HistogramWorker *Worker;
public:
    vtkHistogram()
      {
      this->Range[0] = VTK_DOUBLE_MAX;
      this->Range[1] = VTK_DOUBLE_MIN;
      this->Worker = NULL;
      }
    ~vtkHistogram()
      {
      delete this->Worker;
      }

    void AddRange(vtkDataArray* da)
      {
      if (da)
        {
        double crange[2];
        da->GetRange(crange);
        this->Range[0] = std::min(this->Range[0], crange[0]);
        this->Range[1] = std::max(this->Range[1], crange[1]);
        }
      }
    void PreCompute(MPI_Comm comm, int bins)
      {
      double g_range[2];
      // Find the global max/min
      MPI_Allreduce(&this->Range[0], &g_range[0], 1, MPI_DOUBLE, MPI_MIN, comm);
      MPI_Allreduce(&this->Range[1], &g_range[1], 1, MPI_DOUBLE, MPI_MAX, comm);
      this->Range[0] = g_range[0];
      this->Range[1] = g_range[1];
      this->Worker = new HistogramWorker(g_range, bins);
      }

    void Compute(vtkDataArray* da)
      {
      if (da)
        {
        vtkArrayDispatch::Dispatch::Execute(da, *this->Worker);
        }
      }

    void PostCompute(MPI_Comm comm, int bins, const char* name)
      {
      std::vector<unsigned int> gHist(bins, 0);
      MPI_Reduce(&this->Worker->Histogram[0], &gHist[0], bins, MPI_UNSIGNED, MPI_SUM, 0, comm);

      int rank;
      MPI_Comm_rank(comm, &rank);
      if (rank == 0)
        {
        std::cout << "Histogram '" << name << "' (VTK):\n";
        double width = (this->Range[1] - this->Range[0]) / bins;
        for (int i = 0; i < bins; ++i)
          {
          printf("  %f-%f: %d\n", this->Range[0] + i*width, this->Range[0] + (i+1)*width, gHist[i]);
          }
        }
      }
};

}
127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154

vtkStandardNewMacro(HistogramAnalysisAdaptor);
//-----------------------------------------------------------------------------
HistogramAnalysisAdaptor::HistogramAnalysisAdaptor() :
  Communicator(MPI_COMM_WORLD),
  Bins(0),
  Association(vtkDataObject::FIELD_ASSOCIATION_POINTS)
{
}

//-----------------------------------------------------------------------------
HistogramAnalysisAdaptor::~HistogramAnalysisAdaptor()
{
}

//-----------------------------------------------------------------------------
void HistogramAnalysisAdaptor::Initialize(
  MPI_Comm comm, int bins, int association, const std::string& arrayname)
{
  this->Communicator = comm;
  this->Bins = bins;
  this->ArrayName = arrayname;
  this->Association = association;
}

//-----------------------------------------------------------------------------
bool HistogramAnalysisAdaptor::Execute(vtkInsituDataAdaptor* data)
{
155
  vtkHistogram histogram;
156
  vtkDataObject* mesh = data->GetMesh(/*structure_only*/true);
157
  if (mesh == NULL || !data->AddArray(mesh, this->Association, this->ArrayName.c_str()))
158
    {
159 160 161 162 163 164 165 166 167 168
    histogram.PreCompute(this->Communicator, this->Bins);
    histogram.PostCompute(this->Communicator, this->Bins, this->ArrayName.c_str());
    return true;
    }

  if (vtkCompositeDataSet* cd = vtkCompositeDataSet::SafeDownCast(mesh))
    {
    vtkSmartPointer<vtkCompositeDataIterator> iter;
    iter.TakeReference(cd->NewIterator());
    for (iter->InitTraversal(); !iter->IsDoneWithTraversal(); iter->GoToNextItem())
169
      {
170 171
      vtkDataArray* array = this->GetArray(iter->GetCurrentDataObject());
      histogram.AddRange(array);
172
      }
173 174
    histogram.PreCompute(this->Communicator, this->Bins);
    for (iter->InitTraversal(); !iter->IsDoneWithTraversal(); iter->GoToNextItem())
175
      {
176 177
      vtkDataArray* array = this->GetArray(iter->GetCurrentDataObject());
      histogram.Compute(array);
178
      }
179
    histogram.PostCompute(this->Communicator, this->Bins, this->ArrayName.c_str());
180
    }
181 182 183 184 185 186 187 188 189
  else
    {
    vtkDataArray* array = this->GetArray(mesh);
    histogram.AddRange(array);
    histogram.PreCompute(this->Communicator, this->Bins);
    histogram.Compute(array);
    histogram.PostCompute(this->Communicator, this->Bins, this->ArrayName.c_str());
    }
  return true;
190 191 192
}

//-----------------------------------------------------------------------------
193
vtkDataArray* HistogramAnalysisAdaptor::GetArray(vtkDataObject* dobj)
194
{
195 196
  assert(dobj != NULL && vtkCompositeDataSet::SafeDownCast(dobj) == NULL);
  if (vtkFieldData* fd = dobj->GetAttributesAsFieldData(this->Association))
197
    {
198
    return fd->GetArray(this->ArrayName.c_str());
199
    }
200
  return NULL;
201
}