裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

家紋シリーズ 三つ兎

2021年09月09日 | ブログラミング

家紋シリーズ 三つ兎

自由曲線をベジェ曲線で近似する

# ベジェ曲線
# https://tajimarobotics.com/basis-spline-interpolation-program-2/

function bezier(P; points=50)
    spline(P, size(P, 1) - 1; points)
end

function spline(P, n; points=50)
    # p = size(P, 1)  # Number of Control Ponits
    # n = 3           # Degree of B-Spline
    # n = p -1        # Bezier curve
    p = size(P, 1)
    m = p + n + 1     # Number of Knots in Knot Vector

    # Define Knot Vector
    u = OpenUniformKnotVector(m, n)
    # t = 0:0.01:u[end]
    # t = range()
    t = range(u[1], stop=u[end], length=points)

    # Calculate B-spline Curve
    S = zeros(length(t), 2);
    S[1, :] = P[1, :]
    for i ∈ 2:length(t)
        for j ∈ 1:p
            b = BasisFunction(u, j, n, t[i])
            S[i, :] += P[j, :] * b
        end
    end
    S
end

# u : Knot Vector
# m : Number of Knots in Knot Vector
# n : Degree of B-Spline
function OpenUniformKnotVector(m, n)
    u = zeros(1, m)
    for i ∈ 1:m
          if i < n+1
              u[i] = 0
          elseif i > m - (n + 1)
              u[i] = m - 1 - 2 * n
          else
              u[i] = i - n - 1
          end
    end
    u ./ u[end]
end

# u : Knot Vector
# j : j-th Control Ponit
# k : k-th Basis Function <-- n-th B-Spline
# t : Time
function BasisFunction(u, j, k, t)
    w1, w2 = 0, 0

    if k == 0
        if u[j] < t <= u[j+1]
            var = 1
        else
            var = 0
        end
    else
        if u[j+k+1] != u[j+1]
            w1 = BasisFunction(u, j+1, k-1, t) *
                              (u[j+k+1] - t)/(u[j+k+1]-u[j+1])
        end
        if u[j+k] != u[j]
            w2 = BasisFunction(u, j, k-1, t) *
                              (t - u[j])/(u[j+k] - u[j])
        end
        var = w1 + w2
    end
    var
end

# plotter.jl を include
# https://blog.goo.ne.jp/r-de-r/e/bd71a52a09801335d56f7c47d879bfe3

include("plotter.jl")

function plotwhitecurve(ox, oy, r, t, xy)
    xy = bezier(xy)
    plotwhiteline(ox, oy, r, t, xy)
end

function plotwhiteline(ox, oy, r, t, xy)
    xy2 = r .* xy * t
    plotline(ox .+ xy2[:,1], oy .+ xy2[:, 2], col=:white, lwd=2)
    xy2 = xy; xy2[:, 1] *= -1
    xy2 = r .* xy * t
    plotline(ox .+ xy2[:,1], oy .+ xy2[:, 2], col=:white, lwd=2)
end

function ploteyes(ox, oy, r, t, xy, radius, color)
    xy2 = r .* xy *t
    for i ∈ 1:3
        plotcircle(ox + xy2[i, 1], oy + xy2[i, 2],
                   radius[i], lwd=0, fcol=color[i])
    end
    xy2 = xy; xy2[:, 1] *= -1
    xy2 = r .* xy *t
    for i ∈ 1:3
        plotcircle(ox + xy2[i, 1], oy + xy2[i, 2],
                   radius[i], lwd=0, fcol=color[i])
    end
end

function plothead(ox, oy, r, t)
    # 頭
    xy = bezier([0 692; 36 661; 67 633; 81 626; 138 597;
                 232 576; 237 580; 290 550; 342 530; 392 440]);
    # 耳(上)
    xy = vcat(xy, bezier([250 400; 388 425; 395 440; 392 445;
                          484 430; 494 400; 589 347]));
    # 耳(下)
    xy = vcat(xy, bezier([589 347; 500 290; 426 260; 250 270]));
    # 頬
    xy = vcat(xy, bezier([427 273; 395 157; 346 104; 272 61;
                          188 50; 159 22; 89 22; 5 5; 0 8]));
    xy2 = r .* xy * t
    plotpolygon(ox .+ xy2[:,1], oy .+ xy2[:, 2],
                lwd=0.5, col=:black, fcol=:black)
    xy2 = xy; xy2[:, 1] *= -1
    xy2 = r .* xy2 * t
    plotpolygon(ox .+ xy2[:,1], oy .+ xy2[:, 2],
                lwd=0.5, col=:black, fcol=:black)
end

function usagi1(R, r, θ)
    ox, oy = R * cosd(θ), R * sind(θ)
    θ += 90
    t = [cosd(θ) sind(θ); -sind(θ) cosd(θ)]
    plothead(ox, oy, r, t)
    # おでこ
    plotwhitecurve(ox, oy, r, t,
        [144 363; 112 394; 84 398; 53 412; 0 412])
    # 耳上白線
    plotwhitecurve(ox, oy, r, t,
        [96 329; 131 348; 159 365; 191 379; 230 394; 263 410;
        310 423; 352 435; 382 437; 422 437; 456 431; 489 418])
    # 耳下白線
    plotwhitecurve(ox, oy, r, t,
        [143 269; 145 269; 177 280; 198 286; 226 286; 254 273;
        297 265; 350 258; 392 265; 424 268; 452 274])
    # 目
    ploteyes(ox, oy, r, t,
        [74 216; 73 215; 72 214],
        [32, 22, 12], [:white, :black, :white])
    # 鬚
    plotwhiteline(ox, oy, r, t,
        [89 174; 191 220; NaN NaN; 81 160; 244 178; NaN NaN;
        74 142; 205 142])
    # 鼻
    plotwhitecurve(ox, oy, r, t,
        [-7 104; 0 100; 21 111; 18 125; 17 139; 36 157])
    # 上唇
    plotwhitecurve(ox, oy, r, t,
        [14 114; 32 100; 43 93; 57 90; 74 93; 89 100; 106 107])
    # 下唇
    plotwhitecurve(ox, oy, r, t,
        [0 68; 14 67; 29 75; 36 79; 43 93])
    # 頬
    plotwhitecurve(ox, oy, r, t,
        [99 8; 106 45; 117 60; 131 72; 149 86; 157 89; 180 97])
    # 足
    plotwhiteline(ox, oy, r, t,
        [134 15; 134 44; NaN NaN; 170 26; 170 44])
end

# メインプログラム
function usagi(; r=1, width=400, height=400)
    plotbegin(w=width, h=height)
    R = 725
    for θ ∈ [-30, 90, 210]
        usagi1(R, r, θ)
    end
    #usagi1(r, -0)
    plotend()
end

usagi()

savefig("fig.png")

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村