kopia lustrzana https://github.com/jaseg/gerbolyze
Work on chain approx
jaseg
rodzic
bbf1c02e79
commit
3386e586ac
|
|
@ -35,6 +35,9 @@ namespace gerbolyze {
|
|||
typedef std::array<double, 2> d2p;
|
||||
typedef std::vector<d2p> Polygon;
|
||||
|
||||
typedef std::array<int64_t, 2> i2p;
|
||||
typedef std::vector<i2p> Polygon_i;
|
||||
|
||||
class xform2d {
|
||||
public:
|
||||
xform2d(double xx, double yx, double xy, double yy, double x0=0.0, double y0=0.0) :
|
||||
|
|
|
|||
|
|
@ -1,5 +1,6 @@
|
|||
|
||||
#include <iostream>
|
||||
#include <iomanip>
|
||||
|
||||
#include "nopencv.hpp"
|
||||
|
||||
|
|
@ -37,17 +38,17 @@ static struct {
|
|||
} dir_to_coords[8] = {{0, -1}, {1, -1}, {1, 0}, {1, 1}, {0, 1}, {-1, 1}, {-1, 0}, {-1, -1}};
|
||||
|
||||
static Direction flip_direction[8] = {
|
||||
D_S,
|
||||
D_SW,
|
||||
D_W,
|
||||
D_NW,
|
||||
D_N,
|
||||
D_NE,
|
||||
D_E,
|
||||
D_SE
|
||||
D_S, /* 0 */
|
||||
D_SW, /* 1 */
|
||||
D_W, /* 2 */
|
||||
D_NW, /* 3 */
|
||||
D_N, /* 4 */
|
||||
D_NE, /* 5 */
|
||||
D_E, /* 6 */
|
||||
D_SE /* 7 */
|
||||
};
|
||||
|
||||
static void follow(gerbolyze::nopencv::Image32 &img, int start_x, int start_y, Direction initial_direction, int nbd, int connectivity, Polygon &poly) {
|
||||
static void follow(gerbolyze::nopencv::Image32 &img, int start_x, int start_y, Direction initial_direction, int nbd, int connectivity, Polygon_i &poly) {
|
||||
|
||||
//cerr << "follow " << start_x << " " << start_y << " | dir=" << dir_str[initial_direction] << " nbd=" << nbd << " conn=" << connectivity << endl;
|
||||
int dir_inc = (connectivity == 4) ? 2 : 1;
|
||||
|
|
@ -69,10 +70,11 @@ static void follow(gerbolyze::nopencv::Image32 &img, int start_x, int start_y, D
|
|||
|
||||
if (!found) { /* No nonzero pixels found. This is a single-pixel contour */
|
||||
img.at(start_x, start_y) = nbd;
|
||||
poly.emplace_back(d2p{(double)start_x, (double)start_y});
|
||||
poly.emplace_back(d2p{(double)start_x+1, (double)start_y});
|
||||
poly.emplace_back(d2p{(double)start_x+1, (double)start_y+1});
|
||||
poly.emplace_back(d2p{(double)start_x, (double)start_y+1});
|
||||
poly.emplace_back(i2p{start_x, start_y});
|
||||
poly.emplace_back(i2p{start_x+1, start_y});
|
||||
poly.emplace_back(i2p{start_x+1, start_y+1});
|
||||
poly.emplace_back(i2p{start_x, start_y+1});
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
|
|
@ -106,10 +108,10 @@ static void follow(gerbolyze::nopencv::Image32 &img, int start_x, int start_y, D
|
|||
|
||||
for (int l = (current_direction + 8 - 2 + 1) / 2 * 2; l > k; l -= dir_inc) {
|
||||
switch (l%8) {
|
||||
case 0: poly.emplace_back(d2p{(double)center_x, (double)center_y}); break;
|
||||
case 2: poly.emplace_back(d2p{(double)center_x+1, (double)center_y}); break;
|
||||
case 4: poly.emplace_back(d2p{(double)center_x+1, (double)center_y+1}); break;
|
||||
case 6: poly.emplace_back(d2p{(double)center_x, (double)center_y+1}); break;
|
||||
case 0: poly.emplace_back(i2p{center_x, center_y}); break;
|
||||
case 2: poly.emplace_back(i2p{center_x+1, center_y}); break;
|
||||
case 4: poly.emplace_back(i2p{center_x+1, center_y+1}); break;
|
||||
case 6: poly.emplace_back(i2p{center_x, center_y+1}); break;
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -123,8 +125,15 @@ static void follow(gerbolyze::nopencv::Image32 &img, int start_x, int start_y, D
|
|||
|
||||
|
||||
void gerbolyze::nopencv::find_blobs(gerbolyze::nopencv::Image32 &img, gerbolyze::nopencv::ContourCallback cb) {
|
||||
/* Implementation of the hierarchical contour finding algorithm from Suzuki and Abe, 1983: Topological Structural
|
||||
* Analysis of Digitized Binary Images by Border Following
|
||||
*
|
||||
* Written with these two resources as reference:
|
||||
* https://theailearner.com/tag/suzuki-contour-algorithm-opencv/
|
||||
* https://github.com/FreshJesh5/Suzuki-Algorithm/blob/master/contoursv1/contoursv1.cpp
|
||||
*/
|
||||
int nbd = 1;
|
||||
Polygon poly;
|
||||
Polygon_i poly;
|
||||
for (int y=0; y<img.rows(); y++) {
|
||||
for (int x=0; x<img.cols(); x++) {
|
||||
int val_xy = img.at(x, y);
|
||||
|
|
@ -147,7 +156,7 @@ void gerbolyze::nopencv::find_blobs(gerbolyze::nopencv::Image32 &img, gerbolyze:
|
|||
|
||||
} else if (val_xy >= 1 && img.at_default(x+1, y) == 0) { /* hole border starting point */
|
||||
nbd += 1;
|
||||
follow(img, x, y, D_E, nbd, 4, poly);
|
||||
follow(img, x, y, D_E, nbd, 8, poly);
|
||||
cb(poly, CP_HOLE);
|
||||
poly.clear();
|
||||
}
|
||||
|
|
@ -155,3 +164,178 @@ void gerbolyze::nopencv::find_blobs(gerbolyze::nopencv::Image32 &img, gerbolyze:
|
|||
}
|
||||
}
|
||||
|
||||
static size_t region_of_support(Polygon_i poly, size_t i) {
|
||||
double x0 = poly[i][0], y0 = poly[i][1];
|
||||
size_t sz = poly.size();
|
||||
double last_l = 0;
|
||||
double last_r = 0;
|
||||
size_t k;
|
||||
for (k=0; k<sz/2; k++) {
|
||||
size_t idx1 = (i + k) % sz;
|
||||
size_t idx2 = (i + sz - k) % sz;
|
||||
double x1 = poly[idx1][0], y1 = poly[idx1][1], x2 = poly[idx2][0], y2 = poly[idx2][1];
|
||||
double l = sqrt(pow(x2-x1, 2) + pow(y2-y1, 2));
|
||||
/* https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line
|
||||
* TODO: Check whether distance-to-line is an ok implementation here, the paper asks for distance to chord.
|
||||
*/
|
||||
double d = ((x2-x1)*(y1-y0) - (x1-x0)*(y2-y1)) / sqrt(pow(x2-x1, 2) + pow(y2-y1, 2));
|
||||
double r = d/l;
|
||||
|
||||
bool cond_a = l < last_l;
|
||||
bool cond_b = (d > 0) ? (r < last_r) : (r > last_r);
|
||||
|
||||
if (cond_a || cond_b)
|
||||
break;
|
||||
|
||||
last_l = l;
|
||||
last_r = r;
|
||||
}
|
||||
k -= 1;
|
||||
return k;
|
||||
}
|
||||
|
||||
int freeman_angle(const Polygon_i &poly, size_t i) {
|
||||
/* f:
|
||||
* 2
|
||||
* 3 1
|
||||
* ^
|
||||
* |
|
||||
* 4 <--- X ---> 0
|
||||
* |
|
||||
* v
|
||||
* 5 7
|
||||
* 6
|
||||
*
|
||||
*/
|
||||
size_t sz = poly.size();
|
||||
|
||||
auto &p_last = poly[(i + sz - 1) % sz];
|
||||
auto &p_now = poly[i];
|
||||
auto dx = p_now[0] - p_last[0];
|
||||
auto dy = p_now[1] - p_last[1];
|
||||
/* both points must be neighbors */
|
||||
assert (-1 <= dx && dx <= 1);
|
||||
assert (-1 <= dy && dy <= 1);
|
||||
assert (!(dx == 0 && dy == 0));
|
||||
|
||||
int lut[3][3] = {{3, 2, 1}, {4, -1, 0}, {5, 6, 7}};
|
||||
return lut[dy+1][dx+1];
|
||||
}
|
||||
|
||||
double k_curvature(const Polygon_i &poly, size_t i, size_t k) {
|
||||
size_t sz = poly.size();
|
||||
double acc = 0;
|
||||
for (size_t idx = 0; idx < k; idx++) {
|
||||
acc += freeman_angle(poly, (i + 2*sz - idx) % sz) - freeman_angle(poly, (i+idx + 1) % sz);
|
||||
}
|
||||
return acc / (2*k);
|
||||
}
|
||||
|
||||
double k_cos(const Polygon_i &poly, size_t i, size_t k) {
|
||||
size_t sz = poly.size();
|
||||
int64_t x0 = poly[i][0], y0 = poly[i][1];
|
||||
int64_t x1 = poly[(i + sz + k) % sz][0], y1 = poly[(i + sz + k) % sz][1];
|
||||
int64_t x2 = poly[(i + sz - k) % sz][0], y2 = poly[(i + sz - k) % sz][1];
|
||||
auto xa = x0 - x1, ya = y0 - y1;
|
||||
auto xb = x0 - x2, yb = y0 - y2;
|
||||
auto dp = xa*yb + ya*xb;
|
||||
auto sq_a = xa*xa + ya*ya;
|
||||
auto sq_b = xb*xb + yb*yb;
|
||||
return dp / (sqrt(sq_a)*sqrt(sq_b));
|
||||
}
|
||||
|
||||
ContourCallback gerbolyze::nopencv::simplify_contours_teh_chin(int kcos, ContourCallback cb) {
|
||||
return [&cb, kcos](Polygon_i &poly, ContourPolarity cpol) {
|
||||
size_t sz = poly.size();
|
||||
vector<size_t> ros(sz);
|
||||
vector<double> sig(sz);
|
||||
vector<bool> retain(sz);
|
||||
for (size_t i=0; i<sz; i++) {
|
||||
ros[i] = region_of_support(poly, i);
|
||||
sig[i] = fabs(k_cos(poly, i, kcos));
|
||||
retain[i] = true;
|
||||
}
|
||||
|
||||
cerr << endl;
|
||||
cerr << "Polarity: " << cpol <<endl;
|
||||
cerr << "Coords:"<<endl;
|
||||
cerr << " x: ";
|
||||
for (size_t i=0; i<sz; i++) {
|
||||
cerr << setfill(' ') << setw(2) << poly[i][0] << " ";
|
||||
}
|
||||
cerr << endl;
|
||||
cerr << " y: ";
|
||||
for (size_t i=0; i<sz; i++) {
|
||||
cerr << setfill(' ') << setw(2) << poly[i][1] << " ";
|
||||
}
|
||||
cerr << endl;
|
||||
cerr << "Metrics:"<<endl;
|
||||
cerr << "ros: ";
|
||||
for (size_t i=0; i<sz; i++) {
|
||||
cerr << setfill(' ') << setw(2) << ros[i] << " ";
|
||||
}
|
||||
cerr << endl;
|
||||
cerr << "sig: ";
|
||||
for (size_t i=0; i<sz; i++) {
|
||||
cerr << setfill(' ') << setw(2) << sig[i] << " ";
|
||||
}
|
||||
cerr << endl;
|
||||
|
||||
/* 3a, Pass 1: Non-maxima suppression */
|
||||
for (size_t i=0; i<sz; i++) {
|
||||
for (size_t j=0; j<ros[i]/2; j++) {
|
||||
if (sig[i] < sig[(i + j) % sz] || sig[i] < sig[(i + sz - j) % sz]) {
|
||||
retain[i] = false;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
cerr << "pass 1: ";
|
||||
for (size_t i=0; i<sz; i++) {
|
||||
cerr << (retain[i] ? "#" : ".");
|
||||
}
|
||||
cerr << endl;
|
||||
|
||||
/* 3b, Pass 2: Zero-curvature suppression */
|
||||
for (size_t i=0; i<sz; i++) {
|
||||
if (retain[i] && k_cos(poly, i, kcos) == 0) {
|
||||
retain[i] = false;
|
||||
}
|
||||
}
|
||||
|
||||
cerr << "pass 2: ";
|
||||
for (size_t i=0; i<sz; i++) {
|
||||
cerr << (retain[i] ? "#" : ".");
|
||||
}
|
||||
cerr << endl;
|
||||
|
||||
/* 3c, Pass 3: Further thinning */
|
||||
for (size_t i=0; i<sz; i++) {
|
||||
if (retain[i]) {
|
||||
if (ros[i] == 1) {
|
||||
if (retain[(i + sz - 1) % sz] || retain[(i + 1)%sz]) {
|
||||
if (sig[i] < sig[(i + sz - 1)%sz] || sig[i] < sig[(i + 1)%sz]) {
|
||||
retain[i] = false;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
cerr << "pass 3: ";
|
||||
for (size_t i=0; i<sz; i++) {
|
||||
cerr << (retain[i] ? "#" : ".");
|
||||
}
|
||||
cerr << endl;
|
||||
|
||||
Polygon_i new_poly;
|
||||
for (size_t i=0; i<sz; i++) {
|
||||
if (retain[i]) {
|
||||
new_poly.push_back(poly[i]);
|
||||
}
|
||||
}
|
||||
cb(new_poly, cpol);
|
||||
};
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -37,7 +37,7 @@ namespace gerbolyze {
|
|||
CP_HOLE
|
||||
};
|
||||
|
||||
typedef std::function<void(Polygon, ContourPolarity)> ContourCallback;
|
||||
typedef std::function<void(Polygon_i&, ContourPolarity)> ContourCallback;
|
||||
|
||||
class Image32 {
|
||||
public:
|
||||
|
|
@ -101,6 +101,7 @@ namespace gerbolyze {
|
|||
};
|
||||
|
||||
void find_blobs(Image32 &img, ContourCallback cb);
|
||||
ContourCallback simplify_contours_teh_chin(int kcos, ContourCallback cb);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -40,7 +40,7 @@ MU_TEST(test_complex_example_from_paper) {
|
|||
};
|
||||
Image32 test_img(9, 6, static_cast<int*>(img_data));
|
||||
|
||||
const Polygon expected_polys[3] = {
|
||||
const Polygon_i expected_polys[3] = {
|
||||
{
|
||||
{1,1}, {1,2}, {1,3}, {1,4}, {1,5},
|
||||
{2,5}, {3,5}, {4,5}, {5,5}, {6,5}, {7,5}, {8,5},
|
||||
|
|
@ -64,18 +64,18 @@ MU_TEST(test_complex_example_from_paper) {
|
|||
const ContourPolarity expected_polarities[3] = {CP_CONTOUR, CP_HOLE, CP_HOLE};
|
||||
|
||||
int invocation_count = 0;
|
||||
gerbolyze::nopencv::find_blobs(test_img, [&invocation_count, &expected_polarities, &expected_polys](Polygon poly, ContourPolarity pol) {
|
||||
gerbolyze::nopencv::find_blobs(test_img, [&invocation_count, &expected_polarities, &expected_polys](Polygon_i &poly, ContourPolarity pol) {
|
||||
invocation_count += 1;
|
||||
mu_assert((invocation_count <= 3), "Too many contours returned");
|
||||
|
||||
mu_assert(poly.size() > 0, "Empty contour returned");
|
||||
mu_assert_int_eq(pol, expected_polarities[invocation_count-1]);
|
||||
|
||||
d2p last;
|
||||
i2p last;
|
||||
bool first = true;
|
||||
Polygon exp = expected_polys[invocation_count-1];
|
||||
Polygon_i exp = expected_polys[invocation_count-1];
|
||||
//cout << "poly: ";
|
||||
for (d2p &p : poly) {
|
||||
for (auto &p : poly) {
|
||||
//cout << "(" << p[0] << ", " << p[1] << "), ";
|
||||
if (!first) {
|
||||
mu_assert((fabs(p[0] - last[0]) + fabs(p[1] - last[1]) == 1), "Subsequent contour points have distance other than one");
|
||||
|
|
@ -149,7 +149,7 @@ static void testdata_roundtrip(const char *fn) {
|
|||
<< "xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">" << endl;
|
||||
svg << "<rect width=\"100%\" height=\"100%\" fill=\"black\"/>" << endl;
|
||||
|
||||
gerbolyze::nopencv::find_blobs(ref_img, [&svg](Polygon poly, ContourPolarity pol) {
|
||||
gerbolyze::nopencv::find_blobs(ref_img, [&svg](Polygon_i &poly, ContourPolarity pol) {
|
||||
mu_assert(poly.size() > 0, "Empty contour returned");
|
||||
mu_assert(poly.size() > 2, "Contour has less than three points, no area");
|
||||
mu_assert(pol == CP_CONTOUR || pol == CP_HOLE, "Contour has invalid polarity");
|
||||
|
|
@ -212,7 +212,45 @@ MU_TEST(test_round_trip_two_px) { testdata_roundtrip("testdata/two-p
|
|||
MU_TEST(test_round_trip_two_px_inv) { testdata_roundtrip("testdata/two-px-inv.png"); }
|
||||
|
||||
|
||||
static void chain_approx_test(const char *fn) {
|
||||
int x, y;
|
||||
uint8_t *data = stbi_load(fn, &x, &y, nullptr, 1);
|
||||
Image32 ref_img(x, y);
|
||||
for (int cy=0; cy<y; cy++) {
|
||||
for (int cx=0; cx<x; cx++) {
|
||||
ref_img.at(cx, cy) = data[cy*x + cx] / 255;
|
||||
}
|
||||
}
|
||||
stbi_image_free(data);
|
||||
|
||||
TempfileHack tmp_svg(".svg");
|
||||
|
||||
//ofstream svg(tmp_svg.c_str());
|
||||
ofstream svg("/tmp/teh-chin-test.svg");
|
||||
|
||||
svg << "<svg width=\"" << x << "px\" height=\"" << y << "px\" viewBox=\"0 0 "
|
||||
<< x << " " << y << "\" "
|
||||
<< "xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">" << endl;
|
||||
svg << "<rect width=\"100%\" height=\"100%\" fill=\"black\"/>" << endl;
|
||||
|
||||
gerbolyze::nopencv::find_blobs(ref_img, simplify_contours_teh_chin(2, [&svg](Polygon_i &poly, ContourPolarity pol) {
|
||||
svg << "<path fill=\"" << ((pol == CP_HOLE) ? "black" : "white") << "\" d=\"";
|
||||
svg << "M " << poly[0][0] << " " << poly[0][1];
|
||||
for (size_t i=1; i<poly.size(); i++) {
|
||||
svg << " L " << poly[i][0] << " " << poly[i][1];
|
||||
}
|
||||
svg << " Z\"/>" << endl;
|
||||
}));
|
||||
svg << "</svg>" << endl;
|
||||
svg.close();
|
||||
}
|
||||
|
||||
|
||||
MU_TEST(chain_approx_test_chromosome) { chain_approx_test("testdata/chain-approx-teh-chin-chromosome.png"); }
|
||||
|
||||
|
||||
MU_TEST_SUITE(nopencv_contours_suite) {
|
||||
/*
|
||||
MU_RUN_TEST(test_complex_example_from_paper);
|
||||
MU_RUN_TEST(test_round_trip_blank);
|
||||
MU_RUN_TEST(test_round_trip_white);
|
||||
|
|
@ -229,6 +267,8 @@ MU_TEST_SUITE(nopencv_contours_suite) {
|
|||
MU_RUN_TEST(test_round_trip_two_blobs);
|
||||
MU_RUN_TEST(test_round_trip_two_px);
|
||||
MU_RUN_TEST(test_round_trip_two_px_inv);
|
||||
*/
|
||||
MU_RUN_TEST(chain_approx_test_chromosome);
|
||||
};
|
||||
|
||||
int main(int argc, char **argv) {
|
||||
|
|
|
|||
Plik binarny nie jest wyświetlany.
|
Po Szerokość: | Wysokość: | Rozmiar: 8.1 KiB |
Ładowanie…
Reference in New Issue