Mugichoko's blog

Mugichoko’s blog

プログラミングを中心としたメモ書き.

スプライン曲線の生成から描画まで

3次スプライン曲線 (Cubic Spline Curve)

2次元座標(図中,赤色)を与えると,その間を,3次スプライン補間により補間した点(図中,緑色)を描画するプログラム.数値計算用及び描画用ライブラリとしてOpenCV 2.Xを利用しています.3次元関数を想定しているため,円などは描けません.

f:id:Mugichoko:20180405184817g:plain

main.cpp

//
// 2015.8.18
// Mugichoko
//
#include <iostream>
using namespace std;

#include "CubicSpline.h"

const int wndSize = 512;

int main()
{
    vector<cv::Point2f> vData;
    vData.push_back(cv::Point2f(0.0f, 0.0f));
    vData.push_back(cv::Point2f(100.0f, 200.0f));
    vData.push_back(cv::Point2f(130.0f, 50.0f));
    vData.push_back(cv::Point2f(255.0f, 255.0f));
    CubicSpline<cv::Point2f> myCubicSpline(vData);

    cv::Mat img = cv::Mat::zeros(wndSize, wndSize, CV_8UC3);

    cv::Point2f shift(wndSize / 4, wndSize / 4);
    // Draw control points
    for (int idx = 0; idx < vData.size(); ++idx)
        cv::circle(img, vData[idx] + shift, 5, cv::Scalar(0, 0, 255), -1);
    cv::imshow("VIZ", img);
    cv::waitKey();
    // Draw interpolated points
    for (int i = 0; i < vData.size() - 1; ++i) {
        for (int j = 0; j < 10; ++j) {
            cv::Point2f v = myCubicSpline.getValue(i, j / 10.0f);
            cv::circle(img, v + shift, 2, cv::Scalar(0, 255, 0), -1);

            cv::imshow("VIZ", img);
            cv::waitKey(100);
        }
    }
    cv::circle(img, myCubicSpline.getValue((unsigned int)vData.size() - 2, 1.0f) + shift, 2, cv::Scalar(0, 255, 0), -1);
    cv::imshow("VIZ", img);
    cv::waitKey(2000);

    return 0;
}

CubicSpline.h

//
// 2015.8.18
// Mugichoko
//
#ifndef CUBICSPLINE_H
#define CUBICSPLINE_H

#include <iostream>
using namespace std;
#include <opencv2/opencv.hpp>

struct CubicSplineCoeff
{
    float a, b, c, d;
};

//
// Ref. 1: http://www.ipc.akita-nct.ac.jp/yamamoto/lecture/2003/5E/lecture_5E/Lagrange_Spline.pdf
//
template<typename vec2D>
class CubicSpline
{
private:
    vector<CubicSplineCoeff> m_vCoeff;
    vector<vec2D> m_vData;

    float calc_h(unsigned int idx)
    {
        return this->m_vData[idx + 1].x - this->m_vData[idx].x;
    }
    float calc_v(unsigned int idx)
    {
        assert(idx > 0);
        return (6.0f * ((this->m_vData[idx + 1].y - this->m_vData[idx].y) / calc_h(idx)
            - (this->m_vData[idx].y - this->m_vData[idx - 1].y) / calc_h(idx - 1)));
    }

public:
    CubicSpline(
        const vector<vec2D> &vData = vector<vec2D>()
        )
    {
        // Error handling
        const unsigned int n = unsigned int(vData.size());
        assert(n >= 2);

        // Copy the data
        this->m_vData.resize(n);
        copy(vData.begin(), vData.end(), this->m_vData.begin());

        // Set data to a matrix A
        cv::Mat A = cv::Mat::zeros(n, n, CV_32F);
        cv::Mat b(n, 1, CV_32F);

        A.at<float>(0, 0) = 1.0f;
        b.at<float>(0, 0) = 0.0f;
        for (unsigned int r = 1; r < n - 1; ++r)
        {
            A.at<float>(r, r - 1) = calc_h(r - 1);
            A.at<float>(r, r + 1) = calc_h(r);
            A.at<float>(r, r) = 2.0f * (A.at<float>(r, r - 1) + A.at<float>(r, r + 1));

            b.at<float>(r, 0) = this->calc_v(r);
        }
        A.at<float>(n - 1, n - 1) = 1.0f;
        b.at<float>(n - 1, 0) = 0.0f;

        // Calculate parameters c_i
        cv::Mat u = A.inv(cv::DECOMP_LU) * b;

        // Calculate coefficients
        this->m_vCoeff.resize(n - 1);
        for (unsigned int idx = 0; idx < this->m_vCoeff.size(); ++idx)
        {
            this->m_vCoeff[idx].b = u.at<float>(idx, 0) / 2.0f;
            this->m_vCoeff[idx].a =
                (u.at<float>(idx + 1, 0) - u.at<float>(idx))
                / (6.0f * (this->m_vData[idx + 1].x - this->m_vData[idx].x));
            this->m_vCoeff[idx].d = this->m_vData[idx].y;
            this->m_vCoeff[idx].c =
                (this->m_vData[idx + 1].y - this->m_vData[idx].y) / (this->m_vData[idx + 1].x - this->m_vData[idx].x)
                - (this->m_vData[idx + 1].x - this->m_vData[idx].x) * (2.0f * u.at<float>(idx, 0) + u.at<float>(idx + 1, 0)) / 6.0f;
        }
    }
    ~CubicSpline()
    {
    }

    vec2D getValue(
        unsigned int idx,
        float ratio)
    {
        // Error handling
        assert(ratio >= 0.0f && ratio <= 1.0f);
        assert(idx + 1 < this->m_vData.size());

        vec2D v2;
        v2.x = ratio * this->m_vData[idx + 1].x + (1.0f - ratio) * this->m_vData[idx].x;
        float tmp = v2.x - this->m_vData[idx].x;
        v2.y = this->m_vCoeff[idx].a * tmp * tmp * tmp
            + this->m_vCoeff[idx].b * tmp * tmp
            + this->m_vCoeff[idx].c * tmp
            + this->m_vCoeff[idx].d;

        return v2;
    }
};

Catmull-Romスプライン曲線 (Catmull-Rom Spline Curve)

2次元座標(図中,赤色)を与えると,その間を,Catmull-Romスプライン補間により補間した点(図中,緑色)を描画するプログラム.数値計算用及び描画用ライブラリとしてOpenCV 2.Xを利用しています.3次スプライン曲線程なめらかではないが,描画に対する制約は少ないです.

f:id:Mugichoko:20180405185123g:plain

main.cpp

//
// 2015.8.19
// Mugichoko
//
#include <iostream>
using namespace std;

#include "Catmull-RomSpline.h"

const int wndSize = 512;

int main()
{
    vector<cv::Point2f> vData;
    vData.push_back(cv::Point2f(0.0f, 0.0f));
    vData.push_back(cv::Point2f(100.0f, 200.0f));
    vData.push_back(cv::Point2f(130.0f, 50.0f));
    vData.push_back(cv::Point2f(255.0f, 255.0f));
    CatmullRomSpline<cv::Point2f> myCRSpline(vData);

    cv::Mat img = cv::Mat::zeros(wndSize, wndSize, CV_8UC3);

    cv::Point2f shift(wndSize / 4, wndSize / 4);
    // Draw control points
    for (int idx = 0; idx < vData.size(); ++idx)
        cv::circle(img, vData[idx] + shift, 5, cv::Scalar(0, 0, 255), -1);
    cv::imshow("VIZ", img);
    cv::waitKey();
    // Draw interpolated points
    for (int i = 0; i < vData.size() - 1; ++i) {
        for (int j = 0; j < 10; ++j) {
            cv::Point2f v = myCRSpline.getValue(i, j / 10.0f);
            cv::circle(img, v + shift, 2, cv::Scalar(0, 255, 0), -1);

            cv::imshow("VIZ", img);
            cv::waitKey(100);
        }
    }
    cv::circle(img, myCRSpline.getValue((unsigned int)vData.size() - 2, 1.0f) + shift, 2, cv::Scalar(0, 255, 0), -1);
    cv::imshow("VIZ", img);
    cv::waitKey(2000);

    return 0;
}

Catmull-RomSpline.h

//
// 2015.8.19
// Mugichoko
//
// Ref.: http://markun.cs.shinshu-u.ac.jp/learn/cg/cg3/index5.html
//

#ifndef CATMULLROMSPLINE_H
#define CATMULLROMSPLINE_H

#include <iostream>
using namespace std;
#include <opencv2/opencv.hpp>

static const float power = 0.5f;

template<typename vec2D>
class CatmullRomSpline
{
private:
    vector<vec2D> m_vData;

    float calcVal(float x0, float x1, float v0, float v1, float t)
    {
        return (2.0f * x0 - 2.0f * x1 + v0 + v1) * t * t * t
            + (-3.0f * x0 + 3.0f * x1 - 2.0f * v0 - v1) * t * t
            + v0 * t + x0;
    }

public:
    CatmullRomSpline(
        const vector<vec2D> &vData = vector<vec2D>()
        )
    {
        // Error handling
        const unsigned int n = unsigned int(vData.size());
        assert(n >= 2);

        // Copy the data
        this->m_vData.resize(n + 2);
        copy(vData.begin(), vData.end(), this->m_vData.begin() + 1);
        this->m_vData[0] = vData[0];
        this->m_vData[this->m_vData.size() - 1] = vData[vData.size() - 1];
    }
    ~CatmullRomSpline()
    {
    }

    vec2D getValue(
        unsigned int idx,
        float t)
    {
        // Error handling
        assert(t >= 0.0f && t <= 1.0f);
        assert(idx < this->m_vData.size());

        vec2D pos;
        vec2D &p1 = this->m_vData[idx];
        vec2D &p2 = this->m_vData[idx + 1];
        vec2D &p3 = this->m_vData[idx + 2];
        vec2D &p4 = this->m_vData[idx + 3];
        vec2D v0 = (p3 - p1) * power;   // (p3 - p1) * power (Usually power is 0.5f)
        vec2D v1 = (p4 - p2) * power;   // (p4 - p2) * power (Usually power is 0.5f)

        pos.x = this->calcVal(p2.x, p3.x, v0.x, v1.x, t);
        pos.y = this->calcVal(p2.y, p3.y, v0.y, v1.y, t);

        return pos;
    }
};

#endif