3
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@kan0u(狩野 宏樹)in株式会社イワタ

3次Bézier曲線の弧長の求め方、計り取り方

3
Last updated at Posted at 2026-05-22

はじめに

本記事は、株式会社イワタのnote記事『「イワタつづり」インク溜まり生成自動化スクリプト設計の裏側とフローを徹底図解』（前半）（後半）の補足、技術的な詳細について説明したものです。

note記事で解説したスクリプト中で頻出する処理で、アルゴリズムの説明を省略した Bézier 曲線 (Bézier curve) の弧長 (arc length) を求める方法と、与えられた長さになる Bézier 曲線の分割位置を求める方法を説明します。

弧長を求める

再帰的分割による方法

一番原始的な方法は、曲線を細かく分割して折れ線近似する方法です。曲線の長さは弦の長さ (chord length, $L_c$)、すなわち開始点と終了点の間の直線距離より長く、端点と制御点を順次繋いだ折れ線の長さ (polygon length, $L_p$)より短いのは明らかなので、この二つの値を求めれば求める値を上下から挟み込むことができます。この両者の値の中間値を取れば、より正確な近似になることが期待できます。

分割方法としては、よく知られた de Casteljau のアルゴリズムにより $t = 0.5$ で等分にする操作を再帰的に行って、$2^n$ 分割することにします。

# 点 p, q は x 座標と y 座標の 2 個の浮動小数からなるタプル
def midpoint(p, q):
    return ((p[0]+q[0]) / 2.0, (p[1]+q[1]) / 2.0)

# 3 次 Bezier 曲線を t = 0.5 で前後に等分する
# 返り値の [0]..[3] が分割曲線の前半、[3]..[6] が後半部分
def bisect(p0, p1, p2, p3):
    p01 = midpoint(p0, p1)
    p12 = midpoint(p1, p2)
    p23 = midpoint(p2, p3)
    p012 = midpoint(p01, p12)
    p123 = midpoint(p12, p23)
    p0123 = midpoint(p012, p123)
    return (p0, p01, p012, p0123, p123, p23, p3)

孤長の下限・上限は折れ線の長さとしてして求められます。

# ベクトルと有向線分 (directed line segment) のルーチンから抜粋
def vecAbs(v): return math.hypot(v[0], v[1])
def vecSub(v, w): return (v[0] - w[0], v[1] - w[1])
def dlsLen(p, q): return vecAbs(vecSub(q, p))

def bezierChordLen(bez):
    return dlsLen(bez[0], bez[3])
    
def bezierPolygonLen(bez):
    return dlsLen(bez[0], bez[1]) + dlsLen(bez[1], bez[2]) + dlsLen(bez[2], bez[3])

# 再帰的に n 回 2 分割して、弧長の下限と上限を返す
def bezierArcLenBounds(bez, numBisectionNeeded):
    if numBisectionNeeded <= 0:
        bounds = (bezierChordLen(bez), bezierPolygonLen(bez))
        return bounds
    else:
        b = bisect(bez[0], bez[1], bez[2], bez[3])
        firstHalf, secondHalf =  (b[0:4], b[3:7])
	firstBounds = bezierArcLenBounds(firstHalf, numBisectionNeeded - 1)
        secondBounds = bezierArcLenBounds(secondHalf, numBisectionNeeded - 1)
        bounds = ( firstBounds[0] + secondBounds[0], firstBounds[1] + secondBounds[1] )
        return bounds

以下のようなテストコードで、分割の段数を変えると精度と実行速度がどう変わるかを調べます。なお、テストに使用した曲線は 2 次 Bézier 曲線（すなわち放物線）を 3 次に昇格したものですので弧長を解析的に求めることができ、

L=240(\sqrt2+\log_2(1+\sqrt2)) \approx 550.9409158542331

となります。

# 2 次 Bézier 曲線 ((0,0), (480, 480), (960, 0)) を 3 次に昇格
seg = ((0,0), (160, 160), (320, 240), (480, 240))

prev = time.perf_counter()
for n in range(20):
    sup, inf = bezierArcLenBounds(seg, n)
    now = time.perf_counter()
    milliseconds = (now - prev) * 1000
    prev = now
    avg = (inf+sup) / 2.0
    print(f"{n}|{inf:.13f}|{sup:.13f}|{avg:.13f}|{milliseconds}")

結果は以下のようになります（上限と下限の中間値としては算術平均を採用しています）。

分割段数 弧長の下限 ($L_c$) 弧長の上限 ($L_p$) 平均 ($(L_c+L_p)/2$) 実行時間
0 565.1596081796783 536.6563145999495 550.9079613898139 5 μs
1 554.4846357021841 547.3863375370596 550.9354866196218 14 μs
2 551.8253162993053 550.0558265761100 550.9405714377076 11 μs
3 551.1619190781427 550.7198694974149 550.9408942877787 17 μs
4 550.9961605943457 550.8856684171158 550.9409145057307 29 μs
5 550.9547266599903 550.9271048798947 550.9409157699425 56 μs
6 550.9443685319657 550.9374631659640 550.9409158489648 102 μs
7 550.9417790221846 550.9400526856232 550.9409158539039 204 μs
8 550.9411316461284 550.9407000622967 550.9409158542126 392 μs
9 550.9409698022012 550.9408619062626 550.9409158542319 767 μs
10 550.9409293412248 550.9409023672413 550.9409158542330 1.51 ms
11 550.9409192259810 550.9409124824851 550.9409158542330 3.12 ms
12 550.9409166971701 550.9409150112962 550.9409158542331 6.07 ms
13 550.9409160649674 550.9409156434990 550.9409158542333 12.1 ms
14 550.9409159069168 550.9409158015496 550.9409158542333 24.4 ms
15 550.9409158674041 550.9409158410623 550.9409158542333 48.4 ms
16 550.9409158575259 550.9409158509404 550.9409158542331 96.6 ms
17 550.9409158550563 550.9409158534099 550.9409158542331 194 ms
18 550.9409158544390 550.9409158540274 550.9409158542333 386 ms
19 550.9409158542846 550.9409158541816 550.9409158542331 776 ms

表中の長さの数値が解析解と一致する桁は太字で表示しています。下限・上限とも 1 回分割を増やしても 1 ビット程度しか精度が上がらないのに対し、両者の平均は遥かに早く収束します（この例では曲線が特に滑らかなため、分割段数 10 で精度が頭打ちになっていますが、一般的な場合でも収束するには単精度で 6 段、倍精度で 15 段程度繰り返すのが無難です）。実際、分割数を増やして行ったとき、3次 Bézier 曲線の $L_c$ と $L_p$ の算術平均は $L$ に収束し、一般に、$n$次 Bézier 曲線の弧長は、$(2L_c+(n-1)L_p)/(n+1)$ になることが証明されています1。丸め誤差の影響のため、12 段目以降では値が振動しており、倍精度の最後の1桁まで正確な値を得るには long double (80bit の拡張精度または 4 倍精度) で演算を行う必要があります。

この方法は、制御点の位置が通常フォントのアウトラインには使わないような極端な配置になっていて、尖点やループを持つような場合でも精度の高い解が得られるロバストなアルゴリズムです。しかしながら、そのような場合でも十分な精度を得るには $2^{15} = 32768$ 回の分割処理と、その 8 倍の距離計算が必要になります。距離を計り取るには二分法などの方法で試し切りした曲線の長さを何度も計る必要があるので、遅すぎて実用に耐えません。

数値積分を用いた方法

より正攻法と言えるのが、速度ベクトルの絶対値を積分して移動距離を求める方法です。Bézier曲線の弧長を求める場合、Gauss求積を用いる方法が好んで用いられます。最適なサンプル点における関数値にそれぞれ異なる重みを掛けて足すだけで、台形法のような素朴な方法よりも圧倒的に少ない個数のサンプル点で高い精度が得られるからです。あらかじめ決められたサンプル点における関数値を求め、それに決まった重みを掛けて足すだけで超高精度な積分値が得られるという魔法のような方法ですが、それがうまくいく理由は「任意の関数を多項式近似できる直交多項式の性質をうまく利用しているから」という説明でごまかしておきます（詳しくは数値計算の教科書2を参照してください）。直交多項式としてはLegendere多項式（またはChebyshev多項式）を用いることが多いです。

そのような実装の中で最も高精度と思われたのが、Raph Levien 氏の Blog3 で紹介されていた、24 次の Gauss-Legendre 求積法を用いた弧長計算法で、Pomax 氏の bezier.js で実装されています4。同氏による FAQ5 も参照し、3 次 Bézier 曲線に限定して単純化した Python コードが以下のものになります。

# Gauss-Legendere求積 (24次) のサンプル点と重み係数 (対称性を用いて半分に圧縮)                     
GL_Abscissae = ( # 24次Legendere多項式の根
    0.0640568928626056, 0.1911188674736163, 0.3150426796961634, 0.4337935076260451,
    0.5454214713888395, 0.6480936519369756, 0.7401241915785544, 0.8200019859739029,
    0.8864155270044010, 0.9382745520027328, 0.9747285559713095, 0.9951872199970214
)
GL_Weights = ( # cf. https://pomax.github.io/BezierInfo-2/index.html#arclength                              
    0.1279381953467522, 0.1258374563468283, 0.1216704729278034, 0.1155056680537256,
    0.1074442701159656, 0.0976186521041139, 0.0861901615319533, 0.0733464814110803,
    0.0592985849154367, 0.0442774388174198, 0.0285313886289336, 0.0123412297999871
)
# https://github.com/Pomax/bezierjs/blob/master/src/utils.js を単純化
def bezierArcLength(p0, p1, p2, p3):
    v01, v12, v23 = (vecSub(p1, p0), vecSub(p2, p1), vecSub(p3, p2))
    sum = 0.0
    rep = len(GL_Abscissae) * 2;
    for i in range(rep):
        sign = -1.0 if i % 2 == 0 else 1.0
        t = 0.5 * sign * GL_Abscissae[i>>1] + 0.5 # Gauss求積の定義域は[-1,1]なのを[0,1]に合わせる
        T = 1.0 - t
        a, b, c = (T * T, 2.0 * T * t, t * t)
    	dx = v01[0] * a + v12[0] * b + v23[0] * c
        dy = v01[1] * a + v12[1] * b + v23[1] * c
        sum += GL_Weights[i>>1] * math.hypot(dx, dy)
    return 1.5 * sum

先ほどの曲線例で相対誤差を求めるコードは以下のようになります。

def relativeError(std, value):
    return abs(value - std) / std

seg = ((0,0),　(160,160),　(320,240),　(480,240))
bezlen = bezierArcLength(*seg)
print("total length:", bezlen)
print("relative error:", relativeError(550.9409158542331, bezlen))

結果:

total length: 550.940915854233
relative error: 2.063503262329768e-16

この例の場合、特に近似しやすい形をしているので倍精度浮動小数点数の限界に近い精度を叩き出していますが、歪んだ形のBézier曲線、例えば 2 回の曲率極大点をもつ Z 字形やループを含むようなベジェ曲線では、誤差ははるかに大きくなります。例えば図 1 に示すような曲線の長さをそのまま求めると、相対誤差が $10^{-4}$ のオーダーになってしまいます。1000 メッシュを超える長い曲線だと 1 メッシュ近くズレるおそれがあり、実用的にも見過ごせないレベルです。

fig1.png
図1: 極度に精度が落ちるケース。
曲率極大点で分割して前後別々に長さを計算すると、精度が劇的に改善する

図1のケースの計算:

seg = ((0, 0), (14, 2), (2, -2), (5, 1))

sup, inf = bezierArcLenBounds(seg, 15)
avg = (sup + inf)/2
print("segment:", seg)
print("length:    ", avg)
approx = bezierArcLength(*seg)
print("approx len:", approx)
print("relative error",relativeError(approx, avg))

結果:

segnent: ((0, 0), (14, 2), (2, -2), (5, 1))
length:     10.773841997177922
approx len: 10.775190330611412
relative error 0.00012513314309263015

曲率極大点が2個以上ある場合は、以下のように、最大の曲率をもつ点（図1の赤丸の箇所）の位置で分割すると、劇的に精度が改善します。

def subdividedTwinBezierCurves(p0, p1, p2, p3, t):
    T = 1.0 - t
    p01   = (p0[0]*T + p1[0]*t, p0[1]*T + p1[1]*t)
    p12   = (p1[0]*T + p2[0]*t, p1[1]*T + p2[1]*t)
    p23   = (p2[0]*T + p3[0]*t, p2[1]*T + p3[1]*t)
    p012  = (p01[0]*T + p12[0]*t, p01[1]*T + p12[1]*t)
    p123  = (p12[0]*T + p23[0]*t, p12[1]*T + p23[1]*t)
    p0123 = (p012[0]*T + p123[0]*t, p012[1]*T + p123[1]*t)
    return (p0, p01, p012, p0123, p123, p23, p3)

t0 = 0.3746995425
print("subdivide curve at t =", t0)
twin = subdividedTwinBezierCurves(*seg, t0)
first, second = (twin[0:4], twin[3:7])
print("length:    ", avg)
approx = bezierArcLength(*first) + bezierArcLength(*second)
print("approx len:", approx)
print("relative error", relativeError(approx, avg))

分割計算結果:

subdivide curve at t = 0.3746995425
length:     10.773841997177922
approx len: 10.773841997177922
relative error 0.0

今回は性質の良い曲線のみを扱うと前提することができたので分割処理は行っていないのですが、汎用ルーチン化するならば、パスの形状チェックと分割処理を行う必要があるでしょうし、メインの計算処理はオーバースペック気味（フォントの処理に使うなら、相対誤差の最大値が単精度浮動小数のマシンイプシロンと同程度、およそ$10^{-7}$に設定するのがバランスが取れていると思います）なので、Legendre多項式の次数を下げて高速化した方が良いかと思います。

長さの計り取り方

3 次 Bézier 曲線の弧長 $L$ より短い長さ $l\ (0\lt l\le L)$ が与えられたとき、$p_0$ から始まる弧長が $l$ に一致するような分割点$t\ (0\lt t\le 1)$を求めるには、解析的には解けないので実際に適当な推定点で切って、断片の長さを $l$ と比較し、順次推定の精度を上げていく必要があります。弧長の変化率を求めて Newton 法を使うことも考えられますが、例外的なケース（途中の推定値がで尖点に一致するとか）の対応が面倒なので、関数の連続性のみに依存する求根アルゴリズムを使用しています。本記事のサンプルでは最も単純な二分法を示し、実際に採用しているアルゴリズム (Brent法などのいくつかの方式を比較し、最も効率的だった 修正 Anderson-Bjorck 法という方式を採用) の紹介は、また機会があれば別の記事で行いたいと思います。

# 二分法の繰り返し回数の最大値 (仮数部のビット数が上限) 
BISECTION_METHOD_REPEAT_MAX = 23  # 単精度なら最大 23, 倍精度なら最大 53

NUM_FUNC_CALLED = 0 # ベンチマーク用に、関数の呼び出し回数を数えるカウンタ

def findRoot(func, initial_low_estimate = 0.0, initial_high_estimate = 1.0):
    global NUM_FUNC_CALLED
    lo, hi = (initial_low_estimate, initial_high_estimate)
    vl, vh = (func(lo), func(hi))
    NUM_FUNC_CALLED += 2    
    if vl * vh > 0:
        return None
    dt = hi - lo

    for i in range(BISECTION_METHOD_REPEAT_MAX):
        dt /= 2.0
        mid = lo + dt
        d = func(mid)
        NUM_FUNC_CALLED += 1
        if dt == 0.0 or d == 0.0:
            return mid
        elif d * vl < 0.0:
            hi = mid
        else:
            lo = mid
    return (lo + hi) / 2.0

# Bezier 曲線の 2 分割 (de Casteljau のアルゴリズムに基づく) 結果の前半
def subdividedBezierCurveBefore(p0, p1, p2, p3, t):
    T = 1.0 - t
    p01   = (p0[0]*T + p1[0]*t, p0[1]*T + p1[1]*t)
    p12   = (p1[0]*T + p2[0]*t, p1[1]*T + p2[1]*t)
    p23   = (p2[0]*T + p3[0]*t, p2[1]*T + p3[1]*t)
    p012  = (p01[0]*T + p12[0]*t, p01[1]*T + p12[1]*t)
    p123  = (p12[0]*T + p23[0]*t, p12[1]*T + p23[1]*t)
    p0123 = (p012[0]*T + p123[0]*t, p012[1]*T + p123[1]*t)
    return (p0, p01, p012, p0123)
    
# 分割点の前の弧長
def bezierArcLengthBefore(p0, p1, p2, p3, t):
    P0, P1, P2, P3 = subdividedBezierCurveBefore(p0, p1, p2, p3, t)
    return bezierArcLength(P0, P1, P2, P3)

# 指定した長さ len で 3 次 Bezier曲線を計り取る
def bezierDistantTimeFromStart(p0,p1,p2,p3, len, precision=1.0e-6):
    distFunc = lambda t: len - bezierArcLengthBefore(p0, p1, p2, p3, t)
    return findRoot(distFunc, 0.0, 1.0, precision)

参照文献

  1. Jens Gravesen, "Adaptive subdivision and the length and energy of Bézier curves", Comput. Geom. 8 (1997) 13-31 ↩

  2. 例えば、水島二郎・柳瀬眞一郎『理工学のための数値計算法』（サイエンス社）2002, pp. 31-40. ↩

  3. Raph Levien, "How long is your Bézier?" (Dec 28, 2018) in Raph Levien’s blog ↩

  4. https://github.com/Pomax/bezierjs/blob/master/src/utils.js ↩

  5. Pomax, "§24. Arc length" in A Primer on Bézier Curves ↩

3
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
3
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?