karbon

vcurvefit.cpp
1/* This file is part of the KDE project
2 Copyright (C) 2001, 2002, 2003 The Karbon Developers
3
4 This library is free software; you can redistribute it and/or
5 modify it under the terms of the GNU Library 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 Library General Public License for more details.
13
14 You should have received a copy of the GNU Library General Public License
15 along with this library; see the file COPYING.LIB. If not, write to
16 the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
17 * Boston, MA 02110-1301, USA.
18*/
19
20#include <karbon_part.h>
21#include <karbon_view.h>
22#include <core/vcolor.h>
23#include <core/vcomposite.h>
24#include <core/vfill.h>
25#include <core/vstroke.h>
26#include <core/vglobal.h>
27#include <render/vpainter.h>
28#include <render/vpainterfactory.h>
29#include <commands/vshapecmd.h>
30
31/*
32 An Algorithm for Automatically Fitting Digitized Curves
33 by Philip J. Schneider
34 from "Graphics Gems", Academic Press, 1990
35
36 http://www.acm.org/pubs/tog/GraphicsGems/gems/FitCurves.c
37 http://www.acm.org/pubs/tog/GraphicsGems/gems/README
38*/
39
40
41#include "vcurvefit.h"
42
43#define MAXPOINTS 1000 /* The most points you can have */
44
45
46class FitVector {
47 public:
48 FitVector(KoPoint &p){
49 m_X=p.x();
50 m_Y=p.y();
51 }
52
53 FitVector(){
54 m_X=0;
55 m_Y=0;
56 }
57
58 FitVector(KoPoint &a,KoPoint &b){
59 m_X=a.x()-b.x();
60 m_Y=a.y()-b.y();
61 }
62
63 void normalize(){
64 double len=length();
65 if(len==0.0f)
66 return;
67 m_X/=len; m_Y/=len;
68 }
69
70 void negate(){
71 m_X = -m_X;
72 m_Y = -m_Y;
73 }
74
75 void scale(double s){
76 double len = length();
77 if(len==0.0f)
78 return;
79 m_X *= s/len;
80 m_Y *= s/len;
81 }
82
83 double dot(FitVector &v){
84 return ((m_X*v.m_X)+(m_Y*v.m_Y));
85 }
86
87 double length(){
88 return (double) sqrt(m_X*m_X+m_Y*m_Y);
89 }
90
91 KoPoint operator+(KoPoint &p){
92 KoPoint b(p.x()+m_X,p.y()+m_Y);
93 return b;
94 }
95
96 public:
97 double m_X,m_Y;
98};
99
100double distance(KoPoint *p1,KoPoint *p2){
101 double dx = (p1->x()-p2->x());
102 double dy = (p1->y()-p2->y());
103 return sqrt( dx*dx + dy*dy );
104}
105
106
107FitVector ComputeLeftTangent(TQPtrList<KoPoint> &points,int end){
108 FitVector tHat1(*points.at(end+1),*points.at(end));
109
110 tHat1.normalize();
111
112 return tHat1;
113}
114
115FitVector ComputeRightTangent(TQPtrList<KoPoint> &points,int end){
116 FitVector tHat1(*points.at(end-1),*points.at(end));
117
118 tHat1.normalize();
119
120 return tHat1;
121}
122
123/*
124 * ChordLengthParameterize :
125 * Assign parameter values to digitized points
126 * using relative distances between points.
127 */
128static double *ChordLengthParameterize(TQPtrList<KoPoint> points,int first,int last)
129{
130 int i;
131 double *u; /* Parameterization */
132
133 u = new double[(last-first+1)];
134
135 u[0] = 0.0;
136 for (i = first+1; i <= last; i++) {
137 u[i-first] = u[i-first-1] +
138 distance(points.at(i), points.at(i-1));
139 }
140
141 for (i = first + 1; i <= last; i++) {
142 u[i-first] = u[i-first] / u[last-first];
143 }
144
145 return(u);
146}
147
148static FitVector VectorAdd(FitVector a,FitVector b)
149{
150 FitVector c;
151 c.m_X = a.m_X + b.m_X; c.m_Y = a.m_Y + b.m_Y;
152 return (c);
153}
154static FitVector VectorScale(FitVector v,double s)
155{
156 FitVector result;
157 result.m_X = v.m_X * s; result.m_Y = v.m_Y * s;
158 return (result);
159}
160
161static FitVector VectorSub(FitVector a,FitVector b)
162{
163 FitVector c;
164 c.m_X = a.m_X - b.m_X; c.m_Y = a.m_Y - b.m_Y;
165 return (c);
166}
167
168static FitVector ComputeCenterTangent(TQPtrList<KoPoint> points,int center)
169{
170 FitVector V1, V2, tHatCenter;
171
172 FitVector cpointb = *points.at(center-1);
173 FitVector cpoint = *points.at(center);
174 FitVector cpointa = *points.at(center+1);
175
176 V1 = VectorSub(cpointb,cpoint);
177 V2 = VectorSub(cpoint,cpointa);
178 tHatCenter.m_X= ((V1.m_X + V2.m_X)/2.0);
179 tHatCenter.m_Y= ((V1.m_Y + V2.m_Y)/2.0);
180 tHatCenter.normalize();
181 return tHatCenter;
182}
183
184/*
185 * B0, B1, B2, B3 :
186 * Bezier multipliers
187 */
188static double B0(double u)
189{
190 double tmp = 1.0 - u;
191 return (tmp * tmp * tmp);
192}
193
194
195static double B1(double u)
196{
197 double tmp = 1.0 - u;
198 return (3 * u * (tmp * tmp));
199}
200
201static double B2(double u)
202{
203 double tmp = 1.0 - u;
204 return (3 * u * u * tmp);
205}
206
207static double B3(double u)
208{
209 return (u * u * u);
210}
211
212/*
213 * GenerateBezier :
214 * Use least-squares method to find Bezier control points for region.
215 *
216 */
217KoPoint* GenerateBezier(TQPtrList<KoPoint> &points, int first, int last, double *uPrime,FitVector tHat1,FitVector tHat2)
218{
219 int i;
220 FitVector A[MAXPOINTS][2]; /* Precomputed rhs for eqn */
221 int nPts; /* Number of pts in sub-curve */
222 double C[2][2]; /* Matrix C */
223 double X[2]; /* Matrix X */
224 double det_C0_C1, /* Determinants of matrices */
225 det_C0_X,
226 det_X_C1;
227 double alpha_l, /* Alpha values, left and right */
228 alpha_r;
229 FitVector tmp; /* Utility variable */
230 KoPoint *curve;
231
232 curve = new KoPoint[4];
233 nPts = last - first + 1;
234
235
236 /* Compute the A's */
237 for (i = 0; i < nPts; i++) {
238 FitVector v1, v2;
239 v1 = tHat1;
240 v2 = tHat2;
241 v1.scale(B1(uPrime[i]));
242 v2.scale(B2(uPrime[i]));
243 A[i][0] = v1;
244 A[i][1] = v2;
245 }
246
247 /* Create the C and X matrices */
248 C[0][0] = 0.0;
249 C[0][1] = 0.0;
250 C[1][0] = 0.0;
251 C[1][1] = 0.0;
252 X[0] = 0.0;
253 X[1] = 0.0;
254
255 for (i = 0; i < nPts; i++) {
256 C[0][0] += (A[i][0]).dot(A[i][0]);
257 C[0][1] += A[i][0].dot(A[i][1]);
258 /* C[1][0] += V2Dot(&A[i][0], &A[i][1]);*/
259 C[1][0] = C[0][1];
260 C[1][1] += A[i][1].dot(A[i][1]);
261
262 FitVector vfirstp1(*points.at(first+i));
263 FitVector vfirst(*points.at(first));
264 FitVector vlast(*points.at(last));
265
266 tmp = VectorSub(vfirstp1,
267 VectorAdd(
268 VectorScale(vfirst, B0(uPrime[i])),
269 VectorAdd(
270 VectorScale(vfirst, B1(uPrime[i])),
271 VectorAdd(
272 VectorScale(vlast, B2(uPrime[i])),
273 VectorScale(vlast, B3(uPrime[i])) ))));
274
275
276 X[0] += A[i][0].dot(tmp);
277 X[1] += A[i][1].dot(tmp);
278 }
279
280 /* Compute the determinants of C and X */
281 det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1];
282 det_C0_X = C[0][0] * X[1] - C[0][1] * X[0];
283 det_X_C1 = X[0] * C[1][1] - X[1] * C[0][1];
284
285 /* Finally, derive alpha values */
286 if (det_C0_C1 == 0.0) {
287 det_C0_C1 = (C[0][0] * C[1][1]) * 10e-12;
288 }
289 alpha_l = det_X_C1 / det_C0_C1;
290 alpha_r = det_C0_X / det_C0_C1;
291
292
293 /* If alpha negative, use the Wu/Barsky heuristic (see text) */
294 /* (if alpha is 0, you get coincident control points that lead to
295 * divide by zero in any subsequent NewtonRaphsonRootFind() call. */
296 if (alpha_l < 1.0e-6 || alpha_r < 1.0e-6) {
297 double dist = distance(points.at(last),points.at(first)) /
298 3.0;
299
300 curve[0] = *points.at(first);
301 curve[3] = *points.at(last);
302
303 tHat1.scale(dist);
304 tHat2.scale(dist);
305
306 curve[1] = tHat1 + curve[0];
307 curve[2] = tHat2 + curve[3];
308 return curve;
309 }
310
311 /* First and last control points of the Bezier curve are */
312 /* positioned exactly at the first and last data points */
313 /* Control points 1 and 2 are positioned an alpha distance out */
314 /* on the tangent vectors, left and right, respectively */
315 curve[0] = *points.at(first);
316 curve[3] = *points.at(last);
317
318 tHat1.scale(alpha_l);
319 tHat2.scale(alpha_r);
320
321 curve[1] = tHat1 + curve[0];
322 curve[2] = tHat2 + curve[3];
323
324 return (curve);
325}
326
327/*
328 * Bezier :
329 * Evaluate a Bezier curve at a particular parameter value
330 *
331 */
332static KoPoint BezierII(int degree,KoPoint *V, double t)
333{
334 int i, j;
335 KoPoint Q; /* Point on curve at parameter t */
336 KoPoint *Vtemp; /* Local copy of control points */
337
338 Vtemp = new KoPoint[degree+1];
339
340 for (i = 0; i <= degree; i++) {
341 Vtemp[i] = V[i];
342 }
343
344 /* Triangle computation */
345 for (i = 1; i <= degree; i++) {
346 for (j = 0; j <= degree-i; j++) {
347 Vtemp[j].setX((1.0 - t) * Vtemp[j].x() + t * Vtemp[j+1].x());
348 Vtemp[j].setY((1.0 - t) * Vtemp[j].y() + t * Vtemp[j+1].y());
349 }
350 }
351
352 Q = Vtemp[0];
353 delete[] Vtemp;
354 return Q;
355}
356
357/*
358 * ComputeMaxError :
359 * Find the maximum squared distance of digitized points
360 * to fitted curve.
361*/
362static double ComputeMaxError(TQPtrList<KoPoint> points,int first,int last,KoPoint *curve,double *u,int *splitPoint)
363{
364 int i;
365 double maxDist; /* Maximum error */
366 double dist; /* Current error */
367 KoPoint P; /* Point on curve */
368 FitVector v; /* Vector from point to curve */
369
370 *splitPoint = (last - first + 1)/2;
371 maxDist = 0.0;
372 for (i = first + 1; i < last; i++) {
373 P = BezierII(3, curve, u[i-first]);
374 v = VectorSub(P, *points.at(i));
375 dist = v.length();
376 if (dist >= maxDist) {
377 maxDist = dist;
378 *splitPoint = i;
379 }
380 }
381 return (maxDist);
382}
383
384
385/*
386 * NewtonRaphsonRootFind :
387 * Use Newton-Raphson iteration to find better root.
388 */
389static double NewtonRaphsonRootFind(KoPoint *Q,KoPoint P,double u)
390{
391 double numerator, denominator;
392 KoPoint Q1[3], Q2[2]; /* Q' and Q'' */
393 KoPoint TQ_u, Q1_u, Q2_u; /*u evaluated at Q, Q', & Q'' */
394 double uPrime; /* Improved u */
395 int i;
396
397 /* Compute Q(u) */
398 TQ_u = BezierII(3,Q, u);
399
400 /* Generate control vertices for Q' */
401 for (i = 0; i <= 2; i++) {
402 Q1[i].setX((Q[i+1].x() - Q[i].x()) * 3.0);
403 Q1[i].setY((Q[i+1].y() - Q[i].y()) * 3.0);
404 }
405
406 /* Generate control vertices for Q'' */
407 for (i = 0; i <= 1; i++) {
408 Q2[i].setX((Q1[i+1].x() - Q1[i].x()) * 2.0);
409 Q2[i].setY((Q1[i+1].y() - Q1[i].y()) * 2.0);
410 }
411
412 /* Compute Q'(u) and Q''(u) */
413 Q1_u = BezierII(2, Q1, u);
414 Q2_u = BezierII(1, Q2, u);
415
416 /* Compute f(u)/f'(u) */
417 numerator = (TQ_u.x() - P.x()) * (Q1_u.x()) + (TQ_u.y() - P.y()) * (Q1_u.y());
418 denominator = (Q1_u.x()) * (Q1_u.x()) + (Q1_u.y()) * (Q1_u.y()) +
419 (TQ_u.x() - P.x()) * (Q2_u.x()) + (TQ_u.y() - P.y()) * (Q2_u.y());
420
421 /* u = u - f(u)/f'(u) */
422 uPrime = u - (numerator/denominator);
423 return (uPrime);
424}
425
426/*
427 * Reparameterize:
428 * Given set of points and their parameterization, try to find
429 * a better parameterization.
430 *
431 */
432static double *Reparameterize(TQPtrList<KoPoint> points,int first,int last,double *u,KoPoint *curve)
433{
434 int nPts = last-first+1;
435 int i;
436 double *uPrime; /* New parameter values */
437
438 uPrime = new double[nPts];
439 for (i = first; i <= last; i++) {
440 uPrime[i-first] = NewtonRaphsonRootFind(curve, *points.at(i), u[i-
441 first]);
442 }
443 return (uPrime);
444}
445
446KoPoint *FitCubic(TQPtrList<KoPoint> &points,int first,int last,FitVector tHat1,FitVector tHat2,float error,int &width){
447 double *u;
448 double *uPrime;
449 double maxError;
450 int splitPoint;
451 int nPts;
452 double iterationError;
453 int maxIterations=4;
454 FitVector tHatCenter;
455 KoPoint *curve;
456 int i;
457
458 width=0;
459
460
461 iterationError=error*error;
462 nPts = last-first+1;
463
464 if(nPts == 2){
465 double dist = distance(points.at(last), points.at(first)) / 3.0;
466
467 curve = new KoPoint[4];
468
469 curve[0] = *points.at(first);
470 curve[3] = *points.at(last);
471
472 tHat1.scale(dist);
473 tHat2.scale(dist);
474
475 curve[1] = tHat1 + curve[0];
476 curve[2] = tHat2 + curve[3];
477
478 width=4;
479 return curve;
480 }
481
482 /* Parameterize points, and attempt to fit curve */
483 u = ChordLengthParameterize(points, first, last);
484 curve = GenerateBezier(points, first, last, u, tHat1, tHat2);
485
486
487 /* Find max deviation of points to fitted curve */
488 maxError = ComputeMaxError(points, first, last, curve, u, &splitPoint);
489 if (maxError < error) {
490 delete[] u;
491 width=4;
492 return curve;
493 }
494
495
496 /* If error not too large, try some reparameterization */
497 /* and iteration */
498 if (maxError < iterationError) {
499 for (i = 0; i < maxIterations; i++) {
500 uPrime = Reparameterize(points, first, last, u, curve);
501 curve = GenerateBezier(points, first, last, uPrime, tHat1, tHat2);
502 maxError = ComputeMaxError(points, first, last,
503 curve, uPrime, &splitPoint);
504 if (maxError < error) {
505 delete[] u;
506 width=4;
507 return curve;
508 }
509 delete[] u;
510 u = uPrime;
511 }
512 }
513
514 /* Fitting failed -- split at max error point and fit recursively */
515 delete[] u;
516 delete[] curve;
517 tHatCenter = ComputeCenterTangent(points, splitPoint);
518
519 int w1,w2;
520 KoPoint *cu1=NULL, *cu2=NULL;
521 cu1 = FitCubic(points, first, splitPoint, tHat1, tHatCenter, error,w1);
522
523 tHatCenter.negate();
524 cu2 = FitCubic(points, splitPoint, last, tHatCenter, tHat2, error,w2);
525
526 KoPoint *newcurve = new KoPoint[w1+w2];
527 for(int i=0;i<w1;i++){
528 newcurve[i]=cu1[i];
529 }
530 for(int i=0;i<w2;i++){
531 newcurve[i+w1]=cu2[i];
532 }
533
534 delete[] cu1;
535 delete[] cu2;
536 width=w1+w2;
537 return newcurve;
538}
539
540
541VPath *bezierFit(TQPtrList<KoPoint> &points,float error){
542 FitVector tHat1, tHat2;
543
544 tHat1 = ComputeLeftTangent(points,0);
545 tHat2 = ComputeRightTangent(points,points.count()-1);
546
547 int width=0;
548 KoPoint *curve;
549 curve = FitCubic(points,0,points.count()-1,tHat1,tHat2,error,width);
550
551 VPath *path = new VPath(NULL);
552
553 if(width>3){
554 path->moveTo(curve[0]);
555 path->curveTo(curve[1],curve[2],curve[3]);
556 for(int i=4;i<width;i+=4){
557 path->curveTo(curve[i+1],curve[i+2],curve[i+3]);
558 }
559 }
560
561
562 delete[] curve;
563 return path;
564}
565
VPath
A composite path consists of one or many subpaths.
Definition vcomposite.h:48