補間曲線 (Interpolation Curves)

SVG Path

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

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

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

この SVG ファイルを生成する Python のコードは

#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)

定義

t を 0 から t の実数とする。4点 P₀、P₁、P₂、P₃ を与え、次の R(t) で定まる曲線が 3 次のベジェ曲線である。
R(t) = (1-t)³ P₀ + 3(1-t)² t P₁ + 3(1-t)t²P₂ + t³P₃

R(t) を t で微分すると
dR(t)/dt = 3t²(P₃-P₂)+6(1-t)t(P₂-P₁)+3(1-t)²(P₁-P₀)
これを簡単に DR(t) とすると
DR(0) = 3(P₁-P₀)
DR(1) = 3(P₃-P₂)
この関係は、Catmull-Rom spline curve で使われる。

例1

この SVG ファイルを生成する Python のコードは

#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()

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

• svg module の grid() で label
  
• table module の T class で、混合演算
  
• svg module の Canvas class に mark() method を追加した
  

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

A.value()

逆ポーランドの演算を使えば、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,...]

例2

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

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

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

例3

t=0.4 とし、path を

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

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

t=0.5 では

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

例4

次に、t=0.4 で

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

例5

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

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

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

放物線による補間

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

考え方

次の図を考えよう。

この曲線は

c.path("M -1 1 Q 0 -1 1 1")

A=(-1,1), B=(1,1), C=(0,-1)

従って放物線上の3点

(x1,f(x1))
(x2,f(x2))
(x0,f(x0))	# x0 = (x1+x2)/2

上の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)

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

#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

結局次のように 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 ファイルを生成する Python のコードは

#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 は備えていないようである。自分で作らなくてはならない。それが結構くせ者である。

反省

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

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

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

P1 = (1-a)*P0 + a*Q
P3 = (1-a)*P3 + a*Q

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

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

この SVG ファイルを生成する Python のコードは

#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

コントロールポイントを P[0], P[1], ..., P[n]
目標: 3次の spline 曲線を繋げて、全体として、C² 連続(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()

図では 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

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

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

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

Geometry of Uniform B-Spline Curve

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

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

Open Uniform B-Spline Curve

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

文献

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

Catmull-Rom/Overhauser Spline Curve

補題

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

P0 = R(0)
P3 = R(1)
DP0 = DR(0)
DP1 = DR(1)

P1 = P0 + DP0/3
P2 = P3 - DP3/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

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」がある。

#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=1./3 を選ぶ根拠は、ベジェ曲線の微分とベジェ曲線を決める4点との関係である。

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

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

#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 とした。

双曲線 (ライブラリ line() 関数の使用例)

ライブラリの line() 関数を使った実践的な例を挙げる。

この例では21個の点を与え、そのうち19個の点を通る双曲線を描いている。
21個の点は黒丸で示されている。
両端の2点は、精度を上げるために捨てられている。
プログラムコードは、譜4に示されている。

#!/usr/bin/env python3
# Hyperbola

import math
import svg

sinh=math.sinh
cosh=math.cosh

c = svg.Canvas(0,0,600,500,map=(-4,4,4,-4))
c.grid(1,1,stroke="cyan")
c.text(0,0,text="(0,0)", anchor="NE")
c.text(0,1,text="(0,1)",anchor="NE")
c.text(1,0,text="(1,0)",anchor="NE")
sym = c.symbol("circle",fill="black")

P=[]
for n in range(-10,11):
  t = n/5
  P.append((cosh(t),sinh(t)))
for n in range(0,len(P)):
  c.use(P[n],ref=sym)
c.line(*P,smooth="I",stroke="red;width:2")

P=[]
for n in range(-10,11):
  t = n/5
  P.append((sinh(t),cosh(t)))
for n in range(0,len(P)):
  c.use(P[n],ref=sym)
c.line(*P,smooth="I",stroke="red;width:2")

c.close()

譜4: 双曲線

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

渦巻き曲線

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

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

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

この SVG ファイルを生成する Python のコードは

#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()

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

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

#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()

譜5. 渦巻き曲線

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

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

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

#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 の位置をマークしてある。

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

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

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

この SVG ファイルを生成する Python のコードは

#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()

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

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

結果を図16に示す。

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

図15と違って、今度は最初の端点が捨てられていない。最後の端点は生かされることも、捨てられる事もある。他方、欠点としては、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()

譜6. 渦巻き曲線3

補足4 (問題提起)

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

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

また、それだけの価値があるのか否かが問題である。補完法に適用した場合に、区分ごとに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次ベジェ曲線による補間)

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

赤丸で示した4点 C0,C1,C2,C3 を与え、その下でベジェ曲線を決定する補助点 P0,P1,P2,P3を求める。補助点の定義によって、P0 = C0、P3 = 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 ファイルを生成する Python のコードは

#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 となる点で接線の連続性が破れることにある。
その例を次ぎに示す。

座標原点を 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()