You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
koffice/karbon/tools/vcurvefit.cc

566 lines
13 KiB

/* This file is part of the KDE project
Copyright (C) 2001, 2002, 2003 The Karbon Developers
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Library General Public
License as published by the Free Software Foundation; either
version 2 of the License, or (at your option) any later version.
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Library General Public License for more details.
You should have received a copy of the GNU Library General Public License
along with this library; see the file COPYING.LIB. If not, write to
the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
* Boston, MA 02110-1301, USA.
*/
#include <karbon_part.h>
#include <karbon_view.h>
#include <core/vcolor.h>
#include <core/vcomposite.h>
#include <core/vfill.h>
#include <core/vstroke.h>
#include <core/vglobal.h>
#include <render/vpainter.h>
#include <render/vpainterfactory.h>
#include <commands/vshapecmd.h>
/*
An Algorithm for Automatically Fitting Digitized Curves
by Philip J. Schneider
from "Graphics Gems", Academic Press, 1990
http://www.acm.org/pubs/tog/GraphicsGems/gems/FitCurves.c
http://www.acm.org/pubs/tog/GraphicsGems/gems/README
*/
#include "vcurvefit.h"
#define MAXPOINTS 1000 /* The most points you can have */
class FitVector {
public:
FitVector(KoPoint &p){
m_X=p.x();
m_Y=p.y();
}
FitVector(){
m_X=0;
m_Y=0;
}
FitVector(KoPoint &a,KoPoint &b){
m_X=a.x()-b.x();
m_Y=a.y()-b.y();
}
void normalize(){
double len=length();
if(len==0.0f)
return;
m_X/=len; m_Y/=len;
}
void negate(){
m_X = -m_X;
m_Y = -m_Y;
}
void scale(double s){
double len = length();
if(len==0.0f)
return;
m_X *= s/len;
m_Y *= s/len;
}
double dot(FitVector &v){
return ((m_X*v.m_X)+(m_Y*v.m_Y));
}
double length(){
return (double) sqrt(m_X*m_X+m_Y*m_Y);
}
KoPoint operator+(KoPoint &p){
KoPoint b(p.x()+m_X,p.y()+m_Y);
return b;
}
public:
double m_X,m_Y;
};
double distance(KoPoint *p1,KoPoint *p2){
double dx = (p1->x()-p2->x());
double dy = (p1->y()-p2->y());
return sqrt( dx*dx + dy*dy );
}
FitVector ComputeLeftTangent(TQPtrList<KoPoint> &points,int end){
FitVector tHat1(*points.at(end+1),*points.at(end));
tHat1.normalize();
return tHat1;
}
FitVector ComputeRightTangent(TQPtrList<KoPoint> &points,int end){
FitVector tHat1(*points.at(end-1),*points.at(end));
tHat1.normalize();
return tHat1;
}
/*
* ChordLengthParameterize :
* Assign parameter values to digitized points
* using relative distances between points.
*/
static double *ChordLengthParameterize(TQPtrList<KoPoint> points,int first,int last)
{
int i;
double *u; /* Parameterization */
u = new double[(last-first+1)];
u[0] = 0.0;
for (i = first+1; i <= last; i++) {
u[i-first] = u[i-first-1] +
distance(points.at(i), points.at(i-1));
}
for (i = first + 1; i <= last; i++) {
u[i-first] = u[i-first] / u[last-first];
}
return(u);
}
static FitVector VectorAdd(FitVector a,FitVector b)
{
FitVector c;
c.m_X = a.m_X + b.m_X; c.m_Y = a.m_Y + b.m_Y;
return (c);
}
static FitVector VectorScale(FitVector v,double s)
{
FitVector result;
result.m_X = v.m_X * s; result.m_Y = v.m_Y * s;
return (result);
}
static FitVector VectorSub(FitVector a,FitVector b)
{
FitVector c;
c.m_X = a.m_X - b.m_X; c.m_Y = a.m_Y - b.m_Y;
return (c);
}
static FitVector ComputeCenterTangent(TQPtrList<KoPoint> points,int center)
{
FitVector V1, V2, tHatCenter;
FitVector cpointb = *points.at(center-1);
FitVector cpoint = *points.at(center);
FitVector cpointa = *points.at(center+1);
V1 = VectorSub(cpointb,cpoint);
V2 = VectorSub(cpoint,cpointa);
tHatCenter.m_X= ((V1.m_X + V2.m_X)/2.0);
tHatCenter.m_Y= ((V1.m_Y + V2.m_Y)/2.0);
tHatCenter.normalize();
return tHatCenter;
}
/*
* B0, B1, B2, B3 :
* Bezier multipliers
*/
static double B0(double u)
{
double tmp = 1.0 - u;
return (tmp * tmp * tmp);
}
static double B1(double u)
{
double tmp = 1.0 - u;
return (3 * u * (tmp * tmp));
}
static double B2(double u)
{
double tmp = 1.0 - u;
return (3 * u * u * tmp);
}
static double B3(double u)
{
return (u * u * u);
}
/*
* GenerateBezier :
* Use least-squares method to find Bezier control points for region.
*
*/
KoPoint* GenerateBezier(TQPtrList<KoPoint> &points, int first, int last, double *uPrime,FitVector tHat1,FitVector tHat2)
{
int i;
FitVector A[MAXPOINTS][2]; /* Precomputed rhs for eqn */
int nPts; /* Number of pts in sub-curve */
double C[2][2]; /* Matrix C */
double X[2]; /* Matrix X */
double det_C0_C1, /* Determinants of matrices */
det_C0_X,
det_X_C1;
double alpha_l, /* Alpha values, left and right */
alpha_r;
FitVector tmp; /* Utility variable */
KoPoint *curve;
curve = new KoPoint[4];
nPts = last - first + 1;
/* Compute the A's */
for (i = 0; i < nPts; i++) {
FitVector v1, v2;
v1 = tHat1;
v2 = tHat2;
v1.scale(B1(uPrime[i]));
v2.scale(B2(uPrime[i]));
A[i][0] = v1;
A[i][1] = v2;
}
/* Create the C and X matrices */
C[0][0] = 0.0;
C[0][1] = 0.0;
C[1][0] = 0.0;
C[1][1] = 0.0;
X[0] = 0.0;
X[1] = 0.0;
for (i = 0; i < nPts; i++) {
C[0][0] += (A[i][0]).dot(A[i][0]);
C[0][1] += A[i][0].dot(A[i][1]);
/* C[1][0] += V2Dot(&A[i][0], &A[i][1]);*/
C[1][0] = C[0][1];
C[1][1] += A[i][1].dot(A[i][1]);
FitVector vfirstp1(*points.at(first+i));
FitVector vfirst(*points.at(first));
FitVector vlast(*points.at(last));
tmp = VectorSub(vfirstp1,
VectorAdd(
VectorScale(vfirst, B0(uPrime[i])),
VectorAdd(
VectorScale(vfirst, B1(uPrime[i])),
VectorAdd(
VectorScale(vlast, B2(uPrime[i])),
VectorScale(vlast, B3(uPrime[i])) ))));
X[0] += A[i][0].dot(tmp);
X[1] += A[i][1].dot(tmp);
}
/* Compute the determinants of C and X */
det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1];
det_C0_X = C[0][0] * X[1] - C[0][1] * X[0];
det_X_C1 = X[0] * C[1][1] - X[1] * C[0][1];
/* Finally, derive alpha values */
if (det_C0_C1 == 0.0) {
det_C0_C1 = (C[0][0] * C[1][1]) * 10e-12;
}
alpha_l = det_X_C1 / det_C0_C1;
alpha_r = det_C0_X / det_C0_C1;
/* If alpha negative, use the Wu/Barsky heuristic (see text) */
/* (if alpha is 0, you get coincident control points that lead to
* divide by zero in any subsequent NewtonRaphsonRootFind() call. */
if (alpha_l < 1.0e-6 || alpha_r < 1.0e-6) {
double dist = distance(points.at(last),points.at(first)) /
3.0;
curve[0] = *points.at(first);
curve[3] = *points.at(last);
tHat1.scale(dist);
tHat2.scale(dist);
curve[1] = tHat1 + curve[0];
curve[2] = tHat2 + curve[3];
return curve;
}
/* First and last control points of the Bezier curve are */
/* positioned exactly at the first and last data points */
/* Control points 1 and 2 are positioned an alpha distance out */
/* on the tangent vectors, left and right, respectively */
curve[0] = *points.at(first);
curve[3] = *points.at(last);
tHat1.scale(alpha_l);
tHat2.scale(alpha_r);
curve[1] = tHat1 + curve[0];
curve[2] = tHat2 + curve[3];
return (curve);
}
/*
* Bezier :
* Evaluate a Bezier curve at a particular parameter value
*
*/
static KoPoint BezierII(int degree,KoPoint *V, double t)
{
int i, j;
KoPoint Q; /* Point on curve at parameter t */
KoPoint *Vtemp; /* Local copy of control points */
Vtemp = new KoPoint[degree+1];
for (i = 0; i <= degree; i++) {
Vtemp[i] = V[i];
}
/* Triangle computation */
for (i = 1; i <= degree; i++) {
for (j = 0; j <= degree-i; j++) {
Vtemp[j].setX((1.0 - t) * Vtemp[j].x() + t * Vtemp[j+1].x());
Vtemp[j].setY((1.0 - t) * Vtemp[j].y() + t * Vtemp[j+1].y());
}
}
Q = Vtemp[0];
delete[] Vtemp;
return Q;
}
/*
* ComputeMaxError :
* Find the maximum squared distance of digitized points
* to fitted curve.
*/
static double ComputeMaxError(TQPtrList<KoPoint> points,int first,int last,KoPoint *curve,double *u,int *splitPoint)
{
int i;
double maxDist; /* Maximum error */
double dist; /* Current error */
KoPoint P; /* Point on curve */
FitVector v; /* Vector from point to curve */
*splitPoint = (last - first + 1)/2;
maxDist = 0.0;
for (i = first + 1; i < last; i++) {
P = BezierII(3, curve, u[i-first]);
v = VectorSub(P, *points.at(i));
dist = v.length();
if (dist >= maxDist) {
maxDist = dist;
*splitPoint = i;
}
}
return (maxDist);
}
/*
* NewtonRaphsonRootFind :
* Use Newton-Raphson iteration to find better root.
*/
static double NewtonRaphsonRootFind(KoPoint *Q,KoPoint P,double u)
{
double numerator, denominator;
KoPoint Q1[3], Q2[2]; /* Q' and Q'' */
KoPoint TQ_u, Q1_u, Q2_u; /*u evaluated at Q, Q', & Q'' */
double uPrime; /* Improved u */
int i;
/* Compute Q(u) */
TQ_u = BezierII(3,Q, u);
/* Generate control vertices for Q' */
for (i = 0; i <= 2; i++) {
Q1[i].setX((Q[i+1].x() - Q[i].x()) * 3.0);
Q1[i].setY((Q[i+1].y() - Q[i].y()) * 3.0);
}
/* Generate control vertices for Q'' */
for (i = 0; i <= 1; i++) {
Q2[i].setX((Q1[i+1].x() - Q1[i].x()) * 2.0);
Q2[i].setY((Q1[i+1].y() - Q1[i].y()) * 2.0);
}
/* Compute Q'(u) and Q''(u) */
Q1_u = BezierII(2, Q1, u);
Q2_u = BezierII(1, Q2, u);
/* Compute f(u)/f'(u) */
numerator = (TQ_u.x() - P.x()) * (Q1_u.x()) + (TQ_u.y() - P.y()) * (Q1_u.y());
denominator = (Q1_u.x()) * (Q1_u.x()) + (Q1_u.y()) * (Q1_u.y()) +
(TQ_u.x() - P.x()) * (Q2_u.x()) + (TQ_u.y() - P.y()) * (Q2_u.y());
/* u = u - f(u)/f'(u) */
uPrime = u - (numerator/denominator);
return (uPrime);
}
/*
* Reparameterize:
* Given set of points and their parameterization, try to tqfind
* a better parameterization.
*
*/
static double *Reparameterize(TQPtrList<KoPoint> points,int first,int last,double *u,KoPoint *curve)
{
int nPts = last-first+1;
int i;
double *uPrime; /* New parameter values */
uPrime = new double[nPts];
for (i = first; i <= last; i++) {
uPrime[i-first] = NewtonRaphsonRootFind(curve, *points.at(i), u[i-
first]);
}
return (uPrime);
}
KoPoint *FitCubic(TQPtrList<KoPoint> &points,int first,int last,FitVector tHat1,FitVector tHat2,float error,int &width){
double *u;
double *uPrime;
double maxError;
int splitPoint;
int nPts;
double iterationError;
int maxIterations=4;
FitVector tHatCenter;
KoPoint *curve;
int i;
width=0;
iterationError=error*error;
nPts = last-first+1;
if(nPts == 2){
double dist = distance(points.at(last), points.at(first)) / 3.0;
curve = new KoPoint[4];
curve[0] = *points.at(first);
curve[3] = *points.at(last);
tHat1.scale(dist);
tHat2.scale(dist);
curve[1] = tHat1 + curve[0];
curve[2] = tHat2 + curve[3];
width=4;
return curve;
}
/* Parameterize points, and attempt to fit curve */
u = ChordLengthParameterize(points, first, last);
curve = GenerateBezier(points, first, last, u, tHat1, tHat2);
/* Find max deviation of points to fitted curve */
maxError = ComputeMaxError(points, first, last, curve, u, &splitPoint);
if (maxError < error) {
delete[] u;
width=4;
return curve;
}
/* If error not too large, try some reparameterization */
/* and iteration */
if (maxError < iterationError) {
for (i = 0; i < maxIterations; i++) {
uPrime = Reparameterize(points, first, last, u, curve);
curve = GenerateBezier(points, first, last, uPrime, tHat1, tHat2);
maxError = ComputeMaxError(points, first, last,
curve, uPrime, &splitPoint);
if (maxError < error) {
delete[] u;
width=4;
return curve;
}
delete[] u;
u = uPrime;
}
}
/* Fitting failed -- split at max error point and fit recursively */
delete[] u;
delete[] curve;
tHatCenter = ComputeCenterTangent(points, splitPoint);
int w1,w2;
KoPoint *cu1=NULL, *cu2=NULL;
cu1 = FitCubic(points, first, splitPoint, tHat1, tHatCenter, error,w1);
tHatCenter.negate();
cu2 = FitCubic(points, splitPoint, last, tHatCenter, tHat2, error,w2);
KoPoint *newcurve = new KoPoint[w1+w2];
for(int i=0;i<w1;i++){
newcurve[i]=cu1[i];
}
for(int i=0;i<w2;i++){
newcurve[i+w1]=cu2[i];
}
delete[] cu1;
delete[] cu2;
width=w1+w2;
return newcurve;
}
VPath *bezierFit(TQPtrList<KoPoint> &points,float error){
FitVector tHat1, tHat2;
tHat1 = ComputeLeftTangent(points,0);
tHat2 = ComputeRightTangent(points,points.count()-1);
int width=0;
KoPoint *curve;
curve = FitCubic(points,0,points.count()-1,tHat1,tHat2,error,width);
VPath *path = new VPath(NULL);
if(width>3){
path->moveTo(curve[0]);
path->curveTo(curve[1],curve[2],curve[3]);
for(int i=4;i<width;i+=4){
path->curveTo(curve[i+1],curve[i+2],curve[i+3]);
}
}
delete[] curve;
return path;
}