2013-03-13 84 views
5

有誰知道如何在R中創建模擬柵格高程數據集 - 即實際高程值的2d矩陣? R的jitter似乎不合適。在Java /加工的noise()功能實現了這一具有Perlin noise算法例如爲:R/Perlin噪聲中的真實模擬高程數據

size(200, 200); 
float ns = 0.03; // for scaling 
for (float i=0; i<200; i++) { 
    for (float j=0; j<200; j++) { 
    stroke(noise(i*ns, j*ns) * 255); 
    point(i, j); 
    } 
} 

enter image description here

但我發現中的R文學柏林噪聲的引用。提前致謝。

+1

看到任一'RandomFields'包(方法用於模擬廣各種高斯跑了dom字段),或者可能是分形曲面:http://www.oceanographerschoice.com/2010/10/fractal-landscapes-in-r-part-two/(我有一些舊的代碼) – 2013-03-13 14:00:33

+0

快速搜索[帶來一些C++代碼](http://www.dreamincode.net/forums/topic/66480-perlin-noise/),這應該不太難調整。您可以使用'Rcpp'直接使用它或將其翻譯爲R. – Roland 2013-03-13 14:05:12

+0

感謝Ben,Roland。當我有時間的時候,oceanogherschoice博客就像是一件有用的事情,可以學習C++。 – geotheory 2013-03-13 15:49:30

回答

10

以下是在R, 以下解釋的實施方式中 http://webstaff.itn.liu.se/~stegu/TNM022-2005/perlinnoiselinks/perlin-noise-math-faq.html

perlin_noise <- function( 
    n = 5, m = 7, # Size of the grid for the vector field 
    N = 100, M = 100 # Dimension of the image 
) { 
    # For each point on this n*m grid, choose a unit 1 vector 
    vector_field <- apply(
    array(rnorm(2 * n * m), dim = c(2,n,m)), 
    2:3, 
    function(u) u/sqrt(sum(u^2)) 
) 
    f <- function(x,y) { 
    # Find the grid cell in which the point (x,y) is 
    i <- floor(x) 
    j <- floor(y) 
    stopifnot(i >= 1 || j >= 1 || i < n || j < m) 
    # The 4 vectors, from the vector field, at the vertices of the square 
    v1 <- vector_field[,i,j] 
    v2 <- vector_field[,i+1,j] 
    v3 <- vector_field[,i,j+1] 
    v4 <- vector_field[,i+1,j+1] 
    # Vectors from the point to the vertices 
    u1 <- c(x,y) - c(i,j) 
    u2 <- c(x,y) - c(i+1,j) 
    u3 <- c(x,y) - c(i,j+1) 
    u4 <- c(x,y) - c(i+1,j+1) 
    # Scalar products 
    a1 <- sum(v1 * u1) 
    a2 <- sum(v2 * u2) 
    a3 <- sum(v3 * u3) 
    a4 <- sum(v4 * u4) 
    # Weighted average of the scalar products 
    s <- function(p) 3 * p^2 - 2 * p^3 
    p <- s(x - i) 
    q <- s(y - j) 
    b1 <- (1-p)*a1 + p*a2 
    b2 <- (1-p)*a3 + p*a4 
    (1-q) * b1 + q * b2 
    } 
    xs <- seq(from = 1, to = n, length = N+1)[-(N+1)] 
    ys <- seq(from = 1, to = m, length = M+1)[-(M+1)] 
    outer(xs, ys, Vectorize(f)) 
} 

image(perlin_noise()) 

Perlin noise

你可以有一個更分形結構通過添加那些矩陣, 與不同的網格尺寸。

a <- .6 
k <- 8 
m <- perlin_noise(2,2,2^k,2^k) 
for(i in 2:k) 
    m <- m + a^i * perlin_noise(2^i,2^i,2^k,2^k) 
image(m) 
m[] <- rank(m) # Histogram equalization 
image(m) 

Perlin noise, Sum with different grid sizes

+0

Vincent :)現在我花了幾個小時涉及其他實現 - 這個代碼比我的結果要優雅得多!如果有必要,你知道誰應該歸屬嗎? – geotheory 2013-03-13 15:46:33

+0

也感謝編輯。我對_high_和_low_ extreme點彼此相鄰的重複blips感興趣。除了這些點外,這個模型看起來很完美。 – geotheory 2013-03-13 16:49:44

+2

這種人爲現象是由於在對矢量進行歸一化時缺少'sqrt';我修復了它並更新了相應的圖表。 – 2013-03-13 17:59:15

0

的另一種方法:

require(geoR) 
sim <- grf(441, grid="reg", cov.pars=c(1, .25)) 
image(sim, col=gray(seq(1, .1, l=30))) 

enter image description here

可以提取對象數據與cbind(sim[[1]], z = sim[[2]])