

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


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.


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: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:X2], j=[0:Y2], 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:Y2], k=[1,0]) k?[i+1, i, i+X*Y]: [i+1, i+X*Y, i+1+X*Y]])
let(B=[for(i=[0:Y2], k=[1,0]) let(m=i+(X1)*Y) k?[m, m+1, m+X*Y]: [m+1, m+1+X*Y, m+X*Y]])
let(C=[for(i=[0:X2], 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:X2], k=[1,0]) let(m=(i+1)*Y1) k?[m, m+X*Y, m+Y]: [m+Y, m+X*Y, m+Y+X*Y]])
concat(A, B, C, D);
}


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:X2], j=[0:Y2], 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:Y2], k=[1,0]) k?[i+1, i, i+X*Y]: [i+1, i+X*Y,
i+1+X*Y]])
let(B=[for(i=[0:Y2], k=[1,0]) let(m=i+(X1)*Y) k?[m, m+1, m+X*Y]:
[m+1, m+1+X*Y, m+X*Y]])
let(C=[for(i=[0:X2], 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:X2], k=[1,0]) let(m=(i+1)*Y1) k?[m, m+X*Y, m+Y]:
[m+Y, m+X*Y, m+Y+X*Y]])
concat(A, B, C, D);
}


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 howtouse 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.


Ronaldo wrote
If A is already transposed and the data plot used as such,
Certainly you are right.


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 righthand 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[i1][k]),2))
k==0 ? t1 : t1 + max(lineT1(S, i, k1),1e12);
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 nonuniform spacing. Anyway, if two subsequent data points are equal, the repetition might/should be eliminated as a nonsense.


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:X2], j=[0:Y2], 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:Y2], k=[1,0]) k?[i+1, i, i+X*Y]: [i+1, i+X*Y,
>> i+1+X*Y]])
>> let(B=[for(i=[0:Y2], k=[1,0]) let(m=i+(X1)*Y) k?[m, m+1, m+X*Y]:
>> [m+1, m+1+X*Y, m+X*Y]])
>> let(C=[for(i=[0:X2], 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:X2], k=[1,0]) let(m=(i+1)*Y1) 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/BezierCourveselevationmodeltp20778p20787.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


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.


@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.
_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org


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


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 Ndimensional 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.


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


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 selfintersections. 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


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

