00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include <karbon_part.h>
00021 #include <karbon_view.h>
00022 #include <core/vcolor.h>
00023 #include <core/vcomposite.h>
00024 #include <core/vfill.h>
00025 #include <core/vstroke.h>
00026 #include <core/vglobal.h>
00027 #include <render/vpainter.h>
00028 #include <render/vpainterfactory.h>
00029 #include <commands/vshapecmd.h>
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041 #include "vcurvefit.h"
00042
00043 #define MAXPOINTS 1000
00044
00045
00046 class FitVector {
00047 public:
00048 FitVector(KoPoint &p){
00049 m_X=p.x();
00050 m_Y=p.y();
00051 }
00052
00053 FitVector(){
00054 m_X=0;
00055 m_Y=0;
00056 }
00057
00058 FitVector(KoPoint &a,KoPoint &b){
00059 m_X=a.x()-b.x();
00060 m_Y=a.y()-b.y();
00061 }
00062
00063 void normalize(){
00064 double len=length();
00065 if(len==0.0f)
00066 return;
00067 m_X/=len; m_Y/=len;
00068 }
00069
00070 void negate(){
00071 m_X = -m_X;
00072 m_Y = -m_Y;
00073 }
00074
00075 void scale(double s){
00076 double len = length();
00077 if(len==0.0f)
00078 return;
00079 m_X *= s/len;
00080 m_Y *= s/len;
00081 }
00082
00083 double dot(FitVector &v){
00084 return ((m_X*v.m_X)+(m_Y*v.m_Y));
00085 }
00086
00087 double length(){
00088 return (double) sqrt(m_X*m_X+m_Y*m_Y);
00089 }
00090
00091 KoPoint operator+(KoPoint &p){
00092 KoPoint b(p.x()+m_X,p.y()+m_Y);
00093 return b;
00094 }
00095
00096 public:
00097 double m_X,m_Y;
00098 };
00099
00100 double distance(KoPoint *p1,KoPoint *p2){
00101 double dx = (p1->x()-p2->x());
00102 double dy = (p1->y()-p2->y());
00103 return sqrt( dx*dx + dy*dy );
00104 }
00105
00106
00107 FitVector ComputeLeftTangent(QPtrList<KoPoint> &points,int end){
00108 FitVector tHat1(*points.at(end+1),*points.at(end));
00109
00110 tHat1.normalize();
00111
00112 return tHat1;
00113 }
00114
00115 FitVector ComputeRightTangent(QPtrList<KoPoint> &points,int end){
00116 FitVector tHat1(*points.at(end-1),*points.at(end));
00117
00118 tHat1.normalize();
00119
00120 return tHat1;
00121 }
00122
00123
00124
00125
00126
00127
00128 static double *ChordLengthParameterize(QPtrList<KoPoint> points,int first,int last)
00129 {
00130 int i;
00131 double *u;
00132
00133 u = new double[(last-first+1)];
00134
00135 u[0] = 0.0;
00136 for (i = first+1; i <= last; i++) {
00137 u[i-first] = u[i-first-1] +
00138 distance(points.at(i), points.at(i-1));
00139 }
00140
00141 for (i = first + 1; i <= last; i++) {
00142 u[i-first] = u[i-first] / u[last-first];
00143 }
00144
00145 return(u);
00146 }
00147
00148 static FitVector VectorAdd(FitVector a,FitVector b)
00149 {
00150 FitVector c;
00151 c.m_X = a.m_X + b.m_X; c.m_Y = a.m_Y + b.m_Y;
00152 return (c);
00153 }
00154 static FitVector VectorScale(FitVector v,double s)
00155 {
00156 FitVector result;
00157 result.m_X = v.m_X * s; result.m_Y = v.m_Y * s;
00158 return (result);
00159 }
00160
00161 static FitVector VectorSub(FitVector a,FitVector b)
00162 {
00163 FitVector c;
00164 c.m_X = a.m_X - b.m_X; c.m_Y = a.m_Y - b.m_Y;
00165 return (c);
00166 }
00167
00168 static FitVector ComputeCenterTangent(QPtrList<KoPoint> points,int center)
00169 {
00170 FitVector V1, V2, tHatCenter;
00171
00172 FitVector cpointb = *points.at(center-1);
00173 FitVector cpoint = *points.at(center);
00174 FitVector cpointa = *points.at(center+1);
00175
00176 V1 = VectorSub(cpointb,cpoint);
00177 V2 = VectorSub(cpoint,cpointa);
00178 tHatCenter.m_X= ((V1.m_X + V2.m_X)/2.0);
00179 tHatCenter.m_Y= ((V1.m_Y + V2.m_Y)/2.0);
00180 tHatCenter.normalize();
00181 return tHatCenter;
00182 }
00183
00184
00185
00186
00187
00188 static double B0(double u)
00189 {
00190 double tmp = 1.0 - u;
00191 return (tmp * tmp * tmp);
00192 }
00193
00194
00195 static double B1(double u)
00196 {
00197 double tmp = 1.0 - u;
00198 return (3 * u * (tmp * tmp));
00199 }
00200
00201 static double B2(double u)
00202 {
00203 double tmp = 1.0 - u;
00204 return (3 * u * u * tmp);
00205 }
00206
00207 static double B3(double u)
00208 {
00209 return (u * u * u);
00210 }
00211
00212
00213
00214
00215
00216
00217 KoPoint* GenerateBezier(QPtrList<KoPoint> &points, int first, int last, double *uPrime,FitVector tHat1,FitVector tHat2)
00218 {
00219 int i;
00220 FitVector A[MAXPOINTS][2];
00221 int nPts;
00222 double C[2][2];
00223 double X[2];
00224 double det_C0_C1,
00225 det_C0_X,
00226 det_X_C1;
00227 double alpha_l,
00228 alpha_r;
00229 FitVector tmp;
00230 KoPoint *curve;
00231
00232 curve = new KoPoint[4];
00233 nPts = last - first + 1;
00234
00235
00236
00237 for (i = 0; i < nPts; i++) {
00238 FitVector v1, v2;
00239 v1 = tHat1;
00240 v2 = tHat2;
00241 v1.scale(B1(uPrime[i]));
00242 v2.scale(B2(uPrime[i]));
00243 A[i][0] = v1;
00244 A[i][1] = v2;
00245 }
00246
00247
00248 C[0][0] = 0.0;
00249 C[0][1] = 0.0;
00250 C[1][0] = 0.0;
00251 C[1][1] = 0.0;
00252 X[0] = 0.0;
00253 X[1] = 0.0;
00254
00255 for (i = 0; i < nPts; i++) {
00256 C[0][0] += (A[i][0]).dot(A[i][0]);
00257 C[0][1] += A[i][0].dot(A[i][1]);
00258
00259 C[1][0] = C[0][1];
00260 C[1][1] += A[i][1].dot(A[i][1]);
00261
00262 FitVector vfirstp1(*points.at(first+i));
00263 FitVector vfirst(*points.at(first));
00264 FitVector vlast(*points.at(last));
00265
00266 tmp = VectorSub(vfirstp1,
00267 VectorAdd(
00268 VectorScale(vfirst, B0(uPrime[i])),
00269 VectorAdd(
00270 VectorScale(vfirst, B1(uPrime[i])),
00271 VectorAdd(
00272 VectorScale(vlast, B2(uPrime[i])),
00273 VectorScale(vlast, B3(uPrime[i])) ))));
00274
00275
00276 X[0] += A[i][0].dot(tmp);
00277 X[1] += A[i][1].dot(tmp);
00278 }
00279
00280
00281 det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1];
00282 det_C0_X = C[0][0] * X[1] - C[0][1] * X[0];
00283 det_X_C1 = X[0] * C[1][1] - X[1] * C[0][1];
00284
00285
00286 if (det_C0_C1 == 0.0) {
00287 det_C0_C1 = (C[0][0] * C[1][1]) * 10e-12;
00288 }
00289 alpha_l = det_X_C1 / det_C0_C1;
00290 alpha_r = det_C0_X / det_C0_C1;
00291
00292
00293
00294
00295
00296 if (alpha_l < 1.0e-6 || alpha_r < 1.0e-6) {
00297 double dist = distance(points.at(last),points.at(first)) /
00298 3.0;
00299
00300 curve[0] = *points.at(first);
00301 curve[3] = *points.at(last);
00302
00303 tHat1.scale(dist);
00304 tHat2.scale(dist);
00305
00306 curve[1] = tHat1 + curve[0];
00307 curve[2] = tHat2 + curve[3];
00308 return curve;
00309 }
00310
00311
00312
00313
00314
00315 curve[0] = *points.at(first);
00316 curve[3] = *points.at(last);
00317
00318 tHat1.scale(alpha_l);
00319 tHat2.scale(alpha_r);
00320
00321 curve[1] = tHat1 + curve[0];
00322 curve[2] = tHat2 + curve[3];
00323
00324 return (curve);
00325 }
00326
00327
00328
00329
00330
00331
00332 static KoPoint BezierII(int degree,KoPoint *V, double t)
00333 {
00334 int i, j;
00335 KoPoint Q;
00336 KoPoint *Vtemp;
00337
00338 Vtemp = new KoPoint[degree+1];
00339
00340 for (i = 0; i <= degree; i++) {
00341 Vtemp[i] = V[i];
00342 }
00343
00344
00345 for (i = 1; i <= degree; i++) {
00346 for (j = 0; j <= degree-i; j++) {
00347 Vtemp[j].setX((1.0 - t) * Vtemp[j].x() + t * Vtemp[j+1].x());
00348 Vtemp[j].setY((1.0 - t) * Vtemp[j].y() + t * Vtemp[j+1].y());
00349 }
00350 }
00351
00352 Q = Vtemp[0];
00353 delete[] Vtemp;
00354 return Q;
00355 }
00356
00357
00358
00359
00360
00361
00362 static double ComputeMaxError(QPtrList<KoPoint> points,int first,int last,KoPoint *curve,double *u,int *splitPoint)
00363 {
00364 int i;
00365 double maxDist;
00366 double dist;
00367 KoPoint P;
00368 FitVector v;
00369
00370 *splitPoint = (last - first + 1)/2;
00371 maxDist = 0.0;
00372 for (i = first + 1; i < last; i++) {
00373 P = BezierII(3, curve, u[i-first]);
00374 v = VectorSub(P, *points.at(i));
00375 dist = v.length();
00376 if (dist >= maxDist) {
00377 maxDist = dist;
00378 *splitPoint = i;
00379 }
00380 }
00381 return (maxDist);
00382 }
00383
00384
00385
00386
00387
00388
00389 static double NewtonRaphsonRootFind(KoPoint *Q,KoPoint P,double u)
00390 {
00391 double numerator, denominator;
00392 KoPoint Q1[3], Q2[2];
00393 KoPoint Q_u, Q1_u, Q2_u;
00394 double uPrime;
00395 int i;
00396
00397
00398 Q_u = BezierII(3,Q, u);
00399
00400
00401 for (i = 0; i <= 2; i++) {
00402 Q1[i].setX((Q[i+1].x() - Q[i].x()) * 3.0);
00403 Q1[i].setY((Q[i+1].y() - Q[i].y()) * 3.0);
00404 }
00405
00406
00407 for (i = 0; i <= 1; i++) {
00408 Q2[i].setX((Q1[i+1].x() - Q1[i].x()) * 2.0);
00409 Q2[i].setY((Q1[i+1].y() - Q1[i].y()) * 2.0);
00410 }
00411
00412
00413 Q1_u = BezierII(2, Q1, u);
00414 Q2_u = BezierII(1, Q2, u);
00415
00416
00417 numerator = (Q_u.x() - P.x()) * (Q1_u.x()) + (Q_u.y() - P.y()) * (Q1_u.y());
00418 denominator = (Q1_u.x()) * (Q1_u.x()) + (Q1_u.y()) * (Q1_u.y()) +
00419 (Q_u.x() - P.x()) * (Q2_u.x()) + (Q_u.y() - P.y()) * (Q2_u.y());
00420
00421
00422 uPrime = u - (numerator/denominator);
00423 return (uPrime);
00424 }
00425
00426
00427
00428
00429
00430
00431
00432 static double *Reparameterize(QPtrList<KoPoint> points,int first,int last,double *u,KoPoint *curve)
00433 {
00434 int nPts = last-first+1;
00435 int i;
00436 double *uPrime;
00437
00438 uPrime = new double[nPts];
00439 for (i = first; i <= last; i++) {
00440 uPrime[i-first] = NewtonRaphsonRootFind(curve, *points.at(i), u[i-
00441 first]);
00442 }
00443 return (uPrime);
00444 }
00445
00446 KoPoint *FitCubic(QPtrList<KoPoint> &points,int first,int last,FitVector tHat1,FitVector tHat2,float error,int &width){
00447 double *u;
00448 double *uPrime;
00449 double maxError;
00450 int splitPoint;
00451 int nPts;
00452 double iterationError;
00453 int maxIterations=4;
00454 FitVector tHatCenter;
00455 KoPoint *curve;
00456 int i;
00457
00458 width=0;
00459
00460
00461 iterationError=error*error;
00462 nPts = last-first+1;
00463
00464 if(nPts == 2){
00465 double dist = distance(points.at(last), points.at(first)) / 3.0;
00466
00467 curve = new KoPoint[4];
00468
00469 curve[0] = *points.at(first);
00470 curve[3] = *points.at(last);
00471
00472 tHat1.scale(dist);
00473 tHat2.scale(dist);
00474
00475 curve[1] = tHat1 + curve[0];
00476 curve[2] = tHat2 + curve[3];
00477
00478 width=4;
00479 return curve;
00480 }
00481
00482
00483 u = ChordLengthParameterize(points, first, last);
00484 curve = GenerateBezier(points, first, last, u, tHat1, tHat2);
00485
00486
00487
00488 maxError = ComputeMaxError(points, first, last, curve, u, &splitPoint);
00489 if (maxError < error) {
00490 delete[] u;
00491 width=4;
00492 return curve;
00493 }
00494
00495
00496
00497
00498 if (maxError < iterationError) {
00499 for (i = 0; i < maxIterations; i++) {
00500 uPrime = Reparameterize(points, first, last, u, curve);
00501 curve = GenerateBezier(points, first, last, uPrime, tHat1, tHat2);
00502 maxError = ComputeMaxError(points, first, last,
00503 curve, uPrime, &splitPoint);
00504 if (maxError < error) {
00505 delete[] u;
00506 width=4;
00507 return curve;
00508 }
00509 delete[] u;
00510 u = uPrime;
00511 }
00512 }
00513
00514
00515 delete[] u;
00516 delete[] curve;
00517 tHatCenter = ComputeCenterTangent(points, splitPoint);
00518
00519 int w1,w2;
00520 KoPoint *cu1=NULL, *cu2=NULL;
00521 cu1 = FitCubic(points, first, splitPoint, tHat1, tHatCenter, error,w1);
00522
00523 tHatCenter.negate();
00524 cu2 = FitCubic(points, splitPoint, last, tHatCenter, tHat2, error,w2);
00525
00526 KoPoint *newcurve = new KoPoint[w1+w2];
00527 for(int i=0;i<w1;i++){
00528 newcurve[i]=cu1[i];
00529 }
00530 for(int i=0;i<w2;i++){
00531 newcurve[i+w1]=cu2[i];
00532 }
00533
00534 delete[] cu1;
00535 delete[] cu2;
00536 width=w1+w2;
00537 return newcurve;
00538 }
00539
00540
00541 VPath *bezierFit(QPtrList<KoPoint> &points,float error){
00542 FitVector tHat1, tHat2;
00543
00544 tHat1 = ComputeLeftTangent(points,0);
00545 tHat2 = ComputeRightTangent(points,points.count()-1);
00546
00547 int width=0;
00548 KoPoint *curve;
00549 curve = FitCubic(points,0,points.count()-1,tHat1,tHat2,error,width);
00550
00551 VPath *path = new VPath(NULL);
00552
00553 if(width>3){
00554 path->moveTo(curve[0]);
00555 path->curveTo(curve[1],curve[2],curve[3]);
00556 for(int i=4;i<width;i+=4){
00557 path->curveTo(curve[i+1],curve[i+2],curve[i+3]);
00558 }
00559 }
00560
00561
00562 delete[] curve;
00563 return path;
00564 }
00565