Work on chain approx

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