スタッフBLOG

近頃、天候の具合が思わしくなく、困っています。
特に、休日に雨となると、布団を干したり、シーツやタオルケットを洗ったりするのが厄介ですよね。
先週末と今週末が雨続きだったので、洗濯待ちのシーツが部屋の隅で山になっていてゲンナリです。

もしも厳しい御姑さんなんかがいようものなら、
「まあ！Ｍ子さん！こんなにシーツを溜め込んだりして！まったく、ご実家ではどういう教育を受けてらっしゃったのかしら！？」
とか、軽く下っ腹に効くローなブローをかまされちゃって、
「よよよょ。申し訳ありません、御義母様。Ｍ子、至らなくて。次からはきちんと洗濯しますから、今日のところはこのシーツを洗濯機にぶち込んで洗いあがったところでベランダに干して乾いた加減を見計らって取り込んでキチンと畳んでしまっておいて下さいね。いえ、アイロンかけやがって下さいましとまでは言いませんわ、御義母様次第ですもの。ではＭ子、お仕事いってきま～す」
とか、危うく嫁姑問題に心を砕くハメになっていたところですよ。
あ～独身でよかった（この頬を伝うのは雨粒だよ！）。

ポリアンナをレスペクトする本日の当番のM.L.Kです（LはLamdaのLです。Lamda&Pi!)。
おはようございます。

さて、メディアンカットでございます。

前回は、主成分分析による分割の実装に届かないところで終わってしまいました。
ですので、今回は実装編です！
と、言いながら当番がまわって来る間に、よくよく考えてみたのですが、いわゆる普通のメディアンカットとは実装についての根本的な部分が違うので、もう少しプログラミング的な事柄について、扱ってみたいと思います。

と言うワケで、実装準備編です！
もしかしたら、実装準備編その１かもしれません！悪しからず！

で、前回、固有値λを２次方程式を求めて、その値を代入した連立方程式を解けば、固有ベクトルが求まりますよ～、的なことを言っていたと思うのですが、冷静に考えると、実際は３チャンネルの分割をしないといけないワケですから、３次方程式を解いて、３元の連立方程式を解かなくてはならないワケですよ！

うわ、煩わし！

つ～ワケで、とりあえず、３次方程式の解を求めるコードを書いてみます。

まず、立方根を求められることが必要条件ですので、これを実装します。
何故これが必要かって～と、３次方程式の解の計算では、立方根を求めなければならないのですが、標準関数にはこれが用意されていないからです。
まあ、pow関数を利用するという手もあるのですが、速度と精度の問題から、余り望ましいとは言えないので、敢えて実装するというワケです。

//立方根（実数）を計算する関数
//引数　：x = 立方根を求める値
//戻り値：求められた立方根
double cbrt( double x )
{
　　double dStart = 1.0;
　　bool bPositive = false;

　　//収束の向きを一致させるため、負数は正数に変換する
　　//（値を求めた後で、元の符号に戻します）
　　if( x > 0.0 )
　　{
　　　　bPositive = true;
　　}
　　else if( x < 0.0 )
　　{
　　　　bPositive = false;
　　　　x = -x;
　　}
　　else
　　{
　　　　return 0.0;
　　}

　　if( x > 1.0 )
　　{
　　　　dStart = x;
　　}

　　double dPrev = 0.0;

　　do
　　{
　　　　dPrev = dStart;
　　　　double dRecSq = 1 / (dStart*dStart);
　　　　dStart = ( 2.0 * dStart + x * dRecSq ) * (1.0/3.0);

　　} while( dStart < dPrev );

　　if( bPositive )
　　{
　　　　return dPrev;
　　}

　　return -dPrev;

}

上は、ニュートン法による求め方です。
主に、do～whileの中がそれですね。

ある数Aについての実数の立方根を求めるとは、すなわち、(X-A)^3の実数解を求めることに他ならないので、

f(X) = X^3 - A

とすると、その１次導関数は、

f'(X) = 3X^2

となり、ニュートン法の式は

Xnew = X - f(X)/f'(X)

すなわち

Xnew = X - (X^3 - A)/(3X^2)
Xnew = X - (1/3)(X - A/X^2)
Xnew = (1/3)(3X - X + A/X^2)
Xnew = (1/3)(2X + A/X^2)

となり、求められたXnewを右辺のXに代入することを繰り返せば、次第にＡの立方根に収束していく、という寸法ですね。

続いて、３次方程式の解の計算です。
立方根をニュートン法で求めるんなら、いっそ初めから３次方程式自体をニュートン法で求めりゃイイじゃん、とか仰るイヤに鋭い方もいらっしゃるかもしれませんが、今回は複数の解が欲しいことと、３次方程式は係数の値によっては、収束の振動が厄介なので、ここでは“カルダノの方法”を用います。

ちょっと脱線ですが、皆さんは“カルダノの方法”と記述する時、
「あれ、“タルタリアの方法”だったっけ？」
とか迷ったりしません？
私はよく、どっちが約束を破ったほうだったけ？ってなります。
“タルタルソース”→“食べもの”→“食いものにされた”→“だまされた”
と憶えているつもりだったのですが、実は、
“カルダモン”→“食べもの”→“食いものにされた”→“だまされた”
という連想もできたりするので、結局ワケが分からなくなります。

って、３次方程式の解法についての逸話を知らないと、それこそ全くワケが分からない話題なんですが。

閑話休題、以下実装です。

#define SQRT3 (1.7320508075688772935274463415059)
#define PI (3.1415926535897932384626433832795)

//３次方程式の解を求める関数
//引数　：pDst = 求められた解を格納する配列のポインタ（３要素以上）
//　　　　dA = ３次項の係数
//　　　　dB = ２次項の係数
//　　　　dC = １次項の係数
//　　　　dD = 定数項
//戻り値：求められた解の個数
int CalcCubicEquation( double* pDst,
　　　　　　　　　　　 double dA,
　　　　　　　　　　　 double dB,
　　　　　　　　　　　 double dC,
　　　　　　　　　　　 double dD )
{
　　if( pDst == NULL ) return 0;

　　double dRecA = 1.0 / dA;
　　double dB0 = dB * dRecA * (1.0/3.0);
　　double dC0 = dC * dRecA;
　　double dD0 = dD * dRecA;

　　double dP0 = dB0 * dB0 - dC0 * (1.0/3.0);
　　double dQ0 = ( dB0 * ( dC0 - 2.0 * dB0 * dB0 ) - dD0 ) * (1.0/2.0);
　　double dA0 = dQ0 * dQ0 - dP0 * dP0 * dP0;

　　if( dA0 > 0.0 )
　　{
　　　　double dA3 = 0.0;
　　　　if( dQ0 > 0 )
　　　　{
　　　　　　dA3 = cbrt( dQ0 + sqrt( dA0 ) );
　　　　}
　　　　else
　　　　{
　　　　　　dA3 = cbrt( dQ0 - sqrt( dA0 ) );
　　　　}

　　　　double dB3 = dP0 * (1 / dA3);
　　　　double dX1 = dA3 + dB3 - dB0;
　　　　double dX2 = -(1.0/2.0) * (dA3 + dB3) - dB0;
　　　　double dX3 = fabs(dA3 - dB3) * SQRT3 * (1.0/2.0);

　　　　pDst[0] = dX1;
　　　　pDst[1] = dX2;
　　　　pDst[2] = dX3;

　　　　return 3;
　　}
　　else if( dA0 < 0.0 )
　　{
　　　　double dA0 = sqrt( dP0 );
　　　　double dT = acos( dQ0 * (1.0 / (dP0*dA0)) );
　　　　double dA1 = dA0 + 2.0;
　　　　double dX1 = dA1 * cos( dT * (1.0/3.0) ) - dB0;
　　　　double dX2 = dA1 * cos( (dT + 2.0 * PI) * (1.0/3.0) ) - dB0;
　　　　double dX3 = dA1 * cos( (dT + 4.0 * PI) * (1.0/3.0) ) - dB0;

　　　　pDst[0] = dX1;
　　　　pDst[1] = dX2;
　　　　pDst[2] = dX3;

　　　　return 3;
　　}
　　else
　　{
　　　　double dQ1 = cbrt( dQ0 );
　　　　double dX1 = 2.0 * dQ1 - dB0;
　　　　double dX2 = -dQ1 - dB0;

　　　　pDst[0] = dX1;
　　　　pDst[1] = dX2;

　　　　return 2;
　　}
}

実装的には、“カルダノの方法”そのまんまです。
特に工夫らしいことはしていません。
なにかしら、精度や速度を上げられるのでは、と思って検討してみるのですが、こういった、元より洗練されている手合いは、下手に粗をほじくると逆に酷い目に遭ったりしがちです。

Ｍ子が掃除を終えた後に姑が目ざとく障子の横組子に人差し指を滑らせて、
「あらあら、Ｍ子さん、これは…」
ペロッ
「青酸カリっ！」

みたいな。

なんか、今回は前回以上にメディアンカットから外れた事柄の扱いとなってしまいました。しかも、まだ準備終わってないっぽいんだな～、これが…。