geom/noise/perlin.go

package noise

import (
	"math"
	"math/rand"
)

const perlinB = 0x100
const perlinBM = 0xff
const perlinN = 0x1000

// Perlin contains a perlin noise function (both 1D and 2D).
type Perlin struct {
	Seed int64
	p    [perlinB + perlinB + 2]int
	g1   [perlinB + perlinB + 2]float64
	g2   [perlinB + perlinB + 2][2]float64
}

// NewPerlin creates a new perlin noise function with seed seed.
func NewPerlin(seed int64) *Perlin {
	p := &Perlin{Seed: seed}
	p.Init()
	return p
}

func perlinSmoothCurve(t float64) float64 { return t * t * (3. - 2.*t) }

func perlinLerp(t, a, b float64) float64 { return a + t*(b-a) }

func perlinAt2D(q [2]float64, rx, ry float64) float64 { return rx*q[0] + ry*q[1] }

func (p *Perlin) noise1d(x float64) float64 {
	tx := x + perlinN
	bx0 := int(tx) & perlinBM
	bx1 := (bx0 + 1) & perlinBM
	rx0 := tx - math.Floor(tx)
	rx1 := rx0 - 1.

	sx := perlinSmoothCurve(rx0)
	u := rx0 * p.g1[p.p[bx0]]
	v := rx1 * p.g1[p.p[bx1]]

	return perlinLerp(sx, u, v)
}

func (p *Perlin) noise2d(x, y float64) float64 {
	tx := x + perlinN
	bx0 := int(tx) & perlinBM
	bx1 := (bx0 + 1) & perlinBM
	rx0 := tx - math.Floor(tx)
	rx1 := rx0 - 1.

	ty := y + perlinN
	by0 := int(ty) & perlinBM
	by1 := (by0 + 1) & perlinBM
	ry0 := ty - math.Floor(ty)
	ry1 := ry0 - 1.

	i := p.p[bx0]
	j := p.p[bx1]

	b00 := p.p[i+by0]
	b10 := p.p[j+by0]
	b01 := p.p[i+by1]
	b11 := p.p[j+by1]

	sx := perlinSmoothCurve(rx0)
	sy := perlinSmoothCurve(ry0)

	u := perlinAt2D(p.g2[b00], rx0, ry0)
	v := perlinAt2D(p.g2[b10], rx1, ry0)
	a := perlinLerp(sx, u, v)

	u = perlinAt2D(p.g2[b01], rx0, ry1)
	v = perlinAt2D(p.g2[b11], rx1, ry1)
	b := perlinLerp(sx, u, v)

	return perlinLerp(sy, a, b)
}

func normalize2D(v [2]float64) [2]float64 {
	var s float64

	s = math.Sqrt(v[0]*v[0] + v[1]*v[1])
	v[0] /= s
	v[1] /= s
	return v
}

// Init initialised the Perlin noise source. Must be called unless struct was created with NewPerlin(...).
func (p *Perlin) Init() {
	r := rand.New(rand.NewSource(p.Seed))

	for i := 0; i < perlinB; i++ {
		p.p[i] = i
		p.g1[i] = float64((r.Int()%(perlinB+perlinB))-perlinB) / perlinB

		for j := 0; j < 2; j++ {
			p.g2[i][j] = float64((r.Int()%(perlinB+perlinB))-perlinB) / perlinB
		}
		p.g2[i] = normalize2D(p.g2[i])
		for j := 0; j < 3; j++ {
			r.Int()
		}
	}

	for i := perlinB; i > 0; i-- {
		j := r.Int() % perlinB
		p.p[i], p.p[j] = p.p[j], p.p[i]
	}

	for i := 0; i < perlinB+2; i++ {
		p.p[perlinB+i] = p.p[i]
		p.g1[perlinB+i] = p.g1[i]
		for j := 0; j < 2; j++ {
			p.g2[perlinB+i][j] = p.g2[i][j]
		}
	}
}

// Noise1D generates the noise value for x. alpha defines the weight when the sum is formed. Typically it is 2 (function is noisier when alpha approaches 1). beta  is the harmonic scaling/spacing (typically 2). n defines the number of harmonics.
func (p *Perlin) Noise1D(x, alpha, beta float64, n int) float64 {
	var val, sum float64
	var scale float64 = 1
	var coord float64 = x

	for i := 0; i < n; i++ {
		val = p.noise1d(coord)
		sum += val / scale
		scale *= alpha
		coord *= beta
	}
	return sum
}

// Noise2D generates the noise value for coordinate (x, y). alpha defines the weight when the sum is formed. Typically it is 2 (function is noisier when alpha approaches 1). beta  is the harmonic scaling/spacing (typically 2). n defines the number of harmonics.
func (p *Perlin) Noise2D(x, y, alpha, beta float64, n int) float64 {
	var val, sum float64
	var scale float64 = 1
	var coord = [2]float64{x, y}

	for i := 0; i < n; i++ {
		val = p.noise2d(coord[0], coord[1])
		sum += val / scale
		scale *= alpha
		coord[0] *= beta
		coord[1] *= beta
	}
	return sum
}