sbasis-geometric.cpp revision 2160101b7d0f877a9d1c4f33c2233fab2ba1ea43
/** 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
\relates Piecewise
*/
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
\relates Piecewise, D2
*/
}
/** 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
\relates D2, SBasis
*/
}
/** 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
\relates Piecewise, D2
*/
}
/** 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)
\relates Piecewise, D2
*/
//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)
\relates Piecewise
*/
}
return result;
}
/** returns a function giving the arclength at each point in M.
\param M the Element.
\param tol the maximum error allowed.
\relates Piecewise
*/
return length;
}
/** returns a function giving the arclength at each point in M.
\param M the Element.
\param tol the maximum error allowed.
\relates Piecewise, D2
*/
}
#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.
\relates Piecewise, D2
\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.
\relates Piecewise
\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
\relates Piecewise, D2
*/
unsigned order,
double tol){
u.push_cut(0);
for (unsigned i=0; i < s.size();i++){
continue;
}
}
}
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
\relates Piecewise
*/
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
\relates D2
*/
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.
\relates D2
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.
\relates Piecewise
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.
\relates Piecewise
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.
*
* \relates D2
*/
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
//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 Geom {
/**
* \brief returns all the parameter values of A whose tangent passes through P.
* \relates D2
*/
}
/**
* \brief returns all the parameter values of A whose normal passes through P.
* \relates D2
*/
}
/**
* \brief returns all the parameter values of A whose normal is parallel to vector V.
* \relates D2
*/
}
/**
* \brief returns all the parameter values of A whose tangent is parallel to vector V.
* \relates D2
*/
}
}
//}; // 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:fileencoding=utf-8:textwidth=99 :