src/cpml/cpml-curve.c

   1 /* CPML - Cairo Path Manipulation Library
   2  * Copyright (C) 2007,2008,2009,2010,2011,2012,2013  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  * <refsect2 id="offset">
  52  * <title>Offseting algorithm</title>
  53  * <para>
  54  * Given a cubic Bézier primitive, it must be found the approximated
  55  * Bézier curve parallel to the original one at a specific distance
  56  * (the so called "offset curve"). The four points needed to build
  57  * the new curve must be returned.
  58  *
  59  * To solve the offset problem, a custom algorithm is used. First, the
  60  * resulting curve MUST have the same slope at the start and end point.
  61  * These constraints are not sufficient to resolve the system, so I let
  62  * the curve pass thought a given point (pm, known and got from the
  63  * original curve) at a given time (m, now hardcoded to 0.5).
  64  *
  65  * Firstly, some useful variables are defined:
  66  *
  67  * |[
  68  * v0 = unitvector(p[1] &minus; p[0]) &times; offset;
  69  * v3 = unitvector(p[3] &minus; p[2]) &times; offset;
  70  * p0 = p[0] + normal v0;
  71  * p3 = p[3] + normal v3.
  72  * ]|
  73  *
  74  * The resulting curve must have the same slopes than the original
  75  * one at the start and end points. Forcing the same slopes means:
  76  *
  77  * |[
  78  * p1 = p0 + k0 v0.
  79  * ]|
  80  *
  81  * where %k0 is an arbitrary factor. Decomposing for %x and %y:
  82  *
  83  * |[
  84  * p1.x = p0.x + k0 v0.x;
  85  * p1.y = p0.y + k0 v0.y.
  86  * ]|
  87  *
  88  * and doing the same for the end point:
  89  *
  90  * |[
  91  * p2.x = p3.x + k3 v3.x;
  92  * p2.y = p3.y + k3 v3.y.
  93  * ]|
  94  *
  95  * This does not give a resolvable system, though. The curve will be
  96  * interpolated by forcing its path to pass throught %pm when
  97  * <varname>time</varname> is %m, where <code>0 &le; m &le; 1</code>.
  98  * Knowing the function of the cubic Bézier:
  99  *
 100  * |[
 101  * C(t) = (1 &minus; t)³p0 + 3t(1 &minus; t)²p1 + 3t²(1 &minus; t)p2 + t³p3.
 102  * ]|
 103  *
 104  * and forcing <code>t = m</code> and <code>C(t) = pm</code>:
 105  *
 106  * |[
 107  * pm = (1 &minus; m)³p0 + 3m(1 &minus; m)²p1 + 3m²(1 &minus; m)p2 + m³p3.
 108  *
 109  * (1 &minus; m) p1 + m p2 = (pm &minus; (1 &minus; m)³p0 &minus; m³p3) / (3m (1 &minus; m)).
 110  * ]|
 111  *
 112  * gives this final system:
 113  *
 114  * |[
 115  * p1.x = p0.x + k0 v0.x;
 116  * p1.y = p0.y + k0 v0.y;
 117  * p2.x = p3.x + k3 v3.x;
 118  * p2.y = p3.y + k3 v3.y;
 119  * (1 &minus; m) p1.x + m p2.x = (pm.x &minus; (1 &minus; m)³p0.x &minus; m³p3.x) / (3m (1 &minus; m));
 120  * (1 &minus; m) p1.y + m p2.y = (pm.y &minus; (1 &minus; m)³p0.y &minus; m³p3.y) / (3m (1 &minus; m)).
 121  * ]|
 122  *
 123  * Substituting and resolving for %k0 and %k3:
 124  *
 125  * |[
 126  * (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;
 127  * (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.
 128  *
 129  * (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));
 130  * (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)).
 131  * ]|
 132  *
 133  * Letting:
 134  *
 135  * |[
 136  * pk = (pm &minus; (1 &minus; m)²(1+2m) p0 &minus; m²(3 &minus; 2m) p3) / (3m (1 &minus; m)).
 137  * ]|
 138  *
 139  * reduces the above to this final equations:
 140  *
 141  * |[
 142  * (1 &minus; m) k0 v0.x + m k3 v3.x = pk.x;
 143  * (1 &minus; m) k0 v0.y + m k3 v3.y = pk.y.
 144  * ]|
 145  *
 146  * If <code>v0.x &ne; 0</code>, the system can be resolved for %k0 and
 147  * %k3 calculated accordingly:
 148  *
 149  * |[
 150  * k0 = (pk.x &minus; m k3 v3.x) / ((1 &minus; m) v0.x);
 151  * (pk.x &minus; m k3 v3.x) v0.y / v0.x + m k3 v3.y = pk.y.
 152  *
 153  * k0 = (pk.x &minus; m k3 v3.x) / ((1 &minus; m) v0.x);
 154  * k3 m (v3.y &minus; v3.x v0.y / v0.x) = pk.y &minus; pk.x v0.y / v0.x.
 155  *
 156  * k3 = (pk.y &minus; pk.x v0.y / v0.x) / (m (v3.y &minus; v3.x v0.y / v0.x));
 157  * k0 = (pk.x &minus; m k3 v3.x) / ((1 &minus; m) v0.x).
 158  * ]|
 159  *
 160  * Otherwise, if <code>v3.x &ne; 0</code>, the system can be solved
 161  * for %k3 and %k0 calculated accordingly:
 162  *
 163  * |[
 164  * k3 = (pk.x &minus; (1 &minus; m) k0 v0.x) / (m v3.x);
 165  * (1 &minus; m) k0 v0.y + (pk.x &minus; (1 &minus; m) k0 v0.x) v3.y / v3.x = pk.y.
 166  *
 167  * k3 = (pk.x &minus; (1 &minus; m) k0 v0.x) / (m v3.x);
 168  * k0 (1 &minus; m) (v0.y &minus; k0 v0.x v3.y / v3.x) = pk.y &minus; pk.x v3.y / v3.x.
 169  *
 170  * k0 = (pk.y &minus; pk.x v3.y / v3.x) / ((1 &minus; m) (v0.y &minus; v0.x v3.y / v3.x));
 171  * k3 = (pk.x &minus; (1 &minus; m) k0 v0.x) / (m v3.x).
 172  * ]|
 173  *
 174  * The whole process must be guarded against division by 0 exceptions.
 175  * If either <code>v0.x</code> and <code>v3.x</code> are %0, the first
 176  * equation will be inconsistent. More in general, the
 177  * <code>v0.x &times; v3.y = v3.x &times; v3.y</code> condition must
 178  * be avoided. This is the first situation to avoid, in which case
 179  * an alternative approach should be used.
 180  *
 181  * </para>
 182  * </refsect2>
 183  *
 184  * Since: 1.0
 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  * Since: 1.0
 238  **/
 239 void
 240 cpml_curve_put_pair_at_time(const CpmlPrimitive *curve, double t,
 241                             CpmlPair *pair)
 242 {
 243     CpmlPair p1, p2, p3, p4;
 244     double t_2, t_3, t1, t1_2, t1_3;
 245
 246     cpml_primitive_put_point(curve, 0, &p1);
 247     cpml_primitive_put_point(curve, 1, &p2);
 248     cpml_primitive_put_point(curve, 2, &p3);
 249     cpml_primitive_put_point(curve, 3, &p4);
 250
 251     t_2 = t * t;
 252     t_3 = t_2 * t;
 253     t1 = 1 - t;
 254     t1_2 = t1 * t1;
 255     t1_3 = t1_2 * t1;
 256
 257     pair->x = t1_3 * p1.x + 3 * t1_2 * t * p2.x
 258               + 3 * t1 * t_2 * p3.x + t_3 * p4.x;
 259     pair->y = t1_3 * p1.y + 3 * t1_2 * t * p2.y
 260               + 3 * t1 * t_2 * p3.y + t_3 * p4.y;
 261 }
 262
 263 /**
 264  * cpml_curve_put_vector_at_time:
 265  * @curve:  the #CpmlPrimitive curve data
 266  * @t:      the "time" value
 267  * @vector: the destination vector
 268  *
 269  * Given the @curve Bézier cubic, finds the slope at time @t
 270  * (where 0 is the start and 1 is the end) and stores the result
 271  * in @vector. Keep in mind @t is not homogeneous, so 0.5
 272  * does not necessarily means the mid point.
 273  *
 274  * Since: 1.0
 275  **/
 276 void
 277 cpml_curve_put_vector_at_time(const CpmlPrimitive *curve,
 278                               double t, CpmlVector *vector)
 279 {
 280     CpmlPair p1, p2, p3, p4;
 281     CpmlPair p21, p32, p43;
 282     double t1, t1_2, t_2;
 283
 284     cpml_primitive_put_point(curve, 0, &p1);
 285     cpml_primitive_put_point(curve, 1, &p2);
 286     cpml_primitive_put_point(curve, 2, &p3);
 287     cpml_primitive_put_point(curve, 3, &p4);
 288
 289     p21.x = p2.x - p1.x;
 290     p21.y = p2.y - p1.y;
 291     p32.x = p3.x - p2.x;
 292     p32.y = p3.y - p2.y;
 293     p43.x = p4.x - p3.x;
 294     p43.y = p4.y - p3.y;
 295
 296     t1 = 1 - t;
 297     t1_2 = t1 * t1;
 298     t_2 = t * t;
 299
 300     vector->x = 3 * t1_2 * p21.x + 6 * t1 * t * p32.x + 3 * t_2 * p43.x;
 301     vector->y = 3 * t1_2 * p21.y + 6 * t1 * t * p32.y + 3 * t_2 * p43.y;
 302 }
 303
 304
 305 static void
 306 put_extents(const CpmlPrimitive *curve, CpmlExtents *extents)
 307 {
 308     CpmlPair p1, p2, p3, p4;
 309
 310     extents->is_defined = 0;
 311
 312     cpml_primitive_put_point(curve, 0, &p1);
 313     cpml_primitive_put_point(curve, 1, &p2);
 314     cpml_primitive_put_point(curve, 2, &p3);
 315     cpml_primitive_put_point(curve, 3, &p4);
 316
 317     cpml_extents_pair_add(extents, &p1);
 318     cpml_extents_pair_add(extents, &p2);
 319     cpml_extents_pair_add(extents, &p3);
 320     cpml_extents_pair_add(extents, &p4);
 321 }
 322
 323 static void
 324 offset(CpmlPrimitive *curve, double offset)
 325 {
 326     double m, mm;
 327     CpmlVector v0, v3, vm, vtmp;
 328     CpmlPair p0, p1, p2, p3, pm;
 329
 330     m = 0.5;
 331     mm = 1-m;
 332
 333     /* Firstly, convert the curve points from cairo format to cpml format
 334      * and store them (temporary) in p0..p3 */
 335     cpml_pair_from_cairo(&p0, curve->org);
 336     cpml_pair_from_cairo(&p1, &curve->data[1]);
 337     cpml_pair_from_cairo(&p2, &curve->data[2]);
 338     cpml_pair_from_cairo(&p3, &curve->data[3]);
 339
 340     /* v0 = p1-p0 */
 341     v0.x = p1.x - p0.x;
 342     v0.y = p1.y - p0.y;
 343
 344     /* v3 = p3-p2 */
 345     v3.x = p3.x - p2.x;
 346     v3.y = p3.y - p2.y;
 347
 348     /* pm = point in C(m) offseted the requested @offset distance */
 349     cpml_curve_put_vector_at_time(curve, m, &vm);
 350     cpml_vector_set_length(&vm, offset);
 351     cpml_vector_normal(&vm);
 352     cpml_curve_put_pair_at_time(curve, m, &pm);
 353     pm.x += vm.x;
 354     pm.y += vm.y;
 355
 356     /* p0 = p0 + normal of v0 of @offset magnitude (exact value) */
 357     cpml_pair_copy(&vtmp, &v0);
 358     cpml_vector_set_length(&vtmp, offset);
 359     cpml_vector_normal(&vtmp);
 360     p0.x += vtmp.x;
 361     p0.y += vtmp.y;
 362
 363     /* p3 = p3 + normal of v3 of @offset magnitude, as done for p0 */
 364     cpml_pair_copy(&vtmp, &v3);
 365     cpml_vector_set_length(&vtmp, offset);
 366     cpml_vector_normal(&vtmp);
 367     p3.x += vtmp.x;
 368     p3.y += vtmp.y;
 369
 370     if (v0.x*v3.y == v3.x*v0.y) {
 371         /* Inconsistent equations: use the alternative approach */
 372         p1.x = p0.x + v0.x + vm.x * 4/3;
 373         p1.y = p0.y + v0.y + vm.y * 4/3;
 374         p2.x = p3.x - v3.x + vm.x * 4/3;
 375         p2.y = p3.y - v3.y + vm.y * 4/3;
 376     } else {
 377         CpmlPair pk;
 378         double k0, k3;
 379
 380         pk.x = (pm.x - mm*mm*(1+m+m)*p0.x - m*m*(1+mm+mm)*p3.x) / (3*m*(1-m));
 381         pk.y = (pm.y - mm*mm*(1+m+m)*p0.y - m*m*(1+mm+mm)*p3.y) / (3*m*(1-m));
 382
 383         if (v0.x != 0) {
 384             k3 = (pk.y - pk.x*v0.y / v0.x) / (m*(v3.y - v3.x*v0.y / v0.x));
 385             k0 = (pk.x - m*k3*v3.x) / (mm*v0.x);
 386         } else {
 387             k0 = (pk.y - pk.x*v3.y / v3.x) / (mm*(v0.y - v0.x*v3.y / v3.x));
 388             k3 = (pk.x - mm*k0*v0.x) / (m*v3.x);
 389         }
 390
 391         p1.x = p0.x + k0*v0.x;
 392         p1.y = p0.y + k0*v0.y;
 393         p2.x = p3.x + k3*v3.x;
 394         p2.y = p3.y + k3*v3.y;
 395     }
 396
 397     /* Return the new curve in the original array */
 398     cpml_pair_to_cairo(&p0, curve->org);
 399     cpml_pair_to_cairo(&p1, &curve->data[1]);
 400     cpml_pair_to_cairo(&p2, &curve->data[2]);
 401     cpml_pair_to_cairo(&p3, &curve->data[3]);
 402 }