Re: [OpenSCAD] Bezier Courves - elevation model
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
https://en.wikipedia.org/wiki/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
https://en.wikipedia.org/wiki/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 https://en.wikipedia.org/wiki/Spline_interpolation . This
last method is implemented in splines.scad of thingi:1208001
http://www.thingiverse.com/thing:1208001 . Study the included example and
come back if you have doubts.
--
View this message in context: http://forum.openscad.org/Bezier-Courves-elevation-model-tp20778p20783.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
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([di, dj]) 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]) [id, jd, A[i][j]]];
base = [for(i=[0:len(A)-1], j=[0:len(A[0])-1]) [id, jd, 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+XY]: [i+1, i+XY,
i+1+XY]])
let(B=[for(i=[0:Y-2], k=[1,0]) let(m=i+(X-1)Y) k?[m, m+1, m+XY]:
[m+1, m+1+XY, m+XY]])
let(C=[for(i=[0:X-2], k=[1,0]) let(m=iY) k?[m, m+Y, m+XY]: [m+Y,
m+Y+XY, m+XY]])
let(D=[for(i=[0:X-2], k=[1,0]) let(m=(i+1)Y-1) k?[m, m+XY, m+Y]:
[m+Y, m+XY, m+Y+X*Y]])
concat(A, B, C, D);
}
--
View this message in context: http://forum.openscad.org/Bezier-Courves-elevation-model-tp20778p20786.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
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+XY]])
let(B=[for(i=[0:Y-2], k=[1,0]) let(m=i+(X-1)Y) k?[m, m+1, m+XY]:
[m+1, m+1+XY, m+XY]])
let(C=[for(i=[0:X-2], k=[1,0]) let(m=iY) k?[m, m+Y, m+XY]: [m+Y,
m+Y+XY, m+XY]])
let(D=[for(i=[0:X-2], k=[1,0]) let(m=(i+1)Y-1) k?[m, m+XY, m+Y]:
[m+Y, m+XY, 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.
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.
--
View this message in context: http://forum.openscad.org/Bezier-Courves-elevation-model-tp20778p20788.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
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" rudolf@parkinbot.com escreveu:
Thanks Ronaldo,
well, matrix A was already defined in a transposed fashion ;-)
Ronaldo wrote
If A is already transposed and the data plot used as such,
Certainly you are right.
--
View this message in context: http://forum.openscad.org/Bezier-Courves-elevation-model-tp20778p20790.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
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.
--
View this message in context: http://forum.openscad.org/Bezier-Courves-elevation-model-tp20778p20791.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
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+XY]])
let(B=[for(i=[0:Y-2], k=[1,0]) let(m=i+(X-1)Y) k?[m, m+1, m+XY]:
[m+1, m+1+XY, m+XY]])
let(C=[for(i=[0:X-2], k=[1,0]) let(m=iY) k?[m, m+Y, m+XY]: [m+Y,
m+Y+XY, m+XY]])
let(D=[for(i=[0:X-2], k=[1,0]) let(m=(i+1)Y-1) k?[m, m+XY, m+Y]:
[m+Y, m+XY, 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
Discuss@lists.openscad.org
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
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.
--
View this message in context: http://forum.openscad.org/Bezier-Courves-elevation-model-tp20778p20796.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
meshalab - free - will do poisson disk construction of a surface from a pile
of vertex cooordinates. You could try that...
--
View this message in context: http://forum.openscad.org/Bezier-Courves-elevation-model-tp20778p20797.html
Sent from the OpenSCAD mailing list archive at Nabble.com.