Catmull-Rom Spline

[usao]

参加しましたのでちょっとしたものを．

Catmull-Rom Spline．
全ての制御点の上を通る曲線(の上の点)が簡単に求まるので便利．
（数式とか曲線具合は検索すると見つかります．）

コード:

//Catmull-Rom Splineの，ある区間(隣接する制御点間)のカーブ（上の離散点位置）を計算する．
//
//・制御点[i]～制御点[i+1]までの区間を，"i番目"の区間とする．
//・両端の制御点位置での傾きは関数内でてきとーに決定されるが，
//  自分で制御したい場合には，両端の外側にもう一個それ用のダミー制御点を用意すればできると思われ．
//
//< template引数 >
//  VEC     : 制御点の座標を表すための型．例えば２次元や３次元のベクトル．
//            この型は，関数内で使う以下のoperator等をサポートしていなくてはならない．
//            ・VEC型同士の代入，コピーコンストラクタ
//            ・VEC型同士の加算(+, +=)，減算(-)
//            ・VEC型 = VEC型 * double型定数  乗算（* : 定数は右からかける）
//
//[戻り値]
//  成功時はtrue,失敗時はfalseを返す．
//  なお，falseを返した場合，結果格納コンテナの内容状態はいじられない．
template < class VEC >
bool CalcCRSplineCurveSegment(
    const std::vector<VEC> &rControlPoints, //制御点列
    size_t SegmentIndex,    //計算する区間を指定．( 0 <= val <= rControlPoints.size()-2 )
    size_t nCurvePoints,    //求める点の個数（区間両端の制御点を含む個数）( 3 <= val )
    std::vector<VEC> &rDstCurvePoints   //結果座標群受取用
    )
{
    if( rControlPoints.empty()  ||  nCurvePoints<3 )return false;
    const size_t nCPs = rControlPoints.size();
    if( SegmentIndex+1 >= nCPs )return false;

    const VEC &P0 = rControlPoints[SegmentIndex];
    const VEC &P1 = rControlPoints[SegmentIndex+1];
    const VEC Slope0 = ( P1 - (   SegmentIndex>0   ?   rControlPoints[SegmentIndex-1]   :   P0   ) ) * 0.5;
    const VEC Slope1 = ( (   SegmentIndex+2 < nCPs   ?   rControlPoints[SegmentIndex+2]   :   P1   ) - P0 ) * 0.5;

    // t = 0.0 ～ 1.0 : 区間内での位置を表す媒介変数
    //
    //                             |  1  0   0   0 | |   P0   |
    // Pos(t) = [ 1  t  t^2  t^3 ] |  0  0   1   0 | |   P1   |
    //                             | -3  3  -2  -1 | | Slope0 |
    //                             |  2 -2   1   1 | | Slope1 |
    const VEC Tmp[3] = {
        //P0,
        Slope0,
        (P0 * -3.0) + (P1 * 3.0)  + (Slope0 * -2.0) - Slope1,
        (P0 * 2.0)  + (P1 * -2.0) + Slope0 + Slope1
    };

    rDstCurvePoints.resize( nCurvePoints );
    const size_t nMinus1 = nCurvePoints - 1;
    rDstCurvePoints[0] = P0;
    rDstCurvePoints[nMinus1] = P1;

    double coeff[3];// = { 1.0, 0,0,0 };

    for( size_t i=1;    i<nMinus1;  i++ )
    {
        coeff[0] = (double)i / nMinus1; //t
        coeff[1] = coeff[0] * coeff[0]; //t^2
        coeff[2] = coeff[1] * coeff[0]; //t^3

        VEC Pos = P0;
        for( int j=0; j<3; j++ )
        {
            Pos += ( Tmp[j] * coeff[j] );
        }
        rDstCurvePoints[i] = Pos;
    }

    return true;
}