Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Reaction Diffusion Simulation in GLSL
/*
* 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.
*/
#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);
}
#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)
#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);
}
#version 330 core
in vec4 ciPosition;
in vec2 ciTexCoord0;
out Vertex {
vec2 tex_coord;
} v;
void main()
{
gl_Position = ciPosition;
v.tex_coord = ciTexCoord0;
}
#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));
}
#version 330 core
in vec4 ciPosition;
in vec2 ciTexCoord0;
out Vertex {
vec2 tex_coord;
} v;
void main() {
gl_Position = ciPosition;
v.tex_coord = ciTexCoord0;
}
@sansumbrella

This comment has been minimized.

Copy link
Owner Author

sansumbrella commented Nov 3, 2016

To run:

  1. Create a new Cinder project with TinderBox (add Watchdog in the blocks section).
  2. Paste the .cpp source into your MyNewApp.cpp file.
  3. Copy the glsl files (*.fs, *.vs, *.glsl) into your new application's assets/ directory.

Build and run.

This has only been tested on OSX, but should work cross-platform.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.