TK105000-Add per unit lineprofile (form only) for circle

huangyhg

I start to working on per unit line profile when the considered feature is circle or ellipse.
[11/11/2008 11:50:17 PM] Dan Zwick says: Ok. There's a function that finds the closest point on an ellipse.
[11/11/2008 11:50:30 PM] yanhua Huang says: do you have a function to find meas point according to the given arc length.
[11/11/2008 11:51:35 PM] Dan Zwick says: Hmmm, could be. Arclength on an ellipse is a special function, an "elliptic function."
[11/11/2008 11:51:54 PM] Dan Zwick says: If I don't already have that I'm sure I could find one somewhere.
[11/11/2008 11:57:48 PM] Dan Zwick says: double _DSFITEXPORT CDSFitEllipse2D::closest_point_on_ellipse(double ra, double rb, double x0, double y0,
double alpha, double z0, double w0,
double *t0, double *xc, double *yc)
[11/12/2008 12:19:07 AM] yanhua Huang says: /*-------------------------------------------------------------------*/
/* Find the closest point to a given point on an ellipse. */
/*-------------------------------------------------------------------*/

double _DSFITEXPORT CDSFitEllipse2D::closest_point_on_ellipse(double ra, double rb, double x0, double y0,
double alpha, double z0, double w0,
double *t0, double *xc, double *yc)
the comment is right?
[11/12/2008 12:20:26 AM] yanhua Huang says: I want to calculate ellipse lenth between two meas points.
[11/12/2008 12:27:30 AM] Dan Zwick says: Well, first you need to find the closest points on the ellipse and their parameters. There is no closed form expression for arc length on an ellipse, other than in terms of certain special functions. It might suffice to use Gaussian quadrature to approximate the arc length, though.
[11/12/2008 12:31:24 AM] yanhua Huang says: approximate computation is OK
[11/12/2008 12:31:59 AM] Dan Zwick says: Let me study this for a bit longer.
[11/12/2008 12:33:44 AM] yanhua Huang says: By the way, recently I have a PR about the low efficiency when compute line profile deviation.
[11/12/2008 12:34:12 AM] Dan Zwick says: For straight lines?
[11/12/2008 12:34:39 AM] yanhua Huang says: when the considered feature is scan....
[11/12/2008 12:34:50 AM] Dan Zwick says: Does that mean it takes too long?
[11/12/2008 12:35:07 AM] yanhua Huang says: try to work on it tomorrow
[11/12/2008 12:35:38 AM] yanhua Huang says: PR256106 - Profile dimension calculation very slow
[11/12/2008 12:35:49 AM] yanhua Huang says: (chuckle)
[11/12/2008 12:35:52 AM] Dan Zwick says: How many points are there in that scan?
[11/12/2008 12:36:01 AM] Dan Zwick says: Some people expect miracles.
[11/12/2008 12:36:19 AM] yanhua Huang says: >900
[11/12/2008 12:36:45 AM] Dan Zwick says: Is this profile per unit length?
[11/12/2008 12:37:11 AM] yanhua Huang says: no, old problem...
[11/12/2008 12:37:48 AM] Dan Zwick says: 2D or 3D?
[11/12/2008 12:37:53 AM] yanhua Huang says: Have not investigated the issue, just verify the issue.
[11/12/2008 12:38:59 AM] yanhua Huang says: Let me trace the code firstly tomorrow. Just let you know....now
[11/12/2008 12:39:17 AM] Dan Zwick says: Ok
[11/12/2008 1:12:02 AM] Dan Zwick says: Take a look at ellipse_arc.cpp in MATHTOOL.
[11/12/2008 1:13:58 AM] Dan Zwick says: This is used in HIT_HELP.CPP.
[11/12/2008 1:16:58 AM] yanhua Huang says: which name is the function called?
[11/12/2008 1:17:50 AM] Dan Zwick says: I suppose you need either ellipse_arc or ellipse_arc_inverse.
[11/12/2008 1:19:48 AM] Dan Zwick says: It looks like I implemented all of that ellipse stuff, but I have completely forgotten it.
[11/12/2008 1:20:08 AM] yanhua Huang says: (chuckle)perfect
[11/12/2008 1:20:18 AM] Dan Zwick says: There are a bunch of other ellipse functions in hit_help.cpp.
[11/12/2008 1:20:36 AM] Dan Zwick says: Which I apparently also wrote.
[11/12/2008 1:21:31 AM] Dan Zwick says: In hit_help.cpp the task was to generate points equally spaced (according to arc length) on an ellipse.
[11/12/2008 1:21:36 AM] yanhua Huang says: t_theta = ellipse_arc_inverse(t_theta, s_start + p * (s_end - s_start), M, N);
how to set t_theta?
[11/12/2008 1:22:24 AM] Dan Zwick says: It's an initial guess.
[11/12/2008 1:23:02 AM] Dan Zwick says: If you are starting at the beginning set it to 0. After that set it to whatever the previous theta was.
[11/12/2008 1:23:17 AM] Dan Zwick says: That is, assuming you are going to do it for all points at once.
[11/12/2008 1:24:43 AM] Dan Zwick says: Keep in mind that these ellipse functions assume that the ellipse is in standard position. That is, it has been rotated to the axes align with the x and y axes and the center is at the origin.
[11/12/2008 1:24:58 AM] yanhua Huang says: Which feature should be created if I want to walk into ellipse_angle?
t:\V44B\HITS\HIT_HELP.CPP
[11/12/2008 1:25:35 AM] Dan Zwick says: Also, it might make a difference where the larger semi axis is (along x or along y).
[11/12/2008 1:27:04 AM] yanhua Huang says: max axis set to x?
[11/12/2008 1:27:20 AM] Dan Zwick says: I'm not sure.
[11/12/2008 1:29:03 AM] yanhua Huang says: thanks again, I will trace these code firstly tomorrow.
[11/12/2008 1:29:31 AM] Dan Zwick says: Ok. M and N are the semi-axes, but I'm not sure which is which.
[11/12/2008 1:30:37 AM] Dan Zwick says: theta = atan2(M * sin(t), N * cos(t))
[11/12/2008 1:31:58 AM] Edited by Dan Zwick, 11/12/2008 1:32:22 AM It looks like the equation of the ellipse here is
x = N cos(t), y = M sin(t), or (x/N)^2 + (y/M)^2 = 1.
[11/12/2008 1:33:33 AM] Edited by Dan Zwick, 11/12/2008 1:32:22 AM So N is along x, M is along y, but that doesn't necessarily mean that N > M.
[11/12/2008 1:34:13 AM] yanhua Huang says: where is the comment?
[11/12/2008 1:34:41 AM] yanhua Huang says: got it.
[11/12/2008 1:35:09 AM] Dan Zwick says: Oh, there's also this:

// The following code is used by gen_ellipse_hits. It computes the polar angle of
// the point on the ellipse a proportion p of the way from the point on the perimeter
// corresponding to the polar angle start to the point on the perimeter corresponding
// to the polar angle end.
//
// NOTES.
// 1. The points on the ellipse are given in terms of the ellipse parameter t by
// (x,y) = (M cos(t), N sin(t)). The polar angle theta corresponding to (x,y) is
// arctan(y/x), which gives the conversion formula theta = arctan(N sin(t) / M cos(t)).
// 2. Arclength on an ellipse is expressed as an elliptic integral (of the second
// kind) with respect to the ellipse parameter t, which may be estimated to a high
// degree of accuracy. Once the arclength from start to end is computed, this integral
// has to be inverted in order to find the parameter t corresponding to the arclength
// of the desired point. This is done with Newton's Method.
[11/12/2008 1:35:42 AM] Dan Zwick says: This is in hit_help.cpp.
[11/12/2008 1:36:20 AM] Dan Zwick says: And that is definitely my English!
[11/12/2008 1:38:09 AM] yanhua Huang says: (party) are there other notes?
[11/12/2008 1:38:37 AM] Dan Zwick says: Just in the code.
[11/12/2008 1:39:21 AM] Dan Zwick says: Funny I don't remember any of this. It was only a year ago. :S
[11/12/2008 1:40:51 AM] yanhua Huang says: I try to look for the function using keyword = "ellipse" this afternoon,
[11/12/2008 1:43:36 AM] yanhua Huang says: got >4000 matching lines results.
[11/12/2008 1:44:59 AM] Dan Zwick says: Oh. I searched for elliptic integral.
[11/12/2008 1:46:14 AM] yanhua Huang says: (F)I will go to |-)bed, good day,bye.
[11/12/2008 1:46:36 AM] Dan Zwick says: Good night!
[12:19:16 AM] yanhua Huang says: hi in?
[12:19:35 AM] Dan Zwick says: Yup
[12:20:30 AM] yanhua Huang says: // Find the length of arc on an ellipse from ellipse parameter phi1 to
// ellipse parameter phi2.
//
double ellipse_arc(double phi1, double phi2, double a, double b)
is this function calculating arc lenth from phi1 to phi2 with CCW direction?
[12:20:33 AM] Dan Zwick says: (that's "cowboy talk" for yes)
[12:21:03 AM] Dan Zwick says: Let me check...
[12:22:00 AM] Dan Zwick says: Well, I would think that the answer would be positive if phi1 < phi2 and negative otherwise.
[12:22:15 AM] Dan Zwick says: But I can check...
[12:24:23 AM] Dan Zwick says: Also, the points on the ellipse are (M cos(phi), N sin(phi)), so they start at (M,0) when phi = 0 and proceed CCW around the ellipse as phi increases.
[12:25:07 AM] Dan Zwick says: So I guess the answer is yes.
[12:25:44 AM] yanhua Huang says: ok
[12:27:17 AM] yanhua Huang says: I did add ELLIPSE_ARC.CPP and \ELLIPSE_ARC.H into PCDSFIT project.
t:\V44B\PCDSFIT\INCLUDE\ELLIPSE_ARC.H
t:\V44B\PCDSFIT\SOURCE\ELLIPSE_ARC.CPP
[12:28:08 AM] Dan Zwick says: Do they need to be in PCDSFit? They're already in PC-DMIS.
[12:29:21 AM] yanhua Huang says: because add a new function - LengthArc into t:\V44B\PCDSFIT\SOURCE\STRAIGHTNESS2D.CPP
inline double LengthArc(double *Angles, int i, int j, double a, double b)
{
double dbSum = 0;
if (Angles && Angles[j])
{
dbSum = ellipse_arc(Angles, Angles[j], a, b);
}
return dbSum;
}
[12:30:27 AM] Edited by Dan Zwick, 12:30:37 AM But ellipse_arc is already accessible in MATHTOOL, isn't it?
[12:31:34 AM] yanhua Huang says: which project is mathtool in?
[12:31:59 AM] Dan Zwick says: pcdlrn
[12:39:40 AM] yanhua Huang says: (sweat)I uploaded my changed code just now.
[12:40:17 AM] Dan Zwick says: Does it seem to be working?
[12:42:07 AM] yanhua Huang says: I have not tested per unit profile with ellipse, because have not a PP including ellipse. Don R is traveling
[12:42:40 AM] Dan Zwick says: Ok. Well, good luck.
[12:46:25 AM] yanhua Huang says: but I used the same method to slide per unit window for per unit profile with circle. I am testing the deviation result....now, because Porter is dealing with some issues with me, I need download his latest code, so I uploaded my code just now.
[12:47:00 AM] yanhua Huang says: set major axis = circle diameter, minor axis = circle diameter
[12:47:31 AM] Dan Zwick says: I hope you're not using the ellipse code for the circle.
[12:48:17 AM] yanhua Huang says: why? just e = 1 for circle
[12:48:57 AM] Dan Zwick says: The ellipse code will be very inefficient for the circle. The circle is a degenerate case for the ellipse, so there could also be precision problems.
[12:50:03 AM] Dan Zwick says: It's trivial to compute arc length on a circle (r * theta) and it is trivial to find closest points on a circle. Also, it is not necessary to transform it to standard position. It would be much cleaner and much simpler to do the circle separately.
[12:51:25 AM] yanhua Huang says: (chuckle)Now I did transform all meas points to standard position, because need sort these meas points.
[12:52:16 AM] Edited by Dan Zwick, 12:52:34 AM Well, it's ok to shift the center to the origin, but what about trying to rotate it? The two semiaxes are equal.
[12:54:46 AM] yanhua Huang says: did not rotate
[12:55:00 AM] Dan Zwick says: There's a function in parallel.cpp that sorts points by central angle.
[12:55:19 AM] yanhua Huang says: I agree to do the circle separately
[12:55:26 AM] Dan Zwick says: Ok
[12:56:16 AM] Dan Zwick says: Well, the only advantage of that sort routine is that it doesn't compute the actual angle. But you need the angles anyway.
[12:57:50 AM] Dan Zwick says: So if the circle is at the origin, then the angle for a point is theta = atan2(y,x) and its closest point is (r cos(theta), r sin(theta).
[12:58:40 AM] Dan Zwick says: And the arc length from (r,0) to the point is just r * theta (in radians).
[12:59:45 AM] yanhua Huang says: I did not calculate closet points, just calculate arc length between two polar angles relative to meas points
[12:59:53 AM] Dan Zwick says: One thing you could do is simply add a case for when a == b (or fabs(a -b) < epsilon).
[1:00:13 AM] Dan Zwick says: Ok. The closest points are implicit.
[1:05:08 AM] yanhua Huang says: got it. Thanks
[1:05:15 AM] Dan Zwick says: yw