Сплайны

Привет,

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

Да, забыл, естественно толщина ирает большую роль

Спасибо

Сообщение отредактировано: the_moon - 03.07.02, 12:56

Задача сведется просто к обычному сглаживанию функции берегов твоей реки. В разделе уже обсуждалась данная тема.
Посмотри, например
http://pascal.sources.ru/cgi-bin/forum/YaB...;num=1023115703

Это на паскале. Могу и на C.

 
{====================}  UNIT SplinUNT; {======================}
 
{  В этот модуль включены подпрограммы SPLINE и SEVAL из книги
Форсайт Дж., Малькольм М., Моулер К., Машинные методы математи-
ческих вычислений: Пер. с анг.- М.: Мир, 1980. с. 91-94.
 
   Для повышения точности интерполяции используется арифметика
с плавающей точкой, определяемая типом
 
                    RealType = double
 
   Если Ваш компьютер не имеет арифметического сопроцессора, и
Вы хотите несколько увеличить  скорость вычислений ценой  ухуд-
шения точности, а также во всех других случаях, когда необходи-
мо сменить базовый вещественный тип (например, в целях экономии
памяти), Вы должны желаемым образом изменить определение вещес-
твенного типа, например: type RealType = real;
}
{-------------------------------------------------------------}
 
                         INTERFACE
{$N+,E+}
 
type
   RealType = double;              { базовый вещественный тип }
   RealTypeArray =
      array [1..(2*MaxInt) div Sizeof(RealType)] of RealType;
 
procedure Spline (N : word; var XI, YI, BI, CI, DI);
{
   Вычисляются коэффициенты B[i], C[i] и D[i], i = 1..N
   для кубического интерполяционного сплайна вида
 
   S(X) = Y[i]+B[i]*(X-X[i])+C[i]*(X-X[i])**2+D[i]*(X-X[i])**3
 
   для X[i] <= X <= X[i+1].
 
   ВХОДНАЯ ИНФОРМАЦИЯ:
 
   N - число заданных точек или узлов (N >= 2);
   X - абсциссы узлов в строго возрастающем порядке ;
   Y - ординаты узлов.
 
   ВЫХОДНАЯ ИНФОРМАЦИЯ:
 
   B, C, D, - массивы коэффициентов сплайна, записанного в указанной
              выше форме.
 
   Если обозначить дифференцирование символом ', то
   Y[i] = S(X[i]);
   B[i] = S'(X[i]);
   C[i] = S''(X[i])/2;
   D[i] = S'''(X[i])/6  (правосторонняя производная).
 
   С помощью сопровождающей функции SEVAL можно вычислять значения
   сплайна.
}
{-------------------------------------------------------------}
 
function Seval (N : word; var U : RealType;
                var XI, YI, BI, CI, DI) : RealType;
{
   Эта функция вычисляет значение кубического сплайна:
 
   Seval = Y[i]+B[i]*(U-X[i])+C[i]*(U-X[i])**2+D[i]*(U-X[i])**3
 
   где X[i] <= U <= X[i+1]. Используется схема Горнера.
   Если U < X[1], то берется значение i = 1;
   если U > X[N], то i = N.
 
   ВХОДНАЯ ИНФОРМАЦИЯ:
 
   N - число заданных точек (N >= 2);
   U - абсцисса, для которой вычисляется значение сплайна;
   X, Y - массивы заданных абсцисс и ординат;
   B, C, D - массивы коэффициентов сплайна, вычисленные
             процедурой Spline.
 
   Если по сравнению с предыдущим вызовом Seval значение U не
   находится в том же интервале, то для отыскания нужного ин-
   тервала используется двоичный поиск.
}
{-------------------------------------------------------------}
 
                     IMPLEMENTATION
 
procedure Spline (N : word; var XI, YI, BI, CI, DI);
 
var
   X : RealTypeArray absolute XI;
   Y : RealTypeArray absolute YI;
   B : RealTypeArray absolute BI;
   C : RealTypeArray absolute CI;
   D : RealTypeArray absolute DI;
   T : RealType;
   i : integer;
 
begin { Spline }
 
   if N < 2 then Exit;
 
   if N = 2 then begin
 
      B[1] := (Y[2] - Y[1]) / (X[2] - X[1]);
      C[1] := 0;
      D[1] := 0;
      B[2] := B[1];
      C[2] := 0;
      D[2] := 0
 
   end
   else begin
 
    { Формирование трехдиагональной системы линейных уравнений:
      B = диагональ, D = наддиагональ, C = правые части. }
 
      D[1] := X[2] - X[1];
      C[2] := (Y[2] - Y[1]) / D[1];
      for i := 2 to N-1 do begin
         D[i] := X[i+1] - X[i];
         B[i] := 2 * (D[i-1] + D[i]);
         C[i+1] := (Y[i+1] - Y[i]) / D[i];
         C[i] := C[i+1] - C[i]
      end;
 
    { Граничные условия. Третьи производные в точках X[1] и X[N]
      оцениваются с помощью разделенных разностей. }
 
      B[1] := -D[1];
      B[N] := -D[N-1];
      if N = 3 then begin
         C[1] := 0;
         C[N] := 0
      end
      else begin
         C[1] := C[3]/(X[4]-X[2]) - C[2]/(X[3]-X[1]);
         C[N] := C[N-1]/(X[N]-X[N-2]) - C[N-2]/(X[N-1]-X[N-3]);
         C[1] :=  C[1] * sqr(D[1])/(X[4]-X[1]);
         C[N] := -C[N] * sqr(D[N-1])/(X[N]-X[N-3])
      end;
 
    { Прямой ход метода Гаусса }
 
      for i := 2 to N do begin
         T := D[i-1]/B[i-1];
         B[i] := B[i] - T*D[i-1];
         C[i] := C[i] - T*C[i-1]
      end;
 
    { Обратная подстановка }
 
      C[N] := C[N]/B[N];
      for i := N-1 downto 1 do C[i] := (C[i] - D[i]*C[i+1]) / B[i];
 
    { В C[i] теперь хранится значение Sigma(i), определяемое в книге:
      Дж.Форсайт, М.Малькольм, К.Моулер, Машинные методы математических
      вычислений. }
 
    { Вычисление коэффициентов сплайна B, C и D }
 
      B[N] := (Y[N]-Y[N-1])/D[N-1] + D[N-1]*(C[N-1]+2*C[N]);
      for i := 1 to N-1 do begin
         B[i] := (Y[i+1]-Y[i])/D[i] - D[i]*(C[i+1]+2*C[i]);
         D[i] := (C[i+1]-C[i])/D[i];
         C[i] := 3*C[i]
      end;
      C[N] := 3*C[N];
      D[N] := D[N-1]
 
   end { N <> 2 };
 
end { Spline };
 
{-------------------------------------------------------------}
 
function Seval (N : word; var U : RealType;
                var XI, YI, BI, CI, DI) : RealType;
 
var
   X : RealTypeArray absolute XI;
   Y : RealTypeArray absolute YI;
   B : RealTypeArray absolute BI;
   C : RealTypeArray absolute CI;
   D : RealTypeArray absolute DI;
   DX : RealType;
   j, k : integer;
 
const
   i : integer = 1;
 
begin { Seval }
 
   if i >= N then i := 1;
 
   if (U < X[i]) or (U > X[i+1]) then begin { Двоичный поиск }
      i := 1;
      j := N+1;
      repeat
         k := (i+j) div 2;
         if U < X[k] then j := k else i := k
      until j <= i+1
   end;
 
 { Вычисление сплайна }
 
   DX := U - X[i];
   Seval := Y[i] + DX*(B[i] + DX*(C[i] + DX*D[i]))
 
end { Seval };
 
{=================} END. { Unit SplinUNT } {==================}

Да! да! на Си пожалуйста.

 
/* spline.h
*/
#ifndef __SPLINE_H
 
double *spline (int n, double *x, double *y, double *c);
double *splcvt (int n, double *x, double *y, double *c);
double  seval  (int n, double z, double *x, double *y, double *c);
double  seval1 (int n, double z, double *x, double *y, double *c);
 
#define __SPLINE_H
#endif /* __SPLINE_H */

 
/* spline.c
*/
#include <float.h>
#include <stddef.h>
 
#include "spline.h"
 
#define TOL        (2.0*DBL_MIN/DBL_EPSILON)
#define square(x)  ((x)*(x))
 
 
double *spline (int n, double x[], double y[], double c[])
/*
 * Вычисление коэффициентов кубического интерполяционного сплайна,
 * записанного в виде:
 *
 * S(x) = w(i)*y(i+1) + w'(i)*y(i) +
 *        h(i)^2 * ((w(i)^3 - w(i))*c(i+1) + (w'(i)^3 - w'(i))*c(i)),
 *
 * где h(i)  = x(i+1) - x(i),
 *     w(i)  = (x - x(i)) / h(i),
 *     w'(i) = 1 - w(i).
 *
 * Аргументы:
 *
 * n - число заданных точек (узлов) (n >= 2),
 * x - массив абсцисс узлов в строго возрастающем порядке,
 * y - массив ординат узлов,
 * c - рабочий массив размером не менее 3*n элементов типа double.
 *     После возврата из функции spline в первых n элементах этого
 *     массива находятся коэффициенты сплайна c(i).
 *
 * Возвращаемое значение:
 *
 * Указатель на массив вычисленных коэффициентов c[] либо NULL, если
 * построить сплайн невозможно (n < 2, абсциссы узлов не упорядочены
 * либо среди них имеются совпадающие значения).
*/
{
   double *b = c + n;
   double *d = b + n;
 
   double t;
   int i;
 
   if (n < 2)
 
      return NULL;
 
   else if (n == 2)
 
      c[1] = c[0] = 0.0;
 
   else  {
      /* Формирование трехдиагональной системы линейных уравнений.
         b = диагональ, d = наддиагональ, c = правые части.
      */
      for (i = 0; i < n-1; i++) {
         if ((d[i] = x[i+1] - x[i]) <= TOL) return (NULL);
         c[i+1] = (y[i+1] - y[i]) / d[i];
         if (i) {
            c[i] = c[i+1] - c[i];
            b[i] = 2.0 * (d[i-1] + d[i]);
         }
      }
      /* Граничные условия. Третьи производные в точках x[0] и x[n-1]
         вычисляются с помощью разделенных разностей.
      */
      b[0] = -d[0];
      b[n-1] = -d[n-2];
      if (n == 3)
         c[0] = c[n-1] = 0.0;
      else {
         c[0]   = c[2]/(x[3]-x[1]) - c[1]/(x[2]-x[0]);
         c[n-1] = c[n-2]/(x[n-1]-x[n-3]) - c[n-3]/(x[n-2]-x[n-4]);
         c[0]  *= square(d[0]) / (x[3]-x[0]);
         c[n-1] = -(c[n-1]*square(d[n-2]) / (x[n-1]-x[n-4]));
      }
 
      /* Прямой ход метода Гаусса */
 
      for (i = 0; i < n-1; i++) {
         t = d[i] / b[i];
         b[i+1] -= t*d[i];
         c[i+1] -= t*c[i];
      }
 
      /* Обратная подстановка */
 
      c[n-1] /= b[n-1];
      for (i = n-2; i >= 0; i--)
         c[i] = (c[i] - d[i]*c[i+1]) / b[i];
   }
   return c;
}

 
/* splcvt.c
*/
#include <stddef.h>
 
#include "spline.h"
 
double *splcvt (int n, double x[], double y[], double c[])
/*
 * Преобразование сплайна, построенного с помощью функции spline,
 * к виду
 *
 * S(x) = y[i] + b[i]*(x-x[i]) + c[i]*(x-x[i])^2 + d[i]*(x-x[i])^3
 *
 * где x[i] <= x <= x[i+1].
 *
 * Аргументы:
 *
 * n - число заданных точек (узлов) (n >= 2),
 * x - массив абсцисс узлов в строго возрастающем порядке,
 * y - массив ординат узлов,
 * c - массив размером не менее 3*n элементов типа double.
 *     В первых n элементах этого массива должны находиться коэффициенты
 *     сплайна, рассчитанные функцией spline. После вызова splcvt массив
 *     содержит результат преобразования, т.е. коэффициенты b[i], c[i], d[i]
 *     в следующем порядке:
 *     - первые n элементов массива: коэффициенты c[i],
 *     - следующие n элементов:      коэффициенты b[i],
 *     - последние n элементов:      коэффициенты d[i].
 *
 * Возвращаемое значение:
 *
 * Указатель на массив вычисленных коэффициентов c[] либо NULL, если
 * построить сплайн невозможно (n < 2).
*/
{
   double *b = c + n;
   double *d = b + n;
   int i;
 
   if (n < 2)
 
      return NULL;
 
   else if (n == 2) {
 
      b[1] = b[0] = (y[1] - y[0]) / (x[1] - x[0]);
      d[1] = d[0] = 0.0;
 
   }
   else {
 
      for (i = 0; i < n-1; i++) {
         d[i] = x[i+1] - x[i];
         b[i] = (y[i+1] - y[i]) / d[i];
      }
 
      b[n-1] = b[n-2] + d[n-2] * (c[n-2] + 2.0*c[n-1]);
 
      for (i = 0; i < n-1; i++) {
         b[i] -= d[i] * (c[i+1] + 2.0*c[i]);
         d[i] = (c[i+1] - c[i]) / d[i];
         c[i] *= 3.0;
      }
      d[n-1] = d[n-2];
      c[n-1] *= 3.0;
 
   }
   return c;
}

 
/* seval1.c
*/
#include "spline.h"
 
double seval1 (int n, double z, double x[], double y[], double c[])
/*
 * Значение кубического сплайна в заданной точке.
 *
 * S(z) = y(i) + b(i)*(z-x(i)) + c(i)*(z-x(i))^2 + d(i)*(z-x(i))^3
 *
 * Номер интервала i выбирается таким образом, чтобы для значений z,
 * лежащих внутри диапазона табличных абсцисс, выполнялось условие:
 *
 *     x(i) <= z <= x(i+1).
 *
 * Если z < x(0), берется i = 0; если z >= x(n-1), берется i = n-1.
 *
 * Если по сравнению с предыдущим вызовом seval значение z не находится
 * в том же интервале [x(i), x(i+1)], тo для отыскания нужного интервала
 * применяется двоичный поиск.
 *
 * Аргументы:
 *
 * n - число заданных точек (узлов) (n >= 2),
 * z - абсцисса, для которой вычисляется значение сплайна,
 * x - массив абсцисс узлов в строго возрастающем порядке,
 * y - массив ординат узлов,
 * c - массив коэффициентов сплайна, рассчитанных функцией splcvt.
*/
{
   static int i = 0;
 
   double *b = c + n;
   double *d = b + n;
   double dx;
   int j, k;
 
   if (i >= n-1) i = 0;
 
   /* Двоичный поиск */
 
   if (z < x[i] || z > x[i+1]) {
      i = 0;
      j = n;
      while (j > i+1) {
         k = (i+j)/2;
         if (z < x[k])
            j = k;
         else
            i = k;
      }
   }
 
   /* Вычисление сплайна */
 
   dx = z - x[i];
   return (y[i] + dx*(b[i] + dx*(c[i] + dx*d[i])));
}

 
/* seval.c
*/
#include "spline.h"
 
#define square(x)  ((x)*(x))
 
double seval (int n, double z, double x[], double y[], double c[])
/*
 * Значение кубического сплайна в заданной точке.
 *
 * S(z) = w(i)*y(i+1) + w'(i)*y(i) +
 *        h(i)^2 * ((w(i)^3 - w(i))*c(i+1) + (w'(i)^3 - w'(i))*c(i)),
 *
 * где h(i)  = x(i+1) - x(i),
 *     w(i)  = (z - x(i)) / h(i),
 *     w'(i) = 1 - w(i).
 *
 * Номер интервала i выбирается таким образом, чтобы для значений z,
 * лежащих внутри диапазона табличных абсцисс, выполнялось условие:
 *
 *     x(i) <= z <= x(i+1).
 *
 * Если z < x(0), берется i = 0; если z > x(n-1), берется i = n-2.
 *
 * Если по сравнению с предыдущим вызовом seval значение z не находится
 * в том же интервале [x(i), x(i+1)], тo для отыскания нужного интервала
 * применяется двоичный поиск.
 *
 * Аргументы:
 *
 * n - число заданных точек (узлов) (n >= 2),
 * z - абсцисса, для которой вычисляется значение сплайна,
 * x - массив абсцисс узлов в строго возрастающем порядке,
 * y - массив ординат узлов,
 * c - массив коэффициентов сплайна, рассчитанных функцией spline.
*/
{
   static int i = 0;
 
   double h, w, w1;
   int j, k;
 
   if (i >= n-1) i = 0;
 
   /* Двоичный поиск */
 
   if (z < x[i] || z > x[i+1]) {
      i = 0;
      j = n-1;
      while (j > i+1) {
         k = (i+j)/2;
         if (z < x[k])
            j = k;
         else
            i = k;
      }
   }
 
   /* Вычисление сплайна */
 
   h = x[i+1] - x[i];
   w1 = 1.0 - (w = (z-x[i])/h);
   return (w*y[i+1] + w1*y[i] +
           square(h) * (w*(square(w)-1.0)*c[i+1] + w1*(square(w1)-1.0)*c[i]));
}

а можно то же самое (на си) для функции 3-х переменных?
заранее благодарен....