bezier.mac

/*

https://github.com/t-o-k/Maxima-bezier/bezier.mac

Copyright (c) 2020 Tor Olav Kristensen, http://subcube.com

Use of this source code is governed by the GNU Lesser General Public License version 3, which can be found in the LICENSE file.

*/

load(bernstein);

matrix_sum_1a(m) :=
    block(
        [ size_m, size_s, s_ones ],
        size_m: matrix_size(m),
        size_s: size_m[2],
        s_ones: matrix(makelist(1, i, 1, size_s)),
        m.transpose(s_ones)
    )
;

matrix_sum_2a(m) :=
    block(
        [ size_m, size_u, size_v, u_ones, v_ones ],
        size_m: matrix_size(m),
        size_u: size_m[2],
        size_v: size_m[1],
        u_ones: matrix(makelist(1, i, 1, size_u)),
        v_ones: matrix(makelist(1, i, 1, size_v)),
        v_ones.m.transpose(u_ones)
    )
;

basis_functions_1a(deg_s, s) :=
    matrix(
       makelist(bernstein_poly(i, deg_s, s), i, 0, deg_s)
    )
;

basis_functions_2a(deg_u, deg_v, u, v) :=
   block(
       [ u_fns, v_fns ],
       u_fns(u) := basis_functions_1a(deg_u, u),
       v_fns(v) := basis_functions_1a(deg_v, v),
       transpose(v_fns(v)).u_fns(u)
   )
;

weighted_basis_functions_1a(w, s) :=
    block(
        [ size_w, deg_s ],
        size_w: matrix_size(w),
        deg_s: size_w[2] - 1,
        w*basis_functions_1a(deg_s, s)
   )
;

weighted_basis_functions_2a(w, u, v) :=
    block(
        [ size_w, deg_u, deg_v ],
        size_w: matrix_size(w),
        deg_u: size_w[2] - 1,
        deg_v: size_w[1] - 1,
        w*basis_functions_2a(deg_u, deg_v, u, v)
    )
;

bezier_function_1a(p, s) :=
    matrix_sum_1a(
        weighted_basis_functions_1a(p, s)
    )
;

bezier_function_2a(p, u, v) :=
    matrix_sum_2a(
        weighted_basis_functions_2a(p, u, v)
    )
;

rational_bezier_function_1a(p, w, s) :=
    block(
        [ wb ] ,
        wb: weighted_basis_functions_1a(w, s),
        matrix_sum_1a(p*wb)/matrix_sum_1a(wb)
    )
;

rational_bezier_function_2a(p, w, u, v) :=
    block(
        [ wb ] ,
        wb: weighted_basis_functions_2a(w, u, v),
        matrix_sum_2a(p*wb)/matrix_sum_2a(wb)
    )
;