1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
|
unit unit_interpolate;
{Copyright (C) 2024 by Han Kleijn, www.hnsky.org
email: han.k.. at...hnsky.org
This Source Code Form is subject to the terms of the Mozilla Public
License, v. 2.0. If a copy of the MPL was not distributed with this
file, You can obtain one at https://mozilla.org/MPL/2.0/. }
interface
uses
Classes, SysUtils, math,astap_main;
type
pixel = array[0..2] of double;
function bilinear_interpolation(const img :image_array; x,y:double; out colour: pixel): boolean; {calculate image pixel value on subpixel level}
function BicubicInterpolate(const Img: Image_Array; X, Y: Double; out colour: pixel): boolean;
implementation
//the following function is not used since calculating the rectangle areas covered by an imaginary square pixel works a little better. See function hfd()
function bilinear_interpolation(const img :image_array; x,y:double; out colour: pixel): boolean; //calculate image pixel value on subpixel level
var
X0, Y0, X1, Y1,width5,height5,col,nrcolours: Integer;
Q11, Q12, Q21, Q22: Double;
R1, R2: Double;
begin
result:=true;
width5:=Length(img[0,0]);{width}
height5:=Length(img[0]); {height}
nrcolours:=Length(img);
X0 := Floor(X);
Y0 := Floor(Y);
X1 := X0 + 1;
Y1 := Y0 + 1;
if ((X0 <=0) or (X1 >= width5-1) or (Y0 <= 0) or (Y1 >= height5-1)) then
begin
for col:=0 to nrcolours do colour[col]:=0;
result:=false;
exit;
end;
for col:=0 to nrcolours-1 do
begin
Q11 := img[col,y0,x0];
Q12 := img[col,y0,x1];
Q21 := img[col,y1,x0];
Q22 := img[col,y1,x1];
colour[col]:= (Y1 - Y) * (X1 - X) * Q11 + (Y1 - Y) * (X - X0) * Q21 + (Y - Y0) * (X1 - X) * Q12 + (Y - Y0) * (X - X0) * Q22;
end;
end;
function CatmullRom(p0, p1, p2, p3, t: Double): Double;
// A Catmull-Rom spline is a type of cubic Hermite spline that passes through a set of control points.
// It is C1 continuous, meaning that it has continuous first derivatives.
// The spline is not C2 continuous. The second derivative is linearly interpolated within each segment, causing the curvature to vary linearly over the length of the segment.
// The factor 0.5 is called tention and typically provides a good balance between smoothness and sharpness.
//
//
// [ -1 3 -3 1 ] | P0 |
// [ 2 -5 4 -1 ] | P1 |
// P(t):= 0.5*[t^3, t^2, t, 1] * [ -1 0 1 0 ] * | P2 |
// [ 0 2 0 0 ] | P3 |
//
// This becomes:
//
// P(t) = 0.5 * ((-p0 + 3*p1 -3*p2 + p3)*t*t*t
// + (2*p0 -5*p1 + 4*p2 - p3)*t*t
// + (-p0+p2)*t
// + 2*p1)
//
// The interpolated value represents a point on the curve defined by the four control points p0,p1,p2 and p3.
// This curve segment is influenced by the surrounding points, providing a smooth transition that maintains the shape defined by these points.
//
// The function takes four control points: p0, p1, p2, and p3.
// p1 and p2 are the primary points between which the interpolation occurs.
// p0 and p3 influence the shape of the curve but are not directly interpolated between.
// The output of the Catmull-rom is the estimated value between p1 and p2.
//
// t is the interpolation parameter that determines the position along the curve between two control points p1 and p2.
// When t=0, the interpolation result corresponds exactly to the position of the first control point p1.
// At t=0.5, the point is halfway along the curve between p1 and p2.
// When t=1, the result corresponds exactly to the position of the second control point p2.
var
t2, t3: Double;
a0, a1, a2, a3: Double;
begin
t2 := t * t;
t3 := t2 * t;
// Coefficients for the cubic polynomial with tension parameter
a0 := 0.5 * (-p0 + 3.0 * p1 - 3.0 * p2 + p3);
a1 := 0.5 * (2.0 * p0 - 5.0 * p1 + 4.0 * p2 - p3);
a2 := 0.5 * (-p0 + p2);
a3 := 0.5 * 2.0 * p1;
Result := a0 * t3 + a1 * t2 + a2 * t + a3;
end;
function BicubicInterpolate(const Img: Image_Array; X, Y: Double; out colour: pixel): boolean;
var
i, j, width5,height5,nrcolours,col: Integer;
Row: array[0..3] of Double;
IntX, IntY: Integer;
DX, DY: Double;
begin
result:=true;
width5:=Length(img[0,0]);{width}
height5:=Length(img[0]); {height}
nrcolours:=Length(img);
// Determine integer coordinates
IntX := Floor(X);
IntY := Floor(Y);
// Calculate fractional part of the coordinates
DX := X - IntX;
DY := Y - IntY;
// Boundary checks (ensure we have a 4x4 grid to interpolate within)
if (IntX <= 1) or (IntX >= width5 - 2) or (IntY <= 1) or (IntY >=height5 - 2) then
begin
for i:=0 to 2 do colour[i]:=0;
Exit(false);
end;
for col:=0 to nrcolours-1 do
begin
// Interpolate along the rows
for i := -1 to 2 do
Row[i + 1] := CatmullRom( Img[col,IntY + i,IntX - 1],
Img[col,IntY + i,IntX],
Img[col,IntY + i,(IntX + 1)],
Img[col,IntY + i,(IntX + 2)],
DX );
// Interpolate the result along the column
colour[col] := CatmullRom(Row[0], Row[1], Row[2], Row[3], DY);
end;
end;
//test data
{ img[0,0,0] := 5; img[0,1,0] := 6; img[0,2,0] := 7; img[0,3,0] := 8; img[0,4,0] := 9;
img[0,0,1] := 10; img[0,1,1] := 11; img[0,2,1] := 12; img[0,3,1] := 13; img[0,4,1] := 14;
img[0,0,2] := 15; img[0,1,2] := 16; img[0,2,2] := 17; img[0,3,2] := 18; img[0,4,2] := 19;
img[0,0,3] := 20; img[0,1,3] := 21; img[0,2,3] := 22; img[0,3,3] := 23; img[0,4,3] := 24;
img[0,0,4] := 25; img[0,1,4] := 26; img[0,2,4] := 27; img[0,3,4] := 28; img[0,4,4] := 29;
img[0,0,4] := 30; img[0,1,4] := 31; img[0,2,4] := 32; img[0,3,4] := 33; img[0,4,4] := 34;
val:=value_bilinear(2.5,2.5);//20
val:=BicubicInterpolate(2.5,2.5);
}
begin
end.
|
