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