2014-11-09

曲率のことその２

曲率曲率半径

曲線を観察したとする。曲線上の点を多数観察していて、その順序はわかるが、その順序パラメタは、特に弧長などとの整合性はない(のが普通)
曲線を関数解析などしてパラメタ・関数表示をしたところで、その弧長関数は単純な関数表示にはならないのであるから(ならないらしい)、数値計算的に力技をすることにする

# 観測点の数
n <- 1000
# 適当に曲線上の点を発生させる
# 周期関数を用いて閉曲線とする
t <- seq(from=0,to=1,length=n+1)

x <- cos(t*2*pi) + cos(t*12*pi) * 0.1 + sin(t*18*pi) * 0.05
y <- sin(t*2*pi) + sin(t*12*pi) * 0.1 + sin(t*16*pi) * 0.03
# 最初の点と最後の点の座標は同じなので、その重複を除く
x <- x[-1]
y <- y[-1]
# こんな曲線(２次元表示とx,y座標表示)
par(mfcol=c(1,2))
plot(x,y,asp=TRUE,type="l")
matplot(cbind(c(x,x),c(y,y)),type="l")
par(mfcol=c(1,2))

差分をとって、差分の変化をプロットする

# 差分をとる(接ベクトル)
dx <- c(diff(x),x[1]-x[n])
dy <- c(diff(y),y[1]-y[n])

par(mfcol=c(1,2))
plot(dx,dy,asp=TRUE,type="l")
# 周期性を確認するべく、２周分をプロットする
matplot(cbind(c(dx,dx),c(dy,dy)),type="l")
par(mfcol=c(1,2))

弧長を考慮する。単位接ベクトルにする

# 弧長を考慮する
ds <- sqrt(dx^2+dy^2)
ds.cum <- cumsum(ds) # これが弧長パラメタ
# 単位接ベクトルにする
dx.st <- dx/ds
dy.st <- dy/ds

par(mfcol=c(1,2))
plot(dx.st^2+dy.st^2,type="l") # 単位接ベクトルのノルムは１
plot(dx.st,dy.st,asp=TRUE,type="l") # 単位接ベクトルなので、それを位置ベクトルとみなせば、単位円が描かれる
par(mfcol=c(1,1))
matplot(cbind(c(dx.st,dx.st),c(dy.st,dy.st)),type="l")

単位接ベクトルの方向を追跡すると、１周した後は、角度が2pi増えている

# (単位)接ベクトルの方向を計算する。外積的計算
theta <- rep(0,length(dy.st))
theta[1] <- atan(dy.st[1]/dx.st[1])
for(i in 2:length(theta)){
	tmp <- dx.st[i-1]*dy.st[i]-dy.st[i-1]*dx.st[i]
	tmp2 <- asin(tmp)
	theta[i] <- theta[i-1]+tmp2
}
theta. <- theta
# 方向は変わりつつ、１周すると、2pi増えている
plot(theta.,type="l")
abline(h=pi*((-100):100)+theta.[1])

曲率を求める

# 曲率っていうのは、単位弧長あたりの角変化なので、次のように推定できる
cv <- c(diff(theta.),theta.[1]+2*pi-theta.[length(theta.)])/ds
# 曲率を、その曲率の値を持つ弧長で重みづけをして足し合わせたものが、「１周分の角度 2pi」
sum(ds*cv)/(2*pi)
plot(cv,type="l")
# 弧長パラメタで表示する
par(mfcol=c(1,2))
plot(ds.cum,cv,type="l")
# y軸を狭めて再表示
plot(ds.cum,cv,type="l",ylim=c(-3,3))
par(mfcol=c(1,1))

さて、この正負両方に無限大になりうる曲率を正の範囲に押し込みたい
- circle packingでは、正の曲率が無限大になると、隣接等サイズ円数は1に収束し、負側に最大化すると、6に収束するから
- $1 + \frac{5}{1+exp(-a (X-\bar(X))}$ なる相対密度を入れるのはどうだろうか。ちなみに、Xは曲率、 $\bar(X)$ は平均曲率(その曲率だったら、その周縁長のときに、正円となるはずであるような曲率
- 問題は、係数aがどうとでもとれること…

# 曲率をシフトして正の数にする
mean.cv <- 2*pi/sum(ds) # 平均曲率

cv.tmp <- cv-mean.cv
a <- 100
cv.posi <- 1 + 5 / (1 + exp(-a*cv.tmp))
par(mfcol=c(1,2))
plot(ds.cum,cv.posi,type="l")
plot(ds.cum,cv.posi,type="l",ylim=c(0,9))
par(mfcol=c(1,1))


as <- 1:100
cv.posis <- matrix(0,length(as),length(cv.tmp))
for(i in 1:length(as)){
	cv.posis[i,] <- 1 + 5 / (1 + exp(-as[i]*cv.tmp))
}
matplot(ds.cum,t(cv.posis),type="l")