Reaction diffusion

Gray scott algorithm for reaction diffusion patterns, then using marching squares to trace blobs to paths.

More info on kill ,feed parameters mrob.com/pub/comp/xmorphia/index.html

See sarah.skyon.be/conte…actiondiffusion.html for the iterative process

Log in to post a comment.

// Global code will be evaluated once.

const scale = 50; //min = 0, max = 100, step = 10
const tileWidth = 3;
const tileHeight = 3;
const tileOutsideArea = true; // set to true to make sure the outside area is still drawn on the edges
const randomSeed = 12345; //min = 1, max = 999999, step = 1


const killRatePercentage = 6.3; //min = 1, max = 10, step = 0.1
const feedRatePercentage = 5; //min = 1, max = 10, step = 0.1

const f = feedRatePercentage / 100;
const k = killRatePercentage / 100;

// the number of starting points to introduce the other gas to diffuse
const numberOfInitialSeeds = 150; //min = 1, max = 1000, step = 1

// the area of the simulation, bigger means slower.
const grayScottSize = 200; //min = 100, max = 500, step = 10
 // the more steps, the longer the simulation will run
const steps = 2000; //min = 100, max = 10000, step = 100

const originOffset = 256/4;

// You can find the Turtle API reference here: https://turtletoy.net/syntax
Canvas.setpenopacity(1);

class LineSegment {
    constructor(x1, y1, x2, y2) {
        this.x1 = x1;
        this.y1 = y1;
        this.x2 = x2;
        this.y2 = y2;
    }
}
class Random {
    constructor(seed) {
        if (typeof seed === "undefined")
            seed = ~~(Math.random() * 10000000);
        this.seed = seed;
        this.val = seed;
    }
    next() {
        // this is in no way uniformly distributed, so it's really a bad rng, but it's fast enough
        // and random enough
        let x = Math.sin(this.val++) * 10000;
        return x - Math.floor(x);
    }
    nextBetween(min, max) {
        return min + this.next() * (max - min);
    }
}
class Cell {
    constructor() {
        this.a = 0;
        this.b = 0;
    }
}
class GrayScott {
    constructor(width, height, k, f, DA = 1, DB = 0.5) {
        this.width = width;
        this.height = height;
        this.k = k;
        this.f = f;
        this.DA = DA;
        this.DB = DB;
        this.dt = 1;
        this._t = 0;
        this.minB = Number.POSITIVE_INFINITY;
        this.maxB = Number.NEGATIVE_INFINITY;
        this.minA = Number.POSITIVE_INFINITY;
        this.maxA = Number.NEGATIVE_INFINITY;
        this.buffer = this.createArray();
        this.nextFrameBuffer = this.createArray();
    }
    get t() {
        return this._t;
    }
    createArray() {
        let arr = new Array(this.width);
        for (let i = 0; i < this.width; i++) {
            arr[i] = new Array(this.height);
            for (let j = 0; j < this.height; j++) {
                let c = new Cell();
                c.a = 1;
                c.b = 0;
                arr[i][j] = c;
            }
        }
        return arr;
    }
    addInitialSeed(x, y, rng) {
        for (let i = 0; i < 2; i++) {
            for (let j = 0; j < 2; j++) {
                this.buffer[Math.floor(x + i)][Math.floor(y + j)].a = 1;
                this.buffer[Math.floor(x + i)][Math.floor(y + j)].b = rng.next();
            }
        }
    }
    step() {
        this.minA = Number.POSITIVE_INFINITY;
        this.minB = Number.POSITIVE_INFINITY;
        this.maxA = Number.NEGATIVE_INFINITY;
        this.maxB = Number.NEGATIVE_INFINITY;
        for (let x = 0; x < this.width; x++) {
            for (let y = 0; y < this.height; y++) {
                let oldA = this.buffer[x][y].a;
                let oldB = this.buffer[x][y].b;
                let middle = 0;
                let adjacent = 0;
                let diag = 0;
                middle = oldA * -1;
                let left = x - 1;
                let right = x + 1;
                let top = y - 1;
                let bottom = y + 1;
                if (left < 0)
                    left = this.width - 1;
                if (right >= this.width)
                    right %= this.width;
                if (top < 0)
                    top = this.height - 1;
                if (bottom >= this.height)
                    bottom %= this.height;
                adjacent += this.buffer[x][bottom].a * 0.2;
                adjacent += this.buffer[x][top].a * 0.2;
                adjacent += this.buffer[left][y].a * 0.2;
                adjacent += this.buffer[right][y].a * 0.2;
                diag += this.buffer[right][bottom].a * 0.05;
                diag += this.buffer[right][top].a * 0.05;
                diag += this.buffer[left][top].a * 0.05;
                diag += this.buffer[left][bottom].a * 0.05;
                let factor2 = oldA * oldB * oldB;
                let laplaceA = middle + adjacent + diag;
                let newA = (oldA + (this.DA * laplaceA - factor2) + this.f * (1 - oldA)) * this.dt;
                middle = oldB * -1;
                adjacent = 0;
                diag = 0;
                adjacent += this.buffer[x][bottom].b * 0.2;
                adjacent += this.buffer[x][top].b * 0.2;
                adjacent += this.buffer[left][y].b * 0.2;
                adjacent += this.buffer[right][y].b * 0.2;
                diag += this.buffer[right][bottom].b * 0.05;
                diag += this.buffer[right][top].b * 0.05;
                diag += this.buffer[left][top].b * 0.05;
                diag += this.buffer[left][bottom].b * 0.05;
                let laplaceB = middle + adjacent + diag;
                let newB = (oldB + (this.DB * laplaceB + factor2) - (this.k + this.f) * oldB) * this.dt;
                if (newB < this.minB)
                    this.minB = newB;
                if (newA < this.minA)
                    this.minA = newA;
                if (newB > this.maxB)
                    this.maxB = newB;
                if (newA > this.maxA)
                    this.maxA = newA;
                this.nextFrameBuffer[x][y].a = newA;
                this.nextFrameBuffer[x][y].b = newB;
                // A' = A + (lapl - A*B²  + F * (1-A))* dt
                // B' = B + (lapl + A*B²  - (k+F) * B)* dt
            }
        }
        this.swapBuffers();
        this._t++;
    }
    swapBuffers() {
        let tmp = this.buffer;
        this.buffer = this.nextFrameBuffer;
        this.nextFrameBuffer = tmp;
    }
    foreachCell(func) {
        for (let x = 0; x < this.width; x++) {
            for (let y = 0; y < this.height; y++) {
                let c = this.buffer[x][y];
                let val = 0;
                // normalize
                //val += Math.floor((c.b - this.minB) / (this.maxB - this.minB) * 255);
                val += this.maxA - this.minA > 0 ? Math.floor(((c.a - this.minA) / (this.maxA - this.minA)) * 255) : 0;
                func(x, y, val);
            }
        }
    }
}
class Point {
    constructor(x, y) {
        this.x = x;
        this.y = y;
    }
}
class MarchingSquareCell {
    constructor(nr, p1, p2) {
        this.nr = nr;
        this.p1 = p1;
        this.p2 = p2;
    }
}
/**
 *
 * https://en.wikipedia.org/wiki/Marching_squares
 */
function marchingSquares(values, threshold, width, height, wrapAround = true) {
    let mask;
    let cells = [];
    for (let j = 0; j < height; j++)
        cells.push([]);
    for (let j = 0; j < height; j++) {
        for (let i = 0; i < width; i++) {
            if (!wrapAround) {
                mask = [
                    j + 1 >= height ? values[j][i] > threshold : values[j + 1][i] > threshold,
                    i + 1 >= width || j + 1 >= height ? values[j][i] > threshold : values[j + 1][i + 1] > threshold,
                    i + 1 >= width ? values[j][i] > threshold : values[j][i + 1] > threshold,
                    values[j][i] > threshold,
                ];
            }
            else {
                let nextX = i + 1 >= width ? i + 1 - width : i + 1;
                let nextY = j + 1 >= height ? j + 1 - height : j + 1;
                mask = [
                    values[nextY][i] > threshold,
                    values[nextY][nextX] > threshold,
                    values[j][nextX] > threshold,
                    values[j][i] > threshold,
                ];
            }
            let nr = ((mask[0] ? 1 : 0) << 3) + ((mask[1] ? 1 : 0) << 2) + ((mask[2] ? 1 : 0) << 1) + ((mask[3] ? 1 : 0) << 0);
            if (nr != 0 && nr != 15) {
                let uv;
                if (!wrapAround) {
                    uv = getUV(nr, threshold, values[j][i], j + 1 >= height ? threshold : values[j + 1][i], i + 1 >= width ? threshold : values[j][i + 1], i + 1 >= width || j + 1 >= height ? threshold : values[j + 1][i + 1]);
                }
                else {
                    let nextX = i + 1 >= width ? i + 1 - width : i + 1;
                    let nextY = j + 1 >= height ? j + 1 - height : j + 1;
                    uv = getUV(nr, threshold, values[j][i], values[nextY][i], values[j][nextX], values[nextY][nextX]);
                }
                if (uv.length > 0) {
                    cells[j].push(new MarchingSquareCell(nr, uv[0], uv[1]));
                }
            }
            else
                cells[j].push(new MarchingSquareCell(nr));
        }
    }
    return cells;
}
function getUV(nr, threshold, vLT = 0, vLB = 2, vRT = 0, vRB = 2) {
    // (1 - min) / (max-min)
    let lInterpol = (threshold - vLT) / (vLB - vLT);
    let rInterpol = (threshold - vRT) / (vRB - vRT);
    let tInterpol = (threshold - vLT) / (vRT - vLT);
    let bInterpol = (threshold - vLB) / (vRB - vLB);
    let l = new Point(0, lInterpol);
    let t = new Point(tInterpol, 0);
    let r = new Point(1, rInterpol);
    let b = new Point(bInterpol, 1);
    switch (nr) {
        case 0:
            return [];
        case 1:
            return [l, t];
        case 2:
            return [t, r];
        case 3:
            return [l, r];
        case 4:
            return [r, b];
        case 5:
            return [b, l];
        case 6:
            return [b, t];
        case 7:
            return [b, l];
        case 8:
            return [b, l];
        case 9:
            return [b, t];
        case 10:
            return [b, r];
        case 11:
            return [b, r];
        case 12:
            return [l, r];
        case 13:
            return [t, r];
        case 14:
            return [t, l];
        case 15:
            return [];
        default:
            return [];
    }
}





function getMarchingSquaresResultFromReactionDiffusion(randomSeed, nrOfSeeds, grayScottSize, k, f, steps) {
    const gs = new GrayScott(grayScottSize, grayScottSize, k, f);
    const rng = new Random(randomSeed);
    for (let k = 0; k < nrOfSeeds; k++) {
        let x = Math.floor(rng.next() * (grayScottSize - 2));
        let y = Math.floor(rng.next() * (grayScottSize - 2));
        gs.addInitialSeed(x, y, rng);
    }
    console.log("running gs");
    for (let i = 0; i < steps; i++) {
        gs.step();
    }
    const cells = [];
    for (let i = 0; i < grayScottSize; i++) {
        cells.push([]);
        for (let j = 0; j < grayScottSize; j++) {
            cells[i].push(0);
        }
    }
    gs.foreachCell((i, j, val) => {
        cells[i][j] = 255 - val;
    });
    console.log("done running gs");
    console.log("running marching squares");
    let result = marchingSquares(cells, 128, grayScottSize, grayScottSize);
    return result;
}
function getReactionDiffusionSegments(scale, randomSeed = 12345, numberOfSeeds = 150, grayScottSize = 250, k = 0.063, f = 0.05, steps = 2000) {
    const cellSize = scale / grayScottSize;
    const result = getMarchingSquaresResultFromReactionDiffusion(randomSeed, numberOfSeeds, grayScottSize, k, f, steps);
    let segments = [];
    for (let j = 0; j < grayScottSize; j++) {
        for (let i = 0; i < grayScottSize; i++) {
            const cell = result[j][i];
            if (cell.nr != 0 && cell.nr != 15) {
                let xOffset = i * cellSize;
                let yOffset = j * cellSize;
                let x1 = xOffset + cellSize * cell.p1.x + cellSize / 2;
                let y1 = yOffset + cellSize * cell.p1.y + cellSize / 2;
                let x2 = xOffset + cellSize * cell.p2.x + cellSize / 2;
                let y2 = yOffset + cellSize * cell.p2.y + cellSize / 2;
                segments.push(new LineSegment(x1, y1, x2, y2));
            }
        }
    }
    return segments;
}
function getReactionDiffusionPath(scale, randomSeed = 12345, numberOfSeeds = 150, grayScottSize = 250, k = 0.060, f = 0.06, steps = 2000) {
    const cellSize = scale / grayScottSize;
    const result = getMarchingSquaresResultFromReactionDiffusion(randomSeed, numberOfSeeds, grayScottSize, k, f, steps);
    console.log("rendering marching squares");
    const polygons = getPathsFromMarchingSquaresResult(result, grayScottSize, grayScottSize, cellSize);
    return polygons;
}

function getPathsFromMarchingSquaresResult(result, width, height, cellSize) {
    let polygons = [];
    const visited = [];
    for (let j = 0; j < height; j++) {
        visited.push([]);
        for (let i = 0; i < width; i++) {
            visited[j].push(false);
        }
    }
    //  let debugCanvas: HTMLCanvasElement = document.createElement("canvas");
    //  debugCanvas.width = width;
    //  debugCanvas.height = height;
    //  document.getElementById("container").parentNode.appendChild(debugCanvas);
    //  let debugCtx = debugCanvas.getContext("2d");
    const log = false;
    for (let j = 0; j < height; j++) {
        for (let i = 0; i < width; i++) {
            if (!visited[j][i]) {
                let nr = result[j][i].nr;
                if (nr != 0 && nr != 15) {
                    let stepX = 0;
                    let stepY = 0;
                    let subSegment = [];
                    if (log)
                        console.log("new subsegment at " + (i + "," + j));
                    let curY = j + stepY;
                    let curX = i + stepX;
                    while (!visited[curY][curX]) {
                        let oldNr = nr;
                        nr = result[curY][curX].nr;
                        if (nr != 0 && nr != 15) {
                            let cellP1 = result[curY][curX].p1;
                            let cellP2 = result[curY][curX].p2;
                            let xOffset = (i + stepX) * cellSize;
                            let yOffset = (j + stepY) * cellSize;
                            let x1 = xOffset + cellSize * cellP1.x + cellSize / 2;
                            let y1 = yOffset + cellSize * cellP1.y + cellSize / 2;
                            let x2 = xOffset + cellSize * cellP2.x + cellSize / 2;
                            let y2 = yOffset + cellSize * cellP2.y + cellSize / 2;
                            if (oldNr == nr || (oldNr != nr && nr != 5 && nr != 10)) {
                                // if we wrapped around don't flag it as visited
                                visited[curY][curX] = true;
                            }
                            // debugCtx.fillStyle = `rgb(255, ${nr}, 255)`;
                            // debugCtx.fillRect(i + stepX, j + stepY, 1, 1);
                            if (subSegment.length == 0) {
                                subSegment.push({ x: x1, y: y1 });
                            }
                            switch (nr) {
                                case 4:
                                case 12:
                                case 13:
                                    //right
                                    //     console.log(nr + " going right");
                                    //if (i + stepX + 1 < width)
                                    stepX++;
                                    if (x1 < x2)
                                        subSegment.push({ x: x2, y: y2 });
                                    else
                                        subSegment.push({ x: x1, y: y1 });
                                    break;
                                case 5:
                                    if (log)
                                        console.log(nr + " saddle point at " + (i + stepX) + ", " + (j + stepY));
                                    if (oldNr == 2 || oldNr == 6 || oldNr == 14) {
                                        // came from the bottom -> to right
                                        if (log)
                                            console.log("came from bottom, go to right");
                                        //if (i + stepX + 1 < width)
                                        stepX++;
                                        if (x1 < x2)
                                            subSegment.push({ x: x2, y: y2 });
                                        else
                                            subSegment.push({ x: x1, y: y1 });
                                    }
                                    else {
                                        // came from the top -> to left
                                        if (log)
                                            console.log("came from top, go to left");
                                        //if (i + stepX - 1 >= 0)
                                        stepX--;
                                        if (x2 < x1)
                                            subSegment.push({ x: x2, y: y2 });
                                        else
                                            subSegment.push({ x: x1, y: y1 });
                                    }
                                    break;
                                case 10:
                                    if (log)
                                        console.log(nr + " saddle point at " + (i + stepX) + ", " + (j + stepY));
                                    if (oldNr == 1 || oldNr == 3 || oldNr == 7) {
                                        // came from the right -> to right
                                        if (log)
                                            console.log("came from bottom, go up");
                                        //if (j + stepY - 1 >= 0)
                                        stepY--;
                                        if (y2 < y1)
                                            subSegment.push({ x: x2, y: y2 });
                                        else
                                            subSegment.push({ x: x1, y: y1 });
                                    }
                                    else {
                                        // came from the left -> to bottom
                                        if (log)
                                            console.log("came from left, go down");
                                        //if (j + stepY + 1 < height)
                                        stepY++;
                                        if (y1 < y2)
                                            subSegment.push({ x: x2, y: y2 });
                                        else
                                            subSegment.push({ x: x1, y: y1 });
                                        break;
                                    }
                                    break;
                                case 2:
                                case 6:
                                case 14:
                                    //up
                                    if (log)
                                        console.log(nr + " going up");
                                    //if (j + stepY - 1 >= 0)
                                    stepY--;
                                    if (y2 < y1)
                                        subSegment.push({ x: x2, y: y2 });
                                    else
                                        subSegment.push({ x: x1, y: y1 });
                                    break;
                                case 9:
                                case 11:
                                case 8:
                                    //down
                                    if (log)
                                        console.log(nr + " going down");
                                    //if (j + stepY + 1 < height)
                                    stepY++;
                                    if (y1 < y2)
                                        subSegment.push({ x: x2, y: y2 });
                                    else
                                        subSegment.push({ x: x1, y: y1 });
                                    break;
                                case 1:
                                case 3:
                                case 7:
                                    // left
                                    if (log)
                                        console.log(nr + " going left");
                                    //if (i + stepX - 1 >= 0)
                                    stepX--;
                                    if (x2 < x1)
                                        subSegment.push({ x: x2, y: y2 });
                                    else
                                        subSegment.push({ x: x1, y: y1 });
                                    break;
                                default:
                            }
                        }
                        else {
                            break;
                        }
                        curY = j + stepY;
                        curX = i + stepX;
                        if (curX < 0)
                            curX += width; // wrap around
                        if (curY < 0)
                            curY += height;
                        if (curX >= width)
                            curX -= width;
                        if (curY >= height)
                            curY -= height;
                    }
                    if (subSegment.length > 0) {
                        if (polygons.length > 0) {
                            if (subSegment.every((p) => inside(p, polygons[polygons.length - 1].segments[0]))) {
                                polygons[polygons.length - 1].segments.push(subSegment);
                            }
                            else {
                                polygons.push({
                                    segments: [subSegment],
                                });
                            }
                        }
                        else
                            polygons.push({
                                segments: [subSegment],
                            });
                    }
                }
            }
        }
    }
    return polygons;
}
function inside(point, vs) {
    // ray-casting algorithm based on
    // https://wrf.ecse.rpi.edu/Research/Short_Notes/pnpoly.html/pnpoly.html
    var x = point.x, y = point.y;
    var inside = false;
    for (var i = 0, j = vs.length - 1; i < vs.length; j = i++) {
        var xi = vs[i].x, yi = vs[i].y;
        var xj = vs[j].x, yj = vs[j].y;
        var intersect = yi > y != yj > y && x < ((xj - xi) * (y - yi)) / (yj - yi) + xi;
        if (intersect)
            inside = !inside;
    }
    return inside;
}




const reactionDiffusionPaths = getReactionDiffusionPath(scale, randomSeed, numberOfInitialSeeds, grayScottSize, k, f, steps);

const turtle = new Turtle();


// The walk function will be called until it returns false.
const iterator = draw();

function walk(i) {
    const val = iterator.next();
    return !val.done;
}

function* draw() {

      
  for (let j = tileOutsideArea ? -1 : 0; j < tileWidth + (tileOutsideArea ? 1 : 0); j++) {
    for (let i = tileOutsideArea ? -1 : 0; i < tileHeight + (tileOutsideArea ? 1 : 0); i++) {
      let offsetX = -originOffset + i * scale;
      let offsetY = -originOffset + j * scale;

      for (let polygon of reactionDiffusionPaths) {
          
        
        turtle.penup();

        for (let segment of polygon.segments) {
          turtle.goto(offsetX + segment[0].x, offsetY + segment[0].y);
          turtle.pendown();
          for (let i = 1; i < segment.length; i++) {
            //ctx.lineTo(offsetX + segment[i].x, offsetY + segment[i].y);
            turtle.goto(offsetX + segment[i].x, offsetY + segment[i].y);
          }
          // close path
           turtle.goto(offsetX + segment[0].x, offsetY + segment[0].y);
        }
        
        if (i >= 0 && j >= 0) {
          yield;
        }
      }
    }
  }
}