Yes the non-planar quads are not concave or convex in isolation but when
you add the diagonal the two triangle are. When they form a concave dimple
the shading alternates between light and dark and makes it visible that
they are dimples in a surface that should be all convex. When split the
correct way the shading changes monotonically, so it looks smoothly convex.
Not only is the shading no longer alternating but it also has smaller
changes between adjacent triangles, so you can get away with less of them.

Since the quads in a sweep are not in arbitrary places, they are arranged
in rings, I think it is probably sufficient to work out which way the quad
is twisting to pick the correct diagonal.

How do you get rid of the twist?

On Thu, 24 Jan 2019 at 00:01, Ronaldo Persiano rcmpersiano@gmail.com
wrote:

Yes the minimum rotation algorithm seems to create the twist, not sure why.
I wonder if a real spring is twisted.

Frenet-Serret  isn't twisted but it gets weird at the discontinuities in
curvature  near the ends.

[image: image.png]

I made a mistake earlier when I said the diagonals were equal. I was
subtracting the indices instead of the coordinates. Fixing that and picking
the shortest diagonal works for clockwise and anti-clockwise springs.
Thanks for pointing that out.

Why doesn't linear_extrude pick the shortest diagonal if that is the
solution?

On Thu, 24 Jan 2019 at 11:20, Parkinbot rudolf@digitaldocument.de wrote:

nophead wrote

Good question. You are right!

linear_extrude(50, twist = 360, slices=20)
translate([5,0]) square(10);

I have not observed it because it was never very important to me. To be
honest I don't know, whether linear_extrude() is implemented to produce
triags or quads. But you are right, the result would be way better if it
used always the shorter diagonal. Looks like an adversity that can easily be
fixed.

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

nop head nop.head@gmail.com wrote:

Just calculate the total section angular rotation about its "center" and
distribute this total angle along the path. It is an aproximation because
the twist rate may vary along the path but it works very fine with
helicoidal paths. You can accesses my version in:

https://github.com/RonaldoCMP/list-comprehension-demos/blob/master/sweep.scad
https://github.com/RonaldoCMP/list-comprehension-demos/blob/master/sweep.scad

The untwist is done by  construct_transform_path().

Coil springs mostly act in torsion - they twist and the spring expands and
contracts. They're usually made by wrapping wire or rod around a mandrel, so
they're not that twisted at rest.

Are you sure it's "the solution?"  Maybe it always works right for rotate
extrude, but I'm pretty certain that "longer diagonal" does not work as a
general answer.

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

No its the shorter diagonal, not the longer.  I am pretty sure it isn't a
general solution for any quad patch on a 3D surface. It might be for a
sweep where the quads are always in rings, but I don't know. Seems to work
with helical springs.

The spring machines I have seen feed wire though a nozzle against an anvil
that forces it to turn into a helix . The anvil moves allowing it to vary
the pitch and diameter continuously if necessary. Then it cuts it off and
starts a new one.

On Thu, 24 Jan 2019 at 17:47, NateTG nate-openscadforum@pedantic.org
wrote:

The answer kinda depends on whether the part of the surface in question
exhibits positive or negative curvature. On the outside of your spring,
which has positive curvature, it's clear which diagonal looks better, but
for the negative curvature inside, the answer is less clear. In some cases,
what looks "right" depends more on which direction you're viewing from than
which diagonal is shorter.

The best solution actually comes from staggering your points. I wrote code
for generating mathematical surfaces in the form of z = f(x, y), and I
recently rewrote it to sample points on an equilateral triangular grid
rather than a square grid to solve this problem. By choosing the shorter
diagonal across a quad, you're soft of approximating this solution.

There can still be spots that don't look perfectly smooth (usually because
the function has a high second derivative), but it generally handles these
things a lot better than the square grid. The next step in complexity would
be to write code to use a variable mesh density (or at least one that
remains constant relative to the face's orientation, rather than just
relative to the x/y plane), but without true variables, that would be very
difficult to implement. Even with true variables, I'm sure thesis papers in
topology have been written on the subject, some of which likely resulted in
some of the re-meshing algorithms in Meshlab.

On January 24, 2019 at 11:13:36, nop head (nop.head@gmail.com) wrote:

No its the shorter diagonal, not the longer.  I am pretty sure it isn't a
general solution for any quad patch on a 3D surface. It might be for a
sweep where the quads are always in rings, but I don't know. Seems to work
with helical springs.

The spring machines I have seen feed wire though a nozzle against an anvil
that forces it to turn into a helix . The anvil moves allowing it to vary
the pitch and diameter continuously if necessary. Then it cuts it off and
starts a new one.

On Thu, 24 Jan 2019 at 17:47, NateTG nate-openscadforum@pedantic.org
wrote:

OpenSCAD mailing list
Discuss@lists.openscad.org
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org

This is an ancient story.  The bottom line is that quads lead to all sorts of problems.

Best (if you can manage it) is to deal with purely triangular meshes.

Second best is to be exceedingly careful with non-planar quads.  Again - a good approach is
often to subdivide by producing a good point on your intended surface at the center of the
quad to produce FOUR triangles instead of trying to make do with only TWO.

Finally, when forced to choose a diagonal,  you can evaluate the error at the center of the
quad - or you can blindly choose the smaller of the two diagonals.

--
Kenneth Sloan
KennethRSloan@gmail.com
Vision is the art of seeing what is invisible to others.

On Wed, Jan 23, 2019 at 10:05:58PM +0000, nop head wrote:

The current rendering uses a single shade that depends on the normal
of the triangle.

If the length of the diagonals are not equal, chosing the shorter
triangle reduces the "aspect ratio" (longest possible line segment
divided by the largest perpendicular segment). As an example: imagine
two equilateral triangles joined to a rhombus.

Split them into the two original triangles and the aspect ratio is
almost one (sqrt(3)/2 iirc). Split it the other way and you get a much
larger aspect ratio. Something like sqrt(3) or even more.

The aspect ratio being small means that the points of that triangle
will be colored the same to "close" points. A pointy triangle has
points far away that will be colored the same due to having the same
normal.

For presentation (i.e. not 3D printing) you can do "phong
shading". This has something to do with averaging colors close to an
edge between two triangles. An object shaded that way will look
"perfectly round" until you start examining the image at the pixel
scale. This would be great for quick rendering realistic
representation of an object where you're working with a low $fn during
development but are going to increase $fn for final rendering to be
printed.

That sounds like a plan.

The inside of the spring is concave (in one direction) (I'm guessing that
this is exaclty the "partially convex" that you say :-) ).

Roger. 

--
** R.E.Wolff@BitWizard.nl ** https://www.BitWizard.nl/ ** +31-15-2049110 **
**    Delftechpark 11 2628 XJ  Delft, The Netherlands.  KVK: 27239233    **
The plan was simple, like my brother-in-law Phil. But unlike
Phil, this plan just might work.

I didn't get that in the sense of the general case with n slices. Can you
show linear_extrusion code that has this problem?

Obviously, when you do a twisted linear_extrusion say with 5 slices and the
result is exactly what you want, you are better of.

Whosawhatsis wrote

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

No what I mean't was no shape can be totally concave on the outside, it
must curve back to meet itself and be convex as it does so.

Even on the inside of the spring, a quad oriented along it should still be
folded to be convex. Locally the surface is more convex around the ring and
less concave along it. So perhaps you would look the four neighbouring
quads on its sides and pick the most curved direction to decide it
should be convex or concave.

On Fri, 25 Jan 2019 at 12:18, Rogier Wolff R.E.Wolff@bitwizard.nl wrote:

Kenneth Sloan wrote

This is not practicable for linear_extrude in a general sense, since it uses
an integral slices variable. And even then it is not trivial as it would
imply to also subdivide the polygon, i.e. getting 8 vertices for a square.
So multiplying slices by 2 will multiply the number of triags by four.

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

Actually linear extrude seems to chose the longest diagonal. It definitely
looks bad for negative twists but less so for positive ones. I wonder what
it would look like if it choose the shortest on positive twists.

On Fri, 25 Jan 2019 at 13:43, Parkinbot rudolf@digitaldocument.de wrote:

I think the asymmetry is due to the lighting direction. Even with positive
twist it looks to have concave dimples.  My guess is it would look better
in both case if it chose the shortest diagonal to make the quads convex.

On Fri, 25 Jan 2019 at 14:01, nop head nop.head@gmail.com wrote:

triangulation is scaled non-isotropically.
Can you explain that please.

With linear_extrude it doesn't matter how many slices you have. It is the
twist direction that determines the shortest diagonal. I.e. with short fat
quads or tall thin quads.

The reason we are start with quads is that we are extruding or sweeping a
polygon, so each slice has the same vertices. Joining the vertices that
match on each layer gives quads and then we have to pick the correct
diagonal to give triangles.

On Fri, 25 Jan 2019 at 15:24, Kenneth Sloan kennethrsloan@gmail.com wrote:

nophead wrote

I don't see much difference, to be honest, left is my simulation of
linear_extrude with sweep:

http://forum.openscad.org/file/t887/rectquad.png

I also changed my springs code to get more lengthy quads.

http://forum.openscad.org/file/t887/springquad.png

Of course, the less crooked the rects are, the less negative effects you
get. So better design your code to produce almost rect quads.
According to Kenneth' proposal, linear_extrude could subdivide a polygon
before applying a twist. I also tried that: If you feed the polygon of a
subdivided square in linear_extrude(), OpenSCAD will discard any vertices
inserted by subdivision before extrusion. So no obviously change to the
right side. So I used sweep to show how this looks.

http://forum.openscad.org/file/t887/rectquad1.png

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

I have been doing all my tests on cylinders.

Yes the square is terrible without subdivision. I do that myself by
generating a square with extra vertices along the sides, slightly
perturbed, so they don't get removed.

I think the polygon should be subdivided depending on height / slices, so
that the quads are roughly square.

On Fri, 25 Jan 2019 at 17:35, Parkinbot rudolf@digitaldocument.de wrote:

I have been experimenting with sweeping and found a curious fact: If I use
the absolute minimum rotation angle method then a sweep that orbits the Z
axis gets twisted. If I use the relative method then it doesn't, my spiral
springs are no longer twisted around the wire's axis.

[image: image.png]

By absolute I mean each slice is rotated from flat to align it with the
tangent. This is the result. It goes wrong when the tangent points straight
down as there is no minimum rotation from up to down, they are all the
same.

[image: image.png]

By relative I mean calculate the rotation from each tangent to the next one
and accumulate those rotations by multiplication. That of course has the
advantage that it doesn't go wrong when the path points straight down as a
smooth path never fully doubles back on itself.

A consequence is that closed loops don't always meet at the same rotation,
see the knot below. It is less twisty but it needs to twist 180 to meet
properly.

[image: image.png]

Looking at Ronaldo's solution he ensures it twists the right amount to meet
properly when the path is closed by distributing a calculated axial
rotation.

I was surprised that the absolute rotation gives a different result to the
accumulated relative rotations. I expected them to be mathematically
equivalent, but obviously they are not. I think it is because the absolute
one sort of takes a short cut when it does the minimum rotation all in one
go so arrives at a different orientation around the tangent than it does if
it follows the path of rotations. On the other hand the absolute method
always lines up where it meets because the absolute rotations for the start
and the end must give the same twist as they are the same rotation.

One thing to note is that it is possible to accumulate the rotations in a C
style for loop now. It no longer requires a recursive function.

function skin_points(profile, path, loop) =
let(len = len(path),
n = len - 1,
tangents = [for(i = [0 : n], b = loop ? (i + n) % len : max(i - 1,
0), a = loop ? (i + 1) % len : min(i + 1, n)) path[a] - path[b]],
rotations = [for(i = 0, m = rotate_from_to([0,0,1], tangents[0]); i
< len; m = rotate_from_to(tangents[i], tangents[(i + 1) % len]) * m, i = i

• 1. m])
     [for(i = [0 : n]) let(t = orientate_r(path[i], rotations[i]))
     for(p = profile)
     transform(p, t)
     ];

I am probably reinventing the wheel but I get to understand it a lot more.

On Thu, 24 Jan 2019 at 16:26, Ronaldo Persiano rcmpersiano@gmail.com
wrote:

nop head nop.head@gmail.com wrote:

Any rotation is a minimum rotation for a pair of directions in a specific
set: those in the plane orthogonal to the rotation axis. So, if you compose
a minimum rotation from a to b with the minimum rotation from b to c, the
resulting rotation will be a minimum rotation from a to c if and only if
vectors a, b and c are coplanar, that is they all have the same rotation
axis. I haven't test it but I would expect that the global (absolute) and
local (relative) methods would have the same result if all path tangent
vectors are in a plane (and no twist adjustment is needed for closed
paths). Otherwise, they would differ.

Yes I appreciate that there are an infinite number of ways to specify the
rotation. I am aiming for the simplest code to draw springs and wires, etc,
so I want to specify just the path and let the code work out the rotations.
Speed is more important to me than numerical accuracy.

I was using a mixture of the absolute minimum angle method and
Frenet-Serret for different sections of this ribbon cable depending on
which broke. I.e. up and down are a problem for the former and
discontinuous curvature for the latter.

By moving to the relative mode it all seems stable and doesn't have wild
twists when the curve changes direction, so I can ditch Frenet-Serret
unless I want to design a roller coaster!

[image: image.png]

So thanks for all your help.

On Tue, 29 Jan 2019 at 18:02, Ronaldo Persiano rcmpersiano@gmail.com
wrote:

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/