SimplexNoise.hpp

/*
 * This is a C++ translation/adaptation of 
 * A fast javascript implementation of simplex noise by Jonas Wagner
 * Copyright 2018 Mauro Joel Schütz <maurojoel@gmail.com>

Based on a speed-improved simplex noise algorithm for 2D, 3D and 4D in Java.
Which is based on example code by Stefan Gustavson (stegu@itn.liu.se).
With Optimisations by Peter Eastman (peastman@drizzle.stanford.edu).
Better rank ordering method by Stefan Gustavson in 2012.


 Copyright (c) 2018 Jonas Wagner

 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal
 in the Software without restriction, including without limitation the rights
 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 copies of the Software, and to permit persons to whom the Software is
 furnished to do so, subject to the following conditions:

 The above copyright notice and this permission notice shall be included in all
 copies or substantial portions of the Software.

 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 SOFTWARE.
 */
#ifndef _SIMPLEXNOISE_HPP
#define _SIMPLEXNOISE_HPP

#include <cmath>
#include <ctime>
#include <string>
#include <sstream>
#include <random>
#include <functional>

class SimplexNoise {

	class masher {
		unsigned int n;
		
		template<typename T>
		inline void mash(std::stringstream& ss, T data) {
			ss << data;
		}
		
		template<typename T, typename... Args>
		inline void mash(std::stringstream& ss, T data, Args... args) {
		  ss << data;
		  mash(args...);
		}
		
	public:
		
		masher(): n(0xefc8249d) {
		}
		
		virtual ~masher() {}

		template<typename... Args>
		double operator()(Args... args) {
		  std::stringstream ss;
		  mash(ss, args...);
		  return operator()(ss.str());
		}

		double operator()(std::string const& data) {
			for (int i = 0; i < data.size(); i++) {
				n += data[i];
				double h = 0.02519603282416938 * (double)n;
				n = (unsigned int)h;
				h -= n;
				h *= n;
				n = (unsigned int)h;
				h -= n;
				n += (h * (double)0x100000000); // 2^32
			}
			return (double)(n) * 2.3283064365386963e-10; // 2^-32
		}
	};

	class alea {
		// Johannes Baagøe <baagoe@baagoe.com>, 2010
		double s0;
		double s1;
		double s2;
		double c;
		
		template<typename T>
		inline void calc(masher& mash, T arg) {
			s0 -= mash(arg);
			if (s0 < 0) {
				s0 += 1;
			}
			s1 -= mash(arg);
			if (s1 < 0) {
				s1 += 1;
			}
			s2 -= mash(arg);
			if (s2 < 0) {
				s2 += 1;
			}
		}

		template<typename T, typename... Args>
		inline void calc(masher& mash, T arg, Args... args) {
		  calc(mash, arg);
		  calc(mash, args...);
		}

	public:
		template<typename... Args>
		alea(Args... args) {
			masher mash;
			s0 = mash(' ');
			s1 = mash(' ');
			s2 = mash(' ');
			c = 1;
			calc(mash, args...);
		}
		
		alea() {
			masher mash;
			s0 = mash(' ');
			s1 = mash(' ');
			s2 = mash(' ');
			c = 1;
			calc(mash, std::time(0));
		}
		
		virtual ~alea() {}

		double operator()() {
			double t = 2091639.0 * s0 + c * 2.3283064365386963e-10; // 2^-32
			s0 = s1;
			s1 = s2;
			return s2 = t - (c = (unsigned int)t | 0);
		}
	};
	
	const double F2 { 0.5 * (std::sqrt(3.0) - 1.0) };
	const double G2 { (3.0 - std::sqrt(3.0)) / 6.0 };
	const double F3 { 1.0 / 3.0 };
	const double G3 { 1.0 / 6.0 };
	const double F4 { (std::sqrt(5.0) - 1.0) / 4.0 };
	const double G4 { (5.0 - std::sqrt(5.0)) / 20.0 };

	const double grad3[36] {
		1, 1, 0,
		-1, 1, 0,
		1, -1, 0,

		-1, -1, 0,
		1, 0, 1,
		-1, 0, 1,

		1, 0, -1,
		-1, 0, -1,
		0, 1, 1,

		0, -1, 1,
		0, 1, -1,
		0, -1, -1
	};

	const double grad4[128] {
		0, 1, 1, 1, 0, 1, 1, -1, 0, 1, -1, 1, 0, 1, -1, -1,
		0, -1, 1, 1, 0, -1, 1, -1, 0, -1, -1, 1, 0, -1, -1, -1,
		1, 0, 1, 1, 1, 0, 1, -1, 1, 0, -1, 1, 1, 0, -1, -1,
		-1, 0, 1, 1, -1, 0, 1, -1, -1, 0, -1, 1, -1, 0, -1, -1,
		1, 1, 0, 1, 1, 1, 0, -1, 1, -1, 0, 1, 1, -1, 0, -1,
		-1, 1, 0, 1, -1, 1, 0, -1, -1, -1, 0, 1, -1, -1, 0, -1,
		1, 1, 1, 0, 1, 1, -1, 0, 1, -1, 1, 0, 1, -1, -1, 0,
		-1, 1, 1, 0, -1, 1, -1, 0, -1, -1, 1, 0, -1, -1, -1, 0
	};
	
	std::function<double()> random;
	unsigned char p[256];
	unsigned char perm[512];
	unsigned char permMod12[512];
	
	inline void buildPermutationTable() {
		int i;
		
		for (i = 0; i < 256; i++) {
			p[i] = i;
		}
		for (i = 0; i < 255; i++) {
			unsigned char r = i + ~~(unsigned char)(random() * (256 - i));
			unsigned char aux = p[i];
			p[i] = p[r];
			p[r] = aux;
		}
	}
	
	inline void init() {
		buildPermutationTable();
		for (int i = 0; i < 512; i++) {
			perm[i] = p[i & 255];
			permMod12[i] = perm[i] % 12;
		}
	}
	
public:
	class RandomNumber {
		std::random_device rd;
		std::mt19937 mt;
		std::uniform_real_distribution<double> dist;

	public:
		RandomNumber(): mt(rd()), dist(0.0, 1.0) {}
		RandomNumber(RandomNumber const& rn): rd(), mt(rd()), dist(rn.dist) {}
		virtual ~RandomNumber() {}
		
		double operator()() { return dist(mt); }
	};
	
	SimplexNoise(std::function<double()> rndf = SimplexNoise::RandomNumber()): random(rndf) {
		init();
	}
	
	template<typename... Args>
	SimplexNoise(Args... args): random(SimplexNoise::alea(args...)) {
		init();
	}
	
	virtual ~SimplexNoise() {}
	
	double noise2D(double xin, double yin) {
		double n0 = 0; // Noise contributions from the three corners
		double n1 = 0;
		double n2 = 0;
		
		// Skew the input space to determine which simplex cell we're in
		double s = (xin + yin) * F2; // Hairy factor for 2D
		double i = std::floor(xin + s);
		double j = std::floor(yin + s);
		double t = (i + j) * G2;
		double X0 = i - t; // Unskew the cell origin back to (x,y) space
		double Y0 = j - t;
		double x0 = xin - X0; // The x,y distances from the cell origin
		double y0 = yin - Y0;
		
		// For the 2D case, the simplex shape is an equilateral triangle.
		// Determine which simplex we are in.
		unsigned int i1, j1; // Offsets for second (middle) corner of simplex in (i,j) coords
		if (x0 > y0) {
			i1 = 1;
			j1 = 0;
		}// lower triangle, XY order: (0,0)->(1,0)->(1,1)
		else {
			i1 = 0;
			j1 = 1;
		}
		
		// upper triangle, YX order: (0,0)->(0,1)->(1,1)
		// A step of (1,0) in (i,j) means a step of (1-c,-c) in (x,y), and
		// a step of (0,1) in (i,j) means a step of (-c,1-c) in (x,y), where
		// c = (3-sqrt(3))/6
		double x1 = x0 - i1 + G2; // Offsets for middle corner in (x,y) unskewed coords
		double y1 = y0 - j1 + G2;
		double x2 = x0 - 1.0 + 2.0 * G2; // Offsets for last corner in (x,y) unskewed coords
		double y2 = y0 - 1.0 + 2.0 * G2;
	
		// Work out the hashed gradient indices of the three simplex corners
		unsigned int ii = (unsigned int)i & 255;
		unsigned int jj = (unsigned int)j & 255;

		// Calculate the contribution from the three corners
		double t0 = 0.5 - x0 * x0 - y0 * y0;
		if (t0 >= 0) {
			unsigned int gi0 = permMod12[ii + perm[jj]] * 3;
			t0 *= t0;
			n0 = t0 * t0 * (grad3[gi0] * x0 + grad3[gi0 + 1] * y0); // (x,y) of grad3 used for 2D gradient
		}
		double t1 = 0.5 - x1 * x1 - y1 * y1;
		if (t1 >= 0) {
			unsigned int gi1 = permMod12[ii + i1 + perm[jj + j1]] * 3;
			t1 *= t1;
			n1 = t1 * t1 * (grad3[gi1] * x1 + grad3[gi1 + 1] * y1);
		}
		double t2 = 0.5 - x2 * x2 - y2 * y2;
		if (t2 >= 0) {
			unsigned int gi2 = permMod12[ii + 1 + perm[jj + 1]] * 3;
			t2 *= t2;
			n2 = t2 * t2 * (grad3[gi2] * x2 + grad3[gi2 + 1] * y2);
		}

		// Add contributions from each corner to get the final noise value.
		// The result is scaled to return values in the interval [-1,1].
		return 70.0 * (n0 + n1 + n2);
	}
	
	// 3D simplex noise
	double noise3D(double xin, double  yin, double zin) {
		double n0, n1, n2, n3; // Noise contributions from the four corners
		// Skew the input space to determine which simplex cell we're in
		double s = (xin + yin + zin) * F3; // Very nice and simple skew factor for 3D
		double i = std::floor(xin + s);
		double j = std::floor(yin + s);
		double k = std::floor(zin + s);
		double t = (i + j + k) * G3;
		double X0 = i - t; // Unskew the cell origin back to (x,y,z) space
		double Y0 = j - t;
		double Z0 = k - t;
		double x0 = xin - X0; // The x,y,z distances from the cell origin
		double y0 = yin - Y0;
		double z0 = zin - Z0;
		// For the 3D case, the simplex shape is a slightly irregular tetrahedron.
		// Determine which simplex we are in.
		unsigned int i1, j1, k1; // Offsets for second corner of simplex in (i,j,k) coords
		unsigned int i2, j2, k2; // Offsets for third corner of simplex in (i,j,k) coords
		if (x0 >= y0) {
			if (y0 >= z0) {
				i1 = 1;
				j1 = 0;
				k1 = 0;
				i2 = 1;
				j2 = 1;
				k2 = 0;
			} // X Y Z order
			else if (x0 >= z0) {
				i1 = 1;
				j1 = 0;
				k1 = 0;
				i2 = 1;
				j2 = 0;
				k2 = 1;
			} // X Z Y order
			else {
				i1 = 0;
				j1 = 0;
				k1 = 1;
				i2 = 1;
				j2 = 0;
				k2 = 1;
			} // Z X Y order
		}
		else { // x0<y0
			if (y0 < z0) {
				i1 = 0;
				j1 = 0;
				k1 = 1;
				i2 = 0;
				j2 = 1;
				k2 = 1;
			} // Z Y X order
			else if (x0 < z0) {
				i1 = 0;
				j1 = 1;
				k1 = 0;
				i2 = 0;
				j2 = 1;
				k2 = 1;
			} // Y Z X order
			else {
				i1 = 0;
				j1 = 1;
				k1 = 0;
				i2 = 1;
				j2 = 1;
				k2 = 0;
			} // Y X Z order
		}
		// A step of (1,0,0) in (i,j,k) means a step of (1-c,-c,-c) in (x,y,z),
		// a step of (0,1,0) in (i,j,k) means a step of (-c,1-c,-c) in (x,y,z), and
		// a step of (0,0,1) in (i,j,k) means a step of (-c,-c,1-c) in (x,y,z), where
		// c = 1/6.
		double x1 = x0 - i1 + G3; // Offsets for second corner in (x,y,z) coords
		double y1 = y0 - j1 + G3;
		double z1 = z0 - k1 + G3;
		double x2 = x0 - i2 + 2.0 * G3; // Offsets for third corner in (x,y,z) coords
		double y2 = y0 - j2 + 2.0 * G3;
		double z2 = z0 - k2 + 2.0 * G3;
		double x3 = x0 - 1.0 + 3.0 * G3; // Offsets for last corner in (x,y,z) coords
		double y3 = y0 - 1.0 + 3.0 * G3;
		double z3 = z0 - 1.0 + 3.0 * G3;
		// Work out the hashed gradient indices of the four simplex corners
		unsigned int ii = (unsigned int)i & 255;
		unsigned int jj = (unsigned int)j & 255;
		unsigned int kk = (unsigned int)k & 255;
		// Calculate the contribution from the four corners
		double t0 = 0.6 - x0 * x0 - y0 * y0 - z0 * z0;
		if (t0 < 0) n0 = 0.0;
		else {
			unsigned int gi0 = permMod12[ii + perm[jj + perm[kk]]] * 3;
			t0 *= t0;
			n0 = t0 * t0 * (grad3[gi0] * x0 + grad3[gi0 + 1] * y0 + grad3[gi0 + 2] * z0);
		}
		double t1 = 0.6 - x1 * x1 - y1 * y1 - z1 * z1;
		if (t1 < 0) n1 = 0.0;
		else {
			unsigned int gi1 = permMod12[ii + i1 + perm[jj + j1 + perm[kk + k1]]] * 3;
			t1 *= t1;
			n1 = t1 * t1 * (grad3[gi1] * x1 + grad3[gi1 + 1] * y1 + grad3[gi1 + 2] * z1);
		}
		double t2 = 0.6 - x2 * x2 - y2 * y2 - z2 * z2;
		if (t2 < 0) n2 = 0.0;
		else {
			unsigned int gi2 = permMod12[ii + i2 + perm[jj + j2 + perm[kk + k2]]] * 3;
			t2 *= t2;
			n2 = t2 * t2 * (grad3[gi2] * x2 + grad3[gi2 + 1] * y2 + grad3[gi2 + 2] * z2);
		}
		double t3 = 0.6 - x3 * x3 - y3 * y3 - z3 * z3;
		if (t3 < 0) n3 = 0.0;
		else {
			unsigned int gi3 = permMod12[ii + 1 + perm[jj + 1 + perm[kk + 1]]] * 3;
			t3 *= t3;
			n3 = t3 * t3 * (grad3[gi3] * x3 + grad3[gi3 + 1] * y3 + grad3[gi3 + 2] * z3);
		}
		// Add contributions from each corner to get the final noise value.
		// The result is scaled to stay just inside [-1,1]
		return 32.0 * (n0 + n1 + n2 + n3);
	}
	
	// 4D simplex noise, better simplex rank ordering method 2012-03-09
	double noise4D(double x, double y, double z, double w) {
		double n0, n1, n2, n3, n4; // Noise contributions from the five corners
		// Skew the (x,y,z,w) space to determine which cell of 24 simplices we're in
		double s = (x + y + z + w) * F4; // Factor for 4D skewing
		double i = std::floor(x + s);
		double j = std::floor(y + s);
		double k = std::floor(z + s);
		double l = std::floor(w + s);
		double t = (i + j + k + l) * G4; // Factor for 4D unskewing
		double X0 = i - t; // Unskew the cell origin back to (x,y,z,w) space
		double Y0 = j - t;
		double Z0 = k - t;
		double W0 = l - t;
		double x0 = x - X0; // The x,y,z,w distances from the cell origin
		double y0 = y - Y0;
		double z0 = z - Z0;
		double w0 = w - W0;
		// For the 4D case, the simplex is a 4D shape I won't even try to describe.
		// To find out which of the 24 possible simplices we're in, we need to
		// determine the magnitude ordering of x0, y0, z0 and w0.
		// Six pair-wise comparisons are performed between each possible pair
		// of the four coordinates, and the results are used to rank the numbers.
		double rankx = 0;
		double ranky = 0;
		double rankz = 0;
		double rankw = 0;
		if (x0 > y0) rankx++;
		else ranky++;
		if (x0 > z0) rankx++;
		else rankz++;
		if (x0 > w0) rankx++;
		else rankw++;
		if (y0 > z0) ranky++;
		else rankz++;
		if (y0 > w0) ranky++;
		else rankw++;
		if (z0 > w0) rankz++;
		else rankw++;
		unsigned int i1, j1, k1, l1; // The integer offsets for the second simplex corner
		unsigned int i2, j2, k2, l2; // The integer offsets for the third simplex corner
		unsigned int i3, j3, k3, l3; // The integer offsets for the fourth simplex corner
		// simplex[c] is a 4-vector with the numbers 0, 1, 2 and 3 in some order.
		// Many values of c will never occur, since e.g. x>y>z>w makes x<z, y<w and x<w
		// impossible. Only the 24 indices which have non-zero entries make any sense.
		// We use a thresholding to set the coordinates in turn from the largest magnitude.
		// Rank 3 denotes the largest coordinate.
		i1 = rankx >= 3 ? 1 : 0;
		j1 = ranky >= 3 ? 1 : 0;
		k1 = rankz >= 3 ? 1 : 0;
		l1 = rankw >= 3 ? 1 : 0;
		// Rank 2 denotes the second largest coordinate.
		i2 = rankx >= 2 ? 1 : 0;
		j2 = ranky >= 2 ? 1 : 0;
		k2 = rankz >= 2 ? 1 : 0;
		l2 = rankw >= 2 ? 1 : 0;
		// Rank 1 denotes the second smallest coordinate.
		i3 = rankx >= 1 ? 1 : 0;
		j3 = ranky >= 1 ? 1 : 0;
		k3 = rankz >= 1 ? 1 : 0;
		l3 = rankw >= 1 ? 1 : 0;
		// The fifth corner has all coordinate offsets = 1, so no need to compute that.
		double x1 = x0 - i1 + G4; // Offsets for second corner in (x,y,z,w) coords
		double y1 = y0 - j1 + G4;
		double z1 = z0 - k1 + G4;
		double w1 = w0 - l1 + G4;
		double x2 = x0 - i2 + 2.0 * G4; // Offsets for third corner in (x,y,z,w) coords
		double y2 = y0 - j2 + 2.0 * G4;
		double z2 = z0 - k2 + 2.0 * G4;
		double w2 = w0 - l2 + 2.0 * G4;
		double x3 = x0 - i3 + 3.0 * G4; // Offsets for fourth corner in (x,y,z,w) coords
		double y3 = y0 - j3 + 3.0 * G4;
		double z3 = z0 - k3 + 3.0 * G4;
		double w3 = w0 - l3 + 3.0 * G4;
		double x4 = x0 - 1.0 + 4.0 * G4; // Offsets for last corner in (x,y,z,w) coords
		double y4 = y0 - 1.0 + 4.0 * G4;
		double z4 = z0 - 1.0 + 4.0 * G4;
		double w4 = w0 - 1.0 + 4.0 * G4;
		// Work out the hashed gradient indices of the five simplex corners
		unsigned int ii = (unsigned int)i & 255;
		unsigned int jj = (unsigned int)j & 255;
		unsigned int kk = (unsigned int)k & 255;
		unsigned int ll = (unsigned int)l & 255;
		// Calculate the contribution from the five corners
		double t0 = 0.6 - x0 * x0 - y0 * y0 - z0 * z0 - w0 * w0;
		if (t0 < 0) n0 = 0.0;
		else {
			unsigned int gi0 = (perm[ii + perm[jj + perm[kk + perm[ll]]]] % 32) * 4;
			t0 *= t0;
			n0 = t0 * t0 * (grad4[gi0] * x0 + grad4[gi0 + 1] * y0 + grad4[gi0 + 2] * z0 + grad4[gi0 + 3] * w0);
		}
		double t1 = 0.6 - x1 * x1 - y1 * y1 - z1 * z1 - w1 * w1;
		if (t1 < 0) n1 = 0.0;
		else {
			unsigned int gi1 = (perm[ii + i1 + perm[jj + j1 + perm[kk + k1 + perm[ll + l1]]]] % 32) * 4;
			t1 *= t1;
			n1 = t1 * t1 * (grad4[gi1] * x1 + grad4[gi1 + 1] * y1 + grad4[gi1 + 2] * z1 + grad4[gi1 + 3] * w1);
		}
		double t2 = 0.6 - x2 * x2 - y2 * y2 - z2 * z2 - w2 * w2;
		if (t2 < 0) n2 = 0.0;
		else {
			unsigned int gi2 = (perm[ii + i2 + perm[jj + j2 + perm[kk + k2 + perm[ll + l2]]]] % 32) * 4;
			t2 *= t2;
			n2 = t2 * t2 * (grad4[gi2] * x2 + grad4[gi2 + 1] * y2 + grad4[gi2 + 2] * z2 + grad4[gi2 + 3] * w2);
		}
		double t3 = 0.6 - x3 * x3 - y3 * y3 - z3 * z3 - w3 * w3;
		if (t3 < 0) n3 = 0.0;
		else {
			unsigned int gi3 = (perm[ii + i3 + perm[jj + j3 + perm[kk + k3 + perm[ll + l3]]]] % 32) * 4;
			t3 *= t3;
			n3 = t3 * t3 * (grad4[gi3] * x3 + grad4[gi3 + 1] * y3 + grad4[gi3 + 2] * z3 + grad4[gi3 + 3] * w3);
		}
		double t4 = 0.6 - x4 * x4 - y4 * y4 - z4 * z4 - w4 * w4;
		if (t4 < 0) n4 = 0.0;
		else {
			unsigned int gi4 = (perm[ii + 1 + perm[jj + 1 + perm[kk + 1 + perm[ll + 1]]]] % 32) * 4;
			t4 *= t4;
			n4 = t4 * t4 * (grad4[gi4] * x4 + grad4[gi4 + 1] * y4 + grad4[gi4 + 2] * z4 + grad4[gi4 + 3] * w4);
		}
		// Sum up and scale the result to cover the range [-1,1]
		return 27.0 * (n0 + n1 + n2 + n3 + n4);
	}
	
};

#endif