In that case, I must have applied the Bezier in a wrong way. These curves are indeed built with Bezier approach, but I applied it pairwise (each pair generates 4 CPs), so points between each pair are Bezier points.
$ Runsun Pan, PhD
$ libs: scadx, doctest, faces(git), offline doc(git), runscad.py(2,git), editor of choice: CudaText ( OpenSCAD lexer); $ Tips; $ Snippets 
20161122 16:39 GMT02:00 runsun <[hidden email]>: Ok, I see what you have done, you did right. It is a C1 interpolation spline, sometimes called Hermite spline. You possibly estimates the tangents at the points and find the Bezier CPs of each polynomial arc based on the tangents and the extremes points of the arc. But then, the BezierCPs are not all interpolated and two of them for each arc control the starting end ending arc tangent. The same happens with Bezier CPs of surfaces: some CPs are interpolated and the others control its partial derivatives.In that case, I must have applied the Bezier in a wrong way. These curves _______________________________________________ OpenSCAD mailing list [hidden email] http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org 
In reply to this post by runsun
runsun,
Perhaps, the following code may help you to understand the relationship between the CPs and the Bezier triangle surface. You don't need to understand all the code. Run it, change the first two or three parameters and examine the results.

@ Ronaldo, Wow, this is cool. Thx man.
$ Runsun Pan, PhD
$ libs: scadx, doctest, faces(git), offline doc(git), runscad.py(2,git), editor of choice: CudaText ( OpenSCAD lexer); $ Tips; $ Snippets 
Looking at the result shapes, I now understand how the Bezier Triangle Patch manages to go from CPs to poly.
I guess it's a method to "generate a new surface", which is different from "smoothing existing surface" that I originally imagined.
$ Runsun Pan, PhD
$ libs: scadx, doctest, faces(git), offline doc(git), runscad.py(2,git), editor of choice: CudaText ( OpenSCAD lexer); $ Tips; $ Snippets 
In some sense, Bezier curves and surfaces are a smoothing process of the control points. Look at the following sequence of operations behind a Bezier quadratic: From left to right, the first figure shows the CPs (red marks) of a quadratic arc: it has one corner point. In the second, we cut the corner point out through the median points of the two edges and get 2 new green corners. In the next one, we cut the two new corners by the same process and get 4 new green corners. If we iterate indefinitely this process of cutting corners, we will get a smooth curve (the blue one): a parabola! Note that the median point of each new edge resulting from the corner cut (the blue marks) is a point of the final parabola. The corner cutting method is a smoothing process. It is not what you originally imagined because it is not an interpolation process. However, interpolation methods sometimes produce unexpected oscillations, that is a waving behavior we don't see in the data points. Even natural spline interpolation method is unstable. On the other hand, Bezier and Bspline methods are not interpolations of their CPs but have no more "oscillations" than their CPs. In CADGD, the Bezier and Bspline representations are invaluable when we understand well the relationship between the shapes and their CPs and realize why CPs are called as such. 
In reply to this post by Ronaldo
Thanks Ronaldo,
nice example. To cherish stuff like that an image of what the code produces is invaluable: 
In reply to this post by Ronaldo
I have found a family of methods appropriated to build organic shapes. They are called subdivision algorithms and are founded in the subdivision methods usual to Bezier and Bspline curve and surfaces. They may be regarded as a generalization of the corner cutting method I illustrated in my previous message.
Imagine you start with a cube. As a first step, we cut out each vertex of the cube by a plane through its incident edges. Twenty four new vertices now substitute the initial twelve and eighth new facets are added. Now, apply the same procedure (vertex cuts) to the resulting polyhedron again and again. The polyhedron facet number will grow exponentially and the polyhedron surface may converge to a smooth surface depending on the chosen cuts. The methods I have found are different from that in the way they cut the vertices but that is the general idea behind them. The methods are very simple and could be well implemented in OpenSCAD if we had more powerful means to build data structures. The iteration steps of the method require that the polyhedron definition be refined. I have started to write a triangulation data structure in OpenSCAD. But it is very hard to do it with just lists and even harder to debug it. Here are some references of papers on the subdivision methods: Interpolating Nets Of Curves By Smooth Subdivision Surface Interpolatory sqrt(3)subdivision sqrt(3)subdivision The first reference have nice figures showing how promising is the method for organic modelling. 
@Ronaldo, Thx for the precious references.
Btw, I saw you mentioned data structure in several posts. What do you have in mind ?
$ Runsun Pan, PhD
$ libs: scadx, doctest, faces(git), offline doc(git), runscad.py(2,git), editor of choice: CudaText ( OpenSCAD lexer); $ Tips; $ Snippets 
In reply to this post by Ronaldo
@ronaldo  don't forget the doosabin resmooth operator which is very useful.
It can be easily seen in Wings3D under the Bodies menu. Where you can also compare it to the regular smooth algorithm. 
20161128 20:00 GMT02:00 Neon22 <[hidden email]>: @ronaldo Â don't forget the doosabin resmooth operator which is very Thank you, Neon22, for the reference. I have found references to DooSabin smoothing method but have not crossed with their paper yet. There are lots of method in this subdivision arena. _______________________________________________ OpenSCAD mailing list [hidden email] http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org 
In reply to this post by runsun
Yes, you are right. I already implemented triangulations in Pascal, C, and Ruby. I have no previous experience with functional languages, never studied Haskell for instance. So, I have no idea how data structures are implemented in a functional language environment. I have found an implementation of quadedge in Haskell. It seems simple but it is a mystery for me. And I have found a paper on persistent triangulation data structure and an algorithm for convex hull based on it. So, there are alternatives. I have no idea how OpenSCAD could help me if something besides lists were available. Parkingbot has already mentioned how hard and lengthy is to have to copy a list when you need to change anything in it even a single element. I am doing my best to write a triangulation data structure (based on quadedge) without concepts like objects, records and pointers. I use integers as references instead of pointers and arrays instead of records. It is a nonmutable FORTRAN code. I don't know what could more expressive but I feel myself constrained. That is all. Yes, I know: OpenSCAD is not a general purpose programming language. 
@Ronaldo, what I meant was, what kind of data structure you want ? Can you describe it ?
$ Runsun Pan, PhD
$ libs: scadx, doctest, faces(git), offline doc(git), runscad.py(2,git), editor of choice: CudaText ( OpenSCAD lexer); $ Tips; $ Snippets 
20161128 23:03 GMT02:00 runsun <[hidden email]>: @Ronaldo, what I meant was, what kind of data structure you want ? Can you Any flexible data structure able to store, modify and travel a triangulation or perhaps a more general subdivision of a manifold. I know two data structure used in Brep and computational geometry: wingededge and quadedge. I have already worked with the later. See also QuadEdge Data Structure and Library. The original paper where it was described ("Primitives for the manipulation of general subdivisions and the computation of Voronoi diagrams") is not freely available in the net. But I should warn you it is a lot abstract. _______________________________________________ OpenSCAD mailing list [hidden email] http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org 
@Ronaldo, Many thanks. Like all previous references you cited, these are very helpful to me. I'll keep them for future readings.
In the mean time, although not familiar with geometry studies in general, I spent quite some time working on how I'd like to use OpenSCAD in geometry, which I believe is more or less on the same path of yours. I'd like to share some of my idea here: 1) hash: This is the basic data structure I use in almost every func/mod: data = ["a",3,"b",4] hash( data, "a" ) ==> 3 data2 = update( data, ["b",5, "c",1] ==> ["a",3,"b",5,"c",1] The core of this hash is extremely simple: function und(x,df)= x==undef?df:x; function hash(h,k, notfound)= ( und( [ for(i=[0:2:len(h)2]) if( h[i]==k ) h[i+1] ][0], notfound ) ); It can be set to include actioninstructions. One example is posted in this thread previously. Another one, with a function growPts, pts = growPts( seed, ["x", 3, "y", 4, "z", 5, "transl", [2,3,4] ] ) we can generate a poly from a seed by, literally, "drawing with the code". It's extremely handy when creating a poly. 2) typechecking: this is essential for me to work on geometry, especially checking geometry type (other than the language typechecking like int, str, array, etc) :  0> gtype(P)= "pt"  1> gtype([P,Q])= "ln"  2> gtype([P,Q,R])= "pl"  3> gtype("a")= undef  4> gtype([P,Q,R,P])= undef  5> gtype([[P,Q],1])= "cir"  6> gtype([2,[P,Q]])= "cir"  7> gtype([P,1])= "ball"  8> gtype([0.5,P])= "ball" 3) Other geometry inspections: isSamePt(P,Q) isonline (pt,line ) // is pt on the line ispts(o,dim=3) // are all items in o a point is90( pq1,pq2 ) // are 2 lines 90deg to each other isOnPlane ( pqr, pt ) // is pt on plane pqr isSameSide ( pts,pqr ) // are all pts on the same side of pqr iscoline ( a,b ) // Chk if pts in a and pts in b are all on the same line. a,b could be either pt or pts isparal ( a,b ) // are a and b parallel. a, b could be line or plane ispl ( x,dim=3 ) // is x a plane ? ispt ( x,dim=3 ) // is x a point ? 4) examples of some core geometry functions: angle( pqr ) // return angle of pqr angleBisectPt( pqr ) // Return a pt that the line ptQ cut angle(pqr) in half. anglePt( pqr, a=undef, a2=0, len=undef ...) // Create a pt at certain angle/angles in respect to pqr * left: create pts on 30 and 60 degree on the pqr plane * middle: create a circle on the pqr plane (varying a with a2=0) * right: create a circle 90degree to pqr plane (a=90 with varying a2) * The combination of a and a2 allows creation of any pt in respect to pqr: dist(a,b) // distance between a, b, where a,b could be pt, line, plane. intsec( a,b ) // Return the intersection between a,b where a,b could be line or plane movePts(pts,to,from) // Move a poly from any location to any location * Fig below: Move a poly (yellow) on the [P,Q,R] plane to a [S,T,U] plane (red shape): onlinePt( pq, len, ratio) // Return a pt that is on the line pq with length=len or ratio of pq projPt(a,b,along) // Return the point where a projects on b or b projects on a. where a, b could be: pt, line or plane projPl(pl1,pl2,along,...) // Return projection of pts on plane pl1 onto pl2 along the line "along" These are just part of my lib that makes my journey of OpenSCAD full of joy. As obvious from the above functions, a "[P,Q,R]based" referencing/coordination system is what I use to drive through the space. I've shown the possibility of moving any poly in any pqr location to any other with the movePts. Although I don't know what exactly the move is formally called, I believe it's some sort of coordination transition. The next step is to allow the moved poly be either squeezed or expanded based on the pqr angle (so the destination is no longer an orthogonal coordinate). That might be one step closer to the organic representation. We'll see.
$ Runsun Pan, PhD
$ libs: scadx, doctest, faces(git), offline doc(git), runscad.py(2,git), editor of choice: CudaText ( OpenSCAD lexer); $ Tips; $ Snippets 
Nice approach, runsun. I have being following your work and solutions. But the data structure I need now are quite more complex.
In your code faces.scad, for instance, you build the face list of polyhedron, that is, the topological structure of the polyhedron. You implicitly know the topology, a fixed one (a cube, a rod, a chain), and, for each one of them, you write the polyhedron face lists. What I need to implement the subdivision methods I referred before is something more general. I need a data structure that is able to store any polyhedron but also it should accept local changes and to be traversed in different ways, all that in a efficient manner. The OpenSCAD polyhedron data is not satisfactory. Suppose I need to swap an edge of a triangulation (that is, switch the diagonal of the quadrilateral adjacent to an edge). I will have to find (by searching all the face list) the two faces that are adjacent to the edge. The subdivision refinement needs to do hundreds or thousands swaps besides vertex insertions. That is not an efficient data structure. In the quadedge structure, the atom element is a directed edge (and its symmetrical  its inverted direction counterpart). Each directed edge has two adjacent vertices and two adjacent faces. From each edge, it is possible to get the following edge with the same left face, the following CCW edge with the same origin vertex or the symmetrical edge. All those access are done in constant time. And edge swap is also done in constant time. Well, not always... In a C, Pascal or Ruby, the atom element  the directed edge  would be an object or a record and the connection between edges is done by pointers. So the swap is done by changing a few (fixed) pointers. But, in OpenSCAD, each time I do a swap I need to rewrite the whole structure. Therefore, a swap  and any changing operations  is not a constant time operation anymore. And, as a consequence, the subdivision method may be O(n3) (where n is the number of edges in the last refinement). May be worst than that. I have not enough experience with OpenSCAD to say what would happens in a step of a subdivision method when we reach thousands of vertices. 
I don't know much about the quadedge structure and whether its an optimisatoin for quad shaped meshes but I am very familiar with the winged edge and it guarantees topological correctness for N sided polygons.
(saved format consists of vert lists, edge lists and face lists. Structure is an edge with a left and right faces and CCW and CW edges leading off each of the two edge vertices) The only problem in using these (better IMHO) structures is when loading in a poorly connected error filled nonmanifold polygonal object. You can't get a bad model properly stored in the winged edge structure and the tools to autofix it have not risen out of the morass yet. So you need to have a way to say if an object is not manifold and the tools need to work on that structure also. 
You are right, Neon22. Both data structures were conceived so that only manifolds may be built and stored. If we try to represent a cube with a missing face, for instance, we will end up with a manifold that is not what was intended. And may be very hard to identify that.
Quadedges allows to represent general manifold subdivisions. That means the faces may be polygons of n sides, not only quads. The name comes from the four incarnations of each edge: two opposed directions of the edge and their dual. If we don't need the dual graph, it is possible to represent just two incarnations of each edge. I am trying now is to explore implicit representation of the triangulation topology for certain kinds of triangulation that may be used to model furcation junctions. This way I won't need an explicit topology data structure to apply Kobbelt subdivision, for instance. And the quadedge implementation has been postponed. The advantage of Kobbelt subdivision over DooSabin, for instance, is that it works for any triangulation and its subdivisions are rather regular triangulation. 
In reply to this post by Ronaldo
I c. It seems that we both have the same "goal of organic representation" in mind, but approach it differently. A good analogy of this difference might be that I started it from the bone (to serve as the core for adding skin later), but you skip the tedious details and tackle the skin directly. I didn't know it is possible to handle it your way until after seeing your code and examples, which is a pleasant surprise. Reading your description about a quadedge data structure over and over, I am thinking that I might be able to come up with a structure for that, by using an "edgebased" list/hash instead of all the "pointbased" lists that I used so far. But, like you say, I can imagine it encounters the speed issue you mentioned quickly.
$ Runsun Pan, PhD
$ libs: scadx, doctest, faces(git), offline doc(git), runscad.py(2,git), editor of choice: CudaText ( OpenSCAD lexer); $ Tips; $ Snippets 
After a long struggle computing indices, I have been able to implement Loop's subdivision algorithm for special cases of triangulation. This is my first result, a very rough one.
To avoid the implementation and test a general triangulation data structure, I decided to approach this first test with implicit triangulation of a class of triangulations that would be able to model furcations. The green grid in the image shows the specific triangulation used with this model. The triangulation is made of a circular set of Tpatches triangulations (angular sectors). In this model the triangulation has 12 sectors with 9 triangles each. Above the green triangulation you will find the control point grid for the hand upper surface. It has a total of 6*12 + 1 = 73 control points. The lower surface is a reflection of the upper one. The hand was generated by a modified version of Loop's method where I managed to interpolate prescribed border Bezier curves of degree 3: the openings of the furcation and the interstices between "fingers". They total 12 curves defined by the border control points. The surface of the image was generated by 2 iterations of Loop's method and has 16 subtriangles for each triangle in the triangulation (a total of 1728 triangles). Each iteration of the method multiplies the triangle count by 4: so it grows exponentially. Two iterations required 2 sec. The model has some artifacts. Although it is C2 except in a few points and has the ability to interpolate border curves, I see no way to make the subdivision method to interpolate tangent planes at the border. So the full model with the two surfaces is not even C1 at the joint lines. I have found methods that are able to interpolate even curves along the surface but not tangent planes. Therefore to get a really C1 (C2 almost everywhere) furcation is necessary to build one triangulation for the entire furcation. 
Free forum by Nabble  Edit this page 