reaction_diffusion_simulation.cpp

#include <opencv2/opencv.hpp>
#include <opencv2/highgui.hpp>

using namespace cv;

const int rows = 216;
const int columns = 384;
const int chemicals = 2;
bool is_drawing = false;
int brush_size = 2;

// This loops a value. If the value passes [max] then it returns the min and vice versa.
int loop(int value, int min, int max) 
{
    if (value == 0) 
    {
        return min;
    }
    
    int looped_value = value % (max - min);
    
    if (looped_value == 0) 
    {
        return max * (value / abs(value));
    }
    
    return looped_value + min;
}

// This constrains a value between [min] and [max].
float constrain(float value, float min, float max) 
{
    if (value > max) 
    {
        return max;
    }

    if (value < min) 
    {
        return min;
    }

    return value;
}

// This function takes in a cell in a grid and returns a value based on the status of the cells around it.
float laplace(bool is_b_or_a, int x, int y, float cell_grid[rows][columns][chemicals]) 
{
    float corner_weight = 0.05f;
    float cross_weight = 0.2f;
    float laplace_value = 0.0f;
    int b_or_a = (int)is_b_or_a;

    // Since I'm checking for the cells around a cell, what would happend if the cell was at 0, 0?
    // Since the x is zero, x - 1 would be negative one but negative one doesn't exist in an array.
    // My solution to this is to loop x. If the width of the array was 100, -1 would be 99(98 because array logic);
    int looped_y = loop(y, 0, rows);
    int looped_x = loop(x, 0, columns);
    int looped_upper_y = loop(y + 1, 0, (rows - 1));
    int looped_lower_y = loop(y - 1, 0, (rows - 1));
    int looped_upper_x = loop(x + 1, 0, (columns - 1));
    int looped_lower_x = loop(x - 1, 0, (columns - 1));

    // Calculate a value from the status of the cells around a cell.
    laplace_value += cell_grid[looped_y][looped_x][b_or_a] * -1;
    laplace_value += cell_grid[looped_upper_y][x][b_or_a] * cross_weight;
    laplace_value += cell_grid[looped_lower_y][x][b_or_a] * cross_weight;
    laplace_value += cell_grid[y][looped_upper_x][b_or_a] * cross_weight;
    laplace_value += cell_grid[y][looped_lower_x][b_or_a] * cross_weight;
    laplace_value += cell_grid[looped_upper_y][looped_upper_x][b_or_a] * corner_weight;
    laplace_value += cell_grid[looped_lower_y][looped_lower_x][b_or_a] * corner_weight;
    laplace_value += cell_grid[looped_upper_y][looped_lower_x][b_or_a] * corner_weight;
    laplace_value += cell_grid[looped_lower_y][looped_upper_x][b_or_a] * corner_weight;

    return laplace_value;  
}

// The cell gird has 2 chemicals in it. This adds chemical A to the whole grid.
void initializeCellGrid(float cell_grid[rows][columns][chemicals]) 
{
    for (int j = 0; j < rows; ++j) 
    {
        for (int i = 0; i < columns; ++i) 
        {
            // Set chemical a ammount to 1
            cell_grid[j][i][0] = 1;
            cell_grid[j][i][1] = 0;
        }
    }
}

// This adds chemical B to the grid in a circle to start chemical A and B reactions.
void addChemicalB(float cell_grid[rows][columns][chemicals], float ammount_to_add, int x_position, int y_position, int radius) 
{
    for (int y = 0; y < rows; ++y)
    {
        for (int x = 0; x < columns; ++x)
        {
            float delta_x = (float)x - x_position;
            float delta_y = (float)y - y_position;
            float distanceSquared = delta_x * delta_x + delta_y * delta_y;

            if (distanceSquared <= (radius * radius))
            {
                cell_grid[y][x][1] += ammount_to_add;
            }
        }
    }
}

// I guess this needs no explaining
void addChemicalBOnClick(int event, int x, int y, int flags, void* cell_grid) 
{
    int chemical_b_radius = (int)(((brush_size*0.01f) * min(columns, rows))*0.5f);
    float chemical_b_to_add = 0.5f;

    switch (event) 
    {
        case EVENT_LBUTTONDOWN:
            is_drawing = true;   
            addChemicalB((float(*) [columns][chemicals])cell_grid, chemical_b_to_add, x, y, chemical_b_radius);
            break;

        case EVENT_MOUSEMOVE: 
            if (is_drawing == true) 
            {
                addChemicalB((float(*) [columns][chemicals])cell_grid, chemical_b_to_add, x, y, chemical_b_radius);
            }
            break;

        case EVENT_LBUTTONUP:
            addChemicalB((float(*) [columns][chemicals])cell_grid, chemical_b_to_add, x, y, chemical_b_radius);
            is_drawing = false;
            break;
    }
}

// Since the cell grid itself isn't formatted as an image, I have to the convert the 0 to 1 values to 0 to 255 values to show it as an image.
void updateCellGridImage(Mat &cell_grid_image, float cell_grid[rows][columns][chemicals]) 
{
    for (int j = 0; j < rows; ++j) 
    {
        for (int i = 0; i < columns; ++i) 
        {
            cell_grid_image.at<uint8_t>(j, i) = (uint8_t)(fdimf(cell_grid[j][i][0], cell_grid[j][i][1]) * 255);
        }
    }
}

// I guess this just updates the cell grid
void updateCellGrid(float cell_grid[rows][columns][chemicals], float next_cell_grid[rows][columns][chemicals], float a_diffusion_rate, float b_diffusion_rate, float feed_rate, float kill_rate) 
{
    // This calculates the new cell grid values,

    for (int j = 0; j < rows; ++j) 
    {
        for (int i = 0; i < columns; ++i) 
        {
            float chemical_a_ammount = cell_grid[j][i][0];
            float chemical_b_ammount = cell_grid[j][i][1];

            // This is where the magic happens. 
            // These equations baiscally change the values of chemical A and chemical B
            next_cell_grid[j][i][0] = chemical_a_ammount + a_diffusion_rate * laplace(false, i, j, cell_grid) - chemical_a_ammount * chemical_b_ammount * chemical_b_ammount - chemical_a_ammount * feed_rate + feed_rate;
            next_cell_grid[j][i][1] = chemical_b_ammount + b_diffusion_rate * laplace(true, i, j, cell_grid) + chemical_a_ammount * chemical_b_ammount * chemical_b_ammount - chemical_b_ammount * feed_rate - chemical_b_ammount * kill_rate;
        }
    }

    // and this sets the old grid's values to the new grid's values
    for (int j = 0; j < rows; ++j) 
    {
        for (int i = 0; i < columns; ++i) 
        {
            cell_grid[j][i][0] = constrain(next_cell_grid[j][i][0], 0.0f, 1.0f);
            cell_grid[j][i][1] = constrain(next_cell_grid[j][i][1], 0.0f, 1.0f);
        }
    }
}

int main() 
{
    float cell_grid[rows][columns][chemicals] = {0};
    float next_cell_grid[rows][columns][chemicals] = {0};

    initializeCellGrid(cell_grid);
    initializeCellGrid(next_cell_grid);
    addChemicalB(cell_grid, 1.0f, columns*0.5f, rows*0.5f, 5);
    addChemicalB(next_cell_grid, 1.0f, columns*0.5f, rows*0.5f, 5);

    Mat cell_grid_image = Mat(rows, columns, CV_8U);

    std::string main_window_name = "Reaction Diffusion Simulation";
    std::string slider_window_name = "Sliders";
    namedWindow(main_window_name, WINDOW_NORMAL);
    namedWindow(slider_window_name, WINDOW_NORMAL);
    setMouseCallback(main_window_name, addChemicalBOnClick, (void*)cell_grid);

    // By default, these values are set to simulate coral growth.
    int feed_rate = 55;
    int kill_rate = 62;

    int window_size = 1;

    // Trackbar = Slider
    createTrackbar("Feed Rate", slider_window_name, &feed_rate, 100);
    createTrackbar("Kill Rate", slider_window_name, &kill_rate, 73);
    createTrackbar("Brush Size", slider_window_name, &brush_size, 100);
    createTrackbar("Win Size", slider_window_name, &window_size, 100);

    while (waitKey(1) < 0)
    {
        updateCellGrid(cell_grid, next_cell_grid, 1.0f, 0.5f, (float)feed_rate*0.001f, (float)kill_rate*0.001f);
        updateCellGridImage(cell_grid_image, cell_grid);
        imshow(main_window_name, cell_grid_image);
        resizeWindow(main_window_name, (columns + (1920-columns) * (window_size*0.01f)), ((rows) + (1080-rows) * (window_size*0.01f)));
    }
    
    destroyAllWindows();
}