estimation.cpp

Go to the documentation of this file.

#include <CGAL/IO/Qt_widget.h>
#include "estimation.hpp"
#include "jet.hpp"

using namespace Arak;

PointSetColorEstimator::PointSetColorEstimator(ArakMarkovChain& chain,
                 const Arak::Util::PropertyMap& props) : chain(chain) {
  using namespace Arak::Util;
  int n = 1;
  if (hasp(props, "arak.mcmc.estimation.grid_size")) {
    parse(getp(props, "arak.mcmc.estimation.grid_size"), n);
    assert(n > 0);
  }
  Geometry::Rectangle r = chain.getState().boundary();
  Geometry::Kernel::FT width = r.xmax() - r.xmin();
  Geometry::Kernel::FT height = r.ymax() - r.ymin();
  double aspect_ratio = CGAL::to_double(width) / CGAL::to_double(height);
  int rows = std::max(1, int(sqrt(double(n) / aspect_ratio)));
  int cols = n / rows;
  index = new QueryPointIndex(r, rows, cols);
  index->addListeners(*this);
  chain.getState().addQueryPoints(*index);
}

void PointSetColorEstimator::visualize(CGAL::Qt_widget& widget) const {
  using namespace CGAL;
  widget << PointSize(6) << PointStyle(BOX);
  typedef std::list<PointColorEstimator*>::const_iterator Iterator;
  for (Iterator it = estimators.begin(); it != estimators.end(); it++) {
    const PointColorEstimator* estimator = *it;
    // Gamma-correct the intensities for an accurate visualization.
    double intensity = 1.0 - estimator->estimate();
    const double gamma = 1.0 / 0.45;
    double gc_intensity = pow(intensity, gamma);
    int gray = int(round(255.0 * gc_intensity));
    CGAL::Color c(gray, gray, gray);
    widget << c << estimator->point;
  }
}

GridColorEstimator::GridColorEstimator(ArakMarkovChain& chain,
               const Arak::Util::PropertyMap& props)
  : chain(chain) {
  using namespace Arak::Util;
  int n = 1;
  if (hasp(props, "arak.mcmc.estimation.grid_size")) {
    parse(getp(props, "arak.mcmc.estimation.grid_size"), n);
    assert(n > 0);
  }
  Geometry::Rectangle r = chain.getState().boundary();
  Geometry::Kernel::FT width = r.xmax() - r.xmin();
  Geometry::Kernel::FT height = r.ymax() - r.ymin();
  double aspect_ratio = CGAL::to_double(width) / CGAL::to_double(height);
  int rows = std::max(1, int(sqrt(double(n) / aspect_ratio)));
  int cols = n / rows;
  grid = new EstimatorGrid(r, rows, cols);
  for (int i = 0; i < rows; i++)
    for (int j = 0; j < cols; j++) {
      EstimatorGrid::Cell& cell = grid->getCell(i, j);
      Geometry::Point p = CGAL::midpoint(cell[0], cell[2]);
      Color color = chain.getState().color(p);
      PointColorEstimator* pce = new PointColorEstimator(chain, p, color);
      estimators.push_back(pce);
      cell.add(pce);
    }
  chain.getState().addListener(*this);
}

void GridColorEstimator::recolored(const Geometry::Point& a,
           const Geometry::Point& b,
           const Geometry::Point& c) {
  Geometry::Triangle t(a, b, c);
  EstimatorGrid::TriangleCellIterator it(*grid, t), end;
  do {
    PointColorEstimator* pce = *(it->getItemList().begin());
    if (t.has_on_bounded_side(pce->point))
      pce->recolor(opposite(pce->color));
  } while (++it != end);
}

void GridColorEstimator::visualize(CGAL::Qt_widget& widget) const {
  using namespace CGAL;
  for (int i = grid->numRows() - 1; i >= 0; i--) {
    for (int j = 0; j < grid->numCols(); j++) {
      const EstimatorGrid::Cell& cell = grid->getCell(i, j);
      const PointColorEstimator* pce = *(cell.getItemList().begin());
      int n = int(double(Arak::jet_size - 1) * pce->estimate());    
      const CGAL::Color& c = Arak::jet_colors[n];
      // Cannot avoid plotting the cell boundary!
      widget.setFilled(true);
      widget.setLineWidth(0);
      widget.setPointSize(0);
      widget << c << CGAL::FillColor(c) << cell; 
    }
  }
}

---