📈

ある点に最も近いベジェ曲線上の点を求める

2023/02/21に公開

JavaScript

graphics

tech

はじめに

3次ベジェ曲線を扱うことになり、ある点に最も近いベジェ曲線上の点を計算する必要があったのでJavaScriptで書いてみた。

ベジェ曲線

ベジェ曲線はWikipediaで説明されている通りで、コンピュータグラフィックスなどで使われるパラメトリック曲線のひとつ。今回は3次ベジェ曲線なのでパラメーターをtとするとベジェ曲線上の点は下記式で表すことができる。

B(t) = (1 - t)^3P_0 + 3(1 - t)^2 tP_1 + 3(1 - t)t^2 P_2 + t^3 P_3

ある点に最も近いベジェ曲線上の点

ベジェ曲線上の点の接線（ベジェ曲線の微分）と、点Pからベジェ曲線上の点へのベクトルが垂直になる点が、点Pに最も近いベジェ曲線上の点になる。

ベジェ曲線の微分もWikipediaで説明されている通り。

B'(t) = 3(1 - t)^2 (P_1 - P_0) + 6(1 - t)t(P_2 - P_1) + 3t^2 (P_3 - P_2)

垂直判定はドット積が0になればOK。

(B(t) - P) \cdot B'(t) = 0

このtを2分探索で求めてみることにする。

ベジェ曲線の分割

ベジェ曲線の分割については、ベジエ曲線の長さを求めるに書かれていたものをそのまま使っている。

$P_1, P_2, P_3, P_4$ で定義された3次ベジェ曲線を $t = 0.5$ で分割する場合はこんな感じ。

前半

分割してできる前半の新しいベジェ曲線の始点、制御点1、制御点2、終点は下記。

P_0

\frac{(P_0 + P_1)}{2}

\frac{P_0 + 2P_1 + P_2}{4}

\frac{P_0 + 3P_1 + 3P_2 + P_3}{8}

後半

分割してできる後半の新しいベジェ曲線の始点、制御点1、制御点2、終点は下記。

\frac{P_0 + 3P_1 + 3P_2 + P_3}{8}

\frac{P_1 + 2P_2 + P_3}{4}

\frac{(P_2 + P_3)}{2}

P_3

探索方法

ベジェ曲線を分割しながら、適当なt（とりあえず0.5）で求めたベジェ曲線上の点の接線とある点との直線が垂直か（ドット積が0か）調べ、ドット積の値から次の分割位置を決めて再帰、ベジェ曲線の長さが小さくなればそのベジェ曲線区間の $t=0.5$ を探索点とする方法にしてみた。

これだとこんなベジェ曲線だとうまく動作しなかった。

別の方法として、 $t=0.5$ でベジェ曲線を分割して、ベジェ曲線を内包する矩形の各辺とある点との距離を求めて、分割したベジェ曲線（を内包する矩形）のどちらがある点に近いかを調べ、ベジェ曲線の長さが短くなるまで再帰する方法にしてみた。

直線と点との距離

矩形の各辺とある点との距離は点と線分の距離を求めるを参考に下記で求めることにした。

d2(\overline{t}) = \frac{(a(y1-y0)-b(x1-x0))^2}{a^2+b^2}

ソースコード

ユーティリティ

ベクトル計算が楽になるようユーティリティ関数を用意する。addは和、subは差、mulは要素をv倍、divは要素を1/v、dotはドット積、norm_squaredはベクトルの長さの二乗。

const add = (a, b) => [a[0] + b[0], a[1] + b[1]];

const sub = (a, b) => [a[0] - b[0], a[1] - b[1]];

const mul = (a, v) => a.map(i => i * v);

const div = (a, v) => a.map(i => i / v);

const dot = (a, b) => a[0] * b[0] + a[1] * b[1];

const norm_squared = a => dot(a, a);

ベジェ曲線関連

曲線上の点

曲線上の点の数式をプログラムにしたもの。引数のtはパラメータ、bzは4要素配列で定義した3次ベジェ曲線。

const completion = (t, bz) => {
  const t_ = 1 - t;
  return add(add(add(mul(bz[0], t_ * t_ * t_),
                     mul(bz[1], 3 * t_ * t_ * t)),
                     mul(bz[2], 3 * t_ * t * t)),
                     mul(bz[3], t * t * t));
};

接線

微分の数式をプログラムにしたもの。

const differential = (t, bz) => {
  const t_ = 1 - t;
  return add(add(mul(sub(bz[1], bz[0]), 3 * t_ * t_),
                 mul(sub(bz[2], bz[1]), 6 * t_ * t)),
                 mul(sub(bz[3], bz[2]), 3 * t * t));
};

分割

分割の数式をプログラムにしたもの。 $t=0.5$ で分割して分割後のベジェ曲線を配列で返す。

const split_bezier = bz => {
  const center = completion(0.5, bz);
  return [
    [
      bz[0],
      div(add(bz[0], bz[1]), 2),
      div(add(add(bz[0], mul(bz[1], 2)), bz[2]), 4),
      center
    ],
    [
      center,
      div(add(add(bz[1], mul(bz[2], 2)), bz[3]), 4),
      div(add(bz[2], bz[3]), 2),
      bz[3]
    ]
  ];
};

直線と点との距離（の二乗）

const diff2_to_line = (line, p) => {
  const ps = sub(line[0], p);
  const [a, b] = sub(line[1], line[0]);
  const n2 = norm_squared([a, b]);
  const tt = -(a * ps[0] + b * ps[1]);
  if (tt < 0)
    return norm_squared(ps);
  else if (tt > n2)
    return norm_squared(sub(line[1], p));
  const f1 = a * ps[1] - b * ps[0];
  return f1 * f1 / n2;
};

ベジェ曲線を内包する矩形と点との距離（の二乗）

3次ベジェ曲線を構成する4点で作られる矩形の各辺について長さを求めそれらの最短を求めている。

const diff2_to_polygon = (bz, p) => {
  return Math.min(
    diff2_to_line([bz[0], bz[1]], p),
    diff2_to_line([bz[1], bz[2]], p),
    diff2_to_line([bz[2], bz[3]], p),
    diff2_to_line([bz[3], bz[0]], p)
  );
};

探索処理

bzはベジェ曲線、pはある点、t0とt1は探索区間を表すパラメータの値。

分割した2つのベジェ曲線について距離が短い方について再帰しながらさらに分割し、ベジェ曲線の長さが短くなると再帰を終える。

const done_or_recursive = (bz, p, t0, t1) => {
  const n2 = norm_squared(sub(bz[3], bz[0]));
  if (n2 < 1 * 1)
    return (t0 + t1) * 0.5;
  return neighbor_bezier(bz, p, t0, t1);
};

const neighbor_bezier = (bz, p, t0, t1) => {
  const splitbz = split_bezier(bz);
  const d0 = diff2_to_polygon(splitbz[0], p);
  const d1 = diff2_to_polygon(splitbz[1], p);
  const tcenter = (t0 + t1) * 0.5;
  return d0 < d1
    ? done_or_recursive(splitbz[0], p, t0, tcenter)
    : done_or_recursive(splitbz[1], p, tcenter, t1);
};

動かしてみる

こんなベジェ曲線で動かしてみた。

const bezier = [
  [ 50,  50],
  [ 50, 300],
  [550, 300],
  [550, 550]
];

CanvasでY軸が逆だけど図にするとこんな感じ。

上に挙げたうまく動作しなかったベジェ曲線もOK。

const bezier = [
  [150, 550],
  [ 50,  50],
  [550,  50],
  [250, 550]
];

リポジトリ

ソースコード一式はこちら。

デモ

Canvas使ってブラウザ上で動作確認できるようにしてみた。