pondelok, 18 jún 2012 09:53 Written by 6322 times
Rate this item
(6 votes)

C# - Spline interpolácia

Spline interpolácia je veľmi presná technika interpolovania známych (zadaných) bodov bez toho, aby bolo nutné poznať aktuálnu charakteristiku (krivku).

 

Interpolácia pozostáva z odhadu hodnot (bodov) funkcie medzi zadanými hodnotami (x0, y0), (x1, y1), ... (xi-1, yi-1), (xi, yi), (xi+1, yi+1), ... (xn, yn). Hodnoty sa interpolujú vždy medzi dvojicou zadaných hôdnot, čize pre nejaký interval. Pričom do výpočtu hodnôt z daného intervalu vtupujú hodboty z predchádzajúceho aj nasledujúcého intervalu, tým sa dosiahne hladši a presnejši priebeh interpolovanej funkcie. Podrobné vysvetlenie a odvodzovanie výrazov pre spline interpoláciu je voľne dostupné na webe. Takže len pre informáciu sú nižsie uvedené už len odvodené vzorce, ktoré sa využivajú v programe. Cieľom je vypočitať koeficienty "a, b, c, d" pre jednotlivé intervaly a pomocou polynómu tretieho rádu získať interpolované body funkcie a teda priebeh interpolovanej funkcie.

 

Polynóm tretieho rádu pre výpočet bodov funkcie:

spline1

Na výpočet jednotlivých koeficientov je potrebné poznat prvú a druhú deriváciu:

spline2

Po úprave (úprava je nutná aby v zdrojovom kóde nenastalo delenie nulov, v prípade ak je y súradnica funkcie rovnaká pre rôzne x):

spline8

Pričom pre krajné (hraničné) body platí:

spline3

Druhé derivácie:

spline4

Koeficienty pre jednotlive intevaly zadaných hodnôt sa vypočítajú podľa vzorcov:

spline5

spline6

Ako príklad je uvedená interpolácia pre 6 intervalov, tzn. pre 7 zadaných súradníc (x,y):

spline7

Graf interpolácie - zadané body (červená farba), vypočítané body (modrá farba):

splineGraf

Zdrojový kód

public bool splineInterpolation(out string errorMessage, double[] xInArray, double[] yInArray, out double[] yOutArray, out double[,] koeffArray)
        {
            errorMessage = "";
            if (xInArray.Length < 3)
            {
                errorMessage = "Array Lenght Error: Vstupne pole (polia) ma nespravnu dlzku.";
                yOutArray = new double[0];
                koeffArray = new double[0, 0];
                return false;
            }
            int pointCount = Convert.ToInt32(xInArray[xInArray.Length - 1]) - Convert.ToInt32(xInArray[0]) + 1;
            yOutArray = new double[pointCount];
            double f1x0_deriv = 0;
            double f1x1_deriv = 0;
            double f2x0_deriv = 0;
            double f2x1_deriv = 0;
            int intervalLenght = xInArray.Length - 1;
            koeffArray = new double[intervalLenght, 4];
            try
            {
                if (xInArray.Length != yInArray.Length)
                {
                    errorMessage = "Array Lenght Error: Vstupne pole (polia) ma nespravnu dlzku.";
                    return false;
                }
                //Vypocet koeficientov pre jednotlive intervaly
                for (int iInterval = 0; iInterval < intervalLenght; iInterval++)
                {
                    if (iInterval == 0)
                    {
                        f1x1_deriv = (2 * (yInArray[iInterval + 2] - yInArray[iInterval + 1])
                                    * (yInArray[iInterval + 1] - yInArray[iInterval]))
                                    / ((xInArray[iInterval + 2] - xInArray[iInterval + 1])
                                    * (yInArray[iInterval + 1] - yInArray[iInterval])
                                    + (xInArray[iInterval + 1] - xInArray[iInterval])
                                    * (yInArray[iInterval + 2] - yInArray[iInterval + 1]));
                        f1x0_deriv = ((3 * (yInArray[iInterval + 1] - yInArray[iInterval]))
                                    / (2 * (xInArray[iInterval + 1] - xInArray[iInterval])))
                                    - (f1x1_deriv / 2);
                        f2x0_deriv = ((-2 * (f1x1_deriv + 2 * f1x0_deriv))
                                    / (xInArray[iInterval + 1] - xInArray[iInterval]))
                                    + ((6 * (yInArray[iInterval + 1] - yInArray[iInterval]))
                                    / Math.Pow(xInArray[iInterval + 1] - xInArray[iInterval], 2));
                        f2x1_deriv = ((2 * (2 * f1x1_deriv + f1x0_deriv))
                                    / (xInArray[iInterval + 1] - xInArray[iInterval]))
                                    - ((6 * (yInArray[iInterval + 1] - yInArray[iInterval]))
                                    / Math.Pow(xInArray[iInterval + 1] - xInArray[iInterval], 2));
                    }
                    else
                    {
                        if (iInterval == intervalLenght - 1)
                        {
                            f1x0_deriv = (2 * (yInArray[iInterval + 1] - yInArray[iInterval])
                                        * (yInArray[iInterval] - yInArray[iInterval - 1]))
                                        / ((xInArray[iInterval + 1] - xInArray[iInterval])
                                        * (yInArray[iInterval] - yInArray[iInterval - 1])
                                        + (xInArray[iInterval] - xInArray[iInterval - 1])
                                        * (yInArray[iInterval + 1] - yInArray[iInterval]));
                            f1x1_deriv = ((3 * (yInArray[iInterval + 1] - yInArray[iInterval]))
                                        / (2 * (xInArray[iInterval + 1] - xInArray[iInterval])))
                                        - (f1x0_deriv / 2);
                            f2x0_deriv = ((-2 * (f1x1_deriv + 2 * f1x0_deriv))
                                        / (xInArray[iInterval + 1] - xInArray[iInterval]))
                                        + ((6 * (yInArray[iInterval + 1] - yInArray[iInterval]))
                                        / Math.Pow(xInArray[iInterval + 1] - xInArray[iInterval], 2));
                            f2x1_deriv = ((2 * (2 * f1x1_deriv + f1x0_deriv))
                                        / (xInArray[iInterval + 1] - xInArray[iInterval]))
                                        - ((6 * (yInArray[iInterval + 1] - yInArray[iInterval]))
                                        / Math.Pow(xInArray[iInterval + 1] - xInArray[iInterval], 2));
                        }
                        else
                        {
                            f1x1_deriv = (2 * (yInArray[iInterval + 2] - yInArray[iInterval + 1])
                                        * (yInArray[iInterval + 1] - yInArray[iInterval]))
                                        / ((xInArray[iInterval + 2] - xInArray[iInterval + 1])
                                        * (yInArray[iInterval + 1] - yInArray[iInterval])
                                        + (xInArray[iInterval + 1] - xInArray[iInterval])
                                        * (yInArray[iInterval + 2] - yInArray[iInterval + 1]));
                            f1x0_deriv = (2 * (yInArray[iInterval + 1] - yInArray[iInterval])
                                        * (yInArray[iInterval] - yInArray[iInterval - 1]))
                                        / ((xInArray[iInterval + 1] - xInArray[iInterval])
                                        * (yInArray[iInterval] - yInArray[iInterval - 1])
                                        + (xInArray[iInterval] - xInArray[iInterval - 1])
                                        * (yInArray[iInterval + 1] - yInArray[iInterval]));
                            f2x0_deriv = ((-2 * (f1x1_deriv + 2 * f1x0_deriv))
                                        / (xInArray[iInterval + 1] - xInArray[iInterval]))
                                        + ((6 * (yInArray[iInterval + 1] - yInArray[iInterval]))
                                        / Math.Pow(xInArray[iInterval + 1] - xInArray[iInterval], 2));
                            f2x1_deriv = ((2 * (2 * f1x1_deriv + f1x0_deriv))
                                        / (xInArray[iInterval + 1] - xInArray[iInterval]))
                                        - ((6 * (yInArray[iInterval + 1] - yInArray[iInterval]))
                                        / Math.Pow(xInArray[iInterval + 1] - xInArray[iInterval], 2));
                        }
                    }
                    //koeficienty a, b, c, d - stlpce, intervaly - riadky
                    koeffArray[iInterval, 0] = (f2x1_deriv - f2x0_deriv)
                                            / (6 * (xInArray[iInterval + 1] - xInArray[iInterval]));
                    koeffArray[iInterval, 1] = (xInArray[iInterval + 1] * f2x0_deriv - xInArray[iInterval]
                                            * f2x1_deriv) / (2 * (xInArray[iInterval + 1] - xInArray[iInterval]));
                    koeffArray[iInterval, 2] = ((yInArray[iInterval + 1] - yInArray[iInterval])
                                            - koeffArray[iInterval, 1] * (Math.Pow(xInArray[iInterval + 1], 2)
                                            - Math.Pow(xInArray[iInterval], 2))
                                            - koeffArray[iInterval, 0] * (Math.Pow(xInArray[iInterval + 1], 3)
                                            - Math.Pow(xInArray[iInterval], 3)))
                                            / (xInArray[iInterval + 1] - xInArray[iInterval]);
                    koeffArray[iInterval, 3] = yInArray[iInterval] - koeffArray[iInterval, 2] * xInArray[iInterval]
                                            - koeffArray[iInterval, 1] * Math.Pow(xInArray[iInterval], 2)
                                            - koeffArray[iInterval, 0] * Math.Pow(xInArray[iInterval], 3);
                }
                int start = Convert.ToInt32(xInArray[0]);
                int intervalNum = 0;
                for (int i = start; i < pointCount + start; i++)
                {
                    while ((Convert.ToDouble(i) >= xInArray[intervalNum] && Convert.ToDouble(i) <= xInArray[intervalNum + 1]) != true)
                    {
                        intervalNum++;
                    }
                    yOutArray[i - start] = koeffArray[intervalNum, 3] + koeffArray[intervalNum, 2] * i
                                         + koeffArray[intervalNum, 1] * Math.Pow(i, 2)
                                         + koeffArray[intervalNum, 0] * Math.Pow(i, 3);
                }
            }
            catch (Exception ex)
            {
                errorMessage = ex.Message;
                return false;
            }
            return true;
        }
Last modified on pondelok, 19 december 2016 21:24
Ing.Peter Šuba

Zakladateľ www.projectik.eu

It's time for another revolution

(Why so serious?)

Website: www.projectik.eu