NumCpp  2.5.1
A Templatized Header Only C++ Implementation of the Python NumPy Library
generateThreshold.hpp
Go to the documentation of this file.
1 
29 #pragma once
30 
34 #include "NumCpp/Core/Types.hpp"
35 #include "NumCpp/NdArray.hpp"
37 
38 #include <cmath>
39 #include <string>
40 
41 namespace nc
42 {
43  namespace imageProcessing
44  {
45  //============================================================================
46  // Method Description:
56  template<typename dtype>
57  dtype generateThreshold(const NdArray<dtype>& inImageArray, double inRate)
58  {
60 
61  if (inRate < 0.0 || inRate > 1.0)
62  {
63  THROW_INVALID_ARGUMENT_ERROR("input rate must be of the range [0, 1]");
64  }
65 
66  // first build a histogram
67  auto minValue = static_cast<int32>(std::floor(inImageArray.min().item()));
68  auto maxValue = static_cast<int32>(std::floor(inImageArray.max().item()));
69 
70  if (utils::essentiallyEqual(inRate, 0.0))
71  {
72  return static_cast<dtype>(maxValue);
73  }
74 
75  if (utils::essentiallyEqual(inRate, 1.0))
76  {
78  {
79  return static_cast<dtype>(minValue - 1);
80  }
81 
82  return dtype{ 0 };
83  }
84 
85  const auto histSize = static_cast<uint32>(maxValue - minValue + 1);
86 
87  NdArray<double> histogram(1, histSize);
88  histogram.zeros();
89  for (auto intensity : inImageArray)
90  {
91  const auto bin = static_cast<uint32>(static_cast<int32>(std::floor(intensity)) - minValue);
92  ++histogram[bin];
93  }
94 
95  // integrate the normalized histogram from right to left to make a survival function (1 - CDF)
96  const auto dNumPixels = static_cast<double>(inImageArray.size());
97  NdArray<double> survivalFunction(1, histSize + 1);
98  survivalFunction[-1] = 0.0;
99  for (int32 i = histSize - 1; i > -1; --i)
100  {
101  double histValue = histogram[i] / dNumPixels;
102  survivalFunction[i] = survivalFunction[i + 1] + histValue;
103  }
104 
105  // binary search through the survival function to find the rate
106  uint32 indexLow = 0;
107  uint32 indexHigh = histSize - 1;
108  uint32 index = indexHigh / 2; // integer division
109 
110  constexpr bool keepGoing = true;
111  while (keepGoing)
112  {
113  const double value = survivalFunction[index];
114  if (value < inRate)
115  {
116  indexHigh = index;
117  }
118  else if (value > inRate)
119  {
120  indexLow = index;
121  }
122  else
123  {
124  const int32 thresh = static_cast<int32>(index) + minValue - 1;
126  {
127  return static_cast<dtype>(thresh);
128  }
129 
130  return thresh < 0 ? 0 : static_cast<dtype>(thresh);
131  }
132 
133  if (indexHigh - indexLow < 2)
134  {
135  return static_cast<dtype>(static_cast<int32>(indexHigh) + minValue - 1);
136  }
137 
138  index = indexLow + (indexHigh - indexLow) / 2;
139  }
140 
141  // shouldn't ever get here but stop the compiler from throwing a warning
142  return static_cast<dtype>(histSize - 1);
143  }
144 
145  } // namespace imageProcessing
146 } // namespace nc
#define THROW_INVALID_ARGUMENT_ERROR(msg)
Definition: Error.hpp:36
#define STATIC_ASSERT_ARITHMETIC(dtype)
Definition: StaticAsserts.hpp:37
Holds info about the dtype.
Definition: DtypeInfo.hpp:41
Holds 1D and 2D arrays, the main work horse of the NumCpp library.
Definition: NdArrayCore.hpp:72
size_type size() const noexcept
Definition: NdArrayCore.hpp:4497
NdArray< dtype > min(Axis inAxis=Axis::NONE) const
Definition: NdArrayCore.hpp:3180
NdArray< dtype > max(Axis inAxis=Axis::NONE) const
Definition: NdArrayCore.hpp:3123
dtype generateThreshold(const NdArray< dtype > &inImageArray, double inRate)
Definition: generateThreshold.hpp:57
bool essentiallyEqual(dtype inValue1, dtype inValue2) noexcept
Definition: essentiallyEqual.hpp:52
Definition: Coordinate.hpp:45
dtype floor(dtype inValue) noexcept
Definition: floor.hpp:48
std::int32_t int32
Definition: Types.hpp:36
NdArray< uint32 > histogram(const NdArray< dtype > &inArray, const NdArray< double > &inBinEdges)
Definition: histogram.hpp:58
std::uint32_t uint32
Definition: Types.hpp:40