在R中计算2D样条曲线

Question

我正在尝试计算一条类似Bezier样条曲线，它通过一系列的x-y坐标。例如，Matlab的函数的输出如下（示例链接）：

我相信（不再维护的）grid包曾经可以做到这一点（grid.xspline函数？），但我无法安装存档版本的包，并且没有找到任何确切符合我的要求的示例。

bezier包也看起来很有希望，但它速度非常慢，而且我也无法完全正确地使用它：

library(bezier)

set.seed(1)
n <- 10
x <- runif(n)
y <- runif(n)
p <- cbind(x,y)
xlim <- c(min(x) - 0.1*diff(range(x)), c(max(x) + 0.1*diff(range(x))))
ylim <- c(min(y) - 0.1*diff(range(y)), c(max(y) + 0.1*diff(range(y))))
plot(p, xlim=xlim, ylim=ylim)
text(p, labels=seq(n), pos=3)

bp <- pointsOnBezier(cbind(x,y), n=100)
lines(bp$points)
arrows(bp$points[nrow(bp$points)-1,1], bp$points[nrow(bp$points)-1,2],
  bp$points[nrow(bp$points),1], bp$points[nrow(bp$points),2]
)

正如您所看到的，除了端点之外，它不经过任何点。

我将非常感激您提供的任何指导！

- Marc in the box

---

从技术角度来看，您要找的是Catmull-Rom样条线。它与Bezier曲线属于同一类曲线，但具有不同的“哪些点定义了什么”属性。在这种情况下，点n和n+1之间的每个线段在点n处具有与线(n-1)--(n+1)平行的切线，在点(n+1)处具有切线n--(n+2)，切线向量的长度由您选择的“张力”确定。在张力为1/2时，它是立方Bezier曲线的等效曲线（是的，这意味着第一个和最后一个点的切线需要“编造”=） - Mike 'Pomax' Kamermans

Answer 1

其实用 grid 没有必要。你可以从 graphics 包中访问 xspline。

跟随你的代码和 @mrflick 的 shape：

set.seed(1)
n <- 10
x <- runif(n)
y <- runif(n)
p <- cbind(x,y)
xlim <- c(min(x) - 0.1*diff(range(x)), c(max(x) + 0.1*diff(range(x))))
ylim <- c(min(y) - 0.1*diff(range(y)), c(max(y) + 0.1*diff(range(y))))
plot(p, xlim=xlim, ylim=ylim)
text(p, labels=seq(n), pos=3)

你只需要添加一行额外的代码：

xspline(x, y, shape = c(0,rep(-1, 10-2),0), border="red")

- Derek McCrae Norton

---

太棒了！谢谢这个。 - Marc in the box

Answer 2

这可能不是最佳的方法，但grid绝对不是无用的。它与R安装包一起作为默认软件包提供。 它是绘图库（如lattice和ggplot）的基础图形引擎。您不需要安装它，只需加载即可。以下是我如何将您的代码翻译为使用grid.xpline

set.seed(1)
n <- 10
x <- runif(n)
y <- runif(n)
xlim <- c(min(x) - 0.1*diff(range(x)), c(max(x) + 0.1*diff(range(x))))
ylim <- c(min(y) - 0.1*diff(range(y)), c(max(y) + 0.1*diff(range(y))))

library(grid)
grid.newpage()
pushViewport(viewport(xscale=xlim, yscale=ylim))
grid.points(x, y, pch=16, size=unit(2, "mm"), 
    default.units="native")
grid.text(seq(n), x,y, just=c("center","bottom"), 
    default.units="native")
grid.xspline(x, y, shape=c(0,rep(-1, 10-2),0), open=TRUE, 
    default.units="native")
popViewport()

这导致

请注意，网格是相当低级的，所以它不是非常容易使用，但它确实允许您更多地控制绘图的内容和位置。

如果你想提取曲线上的点而不是绘制它，请查看?xsplinePoints帮助页面。

- MrFlick

---

您可以使用 extendrange 扩展向量的范围。 - Vincent Guillemot

@VincentGuillemot - 感谢您提供的有用函数。 - Marc in the box

@MrFlick - 像往常一样，非常好的答案 - 非常感谢。这肯定会让我入门的。我只需要习惯使用那些低级函数。干杯。 - Marc in the box

从点6到点7的这一段有些奇怪；在p7处给出与线段p6--p8平行的切线，这一段不应该有任何拐点，但看起来好像有。 - Mike 'Pomax' Kamermans

Answer 3

感谢所有帮助的人。我总结了所学到的教训和其他方面。

Catmull-Rom样条和三次B样条

在xspline函数中，负的形状值返回一个通过x-y点的Catmull-Rom类型样条曲线。正的形状值返回一个三次B类型样条曲线。零形状值返回一个尖锐的角。如果只给出一个形状值，则该值用于所有点。端点的形状始终被视为尖锐的角（shape=0），其他值不会影响端点处的样条曲线。

# Catmull-Rom spline vs. cubic B-spline
plot(p, xlim=extendrange(x, f=0.2), ylim=extendrange(y, f=0.2))
text(p, labels=seq(n), pos=3)
# Catmull-Rom spline (-1)
xspline(p, shape = -1, border="red", lwd=2) 
# Catmull-Rom spline (-0.5)
xspline(p, shape = -0.5, border="orange", lwd=2) 
# cubic B-spline (0.5)
xspline(p, shape = 0.5, border="green", lwd=2) 
# cubic B-spline (1)
xspline(p, shape = 1, border="blue", lwd=2)
legend("bottomright", ncol=2, legend=c(-1,-0.5), title="Catmull-Rom spline", col=c("red", "orange"), lty=1)
legend("topleft", ncol=2, legend=c(1, 0.5), title="cubic B-spline", col=c("blue", "green"), lty=1)

从xspline中提取结果以进行外部绘图

这需要一些搜索，但诀窍是将参数draw=FALSE应用于xspline。

# Extract xy values
plot(p, xlim=extendrange(x, f=0.1), ylim=extendrange(y, f=0.1))
text(p, labels=seq(n), pos=3)
spl <- xspline(x, y, shape = -0.5, draw=FALSE) 
lines(spl)
arrows(x0=(spl$x[length(spl$x)-0.01*length(spl$x)]), y0=(spl$y[length(spl$y)-0.01*length(spl$y)]),
       x1=(spl$x[length(spl$x)]), y1=(spl$y[length(spl$y)])
)

- Marc in the box