1. ## Backdooring Instantaneous Radius of Curvature & Functions

I'm a novice coder who seeks help w/ creating a function and possible clean up of code.

Brief history: I'm currently working on a project that requires the use of a C++ DLL to talk to a PLC machine running in a real time production environment. My background is mathematics, not coding. I'm struggling to make this work!

Currently the parameters for this machine require evaluating data by instantaneous radius of curvature standards. (Recall that the smaller the radius of curvature, the more curvy or less straight a curve becomes. This is not difficult to find by way of calculus, but I'm nowhere near ready to code derivatives and second derivatives! My current curve needs to be adjusted to fit the minimum radius parameters. I've found a way of backdooring this calculus requirement with a little geometry and iterating along the curve, looking at three adjacent points each time. (Points are close enough to be considered instantaneous)

Currently this code works. I had to do some strange internal error handling because it was coming up with strange results when I got close to infinity, 1/0 or something close to it. Maybe this could be a corrected suggestion instead of my if, else stmts.

1) Need to reuse this method of finding the radius several more times than currently using it, so need it to be a function instead of part of the main body of code! Not sure how the variables inside of the function need to be set up to conserve memory

2) Recently read an article warning against function use in DLLs because of inability to handle errors. (An error would stop production Not a good thing! Should I be worried about this?) This DLL called is an exported function---Does this mean I can't call a function within a function? I know nothing

3) Can you suggest any restructuring or changes to be made

Here's the code & structure that's currently used:
Code:
```DLLMain()

//stuff to make the PLC machine & DLL talk to each other

extern "C" __declspec(dllexport) int calculate(					EXT_ROUTINE_CONTROL* pERCtrl,
int npar[] , float fpar[], float xcl[], 			float xdia[], float ycl[], float ydia[], 			float zlen[], int nresults[], float 			fresults[])
{
//code
//code
//more code!
//radius of curvature algorithm formatted 			differently
//code
}```
Changes made to the program will require the radius of curvature algorithm to be repeated in a binary search fashion. Repeating that block of code would be extremely unfriendly to the reader.

Code:
```//function prototype
float FINDRADIUS( float& y1, float& z1, float& y2, float& z2, float& y3, float& z3);

DLLMain()

//stuff to make the PLC machine & DLL talk to each other

extern "C" __declspec(dllexport) int calculate(					EXT_ROUTINE_CONTROL* pERCtrl,
int npar[] , float fpar[], float xcl[], 			float xdia[], float ycl[], float ydia[], 			float zlen[], int nresults[], float 			fresults[])
{
//code
//code
//more code!
//call function several times!!!!!!
//code
}

//???????????????????How & where should I declare variables

float FINDRADIUS( float& y1, float& z1, float& y2, float& z2, float& y3, float& z3)
{
yDiff1=y2-y1;
yDiff2=y3-y2;

absYDiff1 = abs(yDiff1);
absYDiff2 = abs(yDiff2);

/*NOTE:the following if, else if tests are important because
division by zero will be a problem since the two lines constructed
have negative reciprocal slopes relative to the ydiff's*/

if (abs(yDiff2-yDiff1)<=.0002)
{
/*log is straight along this interval (i.e--do NOTHING jump code
set radius[j] really high [not important] and return*/
goto Straight_Log;
}//end if straight interval

if (absYDiff1<=.00015)
{
/*To avoid dividing by zero, set absolute value of
m2 higher than 20000 (upper limit check ok, jf 7/2004)*/
if (yDiff1>0)
m1=20000;
else
m1=-20000;
m2= -(z3-z2)/yDiff2;
}//end if yDiff1 is close to 0

else if (absYDiff2<=.00015)
{
/*To avoid dividing by zero, set absolute value of
m2 higher than 20000.  (identical upper limit as above, jf)*/
if (yDiff2>0)
m2=20000;
else
m2=-20000;
m1= -(z2-z1)/yDiff1;
}//end else if yDiff2 is close to 0

if (absYDiff1 > .00015 && absYDiff2 > .00015)
{//finds the perpendicular slope of the two intervals
m1= -(z2-z1)/(y2-y1);
m2= -(z3-z2)/(y3-y2);
}//end if perpendicular slopes

else if (m1==0 && m2==0)
{//error check    (should never hit this)
m1= -(z2-z1)/(y2-y1);
m2= -(z3-z2)/(y3-y2);
}//end else if error check

Z1Mdpt = (z2+z1)/2;
Y1Mdpt = (y2+y1)/2;
Z2Mdpt = (z2+z3)/2;
Y2Mdpt = (y2+y3)/2;

/*We know the slope of the perpendicular bisecting lines
and a point on each of the lines, the midpoints.  Since
the two lines intersect at the radius, solving the system
Try w/ a compass if you don't believe me ;)  */

Zradius = (m1*Z1Mdpt - Y1Mdpt - m2*Z2Mdpt + Y2Mdpt)/(m1-m2);

/*The coordinate of the radius is known.  By using the
distance formula we can find the length of the radius
passing through the three points*/

2. //???????????????????How & where should I declare variables
Just declare the variables you need in FINDRADIUS, within that function.
You are safe calling another function like that, no need to worry.

3. Some observations:
1. What happens if yDiff2 and yDiff1 are BOTH tiny? Say yDiff1 = 0.00014, yDiff2 = -0.00007. Then
1. in the following code, |yDiff1-yDiff2| = 0.00021 > 0.0002, and you will not get rescued by your "goto Straight_Log".
Code:
```if (abs(yDiff2-yDiff1)<=.0002)
{
/*log is straight along this interval (i.e--do NOTHING jump code
set radius[j] really high [not important] and return*/
goto Straight_Log;
}//end if straight interval```
2. The second point is that the goto has you jumping out of the FINDRADIUS function, apparently never to come back, and so FINDRADIUS will never return a value. If you have called FINDRADIUS in some code like
Code:
`r = FINDRADIUS(parameters happen to fit straight log;`
the value of r will be undefined, and control of the program may have gone somewhere totally undesired.

Better to make a function out of whatever code you have following the Straight_Log label, and then have FINDRADIUS return the result of that function.
Code:
```/* in function FINDRADIUS */
if (condition indicating straight log) {
return Straight_Log(whatever arguments are needed);
}```
3. Now in the next bit, you check to see if yDiff1 is small, but then divide by yDiff2 without checking to see if that is ALSO small:
Code:
```if (absYDiff1<=.00015)
{
/*To avoid dividing by zero, set absolute value of
m2 higher than 20000 (upper limit check ok, jf 7/2004)*/
if (yDiff1>0)
m1=20000;
else
m1=-20000;
m2= -(z3-z2)/yDiff2;
}//end if yDiff1 is close to 0```
With the yDiff1 and yDiff2 given above, you will have m2 = -(z3-z2)/(-0.00007). Not only that, but the code to take care of tiny yDiff2 will get skipped over. If that's ok, then nevermind. But it would still be clearer to write your if/else if to test for BOTH yDiff1 and yDiff2 together:
Code:
```if      (yDiff1 <  TINY && yDiff2 <  TINY) {do something}
else if (yDiff1 >= TINY && yDiff2 <  TINY) {do something else}
else if (yDiff1 <  TINY && yDiff2 >= TINY) {do yet again something else}
else if (yDiff1 >= TINY && yDiff2 >= TINY) {the nice case}```
4. Rather than having magic numbers like 0.00015, 0.0002, 20000 littered all over the place, better to make consts out of them. That way, if you ever want to change the values, you'll only have to do it in one place.

2. In your error check section,
Code:
```else if (m1==0 && m2==0)
{	//error check    (should never hit this)
m1= -(z2-z1)/(y2-y1);
m2= -(z3-z2)/(y3-y2);
}	//end else if error check```
m1 and m2 could both be zero if the z's are all equal, and all of the checks you've done so far have only been about the y's. And even if m1 and m2 are not both zero, they might still be equal, so your Zradius formula could end up dividing by zero.

Now a mathematics question: in your coordinate system, is y the horizontal axis, & z the vertical? In the (y,z) coordinate system, I thought delta_z/delta_y is the slope of the line connecting (y1,z1) & (y2,z2), and -delta_y/delta_z is the slope of the line perpendicular to that? Anyway, just check to be sure the numerators & denominators aren't switched in your formula.

3. If the right way to do this is by calculating derivatives (e.g., if y = f(x), then radius(x) = (1+(f')^2)^(3/2) / f'' = exp(1.5 * log(1+sqr(f')))/f''), then there are some pre-packaged routines for calculating derivatives. Check out the numerical recipes in C site: http://www.nr.com, or http://www.library.cornell.edu/nr/cbookcpdf.html (e.g. chapter 5.7). There might be something that will work for you.

4. I sure do appreciate the time that you've spent analyzing the code.

1.)(a) you will not get rescued by your "goto Straight_Log"

The boat comes to rescue me in the if statement following the "goto Straight_Log" if statement. However, it sails away laughing. It only saves me most of the time. (epiphany came after reading part (c) of your suggestions)

(b) goto has you jumping out of the FINDRADIUS function, apparently never to come back

I've tried to modify an existing routine and somehow forgot that line. Thanks!

I will make changes to incorporate the Straight_Log function

(c) Referring to if yDiff2 and yDiff1 are BOTH tiny

You're correct that this if, then situation will not catch an error.

Hypothetically speaking the yDiff1=0.00015 will be acceptable, but the yDiff2=-0.00007 could give me problems. I ran some tests and with both yDiffs so close to zero the intersecting lines intersect in some far off distance of no importance to my analysis. Now abs(yDiff1-yDiff2)=0.00021 and 0.00021>.0002, which causes the code to miss the Straight_Log parameters. But the very next if (absYDiff1<=.00015) would catch the small yDiff1, throw (+ or -)20,000 into the slope, m1, and then find the slope for m2---------------I'm seeing your point (light bulbs ? ? ?)

You are sooo right, this will not be caught by goto statement, but will be caught by the first if yDiff1<=0.00015. It will pick 20,000 for the first slope, but calculate the slope as if everything were fine and dandy for the second slope. I can catch that with a nested if, (right again)

(d) make consts out of them

Will do.

2) m1 and m2 could both be zero if the z's are all equal

In this situation the z's will never be zero. This is in a scanning environment with everything based upon movement. As the data is scanned it is set up to take a measurement around every three inches. Even if production stops, the scanners will wait untill movement resumes to take a measurement to ensure that no z's are the same. (ie: the measurements taken will only be recorded as the z increases incrementally b/c it's motion based scaning)

even if m1 and m2 are not both zero, they might still be equal, so your Zradius formula could end up dividing by zero

I probably need a check for this!

Theoretically if the m1 and m2 are equal this should have been caught by the straight log scenario. The m1 and m2 being equal should only occur if the data is straight or, mathematically, has a constant slope along the interval. In this scenario, since the delta z's can never be zero, it would indicate that the yDiff1=yDiff2 (approximately). Again, theoretically, this should have been caught by the first If statement, if (abs(yDiff1-yDiff2)<=0.0002). But your right too many theoretical webs weave trouble or potential trouble!

3) I'll check it out!

Problem
with this is that the DLL currently fitting data through process of polynomial regression doesn't provide coefficients of the equation. So f' and f'' had to be skirted for now. I recently found one that yields coeff's and am considering making some changes in the near future to make things somewhat more readable to an outside coder. I might just write my own if I get brave.

Have you ever messed around with a Vadermonde Matrix?

Have a great day! Hope you get to read my response so you'll know how much I appreciate your time!

5. I could be wrong but I believe there is a faster way using spherical coordinates.

Any 2D sphere can be setup as:

x=sin(beta)*sin(alpha)
y=sin(beta)*cos(alpha)
z=cos(beta)

To map from spherical to rectangular.

and

0.0f<x<1.0f
0.0f<y<1.0f

For ovals you simply need to add a coefficient into the mix on the desired axis.

You can find the radius if you know two points on the curve. Find the point at which the line through point 1 intersects the line through point 2. Now you know the center.

The normal N at point P for any curve centered about point C is N=P-C since all points on all curves are equidistant from the center.

Also you can create a right triangle out of two known points on the curve and thus using simple trigonometry can then find the angles you need and/or the radius.

I'm sort of confused as to why all the other math is involved. I'm working with radius and curves in my game and I'm not sure what you are getting at.

6. Originally Posted by Bubba
I could be wrong but I believe there is a faster way using spherical coordinates.

Any 2D sphere can be setup as:

x=sin(beta)*sin(alpha)
y=sin(beta)*cos(alpha)
z=cos(beta)
Why would spherical coordinates be used in a 2d plane with z being the horizontal axis and y being the vertical axis in this case?

Do you mean cylindrical? (or am I just confused)

Originally Posted by Bubba
You can find the radius if you know two points on the curve. Find the point at which the line through point 1 intersects the line through point 2. Now you know the center.
The code above does use 2 points. In fact the two points used are the midpoints of the three adjacent points (z1,y1) (z2,y2) (z3,y3). Why in the world would I use this??? Good question. I am changing the coordinates of the center point based upon my instantaneous radius of curvature parameters. If it passes the curvature requirements, I'll move on. If it fails , I'm flagging the center data point as being unacceptable. If it's flagged then I plan on moving (z2,y2) vertically, flattening the curve until the curvature is acceptable. I have a method in mind that will get the right answer but I'm affraid that the overhead required to do so will be unacceptable

Originally Posted by Bubba
Also you can create a right triangle out of two known points on the curve and thus using simple trigonometry can then find the angles you need and/or the radius.

I'm sort of confused as to why all the other math is involved. I'm working with radius and curves in my game and I'm not sure what you are getting at.
I think that I'm going to look at this from your triangle perspective a little more. Maybe I can restrict the angle (opposite the height when increasing or equal to the height when decreasing) of the right triangle. (I will always view this from Left to Right) but I'm going to need to relate that to a circle or radius.

Sorry hitting a blank. Suggestions?

7. Originally Posted by just2peachy

You are sooo right, this will not be caught by goto statement, but will be caught by the first if yDiff1<=0.00015. It will pick 20,000 for the first slope, but calculate the slope as if everything were fine and dandy for the second slope. I can catch that with a nested if, (right again)
Rather than nested ifs, I think it will read easier to test yDiff1 and yDiff2 together. In other words, I think

Code:
```if (absyDiff1 < TINY) {
if (absyDiff2 < TINY) {
do something 1
blah blah blah
}
else if (absyDiff2 >= TINY) {
do something 2
blah blah blah
}
}
else if (absyDiff1 >= TINY) {
if (absyDiff2 < TINY) {
do something3
blah blah blah
}
else if (absyDiff2 >= TINY) {
do something4
blah blah blah
}
}```

Code:
```if      (absyDiff1 <  TINY   &&   absyDiff2 <  TINY) {
do something 1;
blah blah blah
}
else if (absyDiff1 <  TINY   &&   absyDiff2 >= TINY) {
do something 2;
blah blah blah
}
else if (absyDiff1 >= TINY   &&   absyDiff2 <  TINY) {
do something 3;
blah blah blah
}
else if (absyDiff1 >= TINY   &&   absyDiff2 >= TINY) {
do something 4;
blah blah blah
}```
Anyway, best of luck on the project!

8. Then what you really need to do is look into NURBS or non uniform rational b-splines or bezier curves. They are essentially controlled by control points some distance away from the curve. As you move the control points the curve changes. Very complex curves can be created with only a few control points.

Google for NURBS, splines, quaternion splines (squad), spherical linear interpolation (SLERP), and bezier curves.

I'm not sure how your particular device is retrieving coordinates to compute instantaneous radius of curvature so I'm quite lost as to the process involved. But I do understand curves and beziers....but NURBS get a little complex.

You can also use spherical linear interpolation, but to use that you must dump the euler angles and move to quaternions. There is an excellent article about quaternions on this board.

But much of what you are doing can be accomplished with linear and bilinear interpolation.

I will use the double data type here since precision is important.

Code:
```double LI(double v1,double v2,double f1)
{
return v1+f1*(v2-v1);
}

double BI(double v1,double v2,double v3,double v4,double f1,double f2)
{
double value1=v1+f1*(v2-v1);
double value2=v3+f1*(v4-v3);
return value1+f2*(value2-value1);
}

void SLERP(const Quaternion q0,const Quaternion q1,float slerp_factor,Quaternion &qout)
{
float cosOmega=((q0.w*q1.w) + (q0.x*q1.x)+(q0.y*q1.y)+(q0.z*q1.z));

if (cosOmega<0.0f)
{
//q1.Negate()
q1.w=-q1.w;
q1.x=-q1.x;
q1.y=-q1.y;
q1.z=-q1.z;
cosOmega=-cosOmega;
}

float k0=0.0f,kt=0.0f;

if (cosOmega>.9999f)
{
k0=1.0f-slerp_factor;
k1=slerp_factor;
}
else
{
float sinOmega=sqrtf(1.0f-cosOmega*cosOmega);
float omega=atan2(sinOmega,cosOmega);

float OneOverSinOmega=1.0f/sinOmega;

k0=sin((1.0f-slerp_factor)*omega)*OneOverSinOmega;
k1=sin(slerp_factor*omega)*OneOverSinOmega;
}

qout.w=q0.w*k0+q1.w*k1;
qout.x=q0.x*k0+q1.x*k1;
qout.y=q0.y*k0+q1.y*k1;
qout.z=q0.z*k0+q1.z*k1;
}

void SQUAD(Quaternion q0,Quaternion q1,Quaternion q2,Quaternion s1,Quaternion s2,float factor)
{
float new_factor=(2.0*factor)*(1.0-factor);
return SLERP(SLERP(q0,q1,h),SLERP(s1,s2,h),new_factor);
}```
The last two functions are taken from 3D Math Primer for Graphics and Game Development by Fletcher Dunn and Ian Parberry. For an explanation of how it works or how to implement your own, consult their book or buy a book about 3D mathematics.

SQUAD is quaterion spline interpolation and SLERP is spherical linear interpolation. Given the correct control points any curve can be generated using these. There is also math to learn for finding control points on a curve and/or creating control points that two points on a curve would share.

Video games uses these formulas and one like it to smoothly interpolate around objects at any angle and orientation. Very difficult with euler angles.

You may also want to look into the equations for closest point on a 2D line to some point in space, intersection of two lines, etc.

Here is some math for spheres and circles.

Implicit form: || p-c || = r
Essentially what I said earlier.

Implicit Equation:
((x-cx)*(x-cx))+((y-cy)*(y-cy))*((z-cz)*(z-cz))=r*r;

circumference=PI*diameter;

The derivative of the area of a circle is its circumference.
The derivative of the volume of a sphere is its surface area.

9. Originally Posted by Bubba
Then what you really need to do is look into NURBS or non uniform rational b-splines or bezier curves. They are essentially controlled by control points some distance away from the curve. As you move the control points the curve changes. Very complex curves can be created with only a few control points.

Google for NURBS, splines, quaternion splines (squad), spherical linear interpolation (SLERP), and bezier curves.

I'm not sure how your particular device is retrieving coordinates to compute instantaneous radius of curvature so I'm quite lost as to the process involved. But I do understand curves and beziers....but NURBS get a little complex.
I'm currently reading several articles about the Bezier curves and am working on thier implementation. However, I don't fully understand them. In fact, it's a little scary to me. The articles that I've collected are not getting the job done. Maybe it's just me .

Thanks for recommending a specific book. I plan on ordering it & checking out quaternion board discussion. If you have found any other manuals that "break it down" clearly, please drop some titles my way when you have some time.

Again. . . Thank you.