3次元ベジェ曲線のバウンティボックスを求める

ベジエ曲線バウンディングボックス
http://d.hatena.ne.jp/nishiohirokazu/20090616/1245104751


InkscapeSVG解析ですね、ええわかります。
Inkscape SVGは独自拡張記法で回転中心という項があるんですよね、しかし、バウンティボックスという拡張記法は無いんです。
まったく、気の利かない!(逆切れ)


実は私もInkscapeSVGのツールを作っていまして、3次ベジェ曲線のバウンティボックスを求める関数を作っていましたので、共有できればなと・・・


さて、3次のベジェ曲線のバウンティボックス問題なんですが、
下がベジェの方程式です。

f(t) = (1-t)**3*P1 + 3*t*(1-t)**2*P2 + t**2*3*(1-t)*P3 + t**3*P4

バウンティボックス問題は簡単に言ってしまうと「3次方程式の極値」を求めるということです。

極値とはすなわち接線が0の地点ですので、接線を求めるためにベジェ関数を微分します。

f'(t) = 3*t**2*P4-3*t**2*P3+6*(1-t)*
           t*P3-6*(1-t)*t*P2+3*(1-t)**2*
               P2-3*(1-t)**2*P1

接線が0のときのtを求めればいいのですから、方程式「f'(t)=0」のtについて求めます。
2次方程式なので解が2つあります)

t=( sqrt((P1-P2)*P4+P3**2+(-P2-P1)*P3+P2**2)-P3+2*P2-P1)
          /(P4-3*P3+3*P2-P1)


t=(-sqrt((P1-P2)*P4+P3^2+(-P2-P1)*P3+P2^2)-P3+2*P2-P1)
          /(P4-3*P3+3*P2-P1)

ここで注意すべきは虚数解が出る可能性があるという事です。
虚数解の場合は特別な処理が必要です、以下の部分がPythonのプログラムとなります。

# coding:utf-8
import math

def __getBezierMinMax(P1,P2,P3,P4):
    """ベジェ曲線の極値を求める"""
    a = P4 - 3*P3 + 3*P2 - P1
    b = (P1-P2)*P4 + P3**2 + (-P2-P1)*P3+P2**2
    c = P3 - 2*P2 + P1
    v = [P1,P4]

    if a == 0.0:
        d = float(2*(P3-2*P2+P1))
        if d != 0.0:
            t = -(P2-P1) / d
            if 0.0 < t < 1.0:
                v.append(__calcBezier(t,P1,P2,P3,P4))
        
    elif b >= 0.0:
        bb = math.sqrt(b)
        t  = (bb-c)/float(a)
        if 0.0 < t < 1.0:
            v.append(__calcBezier(t,P1,P2,P3,P4))
        t  = -(bb+c)/float(a)
        if 0.0 < t < 1.0:
            v.append(__calcBezier(t,P1,P2,P3,P4))
    return min(v), max(v)

def __calcBezier(t,P1,P2,P3,P4):
    """ベジェ曲線を求める"""
    return (1-t)**3*P1 + 3*t*(1-t)**2*P2 + t**2*3*(1-t)*P3 + t**3*P4

たぶん、あってるはず・・・

※追記
重解の場合を考えていませんでした、この問題は意外と奥が深かった・・・
指摘をしてくだっさった西尾さん、ありがとうございました。

追記

このブログを見た人は、私が数字に強い人だと思うんでしょうねぇ。
次は、私は中学高校と数学の偏差値は30代を叩き出したほどの猛者です。
テストは殆ど名前を書いただけwwwwwwwwwwwwww


「じゃあ、上の計算はどうやったんだ?」といいますと、
wxMaximaというフリーの数値計算ソフトを使いましたwwフヒヒ、サーセンwwwww
ベジェ曲線の問題程度なら、3行打ち込むだけで終了です。

f(t):=(1-t)**3*P1 + 3*t*(1-t)**2*P2 + t**2*3*(1-t)*P3 + t**3*P4;
diff(f(t)=0, t);
solve([%], [t]);

世の中便利になりました、感謝。


■wxMaxima
http://wxmaxima.sourceforge.net/wiki/index.php/Main_Page