基於已知重疊將地形剖面拼接在一起

Question

我有一系列的地形剖面掃描，我想合併以創建單個連續剖面。唯一的問題是每次掃描可能是或者可能不是從不同的高度進行的，因此儘管不同文件在覆蓋區域方面有相當多的重疊，但不同的數據可能沒有共同的參考點絕對高度。

下面是4個不同的掃描。每次掃描包含大約30次測量，最後幾次測量代表新數據，其餘測量與前一次掃描重疊。第一次掃描包含唯一已知的絕對值，因此第一次掃描是「黃金標準」。第二次掃描恰好取自同一高度，因此重疊完美匹配（幾乎）完美，並且僅在前一次掃描中增加了4個新點。第三次和第四次掃描是從不同的高度拍攝的，所以雖然重疊覆蓋了相同的區域（相對），但我不能簡單地將它拼接到前兩次掃描上。

Scan1<-c(5,6,7,8,15,16,18,20,25,23,20,17,15,10,10,9,8,9,11,10,13,16,17,19,20,25,28,30,29,30) 
    Scan2<-c(15,16,18,20,25,23,20,16,15,10,10,9,8,9,11,10,13,16,17,19,20,25,28,30,29,30,32,35,38,37) 
    Scan3<-c(28,25,23,18,18,17,16,17,19,18,21,23,25,27,26,33,36,37,37,38,40,43,46,45,43,42,40,38,32,30) 
    Scan4<-c(27,30,29,36,39,39,40,41,43,46,49,48,46,45,43,41,35,33,30,29,28,30,31,32,35) 

使用R，有沒有辦法將這4個掃描縫合在一起，使連續的地形輪廓？絕對高度需要以第一次掃描爲基礎，每次連續掃描被拼接到前一次掃描。將IE-Scan2拼接到掃描1上，並添加4個數據點，然後將掃描3中的新數據添加到掃描1和掃描2的組合中，然後將掃描4中的新數據添加到掃描1,2和3的組合上，等等......

我假設有一種方法可以通過匹配掃描之間的大量重疊來標準化所有數據，使用某種模式識別來發現Scan3是大約8個單元不同的來自Scan1，並且Scan4約有11個單元關閉。但是請注意，我的數據中存在一些「噪音」，並且重疊的模式不適合。

最終結果應該包含一個包含所有4個掃描的完整地形輪廓，並對實際數字的不同進行一些調整。沿着線的東西：

5,6,7,8,15,16,18,20,25,23,20,16.5,15,10,10,9,8,9,11,10,13,15.5,17,19,19,25,28,29.5,29,30,32,35,38,37,35,34,32,30,24,22,19,18,17,19,20,21,24  

Answer 1

你可能要考慮序列比對 - DNA比對，基本上這個問題，但與鹼而不是數字。

作爲一次快速的翻轉，這裏有一個快速編寫函數來找到「最佳」移位，這是基於找到掃描時滑動值之間的最小差值標準偏差。該函數在給定兩個序列，並將其與（默認爲-15〜15）給出的轉移比較：

aligner <- function(bestsequence, sequence2, shift = (-15):15){ 
    minsd <- sd(bestsequence[1:min(length(sequence2), length(bestsequence))] - sequence2[1:min(length(sequence2), length(bestsequence))]) 
    bestshift <- 0 
    avgdiff <- mean(bestsequence[1:min(length(sequence2), length(bestsequence))] - sequence2[1:min(length(sequence2), length(bestsequence))]) 
    for(i in shift){ 
    if(i < 0){ 
     worksequence2 <- sequence2[abs(i):length(sequence2)] 
     if(sd(bestsequence[1:min(length(worksequence2), length(bestsequence))] 
      - worksequence2[1:min(length(worksequence2), length(bestsequence))]) < minsd){ 
     minsd <- sd(bestsequence[1:min(length(worksequence2), length(bestsequence))]- 
         worksequence2[1:min(length(worksequence2), length(bestsequence))]) 
     bestshift <- i 
     avgdiff <- mean(bestsequence[1:min(length(worksequence2), length(bestsequence))]- 
          worksequence2[1:min(length(worksequence2), length(bestsequence))]) 
     } 
    } 
    if(i > 0){ 
     workbest <- bestsequence[i:length(bestsequence)] 
     if(sd(workbest[1:min(length(sequence2), length(workbest))] 
      -sequence2[1:min(length(sequence2), length(workbest))]) < minsd){ 
     minsd <- sd(workbest[1:min(length(sequence2), length(workbest))]- 
         sequence2[1:min(length(sequence2), length(workbest))]) 
     bestshift <- i 
     avgdiff <- mean(workbest[1:min(length(sequence2), length(workbest))]- 
         sequence2[1:min(length(sequence2), length(workbest))]) 
     } 
    } 
    } 
    return(list(bestshift = bestshift, avgdiff = avgdiff, minsd = minsd)) 
} 

所以，爲您的數據：

aligner(Scan1, Scan2) 

$bestshift 
[1] 5 

$avgdiff 
[1] 0.03846154 

$minsd 
[1] 0.1961161 

所以，你Scan2s第五元素等於Scan1的第一個元素。從這裏開始，它應該很容易通過avgdiff進行子集化，修正並附加新的數據點，然後重新運行。

編輯：下面是如何讓你的最終序列。首先，我們需要一個包裝器來輸出我們想要的序列。它基本上運行前面的命令，然後檢查是否移位是正還是負，然後輸出正確的順序：

wrappedaligner <- function(bestseq, seq2){ 
    z <- aligner(bestseq, seq2) 
    if(z$bestshift==0){ 
    if(length(bestseq) >= length(seq2)){ 
    return(bestseq) 
    } else {return(c(bestseq, (seq2[(length(bestseq)+1):length(seq2)])-z$avgdiff))} 
    } 
    else if(z$bestshift > 0){ 
    if(length(bestseq)-z$bestshift >= length(seq2)){ 
     return(bestseq) 
    } else {return(c(bestseq, seq2[(length(bestseq) - z$bestshift + 2):length(seq2)] - z$avgdiff))} 
    } 
    else if(z$bestshift <0){ 
    if((length(bestseq) - abs(z$bestshift))>= length(seq2)){ 
     return(bestseq) 
    } else {return(c(seq2[1:abs(z$bestshift) - 1] - z$avgdiff, bestseq))} 
    } 
} 

現在，我們需要您的數據遞歸地運行它 - 幸運的是我們可以使用Reduce：

Reduce(wrappedaligner, list(Scan1, Scan2, Scan3, Scan4)) 

[1] 5.00000 6.00000 7.00000 8.00000 15.00000 16.00000 18.00000 20.00000 
[9] 25.00000 23.00000 20.00000 17.00000 15.00000 10.00000 10.00000 9.00000 
[17] 8.00000 9.00000 11.00000 10.00000 13.00000 16.00000 17.00000 19.00000 
[25] 20.00000 25.00000 28.00000 30.00000 29.00000 30.00000 31.96154 34.96154 
[33] 37.96154 36.96154 50.83974 49.83974 47.83974 45.83974 39.83974 37.83974