Sketcher: Ellipse support

Here's the place for discussion related to coding in FreeCAD, C++ or Python. Design, interfaces and structures.
User avatar
DeepSOIC
Posts: 7249
Joined: Fri Aug 29, 2014 12:45 am
Location: Saint-Petersburg, Russia

Re: Sketcher: Ellipse support

Postby DeepSOIC » Sun Sep 21, 2014 5:48 pm

numerical derivative is very simple:
dF(a,b,c,d)/dc=(approx) ( F(a,b,c+eps,d) - F(a,b,c-eps,d) ) / (2*eps)
The key to getting it right is the choice of a small value eps. It should be low enough so that F is reasonably linear in the c+eps..c-eps, yet large enough so that the difference F-F is not too limited by float precision/F internal computation precision.
ulrich1a
Posts: 1920
Joined: Sun Jul 07, 2013 12:08 pm

Re: Sketcher: Ellipse support

Postby ulrich1a » Sun Sep 21, 2014 7:11 pm

abdullah wrote:Ok. I do not fully follow the advantages of using minor instead of major, but for me it is ok. Could be agree to one or the other, so that we do not have to recalculate the partials next week because of changing from one to the other? Ulrich, as you proposed first, what do you think?
I do also not understand the advantage of using the minor axis. Using the major axis has the following advantages:
The error function for the point on ellipse uses the major axis. So the existing function can be used and no extra calculation step is needed.
The major axis is the one parameter in the pins and string construction method.
The creation of the ellipse is more logical to me with the major axis first and I understand is also expected as the normal case.

Ulrich
wmayer
Site Admin
Posts: 15130
Joined: Thu Feb 19, 2009 10:32 am

Re: Sketcher: Ellipse support

Postby wmayer » Sun Sep 21, 2014 7:28 pm

What DeepSOIC said is that for the given focus points F1 and F2 you cannot choose an arbitrary value for the length "a" of the major axis while you can do this for the length "b" of the minor axis. For "b" it suffices to choose b > 0, however for "a" you must choose that "a > 0.5*dist(F1,F2)" because the focus points must be inside the ellipse.

The distance of a focus point to the center C is called "numerical eccentricity": e = dist(F1,C) = dist(F2,C)=0.5*dist(F1,F2) and the relation between e,a,b is;
a*a = e*e + b*b which makes it clear that a > e.
User avatar
DeepSOIC
Posts: 7249
Joined: Fri Aug 29, 2014 12:45 am
Location: Saint-Petersburg, Russia

Re: Sketcher: Ellipse support

Postby DeepSOIC » Sun Sep 21, 2014 7:32 pm

Hey! A generalization for conics will benefit. Clear advantage for major axis way! (a parabola still suffers.. maybe an extra look at different options is worth?)
EDIT: I was about to say that if major radius is smaller than 0.5*|F1-F2|, it is not an ellipse but a hyperbola, but wmayer was first.
ulrich1a
Posts: 1920
Joined: Sun Jul 07, 2013 12:08 pm

Re: Sketcher: Ellipse support

Postby ulrich1a » Sun Sep 21, 2014 8:59 pm

Here is a small patch to fix the constraint incompatibility to the master branch.

Ulrich
Attachments
0001-fix-Solver-constraint-incompatibility-to-master.patch.zip
(677 Bytes) Downloaded 18 times
User avatar
DeepSOIC
Posts: 7249
Joined: Fri Aug 29, 2014 12:45 am
Location: Saint-Petersburg, Russia

Re: Sketcher: Ellipse support

Postby DeepSOIC » Sun Sep 21, 2014 9:08 pm

Sorry, switching topics again.
abdullah wrote:If the constraints are fully independent (the system matrix formed by the gradients has independent rows for each, i.e. if the rank of the matrix produced by the several separate constraints is equal to the number of constraints.), then you can certainly apply two or more solver constraints to implement an UI constraint, and even more, you are encouraged to do so, as this simplifies the code. OIne simple example is the coincident constraint, implemented as two equality constraints:

int System::addConstraintP2PCoincident(Point &p1, Point &p2, int tagId)
{
addConstraintEqual(p1.x, p2.x, tagId);
return addConstraintEqual(p1.y, p2.y, tagId);
}
I have a more human-readable way to tell this: addConstraintXXX is supposed to remove exactly one degree of freedom. UI constraint "coincident" has to remove two degrees of freedom, so has two addConstraintXXX calls. Correct?

Conclustion: a tangent constraint must be a single new row into the matrix, since it removes only one degree of freedom.
Or no? More rows can be added if we introduce additional degrees of freedom.
User avatar
DevJohan
Posts: 41
Joined: Sun Jul 13, 2014 2:36 pm
Location: Stockholm, Sweden

Re: Sketcher: Ellipse support

Postby DevJohan » Sun Sep 21, 2014 9:23 pm

Numeric derivatives should be avoided whenever possible. And it is possible to get code that is much more efficient. A quartic equation has algebraic solutions you just have to choose the right one. I have actually done some computations these last couple of days to visualize the solutions. Assuming A is the major axis the following code produce the image of the theta-angle (named t by abdullah I think) of the closest point on the ellipse.

Code: Select all

A=3; B=2; base_npts = 100;
nptsA = A*base_npts;
nptsB = B*base_npts;

px = 2*(.5+(-nptsA:nptsA))/base_npts;
Px = ones(2*nptsB+1,1)*px;
py = 2*(.5+(-nptsB:nptsB))/base_npts;
Py = py'*ones(1,2*nptsA+1);

a=-(1i*B^2-1i*A^2)/4;
b=-(B*Py+1i*A*Px)/2;
c=0;
d=-(B*Py-1i*A*Px)/2;
e=(1i*B^2-1i*A^2)/4;

th = pi*(0:360)/180;

p     = -3*b.^2./(8*a.^2);
q     =  (b.^3+8*a.^2.*d)./(8*a.^3);
dta0  =  real(12*a.*e-3*b.*d);
dta1  =  real(27*b.^2.*e+27*a.*d.^2);
QP1   =  (dta1.^2-4*dta0.^3);
Qcube =  (dta1+sqrt(QP1))/2;
Q     =  nthroot(real(Qcube),3).*double(imag(Qcube)==0) + (Qcube).^(1/3).*double(~(imag(Qcube)==0));
S     =  1/2*sqrt(-2*p/3+(Q+dta0./Q)./(3*a));
R1    =  -4*S.^2-2*p+q./S;
R2    =  -4*S.^2-2*p-q./S;

X1  = -b./(4*a) - S + 1/2*sqrt(R1);
X2  = -b./(4*a) + S + 1/2*sqrt(R2);
X3  = -b./(4*a) - S - 1/2*sqrt(R1);
X4  = -b./(4*a) + S - 1/2*sqrt(R2);

compl = ...
    X1.*double(   (Px > 0 & QP1 > 0 ) | (Py >0 & QP1 < 0 )  ) + ...
    X4.*double(~( (Px > 0 & QP1 > 0 ) | (Py >0 & QP1 < 0 ) )) ;
T=-1i*log( compl );

figure(1);
image(px,py,64/(2*pi)*(pi+real(T)));
axis equal; hold on;
plot(A*cos(th),B*sin(th),'g','LineWidth',2);
plot([sqrt(A^2-B^2) -sqrt(A^2-B^2)],[0 0],'g*','LineWidth',2);
Which produce the following image.
elipseangle.png
elipseangle.png (53.35 KiB) Viewed 1066 times
This code doesn't take the phi rotation or translation into account but this is just an angle preserving coordinate transform. Getting the partials from this code is quite straightforward and doesn't require 600kb of code. Still it doesn't address the multiple possible solutions but these could also be handled by checking different conditions. For example there are four solutions if QP1 < 0 for a given coordinate.
User avatar
DeepSOIC
Posts: 7249
Joined: Fri Aug 29, 2014 12:45 am
Location: Saint-Petersburg, Russia

Re: Sketcher: Ellipse support

Postby DeepSOIC » Sun Sep 21, 2014 9:24 pm

DevJohan wrote:To avoid the solution jumping from one point on the ellipse and circle to another point on the ellipse and circle I think some kind of extra data is required. An intermediate line is one possibility where you constrain one end of the line to the ellipse and circle curves and add tangency constraints ellipse-line and circle-line.

This tangency point does not necessarily have to be visible, but it can be quite useful if it is.
This is brilliant! Here's how I see it:
A user selects two curves and a point and hits "tangent" button (if no point is supplied, one is created).
The point is constrained onto two objects (two point-on-object constraints, already done). And the final constraint uses the point in calculation of a new, much simpler error function, assuming that the point is the point of tangency. I can take care for making such function, for line, circle and possibly even second ellipse, I feel it's easy! It's just calculating tangent vectors to the curves at the point and F=length of vector product of those two vectors.
If the user deletes the point, it will automatically kill the constraint since it will be associated with the constraint.
What do you think?
User avatar
marktaff
Posts: 47
Joined: Sun Sep 21, 2014 3:25 pm
Location: Seattle, Washington, USA

Re: Sketcher: Ellipse support

Postby marktaff » Sun Sep 21, 2014 10:12 pm

I'm still working on setting up my dev environment, and haven't even peeked at the code yet, though I have read this thread. So take this with a grain of salt--I don't want to step on anyone's toes. :-)

In reading this thread, I've seen you guys beat yourselves up doing this conic section math in cartesian coordinates instead of in polar coordinates where it is much easier, perhaps in a perifocal frame. Is there a reason we aren't using polar forms, and just converting givens/results to and from cartesian forms as required?

Thanks for the clarification.
Available on chat.freenode.net#freecad while working on FreeCad
User avatar
DeepSOIC
Posts: 7249
Joined: Fri Aug 29, 2014 12:45 am
Location: Saint-Petersburg, Russia

Re: Sketcher: Ellipse support

Postby DeepSOIC » Sun Sep 21, 2014 10:35 pm

marktaff wrote:In reading this thread, I've seen you guys beat yourselves up doing this conic section math in cartesian coordinates instead of in polar coordinates where it is much easier, perhaps in a perifocal frame. Is there a reason we aren't using polar forms, and just converting givens/results to and from cartesian forms as required?
Maybe because I don't know how to describe a circle with a center offset from origin in polar coordinates :oops: ... I'm doing my best! I'll start another thread in my brain for "perifocal frame", I haven't thought about it.
EDIT: I also don't know how an ellipse is defined in polar coordinates! (even when it is centered or focused at the origin)
EDIT2:ellipse in perifocal coordinates (looked up in Wikipedia):
ellipse in focal polar.png
ellipse in focal polar.png (1 KiB) Viewed 1044 times
Quite simple! Looks promising.