Logo address

補間曲線 (Interpolation Curves)

目次

ここでは Bézier 曲線を使った補間曲線を採り上げる。

SVG Path

SVG の Path は守備範囲が広い。万屋である。
折れ線、ベジェ曲線、閉じた折れ線、閉じたベジェ曲線が描ける。

SVG の Path は、3次のベジェ曲線までを描く事ができる。3次のベジェ曲線は一般に 4 点を与えるが、特別の場合として、3点だけを与える曲線もある。

次の図は、3点を与えた時の、3種のSVG Path の振る舞いを示す。

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。
図1

このプログラムのコードは

#encoding: utf-8
import svg
c = svg.Canvas(0,0,600,500,map=(-1.5,1.5,1.5,-1.5))
c.grid(1,1,stroke="cyan")
c.text(1,1,text="(1,1)",anchor="NW")
c.text(-1,1,text="(-1,1)",anchor="NE")
c.text(0,-1,text="(0,-1)",anchor="NW")
c.path("M -1 1 L 0 -1 1 1", stroke="black;dasharray:4 4")
c.path("M -1 1 Q 0 -1 1 1")
c.path("M -1 1 S 0 -1 1 1",stroke="red")
c.path("M -1 1 C 0 -1 0 -1 1 1",stroke="blue")
c.close()

このライブラリの特徴は、座標を論理座標で描けることである。

ここに現れる path 関数は、この SVG Library の中では low level の関数で、SVG の仕様がむき出しになっている。(論理座標への変換は行っている)

ここに示した図は L,Q,S,C の意味を捉えるには良いであろう。このうち、黒の線で示した Q は放物線のセグメントである。

ところで、文字列を書く位置の指定は、このプログラムのように、Python/Tk 流の方が分かりやすくて僕の好みである。

3次ベジェ曲線 (Cubic Bézier Curve)

2012/12/02
2012/12/04 改訂

定義

t を 0 から t の実数とする。4点 P0、P1、P2、P3 を与え、次の R(t) で定まる曲線が 3 次のベジェ曲線である。
R(t) = (1-t)3 P0 + 3(1-t)2 t P1 + 3(1-t)t2P2 + t3P3

R(t) を t で微分すると
dR(t)/dt = 3t2(P3-P2)+6(1-t)t(P2-P1)+3(1-t)2(P1-P0)
これを簡単に DR(t) とすると
DR(0) = 3(P1-P0)
DR(1) = 3(P3-P2)
この関係は、Catmull-Rom spline curve で使われる。

例1

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。
図2. 3次ベジェ曲線

#encoding: utf-8
import svg
from table import *
c = svg.Canvas(0,0,500,500,map=(-1,3,3,-1))
c.grid(1,1,stroke="cyan",range=(-0.5,2.5,2.5,-0.5),label=("x","y"),pad=0.2)

mblack = c.symbol("circle",fill="black")
mgreen = c.symbol("circle",fill="green")
mred = c.symbol("circle",fill="red")
mblue = c.symbol("circle",fill="blue")

P = [(0, 0), (1, 2), (2, 1.5), (2, 0)]
lab=("P[0]","P[1]","P[2]","P[3]")
c.path("M %s %s C %s %s %s %s %s %s"%tuple(svg.expand(P)),stroke="width:2")
anchor = ("NE","SE","W","NE")
c.line(P,stroke="green")
c.mark(P,ref=mblack,label=lab,anchor=anchor,pad=0.5)

t = 0.4
for n in range(0,len(P)):
  P[n] = T(*P[n])
A = (1-t)*P[0]+t*P[1]
B = (1-t)*P[1]+t*P[2]
C = (1-t)*P[2]+t*P[3]
M = (1-t)*A+t*B
N = (1-t)*B+t*C
Q = (1-t)*M+t*N

L = [A.value(),B.value(),C.value()]
c.line(L)
c.mark(L,ref=mgreen,label=("A","B","C"),anchor=("SE","SW","W"),pad=0.2)

L = [M.value(),N.value()]
c.line(L)
c.mark(L,ref=mred,label=("M","N"),anchor="S",dy="-0.5")
c.mark(Q.value(),ref=mblue,label="Q",anchor="NW")

c.close()

今回は、新たに次の機能をサポートした。

grid() の label は、細かなチューニング(文字フォント、文字の大きさ、位置など)をサボートしていない。この部分は好みに依存する部分が大きく、まじめにやる価値は殆どない。必要とあれば、直接 text() method で書けばよいし、その事は面倒ではない。
T class は、table データのクラスであり、データの一次元配列、2次元配列をサポートする。(内部構造は次元依存性を持たないが動作の確認が必要である。) 2項演算は、(vector や matrix と異なり)要素毎の演算として定義されている。
mark() は、指定された座標に symbol と(必要なら) label を付ける。label を付けた場合には、anchordxdy でラベルの位置を指定できる。座標は複数指定でき、labelahchordxdy はリストあるいはタプルで与える事ができる。

T は本当に必要なのかと考えたりもする。新たなクラスを作った場合には、その値を取り出すのに

	A.value()
のように、余計なものを必要とする。(もっとも line()T class を認識していれば、この問題は line() のコードの中で解決するが、今のところ外付けである。)

逆ポーランドの演算を使えば、T class は要らない。その場合、

	for n in range(0,len(P)):
	  P[n] = T(*P[n])
	A = (1-t)*P[0]+t*P[1]
	...
	L = [A.value(),...]
は単に
	A = cal((1-t),P[9],"*",t,P[1],"*","+")
	L = [A,...]
と簡潔に書くことになる。(cal() は既に table module でサポートとされている)
でも世の中、逆ポーランドの演算を知らない人が多いからね...

例2

線分P[0]P[3]が線分P[1]P[2]に平行な場合には、幾何学的な考察を行いやすい。

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。

さて、この場合には、t=0.5 で次の図になる。

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。

このことから、この曲線の最大の高さは、台形P[0],P[1],P[2],P[3]の高さの 3/4 であることが解る。

例3

t=0.4 とし、path を

	P = [(0, 0), (1, 2), (1, 2), (2, 0)]
で与えた場合、

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。

この場合には3つの点、BP[1]P[2] が縮退する。

t=0.5 では

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。

従って頂点の高さは、Bの高さの 3/4 である。

例4

次に、t=0.4

	P = [(0, 0), (1, 2), (1, -2), (2, 0)]
とした場合

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。

例5

t=0.4 で、P[2] = P[3] のケース。つまり、

	P = [(0, 0), (1, 2), (2, 0), (2, 0)]

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。

生成された曲線は、左右対称のように見えるけど、本当か?

放物線による補間

ここでは y = f(x) で与えられる関数の補間曲線を放物線で描く事を考える。

考え方

次の図を考えよう。

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。

この曲線は

	c.path("M -1 1 Q 0 -1 1 1")
で描いたもので、放物線になっている。A,B,C の座標は各々
	A=(-1,1), B=(1,1), C=(0,-1)
で、この3点を与えて、2次のベジェ曲線を描くと図のように N=(0,0) を通る放物線が生成される。
MAB の中点で、N は、MC の中点になっている。
こうした関係は、アフィン変換注1によって変化しない。

注1. 平行移動や、回転、さらに、X方向、Y方向に伸縮させる変換、skew (ねじれ) による変換を言う。

従って放物線上の3点

	(x1,f(x1))
	(x2,f(x2))
	(x0,f(x0))	# x0 = (x1+x2)/2
を与えると、次の図のように2次のベジェ曲線によって放物線を描く事が可能であることが解る。

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。

上の2つの図を描いたプログラムのコードは各々次のようになっている。

#encoding: utf-8
import svg
c = svg.Canvas(0,0,600,500,map=(-1.5,1.5,1.5,-1.5))
c.grid(1,1,stroke="cyan")
c.text(1,1,text="(1,1)",anchor="NW")
c.text(-1,1,text="(-1,1)",anchor="NE")
c.text(0,-1,text="(0,-1)",anchor="NW")
c.path("M -1 1 L 0 -1 1 1", stroke="black;dasharray:4 4")
c.path("M -1 1 Q 0 -1 1 1")
c.line(-1,1,1,1,stroke="red")
c.line(0,-1,0,1,stroke="red")
c.text(-1,1,text="A",anchor="SW")
c.text(1,1,text="B",anchor="SW")
c.text(0,-1,text="C",anchor="SW")
c.text(0,1,text="M",anchor="SW")
c.text(0,0,text="N",anchor="SW")
c.close()

#encoding: utf-8
import svg
import math
co = (2.0,-3.0,1.0)

def f(x):
  (a,b,c) = co
  return a*x*x + b*x + c

c = svg.Canvas(0,0,500,600,map=(-1.5,3.5,3.5,-2.5))
c.grid(1,1,stroke="cyan")
c.text(0,0,text="(0,0)",anchor="NE")

x1 = 0.0; y1 = f(x1)
x2 = 2.0; y2 = f(x2)
x0 = (x1+x2)/2
y0 = 2*f(x0) - (f(x1)+f(x2))/2
c.path("M %s %s L %s %s %s %s"%(x1,y1,x0,y0,x2,y2), stroke="black;dasharray:4 4")
c.path("M %s %s Q %s %s %s %s"%(x1,y1,x0,y0,x2,y2))
c.line(x1,y1,x2,y2,stroke="red")
c.line(x0,y0,x0,(y1+y2)/2,stroke="red")
c.text(x1,y1,text="A",anchor="NW")
c.text(x2,y2,text="B",anchor="NW")
c.text(x0,y0,text="C",anchor="NW")
c.text(x0,(y1+y2)/2,text="M",anchor="NW")
c.text(x0,f(x0),text="N",anchor="NW")
c.close()

放物線による関数の近似曲線

ここで述べたアイデアは、放物線による関数の近似曲線に応用できる。

次の図は3次曲線

	f(x) = x*(x-1)*(x-3)
を描いている。

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。
図3: 関数のプロット

これを生成したプログラムは

#encoding: utf-8
import svg

def f(x):
  return x*(x-1)*(x-3)

def draw(c,x1,x2):
  y1,y2 = f(x1),f(x2)
  x0 = (x1+x2)/2
  y0 = 2*f(x0) - (f(x1)+f(x2))/2
  c.path("M %s %s Q %s %s %s %s"%(x1,y1,x0,y0,x2,y2))

c = svg.Canvas(0,0,500,600,map=(-1.5,3.5,3.5,-2.5))
c.grid(1,1,stroke="cyan")
c.text(0,0,text="(0,0)",anchor="NE")

a = 0.5
for n in range(-1,7,1):
  draw(c,a*n,a*(n+1))

c.close()

区間をもっと細かく採れば、近似精度は上がるが、生成される SVG のファイルは大きくなる。
目的に応じた妥協が必要だ。

func() method

2012/12/01

結局次のように func() method を SVG Library に追加することとした。
これを使うと、図3と同じ結果は次のプログラムで得られる。(区間の分割数は増やしてある)

import svg
def f(x):
  return x*(x-1)*(x-3)
c = svg.Canvas(0,0,500,600,map=(-1.5,3.5,3.5,-2.5))
c.grid(1,1,stroke="cyan")
c.text(0,0,text="(0,0)",anchor="NE")
c.func(f,range=(-1,3.5),n=50)
c.close()

n=50 は区間の分割数である。(default 100)

Plot class を追加しようと思ったが、メリットがないので止めた。

Smooth lines

点の集合を与えて、SVG でスムースな曲線を描くにはどうするか?

ここでは、図1の、放物線と3次ベジェを使って描いてみた。

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。
図4スムーズ曲線

このプログラムのコードは

#encoding: utf-8
import svg
c = svg.Canvas(0,0,600,500,map=(-1,6,5,-1.5))
c.grid(1,1,stroke="cyan")
c.text(0,0,text="(0,0)", anchor="NE")
c.text(0,1,text="(0,1)",anchor="SE")
c.text(1,0,text="(1,0)",anchor="N")
p = c.symbol("circle",fill="black")

t = [0, 0, 1, 1, 2, 0, 3, 2, 3, 4, 1, 4, 0, 2]
for n in range(0,len(t),2):
  c.use(t[n],t[n+1],ref=p)

c.line(*t, smooth=0,stroke="black")
c.line(*t, smooth=1,stroke="red")
c.line(*t, smooth=2,stroke="green")
c.close()

譜1

スムーズな曲線の描画アルゴリズムは Python/Tk 流にしてある。

端点は通るが、それ以外の点は通らない事に注意せよ。
このアルゴリズムでは曲線は、(端点を除いて)2つの点の中点を通り、そこで接する。

与えられた点を通る曲線を描きたい場合にはどうするか?
どうやら、そのような曲線を描くオプションを SVG は備えていないようである。自分で作らなくてはならない。それが結構くせ者である。

反省

2012/12/08

3次のベジェ曲線に関しては、もっと柔軟な考え方がある。

次の図は、図1の上に、いくつかの3次のベジェ曲線を追加している。

ベジェ曲線を、P0,P1,P2,P3 で表すとする。P0 = (-1,1)P3 = (1,1)Q = (0,-1) とする。P1P2は、線分P0P3に平行で、各々、直線P0,Q、直線P3,Qとの交点とする。
すると

	P1 = (1-a)*P0 + a*Q
	P3 = (1-a)*P3 + a*Q
である。図では、a
	0.25, 0.5, 2.0/3, 0.75, 1, 1.25, 1.5
の場合について描かれている。

a = 2.0/3 で、このベジェ曲線は、黒い線で描かれた2次の曲線と完全に重なる。この事はベジェ曲線をC(t):

	C(t) = (1-t)**3 * P0 + 3*(1-t)**2 *t * P1 + 3*(1-t)*t**3 * P2 + t**3 * P3
として
	P1 = (1-a)*P0 + a*Q
	P2 = (1-a)*P3 + a*Q
と置くと、C(t)tについて3次の項が a = 2.0/3 で消える事から確認できる。

つまり、譜1の line() method の中の smooth の仕様は、変更した方が良いのだ。smooth の値は、a から決めるべきである。

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。
図5: a を変化させた時の3次のベジェ曲線

#encoding: utf-8
import svg
from table import *
c = svg.Canvas(0,0,600,500,map=(-1.5,1.5,1.5,-1.5))
c.grid(1,1,stroke="cyan")
c.text(1,1,text="(1,1)",anchor="NW")
c.text(-1,1,text="(-1,1)",anchor="NE")
c.text(0,-1,text="(0,-1)",anchor="NW")
c.path("M -1 1 L 0 -1 1 1", stroke="black;dasharray:4 4")
c.path("M -1 1 Q 0 -1 1 1")
c.path("M -1 1 S 0 -1 1 1",stroke="red")
P0 = (-1,1)
P3 = (1,1)
Q = (0,-1)
# P1 = (1-a)*P0 + a*Q
for a in 0.25,0.5,2.0/3,0.75,1,1.25,1.5:
  P1 = cal(1-a,P0,"*",a,Q,"*","+")
  P2 = cal(1-a,P3,"*",a,Q,"*","+")
  c.path("M -1 1 C %s %s %s %s 1 1"%expand(P1,P2),stroke="blue")
c.close()

B-Spline Curves

Uniform B-Spline Curve

2012/12/20

コントロールポイントを P[0], P[1], ..., P[n]
目標: 3次の spline 曲線を繋げて、全体として、C2 連続(2次微分が連続)な曲線を作る。

#encoding: utf-8
#
#	B-Spline
#
import svg
from table import *

def N(i,k,t):
  # Uniform Cox-deBoor recursion formula
  # k: order (degree k-1)
  if k==1:
    if i <= t and t < i+1:
      return 1
    else:
      return 0
  return ((t-i)*N(i,k-1,t)+(i+k-t)*N(i+1,k-1,t))/(k-1)

def q(k,t,P):
  # k: order of B-spline (degree k-1)
  # t: 0 <= t < len(P)
  # P: controle points
  t = float(t)
  sum = (0,0)
  for i in range(0,len(P)):
    # sum = sum + P[i]*N(i,k,t)
    sum = cal(sum,P[i],N(i,k,t),"*","+")
  return sum

c = svg.Canvas(0,0,600,500,map=(-1,6,5,-1.5))
sym = c.symbol("circle",fill="red")
sym1 = c.symbol("cross",stroke="green")
sym2 = c.symbol("circle",fill="black")
c.grid(1,1,stroke="cyan")
c.text(0,0,text="O", anchor="NE",pad=0.3)
#P = [(0,0),(0,0), (0, 0), (1, 1), (2, 0), (3, 2), (3, 4), (1, 4), (0, 2),(0,2)]
P = [(0, 0), (1, 1), (2, 0), (3, 2), (3, 4), (1, 4), (0, 2)]
m = len(P)

# Cox-deBoor for uniform case
U = []
for n in range(0,100):
  U.append(q(4,m*n/100.0,P))
c.mark(U,ref=sym1)

# Alexandre Hardy
Q = [None]*(m-2)
DQ = [None]*(m-2)
for n in range(0,m -2):
  Q[n] = cal(1.0/6,P[n],"*",2.0/3,P[n+1],"*","+", 1.0/6,P[n+2],"*","+")
  DQ[n] = cal(-0.5,P[n],"*",0.5,P[n+2],"*","+")
c.cplot2(Q,DQ,stroke="width:3")

c.mark(P,ref=sym)
c.mark(Q,ref=sym2)
c.close()

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。
図6: B-spline curve

図では B-spline のコントロールポイントを赤丸で示してある。
緑の×印は、Cox-deBoor の再帰式から求めた曲線上の点。
黒の太線は、Hardy の論文をヒントに得られた区間的ベジェ曲線(piece wise bezier curve)。
太線上の黒丸は、ノットの位置。

Hardyは、ドクター論文(文献[3])の中で、B-spline を、解りやすく解説している。彼によると、B-spline の basis 関数は

	B0(t) = (1-t)**3 /6
	B1(t) = (3*t**2 -6*t**2 + 4)/6
	B2(t) = (-3t**3 + 3*t**2 + 3*t + 1)/6
	B3(t) = t**3 /6
で、この事から容易に、ノット(knot)の位置と、接線ベクトルが求められ、そこから、隣り合うノットを結ぶベジェ曲線を描く事が可能である。

プログラムのコントロールポイント P の端に、重複した点を加える、例えば、

	P = [(0,0),(0,0), (0, 0), (1, 1), (2, 0), (3, 2), (3, 4), (1, 4), (0, 2),(0,2)]
のようにすると、得られる曲線は次のようになる。

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。
図7: B-spline curve

つまり、P[0]P[n] に向かう区間的ベジェ曲線を描く事が可能である。

Geometry of Uniform B-Spline Curve

次の図8に Uniform B-spline のコントロールポイント(P[0],...,P[n]) から曲線がどのように定まるかを示す。

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。
図8: Uniform B-spline curve の幾何学

図では、(簡単のために)P[0],...,P[6]P0,...,P6としている。
ノット位置の曲線状の点と、その接線ベクトルが分かれば、ベシェ曲線が定まる。
MP[3]P[5]の中点で、B-spline 曲線はMP4を2対1に内分する点を通る。この点の接線ベクトルは、ベクトル (P5,P3) と平行で、その1/2の大きさである。

Open Uniform B-Spline Curve

2013/01/05

4個のコントロールポイント (P0,P1,P2,P3) から、1つの3次 spline 曲線が得られる。もちろん、その時のspline曲線はベジェ曲線に一致し、(P0,P1,P2,P3)がそのベジェ曲線のコントロールポイントとなる。(従って P0P3 を通る。)
5個のコントロールポイント (P0,P1,P2,P3,P4) から、2つのベジェ曲線が得られる。もちろん、その時のベジェ曲線は P0P3 を通る。

文献

[1] David F.Rogers, J.Alen Adams, (河合慧監訳)
「コンピュータグラフィックス」(日刊工業新聞社,1999)

Catmull-Rom/Overhauser Spline Curve

2012/11/30
2012/12/09 改訂

補題

2012/12/09

t を実数とし、関数R(t)とその微分DR(t)の値が t=0t=1 で求まっているとせよ。以下

	P0 = R(0)
	P3 = R(1)
	DP0 = DR(0)
	DP1 = DR(1)
とする。その下で、R(t) を近似するベジェ曲線は、(P0,P1,P2,P3)で与えられる。ここに
	P1 = P0 + DP0/3
	P2 = P3 - DP3/3
である。

これは、ベジェ曲線の基本的な性質から得られる。

この性質を使って渦巻き曲線を補間してみる。

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。
図9: Catmull-Rom spline curve の構成法

プログラムを次ぎに示す。

#encoding: utf-8
from math import exp,sin,cos
import svg
from table import *
# === logarithmic spiral curve
r0 = 5
a = 0.2
b = 1

def f(t):
  r = r0*exp(-a*t)
  x = r*cos(b*t)
  y = r*sin(b*t)
  return x,y

def df(t):
  r = r0*exp(-a*t)
  dr = -a*r
  co = cos(b*t)
  si = sin(b*t)
  dx = dr*co - r*b*si
  dy = dr*si + r*b*co
  return dx,dy

def cplot2(c,P,DP):
  path = "M %s %s"%tuple(P[0])
  for n in range(0,len(P)-1):
    P0,Q0 = P[n],cal(P[n],DP[n],3,"/","+")
    P1,Q1 = P[n+1],cal(P[n+1],DP[n+1],-3,"/","+")
    path = path + " C %s %s %s %s %s %s"%(expand(Q0,Q1,P1))
  c.path(path)

m = 40
P,DP = [],[]
for n in range(0,m):
  P.append(f(n))
  DP.append(df(n))

c = svg.Canvas(0,0,600,600,map=(-6,6,6,-6))
c.grid(1,1)
sym1 = c.symbol("circle",fill="red")
c.mark(P,ref=sym1)
cplot2(c,P,DP)
c.close()

Catmull-Rom Spline Curve

与えられた点を通る曲線のアルゴリズムとして、「Catmull-Rom spline Curve」がある。

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。
図10: Catmull-Rom spline curve の構成法

#encoding: utf-8
import svg
from table import *
c = svg.Canvas(0,0,600,500,map=(-1,6,5,-1.5))
c.grid(1,1,stroke="cyan")
c.text(0,0,text="(0,0)", anchor="NE")
c.text(0,1,text="(0,1)",anchor="SE")
c.text(1,0,text="(1,0)",anchor="N")
sym = c.symbol("circle",fill="black")
P = [(0, 0), (1, 1), (2, 0), (3, 2), (3, 4), (1, 4), (0, 2)]
c.mark(P,ref=sym,label=("P0","P1","P2","P3","P4","P5","P6"))
P6 = P[6]
P5 = P[5]
P4 = P[4]
P3 = P[3]
V5 = cal(P4,P6,"-",2,"/")
V4 = cal(P5,P3,"-",2,"/")
a = 1./3
H5 = cal(P5,a,V5,"*","+")
H4 = cal(P4,a,V4,"*","+")
c.line(P4,P6,stroke="blue;dasharray:4 4",)
c.line(P5,H5)
c.line(P5,P3,stroke="green;dasharray:4 4")
c.line(P4,H4)
path = "M %s %s C %s %s %s %s %s %s"%expand(P4,H4,H5,P5)
c.path(path)
c.close()

譜2

この方法は、点 P[n] の接線を、隣り合う2点と平行となっていると考え、3次ベジェ曲線を次の4点

	P[n]
	P[n]+a*(P[n+1]-P[n-1])
	P[n+1]-a*(P[n+2]-P[n+1])
	P[n+1]
から決める。
ここで、a はパラメータであるが、Catmull-Rom は a=1.0/3 としている(文献[3])。

a=1./3 を選ぶ根拠は、ベジェ曲線の微分とベジェ曲線を決める4点との関係である。

aの選び方によっては、次の図のように、補間らしからぬ曲線が出来上がる。補間曲線の問題は難しく、どうやら決定的な方法は見つかっていないらしい。

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。
図11: Catmull-Rom spline curve

次のようにサポートすることとした。

#encoding: utf-8
import svg
c = svg.Canvas(0,0,600,500,map=(-1,6,5,-1.5))
c.grid(1,1,stroke="cyan")
c.text(0,0,text="(0,0)", anchor="NE")
c.text(0,1,text="(0,1)",anchor="SE")
c.text(1,0,text="(1,0)",anchor="N")
sym = c.symbol("circle",fill="black")
p = [(0, 0), (1, 1), (2, 0), (3, 2), (3, 4), (1, 4), (0, 2)]
for n in range(0,len(p)):
  c.use(p[n],ref=sym)
c.line(*p,smooth="I;looseness:0",stroke="green")
c.line(*p,smooth="I",stroke="width:2")
c.line(*p,smooth="I;looseness:1",stroke="red")
c.close()

譜3

smooth は簡単に "I" (Interpolation) とした。ついでに、cplot2()svg module の Canvas class でサポートした。

このアルゴリズムでは、点 P[n] の接線の傾きが P[n-1] と P[n+1] を結ぶ直線に平行になるように補間曲線を決める。両端の点は補助点であって、補間から外される。 しかし、曲線を spline に限定したとしても、そのような曲線は無数に存在する。

接線にどの程度引っ張られるかを決めるパラメータをここでは looseness としている。図には looseness が、0, 1./3, 1 の場合が示されている。looseness の default は 1.0/3 とした。

Overhauser Spline Curve

Overhauser spline curve は、Catmull-Rom の改良版であるが、曲線 R(t) が与えられていた時に、そのベジェ曲線による補間を決める場合に使える。詳しくは文献[3]を見よ。

Shemanarev と言う人が Overhauser と良く似たアイデアをネット上に公開している(文献[5])。こちらの方が実際的で、使いやすいかも知れない。

Overhauser の方法や Shemanarev の方法は、アフィン不変性を満たさないのではないかと思える。僕としては、アフィン不変性は譲れない。

文献

[2] James D.Foley,他「コンピュータグラフィックス -- 理論と実践」(オーム社)

[3] Alexander Hardy,"Harmonic Interpolation for Smooth Curves and Surfaces",
https://ujdigispace.uj.ac.za/bitstream/handle/10210/174/thesis.pdf

[4] E.Catmull and R.Rom,"A Class of local interpolating splines."
(Computer Aided Geometric Design,page 317-326,1974,Academic Press)

[5] Maxim Shemanarev, "Interpolation with Bezier Curves"
http://www.antigrain.com/research/bezier_interpolation/index.html

渦巻き曲線

2012/12/02

ここでは次の極座標で与えられている渦巻き曲線を SVG で描く事を考える。
r = r0exp{-λθ}
この曲線は logarithmic spiral curve と呼ばれているらしい。

この曲線を扱うのは、シンプルで、微分が簡単であり、変極点を持たないからである。以下の議論は、式が与えられており、(従って微分が計算でき)、変極点を持たない曲線に対して適用できる。全領域で変極点を持たない事を要求しない。SVG で表現しようとしている領域で変極点を持たなければ、それでも構わない。

考え方は、SVG で近似したい曲線のセグメントを n 個に分割し、各点の座標 P[n] と、その接線を求める。P[n] と P[n+1] の接線の交点を Q[n] とすると、この3点を使って2次の Bézier 曲線を描く。つまり放物線で近似する。この様子を次の図で示す。

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。

#encoding: utf-8
from math import exp,sin,cos
import svg
from table import *
# == debug for l2solution(a)
x = 2
y = 3
print 3*x+2*y, -2*x+5*y

a = [[3,2,12],[-2,5,11]]
print l2solution(a)	# OK

# === logarithmic spiral curve
r0 = 5
a0 = 0.2
b0 = 1
def f(r0,a,b,t):
  r = r0*exp(-a*t)
  x = r*cos(b*t)
  y = r*sin(b*t)
  return x,y

def df(r0,a,b,t):
  r = r0*exp(-a*t)
  dr = -a*r
  co = cos(b*t)
  si = sin(b*t)
  dx = dr*co - r*b*si
  dy = dr*si + r*b*co
  return dx,dy

c = svg.Canvas(0,0,600,600,map=(-6,6,6,-6))
mred = c.marker("circle",fill="red",scale=2)
sym = c.symbol("circle")
sym1 = c.symbol("cross")
c.grid(1,1)
P = []
for n in range(0,20):
  P = P + [f(r0,a0,b0,0.5*n)]
  c.line(f(r0,a0,b0,0.5*n),ope("+",f(r0,a0,b0,0.5*n),df(r0,a0,b0,0.5*n)),stroke="blue")
c.line(*P,marker="start:%s;mid:%s;end:%s"%(mred,mred,mred))

Q = []
m = 6
for n in range(0,m):
  a1 = df(r0,a0,b0,0.5*n)
  a2 = df(r0,a0,b0,0.5*(n+1))
  b1 = f(r0,a0,b0,0.5*n)
  b2 = f(r0,a0,b0,0.5*(n+1))
  a = [[a1[0],-a2[0],b2[0]-b1[0]],[a1[1],-a2[1],b2[1]-b1[1]]]
  t = l2solution(a)[0]
  x = a1[0]*t + b1[0]
  y = a1[1]*t + b1[1]
  Q = Q + [(x,y)]
  c.use(x,y,ref=sym)
  t = l2solution(a)[1]
  x = a2[0]*t + b2[0]
  y = a2[1]*t + b2[1]
  c.use(x,y,ref=sym1)
path = "M 5 0"
for n in range(0,m-1):
  path = path + " Q %s %s %s %s"%(Q[n][0],Q[n][1],P[n+1][0],P[n+1][1])
c.path(path)
c.close()

上の考えに基づいて、ライブラリを作る事とした。

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。
図12. 渦巻き曲線

この結果は、次のプログラムで生成される。

#encoding: utf-8
from math import exp,sin,cos
import svg
from table import *
# === logarithmic spiral curve
r0 = 5
a = 0.2
b = 1

def f(t):
  r = r0*exp(-a*t)
  x = r*cos(b*t)
  y = r*sin(b*t)
  return x,y

def df(t):
  r = r0*exp(-a*t)
  dr = -a*r
  co = cos(b*t)
  si = sin(b*t)
  dx = dr*co - r*b*si
  dy = dr*si + r*b*co
  return dx,dy

c = svg.Canvas(0,0,600,600,map=(-6,6,6,-6))
c.grid(1,1)
parafunc(c,f,df,range=(0,20),n=20)
c.close()

譜4. 渦巻き曲線

parafunc() は取りあえず、table module に放り込んだ。名前とマッチしないので、そのうちに変更されるだろう。

注: 結局廃止

補足1 (2次ベジェ曲線による補間)

2012/12/04

n 個の点と、そこでの接線ベクトルを与える方が、応用が広いのではないかと思う。
その場合、譜4は次のように書く。

#encoding: utf-8
from math import exp,sin,cos
import svg
from table import *
# === logarithmic spiral curve
r0 = 5
a = 0.2
b = 1

def f(t):
  r = r0*exp(-a*t)
  x = r*cos(b*t)
  y = r*sin(b*t)
  return x,y

def df(t):
  r = r0*exp(-a*t)
  dr = -a*r
  co = cos(b*t)
  si = sin(b*t)
  dx = dr*co - r*b*si
  dy = dr*si + r*b*co
  return dx,dy

F = []
DF = []
for n in range(0,20):
  F.append(f(n))
  DF.append(df(n))

c = svg.Canvas(0,0,600,600,map=(-6,6,6,-6))
sym1 = c.symbol("circle",fill="red")
c.grid(1,1)
c.mark(F,ref=sym1)
qplot2(c,F,DF)
c.close()

今回は、参考のために、F の位置をマークしてある。

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。
図13. 渦巻き曲線

注: 3次のベジェ補間があるので、不要

補足2 (2次ベジェ曲線による補間)

2012/12/04

補足1の方法は、点だけを与えて、補間曲線に応用できる。その場合、接線ベクトルは、隣接する2点から決めることになる。端点の接線ベクトルは決められないので、捨てる事になる。

次の図は、図13と計算数のバランスを保つ為に、点の個数を2倍にしている。

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。
図14. 渦巻き曲線2

#encoding: utf-8
from math import exp,sin,cos
import svg
from table import *
# === logarithmic spiral curve
r0 = 5
a = 0.2
b = 1

def f(t):
  r = r0*exp(-a*t)
  x = r*cos(b*t)
  y = r*sin(b*t)
  return x,y

m = 40
F = []
for n in range(0,m):
  F.append(f(0.5*n))

DF = []
for n in range(1,m-1):
  DF.append(ope("-",F[n+1],F[n-1]))

c = svg.Canvas(0,0,600,600,map=(-6,6,6,-6))
sym1 = c.symbol("circle",fill="red")
c.grid(1,1)
c.mark(F,ref=sym1)
# we must discard two end points F[0] and F[m-1], where m=len(F)
qplot2(c,F[1:-1],DF)
c.close()

注: Catmull-Rom/Overhauser spline 補間があるので不要

補足3 (2次ベジェ曲線による補間)

2012/12/04

補間曲線を得る他の考え方は、与えられた3点を通る放物線を求める事である。
P[0],P[1],P[2] を通る放物線は次のようにして求められる。
P[0]P[2]の中点をMとする。その下で点Qを、MQの中点がP[1]となるように定める。そして P[0],Q,P[2] を使って quadratic Bézier 曲線を描く。

結果を図15に示す。

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。
図15. 渦巻き曲線3

この図では説明の為に補助的な情報が書き込まれている。赤丸が P[n]で、×印が Q の位置である。

図14と違って、今度は最初の端点が捨てられていない。最後の端点は生かされることも、捨てられる事もある。他方、欠点としては、n が奇数の点P[n]で、近似曲線の傾きが不連続になるだろう。

これを生成したプログラムを次に示す。

#encoding: utf-8
from math import exp,sin,cos
import svg
from table import *
# === logarithmic spiral curve
r0 = 5
a = 0.2
b = 1

def f(t):
  r = r0*exp(-a*t)
  x = r*cos(b*t)
  y = r*sin(b*t)
  return x,y

m = 40
P = []
for n in range(0,m):
  P.append(f(0.5*n))

Q = []
for n in range(0,(m-1)/2):
  M = cal(P[2*n+2],P[2*n],"+",2,"/")  # Reverse Polish Notation
  Q.append(cal(P[2*n+1],M,"-",2,"*",M,"+"))

c = svg.Canvas(0,0,600,600,map=(-6,6,6,-6))
c.grid(1,1)

sym1 = c.symbol("circle",fill="red")
sym2= c.symbol("cross")
c.mark(P,ref=sym1)
c.mark(Q,ref=sym2)

path = "M %s %s"%tuple(P[0])
for n in range(1,(m+1)/2):
  path = path + " Q %s %s %s %s"%svg.expand(Q[n-1],P[2*n])
c.path(path)
c.close()

譜5. 渦巻き曲線3

注: 4点を通る3次ベジェ曲線による補間法があるので不要

補足4 (問題提起)

2012/12/07

下図16のように、曲線状の3点と、その接線が与えられているとする。曲線そのものは与えられていない。
その下で、3点を通り、与えられた接線に接する Bézier 曲線を求めないと言う問題設定はありそうである。

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。
図16. 3点とその接線

実際に、この問題を解こうとすると、代数的に解くのは至難の業である事がわかる。Qの接線が線分 P[0],P[3] に平行な場合(応用を考えると、このケースが欲しい)でも、無理そう...

文献

ネットで調べてみると、次の文献がこの問題にアプローチしている。

[6] Muhammad Abbas,et al.
"Bezier Curve Interpolation Constrained by a Line"
(Applied Mathematical Sciences, Vol. 5, 2011, no. 37, 1817 - 1832)
http://www.m-hikari.com/ams/ams-2011/ams-37-40-2011/abbasAMS37-40-2011.pdf
この問題に一番近いが.... 不思議な論文だ。

[7] Călin Caba, S.C. Comau Romania S.R.L., Piaţa Ignatie Darabant
"INTERPOLATION WITH CUBIC BÉZIER SPLINES"
(THE INTERNATIONAL CONFERENCE OF THE CARPATHIAN EURO-REGION SPECIALISTS IN INDUSTRIAL SYSTEMS)
http://www.nordtech.ubm.ro/issues/2008/2008.01.09.pdf
これは、与えられた4点を通る Bezier spline の構成法を扱っている。

[8] Kenneth H. Carpenter, EECE KSU
"An introduction to interpolation and splines"
http://www.eece.ksu.edu/~bala/notes/icg/intsplpdf.pdf
目標は同じだが... この論文を読むと、代数的に答えを出すのは難しそう...

補足5 (3次ベジェ曲線による補間)

2012/12/08

与えられた4点を通る3次のペジェ曲線の描き方は Călin の論文(文献[cite:Călin]])から解る。次の図は、この論文に基づいて描いたベジェ曲線である。

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。
図17. 4点を通る3次ベジェ曲線

赤丸で示した4点 C0,C1,C2,C3 を与え、その下でベジェ曲線を決定する補助点 P0,P1,P2,P3を求める。補助点の定義によって、P0 = C0P3 = C3 となる。結果は図では×で記してある。図では計算して得られた P0,P1,P2,P3を基に、実際のベジェ曲線が描かれている。

これを描いたプログラムは

#encoding: utf-8
import svg
from table import *
c = svg.Canvas(0,0,700,800,map=(-3,7,4,-1))
c.grid(1,1,stroke="cyan")
C = [(-2, 0), (-1, 2), (1, 3), (3, 0)]
M = [
  [6,0,0,0],
  [-5,18,-9,2],
  [2,-9,18,-5],
  [0,0,0,6]
  ]
M = cal(M,1.0/6,"*")
P = [None]*len(M)
Cx,Cy = tr(C)
for n in range(0,len(M)):
  x = cal(M[n],Cx,"*","sum")
  y = cal(M[n],Cy,"*","sum")
  P[n] = x,y
sym1 = c.symbol("circle",fill="red")
sym2 = c.symbol("cross")
c.mark(C,ref=sym1,label=("C0","C1","C2","C3"),pad=0.5,anchor="N")
c.mark(P,ref=sym2,label=("P0","P1","P2","P3"),pad=0.5)
path = "M %s %s C %s %s %s %s %s %s"%(expand(P))
c.path(path)
c.close()

この考え方によるプロット関数を cplot1()table module でサポートした。
使い方は次の通り。渦巻き曲線を例題にしている。

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。
図18. 3次ベジェ曲線による補間

#encoding: utf-8
from math import exp,sin,cos
import svg
from table import *
# === logarithmic spiral curve
r0 = 5
a = 0.2
b = 1

def f(t):
  r = r0*exp(-a*t)
  x = r*cos(b*t)
  y = r*sin(b*t)
  return x,y

m = 40
P = []
for n in range(0,m):
  P.append(f(0.5*n))

c = svg.Canvas(0,0,600,600,map=(-6,6,6,-6))
c.grid(1,1)
sym1 = c.symbol("circle",fill="red")
c.mark(P,ref=sym1)
cplot1(c,P)
c.close()

点列を P[0],...,P[n],... とする。cplot1()の欠点は、 n%4 == 3 となる点で接線の連続性が破れることにある。
その例を次ぎに示す。

あなたのブラウザは SVG をサポートしていません。最新のブラウザをお使いください。
図19. 3次ベジェ曲線による補間

座標原点を O で示してある。この点が補間曲線の出発点である。

#encoding: utf-8
import svg
from table import *
c = svg.Canvas(0,0,600,500,map=(-1,6,5,-1.5))
sym = c.symbol("circle",fill="red")
c.grid(1,1,stroke="cyan")
c.text(0,0,text="O", anchor="NE",pad=0.3)
P = [(0, 0), (1, 1), (2, 0), (3, 2), (3, 4), (1, 4), (0, 2)]
cplot1(c,P)
c.mark(P,ref=sym)
c.close()