Estoy buscando una buena biblioteca de C++ que me dé funciones para resolver grandes splines cúbicos (del orden de 1000 puntos) ¿Alguien conoce uno?¿Existen buenas bibliotecas para resolver splines cúbicos en C++?
Respuesta
Prueba el cúbica biblioteca B-spline:
y ALGLIB:
El segundo enlace era más o menos lo que estaba buscando ... excepto que no puedo encontrar ninguna documentación en ese sitio web o en el archivo descargado ... – Faken
La mayoría de la documentación está indexada aquí: http://www.alglib.net/sitemap.php - No creo que haya ninguna que se pueda descargar (se ha planificado por un tiempo, no creo que haya sucedido todavía). – ars
Parece que ALGLIB tiene algunas restricciones: [No admite la interpolación 3D] (http://forum.alglib.net/viewtopic.php?f=2&t=137); fuerza las splines periódicas; [Los valores de las funciones deben estar en una cuadrícula regular] (http://forum.alglib.net/viewtopic.php?f=2&t=137). – Jeff
Tome un vistazo a David Eberly de GeometricTools.com . Estoy empezando, pero el código y el documento son hasta ahora de una calidad excepcional.
(También tiene libros: herramientas geométricas para gráficos por computadora, diseño de motores de juegos en 3D.)
Escribe el tuyo. Aquí es spline()
función que escribí gracias al excelente wiki algorithm:
#include<iostream>
#include<vector>
#include<algorithm>
#include<cmath>
using namespace std;
typedef vector<double> vec;
struct SplineSet{
double a;
double b;
double c;
double d;
double x;
};
vector<SplineSet> spline(vec &x, vec &y)
{
int n = x.size()-1;
vec a;
a.insert(a.begin(), y.begin(), y.end());
vec b(n);
vec d(n);
vec h;
for(int i = 0; i < n; ++i)
h.push_back(x[i+1]-x[i]);
vec alpha;
for(int i = 0; i < n; ++i)
alpha.push_back(3*(a[i+1]-a[i])/h[i] - 3*(a[i]-a[i-1])/h[i-1] );
vec c(n+1);
vec l(n+1);
vec mu(n+1);
vec z(n+1);
l[0] = 1;
mu[0] = 0;
z[0] = 0;
for(int i = 1; i < n; ++i)
{
l[i] = 2 *(x[i+1]-x[i-1])-h[i-1]*mu[i-1];
mu[i] = h[i]/l[i];
z[i] = (alpha[i]-h[i-1]*z[i-1])/l[i];
}
l[n] = 1;
z[n] = 0;
c[n] = 0;
for(int j = n-1; j >= 0; --j)
{
c[j] = z [j] - mu[j] * c[j+1];
b[j] = (a[j+1]-a[j])/h[j]-h[j]*(c[j+1]+2*c[j])/3;
d[j] = (c[j+1]-c[j])/3/h[j];
}
vector<SplineSet> output_set(n);
for(int i = 0; i < n; ++i)
{
output_set[i].a = a[i];
output_set[i].b = b[i];
output_set[i].c = c[i];
output_set[i].d = d[i];
output_set[i].x = x[i];
}
return output_set;
}
int main()
{
vec x(11);
vec y(11);
for(int i = 0; i < x.size(); ++i)
{
x[i] = i;
y[i] = sin(i);
}
vector<SplineSet> cs = spline(x, y);
for(int i = 0; i < cs.size(); ++i)
cout << cs[i].d << "\t" << cs[i].c << "\t" << cs[i].b << "\t" << cs[i].a << endl;
}
'for (int i = 0; ...' y luego accede a 'a [i-1]' y 'h [i-1]' – helium
Gracias por este código. Muy útil! Upvoted. Falta el código para dibujar las curvas spline y hacer algunas comprobaciones de error. He agregado esto y lo he envuelto todo en una clase. https://github.com/JamesBremner/raven-cSpline Déjame saber si desea cambiar su atribución (p. ej., nombre real). – ravenspoint
que tenía que escribir la rutina spline de una "entidad" que estaba siguiendo un camino (serie de waypoints conectados) en un juego que estoy trabajando.
Creé una clase base para manejar un "SplineInterface" y las dos clases derivadas creadas, una basada en la técnica de spline clásica (por ejemplo, Sedgewick/Algorithms) y una segunda basada en Bezier Splines.
Aquí está el código. Es un archivo de cabecera única, que contiene todas las clases splining:
#ifndef __SplineCommon__
#define __SplineCommon__
#include "CommonSTL.h"
#include "CommonProject.h"
#include "MathUtilities.h"
/* A Spline base class. */
class SplineBase
{
private:
vector<Vec2> _points;
bool _elimColinearPoints;
protected:
protected:
/* OVERRIDE THESE FUNCTIONS */
virtual void ResetDerived() = 0;
enum
{
NOM_SIZE = 32,
};
public:
SplineBase()
{
_points.reserve(NOM_SIZE);
_elimColinearPoints = true;
}
const vector<Vec2>& GetPoints() { return _points; }
bool GetElimColinearPoints() { return _elimColinearPoints; }
void SetElimColinearPoints(bool elim) { _elimColinearPoints = elim; }
/* OVERRIDE THESE FUNCTIONS */
virtual Vec2 Eval(int seg, double t) = 0;
virtual bool ComputeSpline() = 0;
virtual void DumpDerived() {}
/* Clear out all the data.
*/
void Reset()
{
_points.clear();
ResetDerived();
}
void AddPoint(const Vec2& pt)
{
// If this new point is colinear with the two previous points,
// pop off the last point and add this one instead.
if(_elimColinearPoints && _points.size() > 2)
{
int N = _points.size()-1;
Vec2 p0 = _points[N-1] - _points[N-2];
Vec2 p1 = _points[N] - _points[N-1];
Vec2 p2 = pt - _points[N];
// We test for colinearity by comparing the slopes
// of the two lines. If the slopes are the same,
// we assume colinearity.
float32 delta = (p2.y-p1.y)*(p1.x-p0.x)-(p1.y-p0.y)*(p2.x-p1.x);
if(MathUtilities::IsNearZero(delta))
{
_points.pop_back();
}
}
_points.push_back(pt);
}
void Dump(int segments = 5)
{
assert(segments > 1);
cout << "Original Points (" << _points.size() << ")" << endl;
cout << "-----------------------------" << endl;
for(int idx = 0; idx < _points.size(); ++idx)
{
cout << "[" << idx << "]" << " " << _points[idx] << endl;
}
cout << "-----------------------------" << endl;
DumpDerived();
cout << "-----------------------------" << endl;
cout << "Evaluating Spline at " << segments << " points." << endl;
for(int idx = 0; idx < _points.size()-1; idx++)
{
cout << "---------- " << "From " << _points[idx] << " to " << _points[idx+1] << "." << endl;
for(int tIdx = 0; tIdx < segments+1; ++tIdx)
{
double t = tIdx*1.0/segments;
cout << "[" << tIdx << "]" << " ";
cout << "[" << t*100 << "%]" << " ";
cout << " --> " << Eval(idx,t);
cout << endl;
}
}
}
};
class ClassicSpline : public SplineBase
{
private:
/* The system of linear equations found by solving
* for the 3 order spline polynomial is given by:
* A*x = b. The "x" is represented by _xCol and the
* "b" is represented by _bCol in the code.
*
* The "A" is formulated with diagonal elements (_diagElems) and
* symmetric off-diagonal elements (_offDiagElemns). The
* general structure (for six points) looks like:
*
*
* | d1 u1 0 0 0 | | p1 | | w1 |
* | u1 d2 u2 0 0 | | p2 | | w2 |
* | 0 u2 d3 u3 0 | * | p3 | = | w3 |
* | 0 0 u3 d4 u4 | | p4 | | w4 |
* | 0 0 0 u4 d5 | | p5 | | w5 |
*
*
* The general derivation for this can be found
* in Robert Sedgewick's "Algorithms in C++".
*
*/
vector<double> _xCol;
vector<double> _bCol;
vector<double> _diagElems;
vector<double> _offDiagElems;
public:
ClassicSpline()
{
_xCol.reserve(NOM_SIZE);
_bCol.reserve(NOM_SIZE);
_diagElems.reserve(NOM_SIZE);
_offDiagElems.reserve(NOM_SIZE);
}
/* Evaluate the spline for the ith segment
* for parameter. The value of parameter t must
* be between 0 and 1.
*/
inline virtual Vec2 Eval(int seg, double t)
{
const vector<Vec2>& points = GetPoints();
assert(t >= 0);
assert(t <= 1.0);
assert(seg >= 0);
assert(seg < (points.size()-1));
const double ONE_OVER_SIX = 1.0/6.0;
double oneMinust = 1.0 - t;
double t3Minust = t*t*t-t;
double oneMinust3minust = oneMinust*oneMinust*oneMinust-oneMinust;
double deltaX = points[seg+1].x - points[seg].x;
double yValue = t * points[seg + 1].y +
oneMinust*points[seg].y +
ONE_OVER_SIX*deltaX*deltaX*(t3Minust*_xCol[seg+1] - oneMinust3minust*_xCol[seg]);
double xValue = t*(points[seg+1].x-points[seg].x) + points[seg].x;
return Vec2(xValue,yValue);
}
/* Clear out all the data.
*/
virtual void ResetDerived()
{
_diagElems.clear();
_bCol.clear();
_xCol.clear();
_offDiagElems.clear();
}
virtual bool ComputeSpline()
{
const vector<Vec2>& p = GetPoints();
_bCol.resize(p.size());
_xCol.resize(p.size());
_diagElems.resize(p.size());
for(int idx = 1; idx < p.size(); ++idx)
{
_diagElems[idx] = 2*(p[idx+1].x-p[idx-1].x);
}
for(int idx = 0; idx < p.size(); ++idx)
{
_offDiagElems[idx] = p[idx+1].x - p[idx].x;
}
for(int idx = 1; idx < p.size(); ++idx)
{
_bCol[idx] = 6.0*((p[idx+1].y-p[idx].y)/_offDiagElems[idx] -
(p[idx].y-p[idx-1].y)/_offDiagElems[idx-1]);
}
_xCol[0] = 0.0;
_xCol[p.size()-1] = 0.0;
for(int idx = 1; idx < p.size()-1; ++idx)
{
_bCol[idx+1] = _bCol[idx+1] - _bCol[idx]*_offDiagElems[idx]/_diagElems[idx];
_diagElems[idx+1] = _diagElems[idx+1] - _offDiagElems[idx]*_offDiagElems[idx]/_diagElems[idx];
}
for(int idx = (int)p.size()-2; idx > 0; --idx)
{
_xCol[idx] = (_bCol[idx] - _offDiagElems[idx]*_xCol[idx+1])/_diagElems[idx];
}
return true;
}
};
/* Bezier Spline Implementation
* Based on this article:
* http://www.particleincell.com/blog/2012/bezier-splines/
*/
class BezierSpine : public SplineBase
{
private:
vector<Vec2> _p1Points;
vector<Vec2> _p2Points;
public:
BezierSpine()
{
_p1Points.reserve(NOM_SIZE);
_p2Points.reserve(NOM_SIZE);
}
/* Evaluate the spline for the ith segment
* for parameter. The value of parameter t must
* be between 0 and 1.
*/
inline virtual Vec2 Eval(int seg, double t)
{
assert(seg < _p1Points.size());
assert(seg < _p2Points.size());
double omt = 1.0 - t;
Vec2 p0 = GetPoints()[seg];
Vec2 p1 = _p1Points[seg];
Vec2 p2 = _p2Points[seg];
Vec2 p3 = GetPoints()[seg+1];
double xVal = omt*omt*omt*p0.x + 3*omt*omt*t*p1.x +3*omt*t*t*p2.x+t*t*t*p3.x;
double yVal = omt*omt*omt*p0.y + 3*omt*omt*t*p1.y +3*omt*t*t*p2.y+t*t*t*p3.y;
return Vec2(xVal,yVal);
}
/* Clear out all the data.
*/
virtual void ResetDerived()
{
_p1Points.clear();
_p2Points.clear();
}
virtual bool ComputeSpline()
{
const vector<Vec2>& p = GetPoints();
int N = (int)p.size()-1;
_p1Points.resize(N);
_p2Points.resize(N);
if(N == 0)
return false;
if(N == 1)
{ // Only 2 points...just create a straight line.
// Constraint: 3*P1 = 2*P0 + P3
_p1Points[0] = (2.0/3.0*p[0] + 1.0/3.0*p[1]);
// Constraint: P2 = 2*P1 - P0
_p2Points[0] = 2.0*_p1Points[0] - p[0];
return true;
}
/*rhs vector*/
vector<Vec2> a(N);
vector<Vec2> b(N);
vector<Vec2> c(N);
vector<Vec2> r(N);
/*left most segment*/
a[0].x = 0;
b[0].x = 2;
c[0].x = 1;
r[0].x = p[0].x+2*p[1].x;
a[0].y = 0;
b[0].y = 2;
c[0].y = 1;
r[0].y = p[0].y+2*p[1].y;
/*internal segments*/
for (int i = 1; i < N - 1; i++)
{
a[i].x=1;
b[i].x=4;
c[i].x=1;
r[i].x = 4 * p[i].x + 2 * p[i+1].x;
a[i].y=1;
b[i].y=4;
c[i].y=1;
r[i].y = 4 * p[i].y + 2 * p[i+1].y;
}
/*right segment*/
a[N-1].x = 2;
b[N-1].x = 7;
c[N-1].x = 0;
r[N-1].x = 8*p[N-1].x+p[N].x;
a[N-1].y = 2;
b[N-1].y = 7;
c[N-1].y = 0;
r[N-1].y = 8*p[N-1].y+p[N].y;
/*solves Ax=b with the Thomas algorithm (from Wikipedia)*/
for (int i = 1; i < N; i++)
{
double m;
m = a[i].x/b[i-1].x;
b[i].x = b[i].x - m * c[i - 1].x;
r[i].x = r[i].x - m * r[i-1].x;
m = a[i].y/b[i-1].y;
b[i].y = b[i].y - m * c[i - 1].y;
r[i].y = r[i].y - m * r[i-1].y;
}
_p1Points[N-1].x = r[N-1].x/b[N-1].x;
_p1Points[N-1].y = r[N-1].y/b[N-1].y;
for (int i = N - 2; i >= 0; --i)
{
_p1Points[i].x = (r[i].x - c[i].x * _p1Points[i+1].x)/b[i].x;
_p1Points[i].y = (r[i].y - c[i].y * _p1Points[i+1].y)/b[i].y;
}
/*we have p1, now compute p2*/
for (int i=0;i<N-1;i++)
{
_p2Points[i].x=2*p[i+1].x-_p1Points[i+1].x;
_p2Points[i].y=2*p[i+1].y-_p1Points[i+1].y;
}
_p2Points[N-1].x = 0.5 * (p[N].x+_p1Points[N-1].x);
_p2Points[N-1].y = 0.5 * (p[N].y+_p1Points[N-1].y);
return true;
}
virtual void DumpDerived()
{
cout << " Control Points " << endl;
for(int idx = 0; idx < _p1Points.size(); idx++)
{
cout << "[" << idx << "] ";
cout << "P1: " << _p1Points[idx];
cout << " ";
cout << "P2: " << _p2Points[idx];
cout << endl;
}
}
};
#endif /* defined(__SplineCommon__) */
Algunas notas
- El spline clásico se bloqueará si le dan un conjunto vertical de puntos. Es por eso que creé el Bezier ... Tengo muchas líneas/caminos verticales a seguir. Podría ser modificado para simplemente dar una línea recta.
- La clase base tiene una opción para eliminar puntos colineales a medida que agrega ellos. Utiliza una comparación de pendiente simple de dos líneas para averiguar si están en la misma línea. No tiene que hacer esto, pero para las rutas largas que son rectas, reduce los ciclos. Cuando hace una gran cantidad de ruta en un gráfico de espaciado regular, tiende a obtener un lote de segmentos continuos.
Aquí es un ejemplo del uso del Bezier spline:
/* Smooth the points on the path so that turns look
* more natural. We'll only smooth the first few
* points. Most of the time, the full path will not
* be executed anyway...why waste cycles.
*/
void SmoothPath(vector<Vec2>& path, int32 divisions)
{
const int SMOOTH_POINTS = 6;
BezierSpine spline;
if(path.size() < 2)
return;
// Cache off the first point. If the first point is removed,
// the we occasionally run into problems if the collision detection
// says the first node is occupied but the splined point is too
// close, so the FSM "spins" trying to find a sensor cell that is
// not occupied.
// Vec2 firstPoint = path.back();
// path.pop_back();
// Grab the points.
for(int idx = 0; idx < SMOOTH_POINTS && path.size() > 0; idx++)
{
spline.AddPoint(path.back());
path.pop_back();
}
// Smooth them.
spline.ComputeSpline();
// Push them back in.
for(int idx = spline.GetPoints().size()-2; idx >= 0; --idx)
{
for(int division = divisions-1; division >= 0; --division)
{
double t = division*1.0/divisions;
path.push_back(spline.Eval(idx, t));
}
}
// Push back in the original first point.
// path.push_back(firstPoint);
}
Notas
- Si bien todo el camino podría ser suavizada, en esta aplicación, ya que el camino estaba cambiando De vez en cuando, era mejor simplemente alisar los primeros puntos y luego conectarlo.
- Los puntos se cargan en orden "inverso" en el vector de ruta. Este puede o no guardar ciclos (he dormido desde entonces).
Este código es parte de una base de código mucho más grande, but you can download it all on github y see a blog entry about it here.
Gracias por este código, he estado buscando uno que opere en cualquier conjunto de puntos 2D (no solo funciones con valores reales) durante muchas horas. acaba de hacer mi día. – Arshia001
- 1. ¿Qué (buenas) bibliotecas de servidor Java RADIUS existen?
- 2. ¿Qué buenas bibliotecas hay para resolver un sistema de ecuaciones no lineales en C++?
- 3. ¿Existen buenas bibliotecas de funciones para usar con colecciones en Java, como
- 4. ¿Buenas bibliotecas de sonido?
- 5. ¿Existen buenas herramientas para integrar GWT con ASP.Net?
- 6. ¿Existen buenas bibliotecas de terceros creadas en la parte superior de openCL todavía?
- 7. Buenas bibliotecas de red Python para construir un servidor TCP?
- 8. ¿Qué bibliotecas JavaScript multiplataforma existen?
- 9. ¿Existen buenas extensiones/plugins de e-commerce para Umbraco?
- 10. ¿Existen buenas asociaciones profesionales para gerentes/programadores de TI?
- 11. ¿Existen buenas suites de pruebas automáticas para Perl?
- 12. ¿Qué bibliotecas existen para crear UI "moving/living" en delphi?
- 13. Buenas bibliotecas de Python para la sincronización de iPod
- 14. ¿Qué son buenas bibliotecas JS para desarrolladores de juegos? (HTML5)
- 15. ¿Qué son buenas bibliotecas de scala para generar fuentes RSS?
- 16. ¿Existen colecciones documentadas y organizadas de bibliotecas para Common Lisp?
- 17. ¿Qué bibliotecas SOAP existen para Python 3.x?
- 18. Creando splines continuos/Haciendo una transición suave entre splines
- 19. ¿Qué objeto JavaScript existen para mapear bibliotecas de objetos?
- 20. Encontrar ciclos hamiltonianos en gráficos planos cúbicos
- 21. ¿Existen guías buenas, modernas y en línea para la optimización manual del código de ensamblaje?
- 22. Bibliotecas de optimización no lineal secuenciales en C++ WITH constraints
- 23. programa para resolver ecuaciones en C#
- 24. Buenas aplicaciones asp.net (C#)?
- 25. splines Catmull-Rom - ¿cómo funcionan?
- 26. Bibliotecas C++ para manipular imágenes
- 27. bibliotecas C# para la internacionalización?
- 28. Splines Mínimos de Python Natural
- 29. Ejemplos de interfaces C para bibliotecas C++?
- 30. C o C++: ¿Bibliotecas para factorizar enteros?
Esta biblioteca tiene un tiempo de O (n) y la aplicación de memoria para splines cúbicos penalizadas con suavizado automático utilizando la validación cruzada o grados de libertad efectivos similares a smooth.splines de R(). Ver skel__Cspplines.h y skel__TestCspplines.h: https://bitbucket.org/aperezrathke/skel – aprstar