GCS: fundamental changes to simplify derivatives

GCS::Vector2D was morphed into a DeriVector2, a derivative-aware vector.
A bunch of vector math methods were added that implicitly calculate
respective derivatives. Now, there is no need to calculate the partials
- most is done implicitly.
This commit is contained in:
DeepSOIC 2014-12-11 23:26:33 +03:00 committed by wmayer
parent 3838bddcdb
commit 397e37ad1c
4 changed files with 182 additions and 175 deletions

View File

@ -428,8 +428,8 @@ int Sketch::addArcOfEllipse(const Part::GeomArcOfEllipse &ellipseSegment, bool f
a.start = p1;
a.end = p2;
a.center = p3;
a.focus1X = f1X;
a.focus1Y = f1Y;
a.focus1.x = f1X;
a.focus1.y = f1Y;
a.radmin = rmin;
a.startAngle = a1;
a.endAngle = a2;
@ -533,8 +533,8 @@ int Sketch::addEllipse(const Part::GeomEllipse &elip, bool fixed)
// set the ellipse for later constraints
GCS::Ellipse e;
e.focus1X = f1X;
e.focus1Y = f1Y;
e.focus1.x = f1X;
e.focus1.y = f1Y;
e.center = c;
e.radmin = rmin;
@ -1867,7 +1867,7 @@ bool Sketch::updateGeometry()
GeomArcOfEllipse *aoe = dynamic_cast<GeomArcOfEllipse*>(it->geo);
Base::Vector3d center = Vector3d(*Points[it->midPointId].x, *Points[it->midPointId].y, 0.0);
Base::Vector3d f1 = Vector3d(*myArc.focus1X, *myArc.focus1Y, 0.0);
Base::Vector3d f1 = Vector3d(*myArc.focus1.x, *myArc.focus1.y, 0.0);
double radmin = *myArc.radmin;
Base::Vector3d fd=f1-center;
@ -1897,7 +1897,7 @@ bool Sketch::updateGeometry()
GeomEllipse *ellipse = dynamic_cast<GeomEllipse*>(it->geo);
Base::Vector3d center = Vector3d(*Points[it->midPointId].x, *Points[it->midPointId].y, 0.0);
Base::Vector3d f1 = Vector3d(*Ellipses[it->index].focus1X, *Ellipses[it->index].focus1Y, 0.0);
Base::Vector3d f1 = Vector3d(*Ellipses[it->index].focus1.x, *Ellipses[it->index].focus1.y, 0.0);
double radmin = *Ellipses[it->index].radmin;
Base::Vector3d fd=f1-center;

View File

@ -769,8 +769,8 @@ int System::addConstraintInternalAlignmentEllipseMajorDiameter(Ellipse &e, Point
double Y_2=*p2.y;
double X_c=*e.center.x;
double Y_c=*e.center.y;
double X_F1=*e.focus1X;
double Y_F1=*e.focus1Y;
double X_F1=*e.focus1.x;
double Y_F1=*e.focus1.y;
double b=*e.radmin;
// P1=vector([X_1,Y_1])
@ -814,8 +814,8 @@ int System::addConstraintInternalAlignmentEllipseMinorDiameter(Ellipse &e, Point
double Y_2=*p2.y;
double X_c=*e.center.x;
double Y_c=*e.center.y;
double X_F1=*e.focus1X;
double Y_F1=*e.focus1Y;
double X_F1=*e.focus1.x;
double Y_F1=*e.focus1.y;
double b=*e.radmin;
// Same idea as for major above, but for minor
@ -841,8 +841,8 @@ int System::addConstraintInternalAlignmentEllipseMinorDiameter(Ellipse &e, Point
int System::addConstraintInternalAlignmentEllipseFocus1(Ellipse &e, Point &p1, int tagId)
{
addConstraintEqual(e.focus1X, p1.x, tagId);
return addConstraintEqual(e.focus1Y, p1.y, tagId);
addConstraintEqual(e.focus1.x, p1.x, tagId);
return addConstraintEqual(e.focus1.y, p1.y, tagId);
}
int System::addConstraintInternalAlignmentEllipseFocus2(Ellipse &e, Point &p1, int tagId)
@ -866,8 +866,8 @@ int System::addConstraintInternalAlignmentEllipseMajorDiameter(ArcOfEllipse &a,
double Y_2=*p2.y;
double X_c=*a.center.x;
double Y_c=*a.center.y;
double X_F1=*a.focus1X;
double Y_F1=*a.focus1Y;
double X_F1=*a.focus1.x;
double Y_F1=*a.focus1.y;
double b=*a.radmin;
// P1=vector([X_1,Y_1])
@ -911,8 +911,8 @@ int System::addConstraintInternalAlignmentEllipseMinorDiameter(ArcOfEllipse &a,
double Y_2=*p2.y;
double X_c=*a.center.x;
double Y_c=*a.center.y;
double X_F1=*a.focus1X;
double Y_F1=*a.focus1Y;
double X_F1=*a.focus1.x;
double Y_F1=*a.focus1.y;
double b=*a.radmin;
// Same idea as for major above, but for minor
@ -938,8 +938,8 @@ int System::addConstraintInternalAlignmentEllipseMinorDiameter(ArcOfEllipse &a,
int System::addConstraintInternalAlignmentEllipseFocus1(ArcOfEllipse &a, Point &p1, int tagId)
{
addConstraintEqual(a.focus1X, p1.x, tagId);
return addConstraintEqual(a.focus1Y, p1.y, tagId);
addConstraintEqual(a.focus1.x, p1.x, tagId);
return addConstraintEqual(a.focus1.y, p1.y, tagId);
}
int System::addConstraintInternalAlignmentEllipseFocus2(ArcOfEllipse &a, Point &p1, int tagId)
@ -956,14 +956,14 @@ double System::calculateAngleViaPoint(Curve &crv1, Curve &crv2, Point &p)
{return calculateAngleViaPoint(crv1, crv2, p, p);}
double System::calculateAngleViaPoint(Curve &crv1, Curve &crv2, Point &p1, Point &p2)
{
GCS::Vector2D n1 = crv1.CalculateNormal(p1);
GCS::Vector2D n2 = crv2.CalculateNormal(p2);
GCS::DeriVector2 n1 = crv1.CalculateNormal(p1);
GCS::DeriVector2 n2 = crv2.CalculateNormal(p2);
return atan2(-n2.x*n1.y+n2.y*n1.x, n2.x*n1.x + n2.y*n1.y);
}
void System::calculateNormalAtPoint(Curve &crv, Point &p, double &rtnX, double &rtnY)
{
GCS::Vector2D n1 = crv.CalculateNormal(p);
GCS::DeriVector2 n1 = crv.CalculateNormal(p);
rtnX = n1.x;
rtnY = n1.y;
}

View File

@ -29,36 +29,69 @@
#include "Geo.h"
namespace GCS{
Vector2D Line::CalculateNormal(Point &p, double* derivparam)
DeriVector2::DeriVector2(const Point &p, double *derivparam)
{
Vector2D p1v(*p1.x, *p1.y);
Vector2D p2v(*p2.x, *p2.y);
x=*p.x; y=*p.y;
dx=0.0; dy=0.0;
if (derivparam == p.x)
dx = 1.0;
if (derivparam == p.y)
dy = 1.0;
}
Vector2D ret(0.0, 0.0);
if(derivparam){
if(derivparam==this->p1.x){
ret.y += -1.0;
//ret.x += 0;
};
if(derivparam==this->p1.y){
//ret.y += 0;
ret.x += 1.0;
};
if(derivparam==this->p2.x){
ret.y += 1.0;
//ret.x += 0;
};
if(derivparam==this->p2.y){
//ret.y += 0;
ret.x += -1.0;
};
double DeriVector2::length(double &dlength)
{
double l = length();
if(l==0){
dlength = 1.0;
return l;
} else {
ret.x = -(p2v.y - p1v.y);
ret.y = (p2v.x - p1v.x);
};
dlength = (x*dx + y*dy)/l;
return l;
}
}
return ret;
DeriVector2 DeriVector2::getNormalized()
{
double l=length();
if(l==0.0) {
return DeriVector2(0, 0, dx/0.0, dy/0.0);
} else {
DeriVector2 rtn;
rtn.x = x/l;
rtn.y = y/l;
//first, simply scale the derivative accordingly.
rtn.dx = dx/l;
rtn.dy = dy/l;
//next, remove the collinear part of dx,dy (make a projection onto a normal)
double dsc = rtn.dx*rtn.x + rtn.dy*rtn.y;//scalar product d*v
rtn.dx -= dsc*rtn.x;//subtract the projection
rtn.dy -= dsc*rtn.y;
return rtn;
}
}
double DeriVector2::scalarProd(const DeriVector2 &v2, double *dprd)
{
if (dprd) {
*dprd = dx*v2.x + x*v2.dx + dy*v2.y + y*v2.dy;
};
return x*v2.x + y*v2.y;
}
DeriVector2 DeriVector2::divD(double val, double dval){
return DeriVector2(x/val,y/val,
dx/val - x*dval/(val*val),
dy/val - y*dval/(val*val)
);
}
DeriVector2 Line::CalculateNormal(Point &p, double* derivparam)
{
DeriVector2 p1v(p1, derivparam);
DeriVector2 p2v(p2, derivparam);
return p2v.subtr(p1v).rotate90ccw();
}
int Line::PushOwnParams(VEC_pD &pvec)
@ -86,35 +119,12 @@ Line* Line::Copy()
//---------------circle
Vector2D Circle::CalculateNormal(Point &p, double* derivparam)
DeriVector2 Circle::CalculateNormal(Point &p, double* derivparam)
{
Vector2D cv (*center.x, *center.y);
Vector2D pv (*p.x, *p.y);
DeriVector2 cv (center, derivparam);
DeriVector2 pv (p, derivparam);
Vector2D ret(0.0, 0.0);
if(derivparam){
if (derivparam == center.x) {
ret.x += 1;
ret.y += 0;
};
if (derivparam == center.y) {
ret.x += 0;
ret.y += 1;
};
if (derivparam == p.x) {
ret.x += -1;
ret.y += 0;
};
if (derivparam == p.y) {
ret.x += 0;
ret.y += -1;
};
} else {
ret.x = cv.x - pv.x;
ret.y = cv.y - pv.y;
};
return ret;
return cv.subtr(pv);
}
int Circle::PushOwnParams(VEC_pD &pvec)
@ -168,76 +178,43 @@ Arc* Arc::Copy()
//--------------ellipse
Vector2D Ellipse::CalculateNormal(Point &p, double* derivparam)
DeriVector2 Ellipse::CalculateNormal(Point &p, double* derivparam)
{
Vector2D cv (*center.x, *center.y);
Vector2D f1v (*focus1X, *focus1Y);
Vector2D pv (*p.x, *p.y);
//fill some vectors in
DeriVector2 cv (center, derivparam);
DeriVector2 f1v (focus1, derivparam);
DeriVector2 pv (p, derivparam);
Vector2D ret(0.0, 0.0);
//calculation.
//focus2:
DeriVector2 f2v = cv.linCombi(2.0, f1v, -1.0); // 2*cv - f1v
Vector2D f2v ( 2*cv.x - f1v.x, 2*cv.y - f1v.y ); //position of focus2
if(derivparam){
//pf1, pf2 = vectors from p to focus1,focus2
DeriVector2 pf1 = f1v.subtr(pv);
DeriVector2 pf2 = f2v.subtr(pv);
//return sum of normalized pf2, pf2
DeriVector2 ret = pf1.getNormalized().sum(pf2.getNormalized());
Vector2D dc;
Vector2D df1;
Vector2D dp;
if (derivparam == center.x) dc.x = 1.0;
if (derivparam == center.y) dc.y = 1.0;
if (derivparam == focus1X) df1.x = 1.0;
if (derivparam == focus1Y) df1.y = 1.0;
if (derivparam == p.x) dp.x = 1.0;
if (derivparam == p.y) dp.y = 1.0;
//todo: exit if all are zero
Vector2D pf1 = Vector2D(pv.x - f1v.x, pv.y - f1v.y);//same as during error function calculation. I reuse the values during derivative calculation
Vector2D pf2 = Vector2D(pv.x - f2v.x, pv.y - f2v.y);
Vector2D pf1e = pf1.getNormalized();
Vector2D pf2e = pf2.getNormalized();
Vector2D df2 (2*dc.x - df1.x, 2*dc.y - df1.y );
Vector2D dpf1 (dp.x - df1.x, dp.y - df1.y);//derivative before normalization
Vector2D dpf1e (dpf1.x/pf1.length(), dpf1.y/pf1.length());//first portion of normalization derivative (normalized' = unnormalized'/len + unnormalized*(1/len)')
dpf1e.x += -pf1.x/pow(pf1.length(),2)*(dpf1.x*pf1e.x + dpf1.y*pf1e.y);//second part of normalization dreivative
dpf1e.y += -pf1.y/pow(pf1.length(),2)*(dpf1.x*pf1e.x + dpf1.y*pf1e.y);
Vector2D dpf2 (dp.x - df2.x, dp.y - df2.y);//same stuff for pf2
Vector2D dpf2e (dpf2.x/pf2.length(), dpf2.y/pf2.length());//first portion of normalization derivative (normalized' = unnormalized'/len + unnormalized*(1/len)')
dpf2e.x += -pf2.x/pow(pf2.length(),2)*(dpf2.x*pf2e.x + dpf2.y*pf2e.y);//second part of normalization dreivative
dpf2e.y += -pf2.y/pow(pf2.length(),2)*(dpf2.x*pf2e.x + dpf2.y*pf2e.y);
ret.x = -(dpf1e.x + dpf2e.x);
ret.y = -(dpf1e.y + dpf2e.y);//DeepSOIC: derivative calculated manually... error-prone =) Tested, fixed, looks good.
//numeric derivatives for testing
#if 0 //make sure to enable DEBUG_DERIVS when enabling
double const eps = 0.00001;
double oldparam = *derivparam;
Vector2D v0 = this->CalculateNormal(p);
*derivparam += eps;
Vector2D vr = this->CalculateNormal(p);
*derivparam = oldparam - eps;
Vector2D vl = this->CalculateNormal(p);
*derivparam = oldparam;
//If not nasty, real derivative should be between left one and right one
Vector2D numretl ((v0.x-vl.x)/eps, (v0.y-vl.y)/eps);
Vector2D numretr ((vr.x-v0.x)/eps, (vr.y-v0.y)/eps);
assert(ret.x <= std::max(numretl.x,numretr.x) );
assert(ret.x >= std::min(numretl.x,numretr.x) );
assert(ret.y <= std::max(numretl.y,numretr.y) );
assert(ret.y >= std::min(numretl.y,numretr.y) );
#endif
} else {
Vector2D pf1 = Vector2D(pv.x - f1v.x, pv.y - f1v.y);
Vector2D pf2 = Vector2D(pv.x - f2v.x, pv.y - f2v.y);
Vector2D pf1e = pf1.getNormalized();
Vector2D pf2e = pf2.getNormalized();
ret.x = -(pf1e.x + pf2e.x);
ret.y = -(pf1e.y + pf2e.y);
};
//numeric derivatives for testing
#if 0 //make sure to enable DEBUG_DERIVS when enabling
if(derivparam) {
double const eps = 0.00001;
double oldparam = *derivparam;
DeriVector2 v0 = this->CalculateNormal(p);
*derivparam += eps;
DeriVector2 vr = this->CalculateNormal(p);
*derivparam = oldparam - eps;
DeriVector2 vl = this->CalculateNormal(p);
*derivparam = oldparam;
//If not nasty, real derivative should be between left one and right one
DeriVector2 numretl ((v0.x-vl.x)/eps, (v0.y-vl.y)/eps);
DeriVector2 numretr ((vr.x-v0.x)/eps, (vr.y-v0.y)/eps);
assert(ret.dx <= std::max(numretl.x,numretr.x) );
assert(ret.dx >= std::min(numretl.x,numretr.x) );
assert(ret.dy <= std::max(numretl.y,numretr.y) );
assert(ret.dy >= std::min(numretl.y,numretr.y) );
}
#endif
return ret;
}
@ -247,8 +224,8 @@ int Ellipse::PushOwnParams(VEC_pD &pvec)
int cnt=0;
pvec.push_back(center.x); cnt++;
pvec.push_back(center.y); cnt++;
pvec.push_back(focus1X); cnt++;
pvec.push_back(focus1Y); cnt++;
pvec.push_back(focus1.x); cnt++;
pvec.push_back(focus1.y); cnt++;
pvec.push_back(radmin); cnt++;
return cnt;
}
@ -256,8 +233,8 @@ void Ellipse::ReconstructOnNewPvec(VEC_pD &pvec, int &cnt)
{
center.x=pvec[cnt]; cnt++;
center.y=pvec[cnt]; cnt++;
focus1X=pvec[cnt]; cnt++;
focus1Y=pvec[cnt]; cnt++;
focus1.x=pvec[cnt]; cnt++;
focus1.y=pvec[cnt]; cnt++;
radmin=pvec[cnt]; cnt++;
}
Ellipse* Ellipse::Copy()

View File

@ -28,27 +28,9 @@
namespace GCS
{
class Vector2D /* DeepSOIC: I tried to reuse Base::Vector2D by #include <Base/Tools2D.h>,
* but I failed to resolve bullshit compilation errors that arose in the process, so...
* Anyway, the benefit is that solver has less dependencies on FreeCAD and can be
* stripped off easier.
* I could have used Eigen's Vector2f, but I found it overblown and too complex to use.
*/
{
public:
Vector2D(){x=0; y=0;}
Vector2D(double x, double y) {this->x = x; this->y = y;}
double x;
double y;
double length() {return sqrt(x*x + y*y);}
const double NAN = std::numeric_limits<double>::quiet_NaN();
const double INF = std::numeric_limits<double>::infinity();
//unlike other vectors in FreeCAD, this normalization creates a new vector instead of modifying existing one.
Vector2D getNormalized(){double l=length(); if(l==0.0) l=1.0; return Vector2D(x/l,y/l);} //returns zero vector if the original is zero.
};
///////////////////////////////////////
// Geometries
///////////////////////////////////////
class Point
{
public:
@ -57,6 +39,56 @@ namespace GCS
double *y;
};
///Class DeriVector2 holds a vector value and its derivative on the
///parameter that the derivatives are being calculated for now. x,y is the
///actual vector (v). dx,dy is a derivative of the vector by a parameter
///(dv/dp). The easiest way to fill the vector in is by passing a point and
///a derivative parameter pointer to its constructor. x,y are read from the
///pointers in Point, and dx,dy are set to either 0 or 1 depending on what
///pointers of Point match the supplied pointer. The derivatives can be set
///manually as well. The class also provides a bunch of methods to do math
///on it (and derivatives are calculated implicitly).
///
class DeriVector2
{
public:
DeriVector2(){x=0; y=0; dx=0; dy=0;}
DeriVector2(double x, double y) {this->x = x; this->y = y; this->dx = 0; this->dy = 0;}
DeriVector2(double x, double y, double dx, double dy) {this->x = x; this->y = y; this->dx = dx; this->dy = dy;}
DeriVector2(const Point &p, double* derivparam);
double x, dx;
double y, dy;
double length() {return sqrt(x*x + y*y);}
double length(double &dlength); //returns length and writes length deriv into the dlength argument.
//unlike other vectors in FreeCAD, this normalization creates a new vector instead of modifying existing one.
DeriVector2 getNormalized(); //returns zero vector if the original is zero.
double scalarProd(const DeriVector2 &v2, double* dprd=0);//calculates scalar product of two vectors and returns the result. The derivative of the result is written into argument dprd.
DeriVector2 sum(const DeriVector2 &v2){//adds two vectors and returns result
return DeriVector2(x + v2.x, y + v2.y,
dx + v2.dx, dy + v2.dy);}
DeriVector2 subtr(const DeriVector2 &v2){//subtracts two vectors and returns result
return DeriVector2(x - v2.x, y - v2.y,
dx - v2.dx, dy - v2.dy);}
DeriVector2 mult(double val){
return DeriVector2(x*val, y*val, dx*val, dy*val);}//multiplies the vector by a number. Derivatives are scaled.
DeriVector2 multD(double val, double dval){//multiply vector by a variable with a derivative.
return DeriVector2(x*val, y*val, dx*val+x*dval, dy*val+y*dval);}
DeriVector2 divD(double val, double dval);//divide vector by a variable with a derivative
DeriVector2 rotate90ccw(){return DeriVector2(-y,x,-dy,dx);}
DeriVector2 rotate90cw(){return DeriVector2(y,-x,dy,-dx);}
DeriVector2 linCombi(double m1, const DeriVector2 &v2, double m2){//linear combination of two vectors
return DeriVector2(x*m1 + v2.x*m2, y*m1 + v2.y*m2,
dx*m1 + v2.dx*m2, dy*m1 + v2.dy*m2);}
};
///////////////////////////////////////
// Geometries
///////////////////////////////////////
class Curve //a base class for all curve-based objects (line, circle/arc, ellipse/arc)
{
public:
@ -66,10 +98,9 @@ namespace GCS
// assumed to be walked counterclockwise, so the vector should point
// into the shape.
//derivparam is a pointer to a curve parameter (or point coordinate) to
// compute the derivative for. if derivparam is nullptr, the actual
// normal vector is returned, otherwise a derivative of normal vector by
// *derivparam is returned.
virtual Vector2D CalculateNormal(Point &p, double* derivparam = 0) = 0;
// compute the derivative for. The derivative is returned through dx,dy
// fields of DeriVector2.
virtual DeriVector2 CalculateNormal(Point &p, double* derivparam = 0) = 0;
//adds curve's parameters to pvec (used by constraints)
virtual int PushOwnParams(VEC_pD &pvec) = 0;
@ -85,7 +116,7 @@ namespace GCS
Line(){}
Point p1;
Point p2;
Vector2D CalculateNormal(Point &p, double* derivparam = 0);
DeriVector2 CalculateNormal(Point &p, double* derivparam = 0);
virtual int PushOwnParams(VEC_pD &pvec);
virtual void ReconstructOnNewPvec (VEC_pD &pvec, int &cnt);
virtual Line* Copy();
@ -97,7 +128,7 @@ namespace GCS
Circle(){rad = 0;}
Point center;
double *rad;
Vector2D CalculateNormal(Point &p, double* derivparam = 0);
DeriVector2 CalculateNormal(Point &p, double* derivparam = 0);
virtual int PushOwnParams(VEC_pD &pvec);
virtual void ReconstructOnNewPvec (VEC_pD &pvec, int &cnt);
virtual Circle* Copy();
@ -123,10 +154,9 @@ namespace GCS
public:
Ellipse(){ radmin = 0;}
Point center;
double *focus1X;
double *focus1Y;
Point focus1;
double *radmin;
Vector2D CalculateNormal(Point &p, double* derivparam = 0);
DeriVector2 CalculateNormal(Point &p, double* derivparam = 0);
virtual int PushOwnParams(VEC_pD &pvec);
virtual void ReconstructOnNewPvec (VEC_pD &pvec, int &cnt);
virtual Ellipse* Copy();
@ -142,8 +172,8 @@ namespace GCS
Point start;
Point end;
//Point center; //inherited
//double *focus1X; //inherited
//double *focus1Y; //inherited
//double *focus1.x; //inherited
//double *focus1.y; //inherited
virtual int PushOwnParams(VEC_pD &pvec);
virtual void ReconstructOnNewPvec (VEC_pD &pvec, int &cnt);
virtual ArcOfEllipse* Copy();