vbInside - 2D-Interpolation auf 4 Arten

Interpolation ist die Kunst, Zwischenwerte zu berechnen, um z.B. ein 10x10 Array zu einem 100x100 Array zu vergrößern und dabei möglichst die Orignaldaten zu rekonstruieren (angenommen das 100x100 Array wurde vorher auf 10x10 verkleinert). Vor allem bei dem Vergrößern von Bildern sind diese Verfahren interessant, da sie die Bilder unterschiedlich scharf machen.

Der Einfachheit halber rechnen wir nur mit Werten von 0 bis 500 statt im RGB-Modus mit 3x 0 bis 255, die Werte werden nachher durch ein Array in eine Farbe umgewandelt, ähnlich einem Verlauf.

Dieser Code stellt die Verfahren Nearest Neighbour, Bilinear Interpolation, Bicubic Interpolation und Sin(x)/x Interpolation vor (letzteres Verfahren beruht auf der sog. Sinc-Kurve).

Die verschiedenen Verfahren unterscheiden sich in der Anzahl der in die Berechnung einbezogenen Punkten und in ihren Formeln. Am komplexesten ist wohl die Sin(x)/x Interpolation, die für die Berechnung eines Zwischenwertes alles Werte benutzt. Am häufigsten findet man jedoch die bikubische Interpolation, da diese einen angemessen Kompromiss zwischen Geschwindigkeit und Qualität darstellt.

History
27.12.2002 Online gestellt

Autor: Dominik Auras <Dominik_auf_vbinside.de>

Code aus Form1.frm Option Explicit
Dim ColScale() As Long

Private Sub Command1_Click()
Dim a(0 To 64, 0 To 64) As Long, x As Long, y As Long
Dim dw As Long, dh As Long, Raster As Long, b(0 To 2, 0 To 2) As Long
Dim fracx As Double, fracy As Double, nx As Long, ny As Long
Dim kx As Long, ky As Long, rows() As Double
Dim offx As Long, offy As Long, hyp As Double, gk As Double, _
alpha As Double
Dim ak As Double, hyp2 As Double, hyp3 As Double, k() As Long
Dim weight As Double, tw As Double, wc As Double

If Option1(0).Value = True Then
ReDim k(0, 0)
'nearest neighbour
offx = 0: offy = 0
ElseIf Option1(1).Value = True Then
ReDim k(1, 1)
'bilinear
offx = 0: offy = 0
ElseIf Option1(2).Value = True Then
ReDim k(3, 3)
'bicubic
offx = -Int(UBound(k, 1) / 2): offy = -Int(UBound(k, 2) / 2)
ElseIf Option1(3).Value = True Then
ReDim k(3, 3)
'sin(x)/x
offx = -Int(UBound(k, 1) / 2): offy = -Int(UBound(k, 2) / 2)
End If

dw = Int(UBound(a, 1) / UBound(b, 1))
dh = Int(UBound(a, 2) / UBound(b, 2))

ReDim rows(UBound(k, 2))

Raster = 6

Randomize Timer
For x = LBound(b, 1) To UBound(b, 1)
For y = LBound(b, 2) To UBound(b, 2)
b(x, y) = Int(Rnd * 500)
Next y
Next x

For x = LBound(a, 1) To UBound(a, 1)
For y = LBound(a, 2) To UBound(a, 2)
'If x = UBound(a, 1) Then Stop

fracx = x / dw
fracy = y / dh

While fracx > 1: fracx = fracx - 1: Wend
While fracy > 1: fracy = fracy - 1: Wend

'Debug.Print fracx, fracy

For ny = 0 To UBound(k, 2)
For nx = 0 To UBound(k, 1)
kx = Int(x / dw) + nx + offx
ky = Int(y / dh) + ny + offy

If x Mod dw = 0 Then
kx = Int((x - 1) / dw) + nx + offx
End If
If y Mod dh = 0 Then
ky = Int((y - 1) / dh) + ny + offy
End If

If kx < 0 Then
kx = 0
End If
If ky < 0 Then
ky = 0
End If

If kx > UBound(b, 1) Then
kx = UBound(b, 1)
End If

If ky > UBound(b, 2) Then
ky = UBound(b, 2)
End If

k(nx, ny) = b(kx, ky)
Next nx
Next ny

If Option1(0).Value = True Then
a(x, y) = k(0, 0)
'Nearest Neighbour
ElseIf Option1(1).Value = True Then
'Bilinear
For ny = 0 To UBound(rows)
rows(ny) = Bilinear(1 - fracx, k(0, ny), fracx, _
k(1, ny))
Next ny

a(x, y) = Bilinear(1 - fracy, rows(0), fracy, rows(1))
ElseIf Option1(2).Value = True Then
'Bicubic
For ny = 0 To UBound(rows)
rows(ny) = Bicubic(fracx, k(0, ny), k(1, ny), k(2, _
ny), k(3, ny))
Next ny

a(x, y) = Bicubic(fracy, rows(0), rows(1), rows(2), _
rows(3))
ElseIf Option1(3).Value = True Then
'Versuch: sin(x)/x
a(x, y) = 0
tw = 0
For ny = 0 To UBound(k, 2)
For nx = 0 To UBound(k, 1)
hyp3 = Sqr((nx - (-offx + fracx)) ^ 2 + (ny - _
(-offy + fracy)) ^ 2)
'vom subpixel zum kernelpixel
'hyp = Sqr((0.5 - fracx) ^ 2 + (0.5 - fracy) ^ _
2) 'vom mittelpunkt zum subpixel
'hyp2 = Sqr((UBound(k, 1) / 2 - nx) ^ 2 + _
(UBound(k, 2) / 2 - ny) ^ 2) 'vom _
mittelpunkt zum kernelpixel
'gk = UBound(k, 2) / 2 - ny + (fracy - 0.5)
'ak = nx - UBound(k, 1) / 2 + (fracx - 0.5)

'alpha = Atann(ak, gk)
If hyp3 <> 0 Then
weight = Sin(hyp3) / hyp3
Else
weight = 1
End If
'If weight < 0 Or hyp3 < 0 Then Stop

tw = tw + weight

'a(x, y) = a(x, y) + k(nx, ny)
Next nx
Next ny

wc = 1.2 / tw

For ny = 0 To UBound(k, 2)
For nx = 0 To UBound(k, 1)
hyp3 = ((nx - (-offx + fracx)) ^ 2 + (ny - _
(-offy + fracy)) ^ 2) ^ (1 / 2)
'vom subpixel zum kernelpixel
If hyp3 <> 0 Then
weight = Sin(hyp3) / hyp3
Else
weight = 1
End If

a(x, y) = a(x, y) + weight * wc * k(nx, ny)
Next nx
Next ny
End If

If a(x, y) < 0 Then
a(x, y) = 0
End If
If a(x, y) > 500 Then
a(x, y) = 500
End If
Next y
Next x

For y = LBound(a, 2) To UBound(a, 2)
For x = LBound(a, 1) To UBound(a, 1)
'Debug.Print Space(3 - Len(CStr(a(x, y)))) & a(x, y) & " ";
Form1.Line (x + Raster * x, y + Raster * y)-(x + _
Raster * (x + 1), y + Raster * (y + 1)), _
ColScale(a(x, y)), BF
Next x
'Debug.Print
Next y
'Debug.Print
End Sub

Function Bilinear(ByVal d0 As Double, ByVal f0 As Double, _
ByVal d1 As Double, ByVal f1 As Double) As Double
Bilinear = f0 * d0 + f1 * d1
End Function

Function Bicubic(ByVal dist As Double, ByVal f0 As Double, _
ByVal f1 As Double, ByVal f2 As Double, ByVal f3 As _
Double) As Double
Bicubic = (-f0 + f1 - f2 + f3) * dist ^ 3 + (2 * f0 - 2 * _
f1 + f2 - f3) * dist ^ 2 + (-f0 + f2) * dist + f1
End Function

Sub load_color_scale()
Dim fname As String, ff As Long

fname = App.Path + "\color.dat"

If Dir(fname$) = "" Then
Exit Sub
End If

ff = FreeFile
ReDim ColScale(500)
Open fname For Binary As #ff
Get #ff, 1, ColScale()
Close #ff
End Sub

Private Sub Form_Load()
load_color_scale
End Sub

Public Function Atann(x As Double, y As Double)
Dim alpha As Double
Dim atan As Double

Const pi = 3.14159265

If x <> 0 Then
atan = Atn(Abs(y) / Abs(x))
Else
atan = 0
End If

Select Case x
Case Is < 0
If y = 0 Then
alpha = pi
ElseIf y > 0 Then
alpha = pi - atan
ElseIf y < 0 Then
alpha = pi + atan
End If
Case Is > 0
If y = 0 Then
alpha = 0
ElseIf y > 0 Then
alpha = atan
ElseIf y < 0 Then
alpha = -atan
End If
Case Is = 0
Select Case y
Case Is < 0
alpha = pi / 2
' 1st
Case Is > 0
alpha = 3 * pi / 2
' 3rd
Case Is = 0
alpha = 0
End Select
End Select

aus_atann:
Atann = alpha
End Function