VPF::[Алгоритмы] Двумерная интерполяция

Burka

Дата 25.6.2007, 12:43 (ссылка)

(нет голосов)

Загрузка ...

Опытный

Профиль
Группа: Участник
Сообщений: 429
Регистрация: 11.5.2007

Репутация: 1
Всего: 2

У меня сложность в реализации двумерной сплайн интерполяции.
Для одномерной интерполяции все просто(есть исходники smile

).
Например: дан массив известных точек A[N] и для каждого необходимого X мы вызываем функцию допустим y=spline(X1, A, N)
Таким образом мы получаем необходимые значения функции. Все просто.

Но я не могу разобраться с двумерной интерполяцией.
Нашел, на мой взгляд, хорошее описание того как это делается:
http://community.livejournal.com/ru_algorithms/24344.html
Цитата из ЖЖ:

Цитата

Так там и формул почти нет, все совсем просто. Стартуем с одномерного линейного или кубического сплайна. (Одномерный линейный сплайн — просто ломаная, проведенная через узлы мнтерполяции. Кубический сплайн описан почти в любой книжке.) При переходе от одномерного линейного сплайна к билинейному и от одномерного кубического сплайна к бикубическому используется идея т.н. «переразложения коэффициентов». (Кстати, она же превращает одномерное преобразование Фурье в двумерное.) Ищем функцию вида f(x,y). Для всех заданных значений x проводим одномерную интерполяцию по переменной y, используя линейный/кубический сплайн. Для каждого x получается набор найденных коэффициентов сплайна. Каждый коэффициент оказывается зависящим от x т.е. для каждого x получаются свое значение для каждого коэффициента. Теперь остается для каждого коэффициента построить интерполяцию по x. После этого можно подставить формулы для параметров (как функций от x) в формулу для зависимости от y. Совершенно так же идея переразложения работает и при произвольном числе переменных.

С математикой туго, поэтому разобратся как это сделать я не смог.
Прошу помощи в реализации двумерной сплайн интерполяции на основе одномерной. smile

--------------------

Великие умы обсуждают идеи; средние умы обсуждают события; мелкие умы обсуждают людей.

HoTMetaL

Дата 29.6.2007, 16:55 (ссылка)

(нет голосов)

Загрузка ...

Шустрый

Профиль
Группа: Участник
Сообщений: 107
Регистрация: 15.4.2007

Репутация: 1
Всего: 4

Предупреждаю - по вычмату у меня за экзамен "удовл." smile

Тем не менее, если у тебя есть функция одномерной интерполяции, то делаешь следующее:
Значение, полученное на каждом шаге интерполяции, запоминаешь и полученный массив загоняешь снова в функцию одномерной итнерполяции.
y
^
| * * * * *
| * * * * *
| * * * * *
| * * * * *
-------------->x
Например, сначала интерполируешь 5 раз по Х, а потом 1 раз по Y. или наоборот - 4 раза по Y и 1 - по X.

PS Тебя "удовл." мой ответ? smile

Это сообщение отредактировал(а) HoTMetaL - 29.6.2007, 17:00

Burka

Дата 29.6.2007, 19:49 (ссылка)

(нет голосов)

Загрузка ...

Опытный

Профиль
Группа: Участник
Сообщений: 429
Регистрация: 11.5.2007

Репутация: 1
Всего: 2

Цитата

PS Тебя "удовл." мой ответ? smile

Не очень))

Вот смотри. Есть такая функция:
y=spline(X, A, N),
где A - это массив известных точек по Х, и где N это количество этих известных точек.
X - Это X, Y которого мы хотим получить(X может входить в A)

то есть прокрутив эту функцию в цикле, меняя X (допустим от 1 до 10) мы получим массив Y[M], в котором точки более или менее похожи на нужные.
Тут ясно.

А теперь мне так же объясни что делать этой функцией для двух переменных)))
Допустим у нас есть двумерный массив известных точек B[K, F]
Вот дальше что-то никак вникнуть не могу.. что с ним делать.

Это сообщение отредактировал(а) Burka - 29.6.2007, 19:56

--------------------

Великие умы обсуждают идеи; средние умы обсуждают события; мелкие умы обсуждают людей.

HoTMetaL

Дата 30.6.2007, 09:00 (ссылка)

(нет голосов)

Загрузка ...

Шустрый

Профиль
Группа: Участник
Сообщений: 107
Регистрация: 15.4.2007

Репутация: 1
Всего: 4

Двумерная интерполяция означает, что нам дана функция от двух переменных z=f(x,y). Для удобства будем считать, что значения известны на прямоугольной сетке интерполяции - массив B[K][F]. Допустим, что нам нужно найти одно значение z0 в точке (x0,y0), т.е. z0=f(x0,y0). Поступим следующим образом: K раз интерполируем функцию по X, т.е.

1)y=y1;
2)y=y2;
3)y=y3;
...........
K)y=yK;

Затем один раз интерполируем по Y.

Таким образом, в рамках нашей задачи делаем следующее:
1) Создаём одномерный массив Temp[K]
2) В цикле заполняем этот массив значениями-результатами работы функции spline();.
3) На выход выдаём значение функции spline(); при подстановке в неё Temp[K].
Т.е. на шаге 2 в функцию spline(); нужно подставлять строки длиной F (K раз) из массива B[K][F].

Это может выглядеть примерно так:

Код


//
#include "spline.h" //модуль с твоей функцией spline();
#include <stdio.h>
#include <conio.h>

#define K 4
#define F 5
#define x0 4.5 // координата x0
#define y0 3.6 // координата y0

float B[K][F]={/*Инициализация массива*/};

void main()
{
   float Stroka[F];// В этот массив будем загонять строки из B[K][F]
   float Temp[K];// Массив для интерполяции по Y.
   float z0;//результат работы программы

   clrscr();// очистка экрана )))

  /*Интерполяция по X*/
   for(int i=0;i<K;i++){
      for(int j=0;j<F;j++)
         Stroka[j]=B[i][j];
      Temp[i]=spline(x0,Stroka,F);
   }
  /*Интерполяция по Y*/
   z0=spline(y0,Temp,K);
   printf("f(%f,%f)=%f",x0,y0,z0);
   getch();
}

Порядок интерполяции можно менять, т.е. сначала пробегаем по столбцам, а потом интерполируем по X.

А теперь "удовл."? smile

Это сообщение отредактировал(а) HoTMetaL - 30.6.2007, 09:29

Burka

Дата 30.6.2007, 15:56 (ссылка)

(нет голосов)

Загрузка ...

Опытный

Профиль
Группа: Участник
Сообщений: 429
Регистрация: 11.5.2007

Репутация: 1
Всего: 2

Так, вопросик есть.

Код



 var
  StrokaX, //
  StrokaY, // 
  Temp    :TReal1DArray; //Одномерные массивы из Double
  MyCtbl : TReal2DArray; //Двумерный массив из double
  i ,j ,Numbers :integer; //Целые

begin

 for i:=0 to K do
   begin
    Numbers :=0; //Сбрасываем счетчик
      for j:=0 to F do //Перебор строк
       if B[i,j]<>0 then //Формируем массивы известных точек
        begin
         StrokaX[Numbers]:=j; //здесь запоминаем значения X
         StrokaY[Numbers]:=B[i,j]; // А здесь значения Y соответствующего X
       end;

   Spline3BuildTable(Numbers, 1, StrokaX, StrokaY, 1, 1, MyCtbl); //Формируем таблицу коэфицентов
   Temp[i]:=  Spline3Interpolate(Numbers, MyCtbl, i); //И заполняем массив для интерполяции по Y
  end;

Извиняюсь что первый пример немного отличается от реальности, хотел упростить.
Вот сложность у меня в создании цикла для заполнения всей матрицы
То есть, мне нужно заполнить матрицу Z[S][D], в которой S>K, D>F так как матрица B входит в Z.

Код


z0=spline(y0,Temp,K);

По-моему это, как минимум, должно быть в одиночном цикле.
-------------------------------------------------
Вроде в таком роде что-то должно быть. smile

Код


  /*Интерполяция по X*/
   for(int i=0;i<K;i++)
    {
      for(int j=0;j<F;j++) Stroka[j]=B[i][j];

      Temp[i]=spline(x0,Stroka,F);

      for(int t=0; t<k; t++) z[i,t]=spline(t,Temp,k); //Т.Е. тут мы как раз и получаем значения Z, как мне кажется -\
    }

UPD че, то стормозил, последний вариант не подходит, так как массив Temp не сформирован полностью.

Это сообщение отредактировал(а) Burka - 30.6.2007, 16:43

--------------------

Великие умы обсуждают идеи; средние умы обсуждают события; мелкие умы обсуждают людей.

HoTMetaL

Дата 30.6.2007, 20:09 (ссылка)

(нет голосов)

Загрузка ...

Шустрый

Профиль
Группа: Участник
Сообщений: 107
Регистрация: 15.4.2007

Репутация: 1
Всего: 4

К первому листингу

Один сплошной баг. Ты его вообще запускал? При заполнении StrokaX[] и StrokaY[] ошибок разве не возникает?
Spline3BuildTable() - что это за функция?
Number где изменяется?

Цитата

Извиняюсь что первый пример немного отличается от реальности, хотел упростить.

А вот это зря.

Код


if B[i,j]<>0 then //Формируем массивы известных точек

Это ещё чё такое?

Ко второму листингу
Никакого цикла для нахождения одной точки не нужно. Если хочешь много значений - вбивай их в массив.

К третьему листингу

Цитата

Код


for(int t=0; t<k; t++) z[i,t]=spline(t,Temp,k); //Т.Е. тут мы как раз и получаем значения Z, как мне кажется -\

Если честно - эту строку я тоже не понял. В этом случае массив Temp[] будет не до конца сформирован и ты получишь огромные порешности (или даже ошибки) при малых i.

А теперь ВНИМАНИЕ!!!! Больше повторять не буду.

Для определённости:
spline() - на входе получает
-точку, для которой идёт одномерная интерполяция;
-массивы известных аргументов и соответстующих им значений;
-размер этих массивов;
spline() - возвращает значение для входной точки.

Теперь дописываешь мой код, согласно этим поправкам.

Вот здесь ->

Код


  /*Интерполяция по Y*/
   z0=spline(y0,Temp,K);

вместо z0 можешь подставить и массив.

PSЕдинственное что добавлю, так это то, что я не учёл наличие массива известных аргументов для массива изветных значений (как можно увидеть из "определённостей" я эту багу исправил)

Дальше сам. Удачи!

Это сообщение отредактировал(а) HoTMetaL - 30.6.2007, 20:16

Burka

Дата 30.6.2007, 20:40 (ссылка)

(нет голосов)

Загрузка ...

Опытный

Профиль
Группа: Участник
Сообщений: 429
Регистрация: 11.5.2007

Репутация: 1
Всего: 2

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

Код


Var
 X, Y, Z :TReal1DArray;
 MyCtbl : TReal2DArray;
 i :integer;
[quote]

[/quote]
const N=5;

begin

//Массивы опорных точек. То есть X1 = 2, x7 = 4 и т.д.
x[0]:=1; y[0]:=2;
x[1]:=7; y[1]:=4;
x[2]:=10; y[2]:=6;
x[3]:=17; y[3]:=2;
x[4]:=24; y[4]:=5;

Spline3BuildTable(N, 1, X, Y, 1, 1, MyCtbl);
//Вот описание процедуры из модуля (писал ее не я)
(*************************************************************************
Процедура Spline3BuildTable служит  для  построения таблицы  коэффициентов
кубического сплайна по заданным точкам и граничным условиям, накладываемым
на производные.
В дальнейшем постренная таблица используется функцией Spline3Interpolate.
Параметры:
    N       - число точек
    DiffN   - тип граничного условия. "1" соответствует граничным условиям
              накладываемым на первые производные, "2" - на вторые.
    xs      - массив абсцисс опорных точек с номерами от 0 до N-1
    ys      - массив ординат опорных точек с номерами от 0 до N-1
    BoundL  - левое граничное условие. Если DiffN равно 1, то первая произ-
              водная  на  левой  границе  равна BoundL, иначе BoundL равна
              вторая.
    BoundR  - аналогично BoundL
    
Выходные значения:
    ctbl    - в этот массив помещается таблица коэффициентов сплайна.
              Эта таблица используется функцией Spline3Interpolate.
*************************************************************************)


//Далее сама интерполяция
for i:=1 to 30 do
 begin
  Z[i] := Spline3Interpolate(N, MyCtbl, i);
 end;

Цитата

Один сплошной баг. Ты его вообще запускал? При заполнении StrokaX[] и StrokaY[] ошибок разве не возникает?

Забыл добавить просто:
Inc(Numbers); //в цикле где формируеются массивы опорных точек
И перед циклом
SetLength(StrokaX, 50); //с запасом
SetLength(StrokaY, 50);

Цитата

if B[i,j]<>0 then //Формируем массивы известных точек

Это ещё чё такое?

UPD да, пожалуй тут это не нужно.(спать пора). Надо бы представить двумерный массив опорных точек как-нибудь по-другому. в одномерных как-нибудь.

Цитата

вместо z0 можешь подставить и массив.

Массив я не могу потставить, потому что я не понимаю как можно с помощью "z0=spline(y0,Temp,K);" заполнить двумерный массив с количеством элементов K*F.

Это сообщение отредактировал(а) Burka - 30.6.2007, 22:30

--------------------

Великие умы обсуждают идеи; средние умы обсуждают события; мелкие умы обсуждают людей.

teesync

Дата 27.7.2022, 01:52 (ссылка)

(нет голосов)

Загрузка ...

Новичок

Профиль
Группа: Участник
Сообщений: 36
Регистрация: 26.7.2022

Репутация: нет
Всего: нет

Модератор: Сообщение скрыто.

Queuego

Дата 5.11.2022, 12:46 (ссылка)

(нет голосов)

Загрузка ...

Новичок

Профиль
Группа: Участник
Сообщений: 34
Регистрация: 4.11.2022

Репутация: нет
Всего: нет

Модератор: Сообщение скрыто.

1 Пользователей читают эту тему (1 Гостей и 0 Скрытых Пользователей)
0 Пользователей:
« Предыдущая тема \| Центр помощи \| Следующая тема »