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

【问题标题】：Playing perfect Tetris: how to align and scale two curves using scaling and translation?玩完美俄罗斯方块：如何使用缩放和平移对齐和缩放两条曲线？
【发布时间】：2012-06-05 07:09:55
【问题描述】：

给定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"))

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

我一直在探索的足迹：

minimize area between curves

algorithm for superposing two curves

使用形状包进行形状分析

Iterative closest point (ICP) 如果有人在 R 中实现它或者可以移植其中一个 C 实现

使用最大似然法的加分（假设误差为正态分布）。

【问题讨论】：

嗯。一个有趣的问题。需要考虑的一件事是，似乎需要做一些事情来避免病理结果。具体来说，局部最优可能是您大规模放大一个形状的情况，因此将所有内容都拟合到直线的一部分......
对如何执行此操作有想法，但有一个简短的问题。假设这两个形状是计算器风格的“5”（即五个相等的水平/垂直线段连接在一起形成“5”形状 - 你肯定明白我的意思）和垂直线段（又名计算器风格“1”）。在这里，我想检查您是否认为最好的解决方案是 100% 重叠的（即第二个形状与第一个形状的垂直部分完全重叠）。换句话说，第一个形状的大部分都未被覆盖，并且与第二个形状相距甚远，这并不重要，对吧？
所以我设法想出了一个有趣的方法来检测精确的叠加 - 但是指定问题的方式我认为你不可避免地会遇到局部最小值的问题，而你没有想远离。请问您是否认为（根据您要解决的任何问题）整个形状 1 可以叠加到形状 2 上？
@TimP，形状 1 和 2 之间的重叠量是任意的。有时，如您所指出的，第一个形状会完全适合第二个形状。

标签： r geometry curve-fitting curve

【解决方案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 的角色从你要求的地方换了。对不起。）

【讨论】：

有趣的努力，谢谢！然而，这似乎是一个比最小化距离更复杂的问题，因为可以达成共识的区域必须比无法叠加的区域具有更高的权重。
我的建议是在我获得该数字结果之前提出的。也可以尝试根据曲线之间的面积进行最小化，但我也可以在该问题规范中看到模棱两可的地方。我不知道 R 的俄罗斯方块播放包。

【解决方案2】：

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

基本诀窍是：我们光栅 em>两条曲线，然后我们可以通过瓷砖的瓷砖进行曲线比较。这似乎是一种相当合理的方式，可以比较曲线叠加，至少。为了鼓励优化器接近曲线，我们还介绍了损失，以惩罚离另一个的曲线太远。

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

 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)

结果：

不太破旧，也许？

你可能必须努力解决光谱的细节，以适应其他问题。而且你可能想调整如何完成oon。

“更智能的”替代方案是为了注意，对于最佳解决方案，可能是至少一对顶点将能够结合的情况。这可以为您提供更好的搜索策略。与面积曲线相比的光栅化方案的优点是，它可能更灵活的WRT到不同的变换和非图形（特别是，您的第一个曲线中的垂直线是一个问题。）您可以避免避免频闪地位适当的计算，但它让我一个头痛只是为了考虑它。

由于我们最大化标准，因此这是最大似然方法。好奇地，我认为它实际上是不可能的 em>使用正常错误的最大似然问题，因为正常错误意味着基于L2的标准，这将着名不是 em>给你准确叠加。

【讨论】：