Сообщество Engee
Engee Community
/
Примеры и проекты
/
Применение базовых алгоритмов интерполяции

Применение базовых алгоритмов интерполяции

Автор
avatar-bogoslmbogoslm
Notebook

Подготовка и объявление базовых функций

В этой части скрипта объявим все необходимые функции для того чтобы успешно и удобно показать работу алгоритмов интерполяции, а также дадим некоторые пояснения

Библиотеки

Нам понадобится библиотека Plots для того чтобы визуализировать наши данные, и сравнить работу алгоритмов

In [ ]: 
using Plots

Функция генерации сетки координат X, Y

Эта функция генерирует две матрицы, которые представляют собой сетку координат X и Y, исходя из векторов x и y

In [ ]: 
function meshgrid(x, y)
    X = transpose([i for i in x, j in 1:length(y)])
    Y = transpose([j for i in 1:length(x), j in y])
    return Matrix(X), Matrix(Y)
end;

Функции оборачивания данных сетки

Для того чтобы успешно решить задачу интерполяции нужно посчитать значения в граничных точках рассматриваемого диапазона. Один из используемых методов для решения этой части задачи - дублирование значений граничных точек (Z) в соседних точках диапазона (функция frame_points). После дублирования решают задачу интерполяции, и в конце отбрасывают лишние точки (функция unframe_points). Здесь для упрощения мы взяли точки, отстоящие от границ исходного интервала на шаг интерполирования.

In [ ]: 
function frame_points(X0, Y0, Z0, step)     
    Y0 = hcat(Y0[:, 1], Y0)
    Y0 = hcat(Y0, Y0[:, end])
    Y0 = Matrix(transpose(Y0))
    Y0 = hcat(Y0[:, 1] .- step, Y0)
    Y0 = hcat(Y0, Y0[:, end] .+ step)
    Y0 = Matrix(transpose(Y0))

    X0 = Matrix(transpose(X0))
    X0 = hcat(X0[:, 1], X0)
    X0 = hcat(X0, X0[:, end])
    X0 = Matrix(transpose(X0))
    X0 = hcat((X0[:, 1] .- step), X0)
    X0 = hcat(X0, (X0[:, end] .+ step))

    Z0 = hcat(Z0[:, 1], Z0)
    Z0 = hcat(Z0, Z0[:, end])
    Z0 = Matrix(transpose(Z0))
    Z0 = hcat(Z0[:, 1], Z0)
    Z0 = hcat(Z0, Z0[:, end])
    Z0 = Matrix(transpose(Z0))

    return X0, Y0, Z0
end


function unframe_points(X0, Y0, Z0)
    return X0[2:end-1, 2:end-1], Y0[2:end-1, 2:end-1], Z0[2:end-1, 2:end-1];
end;

Функции генерации исходных данных F(x, y)

Это функции, при помощи которых мы будем генерировать сетку данных. Для примера сделали несколько вариантов

In [ ]: 
function f1(x, y)
    return x^2 + y^2
end

function f2(x,  y)
    absz = (2*x^2 + y^2)^0.5
    return [absz, -1*absz]
end 

function f3(x, y)
    return abs(x) + abs(y)
end 

function f4(x, y)
    return abs(100 - x^2 + y^2)^0.5
end 

function f5(x, y)
    return y * sin(x) - x * cos(y);
end

function f6(x, y)
    return f3(x, y) * sin(x)
end;

Это функция генерации сетки данных на основании заданных интервалов x и y, а также функции рассчета кординаты z

In [ ]: 
function get_surfe(base_function, x_vector, y_vetor)
    X, Y = meshgrid(x_vector, y_vetor)
    Z = zeros(size(X))
    for i in 1:length(X)
        x = X[i]
        y = Y[i]
        Z[i] = base_function(x, y)
    end
    return X, Y, Z
end;

Решение задачи интерполяции

Функции интерполяции

Рассмотрим два простейших метода интерполяции - линейный и кубический

Билинейный и бикубический методы от линейного и кубического отличаются лишь тем, что приминяются дважды - первый вдоль одной оси, а второй вдоль другой для уже полученных данных.
Функция билинейной интерполяции

Линейный метод интерполяции — это метод, который используется для нахождения промежуточных значений функции между двумя известными точками. Он основан на предположении, что функция изменяется линейно между этими точками.

Формула


Если у нас есть две точки и , которые лежат на графике функции , то мы можем использовать следующую формулу для нахождения значения функции в любой точке между и :

где — искомое значение функции, — точка, в которой мы хотим найти значение функции.

Этот метод широко используется в математике, физике, инженерии и других областях науки и техники для аппроксимации функций и нахождения промежуточных значений.

Билинейная интерполяция

Для случая функции от двух переменных используют билинейную интерполяцию. Ее алгоритм работает на основании принципа линейной интерполяции, но применяется дважды (вдоль осей X и Y)

In [ ]: 
function bilinear_interp(X0, Y0, Z0, step)

    X0, Y0, Z0 = frame_points(X0, Y0, Z0, step)
    
    s = size(X0)
    
    sy = s[1]
    sx = s[2]

    x_line = [X0[1]:step:X0[end];]
    y_line = [Y0[1]:step:Y0[end];]

    nx = length(x_line)
    ny = length(y_line)


    X, Y = meshgrid(x_line, y_line)
    Z = zeros(nx, ny)
    
    ix0 = 2
    for i in 1:nx
        if x_line[i] > X0[1, ix0] || abs(x_line[i] - X0[1, ix0]) < 1e-14
            ix0 = ix0 + 1;
            if ix0 > sx
                break;
            end
        end
        iy0 = 2
        for j in 1:ny
            if (y_line[j] > Y0[iy0, 1] || abs(y_line[j] - Y0[iy0, 1]) < 1e-14)
                iy0 = iy0 + 1;
                if iy0 > sy
                    break;
                end
            end
            
            curr_x = x_line[i]
            curr_y = y_line[j]
            k1 = (Z0[iy0, ix0 - 1] - Z0[iy0 - 1, ix0 - 1]) / (Y0[iy0, ix0 - 1] - Y0[iy0 - 1, ix0 - 1])
            b1 = Z0[iy0 - 1, ix0 - 1] - k1 * Y0[iy0 - 1, ix0 - 1]

            k2 = (Z0[iy0, ix0] - Z0[iy0 - 1, ix0]) / (Y0[iy0, ix0] - Y0[iy0 - 1, ix0])
            b2 = Z0[iy0 - 1, ix0] - k2 * Y0[iy0 - 1, ix0]

            z1 = k1 * curr_y + b1
            z2 = k2 * curr_y + b2

            k3 = (z2 - z1) / (X0[iy0, ix0] - X0[iy0, ix0 - 1])
            b3 = z1 - k3 * X0[iy0, ix0 - 1]

            res = k3 * curr_x + b3

            Z[j, i] = res
        end
    end 

    X, Y, Z = unframe_points(X, Y, Z)
    return X, Y, Z
end;
Функция бикубической интерполяции

Кубический метод интерполяции — это метод, который используется для нахождения промежуточных значений функции между четырьмя известными точками. Он основан на предположении, что функция может быть аппроксимирована кубическим полиномом.

Формула
Если у нас есть четыре точки , , и , которые лежат на графике функции , то мы можем использовать следующую формулу для нахождения значения функции в любой точке между и :

где , , и — коэффициенты кубического полинома, которые можно найти из системы уравнений, составленных на основе известных точек.

Этот метод позволяет получить более гладкую аппроксимацию функции по сравнению с линейной интерполяцией, особенно если функция имеет сложную форму. Однако он также требует больше вычислений и может быть менее точным вблизи границ интервала интерполяции. Заметим что формула представлена полиномом ньютона, однако высчитывать коэфициенты - не такая простая задача. Воспользуемся фактом, что через четыре заданные точки проходит только одна кубическая функция, и апроксимируем эту функцию при помощи палинома лагранжа. Это и будет искомая кубическая функция В нашем примере, для упрощения функция `lagrange` сразу возвращает значение палинома в точке, исходя из соседних точек и заданного x, вместо того чтобы возвращать коэфициенты.

Бикубическая интерполяция

Для случая функции от двух переменных используют бикубическую интерполяцию. Ее алгоритм работает на основании принципа кубической интерполяции, но применяется дважды (вдоль осей X и Y)

In [ ]: 
function lagrange(x, x1, y1, x2, y2, x3, y3, x4, y4)
    l1 = y1 * (x - x2) * (x - x3) * (x - x4) / ((x1 - x2) * (x1 - x3) * (x1 - x4))
    l2 = y2 * (x - x1) * (x - x3) * (x - x4) / ((x2 - x1) * (x2 - x3) * (x2 - x4))
    l3 = y3 * (x - x1) * (x - x2) * (x - x4) / ((x3 - x1) * (x3 - x2) * (x3 - x4))
    l4 = y4 * (x - x1) * (x - x2) * (x - x3) / ((x4 - x1) * (x4 - x2) * (x4 - x3))

    val = l1 + l2 + l3 + l4
    return val
end


function bicubic_interp(X0, Y0, Z0, step)
    for _ in 1:2
        X0, Y0, Z0 = frame_points(X0, Y0, Z0, step)
    end

    s = size(X0)
    sy = s[1]
    sx = s[2]


    x_line = X0[1]:step:X0[end]
    y_line = Y0[1]:step:Y0[end]

    nx = length(x_line)
    ny = length(y_line)

    X, Y = meshgrid(x_line, y_line)
    Z = zeros(nx, ny)
    
    ix0 = 3
    for i in 1:nx
        if x_line[i] > X0[1, ix0] || isequal(x_line[i], X0[1, ix0])
            ix0 = ix0 + 1
            if ix0 > sx - 1
                break
            end
        end
        iy0 = 3
        for j = 1:ny
            if y_line[j] >= Y0[iy0, 1] || isequal(y_line[j], Y0[iy0, 1])
                iy0 = iy0 + 1
                if iy0 > sy - 1
                    break;
                end
            end
            
            curr_x = x_line[i]
            curr_y = y_line[j]

            zz = zeros(1, 4);
            for ind = 1:4
                k = ind - 3
                z = lagrange(
                    curr_x,
                    X0[iy0 + k, ix0 - 2], Z0[iy0 + k, ix0 - 2],
                    X0[iy0 + k, ix0 - 1], Z0[iy0 + k, ix0 - 1],
                    X0[iy0 + k, ix0], Z0[iy0 + k, ix0],
                    X0[iy0 + k, ix0 + 1], Z0[iy0 + k, ix0 + 1])
                zz[ind] = z
            end
            
            res = lagrange(
                curr_y,
                Y0[iy0 - 2, ix0], zz[1],
                Y0[iy0 - 1, ix0], zz[2],
                Y0[iy0, ix0], zz[3],
                Y0[iy0 + 1, ix0], zz[4])
            Z[j, i] = res
        end
    end

    for _ in 1:2
        X, Y, Z = unframe_points(X, Y, Z)
    end

    return X, Y, Z
end;

Запуск генерации данных

Зададим интервал интерполирования и шаг интерполяции. Получим сетку со значениями. С этого момента предположим, что мы не знаем истинных значений в точках (не используем исходную функцию)

In [ ]: 
x_vector = [-5:5;];
y_vector = [-5:5;];

interpolated_data_step = 0.05;
In [ ]: 
X, Y, Z = get_surfe(f5, x_vector, y_vector);

Запуск интерполяции данных

Запустим интерполяцию и получим предположительные данные для промежуточных точек...

In [ ]: 
Xl, Yl, Zl = bilinear_interp(X, Y, Z, interpolated_data_step);
In [ ]: 
Xb, Yb, Zb = bicubic_interp(X, Y, Z, interpolated_data_step);

Построим графики

3D поверхности

так выглядят наши поверхности. Заметим, исходя из формы графика между точками, что для построения поверхности, встроенный интерфейс библиотеки Plots использует линейный метод интерполяции (что в целом логично...)

In [ ]: 
s1 = surface(X, Y, Z, title="input data", c=:grayC);
s2 = surface(Xl, Yl, Zl, title="bilinear interpolated", c=:grayC);
s3 = surface(Xb, Yb, Zb, title="bicubic interpolated", c=:grayC);
In [ ]: 
s = plot(s1, s2, s3, layout = (1, 3))
Out[0]:
2d диаграммы (тепловые карты)

Для большей наглядности воспользуемся интерфейсом построения тепловой карты.

Можно заметить, что бикубическая интерполяция данных выглядит более гладко, чем билинейная. Именно по этому данный метод интерполяции находит широкое применение в области графики и обработки изображений

In [ ]: 
h1 = heatmap(Z, title="input data");
h2 = heatmap(Zl, title="bilinear interpolation");
h3 = heatmap(Zb, title="bicubic interpolation");
In [ ]: 
plot(h1, h2, h3, layout = (1, 3))
Out[0]:
Попробуем интерполировать рисунок сердечка...
In [ ]: 
z = [
    255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255;
    255 255 255 255 255 240 230 240 255 255 240 230 240 255 255 255 255 255 255 255;
    255 255 255 255 230 0   0   0   230 240 0   0   0   230 255 255 255 255 255 255;
    255 255 255 240 0   0   0   0   0   0   0   0   0   0   240 255 255 255 255 255;
    255 255 240 0   0   0   0   0   0   0   0   0   125 0   0   230 255 255 255 255;
    255 255 125 0   0   0   0   0   0   0   0   125 255 0   0   240 255 255 255 255;
    255 255 125 0   0   0   0   0   0   0   0   0   0   0   0   220 255 255 255 255;
    255 255 125 0   0   0   0   0   0   0   0   0   0   0   0   200 255 255 255 255;
    255 255 125 0   0   0   0   0   0   0   0   0   0   0   0   240 255 255 255 255;
    255 255 125 0   0   0   0   0   0   0   0   0   0   0   0   230 255 255 255 255;
    255 255 255 125 0   0   0   0   0   0   0   7   0   0   240 255 255 255 255 255;
    255 255 255 125 0   0   0   0   0   0   0   0   0   0   230 255 255 255 255 255;
    255 230 255 255 125 0   0   0   0   0   0   0   0   240 255 255 255 255 255 255;
    255 255 255 255 255 125 0   0   0   0   0   0   240 255 255 255 255 200 255 255;
    255 255 255 255 255 255 125 0   0   0   0   218 255 255 255 255 255 255 255 255;
    255 240 255 255 255 255 255 125 0   0   250 255 255 255 255 255 255 255 255 255;
    255 255 255 255 220 255 255 255 125 125 255 255 255 255 255 255 255 255 255 255;
    255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255;
    255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255;
    255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
]

zdog = reverse(z')';

x_line = [1:size(zdog)[1];];
y_line = [1:size(zdog)[2];];

xdog, ydog = meshgrid(x_line, y_line);
In [ ]: 
xldog, yldog, zldog = bilinear_interp(xdog, ydog, zdog, interpolated_data_step);
xbdog, ybdog, zbdog = bicubic_interp(xdog, ydog, zdog, interpolated_data_step);
In [ ]: 
h1 = heatmap(zdog, title="input data", c=:Paired_10);
h2 = heatmap(zldog, title="bilinear interpolation", c=:Paired_10);
h3 = heatmap(zbdog, title="bicubic interpolation", c=:Paired_10);
In [ ]: 
plot(h1, h2, h3, layout=(1, 3))
Out[0]:
Группа видимости
ПубличноПросматривать контент могут все пользователи сообщества

Тип

Интерактивные скрипты

Краткое описание

Описали, как работает базовые алгоритмы интерполяции данных для функции двух переменных

Связанные материалы

Пусто

Теги

Категории

МероприятияКонкурс работ Engee 2024

Языки

    Julia

Форматы

    ngscript

Уровень

    Базовый

Статус

Опубликовано

Обновлено

Источник

Сообщество