xarray-spatial/xrspatial/perlin.py at master · Geospatial-Data-Science/xarray-spatial · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
import numpy as np

import xarray as xr
from xarray import DataArray

from xrspatial.utils import ngjit


# TODO: change parameters to take agg instead of height / width
def perlin(width: int,
           height: int,
           freq: tuple = (1, 1),
           seed: int = 5) -> xr.DataArray:
    """
    Generate perlin noise aggregate.

    Parameters:
    ----------
    width: int
        Width of output aggregate array.
    height: int
        Height of output aggregate array.
    freq: tuple (default = (1,1))
        (x, y) frequency multipliers.
    seed: int (default = 5)
        Seed for random number generator.

    Returns:
    ----------
    xarray.DataArray
        2D array, of the same type as the input, of calculated perlin noise values.

    Notes:
    ----------
    Algorithm References:
        numba-ized from Paul Panzer example available here:
        https://stackoverflow.com/questions/42147776/producing-2d-perlin-noise-with-numpy
        http://www.mountaincartography.org/mt_hood/pdfs/nighbert_bump1.pdf

    Examples:
    ----------
    Imports
    >>> import numpy as np
    >>> import xarray as xr
    >>> from xrspatial import perlin

    Generate Perlin Aggregate
    >>> print(perlin(5, 5))
    <xarray.DataArray (y: 5, x: 5)>
    array([[0.38502038, 0.3235394 , 0.13230299, 0.02275815, 0.13502038],
           [0.69650136, 0.6169794 , 0.34002832, 0.10065245, 0.10962136],
           [1.        , 0.90777853, 0.55206348, 0.16047902, 0.01192   ],
           [1.        , 0.92388163, 0.60956174, 0.21797728, 0.        ],
           [0.69650136, 0.6604194 , 0.50362745, 0.29227467, 0.15306136]])
    Dimensions without coordinates: y, x
    Attributes:
        res:      1
    """

    linx = range(width)
    liny = range(height)
    linx = np.linspace(0, 1, width, endpoint=False)
    liny = np.linspace(0, 1, height, endpoint=False)
    x, y = np.meshgrid(linx, liny)
    data = _perlin(x * freq[0], y * freq[1], seed=seed)
    data = (data - np.min(data))/np.ptp(data)
    return DataArray(data, dims=['y', 'x'], attrs=dict(res=1))


@ngjit
def _lerp(a, b, x):
    return a + x * (b-a)


@ngjit
def _fade(t):
    return 6 * t**5 - 15 * t**4 + 10 * t**3


@ngjit
def _gradient(h, x, y):
    vectors = np.array([[0, 1], [0, -1], [1, 0], [-1, 0]])
    dim_ = h.shape
    out = np.zeros(dim_)
    for j in range(dim_[1]):
        for i in range(dim_[0]):
            f = np.mod(h[i, j], 4)
            g = vectors[f]
            out[i, j] = g[0] * x[i, j] + g[1] * y[i, j]
    return out


def _perlin(x, y, seed=0):
    np.random.seed(seed)
    p = np.arange(2**20, dtype=int)
    np.random.shuffle(p)
    p = np.stack([p, p]).flatten()

    # coordinates of the top-left
    xi = x.astype(int)
    yi = y.astype(int)

    # internal coordinates
    xf = x - xi
    yf = y - yi

    # fade factors
    u = _fade(xf)
    v = _fade(yf)

    # noise components
    n00 = _gradient(p[p[xi]+yi], xf, yf)
    n01 = _gradient(p[p[xi]+yi+1], xf, yf-1)
    n11 = _gradient(p[p[xi+1]+yi+1], xf-1, yf-1)
    n10 = _gradient(p[p[xi+1]+yi], xf-1, yf)

    # combine noises
    x1 = _lerp(n00, n10, u)
    x2 = _lerp(n01, n11, u)
    a = _lerp(x1, x2, v)
    return a