フラクタル風景 – 手続き的世界の生成の基本的なアプローチ

私はフラクタル風景を生成し、LaTeXなどで利用可能な機能を活用したいと考えています。

アプローチ:

• Luaを使用してダイヤモンド・スクエア・アルゴリズムによる風景の計算
• pgfplotsを使用して出力することで、多様な表示オプションが利用できる
• メッシュを使用したサーフェスプロット
• Colormapを使用して色付けを行う：海面下は青、山は緑、雪は白、高度によって色が変わる

この解決策は、Marc Lepageによる実装でダイヤモンド・スクエア・アルゴリズムを使用しています。Luaコードはもちろん外部に配置することもでき、それが推奨されています。ここでは扱いやすさのためにドキュメント内に残しています。

• Luaでの計算
• pgfplots Surface-Plotで出力し、colormapを使って色付け（海、山、雪）、viewで視点設定

シード値を変更してバリエーションを加えることができます。これはTerrain関数を呼び出すときの最初のパラメータです。viewだけを変更したい場合はこれを保持します。shader=interpは色を補間しますが、完璧には見えないかもしれません。

計算には時間がかかります。テスト時にはマトリックスの次元やメッシュ行（次元+1）の値を小さくすることがよいでしょう。

\documentclass[border=10pt]{standalone}
\usepackage{pgfplots}
\pgfplotsset{compat=1.15}
\usepackage{luacode}
\begin{luacode*}
  function terrain(seed,dimension,options)
    -- inner functions come from the Heightmap module
    -- Module Copyright (C) 2011 Marc Lepage

    local max, random = math.max, math.random

    -- Find power of two sufficient for size
    local function pot(size)
      local pot = 2
      while true do
        if size <= pot then return pot end
        pot = 2*pot
      end
    end

    -- Create a table with 0 to n zero values
    local function tcreate(n)
      local t = {}
      for i = 0, n do t[i] = 0 end
      return t
    end

    -- Square step
    -- Sets map[x][y] from square of radius d using height function f
    local function square(map, x, y, d, f)
      local sum, num = 0, 0
      if 0 <= x-d then
        if   0 <= y-d   then sum, num = sum + map[x-d][y-d], num + 1 end
        if y+d <= map.h then sum, num = sum + map[x-d][y+d], num + 1 end
      end
      if x+d <= map.w then
        if   0 <= y-d   then sum, num = sum + map[x+d][y-d], num + 1 end
        if y+d <= map.h then sum, num = sum + map[x+d][y+d], num + 1 end
      end
      map[x][y] = f(map, x, y, d, sum/num)
    end

    -- Diamond step
    -- Sets map[x][y] from diamond of radius d using height function f
    local function diamond(map, x, y, d, f)
      local sum, num = 0, 0
      if   0 <= x-d   then sum, num = sum + map[x-d][y], num + 1 end
      if x+d <= map.w then sum, num = sum + map[x+d][y], num + 1 end
      if   0 <= y-d   then sum, num = sum + map[x][y-d], num + 1 end
      if y+d <= map.h then sum, num = sum + map[x][y+d], num + 1 end
      map[x][y] = f(map, x, y, d, sum/num)
    end

    -- Diamond square algorithm generates cloud/plasma fractal heightmap
    -- http://en.wikipedia.org/wiki/Diamond-square_algorithm
    -- Size must be power of two
    -- Height function f must look like f(map, x, y, d, h) and return h'
    local function diamondsquare(size, f)
      -- create map
      local map = { w = size, h = size }
      for c = 0, size do map[c] = tcreate(size) end
      -- seed four corners
      local d = size
      map[0][0] = f(map, 0, 0, d, 0)
      map[0][d] = f(map, 0, d, d, 0)
      map[d][0] = f(map, d, 0, d, 0)
      map[d][d] = f(map, d, d, d, 0)
      d = d/2
      -- perform square and diamond steps
      while 1 <= d do
        for x = d, map.w-1, 2*d do
          for y = d, map.h-1, 2*d do
            square(map, x, y, d, f)
          end
        end
        for x = d, map.w-1, 2*d do
          for y = 0, map.h, 2*d do
            diamond(map, x, y, d, f)
          end
        end
        for x = 0, map.w, 2*d do
          for y = d, map.h-1, 2*d do
            diamond(map, x, y, d, f)
          end
        end
        d = d/2
      end
      return map
    end

    -- Default height function
    -- d is depth (from size to 1 by powers of two)
    -- h is mean height at map[x][y] (from square/diamond of radius d)
    -- returns h' which is used to set map[x][y]
    function defaultf(map, x, y, d, h)
      return h + (random()-0.5)*d
    end

    -- Create a heightmap using the specified height function (or default)
    -- map[x][y] where x from 0 to map.w and y from 0 to map.h
    function create(width, height, f)
      f = f and f or defaultf
      -- make heightmap
      local map = diamondsquare(pot(max(width, height)), f)
      -- clip heightmap to desired size
      for x = 0, map.w do for y = height+1, map.h do map[x][y] = nil end end
      for x = width+1, map.w do map[x] = nil end
      map.w, map.h = width, height
      return map
    end

    -- Initialize pseudo random number generator with seed, to be able to reproduce
    math.randomseed(seed)
    map = create(dimension, dimension)
    if options ~= [[]] then
       tex.sprint("\\addplot3["
         .. options .. "] coordinates{")
    else
      tex.sprint("\\addplot3 coordinates{")
    end
    for x = 0, map.w do
      for y = 0, map.h do
         tex.sprint("("..x..","..y..","..map[x][y]..")")
      end
    end
    tex.sprint("};")
  end
\end{luacode*}
\begin{document}
\begin{tikzpicture}
  \begin{axis}[colormap={terrain}{color(0cm)=(blue!40!black);
    color(1cm)=(blue); color(2cm)=(green!40!black);
    color(4cm)=(green!60!white);color(5cm)=(white!95!black);
    color(8cm)=(white); color(9cm)=(white)},
    hide axis, view={10}{55}]
    \directlua{terrain(14,128,[[surf,mesh/rows=129,mesh/check=false]])}
  \end{axis}
\end{tikzpicture}
\end{document}