cloudy-raytracer/common/noise/worleynoise.cpp

#include <chrono>
#include <iostream>
#include <algorithm>
#include <set>
#include <thread>
#include "worleynoise.h"
#include "common/vector3d.h"

void WorleyNoise::runWorleyNoiseInThread(int xOffset, int xSize, WorleyNoise *noise)
{
    noise->renderNoiseThread(xOffset, xSize);
}

void WorleyNoise::generateNoise()
{
    auto start = std::chrono::high_resolution_clock::now();

    // Generate random points
    points.clear();
    points.reserve(pow(numberOfPoints, 3));
    for (int x = 0; x < numberOfPoints; x++)
    {
        for (int y = 0; y < numberOfPoints; y++)
        {
            for (int z = 0; z < numberOfPoints; z++)
            {
                Vector3d point = getRandomPoint();

                // Scale and translate into subcell
                point /= (float) numberOfPoints;

                Vector3d offset = Vector3d(x, y, z);
                offset /= (float) numberOfPoints;
                point += offset;

                points[x + y * numberOfPoints + z * numberOfPoints * numberOfPoints] = point;
            }
        }
    }

    // Generate getNoise map
    noiseMap.clear();
    noiseMap.resize(size * size * size);

    int const nThreads = (int) std::thread::hardware_concurrency() - 1;
    int threadSize = std::floor((float) size / (float) nThreads);
    int remaining = size - nThreads * threadSize;
    std::vector<std::thread> threads;
    for (int n = 0; n < nThreads; n++)
    {
        threads.emplace_back(runWorleyNoiseInThread, n * threadSize, threadSize, this);
    }

    renderNoiseThread(nThreads * threadSize, remaining);

    // Rejoin the threads
    for (int t = 0; t < nThreads; ++t)
    {
        threads[t].join();
    }

    // Normalize cloudNoise map to [0, 1]
    float min = *std::min_element(noiseMap.begin(), noiseMap.end());
    float max = *std::max_element(noiseMap.begin(), noiseMap.end());

    std::for_each(noiseMap.begin(), noiseMap.end(), [min, max](float &value)
    {
        value = (value - min) / (max - min);
    });

    // Duration of computation
    auto stop = std::chrono::high_resolution_clock::now();
    auto duration = std::chrono::duration_cast<std::chrono::seconds>(stop - start);

    std::cout << "Finished computing Worley noise for size " << size << " and " << numberOfPoints << " points in "
              << duration.count() << " seconds" << std::endl;
}

void WorleyNoise::renderNoiseThread(int xOffset, int xSize)
{
    for (int x = xOffset; x < xOffset + xSize; x++)
    {
        for (int y = 0; y < size; y++)
        {
            for (int z = 0; z < size; z++)
            {
                Vector3d point = Vector3d(x, y, z);
                point /= (float) size;
                setNoise(x, y, z, distanceToClosestPoint(point));
            }
        }
    }
}

WorleyNoise::WorleyNoise(int size, int numberOfPoints, unsigned int seed) : numberOfPoints(numberOfPoints), Noise(size)
{
    // Init random
    if (seed == 0)
    {
        std::random_device rd;
        seed = rd();
    }
    this->generator = std::mt19937(seed);
    this->distribution = std::uniform_real_distribution<float>(0.0f, 1.0f);

    generateNoise();
}

Vector3d WorleyNoise::getRandomPoint()
{
    return Vector3d(distribution(generator), distribution(generator), distribution(generator));
}

float WorleyNoise::distanceToClosestPoint(Vector3d point)
{
    std::vector<Vector3d> closePoints = getSubcellPoints(point);

    // Iterate over all subcells
    float minDistance = INFINITY;
    for (Vector3d p: closePoints)
    {
        float distance = length(p - point);
        if (distance < minDistance)
        {
            minDistance = distance;
        }
    }
    return minDistance;
}

std::vector<Vector3d> WorleyNoise::getSubcellPoints(Vector3d point)
{
    std::vector<Vector3d> closePoints;

    point *= numberOfPoints;
    point.x = std::floor(point.x);
    point.y = std::floor(point.y);
    point.z = std::floor(point.z);

    // Iterate over all subcells
    for (int x = -1; x <= 1; x++)
    {
        for (int y = -1; y <= 1; y++)
        {
            for (int z = -1; z <= 1; z++)
            {
                Vector3d offset = Vector3d(x, y, z);
                Vector3d subcell = point + offset;

                Vector3d cellOffsets = Vector3d(0, 0, 0);

                if (subcell.x < 0)
                {
                    cellOffsets.x = -1;
                }
                if (subcell.y < 0)
                {
                    cellOffsets.y = -1;
                }
                if (subcell.z < 0)
                {
                    cellOffsets.z = -1;
                }
                if (subcell.x >= numberOfPoints)
                {
                    cellOffsets.x = 1;
                }
                if (subcell.y >= numberOfPoints)
                {
                    cellOffsets.y = 1;
                }
                if (subcell.z >= numberOfPoints)
                {
                    cellOffsets.z = 1;
                }

                // Wrap around
                subcell.x = std::fmod(subcell.x + numberOfPoints, numberOfPoints);
                subcell.y = std::fmod(subcell.y + numberOfPoints, numberOfPoints);
                subcell.z = std::fmod(subcell.z + numberOfPoints, numberOfPoints);

                // Get points in subcell
                int index = subcell.x + subcell.y * numberOfPoints + subcell.z * numberOfPoints * numberOfPoints;
                closePoints.push_back(points[index] + cellOffsets);
            }
        }
    }
    return closePoints;
}