Implement Bezier3 simplification

contour_of.c

@@ -210,10 +210,6 @@ int main(int argc, char** argv) {
         }
     }
 
-    if (type_dp == DPBezier3) {
-        ERREUR_FATALE("Bezier3 non implémenté");
-    }
-
     Image i = lire_fichier_image(image_name);
     Mask mask = contour_init_mask(i);
 
@@ -288,8 +284,11 @@ int main(int argc, char** argv) {
             }
             case DPBezier2:
             case DPBezier3: {
-                // TODO gérer Bezier3
-                cb = simplification_douglas_peucker_bezier2(c_tab, 0, c_tab.len - 1, seuil_dp);
+                if (type_dp == DPBezier2) {
+                    cb = simplification_douglas_peucker_bezier2(c_tab, 0, c_tab.len - 1, seuil_dp);
+                } else {
+                    cb = simplification_douglas_peucker_bezier3(c_tab, 0, c_tab.len - 1, seuil_dp);
+                }
                 generic_contour = (GenericContour) {
                     .tag = DPBezier3,
                     .val = { .bezier3 = cb },

geom2d.c

@@ -149,19 +149,61 @@ Bezier2 approx_bezier2(Point *start, unsigned int len) {
         curve.c1 = mul_reel_point(add_point(start[0], start[1]), .5);
         curve.c2 = start[1];
     } else {
-        float n = len - 1;
+        double n = len - 1;
 
         Point sum = set_point(0, 0);
         for (int i = 1; i < (len-1); i++) {
             sum = add_point(sum, start[i]);
         }
 
-        float alpha = (3. * n) / (n*n - 1.);
-        float beta = (1. - (2. * n)) / (2. * (n + 1.));
+        double alpha = (3. * n) / (n*n - 1.);
+        double beta = (1. - (2. * n)) / (2. * (n + 1.));
 
         curve.c2 = start[len-1];
         curve.c1 = add_point(mul_reel_point(sum, alpha), mul_reel_point(add_point(curve.c0, curve.c2), beta));
     }
 
     return curve;
 }
+
+
+double approx_bezier3_gamma(double k, double n) {
+    return (6. * pow(k, 4.)) - (8. * n * pow(k, 3.)) + (6. * k * k) - (4. * n * k) + pow(n, 4.) - (n * n);
+}
+
+
+Bezier3 approx_bezier3(Point* start, unsigned int len) {
+    assert(len > 1);
+
+    if (len <= 3) {
+        Bezier2 b2 = approx_bezier2(start, len);
+        return bezier2_to_bezier3(&b2);
+    } else {
+        Bezier3 curve = {
+                .c0 = start[0],
+                .c3 = start[len-1],
+        };
+
+        double n = len - 1;
+
+        Point sum1 = set_point(0, 0);
+        for (int i = 1; i < (len-1); i++) {
+            sum1 = add_point(sum1, mul_reel_point(start[i], approx_bezier3_gamma(i, n)));
+        }
+
+        Point sum2 = set_point(0, 0);
+        for (int i = 1; i < (len-1); i++) {
+            sum2 = add_point(sum2, mul_reel_point(start[i], approx_bezier3_gamma(n - (double) i, n)));
+        }
+
+        double d = 3. * (n + 2.) * (3. * n * n + 1.);
+        double alpha = ((-15. * n * n * n) + (5. * n * n) + (2. * n) + 4.) / d;
+        double beta = ((10. * n * n * n) - (15. * n * n) + n + 2.) / d;
+        double lambda = (70. * n) / (3. * (n * n - 1.) * (n * n - 4.) * (3. * n * n + 1.));
+
+        curve.c1 = add_point(mul_reel_point(curve.c0, alpha), add_point(mul_reel_point(sum1, lambda), mul_reel_point(curve.c3, beta)));
+        curve.c2 = add_point(mul_reel_point(curve.c0, beta), add_point(mul_reel_point(sum2, lambda), mul_reel_point(curve.c3, alpha)));
+
+        return curve;
+    }
+}

geom2d.h

@@ -117,4 +117,11 @@ Bezier3 bezier2_to_bezier3(Bezier2* self);
  */
 Bezier2 approx_bezier2(Point* start, unsigned int len);
 
+
+/**
+ * Approxime la séquence de points start[0]..start[len-1] par une courbe de Bézier de degré 3 dont les points de
+ * contrôle extrêmes sont le début et la fin de la séquence.
+ */
+Bezier3 approx_bezier3(Point* start, unsigned int len);
+
 #endif

simplification.c

@@ -37,6 +37,39 @@ ListeBezier3 simplification_douglas_peucker_bezier2(TableauPoints c, unsigned in
 }
 
 
+ListeBezier3 simplification_douglas_peucker_bezier3(TableauPoints c, unsigned int i1, unsigned int i2, double seuil) {
+    assert(i1 < i2);
+
+    unsigned int n = i2 - i1;
+    unsigned int np1 = n + 1;
+    Bezier3 bezier = approx_bezier3(&c.inner[i1], np1);
+
+    double distance_max = 0;
+    int index_distance_max = i1;
+    for (int i = i1 + 1; i <= i2; i++) {
+        int j = i - i1;
+        double distance = distance_point_bezier3(&bezier, c.inner[i], (float) j / (float) n);
+        if (distance_max < distance) {
+            distance_max = distance;
+            index_distance_max = i;
+        }
+    }
+
+    ListeBezier3 L;
+    if (distance_max <= seuil) {
+        L = liste_bezier3_new();
+        liste_bezier3_push(&L, bezier);
+    } else {
+        L = simplification_douglas_peucker_bezier3(c, i1, index_distance_max, seuil);
+        ListeBezier3 L2 = simplification_douglas_peucker_bezier3(c, index_distance_max, i2, seuil);
+
+        liste_bezier3_concat(&L, L2);
+    }
+
+    return L;
+}
+
+
 ListePoint simplification_douglas_peucker(TableauPoints c, unsigned int i1, unsigned int i2, double seuil) {
     double distance_max = 0;
     int index_distance_max = i1;

simplification.h

@@ -6,16 +6,18 @@
 
 
 /**
- * Renvoie une courbe de bézier approximant le polygone définit par les `len` points à partir de `start`
+ * Simplifie une séquence de points avec l'algo. de Douglas-Peucker, entre un intervalle i1..i2, pour un seuil donné, en
+ * utilisant des courbes de Bézier de degré 2 pour remplacer les sous-séquences de points à simplifier.
  */
-Bezier2 approx_bezier2(Point *start, unsigned int len);
+ListeBezier3 simplification_douglas_peucker_bezier2(TableauPoints c, unsigned int i1, unsigned int i2, double seuil);
 
 
 /**
  * Simplifie une séquence de points avec l'algo. de Douglas-Peucker, entre un intervalle i1..i2, pour un seuil donné, en
  * utilisant des courbes de Bézier de degré 2 pour remplacer les sous-séquences de points à simplifier.
  */
-ListeBezier3 simplification_douglas_peucker_bezier2(TableauPoints c, unsigned int i1, unsigned int i2, double seuil);
+ListeBezier3 simplification_douglas_peucker_bezier3(TableauPoints c, unsigned int i1, unsigned int i2, double seuil);
+
 
 /**
  * Simplifie une séquence de points avec l'algo. de Douglas-Peucker, entre un intervalle i1..i2, pour un seuil donné, en

test_geom2d.c

@@ -80,7 +80,6 @@ int main() {
             bezier2_C(&bezier2base, 1),
     };
     bezier2 = approx_bezier2(points, 5);
-    printf("(%f;%f)\n", bezier2.c1.x, bezier2.c1.y);
     assert(egaux_points(bezier2base.c0, bezier2.c0));
     assert(egaux_points(bezier2base.c1, bezier2.c1));
     assert(egaux_points(bezier2base.c2, bezier2.c2));
@@ -103,6 +102,16 @@ int main() {
     int ye5 = (2.452381 - bezier2.c1.y) * 100000;
     assert(xe5 == 0 && ye5 == 0);
 
+    // approx_bezier3
+
+    // n=9 (exemple du cours)
+    bezier3 = approx_bezier3(points2, 9);
+    int x1e5 = (1.737287 - bezier3.c1.x) * 100000;
+    int y1e5 = (0.929380 - bezier3.c1.y) * 100000;
+    int x2e5 = (1.844176 - bezier3.c2.x) * 100000;
+    int y2e5 = (3.489158 - bezier3.c2.y) * 100000;
+    assert(x1e5 == 0 && y1e5 == 0 && x2e5 == 0 && y2e5 == 0);
+
     printf("\e[32mtest_geom2d:bezier passé avec succès !\e[0m\n");
 
     printf("\e[32mtest_geom2d:* passé avec succès !\e[0m\n");