cpml/cpml-curve.c

   1 /* CPML - Cairo Path Manipulation Library
   2  * Copyright (C) 2008, 2009  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  * SECTION:curve
  22  * @title: Bézier curves
  23  * @short_description: Bézier cubic curve primitive management
  24  *
  25  * The following functions manipulate %CAIRO_PATH_CURVE_TO #CpmlPrimitive.
  26  * No check is made on the primitive struct, so be sure
  27  * <structname>CpmlPrimitive</structname> is effectively a Bézier curve
  28  * before calling these APIs.
  29  **/
  30
  31 #include "cpml-curve.h"
  32 #include "cpml-pair.h"
  33
  34 /**
  35  * cpml_curve_type_get_npoints:
  36  *
  37  * Returns the number of point needed to properly specify a curve primitive.
  38  *
  39  * Return value: 4
  40  **/
  41 int
  42 cpml_curve_type_get_npoints(void)
  43 {
  44     return 4;
  45 }
  46
  47 /**
  48  * cpml_curve_pair_at_time:
  49  * @curve: the #CpmlPrimitive curve data
  50  * @pair:  the destination pair
  51  * @t:     the "time" value
  52  *
  53  * Given the @curve Bézier cubic, finds the coordinates at time @t
  54  * (where 0 is the start and 1 is the end) and stores the result
  55  * in @pair. Keep in mind @t is not homogeneous, so 0.5 does not
  56  * necessarily means the mid point.
  57  *
  58  * The relation 0 < @t < 1 must be satisfied, as interpolating on
  59  * cubic curves is not allowed.
  60  **/
  61 void
  62 cpml_curve_pair_at_time(const CpmlPrimitive *curve, CpmlPair *pair, double t)
  63 {
  64     cairo_path_data_t *p1, *p2, *p3, *p4;
  65     double t_2, t_3, t1, t1_2, t1_3;
  66
  67     p1 = cpml_primitive_get_point(curve, 0);
  68     p2 = cpml_primitive_get_point(curve, 1);
  69     p3 = cpml_primitive_get_point(curve, 2);
  70     p4 = cpml_primitive_get_point(curve, 3);
  71
  72     t_2 = t * t;
  73     t_3 = t_2 * t;
  74     t1 = 1 - t;
  75     t1_2 = t1 * t1;
  76     t1_3 = t1_2 * t1;
  77
  78     pair->x = t1_3 * p1->point.x + 3 * t1_2 * t * p2->point.x
  79               + 3 * t1 * t_2 * p3->point.x + t_3 * p4->point.x;
  80     pair->y = t1_3 * p1->point.y + 3 * t1_2 * t * p2->point.y
  81               + 3 * t1 * t_2 * p3->point.y + t_3 * p4->point.y;
  82 }
  83
  84 /**
  85  * cpml_curve_pair_at:
  86  * @curve: the #CpmlPrimitive curve data
  87  * @pair:  the destination #AdgPair
  88  * @pos:   the position value
  89  *
  90  * Given the @curve Bézier cubic, finds the coordinates at position
  91  * @pos (where 0 is the start and 1 is the end) and stores the result
  92  * in @pair. It is similar to cpml_curve_pair_at_time() but the @pos
  93  * value is evenly distribuited, that is 0.5 is exactly the mid point.
  94  * If you do not need this feature, use cpml_curve_pair_at_time()
  95  * as it is considerable faster.
  96  *
  97  * The relation 0 < @pos < 1 must be satisfied, as interpolating on
  98  * cubic curves is not allowed.
  99  *
 100  * <important>
 101  * <title>TODO</title>
 102  * <itemizedlist>
 103  * <listitem>To be implemented...</listitem>
 104  * </itemizedlist>
 105  * </important>
 106  **/
 107 void
 108 cpml_curve_pair_at(const CpmlPrimitive *curve, CpmlPair *pair, double pos)
 109 {
 110 }
 111
 112 /**
 113  * cpml_curve_vector_at_time:
 114  * @curve:  the #CpmlPrimitive curve data
 115  * @vector: the destination vector
 116  * @t:      the "time" value
 117  *
 118  * Given the @curve Bézier cubic, finds the slope at time @t
 119  * (where 0 is the start and 1 is the end) and stores the result
 120  * in @vector. Keep in mind @t is not homogeneous, so 0.5
 121  * does not necessarily means the mid point.
 122  *
 123  * @t must be inside the range 0 .. 1, as interpolating is not
 124  * allowed.
 125  **/
 126 void
 127 cpml_curve_vector_at_time(const CpmlPrimitive *curve,
 128                           CpmlVector *vector, double t)
 129 {
 130     cairo_path_data_t *p1, *p2, *p3, *p4;
 131     CpmlPair p21, p32, p43;
 132     double t1, t1_2, t_2;
 133
 134     p1 = cpml_primitive_get_point(curve, 0);
 135     p2 = cpml_primitive_get_point(curve, 1);
 136     p3 = cpml_primitive_get_point(curve, 2);
 137     p4 = cpml_primitive_get_point(curve, 3);
 138
 139     p21.x = p2->point.x - p1->point.x;
 140     p21.y = p2->point.y - p1->point.y;
 141     p32.x = p3->point.x - p2->point.x;
 142     p32.y = p3->point.y - p2->point.y;
 143     p43.x = p4->point.x - p3->point.x;
 144     p43.y = p4->point.y - p3->point.y;
 145
 146     t1 = 1 - t;
 147     t1_2 = t1 * t1;
 148     t_2 = t * t;
 149
 150     vector->x = 3 * t1_2 * p21.x + 6 * t1 * t * p32.x + 3 * t_2 * p43.x;
 151     vector->y = 3 * t1_2 * p21.y + 6 * t1 * t * p32.y + 3 * t_2 * p43.y;
 152 }
 153
 154 /**
 155  * cpml_curve_vector_at:
 156  * @curve:  the #CpmlPrimitive curve data
 157  * @vector: the destination vector
 158  * @pos:    the position value
 159  *
 160  * Given the @curve Bézier cubic, finds the slope at position @pos
 161  * (where 0 is the start and 1 is the end) and stores the result
 162  * in @vector. It is similar to cpml_curve_vector_at_time() but the
 163  * @pos value is evenly distribuited, that is 0.5 is exactly the
 164  * mid point. If you do not need this feature, use
 165  * cpml_curve_vector_at_time() as it is considerable faster.
 166  *
 167  * @pos must be inside the range 0 .. 1, as interpolating is not
 168  * allowed.
 169  *
 170  * <important>
 171  * <title>TODO</title>
 172  * <itemizedlist>
 173  * <listitem>To be implemented...</listitem>
 174  * </itemizedlist>
 175  * </important>
 176  **/
 177 void
 178 cpml_curve_vector_at(const CpmlPrimitive *curve,
 179                      CpmlVector *vector, double pos)
 180 {
 181 }
 182
 183 /**
 184  * cpml_curve_intersection:
 185  * @curve:  the first curve
 186  * @curve2: the second curve
 187  * @dest:   a vector of at least 4 #CpmlPair
 188  *
 189  * Given two Bézier cubic curves (@curve and @curve2), gets their
 190  * intersection points and store the result in @dest. Because two
 191  * curves can have 4 intersections, @dest MUST be at least an array
 192  * of 4 #CpmlPair.
 193  *
 194  * <important>
 195  * <title>TODO</title>
 196  * <itemizedlist>
 197  * <listitem>To be implemented...</listitem>
 198  * </itemizedlist>
 199  * </important>
 200  *
 201  * Return value: the number of intersections (max 4)
 202  **/
 203 int
 204 cpml_curve_intersection(const CpmlPrimitive *curve,
 205                         const CpmlPrimitive *curve2, CpmlPair *dest)
 206 {
 207     return 0;
 208 }
 209
 210 /**
 211  * cpml_curve_intersection_with_arc:
 212  * @curve: a curve
 213  * @arc:   an arc
 214  * @dest:  a vector of at least 4 #CpmlPair
 215  *
 216  * Given a Bézier cubic @curve and an @arc, gets their intersection
 217  * points and store the result in @dest. Because an arc and a cubic
 218  * curve can have up to 4 intersections, @dest MUST be at least an
 219  * array of 4 #CpmlPair.
 220  *
 221  * <important>
 222  * <title>TODO</title>
 223  * <itemizedlist>
 224  * <listitem>To be implemented...</listitem>
 225  * </itemizedlist>
 226  * </important>
 227  *
 228  * Return value: the number of intersections (max 4)
 229  **/
 230 int
 231 cpml_curve_intersection_with_arc(const CpmlPrimitive *curve,
 232                                  const CpmlPrimitive *arc, CpmlPair *dest)
 233 {
 234     return 0;
 235 }
 236
 237 /**
 238  * cpml_curve_intersection_with_line:
 239  * @curve: a curve
 240  * @line:  a line
 241  * @dest:  a vector of at least 4 #CpmlPair
 242  *
 243  * Given a Bézier cubic @curve and a @line, gets their intersection
 244  * points and store the result in @dest. Because a line and a cubic
 245  * curve can have up to 4 intersections, @dest MUST be at least an
 246  * array of 4 #CpmlPair.
 247  *
 248  * <important>
 249  * <title>TODO</title>
 250  * <itemizedlist>
 251  * <listitem>To be implemented...</listitem>
 252  * </itemizedlist>
 253  * </important>
 254  *
 255  * Return value: the number of intersections (max 4)
 256  **/
 257 int
 258 cpml_curve_intersection_with_line(const CpmlPrimitive *curve,
 259                                   const CpmlPrimitive *line, CpmlPair *dest)
 260 {
 261     return 0;
 262 }
 263
 264 /**
 265  * cpml_curve_offset:
 266  * @curve:  the #CpmlPrimitive curve data
 267  * @offset: distance for the computed parallel curve
 268  *
 269  * Given a cubic Bézier primitive in @curve, this function finds
 270  * the approximated Bézier curve parallel to @curve at distance
 271  * @offset (an offset curve). The four points needed to build the
 272  * new curve are returned in the @curve struct.
 273  *
 274  * To solve the offset problem, a custom algorithm is used. First, the
 275  * resulting curve MUST have the same slope at the start and end point.
 276  * These constraints are not sufficient to resolve the system, so I let
 277  * the curve pass thought a given point (pm, known and got from the
 278  * original curve) at a given time (m, now hardcoded to 0.5).
 279  *
 280  * Firstly, I define some useful variables:
 281  *
 282  * v0 = unitvector(p[1]-p[0]) * offset;
 283  * v3 = unitvector(p[3]-p[2]) * offset;
 284  * p0 = p[0] + normal v0;
 285  * p3 = p[3] + normal v3.
 286  *
 287  * Now I want the curve to have the specified slopes at the start
 288  * and end point. Forcing the same slope at the start point means:
 289  *
 290  * p1 = p0 + k0 v0.
 291  *
 292  * where k0 is an arbitrary factor. Decomposing for x and y components:
 293  *
 294  * p1.x = p0.x + k0 v0.x;
 295  * p1.y = p0.y + k0 v0.y.
 296  *
 297  * Doing the same for the end point gives:
 298  *
 299  * p2.x = p3.x + k3 v3.x;
 300  * p2.y = p3.y + k3 v3.y.
 301  *
 302  * Now I interpolate the curve by forcing it to pass throught pm
 303  * when "time" is m, where 0 < m < 1. The cubic Bézier function is:
 304  *
 305  * C(t) = (1-t)³p0 + 3t(1-t)²p1 + 3t²(1-t)p2 + t³p3.
 306  *
 307  * and forcing t=m and C(t) = pm:
 308  *
 309  * pm = (1-m)³p0 + 3m(1-m)²p1 + 3m²(1-m)p2 + m³p3.
 310  *
 311  * (1-m) p1 + m p2 = (pm - (1-m)³p0 - m³p3) / (3m (1-m)).
 312  *
 313  * So the final system is:
 314  *
 315  * p1.x = p0.x + k0 v0.x;
 316  * p1.y = p0.y + k0 v0.y;
 317  * p2.x = p3.x + k3 v3.x;
 318  * p2.y = p3.y + k3 v3.y;
 319  * (1-m) p1.x + m p2.x = (pm.x - (1-m)³p0.x - m³p3.x) / (3m (1-m));
 320  * (1-m) p1.y + m p2.y = (pm.y - (1-m)³p0.y - m³p3.y) / (3m (1-m)).
 321  *
 322  * Substituting and resolving for k0 and k3:
 323  *
 324  * (1-m) k0 v0.x + m k3 v3.x =
 325  *     (pm.x - (1-m)³p0.x - m³p3.x) / (3m (1-m)) - (1-m) p0.x - m p3.x;
 326  * (1-m) k0 v0.y + m k3 v3.y =
 327  *     (pm.y - (1-m)³p0.y - m³p3.y) / (3m (1-m)) - (1-m) p0.y - m p3.y.
 328  *
 329  * (1-m) k0 v0.x + m k3 v3.x =
 330  *     (pm.x - (1-m)²(1+2m) p0.x - m²(3-2m) p3.x) / (3m (1-m));
 331  * (1-m) k0 v0.y + m k3 v3.y =
 332  *     (pm.y - (1-m)²(1+2m) p0.y - m²(3-2m) p3.y) / (3m (1-m)).
 333  *
 334  * Let:
 335  *
 336  * pk = (pm - (1-m)²(1+2m) p0 - m²(3-2m) p3) / (3m (1-m)).
 337  *
 338  * gives the following system:
 339  *
 340  * (1-m) k0 v0.x + m k3 v3.x = pk.x;
 341  * (1-m) k0 v0.y + m k3 v3.y = pk.y.
 342  *
 343  * Now I should avoid division by 0 troubles. If either v0.x and v3.x
 344  * are 0, the first equation will be inconsistent. More in general the
 345  * v0.x*v3.y = v3.x*v3.y condition should be avoided. This is the first
 346  * case to check, in which case an alternative approach is used. In the
 347  * other cases the above system can be used.
 348  *
 349  * If v0.x != 0 I can resolve for k0 and then find k3:
 350  *
 351  * k0 = (pk.x - m k3 v3.x) / ((1-m) v0.x);
 352  * (pk.x - m k3 v3.x) v0.y / v0.x + m k3 v3.y = pk.y.
 353  *
 354  * k0 = (pk.x - m k3 v3.x) / ((1-m) v0.x);
 355  * k3 m (v3.y - v3.x v0.y / v0.x) = pk.y - pk.x v0.y / v0.x.
 356  *
 357  * k3 = (pk.y - pk.x v0.y / v0.x) / (m (v3.y - v3.x v0.y / v0.x));
 358  * k0 = (pk.x - m k3 v3.x) / ((1-m) v0.x).
 359  *
 360  * If v3.x != 0 I can resolve for k3 and then find k0:
 361  *
 362  * k3 = (pk.x - (1-m) k0 v0.x) / (m v3.x);
 363  * (1-m) k0 v0.y + (pk.x - (1-m) k0 v0.x) v3.y / v3.x = pk.y.
 364  *
 365  * k3 = (pk.x - (1-m) k0 v0.x) / (m v3.x);
 366  * k0 (1-m) (v0.y - k0 v0.x v3.y / v3.x) = pk.y - pk.x v3.y / v3.x.
 367  *
 368  * k0 = (pk.y - pk.x v3.y / v3.x) / ((1-m) (v0.y - v0.x v3.y / v3.x));
 369  * k3 = (pk.x - (1-m) k0 v0.x) / (m v3.x).
 370  *
 371  * <important>
 372  * <title>TODO</title>
 373  * <itemizedlist>
 374  * <listitem>By default, interpolation of the new curve is made by offseting
 375  *           the mid point: use a better candidate.</listitem>
 376  * <listitem>When the equations are inconsistent, the alternative approach
 377  *           performs very bad if <varname>v0</varname> and
 378  *           <varname>v3</varname> are opposite or staggered.</listitem>
 379  * </itemizedlist>
 380  * </important>
 381  **/
 382 void
 383 cpml_curve_offset(CpmlPrimitive *curve, double offset)
 384 {
 385     double m, mm;
 386     CpmlVector v0, v3, vm, vtmp;
 387     CpmlPair p0, p1, p2, p3, pm;
 388
 389     m = 0.5;
 390     mm = 1-m;
 391
 392     /* Firstly, convert the curve points from cairo format to cpml format
 393      * and store them (temporary) in p0..p3 */
 394     cpml_pair_from_cairo(&p0, curve->org);
 395     cpml_pair_from_cairo(&p1, &curve->data[1]);
 396     cpml_pair_from_cairo(&p2, &curve->data[2]);
 397     cpml_pair_from_cairo(&p3, &curve->data[3]);
 398
 399     /* v0 = p1-p0 */
 400     cpml_pair_sub(cpml_pair_copy(&v0, &p1), &p0);
 401
 402     /* v3 = p3-p2 */
 403     cpml_pair_sub(cpml_pair_copy(&v3, &p3), &p2);
 404
 405     /* pm = point in C(m) offseted the requested @offset distance */
 406     cpml_curve_vector_at_time(curve, &vm, m);
 407     cpml_vector_set_length(&vm, offset);
 408     cpml_vector_normal(&vm);
 409     cpml_curve_pair_at_time(curve, &pm, m);
 410     cpml_pair_add(&pm, &vm);
 411
 412     /* p0 = p0 + normal of v0 of @offset magnitude (exact value) */
 413     cpml_vector_set_length(cpml_pair_copy(&vtmp, &v0), offset);
 414     cpml_vector_normal(&vtmp);
 415     cpml_pair_add(&p0, &vtmp);
 416
 417     /* p3 = p3 + normal of v3 of @offset magnitude, as done for p0 */
 418     cpml_vector_set_length(cpml_pair_copy(&vtmp, &v3), offset);
 419     cpml_vector_normal(&vtmp);
 420     cpml_pair_add(&p3, &vtmp);
 421
 422     if (v0.x*v3.y == v3.x*v0.y) {
 423         /* Inconsistent equations: use the alternative approach */
 424         p1.x = p0.x + v0.x + vm.x * 4/3;
 425         p1.y = p0.y + v0.y + vm.y * 4/3;
 426         p2.x = p3.x - v3.x + vm.x * 4/3;
 427         p2.y = p3.y - v3.y + vm.y * 4/3;
 428     } else {
 429         CpmlPair pk;
 430         double k0, k3;
 431
 432         pk.x = (pm.x - mm*mm*(1+m+m)*p0.x - m*m*(1+mm+mm)*p3.x) / (3*m*(1-m));
 433         pk.y = (pm.y - mm*mm*(1+m+m)*p0.y - m*m*(1+mm+mm)*p3.y) / (3*m*(1-m));
 434
 435         if (v0.x != 0) {
 436             k3 = (pk.y - pk.x*v0.y / v0.x) / (m*(v3.y - v3.x*v0.y / v0.x));
 437             k0 = (pk.x - m*k3*v3.x) / (mm*v0.x);
 438         } else {
 439             k0 = (pk.y - pk.x*v3.y / v3.x) / (mm*(v0.y - v0.x*v3.y / v3.x));
 440             k3 = (pk.x - mm*k0*v0.x) / (m*v3.x);
 441         }
 442
 443         p1.x = p0.x + k0*v0.x;
 444         p1.y = p0.y + k0*v0.y;
 445         p2.x = p3.x + k3*v3.x;
 446         p2.y = p3.y + k3*v3.y;
 447     }
 448
 449     /* Return the new curve in the original array */
 450     cpml_pair_to_cairo(&p0, curve->org);
 451     cpml_pair_to_cairo(&p1, &curve->data[1]);
 452     cpml_pair_to_cairo(&p2, &curve->data[2]);
 453     cpml_pair_to_cairo(&p3, &curve->data[3]);
 454 }