曲率のこと その３

• 正負の無限大の曲率を持つ閉曲線としてハート関数を使ってみる

# ハート
n <- 100000
t <- seq(from=0,to=1,length=n+1)

x <- cos(t*2*pi)
y <- (sin(t*2*pi)+(x^2)^(1/3))
x <- x[-1]
y <- y[-1]
plot(x,y)

# 最初の点と最後の点の座標は同じなので、その重複を除く
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")

# (単位)接ベクトルの方向を計算する
theta <- atan(dy.st/dx.st)
theta. <- theta
for(i in 2:length(theta.)){
	tmp1 <- theta.[i]
	tmp2 <- theta.[i] + pi
	if(abs(theta.[i-1]-tmp1) > abs(theta.[i-1]-tmp2)){
		theta.[i:length(theta.)] <- theta.[i:length(theta.)] + pi
	}
}
theta <- rep(0,length(dy.st))
theta.diff <- 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]
	theta.diff[i] <- asin(tmp)
	theta[i] <- theta[i-1]+theta.diff[i]
}
theta.diff[length(theta.diff)] <- theta[length(theta)]-theta[1]-2*pi
theta. <- theta
plot(theta.,type="l")
abline(h=pi*((-100):100)+theta.[1])

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

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

cv.tmp <- cv-mean.cv
a <- 1
cv.posi <- 1 + 5 / (1 + exp(-a*cv.tmp))
par(mfcol=c(1,2))
plot(ds.cum[s],cv.posi[s],type="l")
plot(ds.cum[s],cv.posi[s],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[s],t(cv.posis[,s]),type="l")