R / Perlin噪声中的真实模拟高程数据答案

【问题标题】：Realistic simulated elevation data in R / Perlin noiseR / Perlin噪声中的真实模拟高程数据
【发布时间】：2013-03-13 13:55:26
【问题描述】：

有谁知道如何在 R 中创建模拟栅格高程数据集（即现实高程值的二维矩阵）？ R 的jitter 似乎不合适。在 Java/Processing 中，noise() 函数通过 Perlin noise 算法实现这一点，例如：

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);
  }
}

但我在 R 文献中没有发现对 Perlin 噪声的引用。提前致谢。

【问题讨论】：

查看RandomFields 包（模拟各种高斯随机场的方法），或者分形表面：oceanographerschoice.com/2010/10/…（我有一些旧代码）
快速搜索brings up some C++ code，应该不会太难调整。您可以使用Rcpp 直接使用它或将其翻译为R。
谢谢本，罗兰。 Oceanograherschoice 博客就像我有时间可以尝试的有用的东西，就像学习 C++ 一样。

标签： r perlin-noise

【解决方案1】：

这是 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() )

您可以通过添加这些矩阵来获得更分形的结构，具有不同的网格大小。

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)

【讨论】：

发现 Vincent :) 我花了几个小时涉足其他实现 - 这段代码比我的代码要优雅得多！如果有必要，你知道谁应该归属吗？
也感谢您的编辑。我对 high 和 low 极值点彼此相邻的重复光点很感兴趣。除了这些点之外，该模型似乎很完美。
该伪影是由于在规范化向量时缺少sqrt；我已经修复它并相应地更新了图表。

【解决方案2】：

另一种方法：

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

可以用cbind(sim[[1]], z = sim[[2]])提取对象数据

【讨论】：

漂亮而简单。谢谢！

【解决方案3】：

现在 {ambient} 包中还有一些功能。

【讨论】：