Кривые безье, сплайны...

вашу мать..чё никто не знает ? тут блин, что, совсем никто-никто ????????

мне стыдно..мне в натуре больно и обидно, что НИКТО не может провести кривую через N точек...Вы вообще где ???? Тут вообще кто-нить есть ??

никто не знал полиномов ?

куды мир катится....

С кривыми Безье не сталкивался.

Если нужно провести гладкую кривую(не обязательно Безье) могу дать исходники проги интерполятора полиномом/кубическими сплайнами.

2MKuznetsov:

Пъян что ли :-) ?

Разумеется, тут есть перцЫ. Но прямо так ответить на вопрос...

2drF_ckoff:

На секунду забудь умные слова.

У тебя есть n штук точек, Тебе надо провести через них кривую (я почти уверен, что тебе не это надо, BTW),

Вспоминаем пятый класс.

y=a*x^n+b*x^(n-1)+...+c

Подставлям, приравнивам... Если не помогает, лечимся эвтаназией (больше уже ничто не поможет)...

> С кривыми Безье не сталкивался.

Откровение, честно говоря... IMHO слово "сплайн" нынче менее популярно, нежели "кривая Безье"...

Фаил function.h

// function.h: interface for the function class. // ////////////////////////////////////////////////////////////////////// #include <iostream> #include <stdlib.h>

#if !defined(AFX_FUNCTION_H__AC2EA9A3_327E_11DA_BC7C_009027A7118F__INCLUDED_) #define AFX_FUNCTION_H__AC2EA9A3_327E_11DA_BC7C_009027A7118F__INCLUDED_

#if _MSC_VER > 1000 #pragma once #endif // _MSC_VER > 1000

inline double max(long double a,long double b) { if(a>b) return a; return b; }

inline double min(long double a,long double b) { if(a<b) return a; return b; }

inline double abs_v(long double a) { if(a>0) return a; return -a; }

class function { public: virtual long double operator()(long double) = 0; virtual function &diff() = 0; friend std::ostream& operator<<(std::ostream&,function&); virtual void print() = 0; };

#endif // !defined(AFX_FUNCTION_H__AC2EA9A3_327E_11DA_BC7C_009027A7118F__INCLUDED_)

// polynome.h: interface for the polynome class. // //////////////////////////////////////////////////////////////////////

#include "function.h" #include <list> //#include <xutility> #include <utility> #include <iterator>

#if !defined(AFX_POLYNOME_H__AC2EA9A9_327E_11DA_BC7C_009027A7118F__INCLUDED_) #define AFX_POLYNOME_H__AC2EA9A9_327E_11DA_BC7C_009027A7118F__INCLUDED_

static long double nocoeff = 0;

typedef std::pair<long double,long double> interval; typedef std::list<interval> intervals; #define points std::list<long double>

class polynome : public function { int power; long double* coeff; public: polynome(int); virtual ~polynome(); inline long double &operator[](int i){ if((i > power)||(i < 0))return nocoeff; return coeff[i]; }; long double operator()(long double); virtual function &diff(); inline int get_power(){return power;}; virtual void print();

long double* begin() {return coeff;} long double* end() {return coeff + power + 1;} polynome &operator+=(polynome&); polynome &operator+=(long double &c); polynome &operator*=(polynome&); polynome &operator*=(long double); typedef long double* iterator; };

points solve(polynome&,double);

#endif // !defined(AFX_POLYNOME_H__AC2EA9A9_327E_11DA_BC7C_009027A7118F__INCLUDED_)

фаил polynome.cpp

#include "polynome.h"

polynome::polynome(int pow = 0){ power = pow; coeff = (long double*)malloc(sizeof(long double)*(power+1)); for(int i = 0;i <= power;i++)coeff[i] = 0; }

polynome::~polynome(){ free(coeff); }

long double polynome::operator()(long double x){ long double summ = 0; long double cur = 1; for(int i = power;i >= 0;i--){ summ += cur * coeff[i]; cur *= x; } return summ; }

function &polynome::diff(){ polynome* dif = NULL; int i; i = power - 1; if(i < 0)i = 0; dif = new polynome(i); for(int j = 0;j < i + 1;j++){ (*dif)[j] = ((*this)[j])*(i + 1 - j); } return *dif; }

void polynome::print(){ for(int i = 0;i < power;i++) std::cout << (*this)[i] <<"*x^" << power - i << " + "; std::cout << (*this)[power]; }

std::ostream& operator<<(std::ostream& os, const interval& p) { os << '[' << p.first << ','<< p.second << ']'; return os; }

polynome &polynome::operator+=(polynome &p){ int maxpow = (int)max(power,p.power); int i; long double *ad = (long double *)malloc(sizeof(long double)*(1 + maxpow)); for(i = 0;i <=maxpow;i++) ad[i] = 0; for(i = 0;i <= p.power;i++) ad[i + maxpow - p.power] += p[i]; for(i = 0;i <= power;i++) ad[i + maxpow - power] += (*this)[i]; power = maxpow; free(coeff); coeff = ad; return (*this); }

polynome &polynome::operator+=(long double &c){ (*this)[power] += c; return (*this); }

polynome &polynome::operator*=(polynome &p){ long double *ad = (long double *)malloc(sizeof(long double)*(power + p.power + 2)); int i; int j; memset(ad,0,sizeof(long double)*(power + p.power + 2)); for(i = 0; i <= power;i++) for(j = 0; j <= p.power;j++) ad[i+j] += p[j] * (*this)[i]; free(coeff); coeff = ad; power = power + p.power; return (*this); }

polynome &polynome::operator*=(long double c){ for(int i = 0;i <= power;i++) (*this)[i] *= c; return (*this); }

фаил interpolator.h

#include "polynome.h" #include <vector>

#if !defined(AFX_INTERPOLATOR_H__D1A1A9C7_5A0E_41FC_962B_EB293E0BF9D0__INCLUDED_) #define AFX_INTERPOLATOR_H__D1A1A9C7_5A0E_41FC_962B_EB293E0BF9D0__INCLUDED_

class interpolator{ protected: long double *x; long double *f; public: polynome *P; };

class newthon_interpolator: public interpolator{ long double *b; inline long double calc_b(int n,int min,int max); public: newthon_interpolator(int n,long double *x_orig,long double *f_orig); virtual ~newthon_interpolator(); };

class lagranj_interpolator: public interpolator{ inline long double calc_k(int n); int m; public: lagranj_interpolator(int n,long double *x_orig,long double *f_orig); virtual ~lagranj_interpolator(); };

class qube_spliner: public function{ std::vector<polynome*> aprox; long double *x; long double *f; int m; long double *a; public: virtual long double operator()(long double); virtual function &diff(); friend std::ostream& operator<<(std::ostream&,qube_spliner&); virtual void print(); qube_spliner(interpolator *,int ,long double *,long double *); virtual ~qube_spliner(); };

#endif

фаил interpolator.cpp

#include "interpolator.h"

inline long double lagranj_interpolator::calc_k(int n){ long double a = 1; int i; for(i = 0;i < m;i++) if(i != n)a *= x[n] - x[i]; return f[n] / a; }

lagranj_interpolator::lagranj_interpolator(int n,long double *x_orig,long double *f_orig){ x = (long double*)malloc(n * sizeof(long double)); f = (long double*)malloc(n * sizeof(long double)); int i; int j; m = n; long double tmp; P = new polynome(n - 1); memcpy(x,x_orig,n * sizeof(long double)); memcpy(f,f_orig,n * sizeof(long double)); for(i = 0;i < m;i++){ polynome g(0); g[0] = 1; for(j = 0;j < m;j++){ if(i != j){ polynome h(1); h[0] = 1; h[1] = -1 * x[j]; g *= h; } } g *= calc_k(i); *P += g; } }

lagranj_interpolator::~lagranj_interpolator(){ free(x); free(f); delete P; }

inline long double newthon_interpolator::calc_b(int n,int min,int max){ if(n == 0)return f[min]; return (calc_b(n - 1, min + 1, max) - calc_b(n - 1, min, max - 1))/(x[max] - x[min]); }

newthon_interpolator::newthon_interpolator(int n,long double *x_orig,long double *f_orig){ x = (long double*)malloc(n * sizeof(long double)); f = (long double*)malloc(n * sizeof(long double)); b = (long double*)malloc(n * sizeof(long double)); int i; int j; P = new polynome(n - 1); memcpy(x,x_orig,n * sizeof(long double)); memcpy(f,f_orig,n * sizeof(long double)); for(i = 0; i < n; i++){ b[i] = calc_b(i,0,i); } for(i = 1; i <= n - 1; i++){ polynome g(1); g[0] = 1; g[1] = -1 * x[0]; for(j = 1;j < i;j++){ polynome h(1); h[0] = 1; h[1] = -1 * x[j]; g *= h; } g *= b[i]; *P += g; } *P += b[0]; }

newthon_interpolator::~newthon_interpolator(){ free(x); free(f); free(b); delete P; }

qube_spliner::qube_spliner(interpolator* j,int n,long double *x_orig,long double *f_orig){ polynome *P1; P1 = (polynome*)(&(j->P->diff())); int i; x = (long double*)malloc(n * sizeof(long double)); f = (long double*)malloc(n * sizeof(long double)); a = (long double*)malloc(4 * (n - 1) * sizeof(long double)); memcpy(x,x_orig,n * sizeof(long double)); memcpy(f,f_orig,n * sizeof(long double)); m = n; for(i = 0; i < n - 1; i++){ polynome *APROX = new polynome(3); int k; long double dx = x[i+1]-x[i]; long double df = f[i+1]-f[i]; long double dx3 = dx*dx*dx; long double p1 = (*P1)(x[i]); long double p2 = (*P1)(x[i+1]);

a[4 * i + 0] = ((-1)*p2*x[i]*x[i]*x[i+1]*dx + (f[i+1]*x[i]*x[i]*(3*x[i + 1] - x[i])) + (f[i]*x[i+1]*x[i+1]*(x[i+1] - 3*x[i])) - (p1*x[i]*x[i+1]*x[i+1]*dx)) / dx3;

a[4 * i + 1] = (p2*x[i]*(2*x[i+1] + x[i])*dx - 6*df*x[i]*x[i+1] + p1*x[i+1]*(x[i+1]+2*x[i])*dx)/dx3;

a[4 * i + 2] = ((-1)*p2*dx*(x[i+1] + 2 * x[i]) + 3*df*(x[i+1]+x[i]) - p1*(x[i+1] - x[i])*(x[i]+2*x[i+1]))/dx3;

a[4 * i + 3] = (p2*dx - 2*df + p1*dx)/ dx3;

for(k = 0;k < 4;k++) (*APROX)[k] = a[4 * i + 3 - k]; aprox.push_back(APROX); } free(a); delete P1; }

long double qube_spliner::operator()(long double xx){ for(int i = 0; i < m - 1;i++) if(x[i] <= xx && x[i+1] >= xx)return (*aprox[i])(xx); return 0; }

function &qube_spliner::diff(){ qube_spliner* dif = NULL; //TODO: diff return *dif; }

void qube_spliner::print(){ for(int i = 0;i < m - 1;i++) std::cout << x[i] << " < X < " << x[i+1] << ": F(x) = " << (*aprox[i]) << std::endl; }

std::ostream& operator<<(std::ostream& os,qube_spliner& p) { p.print(); return os; }

qube_spliner::~qube_spliner(){ int i; for(i = 0;i < aprox.size();i++) delete aprox[i]; free(x); free(f); }

применение этого кода

#include "interpolator.h" #include <math.h> #include <vector>

using namespace std;

int main(){ #define MAX 6 long double x[MAX] = {0.52360, 0.87267, 1.22173, 1.57080, 1.91986, 2.26893}; long double f[MAX] = {0.000014, 0.00037, 0.00352, 0.01971, 0.08033, 0.26380};

int i;

for(i = 0;i < MAX;i++) cout << "X["<< i << "] = " << x[i] << endl; cout << endl; for(i = 0;i < MAX;i++) cout << "F["<< i << "] = " << f[i] << endl; cout << endl;

newthon_interpolator inter(MAX,x,f);

cout << "Interpolation: " << *(inter.P) << endl; for(i = 0;i < MAX;i++) cout << "F(x[" << i << "]) = " << (*(inter.P))(x[i]) << endl;

qube_spliner spl(&inter,MAX,x,f);

cout << endl << spl << endl;

/* long double zzz;

#define SPLIT 100

for(i = 0;i <= SPLIT;i++){ zzz = x[0] + i*(x[MAX - 1] - x[0])/SPLIT; cout << "X = " << zzz << " Int = "<< (*(inter.P))(zzz) << " Spline = "<< spl(zzz) << endl; }*/

return 0; }

вроде усё 8))

извиняюсь за флуд.

YesSSS:

> извиняюсь за флуд.

Должен я все это убрать? :-)

(или сам понимаешь, что вряд ли тебе по существу ответят)

2Die-hard: полином n-ой степени может выглядеть "хреново", мягко говоря. Вопрос еще в том, что значит "провести". Восстановить значение, между точками - называется интерполяция. У сплайнов фишка в том, что разные точки соединяешь отдельными полиномами, но с такими коэффициентами, чтобы суммарная кривая была дважды гладкой.

anonymous (*) (27.11.2005 18:06:47):

> 2Die-hard: полином n-ой степени может выглядеть "хреново", мягко говоря.

Разумеется!

Вопрос, для чего нужна интерполяция. Если для интегрирования, то обычно достаточно иметь просто негладкие кубические сплайны.

Есть целая наука про сплайны, но все идеи там довольно примитивны. Основное достижение мысли в этой области, про дважды гладкую кривую, ты совершенно правильно озвучил, но это и есть почти кривые Безье (те чуть послабже, дважды непрерывны). Все всегда сводится к сопряжению полиномов, проходящих через n>=2 точек, со сшивкой производных достаточно высокого порядка.

Проблемы начинаются тогда, когда тебе нужна достаточно многократная гладкость -- тогда полиномы обычно пасуют. Но тогда и интерполяция обычно пасует (спасает близкая аппроксимация аналитическими фукциями), что уже достаточно нетривиально.

Вобщем, в приведенной постановке вопрос не совсем корректен. Для чего нужно сглаживание?

2drF_ckoff:

Извиняюсь за свои предыдущие посты, задача поставлена вполне здраво: как провести интерполяцию через n>2 точек кривыми Безье.

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

Это мне надо, это. "Подставляем, приравниваем". Перед этим не забываем "повернуть" систему координат (предварительно вспомнив геометрию, которую я не помню напрочь, как написал в disclaimer'е). Потом все "линиями" рисуем. И всё это на питоне. И будет оно очень быстро.

btw, одна из причин, почему я этот вопрос задаю - это то, что со сплайнами я, худо-бедно, знаком. До них я доучился. А вот что за хрень эти "кривые безье" понимаю с трудом. =)

drF_ckoff:

> ...со сплайнами я, худо-бедно, знаком. ... А вот что за хрень эти "кривые безье" понимаю с трудом. =)

Не могет такого быть, ибо Б-сплайны основаны на Безье.

Тебе нужны не просто Б-сплайны, а "итерполяционные" Б-сплайны, можешь погуглить, наверняка куча вывалится (interpolating B-spline).

Если проблема осталась, открой новый топик в этом форуме -- обещаю час потратить на гугление; мне, скорее всего, скоро самому пригодится.