MedianFilter.h

Go to the documentation of this file.

 
#ifndef SOURCEXTRACTORPLUSPLUS_MEDIANFILTER_H
#define SOURCEXTRACTORPLUSPLUS_MEDIANFILTER_H
 
#include "SEFramework/Image/VectorImage.h"
#include "SEFramework/Image/WriteableImage.h"
#include <array>
 
namespace SourceXtractor {

template<typename T>
class MedianFilter {
public:
  MedianFilter(int box_width, int box_height) : m_box_width(box_width), m_box_height(box_height) {
  }

  explicit MedianFilter(const std::array<int, 2>& box) : MedianFilter(box[0], box[1]) {
  }

  auto operator()(const VectorImage<T>& image, const VectorImage<T>& variance,
                  T threshold = 0) const -> std::pair<std::shared_ptr<VectorImage<T>>, std::shared_ptr<VectorImage<T>>> {
    assert(image.getWidth() == variance.getWidth());
    assert(image.getHeight() == variance.getHeight());
 
    auto out_img = VectorImage<T>::create(image.getWidth(), image.getHeight());
    auto out_var = VectorImage<T>::create(image.getWidth(), image.getHeight());
 
    for (int y = 0; y < image.getHeight(); ++y) {
      for (int x = 0; x < image.getWidth(); ++x) {
        auto box = getBox(image, x, y);
        auto median = getMedian(box);
        auto value = image.getValue(x, y);
        if (std::abs(median - value) >= threshold) {
          out_img->setValue(x, y, median);
          auto var_box = getBox(variance, x, y);
          out_var->setValue(x, y, getMedian(var_box));
        }
        else {
          out_img->setValue(x, y, value);
          out_var->setValue(x, y, variance.getValue(x, y));
        }
      }
    }
 
    return std::make_pair(out_img, out_var);
  }
 
private:
  int m_box_width, m_box_height;

  static T getMedian(std::vector<T>& data) {
    std::sort(data.begin(), data.end());
    auto nitems = data.size();
    if (nitems % 2 == 1)
      return data[nitems / 2];
    return (data[nitems / 2] + data[nitems / 2 - 1]) / 2;
  }

  static int clip(int position, int box_size, int image_size) {
    box_size /= 2;
    if (box_size > position)
      return position;
    if (box_size > image_size - position - 1)
      return image_size - position - 1;
    return box_size;
  }

  std::vector<T> getBox(const VectorImage<T>& img, int x, int y) const {
    int hw = clip(x, m_box_width, img.getWidth());
    int hh = clip(y, m_box_height, img.getHeight());
    std::vector<T> data;
    auto inserter = std::back_inserter(data);
    for (int iy = -hh; iy < hh + 1; ++iy) {
      for (int ix = -hw; ix < hw + 1; ++ix) {
        *inserter = img.getValue(x + ix, y + iy);
        ++inserter;
      }
    }
    return data;
  }
};
 
} // end of namespace SourceXtractor
 
#endif //SOURCEXTRACTORPLUSPLUS_MEDIANFILTER_H