Script for Curvature Calculation

Next
 From:  Karsten (KMRQUS)
6634.1 
Hello,

I want to share my first script. It's dirty but not quick. I hope that someone will test it and compare the results with another software for a validation. Straight line segments can lead to problems, especially in radius-calculations. The performance is very slow (wait - it's not crashed), so I would like to get tipps for an acceleration. Depending on the dimensions of your geometry, you should adjust the scalefactor - sometimes the spikes are so small you can't see them. The curvature at start- and endpoint won't be calculated. So my question is: Can I get a tangent-vector from a curve on a special point or points of curve to calculate it approx by script, without arraycurve-factories?

Kind regard
and thanks to all that have shared they scripts:)
Karsten
Attachments:

Image Attachments:
Size: 432.5 KB, Downloaded: 250 times, Dimensions: 1366x768px
  Reply Reply More Options
Post Options
Reply as PM Reply as PM
Print Print
Mark as unread Mark as unread
Relationship Relationship
IP Logged

Previous
Next
 From:  bemfarmer
6634.2 In reply to 6634.1 
Hi Karsten

By putting a lot of the mathematics in the Script section of the .htm file, a very high increase in speed is / may be possible.
In order to work, it has to be done in certain ways. Not everything works out.
There are several "recent" scripts with examples of how it can be done.
There is recent discussion about it. If something goes wrong, recovery is limited.

Do you have any mathematic papers for reference on the equations?

(I googled a few "porcupine" and "comb" curvature analysis searches last week.)

- Brian


Now to try out your script.
  Reply Reply More Options
Post Options
Reply as PM Reply as PM
Print Print
Mark as unread Mark as unread
Relationship Relationship
IP Logged

Previous
Next
 From:  Karsten (KMRQUS)
6634.3 In reply to 6634.2 
Hello bemfarmer,

no, I haven't any papers of the calculation. The idea is to calculate the Center of an arc given by 3 Points - the length is the radius - 1/R is the curvature. So I try to explain what the script should do:
In a first step it gets an Array of Points lying on the curve. Then it calculates the direction-vectors from the first to the second an from second to third Point. Half them and add the startpoint, you will have a Point for a plane. The direction-vector is the perpendicular vector of the plane. Then you can create a third plane by make a cross-product (axb) of the two normal-vectors (direction vectors). The intersection of the 3 planes is the centerpoint of the arc. You can make an Linear Equations System from the components of the 3 given planes. It is solved with cramers rule (implementin a gauss algorithm in Java-script isn't possible for me at this time). The determinants are solved with saurrus rule. If you have the Center you can calculate the Spikes with some vectormath.

Let me know how it works!

Kind regard
Karsten
  Reply Reply More Options
Post Options
Reply as PM Reply as PM
Print Print
Mark as unread Mark as unread
Relationship Relationship
IP Logged

Previous
Next
 From:  Michael Gibson
6634.4 In reply to 6634.1 
Hi Karsten, great job so far!

> Can I get a tangent-vector from a curve on a special point or points of curve
> to calculate it approx by script, without arraycurve-factories?

Sorry no unfortunately there isn't any methods set up currently for scripts to evaluate properties of curves like tangent vectors at arbitrary points.

That is something that I want to add in the future - right now the scripting in MoI is not focused on doing heavy duty calculations inside of a script, the scripting is primarily used at this time as a kind of control flow mechanism that handles hiding/showing UI, creating picking mechanisms, things like that. So it's currently missing a lot of stuff that would be useful for more detailed calculations like you're trying to do here.

In the future I definitely want to make that sort of thing more feasible to do in script.

- Michael
  Reply Reply More Options
Post Options
Reply as PM Reply as PM
Print Print
Mark as unread Mark as unread
Relationship Relationship
IP Logged

Previous
Next
 From:  Karsten (KMRQUS)
6634.5 In reply to 6634.2 
Hi Brian,

thank you for your tipps about increasing performance. I've took a look in the discussion about inline programming - but it is to complicated for me. I don't know which functions can be put in and which not - and how they interchange data with a mainsript. (Sorry about my englisch and my skills in java script). The only performance killer I found, is that i Calculate every plane twice. I have to check, how much work it is, to change that.

Futher more I think my explanation about the underling idea wasn't the best, so I've draw a model (in Attachment) - a picture is better then thousands of words.

Greetings
Karsten

EDITED: 15 Jun by KMRQUS

  Reply Reply More Options
Post Options
Reply as PM Reply as PM
Print Print
Mark as unread Mark as unread
Relationship Relationship
IP Logged

Previous
Next
 From:  bemfarmer
6634.6 In reply to 6634.5 
Hi Karsten

Maybe I'll have a look this weekend.

(I could not get PDF save of your ccc.3dm to print. It just did outline(?)
So I just get a screensave.)

I'm wondering, if the curve is planar, why not set it as cplane, and intersect the perpendicular bisector lines through the
two centerpoints of P1P2 and P2P3, to find the center of the circle from AP1, P2, P3? Is it necessary to use 3 planes?

- Brian
  Reply Reply More Options
Post Options
Reply as PM Reply as PM
Print Print
Mark as unread Mark as unread
Relationship Relationship
IP Logged

Previous
Next
 From:  Karsten (KMRQUS)
6634.7 In reply to 6634.6 
Hi Brian,

maybe we're talking about different things. The script don't create planes in the geometry database and intersect them there. Also it calculates not the curvature of one single point. (The easiest way to do this, is to draw an arc given by 3points an the secelct it or snap to the center) It calculates the curvature of all samplepoints and displays the "flow" of the curvature as a visualisation of "spikes" with an envelope curve. The curve is not limited to 2D. I don't know if it is possible to adjust a cplane in script, nevertheless you have to calculate the orientation (the third plane!). Because the script calculates the center complete in a mathematical way, a flat curve in the xy-plane will be automatic reduced to a 2D problem. The planes are given as:

A1x+B1y+C1z+D1=0
A2x+B2y+C2z+D2=0
A3x+B3y+C3z+D3=0

for a curve in xy-plane-> z=0->
with C1z=0; C2z=0;A3x=0;B3y=0;D3=0 ->
so you get a linear equation system like that:

A1x+B1y+D=0
A2x+B2y+D=0

a 2D dimensional problem!
Ok in this special configuration it makes some useless calculations, but the atvantage is the general formulation.
The function of the script is to analyse continuity of curves (G2-Stuff). You can find problems of swinging curves, continuity, changes in curvature direction. I used such a feature (with some other) in CATIA V4/V5 and UG NX to analyse, design/reconstuct and repair CLASS A surfaces. And I missed it in Moi. The script is tested with Moi V3 Beta. It doesn't work with V2. It seems there is a problem in the HTML-Call
<moi:NumericInput id="scalefactor" default="0.1" style="persist:true">

Have an nice day!

Kind regards
Karsten
Image Attachments:
Size: 198.6 KB, Downloaded: 140 times, Dimensions: 1366x768px
  Reply Reply More Options
Post Options
Reply as PM Reply as PM
Print Print
Mark as unread Mark as unread
Relationship Relationship
IP Logged

Previous
Next
 From:  bemfarmer
6634.8 In reply to 6634.7 
Hi Karsten

Thank you for the information. I'll study on it after work :-)
I need to understand it better, as well as your code.

Your script makes some very nice curves, even 3D.
I did a loft from a ccc of a clothoid curve to the moved clothoid, for a nice surface.

Yesterday I worked on moving some of the code to htm, for a speedup.
I tried switching from factory/factories to object lists, and even though
I used some of Max's type code, I still got a "type mismatch," so I am doing something wrong.
It is puzzling.

- Brian

Only one thing gets returned. A factory, an array, or an object list.
So I hope to separate out lFactories, acFactory, and pLFactory.

lFactories holds the lines. (?)
acFactory is the array of points. (?)
plFactory is the interpcurve made up of line end points. (?)

EDITED: 21 Apr 2014 by BEMFARMER

  Reply Reply More Options
Post Options
Reply as PM Reply as PM
Print Print
Mark as unread Mark as unread
Relationship Relationship
IP Logged

Previous
Next
 From:  bemfarmer
6634.9 
Outline of proposed modifications to Karsten's CCC script, for faster htm processing.

(All constructive criticism and comments welcome :-)

Note that the .htm script runs in a different (lesser / less escape-able/ ...) processing environment than .js script code.

By bemfarmer April 21, 2014.

After numerous trials, as far as I can tell, only one objectList, or factory, or array,
can be returned by a function in .htm script, to the .js script. (True for javascript in general...)
If additional such items or geometric objects are created in the .htm script, they are not recognized in the .js script, for update purposes. (Not sure about global variables?)

So the three factories, (or their results), acFactory, lFactory, and plFactory, must have their associated .htm functions called separately. Four main .htm functions are proposed in the following outline, with the first two returning point arrays, the third returning lines object list, and the last an interpcurve curve.

1. function getacpoints( Curve, NumItems )
returns ptx, an array of points on the curve, using acFactory.
These points will be use later, in the calculations.
Also, excepting the first and last point, these points will be used later, as the start points of the lines which represent curvature K, or radius.

2. function getplpoints( ptx, analyse_as, scalefactor )
returns plpoints, an array of the scp points, by calculations.
These points will be used later, as the end points of the lines, and also for an interpcurve.

3. function getLines( ptx, plpoints )
returns linesList, an object list of the lines. (or else return the line array.)

4. function getplcurve( plpoints )
returns plcurve, an interpcurve. (Or objectList...)

Note that the clone of ptx, the pts array, is not being used, as far as I can tell.

- Brian

EDITED: 22 Apr 2014 by BEMFARMER

  Reply Reply More Options
Post Options
Reply as PM Reply as PM
Print Print
Mark as unread Mark as unread
Relationship Relationship
IP Logged

Previous
Next
 From:  bemfarmer
6634.10 
For CCC of a square or pentagon, an assert failed error is occurring.

- Brian
  Reply Reply More Options
Post Options
Reply as PM Reply as PM
Print Print
Mark as unread Mark as unread
Relationship Relationship
IP Logged

Previous
Next
 From:  Karsten (KMRQUS)
6634.11 In reply to 6634.10 
Hello Brian,

I'm happy to hear that the script is interesting for you and that you think it is worth, to invest your time to improve it:-)
About the factories You are completly right. May I have a look in your actual code to understand what you are doing now?
The problem with square and n-gons may caused by the straight lines. Straight lines leads to colinear sampling points. That leads to two problems. the bisection planes between them are parallel. The intersection line is a tangent to our known universe;-) The second problem is that the normal vector (Hesse Normal Form) will be created by the cross product from these colinear (parallel) vectors. The result is a Null-Vector (magnitude and all components are Zero). That leads to a singulary matrix. I tried to catch that in the script with:
if(DetA!=0)
{
sp.x=DetX/DetA;
sp.y=DetY/DetA;
sp.z=DetZ/DetA;
} else {
sp.x=p1.x;
sp.y=p1.y;
sp.z=p1.z;
}
return sp;

The problem is, that I tried to return a wrong centerpoint as a workaround for straight line sections in the curve. (Something shoud be returned:-()
Later the length for the curvature vector will be calculated. Normally the vector to the Radius is diffrent from zero. Because I return p1, in that situation it is zero - an so it's magnitude, and what have I done - shit:
var d_vec=moi.vectorMath.createPoint(ptx[n+1].x-cp.x,ptx[n+1].y-cp.y,ptx[n+1].z-cp.z);
var magn=magnitude(d_vec);
var sf=scalefactor*1000.0;
var sd_vec=moi.vectorMath.createPoint((d_vec.x*sf/(magn*magn)),(d_vec.y*sf/(magn*magn)),(d_vec.z*sf/(magn*magn)));

d_vec is zero -> magn is zero -> and then I made a divide by zero op - maybe that kills your script! But sometimes it works, I can't retrace the problem - maybe caused by floating point precision(Attachment). A quick and dirty workaround would be to add a small value to sp before returning. Better catch the divide by zero or catch the whole calculation of collinar points - I don't know which way is the best. (For a better understanding the expression (magn*magn): the first scales the vector to a normalized one (lenght 1) - the second magn scales it to the curvature.)

Would you post your file with the killer square? Then I can make some experiments at the weekend:-)

I hope my speculations are right!???

Kind regards
Karsten
Image Attachments:
Size: 158.2 KB, Downloaded: 77 times, Dimensions: 1366x768px
  Reply Reply More Options
Post Options
Reply as PM Reply as PM
Print Print
Mark as unread Mark as unread
Relationship Relationship
IP Logged

Previous
Next
 From:  bemfarmer
6634.12 In reply to 6634.11 
Hi Karsten

Most any square, or rectangle will do.
Here is a basic rounded square, which also does "assert failed."

I cannot comment on how or where to add exception handling for problems due to straight line segments,
because I have not studied the math (yet) :-)

- Brian
  Reply Reply More Options
Post Options
Reply as PM Reply as PM
Print Print
Mark as unread Mark as unread
Relationship Relationship
IP Logged

Previous
Next
 From:  Karsten (KMRQUS)
6634.13 In reply to 6634.12 
Hi Brian,

with rounded corners I get the message too. At the weekend I will try to fix it!

Greetings
Karsten
  Reply Reply More Options
Post Options
Reply as PM Reply as PM
Print Print
Mark as unread Mark as unread
Relationship Relationship
IP Logged

Previous
Next
 From:  Karsten (KMRQUS)
6634.14 In reply to 6634.12 
Hi Brian,

it seems that the DetA sometimes isn't zero for linesegments. So, caused by a division near zero, it produces very high values for the centerpoint. That crashes the factories, esspecially for interpcurve. I've try to fix it by limitation of a maximum radius of 10000 and a minimum radius for curvature calculation of 0.002 (csp returns for lines p1 - that leads to a radius of zero -> divide by zero - > high values for the factories). Further more the changes supresses the generation of degenerated lines (startpoint equals endpoint). I've marked the changes in the script with //############. Also I made some test with the results with another tool. You may have a look at the pictures. It seem that other tools seperates connected curves for the analyse. So the precission at the connection points, especially with lines, isn't good.

Have a nice day!

Kind regards
Karsten
Attachments:

Image Attachments:
Size: 119.9 KB, Downloaded: 63 times, Dimensions: 815x814px
Size: 178.5 KB, Downloaded: 36 times, Dimensions: 815x814px
Size: 160.5 KB, Downloaded: 42 times, Dimensions: 617x904px
Size: 155.8 KB, Downloaded: 19 times, Dimensions: 617x904px
Size: 187.6 KB, Downloaded: 43 times, Dimensions: 876x722px
Size: 182.4 KB, Downloaded: 68 times, Dimensions: 876x722px
  Reply Reply More Options
Post Options
Reply as PM Reply as PM
Print Print
Mark as unread Mark as unread
Relationship Relationship
IP Logged

Previous
Next
 From:  bemfarmer
6634.15 
After many difficulties, I finally made a version using some htm code, which may be 10% faster.
It is not much of a speed up. It does show the ongoing point generation, and leaves the start point.
Some code was separated. Some code was moved to htm script.

I'm going to work on a helix script now.

- Brian

EDITED: 27 Apr 2014 by BEMFARMER


  Reply Reply More Options
Post Options
Reply as PM Reply as PM
Print Print
Mark as unread Mark as unread
Relationship Relationship
IP Logged

Previous
Next
 From:  hairykiwi
6634.16 
Hi Karsten,

Thanks for your script - and to Brian too for giving it a speed tweak.

The math of your script looks really interesting. Unfortunately, I don't have enough time right now to get back up to speed on matrices and work my way through your code. If you think the Gaussian elimination method for multiple linear equations would give a significant speed advantage over the Cramer method - and I can see how it might - then the following page offers a ready to use javascript function and may be worth a look:
http://martin-thoma.com/solving-linear-equations-with-gaussian-elimination/

When I ran your script run on the output from my NACA airfoil generating script, it looks very close to what Rhino's own curvature analysis was giving, (I originally started developing my script in RhinoPython). Sorry I don't have consistent screen captures of identical airfoil series to show from MoI and Rhino, but you get the idea:








Cheers,
Hamish

  Reply Reply More Options
Post Options
Reply as PM Reply as PM
Print Print
Mark as unread Mark as unread
Relationship Relationship
IP Logged

Previous
 From:  Karsten (KMRQUS)
6634.17 In reply to 6634.16 
Hi Hamish,

Thank You for the interessting link and I am glad to here, that the script is interessting for you.
yes, I think that a gaussian method could be faster, also for this small matricies. But, I have some Problems at the Moment: Frist - Time(as we all;-). The second enlarge the first Problem. The script has a serious design error, that kills also the performance. The function for calculating the centerpoints (csp) has 3 points as Parameters. So every call, biscectional planes were calculated - for the next call the second plane of the previous calculation is the first of the actual calculation. So I have a lot of redundant plane calcs. To improve the Performance a lot of changes have to be done.

Kind regards
Karsten
  Reply Reply More Options
Post Options
Reply as PM Reply as PM
Print Print
Mark as unread Mark as unread
Relationship Relationship
IP Logged
 

Reply to All Reply to All