This project contains the following files (right-click files you'd like to download):
perlinnoise.cpp// [header]
// A simple implementation of Perlin and Improved Perlin Noise
// [/header]
// [compile]
// c++ -o perlinnoise -O3 -Wall perlinnoise.cpp
// [/compile]
// [ignore]
// Copyright (C) 2016 www.scratchapixel.com
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
// [/ignore]
#define _USE_MATH_DEFINES
#include <cmath>
#include <cstdio>
#include <random>
#include <functional>
#include <iostream>
#include <fstream>
#include <algorithm>
template<typename T>
class Vec2
{
public:
Vec2() : x(T(0)), y(T(0)) {}
Vec2(T xx, T yy) : x(xx), y(yy) {}
Vec2 operator * (const T &r) const { return Vec2(x * r, y * r); }
Vec2& operator *= (const T &r) { x *= r, y *= r; return *this; }
T x, y;
};
template<typename T>
class Vec3
{
public:
Vec3() : x(T(0)), y(T(0)), z(T(0)) {}
Vec3(T xx, T yy, T zz) : x(xx), y(yy), z(zz) {}
Vec3 operator * (const T &r) const { return Vec3(x * r, y * r, z * r); }
Vec3 operator - (const Vec3 &v) const { return Vec3(x - v.x, y - v.y, z - v.z); }
Vec3& operator *= (const T &r) { x *= r, y *= r, z *= r; return *this; }
T length2() const { return x * x + y * y + z * z; }
Vec3& operator /= (const T &r) { x /= r, y /= r, z /= r; return *this; }
Vec3 cross(const Vec3 &v) const
{
return Vec3(
y * v.z - z * v.y,
z * v.x - x * v.z,
x * v.y - y * v.x
);
}
Vec3& normalize()
{
T len2 = length2();
if (len2 > 0) {
T invLen = T(1) / sqrt(len2);
x *= invLen, y *= invLen, z *= invLen;
}
return *this;
}
friend std::ostream & operator << (std::ostream &os, const Vec3<T> &v)
{
os << v.x << ", " << v.y << ", " << v.z;
return os;
}
T x, y, z;
};
typedef Vec2<float> Vec2f;
typedef Vec3<float> Vec3f;
template<typename T = float>
inline T dot(const Vec3<T> &a, const Vec3<T> &b)
{ return a.x * b.x + a.y * b.y + a.z * b.z; }
template<typename T = float>
inline T lerp(const T &lo, const T &hi, const T &t)
{
return lo * (1 - t) + hi * t;
}
inline
float smoothstep(const float &t)
{
return t * t * (3 - 2 * t);
}
inline
float quintic(const float &t)
{
return t * t * t * (t * (t * 6 - 15) + 10);
}
inline
float smoothstepDeriv(const float &t)
{
return t * (6 - 6 * t);
}
inline
float quinticDeriv(const float &t)
{
return 30 * t * t * (t * (t - 2) + 1);
}
class PerlinNoise
{
public:
PerlinNoise(const unsigned &seed = 2016)
{
std::mt19937 generator(seed);
std::uniform_real_distribution<float> distribution;
auto dice = std::bind(distribution, generator);
for (unsigned i = 0; i < tableSize; ++i) {
#if 0
// bad
float gradientLen2;
do {
gradients[i] = Vec3f(2 * dice() - 1, 2 * dice() - 1, 2 * dice() - 1);
gradientLen2 = gradients[i].length2();
} while (gradientLen2 > 1);
gradients[i].normalize();
#else
// better
float theta = acos(2 * dice() - 1);
float phi = 2 * dice() * M_PI;
float x = cos(phi) * sin(theta);
float y = sin(phi) * sin(theta);
float z = cos(theta);
gradients[i] = Vec3f(x, y, z);
#endif
permutationTable[i] = i;
}
std::uniform_int_distribution<unsigned> distributionInt;
auto diceInt = std::bind(distributionInt, generator);
// create permutation table
for (unsigned i = 0; i < tableSize; ++i)
std::swap(permutationTable[i], permutationTable[diceInt() & tableSizeMask]);
// extend the permutation table in the index range [256:512]
for (unsigned i = 0; i < tableSize; ++i) {
permutationTable[tableSize + i] = permutationTable[i];
}
}
virtual ~PerlinNoise() {}
//[comment]
// Improved Noise implementation (2002)
// This version compute the derivative of the noise function as well
//[/comment]
float eval(const Vec3f &p, Vec3f& derivs) const
{
int xi0 = ((int)std::floor(p.x)) & tableSizeMask;
int yi0 = ((int)std::floor(p.y)) & tableSizeMask;
int zi0 = ((int)std::floor(p.z)) & tableSizeMask;
int xi1 = (xi0 + 1) & tableSizeMask;
int yi1 = (yi0 + 1) & tableSizeMask;
int zi1 = (zi0 + 1) & tableSizeMask;
float tx = p.x - ((int)std::floor(p.x));
float ty = p.y - ((int)std::floor(p.y));
float tz = p.z - ((int)std::floor(p.z));
float u = quintic(tx);
float v = quintic(ty);
float w = quintic(tz);
// generate vectors going from the grid points to p
float x0 = tx, x1 = tx - 1;
float y0 = ty, y1 = ty - 1;
float z0 = tz, z1 = tz - 1;
float a = gradientDotV(hash(xi0, yi0, zi0), x0, y0, z0);
float b = gradientDotV(hash(xi1, yi0, zi0), x1, y0, z0);
float c = gradientDotV(hash(xi0, yi1, zi0), x0, y1, z0);
float d = gradientDotV(hash(xi1, yi1, zi0), x1, y1, z0);
float e = gradientDotV(hash(xi0, yi0, zi1), x0, y0, z1);
float f = gradientDotV(hash(xi1, yi0, zi1), x1, y0, z1);
float g = gradientDotV(hash(xi0, yi1, zi1), x0, y1, z1);
float h = gradientDotV(hash(xi1, yi1, zi1), x1, y1, z1);
float du = quinticDeriv(tx);
float dv = quinticDeriv(ty);
float dw = quinticDeriv(tz);
float k0 = a;
float k1 = (b - a);
float k2 = (c - a);
float k3 = (e - a);
float k4 = (a + d - b - c);
float k5 = (a + f - b - e);
float k6 = (a + g - c - e);
float k7 = (b + c + e + h - a - d - f - g);
derivs.x = du *(k1 + k4 * v + k5 * w + k7 * v * w);
derivs.y = dv *(k2 + k4 * u + k6 * w + k7 * v * w);
derivs.z = dw *(k3 + k5 * u + k6 * v + k7 * v * w);
return k0 + k1 * u + k2 * v + k3 * w + k4 * u * v + k5 * u * w + k6 * v * w + k7 * u * v * w;
}
//[comment]
// classic/original Perlin noise implementation (1985)
//[/comment]
float eval(const Vec3f &p) const
{
int xi0 = ((int)std::floor(p.x)) & tableSizeMask;
int yi0 = ((int)std::floor(p.y)) & tableSizeMask;
int zi0 = ((int)std::floor(p.z)) & tableSizeMask;
int xi1 = (xi0 + 1) & tableSizeMask;
int yi1 = (yi0 + 1) & tableSizeMask;
int zi1 = (zi0 + 1) & tableSizeMask;
float tx = p.x - ((int)std::floor(p.x));
float ty = p.y - ((int)std::floor(p.y));
float tz = p.z - ((int)std::floor(p.z));
float u = smoothstep(tx);
float v = smoothstep(ty);
float w = smoothstep(tz);
// gradients at the corner of the cell
const Vec3f &c000 = gradients[hash(xi0, yi0, zi0)];
const Vec3f &c100 = gradients[hash(xi1, yi0, zi0)];
const Vec3f &c010 = gradients[hash(xi0, yi1, zi0)];
const Vec3f &c110 = gradients[hash(xi1, yi1, zi0)];
const Vec3f &c001 = gradients[hash(xi0, yi0, zi1)];
const Vec3f &c101 = gradients[hash(xi1, yi0, zi1)];
const Vec3f &c011 = gradients[hash(xi0, yi1, zi1)];
const Vec3f &c111 = gradients[hash(xi1, yi1, zi1)];
// generate vectors going from the grid points to p
float x0 = tx, x1 = tx - 1;
float y0 = ty, y1 = ty - 1;
float z0 = tz, z1 = tz - 1;
Vec3f p000 = Vec3f(x0, y0, z0);
Vec3f p100 = Vec3f(x1, y0, z0);
Vec3f p010 = Vec3f(x0, y1, z0);
Vec3f p110 = Vec3f(x1, y1, z0);
Vec3f p001 = Vec3f(x0, y0, z1);
Vec3f p101 = Vec3f(x1, y0, z1);
Vec3f p011 = Vec3f(x0, y1, z1);
Vec3f p111 = Vec3f(x1, y1, z1);
// linear interpolation
float a = lerp(dot(c000, p000), dot(c100, p100), u);
float b = lerp(dot(c010, p010), dot(c110, p110), u);
float c = lerp(dot(c001, p001), dot(c101, p101), u);
float d = lerp(dot(c011, p011), dot(c111, p111), u);
float e = lerp(a, b, v);
float f = lerp(c, d, v);
return lerp(e, f, w); // g
}
private:
/* inline */
uint8_t hash(const int &x, const int &y, const int &z) const
{
return permutationTable[permutationTable[permutationTable[x] + y] + z];
}
//[comment]
// Compute dot product between vector from cell corners to P with predefined gradient directions
//
// perm: a value between 0 and 255
//
// float x, float y, float z: coordinates of vector from cell corner to shaded point
//[/comment]
float gradientDotV(
uint8_t perm, // a value between 0 and 255
float x, float y, float z) const
{
switch (perm & 15) {
case 0: return x + y; // (1,1,0)
case 1: return -x + y; // (-1,1,0)
case 2: return x - y; // (1,-1,0)
case 3: return -x - y; // (-1,-1,0)
case 4: return x + z; // (1,0,1)
case 5: return -x + z; // (-1,0,1)
case 6: return x - z; // (1,0,-1)
case 7: return -x - z; // (-1,0,-1)
case 8: return y + z; // (0,1,1),
case 9: return -y + z; // (0,-1,1),
case 10: return y - z; // (0,1,-1),
case 11: return -y - z; // (0,-1,-1)
case 12: return y + x; // (1,1,0)
case 13: return -x + y; // (-1,1,0)
case 14: return -y + z; // (0,-1,1)
case 15: return -y - z; // (0,-1,-1)
}
}
static const unsigned tableSize = 256;
static const unsigned tableSizeMask = tableSize - 1;
Vec3f gradients[tableSize];
unsigned permutationTable[tableSize * 2];
};
//[comment]
// Simple class to define a polygonal mesh
//[/comment]
class PolyMesh
{
public:
PolyMesh() : vertices(nullptr), st(nullptr), normals(nullptr) {}
~PolyMesh()
{
if (vertices) delete[] vertices;
if (st) delete[] st;
if (normals) delete[] normals;
}
Vec3f *vertices;
Vec2f *st;
Vec3f *normals;
uint32_t *faceArray;
uint32_t *verticesArray;
uint32_t numVertices;
uint32_t numFaces;
void exportToObj();
};
//[comment]
// Export polygonal mesh to OBJ file (vertex positions, texture coordinates and vertex normals)
//[/comment]
void PolyMesh::exportToObj()
{
std::ofstream ofs;
ofs.open("./polyMesh.obj", std::ios_base::out);
for (uint32_t i = 0; i < numVertices; ++i) {
ofs << "v " << vertices[i].x << " " << vertices[i].y << " " << vertices[i].z << std::endl;
}
for (uint32_t i = 0; i < numVertices; ++i) {
ofs << "vt " << st[i].x << " " << st[i].y << std::endl;
}
for (uint32_t i = 0; i < numVertices; ++i) {
ofs << "vn " << normals[i].x << " " << normals[i].y << " " << normals[i].z << std::endl;
}
for (uint32_t i = 0, k = 0; i < numFaces; ++i) {
ofs << "f ";
for (uint32_t j = 0; j < faceArray[i]; ++j) {
uint32_t objIndex = verticesArray[k + j] + 1;
ofs << objIndex << "/" << objIndex << "/" << objIndex << ((j == (faceArray[i] - 1)) ? "" : " ");
}
ofs << std::endl;
k += faceArray[i];
}
ofs.close();
}
//[comment]
// Simple function to create a polygonal grid centred around the origin
//[/comment]
PolyMesh* createPolyMesh(
uint32_t width = 1,
uint32_t height = 1,
uint32_t subdivisionWidth = 40,
uint32_t subdivisionHeight = 40)
{
PolyMesh *poly = new PolyMesh;
poly->numVertices = (subdivisionWidth + 1) * (subdivisionHeight + 1);
std::cerr << poly->numVertices << std::endl;
poly->vertices = new Vec3f[poly->numVertices];
poly->normals = new Vec3f[poly->numVertices];
poly->st = new Vec2f[poly->numVertices];
float invSubdivisionWidth = 1.f / subdivisionWidth;
float invSubdivisionHeight = 1.f / subdivisionHeight;
for (uint32_t j = 0; j <= subdivisionHeight; ++j) {
for (uint32_t i = 0; i <= subdivisionWidth; ++i) {
poly->vertices[j * (subdivisionWidth + 1) + i] = Vec3f(width * (i * invSubdivisionWidth - 0.5), 0, height * (j * invSubdivisionHeight - 0.5));
poly->st[j * (subdivisionWidth + 1) + i] = Vec2f(i * invSubdivisionWidth, j * invSubdivisionHeight);
}
std::cerr << std::endl;
}
poly->numFaces = subdivisionWidth * subdivisionHeight;
poly->faceArray = new uint32_t[poly->numFaces];
for (uint32_t i = 0; i < poly->numFaces; ++i)
poly->faceArray[i] = 4;
poly->verticesArray = new uint32_t[4 * poly->numFaces];
for (uint32_t j = 0, k = 0; j < subdivisionHeight; ++j) {
for (uint32_t i = 0; i < subdivisionWidth; ++i) {
poly->verticesArray[k] = j * (subdivisionWidth + 1) + i;
poly->verticesArray[k + 1] = j * (subdivisionWidth + 1) + i + 1;
poly->verticesArray[k + 2] = (j + 1) * (subdivisionWidth + 1) + i + 1;
poly->verticesArray[k + 3] = (j + 1) * (subdivisionWidth + 1) + i;
k += 4;
}
}
return poly;
}
#define ANALYTICAL_NORMALS 1
int main(int argc, char **argv)
{
PerlinNoise noise;
PolyMesh *poly = createPolyMesh(3, 3, 30, 30);
// displace and compute analytical normal using noise function partial derivatives
Vec3f derivs;
for (uint32_t i = 0; i < poly->numVertices; ++i) {
Vec3f p((poly->vertices[i].x + 0.5), 0, (poly->vertices[i].z + 0.5));
poly->vertices[i].y = noise.eval(p, derivs);
#if ANALYTICAL_NORMALS
Vec3f tangent(1, derivs.x, 0); // tangent
Vec3f bitangent(0, derivs.z, 1); // bitangent
// equivalent to bitangent.cross(tangent)
poly->normals[i] = Vec3f(-derivs.x, 1, -derivs.z);
poly->normals[i].normalize();
#endif
}
#if !ANALYTICAL_NORMALS
// compute face normal if you want
for (uint32_t k = 0, off = 0; k < poly->numFaces; ++k) {
uint32_t nverts = poly->faceArray[k];
const Vec3f &va = poly->vertices[poly->verticesArray[off]];
const Vec3f &vb = poly->vertices[poly->verticesArray[off + 1]];
const Vec3f &vc = poly->vertices[poly->verticesArray[off + nverts - 1]];
Vec3f tangent = vb - va;
Vec3f bitangent = vc - va;
poly->normals[poly->verticesArray[off]] = bitangent.cross(tangent);
poly->normals[poly->verticesArray[off]].normalize();
off += nverts;
}
#endif
poly->exportToObj();
delete poly;
// output noise map to PPM
const uint32_t width = 512, height = 512;
float *noiseMap = new float[width * height];
for (uint32_t j = 0; j < height; ++j) {
for (uint32_t i = 0; i < width; ++i) {
noiseMap[j * width + i] = (noise.eval(Vec3f(i, 0, j) * (1 / 64.), derivs) + 1) * 0.5;
}
}
std::ofstream ofs;
ofs.open("./noise2.ppm", std::ios::out | std::ios::binary);
ofs << "P6\n" << width << " " << height << "\n255\n";
for (unsigned k = 0; k < width * height; ++k) {
unsigned char n = static_cast<unsigned char>(noiseMap[k] * 255);
ofs << n << n << n;
}
ofs.close();
delete[] noiseMap;
return 0;
}