suche Mathematik-Python-Genie, biquadratische Interpolation

[kussji]

Hallo
Gibt es hier vielleicht jemanden der mir helfen kann, eine mehr oder weniger fertige Lösung zu präsentieren oder Schritt für Schritt zu helfen.

Es geht um eine Biquadratische-Interpolation welche mit Python gelöst werden soll.
Wozu brauche ich das? Es geht um eine GPS-Anwendung welche sehr genaue Resultate liefern soll (Vermessung).

Für eine exakte Höhenmessung gibt es in einem Raster (3x3 km) bekannte Höhenwerte. Aus diesen soll für einen gesuchten Punkt (aktueller GPS-Standort) die gesuchte Höhe durch biquadratische Interpolation der 9 gegebenen 3D-Punkte gefunden werden.

Mir fehlen dazu leider sowohl die mathematischen, als auch die pythonspezifischen Kenntnisse (denke Sympy wäre da mein Kollege).
In einem Mathe-Forum habe hierzu einige Infos erhalten zur mathematischen Lösung. Verstehe das aber nicht ganz und wenn ich es verstehen würde, fehlt mir immer noch die Umsetzung in Python.

Danke für eure Tipps/Hilfen
Grüsse euch
kussji

[sparrow]

Warum musst du das in Python machen? Gibt es denn Programmiersprachen, in denen du sicher stehst?

[kussji]

Hallo
Python behersche ich die Grundlagen nicht schlecht, habe zwar wieder längeren Unterbruch, sollte da aber schnell wieder drin sein.
Die höheren mathematische Funktionen habe ich in sympy mal angesehen. Da wäre bestimmt was drin für mich - nur mir fehlt aber das mathematische Wissen um die richtigen Funktionen zu finden.
Meine Anwendung besteht schon in Python...

Eine Lösung in "Geogebra" wäre ev. auch etwas das mir weiterhelfen würde. Da bin ich aber auch neu...

[Sirius3]

Sind die Stützstellen frei verteilt, dann suchst Du scipy.interpolate.interp2d?

Wenn die Stützstellen auf einem regulären Gitter liegen, dann brauchst Du 4x4 also 16 Werte für die biquadratische Interpolation.
Angenommen, Du hast die Höhenwerte auf einem Gitter `heights` gegeben, wie in einem Bild mit Pixelwerten. Und Du möchtest den Wert eines Zwischenpixels biquadratisch interpolieren, dann ist das einfach:

Code: Alles auswählen

pxi, pyi = int(numpy.floor(px)), int(numpy.floor(py))
dx=px-pxi; dxx=dx*dx; dxxx=dxx*dx
dy=py-pyi; dyy=dy*dy; dyyy=dyy*dy

wx=(0.5 * (    - dx + 2.0*dxx -       dxxx),
    0.5 * (2.0      - 5.0*dxx + 3.0 * dxxx),
    0.5 * (      dx + 4.0*dxx - 3.0 * dxxx),
    0.5 * (         -     dxx +       dxxx))

wy=(0.5 * (    - dy + 2.0*dyy -       dyyy),
    0.5 * (2.0      - 5.0*dyy + 3.0 * dyyy),
    0.5 * (      dy + 4.0*dyy - 3.0 * dyyy),
    0.5 * (         -     dyy +       dyyy))

result = 0
for a, wa in enumerate(wx, pxi-1):
    for b, wb in enumerate(wy, pyi-1):
        result += heights[b, a] * wa * wb

Mehr als Mal und Plus ist da nicht dabei.

[kussji]

Hallo und Danke
Das sieht schon mal nicht so gefürchig aus, wie ich vermutet hatte.
Die Diskusion im Mathe-Forum war von 9 Gleichungen mit 9 Unbekannten, Ableitung etc. - vielleicht viel zu weit gedacht?!

Es ist korrekt, dass das Gitter auf welchem die bekannten zu interpolierenden Korrekturhöhen liegen, ein reguläres Gitter ist. (Je nach System entweder x,y 1x1km-Raster, oder 30 Bogensekunden für Longitude,Latitude)

4x4 Punkte erscheint mir logischer. Von einem profesionellen Onlinerechendienst, welcher ebenfalls die biquadratische Interpolation nutzt, weiss ich, dass er mit den nächsten 9 Punkten (3x3) interpoliert!??? Gibt es hier eine Erklärung?

Die korrekte Verwendung habe ich noch nicht ganz begriffen (habe mal versucht, so weit ich es verstanden habe...)
Das ganze mal in eine Funktion gepackt:

Code: Alles auswählen

import numpy

# Biquadratische Interpolation
def biQuadInt(px, py, heights):
    pxi, pyi = int(numpy.floor(px)), int(numpy.floor(py))
    dx=px-pxi; dxx=dx*dx; dxxx=dxx*dx
    dy=py-pyi; dyy=dy*dy; dyyy=dyy*dy

    wx=(0.5 * (    - dx + 2.0*dxx -       dxxx),
        0.5 * (2.0      - 5.0*dxx + 3.0 * dxxx),
        0.5 * (      dx + 4.0*dxx - 3.0 * dxxx),
        0.5 * (         -     dxx +       dxxx))

    wy=(0.5 * (    - dy + 2.0*dyy -       dyyy),
        0.5 * (2.0      - 5.0*dyy + 3.0 * dyyy),
        0.5 * (      dy + 4.0*dyy - 3.0 * dyyy),
        0.5 * (         -     dyy +       dyyy))

    result = 0
    for a, wa in enumerate(wx, pxi-1):
        for b, wb in enumerate(wy, pyi-1):
            result += heights[b, a] * wa * wb

#------------------------            

#Grundlage:
# x,y-Raster regulär in 1x1 km (x:0-300 km, y:0-200 km)
heights = ??????	# hhmm? wie muss das aussehen? Array[0..15] of Float???
px, py = 258.359, 189.678	#[km] für diesen GPS-Punkt die gesuchte interpolierte Korrekturhöhe
print(biQuadInt(py,py,heights))

[kussji]

... so ist der Aufruf wohl besser - aber auch nicht richtig

Code: Alles auswählen

#Grundlage: x,y-Raster regulär in 1x1 km (x:0-300 km, y:0-200km)
heights = [ [53.2, 53.6, 53.7, 53.5],
            [53.1, 53.3, 53.0, 52.9],
            [52.4, 52.7, 52.1, 52.0],
            [54.6, 54.2, 54.0, 54.5]] # Korrekturhöhe in [m] hhmm? wie muss das aussehen? Array[Float-Werte]???
px, py = 258.3, 189.6 # für diesen GPS-Punkt die gesuchte Interpolierte Höhe
print(biQuadInt(px, py, heights))

[Sirius3]

Wie schon geschrieben, kommt das von der Bildverarbeitung. Pixel-Koordinaten gehen von 0, 1, 2, ...
Du mußt also Deine x-, y-Werte er Rasterpunkte auf diese ganzen Zahlen umrechnen. In Deinem 4x4-Bildchen sind, weil die Randpunkte schlecht quadratisch interpolierbar sind, nur werte zwischen 1 und 2 sinnvoll.

Code: Alles auswählen

heights = numpy.array([ [53.2, 53.6, 53.7, 53.5],
            [53.1, 53.3, 53.0, 52.9],
            [52.4, 52.7, 52.1, 52.0],
            [54.6, 54.2, 54.0, 54.5]])
px, py = 1.3, 1.6 # für diesen Rasterpunkt die gesuchte Interpolierte Höhe
print(biQuadInt(px, py, heights))

Normalerweise hat man ja ein großes Gitter, wo man weiß was für Koordinaten die linke obere oder die rechte untere Ecke haben und kann dann entsprechend in Pixel-Koordinaten umrechnen. Woher bekommst Du Deine Höhen?

Die Gewichtung gilt jetzt nur für den Fall, dass die Erde zwischen den Punkten flach ist, dürfte bei einem 1km-Raster jetzt noch relativ gut passen, zumal Höhen sich innerhalb eines Kilometers deutlich ändern könnten.

[kussji]

Guten Morgen
Danke für die Inputs. Entweder bin ich noch nicht ganz wach oder es liegen andere Probleme vor .

Die Funktion liefert mit Deinem genannten Aufruf als Resultat "None"

Pixel-Koordinaten gehen von 0, 1, 2, ...

...okey - wenn ich aber ein 1km - Raster habe würde das ja passen? Dem entsprechend müssten "heights" auch in [km] stehen. Also:

Code: Alles auswählen

heights = numpy.array([  [0.0532, 0.0536, 0.0537, 0.0538],
                         [0.0531, 0.0533, 0.0530, 0.0529],
                         [0.0524, 0.0527, 0.0521, 0.0520], 
                         [0.0546, 0.0542, 0.0540, 0.0545]])

Woher bekommst Du Deine Höhen?

Die Werte stammen aus einer Tabelle mit xyz-Werten diese beschreibt das Geoid.
Es geht hier nicht um absolute (Meereshöhen/Gebrauchshöhen) sondern um Korrekturwerte (Geoidundulation). Deswegen sind die Änderungen klein (auch in hügeligem/bergigem Gelände). Für meine Fläche von einigen Quadratkilometern ist die maximale Änderung nur 10 Meter. Innerhalb meines gesuchten Punktes gar nur wenige Centimeter)

[kussji]

Entweder bin ich noch nicht ganz wach oder es liegen andere Probleme vor .
Die Funktion liefert mit Deinem genannten Aufruf als Resultat "None"

... ich glaube es lag am noch nicht "wach" bzw. als ich die Funktion geschrieben habe, wohl zu spät in der Nacht...
Ich hatte das "return result" am Ende der Funktion vergessen.

Hier eine Zusammenfassung meiner offenen Fragen:

1. Zitat Sirius3: "Pixel-Koordinaten gehen von 0, 1, 2, ..."
   • ...okey - wenn ich aber ein 1km - Raster habe würde das ja passen?
2. Zitat kussji: "Dem entsprechend müssten "heights" auch in [km] stehen, statt [m]?"
   • denke inzwischen, dem ist nicht zwingend so. Meter passt schon gut, weniger Kommafehler?!
3. Zitat kussji: "4x4 Punkte erscheint mir logischer. Von einem profesionellen Onlinerechendienst, welcher ebenfalls die biquadratische Interpolation nutzt, weiss ich, dass er mit den nächsten 9 Punkten (3x3) interpoliert!??? Gibt es hier eine Erklärung?"

...hier mal die aktuelle Version:
[Output] Resultat: 52.6917632

Code: Alles auswählen

import numpy

def biQuadInt(px, py, heights):
    pxi, pyi = int(numpy.floor(px)), int(numpy.floor(py))

    dx=px-pxi; dxx=dx*dx; dxxx=dxx*dx
    dy=py-pyi; dyy=dy*dy; dyyy=dyy*dy

    wx=(0.5 * (    - dx + 2.0 * dxx -       dxxx),
        0.5 * (2.0      - 5.0 * dxx + 3.0 * dxxx),
        0.5 * (      dx + 4.0 * dxx - 3.0 * dxxx),
        0.5 * (         -       dxx +       dxxx))
    
    wy=(0.5 * (    - dy + 2.0 * dyy -       dyyy),
        0.5 * (2.0      - 5.0 * dyy + 3.0 * dyyy),
        0.5 * (      dy + 4.0 * dyy - 3.0 * dyyy),
        0.5 * (         -       dyy +       dyyy))

    result = 0
    for a, wa in enumerate(wx, pxi-1):
        for b, wb in enumerate(wy, pyi-1):
            result += heights[b, a] * wa * wb
    return result

#------------------------            

heights = numpy.array([ [53.2, 53.6, 53.7, 53.5], [53.1, 53.3, 53.0, 52.9], [52.4, 52.7, 52.1, 52.0], [54.6, 54.2, 54.0, 54.5]])
#heights = numpy.array([  [0.0532, 0.0536, 0.0537, 0.0538], [0.0531, 0.0533, 0.0530, 0.0529], [0.0524, 0.0527, 0.0521, 0.0520], [0.0546, 0.0542, 0.0540, 0.0545]])
px, py = 1.300, 1.600 # für diesen Rasterpunkt die gesuchte Interpolierte Höhe
print("Resultat:",biQuadInt(px, py, heights))

[Sirius3]

Woher soll der Rechner wissen, was ein Kilometer oder was ein Meter ist. Bei den Werten ist die Einheit, die man hineingibt auch die Einheit, die herauskommt. Wenn Du sie mit Metern fütterst, ist das Ergebnis in Metern. Da float intern mit Exponenten arbeitet, ist die Genauigkeit die selbe, egal ob alles in Atomdurchmesser oder Lichtjahren angegeben ist.
Bei den Koordinaten ist das etwas anderes. Ein Array hat nunmal nur diskrete Indizes, [0,0], [0,1] ... [1231, 423]. Wenn Du die Höhen als (x,y,z) gegeben hast, dann mußtest Du x und y ja schon in diese Indizes umgerechnet haben, und das selbe mußt Du auch mit den gesuchten Koordinaten machen, nur dass es hier halt auch Kommazahlen erlaubt sind.

Ob andere Algorithmen mit 9 Stützstellen auskommen mag sein, dieser hier garantiert eine gewisse Symmetrie, weil immer zwei Stützstellen links und zwei rechts liegen, und nicht mal eine links und zwei rechts und mal umgekehrt.

[kussji]

Danke!

Bei den Koordinaten ist das etwas anderes. Ein Array hat nunmal nur diskrete Indizes, [0,0], [0,1] ... [1231, 423]. Wenn Du die Höhen als (x,y,z) gegeben hast, dann mußtest Du x und y ja schon in diese Indizes umgerechnet haben, und das selbe mußt Du auch mit den gesuchten Koordinaten machen, nur dass es hier halt auch Kommazahlen erlaubt sind.

Da stehe ich noch auf dem Schlauch .
Für die gesuchten Koordinaten (x,y) leuchtet mir das noch ein - die Umrechung. Aber für die x,y Koordinaten der Höhenwerte (z) sehe ich das noch nicht? Weil es ja ein reguläres x-y-Gitter ist. Dann müsste ich ja für x,y auch je ein Wertearray haben. Welches in der Funktion berücksichtig werden müsste?

Meine Überlegung ist, dass ich bei jeder Berechnung diese Gitter (4x4-Punkte) als relative Position betrachte und somit x,y immer von [[0,0], bis [3,3]] geht? Die gesuchte x,y-Koordinate müsste dann eine Kommazahl zwischen 1.0 und 2.0 sein, sowohl für x als auch für y

Wo ist mein Denkfehler bzw. denke ich zu weit?
Mit einem Zahlenbeispiel käme ich vielleicht weiter...

[Sirius3]

Genau, das hast Du richtig verstanden. Und wenn Du die selbe Rechenregel für die x- und y-Werte der Höhen anwendest bekommst Du 0 bis 3 als Ergebnis, weshalb man die ja auch nicht übergeben muß, weil es immer die selben Zahlen sind.

[kussji]

Dann habe ich ja Glück gehabt. Vielen Dank für die nette Hilfe.