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