diff options
Diffstat (limited to 'svg-flatten/src/nopencv.cpp')
-rw-r--r-- | svg-flatten/src/nopencv.cpp | 222 |
1 files changed, 203 insertions, 19 deletions
diff --git a/svg-flatten/src/nopencv.cpp b/svg-flatten/src/nopencv.cpp index 22c3fff..f42b6ad 100644 --- a/svg-flatten/src/nopencv.cpp +++ b/svg-flatten/src/nopencv.cpp @@ -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); + }; +} + |