Monotone Cubic Interpolation

[DragonT]

Hallo Community,

bräuchte mal die Hilfe von einem anderen Java Entwickler oder besser von einem Mathematiker.
Möchte für ein vielseitiges Projekt eine Monotone Kubische Interpolation vornehmen weil es zumindest laut Wiki Artikel, exakt das Verhalten ermöglicht das ich brauche (im Gegensatz zur standardmäßigen cubic Interpolation): Monotone cubic interpolation - Wikipedia

Hab dafür diesen Code hier gefunden: Spline interpolation in java. · GitHub
Angeblich soll es aus dem Sourcecode von Android stammen (oder wie ist die erste Komentar-Zeile zu interpretieren?).

Nun ergibt das für folgende Punkte:

X : Y
0.0 : 0.0
30.0 : 100
40.0 : 0.0

Folgenden Graphen ("manuell" gerendert, nicht mit Plotting Funktionen von JFX o.ä.) :

Die interpolierten Punkte schießen über das Maximum von 100 hinaus obwohl das Verfahren ja grade das verhindern soll.
Dies bestätigen auch die Werte. Der Scheitel ist bei 106.69:

Interpolated points: 0 : 0.0
Interpolated points: 1 : 3.5481482
Interpolated points: 2 : 7.4962964
Interpolated points: 3 : 11.799999
Interpolated points: 4 : 16.414816
Interpolated points: 5 : 21.296297
Interpolated points: 6 : 26.400002
Interpolated points: 7 : 31.681482
Interpolated points: 8 : 37.0963
Interpolated points: 9 : 42.600006
Interpolated points: 10 : 48.14815
Interpolated points: 11 : 53.6963
Interpolated points: 12 : 59.2
Interpolated points: 13 : 64.614815
Interpolated points: 14 : 69.89629
Interpolated points: 15 : 75.0
Interpolated points: 16 : 79.881485
Interpolated points: 17 : 84.49629
Interpolated points: 18 : 88.8
Interpolated points: 19 : 92.748146
Interpolated points: 20 : 96.296295
Interpolated points: 21 : 99.4
Interpolated points: 22 : 102.01482
Interpolated points: 23 : 104.09629
Interpolated points: 24 : 105.6
Interpolated points: 25 : 106.48148
Interpolated points: 26 : 106.6963
Interpolated points: 27 : 106.200005
Interpolated points: 28 : 104.94815
Interpolated points: 29 : 102.89629
Interpolated points: 30 : 100.0
Interpolated points: 31 : 95.4
Interpolated points: 32 : 88.53333
Interpolated points: 33 : 79.8
Interpolated points: 34 : 69.600006
Interpolated points: 35 : 58.333332
Interpolated points: 36 : 46.399994
Interpolated points: 37 : 34.200005
Interpolated points: 38 : 22.133331
Interpolated points: 39 : 10.600002
Interpolated points: 40 : 0.0

Wenn ich stattdessen diese Punkte nehme:

X : Y
0.0 : 0.0
20.0 : 100
40.0 : 0.0

Womit der höchste Punkt exakt in der Mitte ist, funktioniert es wie es soll. Der Scheitel ist dann bei 100.

Ist der Algorithmus also nicht für unterschiedliche Abstände auf der X-Achse geeignet?
Kann jemand den Code eventuell nachvollziehen?
Oder kennt ihr eine andere Open-Source Bibliothek mit dieser Funktionalität?
Hab bislang nur klassische Cubic-Interpolation gefunden.

Schon mal riesigen Dank im Vorraus!

 

[Thallius]

Also ich verstehe den Wiki Artikel anders. Das monoton ist zwar glatter als der Standard aber es ist immer noch ein Spline und dami5 ist es logisch das die Werte über 100 gehen denn sonstmwürde der Graph ja zwangsweise platt sein müssen an der Spitze und das wäre dann kein Spline mehr.

 

[DragonT]

Vielen dank für die schnelle Antwort!
Weshalb müsste der Spline denn platt sein? ich habe 3 Punkte durch die der Graph durch soll. Der Scheitel kann bzw. soll doch einfach entsprechend an dem höchsten Punkt sein.
Was ich hier habe scheint ein "Überschwingen" zu sein.
Auch ist die Regel der Monotonie gebrochen:
Zwischen 2 Punkten steigt der Wert aber sinkt auch wieder. Monoton würde heissen dass er nur steigt oder nur sinkt.

Oder meinst du dass es eine zwingende Eigenschaft eiens Splines wäre, dass der Punktabstand immer identisch ist?

 

[JCODA]

Hey,

habe mir das angeschaut und ist ziemlich witzig:
Der Wikipedia-Artikel wird 1 zu 1 in Code umgesetzt bis auf den Schritt 4. Dieser fehlt und ist genau das was du benötigst. Ich hab es mal ergänzt und kommentiert:

/*
* Copyright (C) 2012 The Android Open Source Project
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/

import java.util.Arrays;
import java.util.List;

/**
* Performs spline interpolation given a set of control points.
*
*/
public class SplineInterpolator {

    private final List<Float> mX;
    private final List<Float> mY;
    private final float[] mM;

    private SplineInterpolator(List<Float> x, List<Float> y, float[] m) {
        mX = x;
        mY = y;
        mM = m;
    }

    /**
     * Creates a monotone cubic spline from a given set of control points.
     *
     * The spline is guaranteed to pass through each control point exactly.
     * Moreover, assuming the control points are monotonic (Y is non-decreasing
     * or non-increasing) then the interpolated values will also be monotonic.
     *
     * This function uses the Fritsch-Carlson method for computing the spline
     * parameters. http://en.wikipedia.org/wiki/Monotone_cubic_interpolation
     *
     * @param x
     *            The X component of the control points, strictly increasing.
     * @param y
     *            The Y component of the control points
     * @return
     *
     * @throws IllegalArgumentException
     *             if the X or Y arrays are null, have different lengths or have
     *             fewer than 2 values.
     */
    public static SplineInterpolator createMonotoneCubicSpline(List<Float> x, List<Float> y) {
        if (x == null || y == null || x.size() != y.size() || x.size() < 2) {
            throw new IllegalArgumentException(
                    "There must be at least two control " + "points and the arrays must be of equal length.");
        }

        final int n = x.size();
        float[] d = new float[n - 1]; // could optimize this out
        float[] m = new float[n];

        // Compute slopes of secant lines between successive points.
        for (int i = 0; i < n - 1; i++) {
            float h = x.get(i + 1) - x.get(i);
            if (h <= 0f) {
                throw new IllegalArgumentException(
                        "The control points must all " + "have strictly increasing X values.");
            }
            d[i] = (y.get(i + 1) - y.get(i)) / h;
           
        }

        // Initialize the tangents as the average of the secants.
        m[0] = d[0];
        for (int i = 1; i < n - 1; i++) {
            m[i] = (d[i - 1] + d[i]) * 0.5f;
        }
        m[n - 1] = d[n - 2];

        // Update the tangents to preserve monotonicity.
        for (int i = 0; i < n - 1; i++) {
            if (d[i] == 0f) { // successive Y values are equal
                m[i] = 0f;
                m[i + 1] = 0f;
            } else {
                float a = m[i] / d[i];
                float b = m[i + 1] / d[i];               
                // Here was the missing step 4 of the Wikipedia article.
                // Possible mistake in wiki, because wiki wants to set m[i] = 0, but we need m[i+1] = 0.
                if (a < 0 || b < 0) {
                    m[i+1] = 0.0f;
                } else {
                    float h = (float) Math.hypot(a, b);
                    if (h > 3) {
                        float t = 3f / h;
                        m[i] = t * a * d[i];
                        m[i + 1] = t * b * d[i];
                    }
                }
            }
        }
        return new SplineInterpolator(x, y, m);
    }

    /**
     * Interpolates the value of Y = f(X) for given X. Clamps X to the domain of
     * the spline.
     *
     * @param x
     *            The X value.
     * @return The interpolated Y = f(X) value.
     */
    public float interpolate(float x) {
        // Handle the boundary cases.
        final int n = mX.size();
        if (Float.isNaN(x)) {
            return x;
        }
        if (x <= mX.get(0)) {
            return mY.get(0);
        }
        if (x >= mX.get(n - 1)) {
            return mY.get(n - 1);
        }

        // Find the index 'i' of the last point with smaller X.
        // We know this will be within the spline due to the boundary tests.
        int i = 0;
        while (x >= mX.get(i + 1)) {
            i += 1;
            if (x == mX.get(i)) {
                return mY.get(i);
            }
        }

        // Perform cubic Hermite spline interpolation.
        float h = mX.get(i + 1) - mX.get(i);
        float t = (x - mX.get(i)) / h;
        return (mY.get(i) * (1 + 2 * t) + h * mM[i] * t) * (1 - t) * (1 - t)
                + (mY.get(i + 1) * (3 - 2 * t) + h * mM[i + 1] * (t - 1)) * t * t;
    }

    // For debugging.
    @Override
    public String toString() {
        StringBuilder str = new StringBuilder();
        final int n = mX.size();
        str.append("[");
        for (int i = 0; i < n; i++) {
            if (i != 0) {
                str.append(", ");
            }
            str.append("(").append(mX.get(i));
            str.append(", ").append(mY.get(i));
            str.append(": ").append(mM[i]).append(")");
        }
        str.append("]");
        return str.toString();
    }

    public static void main(String[] args) {
        Float[] x = { 0F, 30F, 40F };
        Float[] y = { 0F, 100F, 0F };

        SplineInterpolator si = SplineInterpolator.createMonotoneCubicSpline(Arrays.asList(x), Arrays.asList(y));

        for (int i = 0; i < 42; i++) {
            System.out.println(i + " " + si.interpolate(i));
        }
    }

}

 

[DragonT]

Hallo, tut mir leid dass ich erst jetzt antworte. Die Hochschule und anderes Zeug haben mich eingeholt und ich hatte mich zeitweise mit der imperfekten Lösung abgefunden.

@JCODA
Deine Lösung wie sie ist erzeugt leider Überschwinger unter anderen Umständen. Sie hat mich aber auf auf die richtige Spur gebracht!

Den Else-Zweig in der zwiten For-Schleife ersetzen durch folgenden Code:

float a = m[i] / d[i];
float b = m[i + 1] / d[i];
           
if (a <= 0)
     m[i] = 0.0f;
           
if (b <= 0)
     m[i+1] = 0.0f;
           
float h = (float) Math.hypot(a, b);
if (h > 3)
{
     float t = 3f / h;
     m[i] = t * a * d[i];
     m[i + 1] = t * b * d[i];
}

Dann funktioniert es perfekt!

Daher Danke!