r - 玩完美俄罗斯方块：如何使用缩放和平移对齐和缩放两条曲线？

Question

给定y轴（s）上的缩放参数和x轴（t）上的平移参数，当目的是最大化曲线叠加（而不是最小化距离）时，如何缩放和对齐两条不重合的曲线？

正如@DWin 所指出的，这可以重新命名为“如何用 R 完美地玩俄罗斯方块”，尽管它的应用程序远远超出了赢得俄罗斯方块游戏的范围。

这个问题的一个变体可能涉及任意数量的刚体变换（旋转、平移和缩放）。

给定曲线 1

curve1<-data.frame(x=c(1,1,2,2,3),
                   y=c(9,6,6,3,3))

with(curve1, plot(x=x, y=y, type="l", xlim=c(0,10), ylim=c(0,10)))

和曲线 2

curve2<-data.frame(x=c(4,5,5,6,6,7),
                   y=c(2,2,1,1,2,3))

with(curve2, plot(x=x, y=y, type="l", xlim=c(0,10), ylim=c(0,10)))

我希望找到使两条曲线之间的叠加最大化的 s 和 t。

理想情况下，该方法将在 R 中使用 optim。

在这个组成的例子中，t=3 和 s=1/3，所以

t=3
s=1/3

with(curve2, plot(x=x, y=y, type="l", xlim=c(0,10), ylim=c(0,10)))

with(curve1, lines(x=x+t, y=y*s, col="red"))

请注意，要获得这样的拟合，可以具有共识的区域必须具有比不能叠加的区域更高的参数化权重，并且共识区域越大，权重越高。

我一直在探索的路径：

最小化曲线之间的面积

叠加两条曲线的算法

使用形状包进行形状分析

迭代最近点 (ICP)，如果有人在 R 中实现它或可以移植 C 实现之一

使用最大似然的方法的奖励分（假设误差的正态分布）。

Answer 1

当curve1在y轴上缩放一个因子“tfac”并在x轴上移动一个量“s”时，这将返回点之间的距离：

 as.matrix( dist( rbind(as.matrix(curve2), 
                 ( matrix(c(rep(s, 5), rep(1,5)), ncol=2) +   # x-shift matrix
                                       as.matrix(curve1) ) %*% 
                   matrix(c( 1, 0, 0, tfac),ncol=2) ) ) # the y-scaling matrix
                )[ 
                           # better not to use 't' as a variable name
      -(1:5), -(6:11)]  # easier to return the relevant distances when in matrix

将其放入要最小化的函数中应该是一件简单的事情：

 dfunc <- function(C1, C2, s, tfac) { sum( .... ) }

我不确定这是否会返回您期望的结果，因为您暗示的目标函数可能不是距离的总和。您可能需要求助于整数编程方法。优化 CRAN 任务视图将是在 R 中查找这些方法的好地方。我想如果出现此问题，另一种选择可能是舍入“s”值并仅缩放到最接近的 2 次方。

dfunc <- function(param, D1=curve1, D2=curve2) { 
           sum( as.matrix(dist( 
                   rbind(as.matrix(D2), 
                         ( matrix(c(rep(param[1], 5), rep(1,5)), ncol=2) + 
                           as.matrix(D1) ) %*% 
                   matrix(c(1,0,0,param[2]),ncol=2) ) ) )[-(1:5), -(6:11)])}
optim(c(1,1), dfunc)
#$par
$[1] 3.3733977 0.2243866
# trimmed output

使用这些值可以获得以下叠加：

因此，我可能会尝试 s=3, tfac=0.25。（回头看，我从您要求的内容中切换了 t 和 s 的角色。对不起。）

Answer 2

好的，这是一个解决方案的尝试。

基本技巧是这样的：我们光栅化两条曲线，然后我们可以逐块地进行曲线比较。至少，这似乎是比较曲线叠加的一种相当合理的方式。为了鼓励优化器接近曲线，我们还引入了一个损失来惩罚离另一条曲线太远的曲线。

不能保证这适用于更复杂的曲线和变换，但至少这是一个想法。

 curve2<-data.frame(x=c(4,5,5,6,6,7),
                    y=c(2,2,1,1,2,3))
 fillin <- function(ax, ay, bx, by, scaling= 10, steps= 100) floor(cbind(seq(from = ax, to = bx, len = steps), seq(from = ay, to = by, len = steps)) * scaling)
 Bmat <- matrix(0, 100, 100)
 for (i in 2:nrow(curve2)){
 Bmat[fillin (curve2[i-1,1], curve2[i-1,2], curve2[i,1], curve2[i,2])] =1
 }
 Bmat.orig = Bmat

 Bmat = Bmat.orig
 #construct utility function based on 
 #manhattan distances to closest point?
 shift = function(mat, offset){
 mat0 = array(0, dim = dim(mat)+2)
 mat0[1:nrow(mat) +1+ offset[1] , 1:ncol(mat) + 1+offset[2]] = mat
 return(mat0[1:nrow(mat) + 1, 1:ncol(mat) + 1])
 }

 for (i in 1:100){
 Bm = (Bmat != 0)
 Btmp1 = shift(Bm, c(1,0))
 Btmp2 = shift(Bm, c(-1,0))
 Btmp3 = shift(Bm, c(0,1))
 Btmp4 = shift(Bm, c(0,-1))

 Bmat = Bmat + pmax(Bm ,Btmp1, Btmp2, Btmp3, Btmp4)/i
 }

 Bmat2 = replace(Bmat, Bmat == max(Bmat), max(Bmat) + 10)

 #construct and compare rasterised versions
 getcurve = function(trans = c(0,1),  curve=data.frame(x=c(1,1,2,2,3) ,
                    y=c(9,6,6,3,3) ), Bmat = Bmat2){
 Amat = array(0, dim = dim(Bmat))
 curve[,1] = curve[,1] + trans[1]
 curve[,2] = curve[,2] * trans[2]
 for (i in 2:nrow(curve)){
 fillin (curve[i-1,1], curve[i-1,2], curve[i,1], curve[i,2]) -> ind
 if (min(ind) < 1 || max(ind) > nrow(Bmat)) return( array(-1, dim= dim(Bmat)))
 Amat[ind] =1
 }
 Amat = (Amat - mean(Amat))/sd(as.vector(Amat))
 Amat
 }
 compcurve = function(trans = c(0,1), curve=data.frame(x=c(1,1,2,2,3) ,
                    y=c(9,6,6,3,3) ) , Bmat = Bmat2){
 Amat = getcurve(trans, curve, Bmat)
 -sum(Amat * Bmat)
 }
 #SANN seems to work for this, but is slow. Beware of finite differencing
 # - criterion is non-smooth! 
 optim(c(0,1), compcurve, method = "SANN", Bmat = Bmat2) -> output
 image(Bmat)
 contour(getcurve(output$par), add = T)

结果：

不会太破旧吧？

您可能不得不捏造光栅化的细节来解决其他问题。你可能想调整优化的完成方式。

一个“更聪明”的替代方法是注意，对于最佳解决方案，可能至少有一对顶点会重合。这可以为您提供更好的搜索策略。与曲线之间的区域相比，光栅化方案的优势在于它可能更灵活地处理不同的转换和非图形（特别是，第一条曲线中的垂直线是一个问题。）您可以潜在地避免光栅化适当的计算，但只是想想就让我头疼。

由于我们正在最大化标准，因此这是一种最大似然方法。奇怪的是，我认为实际上不可能使用正态错误将这个问题表述为最大似然问题，因为正态错误意味着基于 L2 的标准，众所周知，这不会给你精确的叠加。