# Bezier Courves - elevation model

17 messages
Open this post in threaded view
|

## Bezier Courves - elevation model

 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.htmlMy 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
Open this post in threaded view
|

## Re: Bezier Courves - elevation model

 Would 'Surface' do the trick? https://en.wikibooks.org/wiki/OpenSCAD_User_Manual/Other_Language_Features#SurfaceIt'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
Open this post in threaded view
|

## Re: Bezier Courves - elevation model

 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.
Open this post in threaded view
|

## Re: Bezier Courves - elevation model

 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  // http://www.thingiverse.com/thing:1208001A = [  // 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);         }
Open this post in threaded view
|

## Re: Bezier Courves - elevation model

 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);  }
Open this post in threaded view
|

## Re: Bezier Courves - elevation model

 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.
Open this post in threaded view
|

## Re: Bezier Courves - elevation model

 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
Open this post in threaded view
|

## Re: Bezier Courves - elevation model

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

## Re: Bezier Courves - elevation model

 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.
Open this post in threaded view
|

## Re: Bezier Courves - elevation model

Open this post in threaded view
|

## Re: Bezier Courves - elevation model

 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.
Open this post in threaded view
|

## Re: Bezier Courves - elevation model

Open this post in threaded view
|

## Re: Bezier Courves - elevation model

Open this post in threaded view
|

## Re: Bezier Courves - elevation model

 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.
Open this post in threaded view
|

## Re: Bezier Courves - elevation model

 meshalab - free - will do poisson disk construction of a surface from a pile of vertex cooordinates. You could try that...