曲线拟合问题，图形图像高手请进(200分)

如何用圆弧和直线（圆弧优先）来拟合样条曲线？？请高手赐教！要求拟合的圆弧和直线
和标准的样条曲线在允许的误差范围内

三次样条插值

我不知道你为什么需要用圆弧或直线来拟合样条曲线?
样条曲线用于曲线拟合或插值的方法已经相当成熟,绝对比用圆弧或直线拟合出来的效果要好.
而且你使用的样条是几次的样条呢?二次的或三次的,就有很大的不同.

To JohnsonGuo:
如果只是在电脑屏幕上显示或打印出来，当然不需要用圆弧或直线拟合样条曲线，但我
现在是将电脑所绘图形在机床上加工出来，而机床只能有圆弧或直线插补，所以必须要将
样条曲线转化成圆弧或直线（优先圆弧）才能在数控机床上加工出电脑所绘的零件，至于
样条曲线是是二次还是三次，实际上都需要的，不知大侠有何办法解决此问题，同时也希
望其他高手不吝赐教！！！

算法思路:
假设原曲线没有高频振荡(如sin(1000t))
1. 首先使用直线段进行拟合(使用一定步长)
2. 在直线段拟合的基础上,尝试把若干段直线段组合起来使用圆弧拟合,看一下拟合效果是否
提高,如果是则用圆弧段代替该直线段组,否则保留直线段.

程序如下(有三个例子)
如果发现有哪个函数的效果不如理想,请把函数给我.

uses
Windows, SysUtils;

const
InitLineNum = 20; //直线拟合的初始数量
DeltaLineNum = 10; //直线拟合的递增数量

InitSegNum = 20; //多曲线段拟合的初始数量
DeltaSegNum = 10; //多曲线段拟合的递增数量

MinT = 0; //待拟合曲线的时间变化范围
MaxT = 6.28;

MaxLineEps = 5; //直线拟合平均误差的最大误差限
LineCheckCount = 4; //检测直线拟合效果的插值点数量 + 1
InvLineCheckCount = 1 / LineCheckCount;

ArcCheckCount = 10; //检测圆弧拟合效果的插值点数量 + 1
InvArcCheckCount = 1 / ArcCheckCount;

type
//实点结构
TRealPoint = packed record
x, y: Extended;
end;

//多曲线段拟合类型
TSegmentStyle = (ssLine, ssArc);

//多曲线段拟合结构
TSegment = packed record
case Style: TSegmentStyle of
//直线段记下两端点
ssLine: (Pt1, Pt2: TRealPoint);
//圆弧段记下圆心坐标,圆半径,起始角度和结束角度
ssArc: (Org: TRealPoint; Radius, StartAngle, EndAngle: Extended);
end;

var
Pts: array of TRealPoint; //直线拟合的点坐标
Times: array of Extended; //直线拟合的点对应t值
Segs: array of TSegment; //多曲线段拟合的结构数组

function Curve(t: Extended): TRealPoint;
procedure CalcLinePoints;
procedure CalcArcs;

implementation

uses
Math;

var
Pi: Extended;

//已知曲线函数,输入自变量t,返回质点坐标
//Eg.1. 螺旋线 0 < t < Pi, x(t) = 200 cos(t), y(t) = 100 sin(t);
// Pi < t < 2 Pi, x(t) = 140 cos(t) - 60, y(t) = 200 sin(t)
//Eg.2. 正弦曲线 0 < t < 2 Pi, x(t) = t, y(t) = sin(t)
//Eg.3. 星型线 0 < t < 2 Pi, x(t) = 200 cos(t)^3, y(t) = 200 * sin(t)^3
function Curve(t: Extended): TRealPoint;
var
c, s: Extended;
begin
SinCos(t, s, c);

if t > Pi then begin
Result.x := 140 * c - 60;
Result.y := 200 * s;
end else begin
Result.x := 200 * c;
Result.y := 100 * s;
end;

{Result.x := t * 60;
Result.y := Sin(t) * 200; }

{SinCos(t, s, c);
Result.x := 200 * Sqr(c) * c;
Result.y := 200 * Sqr(s) * s; }
end;

//两点距离函数
function Distance(const Pt1, Pt2: TRealPoint): Extended;
begin
Result := Sqrt(Sqr(Pt1.x - Pt2.x) + Sqr(Pt1.y - Pt2.y));
end;

//计算直线段与其上原曲线段的拟合误差估计
function CalcLineEps(const LastT: Extended; const LastPt: TRealPoint;
const ThisT: Extended; const ThisPt: TRealPoint): Extended;
var
DeltaT, t: Extended;
i: Integer;
Pt, LPt: TRealPoint;
begin
DeltaT := (ThisT - LastT) * InvLineCheckCount;
Result := 0;

//把直线分为 LineCheckCount 段进行分别计算其拟合误差取算术平均
for i := 1 to LineCheckCount - 1 do begin
t := LastT + i * DeltaT;
Pt := Curve(t);

LPt.x := (i * ThisPt.x + (LineCheckCount - i) * LastPt.x) * InvLineCheckCount;
LPt.y:= (i * ThisPt.y + (LineCheckCount - i) * LastPt.y) * InvLineCheckCount;

Result := Result + Distance(Pt, LPt);
end;

Result := Result * InvLineCheckCount;
end;

//直线段拟合
procedure CalcLinePoints;
var
DeltaT, t, LastT, ThisDeltaT: Extended;
CurCount, p: Integer;
Pt, LastPt: TRealPoint;
begin
CurCount := InitLineNum + 1;
SetLength(Pts, CurCount);
SetLength(Times, CurCount);
DeltaT := (MaxT - MinT) / InitLineNum; //直线拟合的初始步长

t := MinT;
Pt := Curve(t);

Pts[0] := Pt;
Times[0] := MinT;
p := 1;

LastT := t;
LastPt := Pt;
repeat
ThisDeltaT := DeltaT;
repeat
t := LastT + ThisDeltaT;
if t > MaxT then t := MaxT;

Pt := Curve(t);
//测试选定的直线段是否达到拟合要求,尚未达到则缩小步长重来
if CalcLineEps(LastT, LastPt, t, Pt) <= MaxLineEps then Break;
ThisDeltaT := ThisDeltaT * 0.5;
until t >= MaxT;

//动态分配存储空间
if p >= CurCount then begin
Inc(CurCount, DeltaLineNum);
SetLength(Pts, CurCount);
SetLength(Times, CurCount);
end;

Pts[p] := Pt;
Times[p] := t;
Inc(p);

LastT := t;
LastPt := Pt;
until t >= MaxT;

//保留实际用到的存储空间
SetLength(Pts, p);
SetLength(Times, p);
end;

//用圆弧去拟合原曲线上的一段,返回圆弧相关参数
procedure CalcArcParam(const t1, t2: Extended; var Org: TRealPoint;
var r, StAngle, EdAngle: Extended);
var
Pt1, Pt2, Pt3: TRealPoint;
divd, a, b: Extended;
begin
Pt1 := Curve(t1);
Pt2 := Curve((t1 + t2) * 0.5);
Pt3 := Curve(t2);

divd := 4 * ((Pt2.x - Pt1.x) * (Pt3.y - Pt1.y)
- (Pt2.y - Pt1.y) * (Pt3.x - Pt1.x));
a := Sqr(Pt2.y) - Sqr(Pt1.y) + Sqr(Pt2.x) - Sqr(Pt1.x);
b := Sqr(Pt3.y) - Sqr(Pt1.y) + Sqr(Pt3.x) - Sqr(Pt1.x);
Org.x := 2 * (a * (Pt3.y - Pt1.y) - b * (Pt2.y - Pt1.y)) / divd;
Org.y := 2 * (b * (Pt2.x - Pt1.x) - a * (Pt3.x - Pt1.x)) / divd;
r := Distance(Pt1, Org);
StAngle := ArcTan2(Pt1.y - Org.y, Pt1.x - Org.x);
EdAngle := ArcTan2(Pt3.y - Org.y, Pt3.x - Org.x);
if Abs(EdAngle - StAngle) > Pi then
EdAngle := EdAngle + 2 * Pi;
end;

//计算圆弧与原曲线的拟合误差估计
function CalcArcEps(const t1, t2: Extended; const Org: TRealPoint;
const r, StAngle, EdAngle: Extended): Extended;
var
DeltaT, t, c, s, Angle: Extended;
i: Integer;
Pt, APt: TRealPoint;
begin
DeltaT := (t2 - t1) * InvArcCheckCount;
Result := 0;

//把圆弧分成 ArcCheckCount 段分别计算误差,并取算术平均
for i := 1 to ArcCheckCount - 1 do begin
t := t1 + i * DeltaT;
Pt := Curve(t);

Angle := (i * EdAngle + (ArcCheckCount - i) * StAngle) * InvArcCheckCount;
SinCos(Angle, s, c);
APt.x := Org.x + r * c;
APt.y := Org.y + r * s;

Result := Result + Distance(Pt, APt);
end;

Result := Result * InvArcCheckCount;
end;

//通过直线段拟合的结果进行多曲线段拟合
procedure CalcArcs;
var
CurCount, p1, p2, PtCount, SegCnt, i: Integer;
r, StAngle, EdAngle, LineEps, ArcEps, LastR, LastStAngle, LastEdAngle: Extended;
OrgPt, LastOrgPt: TRealPoint;

//动态分配多曲线段的存储空间
procedure CheckMemory;
begin
if SegCnt >= CurCount then begin
Inc(CurCount, DeltaSegNum);
SetLength(Segs, CurCount);
end;
end;

//添加一条新直线段
procedure AddLine(Point1, Point2: TRealPoint);
begin
with Segs[SegCnt] do begin
Style := ssLine;
Pt1 := Point1;
Pt2 := Point2;
end;
Inc(SegCnt);
end;

//添加一条新圆弧段
procedure AddArc(OrgPt: TRealPoint; r, StAngle, EdAngle: Extended);
begin
with Segs[SegCnt] do begin
Style := ssArc;
Org := OrgPt;
Radius := r;
StartAngle := StAngle;
EndAngle := EdAngle;
end;
Inc(SegCnt);
end;

begin
CurCount := InitSegNum;
SetLength(Segs, CurCount);

p1 := 0;
p2 := 1;
PtCount := Length(Pts);
SegCnt := 0;

LastOrgPt.x := 0;
LastOrgPt.y := 0;
LastR := 0;
LastStAngle := 0;
LastEdAngle := 0;
while p2 < PtCount do begin
//计算使用直线段拟合的误差
LineEps := 0;
for i := p1 + 1 to p2 do
LineEps := LineEps + CalcLineEps(Times[i - 1], Pts[i - 1],
Times, Pts);
LineEps := LineEps / (p2 - p1);

//计算使用圆弧段拟合的误差
CalcArcParam(Times[p1], Times[p2], OrgPt, r, StAngle, EdAngle);
ArcEps := CalcArcEps(Times[p1], Times[p2], OrgPt, r, StAngle, EdAngle);

//测试使用直线拟合和圆弧拟合的拟合效果哪一个好
if LineEps < ArcEps then begin
CheckMemory;

if p2 - p1 = 1 then begin
AddLine(Pts[p1], Pts[p2]);
p1 := p2;
p2 := p1 + 1;
end else begin
AddArc(LastOrgPt, LastR, LastStAngle, LastEdAngle);
p1 := p2 - 1;
end;
end else if p2 = PtCount - 1 then begin
CheckMemory;

if p2 - p1 = 1 then
AddLine(Pts[p1], Pts[p2])
else
AddArc(OrgPt, r, StAngle, EdAngle);

Break;
end else begin
LastOrgPt := OrgPt;
LastStAngle := StAngle;
LastEdAngle := EdAngle;
LastR := r;
Inc(p2);
end;
end;

//保留实际需要的存储空间
SetLength(Segs, SegCnt);
end;

initialization
Pi := System.Pi;

end.

TO JohnsonGuo:
在用圆弧拟合曲线时，必须要保证所拟合的圆弧全部都是相切的，只有这样才能保证
曲线的光顺性要求，你的程序好象并没有考虑到这一点，而这一点是非常非常重要的！
恳请继续赐教！！！！

敬请高手赐教，如嫌分少，我可以再加！！！

gz

贴子提前，请高手赐教！

请高手赐教啊

数学高手在哪里?该问题非常具有挑战性,分不够,我可以再加!!