src/cpml/cpml-curve.c

   1 /* CPML - Cairo Path Manipulation Library
   2  * Copyright (C) 2008,2009,2010,2011  Nicola Fontana <ntd at entidi.it>
   3  *
   4  * This library is free software; you can redistribute it and/or
   5  * modify it under the terms of the GNU Lesser General Public
   6  * License as published by the Free Software Foundation; either
   7  * version 2 of the License, or (at your option) any later version.
   8  *
   9  * This library is distributed in the hope that it will be useful,
  10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12  * Lesser General Public License for more details.
  13  *
  14  * You should have received a copy of the GNU Lesser General Public
  15  * License along with this library; if not, write to the
  16  * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
  17  * Boston, MA  02110-1301, USA.
  18  */
  19
  20
  21 /**
  22  * SECTION:cpml-curve
  23  * @Section_Id:CpmlCurve
  24  * @title: CpmlCurve
  25  * @short_description: Bézier cubic curve primitive management
  26  *
  27  * The following functions manipulate %CPML_CURVE #CpmlPrimitive.
  28  * No validation is made on the input so use the following methods
  29  * only when you are sure the <varname>primitive</varname> argument
  30  * is effectively a cubic Bézier curve.
  31  *
  32  * <important>
  33  * <title>TODO</title>
  34  * <itemizedlist>
  35  * <listitem>the get_length() method must be implemented;</listitem>
  36  * <listitem>actually the put_extents() method is implemented by computing
  37  *           the bounding box of the control polygon and this will likely
  38  *           include some empty space: there is room for improvements;</listitem>
  39  * <listitem>the put_pair_at() method must be implemented;</listitem>
  40  * <listitem>the put_vector_at() method must be implemented;</listitem>
  41  * <listitem>the get_closest_pos() method must be implemented;</listitem>
  42  * <listitem>the put_intersections() method must be implemented;</listitem>
  43  * <listitem>by default, the offset curve is calculated by using the point
  44  *           at t=0.5 as reference: use a better candidate;</listitem>
  45  * <listitem>in the offset() implementation, when the equations are
  46  *           inconsistent, the alternative approach performs very bad
  47  *           if <varname>v0</varname> and <varname>v3</varname> are
  48  *           opposite or staggered.</listitem>
  49  * </itemizedlist>
  50  * </important>
  51  * </para>
  52  * <refsect2 id="offset">
  53  * <title>Offseting algorithm</title>
  54  * <para>
  55  * Given a cubic Bézier primitive, it must be found the approximated
  56  * Bézier curve parallel to the original one at a specific distance
  57  * (the so called "offset curve"). The four points needed to build
  58  * the new curve must be returned.
  59  *
  60  * To solve the offset problem, a custom algorithm is used. First, the
  61  * resulting curve MUST have the same slope at the start and end point.
  62  * These constraints are not sufficient to resolve the system, so I let
  63  * the curve pass thought a given point (pm, known and got from the
  64  * original curve) at a given time (m, now hardcoded to 0.5).
  65  *
  66  * Firstly, some useful variables are defined:
  67  *
  68  * |[
  69  * v0 = unitvector(p[1] &minus; p[0]) &times; offset;
  70  * v3 = unitvector(p[3] &minus; p[2]) &times; offset;
  71  * p0 = p[0] + normal v0;
  72  * p3 = p[3] + normal v3.
  73  * ]|
  74  *
  75  * The resulting curve must have the same slopes than the original
  76  * one at the start and end points. Forcing the same slopes means:
  77  *
  78  * |[
  79  * p1 = p0 + k0 v0.
  80  * ]|
  81  *
  82  * where %k0 is an arbitrary factor. Decomposing for %x and %y:
  83  *
  84  * |[
  85  * p1.x = p0.x + k0 v0.x;
  86  * p1.y = p0.y + k0 v0.y.
  87  * ]|
  88  *
  89  * and doing the same for the end point:
  90  *
  91  * |[
  92  * p2.x = p3.x + k3 v3.x;
  93  * p2.y = p3.y + k3 v3.y.
  94  * ]|
  95  *
  96  * This does not give a resolvable system, though. The curve will be
  97  * interpolated by forcing its path to pass throught %pm when
  98  * <varname>time</varname> is %m, where <code>0 &le; m &le; 1</code>.
  99  * Knowing the function of the cubic Bézier:
 100  *
 101  * |[
 102  * C(t) = (1 &minus; t)³p0 + 3t(1 &minus; t)²p1 + 3t²(1 &minus; t)p2 + t³p3.
 103  * ]|
 104  *
 105  * and forcing <code>t = m</code> and <code>C(t) = pm</code>:
 106  *
 107  * |[
 108  * pm = (1 &minus; m)³p0 + 3m(1 &minus; m)²p1 + 3m²(1 &minus; m)p2 + m³p3.
 109  *
 110  * (1 &minus; m) p1 + m p2 = (pm &minus; (1 &minus; m)³p0 &minus; m³p3) / (3m (1 &minus; m)).
 111  * ]|
 112  *
 113  * gives this final system:
 114  *
 115  * |[
 116  * p1.x = p0.x + k0 v0.x;
 117  * p1.y = p0.y + k0 v0.y;
 118  * p2.x = p3.x + k3 v3.x;
 119  * p2.y = p3.y + k3 v3.y;
 120  * (1 &minus; m) p1.x + m p2.x = (pm.x &minus; (1 &minus; m)³p0.x &minus; m³p3.x) / (3m (1 &minus; m));
 121  * (1 &minus; m) p1.y + m p2.y = (pm.y &minus; (1 &minus; m)³p0.y &minus; m³p3.y) / (3m (1 &minus; m)).
 122  * ]|
 123  *
 124  * Substituting and resolving for %k0 and %k3:
 125  *
 126  * |[
 127  * (1 &minus; m) k0 v0.x + m k3 v3.x = (pm.x &minus; (1 &minus; m)³p0.x &minus; m³p3.x) / (3m (1 &minus; m)) &minus; (1 &minus; m) p0.x &minus; m p3.x;
 128  * (1 &minus; m) k0 v0.y + m k3 v3.y = (pm.y &minus; (1 &minus; m)³p0.y &minus; m³p3.y) / (3m (1 &minus; m)) &minus; (1 &minus; m) p0.y &minus; m p3.y.
 129  *
 130  * (1 &minus; m) k0 v0.x + m k3 v3.x = (pm.x &minus; (1 &minus; m)²(1+2m) p0.x &minus; m²(3 &minus; 2m) p3.x) / (3m (1 &minus; m));
 131  * (1 &minus; m) k0 v0.y + m k3 v3.y = (pm.y &minus; (1 &minus; m)²(1+2m) p0.y &minus; m²(3 &minus; 2m) p3.y) / (3m (1 &minus; m)).
 132  * ]|
 133  *
 134  * Letting:
 135  *
 136  * |[
 137  * pk = (pm &minus; (1 &minus; m)²(1+2m) p0 &minus; m²(3 &minus; 2m) p3) / (3m (1 &minus; m)).
 138  * ]|
 139  *
 140  * reduces the above to this final equations:
 141  *
 142  * |[
 143  * (1 &minus; m) k0 v0.x + m k3 v3.x = pk.x;
 144  * (1 &minus; m) k0 v0.y + m k3 v3.y = pk.y.
 145  * ]|
 146  *
 147  * If <code>v0.x &ne; 0</code>, the system can be resolved for %k0 and
 148  * %k3 calculated accordingly:
 149  *
 150  * |[
 151  * k0 = (pk.x &minus; m k3 v3.x) / ((1 &minus; m) v0.x);
 152  * (pk.x &minus; m k3 v3.x) v0.y / v0.x + m k3 v3.y = pk.y.
 153  *
 154  * k0 = (pk.x &minus; m k3 v3.x) / ((1 &minus; m) v0.x);
 155  * k3 m (v3.y &minus; v3.x v0.y / v0.x) = pk.y &minus; pk.x v0.y / v0.x.
 156  *
 157  * k3 = (pk.y &minus; pk.x v0.y / v0.x) / (m (v3.y &minus; v3.x v0.y / v0.x));
 158  * k0 = (pk.x &minus; m k3 v3.x) / ((1 &minus; m) v0.x).
 159  * ]|
 160  *
 161  * Otherwise, if <code>v3.x &ne; 0</code>, the system can be solved
 162  * for %k3 and %k0 calculated accordingly:
 163  *
 164  * |[
 165  * k3 = (pk.x &minus; (1 &minus; m) k0 v0.x) / (m v3.x);
 166  * (1 &minus; m) k0 v0.y + (pk.x &minus; (1 &minus; m) k0 v0.x) v3.y / v3.x = pk.y.
 167  *
 168  * k3 = (pk.x &minus; (1 &minus; m) k0 v0.x) / (m v3.x);
 169  * k0 (1 &minus; m) (v0.y &minus; k0 v0.x v3.y / v3.x) = pk.y &minus; pk.x v3.y / v3.x.
 170  *
 171  * k0 = (pk.y &minus; pk.x v3.y / v3.x) / ((1 &minus; m) (v0.y &minus; v0.x v3.y / v3.x));
 172  * k3 = (pk.x &minus; (1 &minus; m) k0 v0.x) / (m v3.x).
 173  * ]|
 174  *
 175  * The whole process must be guarded against division by 0 exceptions.
 176  * If either <code>v0.x</code> and <code>v3.x</code> are %0, the first
 177  * equation will be inconsistent. More in general, the
 178  * <code>v0.x &times; v3.y = v3.x &times; v3.y</code> condition must
 179  * be avoided. This is the first situation to avoid, in which case
 180  * an alternative approach should be used.
 181  *
 182  * </para>
 183  * </refsect2>
 184  * <para>
 185  **/
 186
 187
 188 #include "cpml-internal.h"
 189 #include "cpml-extents.h"
 190 #include "cpml-segment.h"
 191 #include "cpml-primitive.h"
 192 #include "cpml-primitive-private.h"
 193 #include "cpml-curve.h"
 194
 195
 196 static void     put_extents             (const CpmlPrimitive    *curve,
 197                                          CpmlExtents            *extents);
 198 static void     offset                  (CpmlPrimitive          *curve,
 199                                          double                  offset);
 200
 201
 202 const _CpmlPrimitiveClass *
 203 _cpml_curve_get_class(void)
 204 {
 205     static _CpmlPrimitiveClass *p_class = NULL;
 206
 207     if (p_class == NULL) {
 208         static _CpmlPrimitiveClass class_data = {
 209             "curve to", 4,
 210             NULL,
 211             put_extents,
 212             NULL,
 213             NULL,
 214             NULL,
 215             NULL,
 216             offset,
 217             NULL
 218         };
 219         p_class = &class_data;
 220     }
 221
 222     return p_class;
 223 }
 224
 225
 226 /**
 227  * cpml_curve_put_pair_at_time:
 228  * @curve: the #CpmlPrimitive curve data
 229  * @t:     the "time" value
 230  * @pair:  the destination pair
 231  *
 232  * Given the @curve Bézier cubic, finds the coordinates at time @t
 233  * (where 0 is the start and 1 is the end) and stores the result
 234  * in @pair. Keep in mind @t is not homogeneous, so 0.5 does not
 235  * necessarily means the mid point.
 236  **/
 237 void
 238 cpml_curve_put_pair_at_time(const CpmlPrimitive *curve, double t,
 239                             CpmlPair *pair)
 240 {
 241     cairo_path_data_t *p1, *p2, *p3, *p4;
 242     double t_2, t_3, t1, t1_2, t1_3;
 243
 244     p1 = cpml_primitive_get_point(curve, 0);
 245     p2 = cpml_primitive_get_point(curve, 1);
 246     p3 = cpml_primitive_get_point(curve, 2);
 247     p4 = cpml_primitive_get_point(curve, 3);
 248
 249     t_2 = t * t;
 250     t_3 = t_2 * t;
 251     t1 = 1 - t;
 252     t1_2 = t1 * t1;
 253     t1_3 = t1_2 * t1;
 254
 255     pair->x = t1_3 * p1->point.x + 3 * t1_2 * t * p2->point.x
 256               + 3 * t1 * t_2 * p3->point.x + t_3 * p4->point.x;
 257     pair->y = t1_3 * p1->point.y + 3 * t1_2 * t * p2->point.y
 258               + 3 * t1 * t_2 * p3->point.y + t_3 * p4->point.y;
 259 }
 260
 261 /**
 262  * cpml_curve_put_vector_at_time:
 263  * @curve:  the #CpmlPrimitive curve data
 264  * @t:      the "time" value
 265  * @vector: the destination vector
 266  *
 267  * Given the @curve Bézier cubic, finds the slope at time @t
 268  * (where 0 is the start and 1 is the end) and stores the result
 269  * in @vector. Keep in mind @t is not homogeneous, so 0.5
 270  * does not necessarily means the mid point.
 271  **/
 272 void
 273 cpml_curve_put_vector_at_time(const CpmlPrimitive *curve,
 274                               double t, CpmlVector *vector)
 275 {
 276     cairo_path_data_t *p1, *p2, *p3, *p4;
 277     CpmlPair p21, p32, p43;
 278     double t1, t1_2, t_2;
 279
 280     p1 = cpml_primitive_get_point(curve, 0);
 281     p2 = cpml_primitive_get_point(curve, 1);
 282     p3 = cpml_primitive_get_point(curve, 2);
 283     p4 = cpml_primitive_get_point(curve, 3);
 284
 285     p21.x = p2->point.x - p1->point.x;
 286     p21.y = p2->point.y - p1->point.y;
 287     p32.x = p3->point.x - p2->point.x;
 288     p32.y = p3->point.y - p2->point.y;
 289     p43.x = p4->point.x - p3->point.x;
 290     p43.y = p4->point.y - p3->point.y;
 291
 292     t1 = 1 - t;
 293     t1_2 = t1 * t1;
 294     t_2 = t * t;
 295
 296     vector->x = 3 * t1_2 * p21.x + 6 * t1 * t * p32.x + 3 * t_2 * p43.x;
 297     vector->y = 3 * t1_2 * p21.y + 6 * t1 * t * p32.y + 3 * t_2 * p43.y;
 298 }
 299
 300
 301 static void
 302 put_extents(const CpmlPrimitive *curve, CpmlExtents *extents)
 303 {
 304     CpmlPair p1, p2, p3, p4;
 305
 306     extents->is_defined = 0;
 307
 308     cpml_pair_from_cairo(&p1, cpml_primitive_get_point(curve, 0));
 309     cpml_pair_from_cairo(&p2, cpml_primitive_get_point(curve, 1));
 310     cpml_pair_from_cairo(&p3, cpml_primitive_get_point(curve, 2));
 311     cpml_pair_from_cairo(&p4, cpml_primitive_get_point(curve, 3));
 312
 313     cpml_extents_pair_add(extents, &p1);
 314     cpml_extents_pair_add(extents, &p2);
 315     cpml_extents_pair_add(extents, &p3);
 316     cpml_extents_pair_add(extents, &p4);
 317 }
 318
 319 static void
 320 offset(CpmlPrimitive *curve, double offset)
 321 {
 322     double m, mm;
 323     CpmlVector v0, v3, vm, vtmp;
 324     CpmlPair p0, p1, p2, p3, pm;
 325
 326     m = 0.5;
 327     mm = 1-m;
 328
 329     /* Firstly, convert the curve points from cairo format to cpml format
 330      * and store them (temporary) in p0..p3 */
 331     cpml_pair_from_cairo(&p0, curve->org);
 332     cpml_pair_from_cairo(&p1, &curve->data[1]);
 333     cpml_pair_from_cairo(&p2, &curve->data[2]);
 334     cpml_pair_from_cairo(&p3, &curve->data[3]);
 335
 336     /* v0 = p1-p0 */
 337     v0.x = p1.x - p0.x;
 338     v0.y = p1.y - p0.y;
 339
 340     /* v3 = p3-p2 */
 341     v3.x = p3.x - p2.x;
 342     v3.y = p3.y - p2.y;
 343
 344     /* pm = point in C(m) offseted the requested @offset distance */
 345     cpml_curve_put_vector_at_time(curve, m, &vm);
 346     cpml_vector_set_length(&vm, offset);
 347     cpml_vector_normal(&vm);
 348     cpml_curve_put_pair_at_time(curve, m, &pm);
 349     pm.x += vm.x;
 350     pm.y += vm.y;
 351
 352     /* p0 = p0 + normal of v0 of @offset magnitude (exact value) */
 353     cpml_pair_copy(&vtmp, &v0);
 354     cpml_vector_set_length(&vtmp, offset);
 355     cpml_vector_normal(&vtmp);
 356     p0.x += vtmp.x;
 357     p0.y += vtmp.y;
 358
 359     /* p3 = p3 + normal of v3 of @offset magnitude, as done for p0 */
 360     cpml_pair_copy(&vtmp, &v3);
 361     cpml_vector_set_length(&vtmp, offset);
 362     cpml_vector_normal(&vtmp);
 363     p3.x += vtmp.x;
 364     p3.y += vtmp.y;
 365
 366     if (v0.x*v3.y == v3.x*v0.y) {
 367         /* Inconsistent equations: use the alternative approach */
 368         p1.x = p0.x + v0.x + vm.x * 4/3;
 369         p1.y = p0.y + v0.y + vm.y * 4/3;
 370         p2.x = p3.x - v3.x + vm.x * 4/3;
 371         p2.y = p3.y - v3.y + vm.y * 4/3;
 372     } else {
 373         CpmlPair pk;
 374         double k0, k3;
 375
 376         pk.x = (pm.x - mm*mm*(1+m+m)*p0.x - m*m*(1+mm+mm)*p3.x) / (3*m*(1-m));
 377         pk.y = (pm.y - mm*mm*(1+m+m)*p0.y - m*m*(1+mm+mm)*p3.y) / (3*m*(1-m));
 378
 379         if (v0.x != 0) {
 380             k3 = (pk.y - pk.x*v0.y / v0.x) / (m*(v3.y - v3.x*v0.y / v0.x));
 381             k0 = (pk.x - m*k3*v3.x) / (mm*v0.x);
 382         } else {
 383             k0 = (pk.y - pk.x*v3.y / v3.x) / (mm*(v0.y - v0.x*v3.y / v3.x));
 384             k3 = (pk.x - mm*k0*v0.x) / (m*v3.x);
 385         }
 386
 387         p1.x = p0.x + k0*v0.x;
 388         p1.y = p0.y + k0*v0.y;
 389         p2.x = p3.x + k3*v3.x;
 390         p2.y = p3.y + k3*v3.y;
 391     }
 392
 393     /* Return the new curve in the original array */
 394     cpml_pair_to_cairo(&p0, curve->org);
 395     cpml_pair_to_cairo(&p1, &curve->data[1]);
 396     cpml_pair_to_cairo(&p2, &curve->data[2]);
 397     cpml_pair_to_cairo(&p3, &curve->data[3]);
 398 }