1 //============================================================================
2 // Copyright (c) Kitware, Inc.
3 // All rights reserved.
4 // See LICENSE.txt for details.
5 //
6 // This software is distributed WITHOUT ANY WARRANTY; without even
7 // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
8 // PURPOSE. See the above copyright notice for more information.
9 //============================================================================
10
11 #ifndef vtk_m_worklet_FieldHistogram_h
12 #define vtk_m_worklet_FieldHistogram_h
13
14 #include <vtkm/Math.h>
15 #include <vtkm/cont/Algorithm.h>
16 #include <vtkm/cont/ArrayGetValues.h>
17 #include <vtkm/cont/ArrayHandle.h>
18 #include <vtkm/cont/ArrayHandleCounting.h>
19 #include <vtkm/worklet/DispatcherMapField.h>
20 #include <vtkm/worklet/WorkletMapField.h>
21
22 #include <vtkm/cont/Field.h>
23
24 namespace
25 {
26 // GCC creates false positive warnings for signed/unsigned char* operations.
27 // This occurs because the values are implicitly casted up to int's for the
28 // operation, and than casted back down to char's when return.
29 // This causes a false positive warning, even when the values is within
30 // the value types range
31 #if defined(VTKM_GCC)
32 #pragma GCC diagnostic push
33 #pragma GCC diagnostic ignored "-Wconversion"
34 #endif // gcc
35 template <typename T>
compute_delta(T fieldMinValue,T fieldMaxValue,vtkm::Id num)36 T compute_delta(T fieldMinValue, T fieldMaxValue, vtkm::Id num)
37 {
38 using VecType = vtkm::VecTraits<T>;
39 const T fieldRange = fieldMaxValue - fieldMinValue;
40 return fieldRange / static_cast<typename VecType::ComponentType>(num);
41 }
42 #if defined(VTKM_GCC)
43 #pragma GCC diagnostic pop
44 #endif // gcc
45 }
46
47 namespace vtkm
48 {
49 namespace worklet
50 {
51
52 //simple functor that prints basic statistics
53 class FieldHistogram
54 {
55 public:
56 // For each value set the bin it should be in
57 template <typename FieldType>
58 class SetHistogramBin : public vtkm::worklet::WorkletMapField
59 {
60 public:
61 using ControlSignature = void(FieldIn value, FieldOut binIndex);
62 using ExecutionSignature = void(_1, _2);
63 using InputDomain = _1;
64
65 vtkm::Id numberOfBins;
66 FieldType minValue;
67 FieldType delta;
68
69 VTKM_CONT
SetHistogramBin(vtkm::Id numberOfBins0,FieldType minValue0,FieldType delta0)70 SetHistogramBin(vtkm::Id numberOfBins0, FieldType minValue0, FieldType delta0)
71 : numberOfBins(numberOfBins0)
72 , minValue(minValue0)
73 , delta(delta0)
74 {
75 }
76
77 VTKM_EXEC
operator()78 void operator()(const FieldType& value, vtkm::Id& binIndex) const
79 {
80 binIndex = static_cast<vtkm::Id>((value - minValue) / delta);
81 if (binIndex < 0)
82 binIndex = 0;
83 else if (binIndex >= numberOfBins)
84 binIndex = numberOfBins - 1;
85 }
86 };
87
88 // Calculate the adjacent difference between values in ArrayHandle
89 class AdjacentDifference : public vtkm::worklet::WorkletMapField
90 {
91 public:
92 using ControlSignature = void(FieldIn inputIndex, WholeArrayIn counts, FieldOut outputCount);
93 using ExecutionSignature = void(_1, _2, _3);
94 using InputDomain = _1;
95
96 template <typename WholeArrayType>
operator()97 VTKM_EXEC void operator()(const vtkm::Id& index,
98 const WholeArrayType& counts,
99 vtkm::Id& difference) const
100 {
101 if (index == 0)
102 difference = counts.Get(index);
103 else
104 difference = counts.Get(index) - counts.Get(index - 1);
105 }
106 };
107
108 // Execute the histogram binning filter given data and number of bins
109 // Returns:
110 // min value of the bins
111 // delta/range of each bin
112 // number of values in each bin
113 template <typename FieldType, typename Storage>
Run(vtkm::cont::ArrayHandle<FieldType,Storage> fieldArray,vtkm::Id numberOfBins,vtkm::Range & rangeOfValues,FieldType & binDelta,vtkm::cont::ArrayHandle<vtkm::Id> & binArray)114 void Run(vtkm::cont::ArrayHandle<FieldType, Storage> fieldArray,
115 vtkm::Id numberOfBins,
116 vtkm::Range& rangeOfValues,
117 FieldType& binDelta,
118 vtkm::cont::ArrayHandle<vtkm::Id>& binArray)
119 {
120 const vtkm::Vec<FieldType, 2> initValue{ vtkm::cont::ArrayGetValue(0, fieldArray) };
121
122 vtkm::Vec<FieldType, 2> result =
123 vtkm::cont::Algorithm::Reduce(fieldArray, initValue, vtkm::MinAndMax<FieldType>());
124
125 this->Run(fieldArray, numberOfBins, result[0], result[1], binDelta, binArray);
126
127 //update the users data
128 rangeOfValues = vtkm::Range(result[0], result[1]);
129 }
130
131 // Execute the histogram binning filter given data and number of bins, min,
132 // max values.
133 // Returns:
134 // number of values in each bin
135 template <typename FieldType, typename Storage>
Run(vtkm::cont::ArrayHandle<FieldType,Storage> fieldArray,vtkm::Id numberOfBins,FieldType fieldMinValue,FieldType fieldMaxValue,FieldType & binDelta,vtkm::cont::ArrayHandle<vtkm::Id> & binArray)136 void Run(vtkm::cont::ArrayHandle<FieldType, Storage> fieldArray,
137 vtkm::Id numberOfBins,
138 FieldType fieldMinValue,
139 FieldType fieldMaxValue,
140 FieldType& binDelta,
141 vtkm::cont::ArrayHandle<vtkm::Id>& binArray)
142 {
143 const vtkm::Id numberOfValues = fieldArray.GetNumberOfValues();
144
145 const FieldType fieldDelta = compute_delta(fieldMinValue, fieldMaxValue, numberOfBins);
146
147 // Worklet fills in the bin belonging to each value
148 vtkm::cont::ArrayHandle<vtkm::Id> binIndex;
149 binIndex.Allocate(numberOfValues);
150
151 // Worklet to set the bin number for each data value
152 SetHistogramBin<FieldType> binWorklet(numberOfBins, fieldMinValue, fieldDelta);
153 vtkm::worklet::DispatcherMapField<SetHistogramBin<FieldType>> setHistogramBinDispatcher(
154 binWorklet);
155 setHistogramBinDispatcher.Invoke(fieldArray, binIndex);
156
157 // Sort the resulting bin array for counting
158 vtkm::cont::Algorithm::Sort(binIndex);
159
160 // Get the upper bound of each bin number
161 vtkm::cont::ArrayHandle<vtkm::Id> totalCount;
162 vtkm::cont::ArrayHandleCounting<vtkm::Id> binCounter(0, 1, numberOfBins);
163 vtkm::cont::Algorithm::UpperBounds(binIndex, binCounter, totalCount);
164
165 // Difference between adjacent items is the bin count
166 vtkm::worklet::DispatcherMapField<AdjacentDifference> dispatcher;
167 dispatcher.Invoke(binCounter, totalCount, binArray);
168
169 //update the users data
170 binDelta = fieldDelta;
171 }
172 };
173 }
174 } // namespace vtkm::worklet
175
176 #endif // vtk_m_worklet_FieldHistogram_h
177