getquad.cpp 2.97 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
#include "getquad.hpp"
#include <mln/accu/accumulators/covariance.hpp>
#include <mln/colors/literal.hpp>
#include <mln/core/algorithm/fill.hpp>
#include <mln/io/imsave.hpp>
#include <mln/accu/accumulators/infsup.hpp>
#include <mln/draw/range/line2d.hpp>
#include <Eigen/Dense>

pquad_t
getquad(const mln::image2d<bool>& bin,
        const char* debug_filename)
{
  using namespace mln;
  Eigen::Array<double, 2, 2> COV;
  Eigen::Array<double, 2, 2> P, Pinv;

  // 1. Computee covariance matrix.
  {
    accu::accumulators::covariance<point2d, double, double> accu;

    mln_foreach(auto px, bin.pixels())
      if (px.val())
        accu.take(px.point());

    COV = accu.to_result();
  }



  // 2. Compute the eigen vectors
  {
    Eigen::SelfAdjointEigenSolver<Eigen::Matrix2d> eigensolver(COV.matrix());
    if (eigensolver.info() != Eigen::Success)
      std::abort();
    P = eigensolver.eigenvectors();
    Pinv = P.matrix().inverse();

    // Debug
    std::cout << "COV: \n" << COV << std::endl;
    std::cout << "Lambdas: \n" <<  eigensolver.eigenvalues() << std::endl;
    std::cout << "P: \n" << P << std::endl;
    std::cout << "Pinv: \n" << Pinv << std::endl;
  }


  // 3. Compute the bbox
  point2df pmin, pmax;
  {
    accu::accumulators::infsup<point2df> accu;

    mln_foreach(auto px, bin.pixels())
      if (px.val())
        {
          point2df p;
          p[0] = P(0,0) * px.point()[0] + P(1,0) * px.point()[1];
          p[1] = P(0,1) * px.point()[0] + P(1,1) * px.point()[1];
          accu.take(p);
        }

    pmin = accu::extractor::inf(accu);
    pmax = accu::extractor::sup(accu);
  }

  pquad_t res;
  res.tl[0] = Pinv(0,0) * pmin[0] + Pinv(1,0) * pmin[1];
  res.tl[1] = Pinv(0,1) * pmin[0] + Pinv(1,1) * pmin[1];

  res.tr[0] = Pinv(0,0) * pmin[0] + Pinv(1,0) * pmax[1];
  res.tr[1] = Pinv(0,1) * pmin[0] + Pinv(1,1) * pmax[1];

  res.br[0] = Pinv(0,0) * pmax[0] + Pinv(1,0) * pmax[1];
  res.br[1] = Pinv(0,1) * pmax[0] + Pinv(1,1) * pmax[1];

  res.bl[0] = Pinv(0,0) * pmax[0] + Pinv(1,0) * pmin[1];
  res.bl[1] = Pinv(0,1) * pmax[0] + Pinv(1,1) * pmin[1];

  std::cout << res.tl << ","
            << res.tr << ","
            << res.br << ","
            << res.bl << std::endl;


  if (debug_filename != NULL)
    {
      accu::accumulators::infsup<point2df> accu;
      accu.take(bin.domain().pmin);
      accu.take(bin.domain().pmax);
      accu.take(res.tl);
      accu.take(res.tr);
      accu.take(res.bl);
      accu.take(res.br);

      box2d newdomain = { (point2d) accu::extractor::inf(accu),
                          (point2d) accu::extractor::sup(accu) };


      image2d<rgb8> out(newdomain);
      fill(out | where(bin), colors::literal::white);
      fill(out | draw::line2d(res.tl, res.tr), colors::literal::red);
      fill(out | draw::line2d(res.tr, res.br), colors::literal::red);
      fill(out | draw::line2d(res.br, res.bl), colors::literal::red);
      fill(out | draw::line2d(res.bl, res.tl), colors::literal::red);

      io::imsave(out, debug_filename);
    }


  return res;
}