src/cpml/cpml-curve.c

   1 /* CPML - Cairo Path Manipulation Library
   2  * Copyright (C) 2008,2009,2010  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, I define some useful variables:
  67  *
  68  * |[
  69  * v0 = unitvector(p[1]-p[0]) * offset;
  70  * v3 = unitvector(p[3]-p[2]) * offset;
  71  * p0 = p[0] + normal v0;
  72  * p3 = p[3] + normal v3.
  73  * ]|
  74  *
  75  * Now I want the curve to have the specified slopes at the start
  76  * and end point. Forcing the same slope at the start point means:
  77  *
  78  * |[
  79  * p1 = p0 + k0 v0.
  80  * ]|
  81  *
  82  * where k0 is an arbitrary factor. Decomposing for x and y components:
  83  *
  84  * |[
  85  * p1.x = p0.x + k0 v0.x;
  86  * p1.y = p0.y + k0 v0.y.
  87  * ]|
  88  *
  89  * Doing the same for the end point gives:
  90  *
  91  * |[
  92  * p2.x = p3.x + k3 v3.x;
  93  * p2.y = p3.y + k3 v3.y.
  94  * ]|
  95  *
  96  * Now I interpolate the curve by forcing it to pass throught pm
  97  * when "time" is m, where 0 < m < 1. The cubic Bézier function is:
  98  *
  99  * |[
 100  * C(t) = (1-t)³p0 + 3t(1-t)²p1 + 3t²(1-t)p2 + t³p3.
 101  * ]|
 102  *
 103  * and forcing t=m and C(t) = pm:
 104  *
 105  * |[
 106  * pm = (1-m)³p0 + 3m(1-m)²p1 + 3m²(1-m)p2 + m³p3.
 107  *
 108  * (1-m) p1 + m p2 = (pm - (1-m)³p0 - m³p3) / (3m (1-m)).
 109  * ]|
 110  *
 111  * So the final system is:
 112  *
 113  * |[
 114  * p1.x = p0.x + k0 v0.x;
 115  * p1.y = p0.y + k0 v0.y;
 116  * p2.x = p3.x + k3 v3.x;
 117  * p2.y = p3.y + k3 v3.y;
 118  * (1-m) p1.x + m p2.x = (pm.x - (1-m)³p0.x - m³p3.x) / (3m (1-m));
 119  * (1-m) p1.y + m p2.y = (pm.y - (1-m)³p0.y - m³p3.y) / (3m (1-m)).
 120  * ]|
 121  *
 122  * Substituting and resolving for k0 and k3:
 123  *
 124  * |[
 125  * (1-m) k0 v0.x + m k3 v3.x =
 126  *     (pm.x - (1-m)³p0.x - m³p3.x) / (3m (1-m)) - (1-m) p0.x - m p3.x;
 127  * (1-m) k0 v0.y + m k3 v3.y =
 128  *     (pm.y - (1-m)³p0.y - m³p3.y) / (3m (1-m)) - (1-m) p0.y - m p3.y.
 129  *
 130  * (1-m) k0 v0.x + m k3 v3.x =
 131  *     (pm.x - (1-m)²(1+2m) p0.x - m²(3-2m) p3.x) / (3m (1-m));
 132  * (1-m) k0 v0.y + m k3 v3.y =
 133  *     (pm.y - (1-m)²(1+2m) p0.y - m²(3-2m) p3.y) / (3m (1-m)).
 134  * ]|
 135  *
 136  * Let:
 137  *
 138  * |[
 139  * pk = (pm - (1-m)²(1+2m) p0 - m²(3-2m) p3) / (3m (1-m)).
 140  * ]|
 141  *
 142  * gives the following system:
 143  *
 144  * |[
 145  * (1-m) k0 v0.x + m k3 v3.x = pk.x;
 146  * (1-m) k0 v0.y + m k3 v3.y = pk.y.
 147  * ]|
 148  *
 149  * Now I should avoid division by 0 troubles. If either v0.x and v3.x
 150  * are 0, the first equation will be inconsistent. More in general the
 151  * v0.x*v3.y = v3.x*v3.y condition should be avoided. This is the first
 152  * case to check, in which case an alternative approach is used. In the
 153  * other cases the above system can be used.
 154  *
 155  * If v0.x != 0 I can resolve for k0 and then find k3:
 156  *
 157  * |[
 158  * k0 = (pk.x - m k3 v3.x) / ((1-m) v0.x);
 159  * (pk.x - m k3 v3.x) v0.y / v0.x + m k3 v3.y = pk.y.
 160  *
 161  * k0 = (pk.x - m k3 v3.x) / ((1-m) v0.x);
 162  * k3 m (v3.y - v3.x v0.y / v0.x) = pk.y - pk.x v0.y / v0.x.
 163  *
 164  * k3 = (pk.y - pk.x v0.y / v0.x) / (m (v3.y - v3.x v0.y / v0.x));
 165  * k0 = (pk.x - m k3 v3.x) / ((1-m) v0.x).
 166  * ]|
 167  *
 168  * If v3.x != 0 I can resolve for k3 and then find k0:
 169  *
 170  * |[
 171  * k3 = (pk.x - (1-m) k0 v0.x) / (m v3.x);
 172  * (1-m) k0 v0.y + (pk.x - (1-m) k0 v0.x) v3.y / v3.x = pk.y.
 173  *
 174  * k3 = (pk.x - (1-m) k0 v0.x) / (m v3.x);
 175  * k0 (1-m) (v0.y - k0 v0.x v3.y / v3.x) = pk.y - pk.x v3.y / v3.x.
 176  *
 177  * k0 = (pk.y - pk.x v3.y / v3.x) / ((1-m) (v0.y - v0.x v3.y / v3.x));
 178  * k3 = (pk.x - (1-m) k0 v0.x) / (m v3.x).
 179  * ]|
 180  * </para>
 181  * </refsect2>
 182  * <para>
 183  **/
 184
 185 /**
 186  * CPML_CURVE:
 187  *
 188  * The type code used to identify "curve-to" primitives.
 189  * It is equivalent to the %CAIRO_PATH_CURVE_TO cairo constant.
 190  **/
 191
 192
 193 #include "cpml-internal.h"
 194 #include "cpml-extents.h"
 195 #include "cpml-segment.h"
 196 #include "cpml-primitive.h"
 197 #include "cpml-primitive-private.h"
 198 #include "cpml-curve.h"
 199
 200
 201 static void     put_extents             (const CpmlPrimitive    *curve,
 202                                          CpmlExtents            *extents);
 203 static void     offset                  (CpmlPrimitive          *curve,
 204                                          double                  offset);
 205
 206
 207 const _CpmlPrimitiveClass *
 208 _cpml_curve_get_class(void)
 209 {
 210     static _CpmlPrimitiveClass *p_class = NULL;
 211
 212     if (p_class == NULL) {
 213         static _CpmlPrimitiveClass class_data = {
 214             "curve to", 4,
 215             NULL,
 216             put_extents,
 217             NULL,
 218             NULL,
 219             NULL,
 220             NULL,
 221             offset,
 222             NULL
 223         };
 224         p_class = &class_data;
 225     }
 226
 227     return p_class;
 228 }
 229
 230
 231 /**
 232  * cpml_curve_put_pair_at_time:
 233  * @curve: the #CpmlPrimitive curve data
 234  * @t:     the "time" value
 235  * @pair:  the destination pair
 236  *
 237  * Given the @curve Bézier cubic, finds the coordinates at time @t
 238  * (where 0 is the start and 1 is the end) and stores the result
 239  * in @pair. Keep in mind @t is not homogeneous, so 0.5 does not
 240  * necessarily means the mid point.
 241  **/
 242 void
 243 cpml_curve_put_pair_at_time(const CpmlPrimitive *curve, double t,
 244                             CpmlPair *pair)
 245 {
 246     cairo_path_data_t *p1, *p2, *p3, *p4;
 247     double t_2, t_3, t1, t1_2, t1_3;
 248
 249     p1 = cpml_primitive_get_point(curve, 0);
 250     p2 = cpml_primitive_get_point(curve, 1);
 251     p3 = cpml_primitive_get_point(curve, 2);
 252     p4 = cpml_primitive_get_point(curve, 3);
 253
 254     t_2 = t * t;
 255     t_3 = t_2 * t;
 256     t1 = 1 - t;
 257     t1_2 = t1 * t1;
 258     t1_3 = t1_2 * t1;
 259
 260     pair->x = t1_3 * p1->point.x + 3 * t1_2 * t * p2->point.x
 261               + 3 * t1 * t_2 * p3->point.x + t_3 * p4->point.x;
 262     pair->y = t1_3 * p1->point.y + 3 * t1_2 * t * p2->point.y
 263               + 3 * t1 * t_2 * p3->point.y + t_3 * p4->point.y;
 264 }
 265
 266 /**
 267  * cpml_curve_put_vector_at_time:
 268  * @curve:  the #CpmlPrimitive curve data
 269  * @t:      the "time" value
 270  * @vector: the destination vector
 271  *
 272  * Given the @curve Bézier cubic, finds the slope at time @t
 273  * (where 0 is the start and 1 is the end) and stores the result
 274  * in @vector. Keep in mind @t is not homogeneous, so 0.5
 275  * does not necessarily means the mid point.
 276  **/
 277 void
 278 cpml_curve_put_vector_at_time(const CpmlPrimitive *curve,
 279                               double t, CpmlVector *vector)
 280 {
 281     cairo_path_data_t *p1, *p2, *p3, *p4;
 282     CpmlPair p21, p32, p43;
 283     double t1, t1_2, t_2;
 284
 285     p1 = cpml_primitive_get_point(curve, 0);
 286     p2 = cpml_primitive_get_point(curve, 1);
 287     p3 = cpml_primitive_get_point(curve, 2);
 288     p4 = cpml_primitive_get_point(curve, 3);
 289
 290     p21.x = p2->point.x - p1->point.x;
 291     p21.y = p2->point.y - p1->point.y;
 292     p32.x = p3->point.x - p2->point.x;
 293     p32.y = p3->point.y - p2->point.y;
 294     p43.x = p4->point.x - p3->point.x;
 295     p43.y = p4->point.y - p3->point.y;
 296
 297     t1 = 1 - t;
 298     t1_2 = t1 * t1;
 299     t_2 = t * t;
 300
 301     vector->x = 3 * t1_2 * p21.x + 6 * t1 * t * p32.x + 3 * t_2 * p43.x;
 302     vector->y = 3 * t1_2 * p21.y + 6 * t1 * t * p32.y + 3 * t_2 * p43.y;
 303 }
 304
 305
 306 static void
 307 put_extents(const CpmlPrimitive *curve, CpmlExtents *extents)
 308 {
 309     CpmlPair p1, p2, p3, p4;
 310
 311     extents->is_defined = 0;
 312
 313     cpml_pair_from_cairo(&p1, cpml_primitive_get_point(curve, 0));
 314     cpml_pair_from_cairo(&p2, cpml_primitive_get_point(curve, 1));
 315     cpml_pair_from_cairo(&p3, cpml_primitive_get_point(curve, 2));
 316     cpml_pair_from_cairo(&p4, cpml_primitive_get_point(curve, 3));
 317
 318     cpml_extents_pair_add(extents, &p1);
 319     cpml_extents_pair_add(extents, &p2);
 320     cpml_extents_pair_add(extents, &p3);
 321     cpml_extents_pair_add(extents, &p4);
 322 }
 323
 324 static void
 325 offset(CpmlPrimitive *curve, double offset)
 326 {
 327     double m, mm;
 328     CpmlVector v0, v3, vm, vtmp;
 329     CpmlPair p0, p1, p2, p3, pm;
 330
 331     m = 0.5;
 332     mm = 1-m;
 333
 334     /* Firstly, convert the curve points from cairo format to cpml format
 335      * and store them (temporary) in p0..p3 */
 336     cpml_pair_from_cairo(&p0, curve->org);
 337     cpml_pair_from_cairo(&p1, &curve->data[1]);
 338     cpml_pair_from_cairo(&p2, &curve->data[2]);
 339     cpml_pair_from_cairo(&p3, &curve->data[3]);
 340
 341     /* v0 = p1-p0 */
 342     v0.x = p1.x - p0.x;
 343     v0.y = p1.y - p0.y;
 344
 345     /* v3 = p3-p2 */
 346     v3.x = p3.x - p2.x;
 347     v3.y = p3.y - p2.y;
 348
 349     /* pm = point in C(m) offseted the requested @offset distance */
 350     cpml_curve_put_vector_at_time(curve, m, &vm);
 351     cpml_vector_set_length(&vm, offset);
 352     cpml_vector_normal(&vm);
 353     cpml_curve_put_pair_at_time(curve, m, &pm);
 354     pm.x += vm.x;
 355     pm.y += vm.y;
 356
 357     /* p0 = p0 + normal of v0 of @offset magnitude (exact value) */
 358     cpml_pair_copy(&vtmp, &v0);
 359     cpml_vector_set_length(&vtmp, offset);
 360     cpml_vector_normal(&vtmp);
 361     p0.x += vtmp.x;
 362     p0.y += vtmp.y;
 363
 364     /* p3 = p3 + normal of v3 of @offset magnitude, as done for p0 */
 365     cpml_pair_copy(&vtmp, &v3);
 366     cpml_vector_set_length(&vtmp, offset);
 367     cpml_vector_normal(&vtmp);
 368     p3.x += vtmp.x;
 369     p3.y += vtmp.y;
 370
 371     if (v0.x*v3.y == v3.x*v0.y) {
 372         /* Inconsistent equations: use the alternative approach */
 373         p1.x = p0.x + v0.x + vm.x * 4/3;
 374         p1.y = p0.y + v0.y + vm.y * 4/3;
 375         p2.x = p3.x - v3.x + vm.x * 4/3;
 376         p2.y = p3.y - v3.y + vm.y * 4/3;
 377     } else {
 378         CpmlPair pk;
 379         double k0, k3;
 380
 381         pk.x = (pm.x - mm*mm*(1+m+m)*p0.x - m*m*(1+mm+mm)*p3.x) / (3*m*(1-m));
 382         pk.y = (pm.y - mm*mm*(1+m+m)*p0.y - m*m*(1+mm+mm)*p3.y) / (3*m*(1-m));
 383
 384         if (v0.x != 0) {
 385             k3 = (pk.y - pk.x*v0.y / v0.x) / (m*(v3.y - v3.x*v0.y / v0.x));
 386             k0 = (pk.x - m*k3*v3.x) / (mm*v0.x);
 387         } else {
 388             k0 = (pk.y - pk.x*v3.y / v3.x) / (mm*(v0.y - v0.x*v3.y / v3.x));
 389             k3 = (pk.x - mm*k0*v0.x) / (m*v3.x);
 390         }
 391
 392         p1.x = p0.x + k0*v0.x;
 393         p1.y = p0.y + k0*v0.y;
 394         p2.x = p3.x + k3*v3.x;
 395         p2.y = p3.y + k3*v3.y;
 396     }
 397
 398     /* Return the new curve in the original array */
 399     cpml_pair_to_cairo(&p0, curve->org);
 400     cpml_pair_to_cairo(&p1, &curve->data[1]);
 401     cpml_pair_to_cairo(&p2, &curve->data[2]);
 402     cpml_pair_to_cairo(&p3, &curve->data[3]);
 403 }