1

I have two curves of same count of points filled using approx function, for both x and y values separately for each curve. Both x and y axis values are logarithmic, so I convert back to normal decimal scale when approximating and interpolating. Black and blue lines are original lines and the red one is interpolated in between. As you can see the red line doesn't mimic the bend on the right side, since interpolation is performed based on assumption that each x and y pair are the closest.

Is there any way how to perform interpolation between curves in R based on the real closest points in between? Maybe there exists algorithms for that? Anything would be useful as I am not sure how it is called in mathematics.

    base="ftp://cdsarc.u-strasbg.fr/pub/cats/J/A+A/508/355/ms/"
    setwd("~/Desktop")
    file1=paste(base,"z001y23_1.60.dat",sep="")
    file2=paste(base,"z001y23_1.70.dat",sep="")

    cols=c("no","age","logL","logTef", "grav","stage")
    ncol <- length(count.fields(file=file1, sep = ","))
    second=read.table(file=file1,fill=T, blank.lines.skip=F, skip=2, header=F, strip.white=T, col.names = paste("V", seq_len(ncol)))
    second$V.6<-second$V.23
    colnames(second) <-cols
    second$logL=as.numeric(second$logL)
    #performing some filtering of data here
    pos1=which(second$stage == "trgb")[1]
    second=second[1:pos1,]

    ncol <- length(count.fields(file=file2, sep = ","))
    first=read.table(file=file2,fill=T, blank.lines.skip=F, skip=2, header=F, strip.white=T, col.names = paste("V", seq_len(ncol)))
    first$V.6<-first$V.23
    colnames(first) <-cols
    #performing some filtering of data here
    pos2=which(first$stage == "trgb")[1]
    first=first[1:pos2,]

    #plotting data
    len=max(c(min(first[[4]]),min(second[[4]])))
    first=first[first[[4]]>len,]
    second=second[second[[4]]>len,]

    plot(second[[4]],second[[3]],t="l",xlim=rev(range(second[[4]])),xlab="x",ylab="y")
    lines(first[[4]],first[[3]],t="l",col="blue")
    n=max(c(length(second[[4]]),length(first[[4]])))
    #approximating missing points
    xf1 <- approx(10^second[[4]],n=n)
    yf1 <- approx(10^second[[3]],n=n)

    xf2 <- approx(10^first[[4]],n=n)
    yf2 <- approx(10^first[[3]],n=n)

    #calculating interpolated line
    ratio=2
    s1<-log10((xf1$y-xf2$y)/ratio+xf2$y)
    s2<-log10((yf1$y-yf2$y)/ratio+yf2$y)
    lines(s1,s2, col ="red")

enter image description here

4

1 回答 1

3

虽然不是最终答案，这是从我前段时间为流通道迁移所做的改编而来的。请注意，这些通常不是自穿越，因此您的里程可能会有所不同。整个想法是计算曲率并使用动态时间扭曲来匹配极值。

大致可以这样概括：

  1. 对两条曲线进行参数化，因此 L1 和 L2 是表示从曲线开始到相关索引的长度的向量。
  2. 使用 L1 和 L2计算smooth.spline每条曲线的 x 和 y 的 xsp1、ysp1、xsp2、ysp2。注意平滑参数，因为您的曲线有时看起来很尖锐。
  3. 显式获取每条平滑线的有符号曲率
  4. 使用 dtw 匹配每条平滑线的曲率峰值
  5. 使用 dtw 返回的索引建立曲线之间的映射
  6. ...
  7. 利润！！！

请注意，dtw 不会创造奇迹，因此需要进行一些实验。

PS 为了节省您的时间，我尝试在没有曲率的情况下直接在 x 和 y 上使用 dtw，但结果并不好，因为我们希望同时映射两个坐标。

编辑

library(dtw)
df1 <- data.frame(x=first[[4]], y=first[[3]])
df2 <- data.frame(x=second[[4]], y=second[[3]])
measure <- function(df)
  within(df, m <- c(0, cumsum(diff(x)^2 + diff(y)^2)))
df1 <- measure(df1)
df2 <- measure(df2)

curvify <- function(df) {
  xsp <- with(df, smooth.spline(m, x))
  ysp <- with(df, smooth.spline(m, y))
  xx <- predict(xsp, df$m)$y
  yy <- predict(ysp, df$m)$y
  xp <- predict(xsp, df$m, deriv=1)$y
  xpp <- predict(xsp, df$m, deriv=2)$y
  yp <- predict(ysp, df$m, deriv=1)$y
  ypp <- predict(ysp, df$m, deriv=2)$y
  # http://en.wikipedia.org/wiki/Curvature#Signed_curvature
  within(df, c <- (xp*ypp - yp*xpp)/(xp^2 + yp^2)^1.5)
}

df1 <- curvify(df1)
df2 <- curvify(df2)

d <- dtw(df1$c, df2$c, keep=TRUE)
# plot(d, type='three')

xx <- ( df1$x[d$index1] + df2$x[d$index2] ) /2
yy <- ( df1$y[d$index1] + df2$y[d$index2] ) /2

lines(xx, yy, col="green")

在此处输入图像描述

在此处输入图像描述

编辑

使用 1/2 以外的权重进行插值

fr <- 1/3
xx <- df1$x[d$index1] * fr + df2$x[d$index2] * (1-fr)
yy <- df1$y[d$index1] * fr + df2$y[d$index2] * (1-fr)
lines(xx, yy, col="yellow")

fr <- 2/3
xx <- df1$x[d$index1] * fr + df2$x[d$index2] * (1-fr)
yy <- df1$y[d$index1] * fr + df2$y[d$index2] * (1-fr)
lines(xx, yy, col="brown")
于 2014-06-13T01:59:51.610 回答