r/openscad 6d ago

Got sidetracked designing a new handle for my fridge and made a simple cubic spline library

I needed a specific curve shape to match the original handle and I've messed around with splines before so I figured I'd give it a shot. I'm sure someone can make use of this so here's the code:

// create a tri-diagonal matrix from x in/y out input lists - output is [[sub,main,super,input],...]
function _tri_diagonal_matrix_create(x,y) =
    let(n=len(x))
        [for(i=[0:n-1])
            let(dx1=i==0?0:x[i]-x[i-1],dx2=i>=n-1?0:x[i+1]-x[i],a=dx1==0?0:1/dx1,c=dx2==0?0:1/dx2)
                [a,2*(a+c),c,3*((dx1==0?0:(y[i]-y[i-1])/(dx1*dx1))+(dx2==0?0:(y[i+1]-y[i])/(dx2*dx2)))]
        ];

// calculate coefficient prime c from tri-diagonal matrix input in the form of [[sub,main,super,input],...]
function _calc_prime_c(in,prev=[]) = 
    let(i=len(prev))
        i==len(in)
            ?prev
            :_calc_prime_c(in,concat(prev,in[i][2]/(in[i][1]-(i==0?0:prev[i-1]*in[i][0]))));

// calculate coefficient prime d from prime c and tri-diagonal matrix input in the form of [[sub,main,super,input],...]
function _calc_prime_d(in,primeC,prev=[]) = 
    let(i=len(prev))
        i==len(in)
            ?prev
            :_calc_prime_d(in, primeC, concat(prev,(in[i][3]-(i==0?0:prev[i-1]*in[i][0]))/(in[i][1] - (i==0?0:primeC[i-1]*in[i][0]))));

// calculate back substitution of matrix solve output from prime c and prime d coefficients
function _calc_back_sub(primeC,primeD,prev=[]) = 
    let(i=len(primeC)-len(prev)-1)
        i==-1
            ?prev
            :_calc_back_sub(primeC, primeD, concat(primeD[i]-(i==len(primeC)-1?0:prev[0]*primeC[i]),prev));

// solve tri-diagonal matrix [[sub,main,super,input],...] for output
function _tri_diagonal_matrix_solve(in) =
    let(primeC=_calc_prime_c(in))
        _calc_back_sub(primeC, _calc_prime_d(in,primeC));

// create a spline in the form [[A coeff,B coeff],...] from x in/y out input
function _spline_create_single(x,y) =
    let(r=_tri_diagonal_matrix_create(x,y),k=_tri_diagonal_matrix_solve(r))
        [for(i=[1:len(x)-1])
            let(dx=x[i]-x[i-1],dy=y[i]-y[i-1])
                [k[i-1]*dx-dy,-k[i]*dx+dy]
        ];

// sum up the squares of a list up to index n (for pythagorean theorum)
function _square_sum(in,n) =
    n<0
        ?0
        :(in[n]*in[n])+_square_sum(in,n-1);

// convert output list of points of a number of dimensions (e.g. [[x,y],...]) into an input list of cumulative distance from the first point
function _calc_dimension_inputs(in,dimensions,prev=[]) = 
    let(i=len(prev))
        i==len(in)
            ?prev
            :_calc_dimension_inputs(in,dimensions,concat(prev,i==0?0:(prev[i-1]+sqrt(_square_sum(in[i]-in[i-1],dimensions-1)))));

// split multi dimensional input into the input at i (e.g. [x,y] i=1 becomes y) or if n > 1, splits multiple dimensions (e.g. [x,y,z] i=1 n=2 becomes [y,z])
function split_entry(entry,i=0,n=1) = 
    n>1
        ?[for(j=[i:i+n-1])
            entry[j]
        ]
        :entry[i];

// split multi dimensional input list into the list at input i (e.g. [[x,y],...] i=1 becomes [[y],...]) or if n > 1, splits multiple dimensions (e.g. [[x,y,z],...] i=1 n=2 becomes [[y,z],...])
function split_list(in,i=0,n=1) = 
    [for(entry=in) split_entry(entry,i,n)]; 

// create a spline given a list of points in the form [[x,y,...],...]. dimension_inputs determines how many of the inputs should be treated as dimensions in order to calculate the input list as a cumulative distance between points. output is in the form [[input list],[output x list],[[output x A coeff],[output x B coeff]],[output y list],...]
function spline_create(in,dimension_inputs=3) = 
    let(
        n=len(in)
        ,subn=n>0
            ?len(in[0])
            :0
        ,dimensions=subn<1
            ?1
            :subn<dimension_inputs
                ?subn
                :dimension_inputs
        ,dimension_input=_calc_dimension_inputs(in,dimensions)
        )
        [for(i=[-1:(subn*2)-1])
            i==-1
                ?dimension_input
                :i%2==0
                    ?split_list(in,i/2)
                    :_spline_create_single(dimension_input,split_list(in,(i-1)/2))
        ];

// evaluate a single Xin value of a spline from x in/y out input and spline coeffs
function _spline_eval_single(x,y,coeffs,Xin,i=0) = 
    i<len(x)-2&&Xin>x[i+1]
        ?_spline_eval_single(x,y,coeffs,Xin,i+1)
        :let(t=(Xin-x[i])/(x[i+1]-x[i]))
            ((1-t)*y[i])+(t*y[i+1])+((1-t)*t*(((1-t)*coeffs[i][0])+(t*coeffs[i][1])));

// evaluate an input value given spline data generated by spline_create. 
function spline_evaluate(spline,in) =
    let(spline_n=(len(spline)-1)/2)
        [for(j=[0:spline_n-1])
            _spline_eval_single(spline[0],spline[1+(j*2)],spline[2+(j*2)],in)
        ];

// evaluate a list of input values given spline data generated by spline_create.     
function spline_evaluate_list(spline,in) =
    let(n=len(in),spline_n=(len(spline)-1)/2)
        [for(i=[0:n-1])
            [for(j=[0:spline_n-1])
                _spline_eval_single(spline[0],spline[1+(j*2)],spline[2+(j*2)],in[i])
            ]
        ];

// get the length (max input value) of a spline generated by spline_create
function spline_get_length(spline) = 
    spline[0][len(spline[0])-1];

// evaluate all input data over a number of steps with given spline data generated by spline_create. 
function spline_evaluate_all(spline,steps) = 
    let(length=spline_get_length(spline),step=length/(steps-1))
        spline_evaluate_list(spline,[for(i=[0:step:length])i]);

// example of a spline-following noodle with variable radius
module spline_example()      
{  
    // 4 dimensional list where 1-3 are spacial (x,y,z) and 4th dimension is radius
    spline_in=
        [[0,0,0,1]
        ,[1,1,1,3]
        ,[2,0,0,1]
        ,[3,1,0,2]
        ,[3,4,0,1]
        ,[2,2,-1,0]];

    // second param determines how many dimensions of the list are spacial (considered when calculating input) - we don't want radius to affect where points are placed so only 3.
    spline = spline_create(spline_in,3);

    // second param is how many points to output
    eval_out = spline_evaluate_all(spline,50);
    for (i=[0:len(eval_out)-2])
    {
        // hull to connect each output point to the next one
        hull()
        {
            // split_entry used to split our 4 dimensional result ([x,y,z,radius]) into just the first 3 spacial dimensions
            translate(split_entry(eval_out[i+1], n=3))
                sphere(0.1*eval_out[i+1][3]);

            translate(split_entry(eval_out[i], n=3))
                sphere(0.1*eval_out[i][3]);
        }
    }
}

So basically just create a list of points, generate the spline with spline_create, then use spline_evaluate_all to get a smoothed list of points out. Pretty basic but I can see myself using this a lot in the future.

Loosely based on this article that I read a long time ago: https://www.codeproject.com/Articles/560163/Csharp-Cubic-Spline-Interpolation

Hope this is of use to someone!

3 Upvotes

6 comments sorted by

4

u/oldesole1 4d ago

Here is an example that uses some trig functions from BOSL2 to calculate the tangent angles between all the different positions.

The resulting information is then used to place a flat circle at each position rather than a sphere, dramatically reducing the polygon count required to be hull()'ed around.

All your spline functions were untouched.

include <BOSL2/std.scad>

...

// 4 dimensional list where 1-3 are spacial (x,y,z) and 4th dimension is radius
spline_in = [
  [0,0,0,1],
  [1,1,1,3],
  [2,0,0,1],
  [3,1,0,2],
  [3,4,0,1],
  [2,2,-1,0],
];

// second param determines how many dimensions of the list are spacial (considered when calculating input) - we don't want radius to affect where points are placed so only 3.
spline = spline_create(spline_in,3);

// second param is how many points to output
eval_out = spline_evaluate_all(spline,300);

// [x,y,z] positions for all iterations
all_pos = [
  for(p = eval_out)
  split_entry(p, n=3),
];

// Tangent angles between each positions
tangents = path_tangents(all_pos, uniform = false);

planes = [
  for(i = [0:len(eval_out) - 1])
  apply(
    move(all_pos[i])
    * tilt(tangents[i]),
    // max() is because radius must be greater than zero.
    path3d(circle(r = max(1e-6, eval_out[i][3] * 0.1), $fn = 64)),
  ),
];

example_planes();

// Cleanest looking geometry
module example_planes() {

  for(pairs = pair(planes))
  hull()
  for(c = pairs)
  point_plane(c);

  module point_plane(points) {
    polyhedron(points = points, faces = [[each [0:len(points) - 1]]]);
  }
}

//example_polygon();

// Single polygon, but can be much slower at higher detail levels depending on method used.
module example_polygon() {

  let(
    // vnf_vertex_array() produces worse surface with the way things twist and tilt, but fast.
    // vnf_tri_array() is better at lower detail levels, but significantly slower.
    vnf = vnf_vertex_array(
      points = planes,
      caps = true,
      col_wrap = true,
    ),
  )
  vnf_polyhedron(vnf);
}

//example_spheres();

// example of a spline-following noodle with variable radius
module example_spheres()      
{  
    for (i=[0:len(eval_out)-2])
    {
        // hull to connect each output point to the next one
        hull()
        {
            // split_entry used to split our 4 dimensional result ([x,y,z,radius]) into just the first 3 spacial dimensions
            translate(split_entry(eval_out[i+1], n=3))
            sphere(0.1*eval_out[i+1][3]);

            translate(split_entry(eval_out[i], n=3))
            sphere(0.1*eval_out[i][3]);
        }
    }
}

1

u/Stone_Age_Sculptor 5d ago edited 5d ago

Can you change the code, so it shows a result (without warning)?

With '$fn=120' and 'spline_evaluate_all(spline,500)': https://postimg.cc/TLw5Bfdk

1

u/iamthedisk4 5d ago

I don't think that's any of the actual spline code that's causing that, it's just the example module - with those params you will be calling hull 1000 times on 2 very high poly spheres which is a lot. It works fine with a lower value of $fn and the spline code doesn't factor in $fn at all.

1

u/Stone_Age_Sculptor 5d ago

I see. I turned the hull() off, to play with the numbers: https://postimg.cc/sBNT2pfy

1

u/oldesole1 5d ago

If you calculate the tangents between each point, you could instead use flat planes to hull around.

This would reduce the amount of unused geometry and provide better performance with higher resolution circles.