Reaction Diffusion Simulation in GLSL

## LICENSE
/*
 * Copyright (c) 2016 David Wicks, sansumbrella.com
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or
 * without modification, are permitted provided that the following
 * conditions are met:
 *
 * Redistributions of source code must retain the above copyright
 * notice, this list of conditions and the following disclaimer.
 * Redistributions in binary form must reproduce the above copyright
 * notice, this list of conditions and the following disclaimer in the
 * documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

## math.glsl
#version 330 core

///
/// Some shared math utilities.
///

float vmax(vec2 v) {
    return max(v.x, v.y);
}

float vmax(vec3 v) {
    return max(max(v.x, v.y), v.z);
}

/// reduce vector values to [0.0, 1.0] while retaining linear relationship
vec3 unbloom(vec3 v) {
    float m = vmax(v);
    if (m > 1.0) {
        return v / m;
    }
    return v;
}

float sdf_square(vec2 center, float half_size, vec2 coord) {
    float d = vmax(abs(center - coord));
    return step(d, half_size);
}

float sdf_circle(vec2 center, float radius, vec2 coord) {
    float l = length(center - coord);
    return step(l, radius);
}

## PingPongApp.cpp
#include "cinder/app/App.h"
#include "cinder/app/RendererGl.h"
#include "cinder/gl/gl.h"
#include "cinder/Log.h"

#include "Watchdog.h"
#include "CinderImGui.h"

using namespace ci;
using namespace ci::app;
using namespace std;

///
/// Demonstrates using Framebuffer Objects as data buffers to run a simulation using GLSL.
/// The two data buffers (source and destination) store floating point data and alternate being read from and written to.
/// The draw calls simply draw fullscreen quads to cover the buffers (in opengl that's just [-1, 1], [-1, 1] without any projection).
///
/// This uses [WatchDog]( https://github.com/simongeilfus/Watchdog ) to automatically reload files.
/// If you don't have/want to download watchdog, you can remove the includes/usage of watchdog and
/// use key presses to reload shaders while the app is running.
///

class PingPongShadersApp : public App {
public:
    void setup() override;
    void keyDown(KeyEvent event) override;
    void update() override;
    void draw() override;
    void mouseDown(MouseEvent event) override;
    void mouseUp(MouseEvent event) override;
    void mouseDrag(MouseEvent event) override;

    void loadSimulationShader();
    void loadRenderShader();
    void seedSimulation();
    void addSquare();
private:
    gl::FboRef   _source, _destination;
    gl::BatchRef _simulation_batch;
    gl::BatchRef _render_batch;
    vec2 _mouse_pos;
    bool _drawing = false;
};

void PingPongShadersApp::loadSimulationShader()
{
    try {
        auto glsl = gl::GlslProg::create(loadAsset("shaders/simulate.vs"), loadAsset("shaders/simulate.fs"));
        glsl->uniform("texture_0", 0);
        glsl->uniform("pixel_size", vec2(1.0f) / vec2(_destination->getSize()));
        _simulation_batch = gl::Batch::create(geom::Rect(Rectf(-1.0f, -1.0f, 1.0f, 1.0f)).texCoords(vec2(0, 0), vec2(1, 0), vec2(1, 1), vec2(0, 1)), glsl);
        CI_LOG_I("Reloaded simulation shader.");
    }
    catch (const std::exception &exc) {
        CI_LOG_W("Unable to load simulation shaders: " << exc.what());
    }
}

void PingPongShadersApp::loadRenderShader()
{
    try {
        auto glsl = gl::GlslProg::create(loadAsset("shaders/render.vs"), loadAsset("shaders/render.fs"));
        _render_batch = gl::Batch::create(geom::Rect(Rectf(-1.0f, -1.0f, 1.0f, 1.0f)).texCoords(vec2(0, 0), vec2(1, 0), vec2(1, 1), vec2(0, 1)), glsl);
        CI_LOG_I("Reloaded render shader.");
    }
    catch (const std::exception &exc) {
        CI_LOG_W("Unable to load render shaders: " << exc.what());
    }
}

void PingPongShadersApp::setup()
{
    ui::initialize();

    gl::Fbo::Format fmt;
    // We only need 2 color channels to simulate chemicals A and B
    // If your GPU doesn't like GL_RG32F, you can try some other values, like: GL_RGB32F, GL_RGB16F
    fmt.setColorTextureFormat(gl::Texture::Format().internalFormat(GL_RG32F).wrap(GL_CLAMP_TO_EDGE));
    _source = gl::Fbo::create(512, 512, fmt);
    _destination = gl::Fbo::create(512, 512, fmt);
    seedSimulation();
    addSquare();

    wd::watch("shaders/simulate.*", [this] (const fs::path &path) {
        loadSimulationShader();
    });

    wd::watch("shaders/render.*", [this] (const fs::path &path) {
        loadRenderShader();
    });
}

void PingPongShadersApp::seedSimulation()
{
    gl::ScopedFramebuffer fbo(_source);
    gl::clear(Color(1, 0, 0)); // fill with chemical "a"
}

void PingPongShadersApp::addSquare()
{
    gl::ScopedFramebuffer fbo(_source);
    gl::ScopedViewport viewport(ivec2(0), _source->getSize());
    gl::ScopedBlendAdditive blend;
    gl::ScopedMatrices matrices;
    gl::setMatricesWindow(_source->getSize());

    gl::ScopedColor color(Color(0, 1, 0)); // draw in chemical "b"
    vec2 center = vec2(_source->getSize()) * 0.5f;
    gl::drawSolidRect(Rectf(center - vec2(5), center + vec2(5)));
}

void PingPongShadersApp::mouseDown(cinder::app::MouseEvent event)
{
    _drawing = true;
    _mouse_pos = event.getPos();
}

void PingPongShadersApp::mouseUp(cinder::app::MouseEvent event)
{
    _drawing = false;
}

void PingPongShadersApp::mouseDrag(cinder::app::MouseEvent event)
{
    _mouse_pos = event.getPos();
}

void PingPongShadersApp::keyDown(cinder::app::KeyEvent event)
{
    switch (event.getCode()) {
        case KeyEvent::KEY_r:
            loadSimulationShader();
            loadRenderShader();
        break;
        case KeyEvent::KEY_c:
            seedSimulation();
        break;
    }
}

void PingPongShadersApp::update()
{
    if (_drawing) {
        // add chemical b by drawing with the mouse
        gl::ScopedFramebuffer fbo(_source);
        gl::ScopedViewport viewport(ivec2(0), _source->getSize());
        gl::ScopedBlendAdditive blend;
        gl::ScopedMatrices matrices;
        gl::setMatricesWindow(getWindowSize());

        gl::ScopedColor color(Color(0, 1, 0)); // add in chemical "b"
        gl::drawSolidCircle(_mouse_pos, 30.0f);
    }

    if (_simulation_batch) {
        // Run the simulation many times each frame so we can watch it unfold faster.
        for (auto i = 0; i < 12; i += 1) {
            gl::ScopedFramebuffer fbo(_destination);
            gl::clear(Color(0, 0, 0));
            gl::ScopedViewport    viewport(ivec2(0), _destination->getSize());
            gl::ScopedTextureBind tex0(_source->getColorTexture(), 0);

            _simulation_batch->draw();
            std::swap(_source, _destination);
        }
    }
}

void PingPongShadersApp::draw()
{
    gl::clear(Color(0, 0, 0));

    if (_render_batch) {
        gl::ScopedTextureBind tex0(_source->getColorTexture(), 0);
        _render_batch->draw();
    }
}

void prepare_settings(App::Settings *settings)
{
    settings->setWindowSize(960, 960);
}

CINDER_APP(PingPongShadersApp, RendererGl, prepare_settings)

## render.fs
#version 330 core
#include "math.glsl"

uniform sampler2D texture_0;
uniform float ciElapsedSeconds;

in Vertex {
    vec2 tex_coord;
} v;

out vec4 f_color;

void main() {
    // Calculate gradient across surface
    vec2 pixel_size = vec2(1.0 / 512.0);
    float dbdx = (texture(texture_0, v.tex_coord + vec2(pixel_size.x, 0)).y - texture(texture_0, v.tex_coord + vec2(-pixel_size.x, 0)).y) / 2.0;
    float dbdy = (texture(texture_0, v.tex_coord + vec2(0, pixel_size.y)).y - texture(texture_0, v.tex_coord + vec2(0, -pixel_size.y)).y) / 2.0;
    float gradient_softness = 0.5 / pixel_size.x;
    vec3 gx = normalize(vec3(pixel_size.x * gradient_softness, 0, dbdx));
    vec3 gy = normalize(vec3(0, pixel_size.y * gradient_softness, dbdy));
    vec3 normal = normalize(cross(gy, gx));

    float simulation_value = smoothstep(0.1, 0.2, texture(texture_0, v.tex_coord).y);

//    float theta = ciElapsedSeconds * 0.5;
//    vec3 light_direction = normalize(vec3(cos(theta), sin(theta), -1.0));
    vec3 light_direction = normalize(vec3(1.0, -0.33, -1.0));
    float diffuse = max(dot(normal, light_direction), 0.0);
    vec3 reflection = reflect(light_direction, normal);
    const vec3 eye_direction = normalize(vec3(0.0, 0.0, 1.0));
    float specular = max(dot(eye_direction, reflection), 0.0);
    specular = pow(specular, 6.0);

    // heightmap-shaded view
    vec3 a = vec3(0.0, 0.0, 0.0);
    vec3 b = vec3(1.1, 1.05, 1.0);
    vec3 specular_color = vec3(0.33, 1.1, 1.4);
    f_color = vec4(mix(a, b, diffuse) * simulation_value + specular_color * specular, 1.0);
    // value-only black & white view
//    f_color = vec4(vec3(simulation_value), 1);
    // normal-map view
//    f_color = vec4((normal + vec3(1.0)) / 2.0, 1.0);
    // raw rgba (ab01) view.
//    f_color = vec4(texture(texture_0, v.tex_coord).xy, 0, 1);
}

## render.vs
#version 330 core

in vec4 ciPosition;
in vec2 ciTexCoord0;

out Vertex {
    vec2 tex_coord;
} v;

void main()
{
    gl_Position = ciPosition;
    v.tex_coord = ciTexCoord0;
}

## simulate.fs
#version 330 core
#include "math.glsl"

uniform sampler2D texture_0;
uniform vec2 pixel_size;

in Vertex {
    vec2 tex_coord;
} v;

out vec4 f_color;

void main() {
    vec2 uv = v.tex_coord;

    // Read the simulation values at our current position.
    vec2 ab = texture(texture_0, uv).xy;
    // The laplace is a measure of the difference between a cell's value and the average of its neighbor values.
    // We need to add everything up to calculate this.
    vec2 laplace = vec2(0);

    laplace += texture(texture_0, uv + pixel_size * vec2(-1, -1)).xy * 0.05;
    laplace += texture(texture_0, uv + pixel_size * vec2(0, -1)).xy * 0.2;
    laplace += texture(texture_0, uv + pixel_size * vec2(1, -1)).xy * 0.05;
    laplace += texture(texture_0, uv + pixel_size * vec2(-1, 0)).xy * 0.2;
    laplace += ab * (-1.0);
    laplace += texture(texture_0, uv + pixel_size * vec2(1, 0)).xy * 0.2;
    laplace += texture(texture_0, uv + pixel_size * vec2(-1, 1)).xy * 0.05;
    laplace += texture(texture_0, uv + pixel_size * vec2(0, 1)).xy * 0.2;
    laplace += texture(texture_0, uv + pixel_size * vec2(1, 1)).xy * 0.05;

    float diffusion_a = mix(0.9, 1.03, uv.x);
    float diffusion_b = mix(0.4, 0.58, uv.y);
    float feed = mix(0.05, 0.06, uv.x);
    float kill = mix(0.060, 0.064, uv.y);
    float dt = 1.0;

    float a = ab.x;
    float b = ab.y;
    // Calculate our new simulation values using the Gray-Scott equation, ignoring timestep.
    float ap = clamp(a + ((diffusion_a * laplace.x) - (a * b * b) + (feed * (1.0 - a))) * dt, 0, 1);
    float bp = clamp(b + ((diffusion_b * laplace.y) + (a * b * b) - ((kill + feed) * b)) * dt, 0, 1);

    f_color = vec4(ap, bp, 0, 1);

    // Optionally seed the simulation with an SDF shape here:
    vec4 food_only = vec4(1, 0, 0, 1);
    vec4 food_and_consumer = vec4(1, 1, 0, 1);

    float square = sdf_square(vec2(0.5, 0.55), 0.1, v.tex_coord);
    float circle = sdf_circle(vec2(0.6, 0.3), 0.2, v.tex_coord);
//    f_color = mix(food_only, food_and_consumer, circle);
//    f_color = mix(food_only, food_and_consumer, max(circle, square));
}

## simulate.vs
#version 330 core

in vec4 ciPosition;
in vec2 ciTexCoord0;

out Vertex {
    vec2 tex_coord;
} v;

void main() {
    gl_Position = ciPosition;
    v.tex_coord = ciTexCoord0;
}