aboutsummaryrefslogtreecommitdiff
path: root/svg-flatten
diff options
context:
space:
mode:
Diffstat (limited to 'svg-flatten')
-rw-r--r--svg-flatten/include/geom2d.hpp3
-rw-r--r--svg-flatten/src/nopencv.cpp222
-rw-r--r--svg-flatten/src/nopencv.hpp3
-rw-r--r--svg-flatten/src/nopencv_test.cpp52
-rw-r--r--svg-flatten/testdata/chain-approx-teh-chin-chromosome.pngbin0 -> 8310 bytes
5 files changed, 254 insertions, 26 deletions
diff --git a/svg-flatten/include/geom2d.hpp b/svg-flatten/include/geom2d.hpp
index fa5c545..b4cccf2 100644
--- a/svg-flatten/include/geom2d.hpp
+++ b/svg-flatten/include/geom2d.hpp
@@ -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) :
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);
+ };
+}
+
diff --git a/svg-flatten/src/nopencv.hpp b/svg-flatten/src/nopencv.hpp
index 98f6f50..48bb82b 100644
--- a/svg-flatten/src/nopencv.hpp
+++ b/svg-flatten/src/nopencv.hpp
@@ -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);
}
}
diff --git a/svg-flatten/src/nopencv_test.cpp b/svg-flatten/src/nopencv_test.cpp
index 21710b9..901c7a8 100644
--- a/svg-flatten/src/nopencv_test.cpp
+++ b/svg-flatten/src/nopencv_test.cpp
@@ -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) {
diff --git a/svg-flatten/testdata/chain-approx-teh-chin-chromosome.png b/svg-flatten/testdata/chain-approx-teh-chin-chromosome.png
new file mode 100644
index 0000000..86ae8a2
--- /dev/null
+++ b/svg-flatten/testdata/chain-approx-teh-chin-chromosome.png
Binary files differ
generated by cgit (git 2.34.1) at 2024-11-05 03:01:29 +0000