SourceXtractorPlusPlus  0.15
Please provide a description of the project.
ImageMode.cpp
Go to the documentation of this file.
1 
18 #include <Histogram/Histogram.h> // From Alexandria
19 
23 
24 
26 
27 namespace SourceXtractor {
28 
29 template<typename T>
31  int cell_w, int cell_h,
32  T invalid_value, T kappa1, T kappa2, T kappa3,
33  T rtol, size_t max_iter): m_image(image),
34  m_cell_w(cell_w), m_cell_h(cell_h),
35  m_invalid(invalid_value),
36  m_kappa1(kappa1), m_kappa2(kappa2), m_kappa3(kappa3),
37  m_rtol(rtol), m_max_iter(max_iter) {
38  auto hist_width = std::div(image->getWidth(), m_cell_w);
39  if (hist_width.rem)
40  ++hist_width.quot;
41  auto hist_height = std::div(image->getHeight(), m_cell_h);
42  if (hist_height.rem)
43  ++hist_height.quot;
44 
45  m_mode = VectorImage<T>::create(hist_width.quot, hist_height.quot);
46  m_sigma = VectorImage<T>::create(hist_width.quot, hist_height.quot);
47 
48  if (variance) {
49  m_var_mode = VectorImage<T>::create(hist_width.quot, hist_height.quot);
50  m_var_sigma = VectorImage<T>::create(hist_width.quot, hist_height.quot);
51 
52  for (int y = 0; y < m_mode->getHeight(); ++y) {
53  for (int x = 0; x < m_mode->getWidth(); ++x) {
54  processCell(*image, x, y, *m_mode, *m_sigma);
55  processCell(*variance, x, y, *m_var_mode, *m_var_sigma);
56  }
57  }
58  }
59  else {
60  for (int y = 0; y < m_mode->getHeight(); ++y) {
61  for (int x = 0; x < m_mode->getWidth(); ++x) {
62  processCell(*image, x, y, *m_mode, *m_sigma);
63  }
64  }
65  }
66 }
67 
68 template<typename T>
70  return m_mode;
71 }
72 
73 template<typename T>
75  return m_sigma;
76 }
77 
78 template<typename T>
80  return m_var_mode;
81 }
82 
83 template<typename T>
85  return m_var_sigma;
86 }
87 
88 template<typename T>
90  Histogram<double, int64_t> histo(data.begin(), data.end(), KappaSigmaBinning<double>(m_kappa1, m_kappa2));
91 
92  auto ref_bin = histo.getBinEdges(0);
93  auto atol = (ref_bin.second - ref_bin.first) * 0.1;
94 
95  T mean, median, sigma;
96  std::tie(mean, median, sigma) = histo.getStats();
97  T prev_sigma = sigma * 10;
98 
99  assert(!std::isnan(mean));
100 
101  for (size_t iter = 0; iter < m_max_iter && sigma > atol && std::abs(sigma / prev_sigma - 1.0) > m_rtol; ++iter) {
102  histo.clip(median - sigma * m_kappa3, median + sigma * m_kappa3);
103  prev_sigma = sigma;
104  std::tie(mean, median, sigma) = histo.getStats();
105  }
106 
107  // Sigma is 0
108  T mode;
109  if (std::abs(sigma) == 0) {
110  mode = mean;
111  }
112  // Not crowded: mean and median do not differ more than 30%
113  else if (std::abs((mean - median) / sigma) < 0.3) {
114  mode = 2.5 * median - 1.5 * mean;
115  }
116  // Crowded case: we use the median
117  else {
118  mode = median;
119  }
120  return std::make_tuple(mode, sigma);
121 }
122 
123 template<typename T>
124 void ImageMode<T>::processCell(const Image<T>& img, int x, int y,
125  VectorImage<T>& out_mode, VectorImage<T>& out_sigma) const {
126  int off_x = x * m_cell_w;
127  int off_y = y * m_cell_h;
128  int w = std::min(m_cell_w, img.getWidth() - off_x);
129  int h = std::min(m_cell_h, img.getHeight() - off_y);
130 
131  auto img_chunk_ptr = img.getChunk(off_x, off_y, w, h);
132  auto& img_chunk = *img_chunk_ptr;
133 
134  std::vector<T> filtered;
135  filtered.reserve(w * h);
136 
137  for (int y = 0; y < h; ++y) {
138  for (int x = 0; x < w; ++x) {
139  auto v = img_chunk.getValue(x, y);
140  if (v != m_invalid)
141  filtered.emplace_back(v);
142  }
143  }
144 
145  if (filtered.size() / static_cast<float>(w * h) < 0.5) {
146  out_mode.setValue(x, y, m_invalid);
147  out_sigma.setValue(x, y, m_invalid);
148  }
149  else {
150  T mode, sigma;
151  std::tie(mode, sigma) = getBackGuess(filtered);
152  out_mode.setValue(x, y, mode);
153  out_sigma.setValue(x, y, sigma);
154  }
155 }
156 
157 template
158 class ImageMode<SeFloat>;
159 
160 } // end of namespace SourceXtractor
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
T atol(T... args)
T begin(T... args)
std::tuple< VarType, VarType, VarType > getStats() const
void clip(VarType min, VarType max)
std::pair< VarType, VarType > getBinEdges(size_t i) const
std::shared_ptr< VectorImage< T > > getVarianceSigmaImage() const
Definition: ImageMode.cpp:84
std::shared_ptr< VectorImage< T > > m_var_mode
Definition: ImageMode.h:107
std::shared_ptr< VectorImage< T > > getVarianceModeImage() const
Definition: ImageMode.cpp:79
void processCell(const Image< T > &img, int x, int y, VectorImage< T > &out_mode, VectorImage< T > &out_sigma) const
Definition: ImageMode.cpp:124
std::shared_ptr< VectorImage< T > > getSigmaImage() const
Definition: ImageMode.cpp:74
std::shared_ptr< VectorImage< T > > m_mode
Definition: ImageMode.h:106
std::shared_ptr< VectorImage< T > > getModeImage() const
Definition: ImageMode.cpp:69
std::shared_ptr< VectorImage< T > > m_var_sigma
Definition: ImageMode.h:107
ImageMode(const std::shared_ptr< Image< T >> &image, const std::shared_ptr< Image< T >> &variance, int cell_w, int cell_h, T invalid_value, T kappa1=2, T kappa2=5, T kappa3=3, T rtol=1e-4, size_t max_iter=100)
Definition: ImageMode.cpp:30
std::shared_ptr< VectorImage< T > > m_sigma
Definition: ImageMode.h:106
std::tuple< T, T > getBackGuess(const std::vector< T > &data) const
Definition: ImageMode.cpp:89
Interface representing an image.
Definition: Image.h:43
virtual std::shared_ptr< ImageChunk< T > > getChunk(int x, int y, int width, int height) const =0
virtual int getHeight() const =0
Returns the height of the image in pixels.
virtual int getWidth() const =0
Returns the width of the image in pixels.
Image implementation which keeps the pixel values in memory.
Definition: VectorImage.h:52
void setValue(int x, int y, T value) final
Definition: VectorImage.h:124
static std::shared_ptr< VectorImage< T > > create(Args &&... args)
Definition: VectorImage.h:100
T div(T... args)
T emplace_back(T... args)
T end(T... args)
T isnan(T... args)
T make_tuple(T... args)
T min(T... args)
T reserve(T... args)
T size(T... args)
T tie(T... args)