1 /* This source file is part of the Avogadro project.
2 It is released under the 3-Clause BSD License, see "LICENSE". */
3
4 #ifndef AVOGADRO_COMPUTE_HISTOGRAM_H
5 #define AVOGADRO_COMPUTE_HISTOGRAM_H
6
7 #include <vtkDoubleArray.h>
8 #include <vtkFloatArray.h>
9 #include <vtkImageData.h>
10 #include <vtkIntArray.h>
11 #include <vtkMath.h>
12 #include <vtkPointData.h>
13 #include <vtkTable.h>
14
15 #include <cmath>
16
17 namespace Avogadro {
18
19 /** Single component integral type specialization. */
20 template <typename T,
21 typename std::enable_if<std::is_integral<T>::value>::type* = nullptr>
calcHistogram(T * values,const vtkIdType numTuples,const float min,const float inv,int * pops,int &)22 void calcHistogram(T* values, const vtkIdType numTuples, const float min,
23 const float inv, int* pops, int&)
24 {
25 for (vtkIdType j = 0; j < numTuples; ++j) {
26 ++pops[static_cast<int>((*values++ - min) * inv)];
27 }
28 }
29
30 /** Needs to be present, should never be compiled. */
31 template <typename T>
calcHistogram(T *,const vtkIdType,int *)32 void calcHistogram(T*, const vtkIdType, int*)
33 {
34 static_assert(!std::is_same<unsigned char, T>::value, "Invalid type");
35 }
36
37 /** Single component unsigned char covering 0 -> 255 range. */
calcHistogram(unsigned char * values,const vtkIdType numTuples,int * pops)38 void calcHistogram(unsigned char* values, const vtkIdType numTuples, int* pops)
39 {
40 for (vtkIdType j = 0; j < numTuples; ++j) {
41 ++pops[*values++];
42 }
43 }
44
45 /** Single component floating point type specialization. */
46 template <typename T,
47 typename std::enable_if<!std::is_integral<T>::value>::type* = nullptr>
calcHistogram(T * values,const vtkIdType numTuples,const float min,const float inv,int * pops,int & invalid)48 void calcHistogram(T* values, const vtkIdType numTuples, const float min,
49 const float inv, int* pops, int& invalid)
50 {
51 for (vtkIdType j = 0; j < numTuples; ++j) {
52 T value = *(values++);
53 if (std::isfinite(value)) {
54 ++pops[static_cast<int>((value - min) * inv)];
55 } else {
56 ++invalid;
57 }
58 }
59 }
60
61 /**
62 * Computes a histogram from an array of values.
63 * \param values The array from which to compute the histogram.
64 * \param numTuples Number of tuples in the array.
65 * \param numComponents Number of components in each tuple.
66 * \param min Minimum value in range
67 * \param max Maximum value in range
68 * \param inv Inverse of bin size, numBins is the number of bins
69 * in the histogram (or length of the pops array), and invalid is a return
70 * parameter indicating how many values in the array had a non-finite value.
71 */
72 template <typename T>
CalculateHistogram(T * values,const vtkIdType numTuples,const vtkIdType numComponents,const float min,const float max,int * pops,const float inv,int & invalid)73 void CalculateHistogram(T* values, const vtkIdType numTuples,
74 const vtkIdType numComponents, const float min,
75 const float max, int* pops, const float inv,
76 int& invalid)
77 {
78 // Single component is a simpler/faster path, let's dispatch separately.
79 if (numComponents == 1) {
80 // Very fast path for unsigned char in 0 -> 255 range, or fast path.
81 if (std::is_same<T, unsigned char>::value && min == 0.f && max == 255.f) {
82 calcHistogram(values, numTuples, pops);
83 } else {
84 calcHistogram(values, numTuples, min, inv, pops, invalid);
85 }
86 } else {
87 // Multicomponent magnitude
88 for (vtkIdType j = 0; j < numTuples; ++j) {
89 // Check that all components are valid.
90 bool valid = true;
91 double squaredSum = 0.0;
92 for (vtkIdType c = 0; c < numComponents; ++c) {
93 T value = *(values + c);
94 if (!vtkMath::IsFinite(value)) {
95 valid = false;
96 break;
97 }
98 squaredSum += (value * value);
99 }
100 if (valid) {
101 int index = static_cast<int>((sqrt(squaredSum) - min) * inv);
102 ++pops[index];
103 } else {
104 ++invalid;
105 }
106 values += numComponents;
107 }
108 }
109 }
110
PopulateHistogram(vtkImageData * input,vtkTable * output)111 void PopulateHistogram(vtkImageData* input, vtkTable* output)
112 {
113 // The output table will have the twice the number of columns, they will be
114 // the x and y for input column. This is the bin centers, and the population.
115 double minmax[2] = { 0.0, 0.0 };
116
117 // This number of bins in the 2D histogram will also be used as the number of
118 // bins in the 2D transfer function for X (scalar value) and Y (gradient mag.)
119 const int numberOfBins = 256;
120
121 // Keep the array we are working on around even if the user shallow copies
122 // over the input image data by incrementing the reference count here.
123 vtkSmartPointer<vtkDataArray> arrayPtr = input->GetPointData()->GetScalars();
124 if (!arrayPtr) {
125 return;
126 }
127
128 // The bin values are the centers, extending +/- half an inc either side
129 arrayPtr->GetFiniteRange(minmax, -1);
130 if (minmax[0] == minmax[1]) {
131 minmax[1] = minmax[0] + 1.0;
132 }
133
134 double inc = (minmax[1] - minmax[0]) / (numberOfBins - 1);
135 double halfInc = inc / 2.0;
136 vtkSmartPointer<vtkFloatArray> extents =
137 vtkFloatArray::SafeDownCast(output->GetColumnByName("image_extents"));
138 if (!extents) {
139 extents = vtkSmartPointer<vtkFloatArray>::New();
140 extents->SetName("image_extents");
141 }
142 extents->SetNumberOfTuples(numberOfBins);
143 double min = minmax[0] + halfInc;
144 for (int j = 0; j < numberOfBins; ++j) {
145 extents->SetValue(j, min + j * inc);
146 }
147 vtkSmartPointer<vtkIntArray> populations =
148 vtkIntArray::SafeDownCast(output->GetColumnByName("image_pops"));
149 if (!populations) {
150 populations = vtkSmartPointer<vtkIntArray>::New();
151 populations->SetName("image_pops");
152 }
153 populations->SetNumberOfTuples(numberOfBins);
154 auto pops = static_cast<int*>(populations->GetVoidPointer(0));
155 for (int k = 0; k < numberOfBins; ++k) {
156 pops[k] = 0;
157 }
158 int invalid = 0;
159
160 switch (arrayPtr->GetDataType()) {
161 vtkTemplateMacro(CalculateHistogram(
162 reinterpret_cast<VTK_TT*>(arrayPtr->GetVoidPointer(0)),
163 arrayPtr->GetNumberOfTuples(), arrayPtr->GetNumberOfComponents(),
164 minmax[0], minmax[1], pops, 1.0 / inc, invalid));
165 default:
166 cout << "UpdateFromFile: Unknown data type" << endl;
167 }
168
169 #ifndef NDEBUG
170 vtkIdType total = invalid;
171 for (int i = 0; i < numberOfBins; ++i)
172 total += pops[i];
173 assert(total == arrayPtr->GetNumberOfTuples());
174 #endif
175 if (invalid) {
176 cout << "Warning: NaN or infinite value in dataset" << endl;
177 }
178
179 output->AddColumn(extents);
180 output->AddColumn(populations);
181 }
182
183 } // namespace Avogadro
184 #endif // AVOGADRO_COMPUTE_HISTOGRAM_H
185