Bezier曲線をC++, Python, Juliaで書いてみた

目次

はじめに

ベジエ曲線を題材にして、C++、Python、Juliaでそれぞれのコードの書き方の違いをメモしておくための記事です。自分用のメモなのでコードの書き方などおかしなところがあるかもしないので、その時は指摘していただけると幸いです。

ベジエ曲線とは

ベジェ曲線（ベジェきょくせん、Bézier Curve）またはベジエ曲線とは、N 個の制御点から得られる N − 1 次曲線である。フランスの自動車メーカー、シトロエン社のド・カステリョ（英語版） とルノー社のピエール・ベジェにより独立に考案された。ド・カステリョの方が先んじていたが、その論文が公知とならなかったためベジェの名が冠されている[1]。コンピューター上で滑らかな曲線を描くのに2次ベジェ曲線 (Quadratic Bézier curve) や3次ベジェ曲線 (Cubic Bézier curve) などが広く利用されている。
出典：「ベジェ曲線 - Wikipedia」

C++

// "Bezier_curves.cpp"
// @Fumihachi

#include <iostream>
#include <random>

#include <opencv2/core.hpp>
#include <opencv2/highgui.hpp>
#include <opencv2/line_descriptor.hpp>

// OpecCVのライブラリの読み込み
#ifdef _DEBUG
#pragma comment(lib,"opencv_core330d.lib")
#pragma comment(lib,"opencv_highgui330d.lib")
#pragma comment(lib,"opencv_line_descriptor330d.lib")
#pragma comment(lib,"opencv_imgproc330d.lib")
#else
#pragma comment(lib,"opencv_core330.lib")
#pragma comment(lib,"opencv_highgui330.lib")
#pragma comment(lib,"opencv_line_descriptor330.lib")
#pragma comment(lib,"opencv_imgproc330.lib")
#endif

#define N_CTRL 4  // 制御点の個数
#define RANDAM 0  // 1->制御点をランダムに生成

// 二項係数nCrの計算
int comb(int n, int r)
{
    int num = n - r + 1, x = 1;

    if (n == r) 
        return 1;
    else if (n < r) 
        std::cout << "Error : please set n >= r" << std::endl;

    for (int i = 1; i <= r; i++)
        x = x*num++ / i;

    return x;
}

// Bernstein係数
float bernstein(int n, int i, float t)
{
    int k = 0;
    float a = 1.0f, b = 1.0f;

    for (k = 0; k < i; k++)
        a *= t;

    for (k = 0; k < n-i; k++) 
        b *= 1 - t;

    return float(comb(n, i)) * a * b;
}

// Bezier curve の算出
cv::Point2f bezier_curve(int n, float t, float *px, float *py)
{
    float x = 0.0f, y = 0.0f, B = 0.0f;

    for (int i = 0; i <= n; i++) {
        B  = bernstein(n, i, t);
        x += B * px[i];
        y += B * py[i];
    }

    return cv::Point2f(x,y);
}

int main(void)
{
    // 制御点の指定
    float px[N_CTRL] = { 0.605222, 0.697092, 2.45667, 5.75014 };
    float py[N_CTRL] = { 0.983481, 4.994620, 2.51148, 8.39767 };

    // ランダムに制御点を生成
#if RANDAM
    std::cout << "Use randam control points." << std::endl;
    std::random_device rnd;
    for (int i = 0; i < N_CTRL; ++i) {
        px[i] = 10.0f * float(rnd()) / std::random_device::max();
        py[i] = 10.0f * float(rnd()) / std::random_device::max();
        std::cout << "[x, y] = " << px[i] << ", " << py[i] << std::endl;
    }
#else
    std::cout << "Use default control points." << std::endl;
    for (int i = 0; i < N_CTRL; ++i)
        std::cout << "[x, y] = " << px[i] << ", " << py[i] << std::endl;
#endif // RANDAM


    // Bezier曲線
    std::vector<cv::Point2f> p;
    for (float t = 0; t <= 1.0; t += 0.01)
        p.push_back(bezier_curve(3, t, px, py));

    // openCVでグラフを描画
    float scale = 50.0;
    float margin = 10.0;

    cv::Mat img(cv::Size(500, 500), CV_8UC3, cv::Scalar(255, 255, 255));

    for (int i = 0; i < 4; i++) 
        cv::circle(img, 
            cv::Point2f(
                px[i] * scale + margin, 
                img.size().height-py[i] * scale - margin), 
            10, cv::Scalar(0, 0, 255));

    for (int i = 0; i < p.size()-1; i++) {
        cv::Point2f p0 = cv::Point2f(
            p[i].x*scale + margin, 
            img.size().height - p[i].y*scale - margin);
        cv::Point2f p1 = cv::Point2f(
            p[i+1].x*scale + margin, 
            img.size().height - p[i+1].y*scale - margin);
        cv::line(img, p0, p1, cv::Scalar(255, 0, 0));
    }

    cv::imshow("image", img);
    cv::waitKey(-1);
}

Bezier curve (cpp)

Python

# "Bezier_curves.py"
# @Fumihachi

import matplotlib.pyplot as plt
import numpy as np
from numpy.random import *
from scipy.special import comb


def bernstein(n, i, t):
    return comb(n, i) * t ** i * (1 - t) ** (n - i)


def bezier_curve(n, t, p):
    return np.dot([bernstein(n, i, t) for i in range(n + 1)], p)


# 制御点をランダムに生成
N = 4  # Number of control points
L = 10  # Max range of X, Y
p = rand(N, 2) * L  # Create control points
p = p[p[:, 0].argsort(), :]  # Sort to x-axis
print(p)

# Calculate
x = np.array([bezier_curve(N - 1, t, p) for t in np.linspace(0, 1, 100)])

# Draw
plt.plot(p[:, 0], p[:, 1], marker='o')
plt.plot(x[:, 0], x[:, 1])
plt.title('Bezier curve')
plt.xlabel('X')
plt.ylabel('Y')
plt.grid('on')
plt.show()

Bezier curve (python)

Julia

# "Bezier_curves.jl"
# @Fumihachi

using PyPlot
using LinearAlgebra

function bernstein(n, i, t)
    return binomial(n, i) * t^i * (1-t)^(n-i)
end

function bezier_curve(n, t, p)
    x = [dot([bernstein(n, i, t) for i = 0:n], p[:,1]) dot([bernstein(n, i, t) for i = 0:n], p[:,2])];
    return x
end

# Create control points at random
N = 4;
L = 10;
p = rand(N, 2) * L;
p = p[sortperm(p[:,1]), :]

# Calculate
nt  = 100
x   = Array{Float64,2}(undef, nt+1, 2)
for k = 1:nt+1
    x[k,:] = bezier_curve(N - 1, (k-1)/100, p)
end

# Draw
plot(p[:,1], p[:,2], marker="o")
plot(x[:,1], x[:,2])
title("Bezier curve")
xlabel("X")
ylabel("Y")
grid("on")
show()

Bezier curve (julia)