sbasis-geometric.cpp revision 57d645121d3dc83796c5fa418a2047b329f74f29
/** Geometric operators on D2<SBasis> (1D->2D).
* Copyright 2007 JF Barraud
* Copyright 2007 N Hurst
*
* The functions defined in this header related to 2d geometric operations such as arc length,
* unit_vector, curvature, and centroid. Most are built on top of unit_vector, which takes an
* arbitrary D2 and returns a D2 with unit length with the same direction.
*
* arclength D2 -> sbasis (giving arclength function)
* does uniform_speed return natural parameterisation?
* integrate sb2d code from normal-bundle
* angle(md<2>) -> sbasis (gives angle from vector - discontinuous?)
* osculating circle center?
*
**/
//namespace Geom{
using namespace Geom;
using namespace std;
//Some utils first.
//TODO: remove this!!
/**
* Return a list of doubles that appear in both a and b to within error tol
* a, b, vector of double
* tol tolerance
*/
unsigned i=0,j=0;
i+=1;
j+=1;
}else if (a[i]<b[j]){
i+=1;
}else if (a[i]>b[j]){
j+=1;
}
}
return inter;
}
//------------------------------------------------------------------------------
if ( k>=(int)a.size()){
//make sure a is 0?
return SBasis();
}
if(k < 0) return shift(a,-k);
SBasis c;
return c;
}
if(k < 0) {
for (int i=2; i<=-k; i++){
c*=c;
}
c*=a;
return(c);
}else{
for (int i=2; i<=k; i++){
c*=c;
}
c*=a;
return(divide_by_sk(c,k));
}
}
if(k < 0) {
for (int i=2; i<=-k; i++){
c*=c;
}
c*=a;
return(c);
}else{
for (int i=2; i<=k; i++){
c*=c;
}
c*=a;
return(divide_by_sk(c,k));
}
}
//TODO: divide by all the s at once!!!
M[0] = divide_by_sk(M[0],1);
}
M[0] = divide_by_t0k(M[0],1);
}
M[0] = divide_by_t1k(M[0],1);
}
return M;
}
/*static D2<SBasis> RescaleForNonVanishing(D2<SBasis> const &MM, double ZERO=1.e-4){
std::vector<double> levels;
levels.push_back(-ZERO);
levels.push_back(ZERO);
//std::vector<std::vector<double> > mr = multi_roots(MM, levels);
}*/
//=================================================================
//TODO: what's this for?!?!
for (unsigned i=0; i<M.size(); i++){
}
}
}
/** Return a function which gives the angle of vect at each point.
\param vect a piecewise parameteric curve.
\param tol the maximum error allowed.
\param order the maximum degree to use for approximation
*/
for (unsigned i=0; i<v.size(); i++){
//TODO: I don't understand this - sign.
//TODO: deal with 2*pi jumps form one seg to the other...
//TODO: not exact at t=1 because of the integral.
//TODO: force continuity?
}
return result;
}
/** Return a function which gives the angle of vect at each point.
\param vect a piecewise parameteric curve.
\param tol the maximum error allowed.
\param order the maximum degree to use for approximation
*/
}
/** tan2 is the pseudo-inverse of atan2. It takes an angle and returns a unit_vector that points in the direction of angle.
\param angle a piecewise function of angle wrt t.
\param tol the maximum error allowed.
\param order the maximum degree to use for approximation
*/
}
/** tan2 is the pseudo-inverse of atan2. It takes an angle and returns a unit_vector that points in the direction of angle.
\param angle a piecewise function of angle wrt t.
\param tol the maximum error allowed.
\param order the maximum degree to use for approximation
*/
}
/** Return a Piecewise<D2<SBasis> > which points in the same direction as V_in, but has unit_length.
\param V_in the original path.
\param tol the maximum error allowed.
\param order the maximum degree to use for approximation
unitVector(x,y) is computed as (b,-a) where a and b are solutions of:
ax+by=0 (eqn1) and a^2+b^2=1 (eqn2)
*/
//TODO: Handle vanishing vectors...
// -This approach is numerically bad. Find a stable way to rescale V_in to have non vanishing ends.
// -This done, unitVector will have jumps at zeros: fill the gaps with arcs of circles.
SBasis x = V[0], y = V[1];
r_eqn1 = -(a*x+b*y);
for (unsigned k=1; k<=order; k++){
//the equations to solve at this point are:
// a0*x(0)+b0*y(0)=r0 & 2*a0*a(0)+2*b0*b(0)=rr0
//and
// a1*x(1)+b1*y(1)=r1 & 2*a1*a(1)+2*b1*b(1)=rr1
//TODO: use "incremental" rather than explicit formulas.
r_eqn1 = -(a*x+b*y);
}
//our candidate is:
unitV[0] = b;
unitV[1] = -a;
//is it good?
//if not: subdivide and concat results.
return(unitV0);
}else{
//if yes: return it as pw.
return result;
}
}
/** Return a Piecewise<D2<SBasis> > which points in the same direction as V_in, but has unit_length.
\param V_in the original path.
\param tol the maximum error allowed.
\param order the maximum degree to use for approximation
unitVector(x,y) is computed as (b,-a) where a and b are solutions of:
ax+by=0 (eqn1) and a^2+b^2=1 (eqn2)
*/
}
return result;
}
/** returns a function giving the arclength at each point in M.
\param M the Element.
\param tol the maximum error allowed.
*/
return length;
}
/** returns a function giving the arclength at each point in M.
\param M the Element.
\param tol the maximum error allowed.
*/
}
#if 0
double
double tol){
}
double
double tol){
}
#endif
/** returns a function giving the curvature at each point in M.
\param M the Element.
\param tol the maximum error allowed.
Todo:
* claimed incomplete. Check.
*/
return(k);
}
/** returns a function giving the curvature at each point in M.
\param M the Element.
\param tol the maximum error allowed.
Todo:
* claimed incomplete. Check.
*/
}
return result;
}
//=================================================================
/** Reparameterise M to have unit speed.
\param M the Element.
\param tol the maximum error allowed.
\param order the maximum degree to use for approximation
*/
unsigned order,
double tol){
u.push_cut(0);
for (unsigned i=0; i < s.size();i++){
}
}
return u;
}
/** Reparameterise M to have unit speed.
\param M the Element.
\param tol the maximum error allowed.
\param order the maximum degree to use for approximation
*/
unsigned order,
double tol){
for (unsigned i=0; i<M.size(); i++ ){
}
return(result);
}
#include <gsl/gsl_integration.h>
static double sb_length_integrating(double t, void* param) {
}
/** Calculates the length of a D2<SBasis> through gsl integration.
\param B the Element.
\param tol the maximum error allowed.
\param result variable to be incremented with the length of the path
\param abs_error variable to be incremented with the estimated error
*/
gsl_function F;
= gsl_integration_workspace_alloc (20);
double quad_result, err;
/* We could probably use the non adaptive code here if we removed any cusps first. */
result += quad_result;
}
/** Calculates the length of a D2<SBasis> through gsl integration.
\param s the Element.
\param tol the maximum error allowed.
If you only want the total length, this routine faster and more accurate than constructing an arcLengthSb.
*/
double
double tol){
double result = 0;
double abs_error = 0;
return result;
}
/** Calculates the length of a Piecewise<D2<SBasis> > through gsl integration.
\param s the Element.
\param tol the maximum error allowed.
If you only want the total length, this routine faster and more accurate than constructing an arcLengthSb.
*/
double
double tol){
double result = 0;
double abs_error = 0;
for (unsigned i=0; i < s.size();i++){
}
return result;
}
/**
* Centroid using sbasis integration.
\param p the Element.
\param centroid on return contains the centroid of the shape
\param area on return contains the signed area of the shape.
This approach uses green's theorem to compute the area and centroid using integrals. For curved shapes this is much faster than converting to polyline. Note that without an uncross operation the output is not the absolute area.
* Returned values:
0 for normal execution;
2 if area is zero, meaning centroid is meaningless.
*/
Point centroid_tmp(0,0);
double atmp = 0;
for(unsigned i = 0; i < p.size(); i++) {
}
// join ends
centroid_tmp *= 2;
if (atmp != 0) {
return 0;
}
return 2;
}
/**
* Find cubics with prescribed curvatures at both ends.
*
* this requires to solve a system of the form
*
* \f[
* \lambda_1 = a_0 \lambda_0^2 + c_0
* \lambda_0 = a_1 \lambda_1^2 + c_1
* \f]
*
* which is a deg 4 equation in lambda 0.
* Below are basic functions dedicated to solving this assuming a0 and a1 !=0.
*/
static OptInterval
int insist_on_speeds_signs){
if ( delta < 0 )
return OptInterval();//return empty interval
if (flip)
result *= -1;
if (insist_on_speeds_signs == 1){
return OptInterval();//return empty interval
}
return result;
}
static
int insist_on_speeds_signs){
if ( !domain )
}
return rts;
}
/**
* \brief returns the cubics fitting direction and curvature of a given
* input curve at two points.
*
* The input can be the
* value, speed, and acceleration
* or
* value, speed, and cross(acceleration,speed)
* of the original curve at the both ends.
* (the second is often technically usefull, as it avoids unnecessary division by |v|^2)
* Recall that K=1/R=cross(acceleration,speed)/|speed|^3.
*
* Moreover, a 7-th argument 'insist_on_speed_signs' can be supplied to select solutions:
* If insist_on_speed_signs == 1, only consider solutions where speeds at both ends are positively
* proportional to the given ones.
* If insist_on_speed_signs == 0, allow speeds to point in the opposite direction (both at the same time)
* If insist_on_speed_signs == -1, allow speeds to point in both direction independantly.
*/
double epsilon){
//speed of cubic bezier will be lambda0*dM0 and lambda1*dM1,
//with lambda0 and lambda1 s.t. curvature at both ends is the same
//as the curvature of the given curve.
return result;
}
return result;
}
}else{
//solve: lambda1 = a0 lambda0^2 + c0
// lambda0 = a1 lambda1^2 + c1
}else{
//find lamda0 by solving a deg 4 equation d0+d1*X+...+d4*X^4=0
double a[5];
a[1] = -1;
a[3] = 0;
//vector<double> solns=solve_poly(a,4);
//is this solution pointing in the + direction at both ends?
}
//is this solution pointing in the - direction at both ends?
}
//ok,this solution is pointing in the + and - directions.
else if (insist_on_speed_signs<0){
}
}
}
}
}
#if 0
#endif
}
return(result);
}
double epsilon){
}
double epsilon){
double length;
}
//}; // namespace
/*
Local Variables:
mode:c++
c-file-style:"stroustrup"
c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
indent-tabs-mode:nil
fill-column:99
End:
*/
// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :