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  * </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  * Since: 1.0
 187  **/
 188
 189
 190 #include "cpml-internal.h"
 191 #include "cpml-extents.h"
 192 #include "cpml-segment.h"
 193 #include "cpml-primitive.h"
 194 #include "cpml-primitive-private.h"
 195 #include "cpml-curve.h"
 196
 197
 198 static void     put_extents             (const CpmlPrimitive    *curve,
 199                                          CpmlExtents            *extents);
 200 static void     offset                  (CpmlPrimitive          *curve,
 201                                          double                  offset);
 202
 203
 204 const _CpmlPrimitiveClass *
 205 _cpml_curve_get_class(void)
 206 {
 207     static _CpmlPrimitiveClass *p_class = NULL;
 208
 209     if (p_class == NULL) {
 210         static _CpmlPrimitiveClass class_data = {
 211             "curve to", 4,
 212             NULL,
 213             put_extents,
 214             NULL,
 215             NULL,
 216             NULL,
 217             NULL,
 218             offset,
 219             NULL
 220         };
 221         p_class = &class_data;
 222     }
 223
 224     return p_class;
 225 }
 226
 227
 228 /**
 229  * cpml_curve_put_pair_at_time:
 230  * @curve: the #CpmlPrimitive curve data
 231  * @t:     the "time" value
 232  * @pair:  the destination pair
 233  *
 234  * Given the @curve Bézier cubic, finds the coordinates at time @t
 235  * (where 0 is the start and 1 is the end) and stores the result
 236  * in @pair. Keep in mind @t is not homogeneous, so 0.5 does not
 237  * necessarily means the mid point.
 238  *
 239  * Since: 1.0
 240  **/
 241 void
 242 cpml_curve_put_pair_at_time(const CpmlPrimitive *curve, double t,
 243                             CpmlPair *pair)
 244 {
 245     CpmlPair p1, p2, p3, p4;
 246     double t_2, t_3, t1, t1_2, t1_3;
 247
 248     cpml_primitive_put_point(curve, 0, &p1);
 249     cpml_primitive_put_point(curve, 1, &p2);
 250     cpml_primitive_put_point(curve, 2, &p3);
 251     cpml_primitive_put_point(curve, 3, &p4);
 252
 253     t_2 = t * t;
 254     t_3 = t_2 * t;
 255     t1 = 1 - t;
 256     t1_2 = t1 * t1;
 257     t1_3 = t1_2 * t1;
 258
 259     pair->x = t1_3 * p1.x + 3 * t1_2 * t * p2.x
 260               + 3 * t1 * t_2 * p3.x + t_3 * p4.x;
 261     pair->y = t1_3 * p1.y + 3 * t1_2 * t * p2.y
 262               + 3 * t1 * t_2 * p3.y + t_3 * p4.y;
 263 }
 264
 265 /**
 266  * cpml_curve_put_vector_at_time:
 267  * @curve:  the #CpmlPrimitive curve data
 268  * @t:      the "time" value
 269  * @vector: the destination vector
 270  *
 271  * Given the @curve Bézier cubic, finds the slope at time @t
 272  * (where 0 is the start and 1 is the end) and stores the result
 273  * in @vector. Keep in mind @t is not homogeneous, so 0.5
 274  * does not necessarily means the mid point.
 275  *
 276  * Since: 1.0
 277  **/
 278 void
 279 cpml_curve_put_vector_at_time(const CpmlPrimitive *curve,
 280                               double t, CpmlVector *vector)
 281 {
 282     CpmlPair p1, p2, p3, p4;
 283     CpmlPair p21, p32, p43;
 284     double t1, t1_2, t_2;
 285
 286     cpml_primitive_put_point(curve, 0, &p1);
 287     cpml_primitive_put_point(curve, 1, &p2);
 288     cpml_primitive_put_point(curve, 2, &p3);
 289     cpml_primitive_put_point(curve, 3, &p4);
 290
 291     p21.x = p2.x - p1.x;
 292     p21.y = p2.y - p1.y;
 293     p32.x = p3.x - p2.x;
 294     p32.y = p3.y - p2.y;
 295     p43.x = p4.x - p3.x;
 296     p43.y = p4.y - p3.y;
 297
 298     t1 = 1 - t;
 299     t1_2 = t1 * t1;
 300     t_2 = t * t;
 301
 302     vector->x = 3 * t1_2 * p21.x + 6 * t1 * t * p32.x + 3 * t_2 * p43.x;
 303     vector->y = 3 * t1_2 * p21.y + 6 * t1 * t * p32.y + 3 * t_2 * p43.y;
 304 }
 305
 306
 307 static void
 308 put_extents(const CpmlPrimitive *curve, CpmlExtents *extents)
 309 {
 310     CpmlPair p1, p2, p3, p4;
 311
 312     extents->is_defined = 0;
 313
 314     cpml_primitive_put_point(curve, 0, &p1);
 315     cpml_primitive_put_point(curve, 1, &p2);
 316     cpml_primitive_put_point(curve, 2, &p3);
 317     cpml_primitive_put_point(curve, 3, &p4);
 318
 319     cpml_extents_pair_add(extents, &p1);
 320     cpml_extents_pair_add(extents, &p2);
 321     cpml_extents_pair_add(extents, &p3);
 322     cpml_extents_pair_add(extents, &p4);
 323 }
 324
 325 static void
 326 offset(CpmlPrimitive *curve, double offset)
 327 {
 328     double m, mm;
 329     CpmlVector v0, v3, vm, vtmp;
 330     CpmlPair p0, p1, p2, p3, pm;
 331
 332     m = 0.5;
 333     mm = 1-m;
 334
 335     /* Firstly, convert the curve points from cairo format to cpml format
 336      * and store them (temporary) in p0..p3 */
 337     cpml_pair_from_cairo(&p0, curve->org);
 338     cpml_pair_from_cairo(&p1, &curve->data[1]);
 339     cpml_pair_from_cairo(&p2, &curve->data[2]);
 340     cpml_pair_from_cairo(&p3, &curve->data[3]);
 341
 342     /* v0 = p1-p0 */
 343     v0.x = p1.x - p0.x;
 344     v0.y = p1.y - p0.y;
 345
 346     /* v3 = p3-p2 */
 347     v3.x = p3.x - p2.x;
 348     v3.y = p3.y - p2.y;
 349
 350     /* pm = point in C(m) offseted the requested @offset distance */
 351     cpml_curve_put_vector_at_time(curve, m, &vm);
 352     cpml_vector_set_length(&vm, offset);
 353     cpml_vector_normal(&vm);
 354     cpml_curve_put_pair_at_time(curve, m, &pm);
 355     pm.x += vm.x;
 356     pm.y += vm.y;
 357
 358     /* p0 = p0 + normal of v0 of @offset magnitude (exact value) */
 359     cpml_pair_copy(&vtmp, &v0);
 360     cpml_vector_set_length(&vtmp, offset);
 361     cpml_vector_normal(&vtmp);
 362     p0.x += vtmp.x;
 363     p0.y += vtmp.y;
 364
 365     /* p3 = p3 + normal of v3 of @offset magnitude, as done for p0 */
 366     cpml_pair_copy(&vtmp, &v3);
 367     cpml_vector_set_length(&vtmp, offset);
 368     cpml_vector_normal(&vtmp);
 369     p3.x += vtmp.x;
 370     p3.y += vtmp.y;
 371
 372     if (v0.x*v3.y == v3.x*v0.y) {
 373         /* Inconsistent equations: use the alternative approach */
 374         p1.x = p0.x + v0.x + vm.x * 4/3;
 375         p1.y = p0.y + v0.y + vm.y * 4/3;
 376         p2.x = p3.x - v3.x + vm.x * 4/3;
 377         p2.y = p3.y - v3.y + vm.y * 4/3;
 378     } else {
 379         CpmlPair pk;
 380         double k0, k3;
 381
 382         pk.x = (pm.x - mm*mm*(1+m+m)*p0.x - m*m*(1+mm+mm)*p3.x) / (3*m*(1-m));
 383         pk.y = (pm.y - mm*mm*(1+m+m)*p0.y - m*m*(1+mm+mm)*p3.y) / (3*m*(1-m));
 384
 385         if (v0.x != 0) {
 386             k3 = (pk.y - pk.x*v0.y / v0.x) / (m*(v3.y - v3.x*v0.y / v0.x));
 387             k0 = (pk.x - m*k3*v3.x) / (mm*v0.x);
 388         } else {
 389             k0 = (pk.y - pk.x*v3.y / v3.x) / (mm*(v0.y - v0.x*v3.y / v3.x));
 390             k3 = (pk.x - mm*k0*v0.x) / (m*v3.x);
 391         }
 392
 393         p1.x = p0.x + k0*v0.x;
 394         p1.y = p0.y + k0*v0.y;
 395         p2.x = p3.x + k3*v3.x;
 396         p2.y = p3.y + k3*v3.y;
 397     }
 398
 399     /* Return the new curve in the original array */
 400     cpml_pair_to_cairo(&p0, curve->org);
 401     cpml_pair_to_cairo(&p1, &curve->data[1]);
 402     cpml_pair_to_cairo(&p2, &curve->data[2]);
 403     cpml_pair_to_cairo(&p3, &curve->data[3]);
 404 }