From 3386e586ac0f0ae07535036efe6293f6118f7e2b Mon Sep 17 00:00:00 2001
From: jaseg <git@jaseg.de>
Date: Tue, 1 Jun 2021 23:36:32 +0200
Subject: [PATCH] Work on chain approx

---
 svg-flatten/include/geom2d.hpp                |   3 +
 svg-flatten/src/nopencv.cpp                   | 222 ++++++++++++++++--
 svg-flatten/src/nopencv.hpp                   |   3 +-
 svg-flatten/src/nopencv_test.cpp              |  52 +++-
 .../chain-approx-teh-chin-chromosome.png      | Bin 0 -> 8310 bytes
 5 files changed, 254 insertions(+), 26 deletions(-)
 create mode 100644 svg-flatten/testdata/chain-approx-teh-chin-chromosome.png

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 0000000000000000000000000000000000000000..86ae8a2bd2e7e39259898d41bdb585ca99367745
GIT binary patch
literal 8310
zcmeAS@N?(olHy`uVBq!ia0y~yU=U?sU=ZVAU|?Y2deQNOfk7~&DkP#LD6w3jpeR2r
zGbdG{q_QAYA+;hije()!*4kN_Ww%dFYWaU%tBTW<|Mbd_bM86n@87R6VUzCZ@3-pw
zw{UQ=WSvk{GD_V1|9yS`Px-n@pKl%s$;z{gu1o4ZTW%Nf`_Gr(?{fCMKQI6C>)hWT
z_Z0s8+WI`?fy6xf*RS_f&o7r}wrM}N{L5X3(tZ0R?CUN*%Kq~!bHQ~NN9|8F5$hXo
zWiNOzU;j>a!@r!<&%#cBkCymtx$j;5nFH%D-maKm`}Y6m`x_JXPMx7%vA>ynTDknA
z%?%6U)196>KeBD9Ry=<0epUIxePK<G`!-Gfn-}_QbGYB1TiWvT?JQNFJm2&1NdK`j
z^2hbhKmMoRo4>Ae|7^=M{WI$=>=%EYzy0$++w{=lb9H6claDR=xpn3Duy^WuGBa%M
z+%{T&DF0FM&h^JuOD!sY^eJ_@?4D)W=N?acyeN9z?#JJ3{mdhizT6h?yXP;nwAfVS
z{vV!e4`#2JZoRKXENn$Z-lKa5BL8m`D}SC_zU!ryu>SF>_Y7w|=F9)`?S94Iq6Nhw
zpU&O<wJCT?{a&7PMo%pF7K`!~K9_ZWDfQ#;uRqHB*=%O4o2f1Rrp(5<=Skegr50S*
zFE(-uA6Hk1ja<)YAXv<7d4q}n#v7-in>QcT7V=q#{kbf;Vv<SS1E%Sd;@4a{E?{(S
zWn6Kwf8qD2`G1O~&Kwfdd=gOL%Pkx4+R*z$bEA;6ZL3JIm*v!uU~k)Ync`EH({(<)
zR-F<NbaLsGsHm5wS0aMXZoQI~wR-PQF^1z;inP~myH=Hc>D-CAyRF5)*!|&9?lbwM
zvAECdmr70AoXlsDX{MXc?)$^VJtv}aYG`<L?QUO<w>vJM+x=$O?RV_~)(?+N_Op9>
z<+9)Z53`KT-`??GfB*Ze>gu!c{x1)GdT5%b`_#5Yc#g_Xp8lUNK2%I(=g82XfAdw+
z7q@kj&*um)Q$4RWdCulJHf`*W%j{&oKir)^@AtucwPXKJ<o|nXseJvfxlQ_$UGWv~
z()ZUs+^V{5rrC_wic1&vWHG-y?s_-;-$Bp)^(Om&B(+Sw*45tF&TYKrT-UX}<!dM3
z?dVP2)NAT>)A*NZ&zU&&l&i~LvxnF&dVQ_*jjot%t^Y0O^=5BYMV9@1Y4!PV>vNm?
zw#s)d7ra_7c;k<FxTK!+>b93m1!o)tUsq3WeiLzi-qR;_P8LVs)<j9Gp0-@hBkdHu
z`0chC$|cGF4{eB-{a?D^!j^*VU+ki19Pbod^?r6qH%md|M~;r{#k;=lOTN8y+o3v{
ziM_lvuam1!uC-p#=ls2xJ+ioFjjGWb4T<S*;`f>9b!?mex;a-e$LzhB<-Y5DXYSO#
zna)tRv+nc#jRMTF=6<SwV)Qa6|62ZSpO?((_o9!o?ruBtX<K&f_W*&}lRb1-Ot|=3
zg!9Q{E3>Y@TcrHz-dw5C{d=tI&B<lcq~&f-J>=DVC3)6FDVaKtdx7FBrw4yk*WdhD
z<ABc!%Y#QVkLo2T{p|R{9R4OGX_DL%-A2`|jGl8BaBEk;>7F6on&_)E_xOwjPA_fq
z9~obL+!$!yc(QF@ac|B3i`GwnvE3*Q7hq+1ohSI>X2cPh+sbxj<<={X?>=bo$+tDK
zV7ZWbZ@O&Wyr`Id=@6Fb=dym6R3Dm=mZh1zeyfv2+437XjZ*z*SswpjUjKLPZu{LQ
z{=8W%^lsZ5iKoe@jcv@AT%Yz?&SmZArQasqDpZvY`Cw&r*H5m6?fWe|z7M@&&c|K4
zLzH*7TL)*p`{uq)=fKmGO<Sh_oW_<b9NcKq$#CxM(vV~BdYV(cg2g81)-!CaZ@;O(
zf@8UmOk-B^{DWSfZmzf}b8bRr!(JhqhvuL7maduE@Pg+g14rgOZsm;gt~GZa3e3|;
zThvpa%)j*}_jXC)6kd*{w&_hwA4Nau&HB#1K|6>^?#c$Ct}}03?9V(Y^Zw+#%9_=?
zYW=%f+XcLnbasat+WReCVYyMO<YW0v%exQOwm*2gx<4jyzEQxT;`t5iZioI92z1T<
z<8##gL$BjxP4jNI1BOnQM435WvWp!&AA012?Yl!|$xnqD=gfW0S#Gvw!>pUew@u5~
zL{FUCQ@(B1DJ?0j&fB6=vpcLe*nD`sN<-t`Wruf8rB;{frHiurv-i~;{#&x;vXj7x
z9Is#Xg=}+D#TT7+xg_7=QnhkP2Sdco)l;6NxtwJ%Oze|>oN_|t(_OEvIXq3>3g#TY
zBWHg+DfjQB-?S^~A2X5}YxU(%ez11`+^v7?P5X+K={X6P8`nGZ-nTxox#$(^N|x^n
zv=|t|1O3@tpH;6Z4qM$+zEVLi?!itEcb4Njvplr>Dt(@2@mNPxY)r_=;1A+|YH&g7
zfbNY)9FHtJMVRuY7ey(uhN-y;o?iTpwQl|o-G+y?)~mw5JeFM|zbRdK^^sQg-`|<0
zyzUEN{4-g&@qu1i*b@IsCG0VYJ+IHd6XM%-Y~H#Tms8Y@x4*TraXRbKr_y!#UYD$h
zw%%#i>|<V2ZnJ!Mw7!?-wAOde(@BmT20k4pu2_Fsw?<Ub<$=bR8-ixd@d5K3{waiV
zPTd%|CDqa5@;TqAw!{esXYNWD{6A&QDJCPO#qN`R-kq{jiJUXt<Jqm4ea{?!s+CI}
z=&?<@a78US$i1*zS@Xogcl+|cZS?dw7V8~g&3)3GZO7_)s{L-u!Y}OFEW4rhaJE3$
zAsvn{JVyE=MVXx$Y&V+MD*89f(Q&`%&9D2?K}Gtf__4ref$HL6F8gkF_IGxyR`h%`
z^X}x(1^3Uc?+Q_sHQXptF8;Dd_hz7U=@AZ<>vaM~Rw^f=YRx_z_Y+Z4446L0u6pvj
zrT-4Tef&`1j-AvWzBUtq)f-gOtC(hZpZu_<+)()2lT}$=vYMqFRp!P=gH@Y?_ME=C
z<NM`(EXwEi{kYY=LvgWIZF%!s^UGf?XUSg@Vf_-P|H4v+>qIlJ^8$CZ)-Iz(D>5gU
z`AkpJ7n4ot&Tjb6SMkQ}Zs-c@9`pJ67gd~?o-^;^w)}K-UfepPHz(d*GV|g~IN{Rb
ze|EACi^_zqS4FF?t@!qX=c?B)FO!LG_gN&|KmYsC?;aNIAn<bY)}onL>H->s+#UC=
zxqsFmb-u&ux0{2foD0>CJ7t!(E#7ge+Kw*?M=i|Ph=fmUdU@<t=f3B&Cm(xXmB^jD
z<8}N8J5>qRn~}?Qm;1JvEO7}mQhYRTQ_hwiDTi+fQ7hEi+f1e%x%2Q^nq7a^(dEZy
zO2=OL%w(~o_lBVRLMP8NLMuANGIu@`b6+TD^k~tX3l0o{(ry#~I9AR4Q80^Hb;*ri
zPLVpM9;r|4EoNP3U-c|zas-c#xOAVz#TRR?nTzlynVKJ$;pn=0;qno_2SwL&-n*Su
zX*q3bKJlkh(x%C=&1+`9jeGne=S9>8A%@fo4DJt%4)<RQw4VJSU|DA1Kjkg+yu5=t
zloam!24@Sd{y(qc+0+aTMTY#}7HV<UJzRgpN;wY(HQ5Pz-CSJi!h5>>y}}`l%#bTe
zdM7pVp7*YMDzNn1$=`pv^e(+U*k+*7BHyP_pyyn`e=NseamQIz&I%R#qs?5s+)74g
zbxkjD%Nc*u@ZIbD=w8D-=Lx*mn9R=z?q+kUe$1qC!Tdw^!LBvh=T32)R9w7l)eC!V
z#civ^!n(Elp7A8J-VpxUcj*4C!gZ}-+_$pk6=omvTDNIIYKVXD(~#V|j~XUe<gfbj
za^b1s%SEb@MRhElS*P1&98}Zgg4K>OFX>mk5Xzuo(b+${-F6il!|hAcUbw!lSWz4E
zeAAQ(&lcDiYd?K`uawQy>ojAEj@^m9e1^MYY+Vn1_iqu4R_T$L7JJa}v6J}Kn79<H
z$LnNf2pN`IE)(+jT+}nU%XCA6;PK=*-3r&Wb^Q!`cCB5sHJ_vLZfg&p@9agBehIxk
zzbvn8a!5wNJV~L+|09xvrT?<Y%!sd@xHLs&de;}_2P$b@f?r((xeY(0Pm0>l(7V4U
z_vOKK_iG_%q#o9PG1$eT=a_e7wQRmp#Hpnd3#}%xP35}kqq8kWah}gs%g~og+t#S?
zmHlm=^l10upP_HLmnhBYbN+i>dIERt#wdq*Z2>2aR2pgtd}IlZWm9uWO=1`MC|1rV
zB(l@&ymIG$fs>(UzkJp()KqwB@Y`kE7L91lpJg`-9y*^|Xv4)jiEow7G^-QSn>*OI
zNL`%vr*K2bhP2gBJF{|6C>b7*@=$eI{k|g6&&)V#a!T@A5t~;V>y(n7$jv(B&~r*s
zC(t?Pbo#~_+T5{47Z|%c8>(2kx}-l%&ba;Nd)Kb_vlX8`&+E7mtS9wOjD1~tiU;S-
zscp9tzxyaMdHvX<`|;}R@(p_~9-hBqm+4=w{SBHeuP3VR^!Te8*2K!c(R0&cnVQun
zP2L`9T<prPH}ND#W!AmR9^E|a(qzK(Zz~<lk)0|#;lIk=TgE9t261OE@``%rth@U{
zFUViwh1*B&gCRTj3i$kB?PI?dn)l|q_USMVk*ggXk9wNaRtgCuHOB<~a{raEeuk8%
zQ*QJ!xoWXLYz`b@2UC1n&pTL|C3>8_JySdR&M}$7GU?*nrEUBCG<{36tzD$<E;_<C
z;a{V|ewKaV-QQp9v8?ve()Mujo+s{H$SSAmQz!O3kVEpu-3PZ7H;3-`KlPqbw~Z|;
z;qdK^AGc5EsaW`bxm1E%#KzDq6`N#K-9_?Dm3D4lscwGhPS>=V>P@ldtGA23QtGpB
z`Rb4}^LATycCdu)+@5mjb-!nJ-zs=@DA4~cYvdh0nUHI5R&yWM_5Xiu#o4(RE`RH|
zbh&QVk?au5sJZW7&aU6Zt}~ruLB!naWqDr9ZiQd>I#y@8?nKenb$O4yCTm>1x^Bfj
zpBG}LcD9~71pUJfMRalNu79qqp1znbuz&4=|7TN^L}e!yUrI8%Bevp)!HGRV!Mu~$
zBj$;9O4--io~_=j;=bdRJQsV{7Q2Zx>SuqiN%cKyFpKB^D)nmS8$2I)e|0}-nPu!?
zev#XwPBMJ@%kbd*(8o5Hq!hbXgiVzVj`~);I>@#zKioW?gQ0ZQ1owIowX&4(!vP@`
z8(oc`hqxL?YtM^xbb7XEmhHc3bN**`n_2|^Xx-X&N;AArz@)hAFYEjTkzB65yh3gU
zLJJDpmwUZxx9$&ETJz)CN!|EQr_F!Anyq~%P*glGt|I%~<?6Z(->!Y}Kk+~F`?)_?
z_5~Y1uNG8z{eN{x?D@@$Z|7E<-@cM{dE36Nx9qMr@yTERW24%1O+dR=z}?J-e}7Tr
zk?pTmZ`3Z;6D^crW!Y%{_h9pl`Cfl>j~s95XPlb7N#UZ=p|j=|hBmei+f0N^YyPf1
zboY;^^SfUAuHFSj$}z59Z|0^1*a>VpRIFL3?B@7hK(=bn>j_qxyMI4;E%xgU+nZlM
zeQg#VVu~wc%lvlD?}g3scGdj1Q<T3|eA1scy(H+#ii}F((0#1Noq2k}yo+3;G&Cx&
z@;X(%>5K6X5psOHh$&X_+wHDx7FRt3HE!PV`lxYuxuv~n!M}C&$$GmNe~h{1AYiq7
zOW*UVTf82x1*Lo+lx`Gub1?oUGSTLC+7vs}O}DzFHHzA-mNHELaCl=tZO}`m@D-me
zuAEm@=P=nHefsk@iyOIf6lPcz{|esA_vvVZ0+R=O)L!ezbNZzn@5Sc3w?#fQG10LL
z(muYlwMJsqitBH;3%UO^i}~o&{!I75vj3M{5(DxSZ#M*9ih6rxfgcmg<D-RZ7?lq!
zn>eF={kEODtxMM_GO(8`CI5Rhv+d{H<1UU@Z#DY}uH7Q{YEsHKktmTR{T%_uYwXU>
zSle#0K00ByeBQ#2n(jm2mE~>QwB9jFnAWPcY9>9>O^T@s_;}jvfV~OhKefzy{oh|5
z0va;@Ue0J`d+>7cSB>)m+p3P#R^Jvm{^G>k&ryeG-dp%tI7jU%XR1n~<hMr;K5uY;
zx=GIJxN1uO^(>CbxBV&v-uyiE;{S!hJKN9y>uE46VCoEAsTp{fqw=J5qpIa(4x6H9
z-{KmwLs(uJi935d_fYzixckNO#Xo;>|K)GKdged#0+oGA8<hB-FfcHxWI8(scse`7
z##a~^D(1|cVC!`xK%{kl&;E|Ee<Dhg4wX4LF5e_6+O$C6F^Ah~pGgM~s<f#be!zF9
zQ!kCnv4B0#;oz;$i&|SVYYHsor&ph^eD8SL@6eVx_u5QWa~*&EYx=*bw%?EO%wDtU
zw4cDV%Il^688-9RN4Uk`I<~s<-Fai_qp2EmS8OuS*V^#y=-k4o(qS8mZhz)G%AaKI
zwb|p3%Dt{jCpO&tdMQQZU&XB3I>&Av-T890S%Jkf5%Jk}%55qkTb9nWaxb|t*-C`v
zPOy-daMIUfvsUg9oLW?5A@J2@m#_1^s=3d#mg@H}+~addcBPH4+_U=q*X_=|^3Bej
z#O@lfeObubgADiN6Hm53O^?`-cVgE6zuVeo`q_RyZ(hfkG$&w@;|)D~1_oZK%#etZ
z2wxwo<osN{#FYG`RK1Ga0?>dJgN=PfL1J=tYKcNpYI<g#t<wA5`;`<jQ*4#OO??e~
z3!GCkGCit_QiH4fLXxYplPyz}?CrR0Dy%AU3vyFS5)~?PbMlI<Dr}X&=33<yTY-f2
zl@!2AO0sR0B76fBob!uP6-@L@^bB-eD>BQBl<c?^Y>HCStb$zJphgs>q}eKEl#~=$
z>Fbx5m+O@q>*W`v>l<2HTIw4Z=^Gj87Nw-=7FXt#Bv$C=6)S^`fSBQuTAW;zSx}Oh
zpQivaGchT@w8U0PiAzC20cv<fZh@~a-1xj=XecD-=jvA^7U&!58R+Nc>LcqYE^#d@
z!LOq@q_QAYKPa_0zqBYh)wL`&uS6Ny2zO840JsStFQsSZ7Jxlfkz0UbSxPcoT|sG4
z4p>b}vVLk#YHn&?NwL16o*|0*lFT$jV4?U1Bm?#ivO_9z3*hFXsD}9+tQZ_dRxbI;
zr65Z^U2K&=?zKwEPtHuS0y9l5O-+r`5>s_e4Gqk6O;U}_brVe!lXVSD4N?sf%@Y%o
zEDVv1^2{qPNz6-51sPS5TcDSjnPO#RYGh$<W?`;tm}Zu!YhqxYtecdUY^s}Tnv|Ms
zW|Wd-VPK48gnv<HdS+f?4zjC2Mx|t?SXri87#pP}ChA%kn;Pkwq#9W0T3Du9=q4K)
zC8t>!ni(aU7=euf#e<b&fTyjJk)8oUBp@fTBrU%v*H+0VGq1D)ArX?9n;Kk_2nsYq
zGb1AtLqjtI10xGdOG`6^qOjDW;>`R!keP-CdPXK-nUrKJxBQ~q#1dPj%-qEERQ-ay
zbg%>{TCDtwGE?(P5<%(IRtX+ZR?bDKi6!|(A^G_^wn`wkDj4Y*8i13xf(;~6JTi+*
z@{20%z$qG>Z-P?`Av{PH$jJmtDJX!GtyN+&#FFB~veXo?0~BCV$r*`x>8UBUO3(xa
z6VJpFS%zjQscELkmbw-x$p*S6W@(AKNg!@=qM?bAQL0g*xluBz>BaeJCFO}lsgCKX
zc_p?=?wPp-VBaWcfWuf5RdsnrDkuOAj0|-ROm&S7LX3>9OpUFK473f5tPBj4^dVlf
z(Ff&fn2&ArF~S6*0I3+T<5GZ#1-ZD{aoOmD%OX&@1Thd)F3=J~<ARn}C@73tLQ?pS
z2G?kCkrV<XDIQH-qrpW|2#}<BG<8uexVR7<p47Y)TcvU(d%L?gR<2`UVBjq9h%9Dc
z;5!1sj8nDwq!}0(*h@TpUD=<q%88na1l0bJWnf^4_jGX#k+__D(UFfqk%Q%++~fcI
zZK@MY8!9$8$u2fJd@EJm;J_@m)mb6`90Gk=kFL1F?Z&ipwx;Pd>7(J?a%&}GeZSXE
znEUp@%UudahYAn;ti9WJze%gDjHlGu<li%Pvq*6xe>X)3Q2)%+)z4*}Q$iB}Ph|8L

literal 0
HcmV?d00001