Parkinbot rudolf@digitaldocument.de wrote:

Parkinbot,

That was unfair! The accumulated rotation should be compared to Rz(100) and
not with the identity matrix. With:

echo(error=matrixNorm(Rz(100)-Rm));

I got the following Euclidean norm error :

ECHO: error = 4.06721e-15

Using the max norm:

function norm(p) = max([for(pi=p) abs(pi)]);

the error decreased somewhat to:

ECHO: error = 2.88658e-15

And using the sum norm:

function norm(p) = sum([for(pi=p) abs(pi)]);

the error increases to 6.57807e-15 .

Parkinbot rudolf@digitaldocument.de wrote:

The harder step is not to compute the path but to compute the section
rotations along the path. The helicoidal spring is relatively simple but
with other more general path it is not. For instance:

p = [[0,0,0], [20,0,0], [10,30,0], [0,0,0], [0,0,20]];
n= 20;
path = [for(i=[0:n]) Bz(i/n,p,4)];

line(path);

function Bz(u, p, n=3, from=0) =
(n == 1)?
u*p[from+1] + (1-u)p[from]:
uBz(u,p,n-1,from+1) + (1-u)*Bz(u,p,n-1,from);

module line(p,w=0.1)
for(i=[0:len(p)-2])
hull(){ translate(p[i])  sphere(w);
translate(p[(i+1)%len(p)]) sphere(w); }

[image: Bz4.PNG]

Ronaldo wrote

Ups, sorry, you are right, I didn't see that in my quick review of the code.
I wasn't expecting such a good result and still have my doubts with respect
to applying rotations to a polygon again and again.

Are you sure that it is enough to calc the error of R^m to get a measure of
the error your code will produce?
Am I wrong that in code you calculate a sequence of relative rotations and
translations and apply them to the polygon by using the result of step n as
operand for step n+1?
That means you apply rotations R_i for a polygon P over and over: P1 =
(((PR1)R2)R3)R4 ... So one should compare P1 with P2 = P(R1R2R3R4
... ) which has the problem that R1R2R3*R4 can have a bad condition as
well as any matrix R1, R2 ...

Presuming that your sweep has many steps, than the relative rotation angle
gets smaller and smaller. This  is good news for the condition of the
relative rotation matrices. So you are mainly left with a possible bad
condition of the first rotation R1, which you could split ...

--
Sent from: http://forum.openscad.org/

Parkinbot rudolf@digitaldocument.de wrote:

That is exactly what the code does.

To consider that I have done another test code:

m = 10000;
// Rz(360/m)^m
Rm = [for(  a = 360/m, i = 0, R = Rz(0);
i<=m;
R = R*Rz(a), i = i+1 )
if(i==m) R ][0];

P = circ();
Pm = PRm;
echo("Errors of P-PR^m");
echo(errEucl= matrixNorm(P-Pm));
echo(errMax = norMax(P-Pm));
echo(errSum = norSum(P-Pm));

// compute the rotations of P like in sweep.scad and keep the last one
PR_ = [for( a = 360/m, i = 0, Pi = P;
i<=m;
Pi = Pi*Rz(a), i = i+1 )
if(i==m) Pi][0];

echo("Errors of P-(((P*R)*R)*R...)");
echo(errEucl= matrixNorm(P-PR_));
echo(errMax = norMax(P-PR_));
echo(errSum = norSum(P-PR_));

echo(Id=Rm*transpose(Rm));

function Rz(a)= let(c=cos(a), s=sin(a)) [[c,s,0],[-s,c,0],[0,0,1]];
function matrixNorm(M) = norm([for(Mi=M) norm(Mi)]);
function norMax(p) = max([for(pi=p) norm(pi)]);
function norSum(p) = [for(pi=p) norm(pi)]*[for(pi=p) 1];

function circ(m=20) = [for(i=[0:m-1]) [cos(360i/m), sin(360i/m),0] ];

where P is a list of points. The code produced the following output:

ECHO: "Errors of P-P*R^m"

ECHO: errEucl = 1.84876e-12

ECHO: errMax = 4.1515e-13

ECHO: errSum = 8.26788e-12

ECHO: "Errors of P-(((P*R)*R)*R...)"

ECHO: errEucl = 1.84209e-12

ECHO: errMax = 4.16634e-13

ECHO: errSum = 8.23773e-12

ECHO: Id = [[1, 1.13101e-14, 0], [1.13101e-14, 1, 0], [0, 0, 1]]

which shows that there is no significant difference among the various ways
of computing errors as I was expecting. Note that Rm multiplied by its
transpose differs from the identity matrix by an error even smaller.

As I said before, if this is critical for some application, the use of
quaternions may decrease the accumulated errors.

As you could read in the last sentence of my previous post, I meanwhile
understand, that it is not the condition of the rotation matrix itself that
is problematic (orthogonal matrices have a perfect condition: 1),
but the condition of the function the determines the rotation matrix from a
given normal vector. A normal defining a rotation with an absolute angle
near 90° has itself a small angle against the plane defined by the polygon
and suffers under the tangens problem, respectively under a problematic
divison, spoiling the condition.
Your relative angle extrusion scheme therefore has indeed no serious
numerical problems exept for the first rotation, given that you don't get
into the pitfall of near 90° angles. In a sufficiently refined extrusion
path, this case will probably happen only with the first rotation, and the
good news is, that one cannot see it, because you build up upon the first
rotation in a relative way.
Therefore the significant error in your calculus will indeed be produced by
determinining the first rotation (and any relative rotation near 90°) and
can be opposed by splitting the rotation into a sequence of two rotations: a
90° rotation and a small 90°-α rotation.

As for the determination of the best z-rotation I could imagine that your
scheme could also be extented to optionally accept user input. Imagine a
user wants to extrude an airfoil or any other non rotation symmetric shape
in a specific twisted way.

To also allow for a parametrisation of the shape (e.g. a wing is a morphing
sequence of airfoils, therefore the shape generator needs at least one
argument) a last comment on

Ronaldo wrote

The analysis has so far shown that it is not the rotations that spoil the
result, but the determination of the rotation matrices. In order to allow
for morphing shapes along the path you can indeed calculate

Pi = PR  with R= R1R2*...*Ri

instead of P(i+1) = Pi*R(i+1)

But I guess you do this anyway.

BTW: if I had to extrude a shape along the line you showed, I would probably
exactly follow this way and implement a preprocessor that generates the
polygon sequence for extrusion with Naca_Sweep. :-)

--
Sent from: http://forum.openscad.org/

Parkinbot rudolf@digitaldocument.de wrote:

You probably meant 180°. And that requires a big account.

You are right saying we need to consider the condition of the process
during a sweep.  At each step in the global method (the one in the code of
sweep.scad in the List Comprehension Demos), it is computed a matrix (of a
minimum rotation) transforming the section from its base in the xy plane to
a plane orthogonal to the path tangent at a given point. At each step in
the local method (in my github repository), the minimum rotation is
computed to transform the previously transformed section to the next
tangent direction. Both method use the same function to compute the
rotation:

function rotate_from_to(a,b) =
let( axis = unit(cross(a,b)) )
axis*axis >= 0.99 ?
transpose([unit(b), axis, cross(axis, unit(b))]) * [unit(a), axis,
cross(axis, unit(a))] :
[[1, 0, 0], [0, 1, 0], [0, 0, 1]];

a minimum rotation that brings a to b.

When a and b have opposing directions (180°), there is an infinite number
of possible minimum rotation. That function solves this singularity by
arbitrarily outputing the identity matrix. This is a bug indeed, a subtle
one! In the context of a sweep, that arbitrary solution may lead to an
inversion of the section vertex circulation from CCW to CW which is the
source of wild twists reported in
http://forum.openscad.org/Twisty-problem-with-scad-utils-quot-sweep-quot-tc9775.html
. In that thread, may be found the original proposal of the local method.

When a path starts with a tangent direction [0,0,-1], all the sweep faces
will have a CW circulation and any boolean operation with the polyhedron
will be rejected by CGAL. So, yes, Parkinbot, you are right: there is a bug
in sweep. But not only at the start point: if the path has three or more
collinear points, the local method will show wild twists and face
circulation reversal!

To solve the bug I rewrote the minimum rotation function:

function rotate_from_to(a,b) =
let( axis = unit(cross(a,b)) )
axis*axis >= 0.99 ?
transpose([unit(b), axis, cross(axis, unit(b))]) * [unit(a), axis,
cross(axis, unit(a))] :
rotate_from_to(a,b+[1e-24,0,0]);

where I arbitrarily perturb the vector b to escape from the singularity.
That is enough to avoid the circulation changes and wild twists. I have not
check it but I guess the original global method, which is simpler, more
precise and faster then the local one, will produce basically the same
results then the local one with no wildness. Note that the rotation of the
first section of the sweep about to its own center is arbitrary but the
following sections will have a consistent rotation with the first one.

The non-diagonal elements of matrix rotate_from_to(a,b+[1e-24,0,0]) are
very small and may be reset to 0 without any drawbacks.

As for the determination of the best z-rotation I could imagine that your

It really accepts an initial rotation angle to cover that.

Yes, that is the way it is done.

In my version of sweep.scad, a function sweep_polyhedron() generates a
polyhedron data ([points,faces]) ready to be sent to polyhedron().  The
sweep() module just call that function and polyhedron() with its output.
The function sweep_polyhedron() may or may not add caps to the sweep ends
allowing additions of special caps or even polygons with holes (elsewhere
pre-processed). So, to do a sweep of a polygon with holes it is just a
matter of:

outer = sweep_polyhedron(shape=outer_section, path_transforms=pt,
caps=false);
hole  = sweep_polyhedron(shape=hole_section,  path_transforms=pt,
caps=false);
ppoly = partition_polygon_with_holes(outer,[hole]); // generate [points,
faces]
caps  = caps_of(ppoly, pt); // transform ppoly points by pt[0] and
pt[len(pt)-1], keep faces
lazzy_union([outer, hole, caps]);

when function partition_polygon_with_holes() be available. :)

simpler, more precise and faster then the local one, will produce basically
the same results then the local one with no wildness.

The problem with the global one is it produces a 360 degree twist for every
orbit of the Z axis. The relative one does not.

To avoid the singularity shouldn't you look at the next tangent and orient
it around the axis the way that indicates? How does a fixed direction
perturbation help?

On Wed, 30 Jan 2019 at 19:02, Ronaldo Persiano rcmpersiano@gmail.com
wrote:

Ronaldo wrote

you are right it is 180° and the critical operation is unit(). I haven't
studied the orginal sweep code extensively and it is a long while ago. But
it is obvious that the code can run into badly conditioned situations that
need special treatment. I will refer back to this discussion, and what I
have learned from it, when it is less academic for me, i.e. when I encounter
the practical problem having to extrude along a "strange" trajectory.

Because you mentioned it: My current implementation of the multi hole sweep
looks like this:

module sweep_multihole(outer, holes, convexity = 5)
difference()
{
sweep(outer, convexity = convexity);
sweeps(holes, convexity = convexity);
}

module sweeps(Dat, convexity = 5) // lazy union sweep
{
soup = sweeps(Dat);
polyhedron(soup[0], soup[1], convexity = convexity);
}

Yours looks straight forward to me.

--
Sent from: http://forum.openscad.org/