插值算法 (100分)

欲寻找Delphi下线性插值和抛物线插值算法。

这本书里有，
书名：《常用数值算法丛书——Delphi常用数值算法集》
作者：何光渝等
ISBN号：703009699
出版社：科学出版社
出版日期：2001年9月

我用c写过,现在不知能不能找到
找一本数值分析的书看一下就可以搞定了

我这里有一个。
二次线性插值。没有处理四边，但是没有关系，
经我测试，放大后，显示效果很好。
function TForm1.MyInsertValue(lWidth, lHeight: integer;
pRGB: pRGBColor;
LineNum, ZoomScale: Single;
Bitmap: TBitmap): pRGBColor;
var
i1, i2, j1, j2, j, i: integer;
f1, f2, f3, f4, f12, f34, fNew: byte;
x, y: single;
pRGBTemp, pRGBNext: pRGBColor;
begin
//行数
x := (LineNum / ZoomScale);
i1 := trunc(x);
i2 := i1 + 1;
if i2>lWidth - 1 then
i2 := lWidth - 1;
Bitmap.PixelFormat := pf24Bit;
pRGBTemp := Bitmap.ScanLine[i1];
pRGBNext := Bitmap.ScanLine[i2];
//列数
for j := 0 to trunc(lWidth * ZoomScale) - 1do
begin
y := (j / ZoomScale);
j1 := trunc;
j2 := j1 + 1;
if j2>lHeight - 1 then
j2 := lHeight - 1;
if (x=0) or (x=(lWidth*ZoomScale-1)/ZoomScale) or (y=0) or (y=(lWidth*ZoomScale-1)/ZoomScale) then
begin
pRGB[j].rgbtBlue := 0;
pRGB[j].rgbtGreen := 0;
pRGB[j].rgbtRed := 0;
end else
begin
{if (i1 = x) then
begin
if (j1 = y) then
begin
f1 := pRGBTemp^[j1].rgbtBlue;
fNew := f1;
end else
begin
f1 := pRGBTemp^[j1].rgbtBlue;
f2 := pRGBTemp^[j2].rgbtBlue;
fNew := trunc(f1 + (y - j1)*(f2 - f1));
end;
end else
begin
}
f1 := pRGBTemp^[j1].rgbtBlue;
f2 := pRGBTemp^[j2].rgbtBlue;
f3 := pRGBNext^[j1].rgbtBlue;
f4 := pRGBNext^[j2].rgbtBlue;
f12 := trunc(f1 + (y - j1)*(f2 - f1));
f34 := trunc(f3 + (y - j1)*(f4 - f3));
if f12 > 255 then
f12 := 255;
if f34 > 255 then
f34 := 255;
fNew := trunc(f12 + (x - i1) * (f34 - f12));
if fNew > 255 then
fNew := 255;
//end;

pRGB[j].rgbtBlue := fNew;
pRGB[j].rgbtGreen := fNew;
pRGB[j].rgbtRed := fNew;
end;
end;
MyInsertValue := pRGB;
end;

pRGBColor的定义是不是：
TRGBColor = record R, G, B: Byte;
end;
PRGBColor = ^TRGBColor;
我运行了，出错呀[]

1) SPLINE interpolation, construct a cubic spline
procedure spline(x, y: glnarray;
n: Integer;
yp1, ypn: Real;
var y2: glnarray);
(* Programs using routine SPLINE must define the type
TYPE
glnarray = ARRAY [1..n] OF real;
in the main routine. *)
var
i, k: Integer;
p, qn, sig, un: Real;
u: glnarray;
begin
if (yp1 > 0.99e30) then
begin
y2[1] := 0.0;
u[1] := 0.0
end
else
begin
y2[1] := -0.5;
u[1] := (3.0 / (x[2] - x[1])) *
((y[2] - y[1]) / (x[2] - x[1]) - yp1)
end;
for i := 2 to n - 1do
begin
sig := (x - x[i - 1]) / (x[i + 1] - x[i - 1]);
p := sig * y2[i - 1] + 2.0;
y2 := (sig - 1.0) / p;
u := (y[i + 1] - y) / (x[i + 1] - x) -
(y - y[i - 1]) / (x - x[i - 1]);
u := (6.0 * u / (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p
end;
if (ypn > 0.99e30) then
begin
qn := 0.0;
un := 0.0
end
else
begin
qn := 0.5;
un := (3.0 / (x[n] - x[n - 1])) *
(ypn - (y[n] - y[n - 1]) / (x[n] - x[n - 1]))
end;
y2[n] := (un - qn * u[n - 1]) / (qn * y2[n - 1] + 1.0);
for k := n - 1do
wnto 1do
begin
y2[k] := y2[k] * y2[k + 1] + u[k]
end
end;

2)SPLINT interpolation, evaluate a cubic spline
procedure splint(xa, ya, y2a: glnarray;
n: Integer;
x: Real;
var y: Real);
(* Programs using routine SPLINT must define the type
TYPE
glnarray = ARRAY [1..n] OF real;
in the main routine. *)
var
klo, khi, k: Integer;
h, b, a: Real;
begin
klo := 1;
khi := n;
while (khi - klo > 1)do
begin
k := (khi + klo) div 2;
if (xa[k] > x) then
khi := k
else
klo := k
end;
h := xa[khi] - xa[klo];
if (h = 0.0) then
begin
writeln('pause in routine SPLINT');
writeln(' ... bad XA input');
readln
end;
a := (xa[khi] - x) / h;
b := (x - xa[klo]) / h;
y := a * ya[klo] + b * ya[khi] + ((a * a * a - a) * y2a[klo] +
(b * b * b - b) * y2a[khi]) * (h * h) / 6.0
end;

谢谢