r - 完璧なテトリスをプレイ: スケーリングと変換を使用して 2 つの曲線を整列およびスケーリングする方法は?

Question

y 軸 (s) のスケーリングと x 軸 (t) の変換のパラメーターが与えられた場合、目的が曲線の重ね合わせを最大化する (距離を最小化しない) ことである場合に、一致しない 2 つの曲線をスケーリングして整列させる方法は?

@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)))

2 つの曲線間の重ね合わせを最大化する s と t を見つけたいと考えています。

理想的には、メソッドは optim を使用して R にあるでしょう。

この構成例では、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"))

このようなフィッティングを得るには、コンセンサスを持つことができる領域は、重ね合わせることができない領域よりもパラメーター化の重みが高くなければならず、コンセンサス領域が大きいほど重みが高くなることに注意してください。

私が探索しているトレイル：

曲線間の領域を最小化

2 つの曲線を重ね合わせるアルゴリズム

形状パッケージを使用した形状分析

誰かが R で実装した場合、または C 実装の 1 つを移植できる場合は、反復最近点 (ICP)

最尤法を使用する方法のボーナス ポイント (エラーの正規分布を仮定)。

Answer 1

これは、curve1 が係数 "tfac" によって y 軸上でスケーリングされ、量 "s" によって x 軸上で移動されたときのポイント間の距離を返します。

 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

さて、ここで解決策を試みます。

基本的なトリックは次のとおりです。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)

結果：

あまりにもみすぼらしくないでしょうか？

他の問題に対処するには、ラスター化の詳細をごまかす必要があるかもしれません。そして、最適化がどのように行われるかを調整したいかもしれません。

「より賢明な」代替案は、最適な解を得るには、おそらく少なくとも 1 対の頂点が一致することに注意することです。これにより、より良い検索戦略が得られます。曲線間の領域と比較したラスター化スキームの利点は、さまざまな変換や非グラフに対してより柔軟である可能性があることです (特に、最初の曲線の垂直線が問題です)。適切な計算ですが、考えるだけで頭が痛くなります。

基準を最大化するため、これは最尤法です。不思議なことに、この問題を法線誤差を使用した最尤問題と表現することは実際には不可能だと思います。法線誤差は L2 ベースの基準を意味するため、正確な重ね合わせが得られないことはよく知られています。