Bezier Courves - elevation model

classic Classic list List threaded Threaded
17 messages Options
Reply | Threaded
Open this post in threaded view
|

Bezier Courves - elevation model

rob63
This post has NOT been accepted by the mailing list yet.
Hi,

i'm new to OpenSCAD and i hav a problem with no idea how to solve it.

I have a raster of 150 points x and y coordinates are each 25 so it is a map of squares.
For each point on this map i have different z measures.
Now i want to create something like this
http://www.ozflux.org.au/monitoringsites/capetribulation/capetrib_dem.html

My idea was to project a bezier through this points, in x and z axis horizontal, vertical and with 45 Degrees.

How can i do this?
A small examlpe would be helpful.

Regards
Rob
Reply | Threaded
Open this post in threaded view
|

Re: Bezier Courves - elevation model

droftarts
Would 'Surface' do the trick? https://en.wikibooks.org/wiki/OpenSCAD_User_Manual/Other_Language_Features#Surface
It's not smooth like a Bezier curve, though. I think that's probably beyond OpenSCAD. You could smooth it more by sub-sampling between the points in the heightmap.

Ian
Reply | Threaded
Open this post in threaded view
|

Re: Bezier Courves - elevation model

Ronaldo
In reply to this post by rob63
This a classical problem of interpolation of points in a grid. There is many solutions for it. If you don't need smoothness of the surface, the simplest is to apply a bilinear interpolation. The surface will have fold lines at the grid edges. I guess that is the method surface() applies.

To get some smoothness (C1), the natural way is to estimate the partial derivatives at the grid points and adjust a cubic polynomial between two adjacent points in the grid (Bicubic interpolation). A evaluations of the cubics at intermediate points are needed.

To get a second order smoothness (C2) the usual way is to choose spline interpolation. This last method is implemented in splines.scad of thingi:1208001. Study the included example and come back if you have doubts.
Reply | Threaded
Open this post in threaded view
|

Re: Bezier Courves - elevation model

Parkinbot
Ronaldo,
thanks pointing that out. The scheme can easily be used to interpolate surfaces. For this nNspline() must be called twice. The second call operates over the transposed result of the first call.  
My quickly composed example code also compares two possible output schemes: cubes and polyhedron



use <splines.scad> // http://www.thingiverse.com/thing:1208001

A = [  // coarse data of surface
  [1,   1, 0.1, 1,  2],
  [1,   2, 4,   2,  1],
  [0.1, 4, 2,   1,  3],
  [1,   1, 1,   1,  4],
];


spline2D();

module spline2D()
{
  B = nSpline(A, 100);            // interpolate X
  C = nSpline(transpose(B), 100); // interpolate Y
 
  // using cubes  and very(!!) slow F6 rendering
  translate([-5, 0, 0])  
    surf(A);                // plot coarse data
  translate([1, 0, 0])
    surf(C, len(A)/len(B)); // plot interpolated data

  // using triangulation with polyhedron fast F6 rendering
  translate([-4.5, -5, 0])  
    surfaceplot(A);                // plot coarse data
  translate([1, -5, 0])
    surfaceplot(C, len(A)/len(B)); // plot interpolated data
}


function transpose(A) =  [for(i=[0:len(A[0])-1]) [for(j=[0:len(A)-1]) A[j][i]]];

module surf(A, d=1) // test
   for(i=[0:len(A)-1])
      for(j=[0:len(A[0])-1])
        translate([d*i, d*j]) cube([d,d,A[i][j]], center = false);

module surfaceplot(A, d=1)
{
  map = [for(i=[0:len(A)-1], j=[0:len(A[0])-1]) [i*d, j*d, A[i][j]]];
  base = [for(i=[0:len(A)-1], j=[0:len(A[0])-1]) [i*d, j*d, 0]];  
 
  faces = concat(faces(A), faces(A, len(A)*len(A[0])), sidefaces(A));
  polyhedron(concat(map, base), faces);
 
  // used for top and bottom faces  
  function faces(A, offs=0) =   let(X = len(A), Y= len(A[0]))
    [for(i=[0:X-2], j=[0:Y-2], k=[1,0]) let(m=i*Y+j+offs)
      offs? k?[m+1, m, m+Y]:[m+1, m+Y, m+Y+1] : k?[m, m+1, m+Y]:[m+1, m+Y+1, m+Y]];

  function sidefaces(A) =
    let(X = len(A), Y= len(A[0]))
    let(A=[for(i=[0:Y-2], k=[1,0]) k?[i+1, i, i+X*Y]: [i+1, i+X*Y, i+1+X*Y]])
    let(B=[for(i=[0:Y-2], k=[1,0]) let(m=i+(X-1)*Y) k?[m, m+1, m+X*Y]: [m+1, m+1+X*Y, m+X*Y]])
    let(C=[for(i=[0:X-2], k=[1,0]) let(m=i*Y) k?[m, m+Y, m+X*Y]: [m+Y, m+Y+X*Y, m+X*Y]])
    let(D=[for(i=[0:X-2], k=[1,0]) let(m=(i+1)*Y-1) k?[m, m+X*Y, m+Y]: [m+Y, m+X*Y, m+Y+X*Y]])
    concat(A, B, C, D);  
     
}
Reply | Threaded
Open this post in threaded view
|

Re: Bezier Courves - elevation model

Ronaldo
Parkinbot,

That was basically the solution I had devised. However, there is two small slips in your code: the spline map has a reflection in relation to the data plot and its scale in x and y are different from the data plot (both noticeable in your image). As the Nspline() just compute the z values, the transposition must be applied twice in order to agree with the x,y assignments in surfaceplot(). The lack of a domain data is responsible for the difference in scales. The following code correct these slips. I include here just the changed modules where I marked with //<<< the changed lines. There is a missing scale/domain in module surf() but I ignored it and dropped it from the code.

module spline2D()
 {
   B = nSpline(transpose(A), 100); // interpolate X   //<<<
   C = nSpline(transpose(B), 100); // interpolate Y

  // using cubes  and very(!!) slow F6 rendering
  translate([-5, 0, 0])  
    surf(A);                // plot coarse data
  translate([1, 0, 0])
    surf(C); // plot interpolated data

   // using triangulation with polyhedron fast F6 rendering
   translate([-4.5, -5, 0])
     surfaceplot(A);                // plot coarse data
   translate([1, -5, 0])
     surfaceplot(C); // plot interpolated data   //<<<
 }

 module surfaceplot(A, d=[0,5,0,5])    //<<<
// d =[ xmin, xmax, ymin, ymax ]
 {
   map = [for(i=[0:len(A)-1], j=[0:len(A[0])-1]) [d[0]+i*(d[1]-d[0])/(len(A)-1), d[2]+j*(d[3]-d[2])/(len(A[0])-1), A[i][j]]];    //<<<
   base = [for(i=[0:len(A)-1], j=[0:len(A[0])-1]) [d[0]+i*(d[1]-d[0])/(len(A)-1), d[2]+j*(d[3]-d[2])/(len(A[0])-1), 0]];   //<<<

   faces = concat(faces(A), faces(A, len(A)*len(A[0])), sidefaces(A));
   polyhedron(concat(map, base), faces);

   // used for top and bottom faces
   function faces(A, offs=0) =   let(X = len(A), Y= len(A[0]))
     [for(i=[0:X-2], j=[0:Y-2], k=[1,0]) let(m=i*Y+j+offs)
       offs? k?[m+1, m, m+Y]:[m+1, m+Y, m+Y+1] : k?[m, m+1, m+Y]:[m+1,
 m+Y+1, m+Y]];

   function sidefaces(A) =
     let(X = len(A), Y= len(A[0]))
     let(A=[for(i=[0:Y-2], k=[1,0]) k?[i+1, i, i+X*Y]: [i+1, i+X*Y,
 i+1+X*Y]])
     let(B=[for(i=[0:Y-2], k=[1,0]) let(m=i+(X-1)*Y) k?[m, m+1, m+X*Y]:
 [m+1, m+1+X*Y, m+X*Y]])
     let(C=[for(i=[0:X-2], k=[1,0]) let(m=i*Y) k?[m, m+Y, m+X*Y]: [m+Y,
 m+Y+X*Y, m+X*Y]])
     let(D=[for(i=[0:X-2], k=[1,0]) let(m=(i+1)*Y-1) k?[m, m+X*Y, m+Y]:
 [m+Y, m+X*Y, m+Y+X*Y]])
     concat(A, B, C, D);

 }
Reply | Threaded
Open this post in threaded view
|

Re: Bezier Courves - elevation model

Parkinbot
Thanks Ronaldo,

well, matrix A was already defined in a transposed fashion ;-) But you are right with the scaling. I more or less sketched the code quickly and without much care, just to point out some how-to-use code example to the thread starter. surfaceplot() should also calcuate the minimum of A and use a ground level below this value to guarantee for a welldefined polyhedron. But this aspects are already beyond the scope example code.

Rob, it would be nice if you could post the data you want to apply the scheme for. Because when playing around, I also found another, more severe data related problem with nSpline(), e.g. when trying to interpolate a series of points with constant values, a matrix full of NaNs results.

Whenever Ronaldo and me will attack our common interpolation library it will be time to revise this (terrible) piece of code.
Reply | Threaded
Open this post in threaded view
|

Re: Bezier Courves - elevation model

Ronaldo
If A is already transposed and the data plot used as such, it will represent its surface reflected. So, you need to transpose A before send it to the plot data module. Look your image to see what I mean. 

I will take a look later in the Nan issue.

Em 7 de mar de 2017 14:05, "Parkinbot" <[hidden email]> escreveu:
Thanks Ronaldo,

well, matrix A was already defined in a transposed fashion ;-) 

_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org
Reply | Threaded
Open this post in threaded view
|

Re: Bezier Courves - elevation model

Parkinbot
Ronaldo wrote
If A is already transposed and the data plot used as such,
Certainly you are right.
Reply | Threaded
Open this post in threaded view
|

Re: Bezier Courves - elevation model

Ronaldo
In reply to this post by Parkinbot
Parkinbot wrote
Because when playing around, I also found another, more severe data related problem with nSpline(), e.g. when trying to interpolate a series of points with constant values, a matrix full of NaNs results.
Rudolf,

I have found the origin of the NaNs. In your spline code, you make the spline knot intervals proportional to the input data segment length. When we have two subsequent equal data points, this will define two subsequent equal knots. However, the distances between subsequent knots are divisors of the right-hand terms of the equations system to solve.

A simple workaround is to change lineT1(), that is part of linear_integral() computation, avoiding 0 length knot intervals:

function lineT1(S, i, k) =
// recursive helper of lineT2()
    let(t1 = pow((S[i][k]-S[i-1][k]),2))
    k==0 ? t1 : t1 + max(lineT1(S, i, k-1),1e-12);
 
BTW, I don't think you need to recursively find the arc length sequence at all. To build the linear system you only need the knot interval lengths and this is just the input data segment lengths.

In my version of general purpose natural splines, I use uniform knot spacing. My comparison tests show not much difference between uniform and non-uniform spacing. Anyway, if two subsequent data points are equal, the repetition might/should be eliminated as a non-sense.

Reply | Threaded
Open this post in threaded view
|

Re: Bezier Courves - elevation model

jon_bondy
In reply to this post by Ronaldo
I replaced the original spline2D() and surfaceplot() routines with the
ones below, and went from 4 small squares to 3 small squares and one
huge one.

Whenever the dust settles, if one of you could provide a final version,
I would appreciate it.  I'm keeping copies of things like this for
reference should I need it.

Thank you for all that you do for us!

Jon


On 3/7/2017 11:06 AM, Ronaldo wrote:

> Parkinbot,
>
> That was basically the solution I had devised. However, there is two small
> slips in your code: the spline map has a reflection in relation to the data
> plot and its scale in x and y are different from the data plot (both
> noticeable in your image). As the Nspline() just compute the z values, the
> transposition must be applied twice in order to agree with the x,y
> assignments in surfaceplot(). The lack of a domain data is responsible for
> the difference in scales. The following code correct these slips. I include
> here just the changed modules where I marked with //<<< the changed lines.
> There is a missing scale/domain in module surf() but I ignored it and
> dropped it from the code.
>
>
>> module spline2D()
>>   {
>>     B = nSpline(transpose(A), 100); // interpolate X   //<<<
>>     C = nSpline(transpose(B), 100); // interpolate Y
>>
>>    // using cubes  and very(!!) slow F6 rendering
>>    translate([-5, 0, 0])
>>      surf(A);                // plot coarse data
>>    translate([1, 0, 0])
>>      surf(C); // plot interpolated data
>>
>>     // using triangulation with polyhedron fast F6 rendering
>>     translate([-4.5, -5, 0])
>>       surfaceplot(A);                // plot coarse data
>>     translate([1, -5, 0])
>>       surfaceplot(C); // plot interpolated data   //<<<
>>   }
>>
>>   module surfaceplot(A, d=[0,5,0,5])    //<<<
>> // d =[ xmin, xmax, ymin, ymax ]
>>   {
>>     map = [for(i=[0:len(A)-1], j=[0:len(A[0])-1])
>> [d[0]+i*(d[1]-d[0])/(len(A)-1), d[2]+j*(d[3]-d[2])/(len(A[0])-1),
>> A[i][j]]];    //<<<
>>     base = [for(i=[0:len(A)-1], j=[0:len(A[0])-1])
>> [d[0]+i*(d[1]-d[0])/(len(A)-1), d[2]+j*(d[3]-d[2])/(len(A[0])-1), 0]];
>> //<<<
>>
>>     faces = concat(faces(A), faces(A, len(A)*len(A[0])), sidefaces(A));
>>     polyhedron(concat(map, base), faces);
>>
>>     // used for top and bottom faces
>>     function faces(A, offs=0) =   let(X = len(A), Y= len(A[0]))
>>       [for(i=[0:X-2], j=[0:Y-2], k=[1,0]) let(m=i*Y+j+offs)
>>         offs? k?[m+1, m, m+Y]:[m+1, m+Y, m+Y+1] : k?[m, m+1, m+Y]:[m+1,
>>   m+Y+1, m+Y]];
>>
>>     function sidefaces(A) =
>>       let(X = len(A), Y= len(A[0]))
>>       let(A=[for(i=[0:Y-2], k=[1,0]) k?[i+1, i, i+X*Y]: [i+1, i+X*Y,
>>   i+1+X*Y]])
>>       let(B=[for(i=[0:Y-2], k=[1,0]) let(m=i+(X-1)*Y) k?[m, m+1, m+X*Y]:
>>   [m+1, m+1+X*Y, m+X*Y]])
>>       let(C=[for(i=[0:X-2], k=[1,0]) let(m=i*Y) k?[m, m+Y, m+X*Y]: [m+Y,
>>   m+Y+X*Y, m+X*Y]])
>>       let(D=[for(i=[0:X-2], k=[1,0]) let(m=(i+1)*Y-1) k?[m, m+X*Y, m+Y]:
>>   [m+Y, m+X*Y, m+Y+X*Y]])
>>       concat(A, B, C, D);
>>
>>   }
>
>
>
>
> --
> View this message in context: http://forum.openscad.org/Bezier-Courves-elevation-model-tp20778p20787.html
> Sent from the OpenSCAD mailing list archive at Nabble.com.
>
> _______________________________________________
> OpenSCAD mailing list
> [hidden email]
> http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org
>
>
>
> -----
> No virus found in this message.
> Checked by AVG - www.avg.com
> Version: 2016.0.7998 / Virus Database: 4756/14073 - Release Date: 03/07/17
>
>


_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org
Reply | Threaded
Open this post in threaded view
|

Re: Bezier Courves - elevation model

cbernhardt
In reply to this post by rob63
rob63:
Since your example looks like a DEM (digital elevation mode) of the surface of the earth, I will tell you how to solve this problem the way land surveyors do it, but you probably will not like the answer.  I was a land surveyor for many years and wrote code in AutoCad to solve this problem.  Basically surveying data is just a group of x,y,z coordinates, not necessarily, and usually not, in a grid pattern.  I used the software to calculate volumes of coal piles and computing earthwork volumes on large grading projects.

The process involves creating a TIN (triangular irregular network) from the x,y,z coordinate data.  If you do a search for TIN,  Renka, and Delaunay you will probably find some sample code that you may be able to adapt to your particular situation.

Personally I would not even attempt to try programming this method in OpenSCAD.  Below is a model of a large coal pile I created in AutoCad.  There were over 26,000 data points and the TIN was created in about 2 seconds on a 32 bit 3.0GHz machine.
Reply | Threaded
Open this post in threaded view
|

Re: Bezier Courves - elevation model

Ronaldo
@cbernhardt

Delaunay triangulation is standard for scattered data. But as stated, droftarts input data are distributed on a regular grid (whose Delaunay triangulation is trivial). I have tried to code a Delaunay triangulation method in OpenSCAD. The hardest part is to implement the triangulation data structure with just simple lists. I gave up.

2017-03-07 17:59 GMT-03:00 cbernhardt <[hidden email]>:
rob63:
Since your example looks like a DEM (digital elevation mode) of the surface
of the earth, I will tell you how to solve this problem the way land
surveyors do it, but you probably will not like the answer.  I was a land
surveyor for many years and wrote code in AutoCad to solve this problem.
Basically surveying data is just a group of x,y,z coordinates, not
necessarily, and usually not, in a grid pattern.  I used the software to
calculate volumes of coal piles and computing earthwork volumes on large
grading projects.

The process involves creating a TIN (triangular irregular network) from the
x,y,z coordinate data.  If you do a search for TIN,  Renka, and Delaunay you
will probably find some sample code that you may be able to adapt to your
particular situation.

Personally I would not even attempt to try programming this method in
OpenSCAD.  Below is a model of a large coal pile I created in AutoCad.
There were over 26,000 data points and the TIN was created in about 2
seconds on a 32 bit 3.0GHz machine.
<http://forum.openscad.org/file/n20793/tin.jpg>



--
View this message in context: http://forum.openscad.org/Bezier-Courves-elevation-model-tp20778p20793.html
Sent from the OpenSCAD mailing list archive at Nabble.com.

_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org


_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org
Reply | Threaded
Open this post in threaded view
|

Re: Bezier Courves - elevation model

cacb
In reply to this post by cbernhardt
On 07. mars 2017 21:59, cbernhardt wrote:

> The process involves creating a TIN (triangular irregular network) from the
> x,y,z coordinate data.  If you do a search for TIN,  Renka, and Delaunay you
> will probably find some sample code that you may be able to adapt to your
> particular situation.
>
> Personally I would not even attempt to try programming this method in
> OpenSCAD.  Below is a model of a large coal pile I created in AutoCad.
> There were over 26,000 data points and the TIN was created in about 2
> seconds on a 32 bit 3.0GHz machine.
> <http://forum.openscad.org/file/n20793/tin.jpg>

Very cool. Does this dataset or similar data exist on the net for download?

Carsten Arnholm

_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org
Reply | Threaded
Open this post in threaded view
|

Re: Bezier Courves - elevation model

Parkinbot
In reply to this post by Ronaldo
Ronaldo wrote
I have found the origin of the NaNs. In your spline code, you make the spline knot intervals proportional to the input data segment length.
Yes that must be the cause. Thanks. Now I remember. I implemented the code to specifically interpolate trajectories of N-dimensional points (for N=3 this is a polyline in 3D). To get a better coupling of the dimensions (and smoother splines) I used a specific version of a line integral that is also used in Matlabs nSpline(). This approach is of course not the best for Rob's use case. But for now, let's see how it works with Rob's data.
Reply | Threaded
Open this post in threaded view
|

Re: Bezier Courves - elevation model

Neon22
meshalab - free - will do poisson disk construction of a surface from a pile of vertex cooordinates. You could try that...
Reply | Threaded
Open this post in threaded view
|

Re: Bezier Courves - elevation model

Ronaldo
Playing with Rudolf's code, it became apparent to me that natural splines may not the best choice for droftarts' problem.
If we take the following input data:

 A = [  // coarse data of surface
   [  1,   2, 0.1,   1,  2],
   [  2,   2,   4,   2,  1],
   [0.2, 0.2,   2,   1,  3],
   [  1,   1,   1,   1,  4],
 ];


the overshoot tendency of natural spline will produce a graph with negative values (although all data values are positive) and the graph volume itself has self-intersections. As this is in the very nature of the natural splines I would not recommend it to general terrain modelling.



There are some spline interpolation methods that keep the solution in the bounds of the input data but I never studied them.

_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org
Reply | Threaded
Open this post in threaded view
|

Re: Bezier Courves - elevation model

Ronaldo
Wrong example, sorry. That is the corrected one:

 A = [  // coarse data of surface
   [  1,   2, 0.1,   1,  2],
   [  4,   3,   4,   2,  1],
   [0.1, 0.1,   3,   2,  3],
   [  1,   2,   1,   1,  4],
 ];





_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org