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