diff --git a/example/reaction_diffusion.cpp b/example/reaction_diffusion.cpp new file mode 100644 index 000000000..77725a6ca --- /dev/null +++ b/example/reaction_diffusion.cpp @@ -0,0 +1,255 @@ +//---------------------------------------------------------------------------// +// Copyright (c) 2020 Adam Wulkiewicz +// +// Distributed under the Boost Software License, Version 1.0 +// See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt +// +// See http://boostorg.github.com/compute for more information. +//---------------------------------------------------------------------------// + +#include + +#ifdef __APPLE__ +#include +#else +#include +#endif + +namespace bc = boost::compute; + +// Inspired by Reaction-Diffusion Tutorial by Karl Sims +// see: http://www.karlsims.com/rd.html +class reaction +{ +public: + reaction(size_t w, size_t h) + : m_device(bc::system::default_device()) + , m_context(m_device) + , m_queue(m_context, m_device) + , m_program(m_context) + , m_kernel(m_program.program, "my_program") + , m_current(0) + , m_w(w) + , m_h(h) + { + std::vector vec(w * h); + for (size_t j = 0; j < h; j++) + { + for (size_t i = 0; i < w; i++) + { + size_t id = j * w + i; + vec[id].x = 1; + vec[id].y = 0; + } + } + + for (size_t j = h / 2 - 5; j < h / 2 + 5; j++) + { + for (size_t i = w / 2 - 5; i < w / 2 + 5; i++) + { + size_t id = j * w + i; + vec[id].y = 1; + } + } + + m_device_vectors[0] = bc::vector(w * h, m_context); + m_device_vectors[1] = bc::vector(w * h, m_context); + m_current = 0; + + bc::copy(vec.begin(), vec.end(), m_device_vectors[0].begin(), m_queue); + bc::copy(m_device_vectors[0].begin(), m_device_vectors[0].end(), m_device_vectors[1].begin(), m_queue); + } + + template + void next(unsigned n, Vector & texture) + { + BOOST_ASSERT(texture.size() == m_w * m_h); + + for (unsigned i = 0 ; i < n ; i++) + calculate_next(); + + bc::vector const& current_vector = m_device_vectors[m_current]; + std::vector vec(m_w * m_h); + bc::copy(current_vector.begin(), current_vector.end(), vec.begin(), m_queue); + for (size_t k = 0; k < vec.size(); k++) + texture[k] = luminosity(vec[k]); + } + +private: + void calculate_next() + { + size_t next = (m_current + 1) % 2; + bc::vector const& current_vector = m_device_vectors[m_current]; + bc::vector const& next_vector = m_device_vectors[next]; + m_current = next; + + m_kernel.set_arg(0, current_vector.get_buffer()); + m_kernel.set_arg(1, next_vector.get_buffer()); + m_kernel.set_arg(2, (bc::uint_)m_w); + m_kernel.set_arg(3, (bc::uint_)m_h); + + size_t origin[2] = { 1, 1 }; + size_t region[2] = { m_w - 2, m_h - 2 }; + + m_queue.enqueue_nd_range_kernel(m_kernel, 2, origin, region, 0); + m_queue.finish(); + } + + static unsigned char luminosity(bc::float2_ const& f2) + { + //float val = bounded((f2.x - f2.y - 0.5f) * 2.0f, -0.5f, 0.5f) + 0.5; + float val = bounded((f2.x - f2.y) * 2.0f, 0.0f, 1.0f); + return (unsigned char)(255 * val); + } + + template + static T bounded(T val, T mi, T ma) + { + return (std::min)((std::max)(val, mi), ma); + } + + struct program_holder + { + program_holder(bc::context const& context) + : program(bc::program::create_with_source(source(), context)) + { + program.build(); + } + + static const char * source() + { + return BOOST_COMPUTE_STRINGIZE_SOURCE( + __kernel void my_program(__global __read_only float2* curr, + __global __write_only float2* next, + uint w, + uint h) + { + uint i = get_global_id(0); + uint j = get_global_id(1); + + // Parameters + float da = 1.0f; + float db = 0.5f; + //float f = 0.04f; + //float k = 0.0649f; + //float f = 0.0545f; + //float k = 0.062f; + //float f = 0.055f; + //float k = 0.062f; + float f = (float)j / h * (0.07f - 0.03f) + 0.01f; + float k = (float)i / w * (0.07f - 0.045f) + 0.045f; + + // 2D Laplacian + // id_0 - 1 | id_0 | id_0 + 1 + // -------------------------- + // id_1 - 1 | id_1 | id_1 + 1 + // -------------------------- + // id_2 - 1 | id_2 | id_2 + 1 + int id_0 = (j - 1) * w + i; + int id_1 = j * w + i; + int id_2 = (j + 1) * w + i; + float2 l = 0.05f * curr[id_0 - 1] + + 0.2f * curr[id_0] + + 0.05f * curr[id_0 + 1] + + 0.2f * curr[id_1 - 1] + - 1 * curr[id_1] + + 0.2f * curr[id_1 + 1] + + 0.05f * curr[id_2 - 1] + + 0.2f * curr[id_2] + + 0.05f * curr[id_2 + 1]; + + // New values + float a = curr[id_1].x; + float b = curr[id_1].y; + float la = l.x; + float lb = l.y; + float abb = a * b * b; + float na = a + (da * la - abb + f * (1 - a)); + float nb = b + (db * lb + abb - (k + f) * b); + + // Normaliziation + if (na < 0.0f) na = 0.0f; + else if (na > 1.0f) na = 1.0f; + if (nb < 0.0f) nb = 0.0f; + else if (nb > 1.0f) nb = 1.0f; + + next[id_1].x = na; + next[id_1].y = nb; + } + ); + } + + bc::program program; + }; + + bc::device m_device; + bc::context m_context; + bc::command_queue m_queue; + program_holder m_program; + bc::kernel m_kernel; + + bc::vector m_device_vectors[2]; + size_t m_current; + + size_t m_w; + size_t m_h; +}; + +int width = 800; +int height = 800; + +reaction react(width, height); + +void init() +{ + GLuint texture_id; + glEnable(GL_TEXTURE_2D); + glGenTextures(1, &texture_id); + glBindTexture(GL_TEXTURE_2D, texture_id); + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST); + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST); + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP); + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP); +} + +void render_scene() +{ + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); + + glBegin(GL_QUADS); + glTexCoord2f(0, 0); + glVertex3f(-1, -1, 0); + glTexCoord2f(1, 0); + glVertex3f(1, -1, 0); + glTexCoord2f(1, 1); + glVertex3f(1, 1, 0); + glTexCoord2f(0, 1); + glVertex3f(-1, 1, 0); + glEnd(); + + glutSwapBuffers(); +} + +void idle() +{ + std::vector texture(width * height); + react.next(64, texture); + + glTexImage2D(GL_TEXTURE_2D, 0, GL_LUMINANCE, width, height, 0, GL_LUMINANCE, GL_UNSIGNED_BYTE, &texture[0]); + + glutPostRedisplay(); +} + +int main(int argc, char **argv) +{ + glutInit(&argc, argv); + glutInitDisplayMode(GLUT_DEPTH | GLUT_DOUBLE | GLUT_RGBA); + glutInitWindowPosition(100, 100); + glutInitWindowSize(width, height); + glutCreateWindow("Boost.Compute - Reaction Diffusion"); + init(); + glutDisplayFunc(render_scene); + glutIdleFunc(idle); + glutMainLoop(); +}