ベジエ曲線の制御点を指定して曲率を求めたい

Question

実現したいこと

ベジエ曲線の曲率を求められるようにしたい。

制御点を指定して作図するプログラムはできたのですが、そこから曲率を求めるプログラムの作成に苦戦してます。
以下のサイトを参考に曲率半径の計算を行おうとしています。
https://qiita.com/HMMNRST/items/39f77e5051df2928eb93

該当のソースコード

import numpy as np
import math
from matplotlib import pyplot as plt
import os
import glob
from PIL import Image

# Bernstein多項式を計算する関数
def bernstein(n, t):
    B = []
    for k in range(n + 1):
        # 二項係数を計算してからBernstein多項式を計算
        nCk = math.factorial(n) / (math.factorial(k) * math.factorial(n - k))
        B.append(nCk * t ** k * (1 - t) ** (n - k))
        print(nCk, k, n - k)
    return B


# ベジェ曲線を描く関数
def bezie_curve(Q):
    n = len(Q) - 1
    dt = 0.01
    t = np.arange(0, 1 + dt, dt)
    B = bernstein(n, t)
    px = 0
    py = 0
    for i in range(len(Q)):
        px += np.dot(B[i], Q[i][0])
        py += np.dot(B[i], Q[i][1])
    return px, py, t


# 点座標を準備

q1 = [0., 0]
q2 = [0.5, 1.]
q3 = [0.75, 0.]
q4 = [1.25, 1.]
q5 = [1.5, 0.]
q6 = [2., 1.]

Q = [q1, q2, q3, q4, q5, q6]

# ベジェ曲線を描く関数を実行
px, py, t = bezie_curve(Q)


for j in range(len(t)):
    print(j)
    # 制御点個数-1（n-1)の線分ができるので、その方程式（傾きと切片）を得る
    tx = []
    ty = []
    for i in range(len(Q) - 1):
        diff_x = Q[i + 1][0] - Q[i][0]
        if diff_x == 0:
            diff_y = Q[i + 1][1] - Q[i][1]
            tx.append(Q[i][0])
            ty.append(Q[i][1] + t[j] * diff_y)
        else:
            grad = (Q[i + 1][1] - Q[i][1]) / diff_x
            intercept = Q[i][1] - grad * Q[i][0]
            tx_temp = Q[i][0] + t[j] * diff_x
            tx.append(tx_temp)
            ty.append(grad * tx_temp + intercept)
    tx2 = []
    ty2 = []
    for i in range(len(tx) - 1):
        diff_x2 = tx[i + 1] - tx[i]
        if diff_x2 == 0:
            diff_y2 = ty[i + 1] - ty[i]
            tx2.append(tx[i])
            ty2.append(ty[i] + t[j] * diff_y2)
        else:
            grad2 = (ty[i + 1] - ty[i]) / diff_x2
            intercept2 = ty[i] - grad2 * tx[i]
            tx_temp2 = tx[i] + t[j] * diff_x2
            tx2.append(tx_temp2)
            ty2.append(grad2 * tx_temp2 + intercept2)

曲率の式をどのように適用すればいいのか三時間くらいネットサーフィンをして探したり考えたりしてますが、いまだにできてません。
どうか教えてください。

補足情報（FW/ツールのバージョンなど）

python 3.10

Answer 1

あらかじめバーンスタイン法によってベジエ曲線を構成する点を計算済み（≒微分済み）であると考えれば、
隣接する3点を通る円の半径の逆数を、その点を通る曲線の曲率と見なせばいいのではないでしょうか。

（略）
def get_curvature(x1, y1, x2, y2, x3, y3):
    # 3点(x1,y1) (x2,y2) (x3,y3) を通る円の半径の逆数を計算
    d = ((y1 - y3) * (x1 - x2) - (y1 - y2) * (x1 - x3)) * 2
    x = ((y1 - y3) * (y1 ** 2 - y2 ** 2 + x1 ** 2 - x2 ** 2) - (y1 - y2) * (y1 ** 2 - y3 ** 2 + x1 ** 2 - x3 ** 2)) / d
    y = ((x1 - x3) * (x1 ** 2 - x2 ** 2 + y1 ** 2 - y2 ** 2) - (x1 - x2) * (x1 ** 2 - x3 ** 2 + y1 ** 2 - y3 ** 2)) / -d
    r = math.sqrt((x - x1) ** 2 + (y - y1) ** 2)
    if r == 0:
        return 0
    return 1 / r

# 最初の端点における曲率＝曲線を構成するために分割して得た点のうち、
# 最初の3点を通る円の半径の逆数
print(get_curvature( px[0],py[0],
                               px[1],py[1],
                               px[2],py[2]))

Comments

[Tyutohannpa_]

ありがとうございます！曲率計算ができるようになりました。
しかし、そのまま適用するのと
以下のように直接、値を指定すると違う値が表示されるのですが
どうしてでしょうか
def get_curvature():
# 3点(x1,y1) (x2,y2) (x3,y3) を通る円の半径の逆数を計算
(x1, y1) = [0., 0]
(x2, y2) = [1.75, 1.]
(x3, y3) = [1.75, 1.]
～以下略～

スキル不足ですみません。

[退会済みユーザー]

get_curvature 関数には制御点の座標ではなく、曲線を構成する、隣接する3点を指定します。

ある曲線の、ある点における曲率＝「その曲線にその点で接する円の半径の逆数」なので。

したがって、曲線のスタートである端点における曲率は、
端点と、端点の1つ隣の点と、そのまた１つ隣の点　この３つの点を通る円の半径の逆数で近似できます。
回答したコードで
print(get_curvature( px[0],py[0],
px[1],py[1],
px[2],py[2]))
としているのですから
get_curvature に渡している引数
px[0],py[0]
px[1],py[1]
px[2],py[2]
の座標をprintしてみてください。
いずれもx座標どうし、y座標どうしは非常に近い値になっていると思います。これは「隣り合っている」点（最初の3つの点）だからです。

[Tyutohannpa_]

では、次の3つの点の曲率を調べたいときは
px[1],py[1]
px[2],py[2]
px[3],py[3]
とすればいいですか？

[退会済みユーザー]

そうですね。「スタートの端点の１つ隣から３つ」の部分での曲率を得るならばそうなります。

ちなみに元質問のコードは t=0.01なので、曲線は101個の点をつなげた100個の線分でできています。
したがって、終端点における曲率を知りたいなら、
get_curvage(px[98],py[98],px[99],py[99],px[100],py[100])
とすればよいです。
（曲線を構成する点は、０番目～100番目の101個なので）

[Tyutohannpa_]

分かりやすい解説ありがとうございました！！
とても満足いたしました！