VapoMotion

// options:

const fluidD = 3;

const mirror = false;

const maxTransparency = 0.8;

​

const fadeSpeed = 0.018f;

const deltaT = 1.2f;

const visc = 0.000001f;

​

const mainWidth = 640;

const mainHeight = 480;

//

​

var video = document.createElement("video");

video.setAttribute("style", "display:none;");

video.setAttribute("id", "videoOutput");

video.setAttribute("width", mainWidth + "px");

video.setAttribute("height", mainHeight + "px");

video.setAttribute("autoplay", "true");

if(document.body!=null) document.body.appendChild(video);

​

async function init() {

  stream = await navigator.mediaDevices.getUserMedia({video: true});

  video.srcObject = stream;

}

​

var ctx;

var currentImageData;

var lastImageData;

​

int nb=20;

​

//if the average of U and V for a zone is lower than this value, it is not rendered.

float UVCutoff = 4;

​

​

void setup()

{

    size(mainWidth,mainHeight);

    ctx = externals.context;

    ellipseMode(CORNER);

    smooth();

    init();

    invWidth = 1.0f/width;

    invHeight = 1.0f/height;

    aspectRatio = width * invHeight;

    aspectRatio2 = aspectRatio * aspectRatio;

​

    // create fluid and set options

    fluidSolver = new MSAFluidSolver2D((int)(FLUID_WIDTH), (int)(FLUID_WIDTH * height/width));

    fluidSolver.enableRGB(true).setFadeSpeed(fadeSpeed).setDeltaT(deltaT).setVisc(visc);

​

    // create image to hold fluid picture

    imgFluid = createImage(fluidSolver.getWidth(), fluidSolver.getHeight(), P2D);

    // create particle system

    //particleSystem = new ParticleSystem();  

}

​

void AdvanceBuffer()

{

  lastImageData = currentImageData;

  currentImageData = ctx.getImageData(0,0,mainWidth,mainHeight).data;

}

​

​

​

void draw()

{

    ctx.drawImage(video, 0, 0, mainWidth,mainHeight); //video is defined outside processing code

    AdvanceBuffer();

    ArrayList<FlowZone> fZones = calculateFlow(lastImageData, currentImageData, mainWidth,mainHeight);

    for(int i=0; i<fluidSolver.getNumCells(); i++) 

    {

      float d = fluidD;

      float r = fluidSolver.r[i]*d, g = fluidSolver.g[i]*d, b = fluidSolver.b[i]*d;

      float a = min(maxTransparency, r+g+b);//min(brightness(color(r,g,b)), 0.8);//saturation?

      imgFluid.pixels[i] = color(r,g,b,a);

    }           

    imgFluid.updatePixels();

    if(mirror){

        pushMatrix();

      scale(-1,1);

      translate(-mainWidth,0);

    }

    image(currentImageData, 0,0, mainWidth,mainHeight );

    image(imgFluid, 0,0, mainWidth,mainHeight );

    if(mirror) popMatrix();         

    for(FlowZone fz : fZones)

    {

      if(fz.ProductAboveLimit())

      {

        addForce(

          (fz.X * invWidth), 

          (fz.Y * invHeight), 

          (fz.U) * invWidth, 

          (fz.V) * invHeight

        );

      }

    }

    fluidSolver.update();

}

​

ArrayList<FlowZone> calculateFlow (oldImage, newImage,width,height) 

{

  ArrayList<FlowZone> zones = new ArrayList();

  if(oldImage == null)

    return zones;  

  var step = 8;

  var winStep = step * 2 + 1;

​

  var A2, A1B2, B1, C1, C2;

  var u, v, uu, vv;

  uu = vv = 0;

  var wMax = width - step - 1;

  var hMax = height - step - 1;

  var globalY, globalX, localY, localX;

​

  for (globalY = step + 1; globalY < hMax; globalY += winStep) 

  {

    for (globalX = step + 1; globalX < wMax; globalX += winStep) 

    {

      A2 = A1B2 = B1 = C1 = C2 = 0;

​

      for (localY = -step; localY <= step; localY++) 

      {

        for (localX = -step; localX <= step; localX++) 

        {

          var address = (globalY + localY) * width + globalX + localX;

​

          var gradX = (newImage[(address - 1) * 4]) - (newImage[(address + 1) * 4]);

          var gradY = (newImage[(address - width) * 4]) - (newImage[(address + width) * 4]);

          var gradT = (oldImage[address * 4]) - (newImage[address * 4]);

​

          A2 += gradX * gradX;

          A1B2 += gradX * gradY;

          B1 += gradY * gradY;

          C2 += gradX * gradT;

          C1 += gradY * gradT;

        }

      }

​

      var delta = (A1B2 * A1B2 - A2 * B1);

​

      if (delta !== 0) {

        /* system is not singular - solving by Kramer method */

        var Idelta = step / delta;

        var deltaX = -(C1 * A1B2 - C2 * B1);

        var deltaY = -(A1B2 * C2 - A2 * C1);

        u = deltaX * Idelta;

        v = deltaY * Idelta;

      }

      else

      {

        /* singular system - find optical flow in gradient direction */

        var norm = (A1B2 + A2) * (A1B2 + A2) + (B1 + A1B2) * (B1 + A1B2);

        if (norm !== 0) 

        {

          var IGradNorm = step / norm;

          var temp = -(C1 + C2) * IGradNorm;

          u = (A1B2 + A2) * temp;

          v = (B1 + A1B2) * temp;

        }

        else

        {

          u = v = 0;

        }

      }

​

      if (-winStep < u && u < winStep &&

          -winStep < v && v < winStep) {

        uu += u;

        vv += v;

        zones.add(new FlowZone(globalX, globalY, u, v));

      }

    }

  }

  return zones;

}

​

​

class FlowZone

{

  public float X, Y, U, V;

​

  public FlowZone(float x, float  y,float  u,float  v)

  {

    X = x;

    Y = y;

    U = u;

    V = v;

  }

​

  public void draw()

  {

    pushStyle();

    stroke(255);

    strokeWeight(2);

    //fill(0,0,0,80);

    line(X,Y,U*3, V*3);

    popStyle();

  }

  public bool ProductAboveLimit()

  {

    return ((abs(U) + abs(V)) / 2.0 > UVCutoff);

  }

}

​

//MSA fluid

final float FLUID_WIDTH = 256;

final static int maxParticles = 500;

​

float invWidth, invHeight;    // inverse of screen dimensions

float aspectRatio, aspectRatio2;

​

MSAFluidSolver2D fluidSolver;

ParticleSystem particleSystem;

​

PImage imgFluid; 

boolean untouched=true;

      // add force and dye to fluid, and create particles

      public void addForce(float x, float y, float dx, float dy) {

          float speed = dx * dx  + dy * dy * aspectRatio2;    // balance the x and y components of speed with the screen aspect ratio

          if(speed > 0) {

              if(x<0) x = 0;

              else if(x>1) x = 1;

              if(y<0) y = 0;

              else if(y>1) y = 1;

              float colorMult = 5;

              float velocityMult = 30.0f;

              int index = fluidSolver.getIndexForNormalizedPosition(x, y);

              int drawColor;

              colorMode(HSB, 360, 1, 1);

              float hue = ((x + y) * 180 + frameCount) % 360;

​

              drawColor = color(hue, 1, 1);

              //int offset = (y * FLUID_WIDTH + x ) * 4;

              //drawColor = color(currentImageData[offset+0], currentImageData[offset+1], currentImageData[offset+2]);

              //float h = hue(drawColor);

              //colorMode(HSB, 360, 1, 1);

              //  drawColor = color(h, 1, 1);

              colorMode(RGB, 1);

              fluidSolver.rOld[index]  += red(drawColor) * colorMult;

              fluidSolver.gOld[index]  += green(drawColor) * colorMult;

              fluidSolver.bOld[index]  += blue(drawColor) * colorMult;

            //i took out the particle system. sorry.

              //particleSystem.addParticles(x * width, y * height, 10);

              fluidSolver.uOld[index] += dx * velocityMult;

              fluidSolver.vOld[index] += dy * velocityMult;

          }

      }

      class Particle {

          final static float MOMENTUM = 0.5f;

          final static float FLUID_FORCE = 0.6f;

          float x, y;

          float vx, vy;

          float radius;       // particle's size

          float alpha;

          float mass;

          public void init(float x, float y) {

              this.x = x;

              this.y = y;

              vx = 0;

              vy = 0;

              radius = 5;

              alpha  = random(0.3f, 1);

              mass = random(0.1f, 1);

          }

          public void update() {

              // only update if particle is visible

              if(alpha == 0) return;

              // read fluid info and add to velocity

              int fluidIndex = fluidSolver.getIndexForNormalizedPosition(x * invWidth, y * invHeight);

              vx = fluidSolver.u[fluidIndex] * width * mass * FLUID_FORCE + vx * MOMENTUM;

              vy = fluidSolver.v[fluidIndex] * height * mass * FLUID_FORCE + vy * MOMENTUM;

              // update position

              x += vx;

              y += vy;

              // bounce of edges

              if(x<0) {

                  x = 0;

                  vx *= -1;

              }

              else if(x > width) {

                  x = width;

                  vx *= -1;

              }

              if(y<0) {

                  y = 0;

                  vy *= -1;

              }

              else if(y > height) {

                  y = height;

                  vy *= -1;

              }

              // hackish way to make particles glitter when the slow down a lot

              if(vx * vx + vy * vy < 1) {

                  vx = random(-1, 1);

                  vy = random(-1, 1);

              }

              // fade out a bit (and kill if alpha == 0);

              alpha *= 0.999f;

              if(alpha < 0.01f) alpha = 0;

          }

          public void drawOldSchool() {

              strokeWeight(alpha*1.5f);

   stroke(alpha, alpha, alpha,alpha);

              line(x-vx, y-vy,x, y);

          }

      }

      class ParticleSystem {

          int curIndex;

          Particle[] particles;

          ParticleSystem() {

              particles = new Particle[maxParticles];

              for(int i=0; i<maxParticles; i++) particles[i] = new Particle();

              curIndex = 0;

          }

          public void updateAndDraw(){

                  for(int i=0; i<maxParticles; i++) {

                      if(particles[i].alpha > 0) {

                          particles[i].update();

                          particles[i].drawOldSchool();    // use oldschool renderng

                      }

                  }

          }

          public void addParticles(float x, float y, int count ){

              for(int i=0; i<count; i++) addParticle(x + random(-15, 15), y + random(-15, 15));

          }

          public void addParticle(float x, float y) {

              particles[curIndex].init(x, y);

              curIndex++;

              if(curIndex >= maxParticles) curIndex = 0;

          }

      }

/***********************************************************************

 Copyright (c) 2008, 2009, Memo Akten, www.memo.tv

 *** The Mega Super Awesome Visuals Company ***

 * 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.

 *     * Neither the name of MSA Visuals nor the names of its contributors

 *       may be used to endorse or promote products derived from this software

 *       without specific prior written permission.

 *

 * 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.

 *

 * ***********************************************************************/

/**

 * this is a class for solving real-time fluid dynamics simulations based on Navier-Stokes equations

 * and code from Jos Stam's paper "Real-Time Fluid Dynamics for Games" http://www.dgp.toronto.edu/people/stam/reality/Research/pdf/GDC03.pdf

 * Other useful resources and implementations I looked at while building this lib:

 * Mike Ash (C) - http://mikeash.com/?page=pyblog/fluid-simulation-for-dummies.html

 * Alexander McKenzie (Java) - http://www.multires.caltech.edu/teaching/demos/java/stablefluids.htm

 * Pierluigi Pesenti (AS3 port of Alexander's) - http://blog.oaxoa.com/2008/01/21/actionscript-3-fluids-simulation/

 * Gustav Taxen (C) - http://www.nada.kth.se/~gustavt/fluids/

 * Dave Wallin (C++) - http://nuigroup.com/touchlib/ (uses portions from Gustav's)

 *

 *

 * @example MSAFluid

 * @author Memo Akten

 *

 */

public class MSAFluidSolver2D {

  public float[]  r;

  public float[]  g;

  public float[]  b;

  public float[]  u;

  public float[]  v;

  public float[]  rOld;

  public float[]  gOld;

  public float[]  bOld;

  public float[]  uOld;

  public float[]  vOld;

  public final String VERSION = "1.3.0";

  /**

   * return the version of the library.

   *

   * @return String

   */

  public String version() {

    return VERSION;

  }

  final static float    FLUID_DEFAULT_NX                    = 100;

  final static float    FLUID_DEFAULT_NY                    = 100;

  final static float    FLUID_DEFAULT_DT                    = 1.0f;

  final static float    FLUID_DEFAULT_VISC                  = 0.0001f;

  final static float    FLUID_DEFAULT_FADESPEED             = 0;

  final static int      FLUID_DEFAULT_SOLVER_ITERATIONS    = 10;

  /**

   * Constructor to initialize solver and setup number of cells

   * @param NX number of cells in X direction

   * @param NY number of cells in Y direction

   */

  public MSAFluidSolver2D(int NX, int NY) {

    r    = null;

    rOld = null;

    g    = null;

    gOld = null;

    b    = null;

    bOld = null;

    u    = null;

    uOld = null;

    v    = null;

    vOld = null;

    _isInited = false;

    setup(NX, NY);

  }

  /**

   * (OPTIONAL SETUP) re-initialize solver and setup number of cells

   * @param NX number of cells in X direction

   * @param NY number of cells in X direction

   * @return instance of MSAFluidSolver2D for further configuration

   */

  public MSAFluidSolver2D setup(int NX, int NY) {

    setDeltaT(FLUID_DEFAULT_DT);

    setFadeSpeed(FLUID_DEFAULT_FADESPEED  );

    setSolverIterations(FLUID_DEFAULT_SOLVER_ITERATIONS);

    _NX = NX;

    _NY = NY;

    _numCells = (_NX + 2) * (_NY + 2);

    _invNumCells = 1.0f/ _numCells;

    //    reset();

    _invNX = 1.0f / _NX;

    _invNY = 1.0f / _NY;

    width    = getWidth();

    height    = getHeight();

    invWidth  = 1.0f/width;

    invHeight  = 1.0f/height;

    reset();

    enableRGB(false);

    return this;

  }

  /**

   * (OPTIONAL SETUP) set timestep

   * @param dt timestep

   * @return instance of MSAFluidSolver2D for further configuration

   */

  public MSAFluidSolver2D setDeltaT(float dt) {

    _dt = dt;

    return this;

  }

  /**

   * (OPTIONAL SETUP) set how quickly the fluid dye dissipates and fades out

   * @param fadeSpeed (0...1)

   * @return instance of MSAFluidSolver2D for further configuration

   */

  public MSAFluidSolver2D setFadeSpeed(float fadeSpeed) {

    _fadeSpeed = fadeSpeed;

    return this;

  }

  /**

   * (OPTIONAL SETUP) set number of iterations for solver (higher is slower but more accurate)

   * @param solverIterations

   * @return instance of MSAFluidSolver2D for further configuration

   */

  public MSAFluidSolver2D setSolverIterations(int solverIterations) {

    _solverIterations = solverIterations;

    return this;

  }

  /**

   * (OPTIONAL SETUP) set whether solver should work with monochrome dye (default) or RGB

   * @param isRGB true or false

   * @return instance of MSAFluidSolver2D for further configuration

   */

  public MSAFluidSolver2D enableRGB(boolean isRGB) {

    _isRGB = isRGB;

    return this;

  }

  /**

   * (OPTIONAL SETUP) set viscosity

   * @param newVisc

   * @return instance of MSAFluidSolver2D for further configuration

   */

  public MSAFluidSolver2D setVisc(float newVisc) {

    visc = newVisc;

    return this;

  }

  /**

   * (OPTIONAL SETUP) randomize dye (useful for debugging)

   */

  public void randomizeColor() {

    for(int i=0; i< getWidth(); i++) {

      for(int j=0; j< getHeight(); j++) {

        int index = FLUID_IX(i, j);

        r[index] = rOld[index] = (float)Math.random();

        if(_isRGB) {

          g[index] = gOld[index] = (float)Math.random();

          b[index] = bOld[index] = (float)Math.random();

        }

      }

    }

  }

  /**

   * destroy solver and release all memory

   */

  public void destroy() {

    _isInited = false;

    r    = null;

    rOld = null;

    g    = null;

    gOld = null;

    b    = null;

    bOld = null;

    u    = null;

    uOld = null;

    v    = null;

    vOld = null;

  }

  /**

   * initialize solver (remove all velocities and dye)

  */

  public void reset() {

    destroy();

    _isInited = true;

    r    = new float[_numCells];

    rOld = new float[_numCells];

    g    = new float[_numCells];

    gOld = new float[_numCells];

    b    = new float[_numCells];

    bOld = new float[_numCells];

    u    = new float[_numCells];

    uOld = new float[_numCells];

    v    = new float[_numCells];

    vOld = new float[_numCells];

    for (int i = 0; i < _numCells; i++) {

      u[i] = uOld[i] = v[i] = vOld[i] = 0.0f;

      r[i] = rOld[i] = g[i] = gOld[i] = b[i] = bOld[i] = 0;

    }

  }

  /**

   * (INFO) get fluid cell index for (i,j) cell coordinates

   * @param i fluid cell index in x direction

   * @param j fluid cell index in y direction

   * @return cell index (to be used in r, g, b, u, v arrays)

   */

  public int getIndexForCellPosition(int i, int j) {

    if(i < 1) i=1; else if(i > _NX) i = _NX;

    if(j < 1) j=1; else if(j > _NY) j = _NY;

    return FLUID_IX(i, j);

  }

  /**

   * (INFO) get fluid cell index for normalized (x, y) coordinates

   * @param x 0...1 normalized x position

   * @param y 0...1 normalized y position

   * @return cell index (to be used in r, g, b, u, v arrays)

   */

  public int getIndexForNormalizedPosition(float x, float y) {

    return getIndexForCellPosition((int)Math.floor(x * (_NX+2)), (int)Math.floor(y * (_NY+2)));

  }

  /**

   * (INFO) whether the solver has been setup or not

   */

  public boolean isInited() {

    return _isInited;

  }

  /**

   * (INFO) return total number of cells (_NX+2) * (_NY+2)

  */

  public int getNumCells() {

    return _numCells;

  }

  /**

   * (INFO) return number of cells in x direction (_NX+2)

  */

  public int getWidth() {

    return _NX + 2;

  }

  /**

   * (INFO) return number of cells in y direction (_NY+2)

  */

  public int getHeight() {

    return _NY + 2;

  }

  /**

   * (INFO) return viscosity

  */

  public float getVisc() {

    return visc;

  }

  /**

   * (INFO) return average density of fluid

  */

  public float getAvgDensity() {

    return _avgDensity;

  }

  /**

   * (INFO) return average uniformity (distribution of densities and dye)

  */

  public float getUniformity() {

    return uniformity;

  }

  /**

   * (INFO) return average speed of fluid

  */

  public float getAvgSpeed() {

    return _avgSpeed;

  }

  /*

    public void getInfoAtPos(float x, float y, PVector vel, PVector col) {

      int i= (int)(x * (_NX+2));

      int j= (int)(y * (_NY+2));

      getInfoAtCell(i, j, vel, col);

    }

    public void getInfoAtCell(int i, int j,  PVector vel, PVector col) {

      if(i<0) i = 0; else if(i > _NX+1) i = _NX+1;

      if(j<0) j = 0; else if(j > _NY+1) j = _NY+1;

      getInfoAtCell(FLUID_IX(i, j), vel, col);

    }

    public void getInfoAtCell(int i,  PVector vel, PVector col) {

      //  if(safeToRun()){

      if(vel != null) vel.set(u[i] * _invNX, v[i] * _invNY, 0);

      if(col != null) {

        if(_isRGB) col.set(r[i], g[i], b[i]);

        else col.set(r[i], r[i], r[i]);

      }

      //    unlock();

      //  }

    }

    public void addForceAtPos(float x, float y, float vx, float vy) {

      int i = (int) (x * _NX + 1);

      int j = (int) (y * _NY + 1);

      if(i<0 || i>_NX+1 || j<0 || j>_NY+1) return;

      addForceAtCell(i, j, vx, vy);

    }

    public void addForceAtCell(int i, int j, float vx, float vy) {

      //  if(safeToRun()){

      int index = FLUID_IX(i, j);

      uOld[index] += vx * _NX;

      vOld[index] += vy * _NY;

      //    unlock();

      //  }

    }

    public void addColorAtPos(float x, float y, float r, float g, float b) {

      int i = (int) (x * _NX + 1);

      int j = (int) (y * _NY + 1);

      if(i<0 || i>_NX+1 || j<0 || j>_NY+1) return;

      addColorAtCell(i, j, r, g, b);

    }

    public void addColorAtCell(int i, int j, float r, float g, float b) {

      //  if(safeToRun()){

      int index = FLUID_IX(i, j);

      rOld[index] += r;

      if(_isRGB) {

        gOld[index] += g;

        bOld[index] += b;

      }

      //    unlock();

      //  }

    }

    public void randomizeColor() {

      for(int i=0; i< getWidth(); i++) {

        for(int j=0; j< getHeight(); j++) {

          int index = FLUID_IX(i, j);

          r[index] = rOld[index] = random(0, 1);

          if(_isRGB) {

            g[index] = gOld[index] = random(0, 1);

            b[index] = bOld[index] = random(0, 1);

          }

        }

      }

    }

  */

  /**

   * this must be called once every frame to move the solver one step forward

   * i.e. in your sketch draw() method

  */

  public void update() {

//    ADD_SOURCE_UV();

    addSourceUV();

    swapU();

    swapV();

//    DIFFUSE_UV();

    diffuseUV(0, visc);

    project(u, v, uOld, vOld);

    swapU();

    swapV();

    advect(1, u, uOld, uOld, vOld);

    advect(2, v, vOld, uOld, vOld);

    project(u, v, uOld, vOld);

    if(_isRGB) {

      //ADD_SOURCE_RGB();

      addSourceRGB();

      swapRGB();

      //DIFFUSE_RGB();

      diffuseRGB(0, 0);

      swapRGB();

      //ADVECT_RGB();

      advectRGB(0, u, v);

      fadeRGB();

    } else {

      addSource(r, rOld);

      swapR();

      diffuse(0, r, rOld, 0);

      swapRGB();

      advect(0, r, rOld, u, v);

      fadeR();

    }

  }

  protected void fadeR() {

    // I want the fluid to gradually fade out so the screen doesn't fill. the amount it fades out depends on how full it is, and how uniform (i.e. boring) the fluid is...

//    float holdAmount = 1 - _avgDensity * _avgDensity * _fadeSpeed;  // this is how fast the density will decay depending on how full the screen currently is

    float holdAmount = 1 - _fadeSpeed;

    _avgDensity = 0;

    _avgSpeed = 0;

    float totalDeviations = 0;

    float currentDeviation;

    //  float uniformityMult = uniformity * 0.05f;

    _avgSpeed = 0;

    for (int i = 0; i < _numCells; i++) {

      // clear old values

      uOld[i] = vOld[i] = 0;

      rOld[i] = 0;

      //    gOld[i] = bOld[i] = 0;

      // calc avg speed

      _avgSpeed += u[i] * u[i] + v[i] * v[i];

      // calc avg density

      r[i] = Math.min(1.0f, r[i]);

      //    g[i] = Math.min(1.0f, g[i]);

      //    b[i] = Math.min(1.0f, b[i]);

      //    float density = Math.max(r[i], Math.max(g[i], b[i]));

      float density = r[i];

      _avgDensity += density;  // add it up

      // calc deviation (for uniformity)

      currentDeviation = density - _avgDensity;

      totalDeviations += currentDeviation * currentDeviation;

      // fade out old

      r[i] *= holdAmount;

    }

    _avgDensity *= _invNumCells;

    //  _avgSpeed *= _invNumCells;

    //  println("%.3f\n", _avgSpeed);

    uniformity = 1.0f / (1 + totalDeviations * _invNumCells);    // 0: very wide distribution, 1: very uniform

  }

  protected void fadeRGB() {

    // I want the fluid to gradually fade out so the screen doesn't fill. the amount it fades out depends on how full it is, and how uniform (i.e. boring) the fluid is...

//    float holdAmount = 1 - _avgDensity * _avgDensity * _fadeSpeed;  // this is how fast the density will decay depending on how full the screen currently is

    float holdAmount = 1 - _fadeSpeed;

    _avgDensity = 0;

    _avgSpeed = 0;

    float totalDeviations = 0;

    float currentDeviation;

    //  float uniformityMult = uniformity * 0.05f;

    _avgSpeed = 0;

    for (int i = 0; i < _numCells; i++) {

      // clear old values

      uOld[i] = vOld[i] = 0;

      rOld[i] = 0;

      gOld[i] = bOld[i] = 0;

      // calc avg speed

      _avgSpeed += u[i] * u[i] + v[i] * v[i];

      // calc avg density

      r[i] = Math.min(1.0f, r[i]);

      g[i] = Math.min(1.0f, g[i]);

      b[i] = Math.min(1.0f, b[i]);

      float density = Math.max(r[i], Math.max(g[i], b[i]));

      //float density = r[i];

      _avgDensity += density;  // add it up

      // calc deviation (for uniformity)

      currentDeviation = density - _avgDensity;

      totalDeviations += currentDeviation * currentDeviation;

      // fade out old

      r[i] *= holdAmount;

      g[i] *= holdAmount;

      b[i] *= holdAmount;

    }

    _avgDensity *= _invNumCells;

    _avgSpeed *= _invNumCells;

    //println("%.3f\n", _avgDensity);

    uniformity = 1.0f / (1 + totalDeviations * _invNumCells);    // 0: very wide distribution, 1: very uniform

  }

  protected void addSourceUV() {

    for (int i = 0; i < _numCells; i++) {

      u[i] += _dt * uOld[i];

      v[i] += _dt * vOld[i];

    }

  }

  protected void addSourceRGB() {

    for (int i = 0; i < _numCells; i++) {

      r[i] += _dt * rOld[i];

      g[i] += _dt * gOld[i];

      b[i] += _dt * bOld[i];

    }

  }

  protected void addSource(float[] x, float[] x0) {

    for (int i = 0; i < _numCells; i++) {

      x[i] += _dt * x0[i];

    }

  }

  protected void advect(int b, float[] _d, float[] d0, float[] du, float[] dv) {

    int i0, j0, i1, j1;

    float x, y, s0, t0, s1, t1, dt0;

    dt0 = _dt * _NX;

    for (int i = 1; i <= _NX; i++) {

      for (int j = 1; j <= _NY; j++) {

        x = i - dt0 * du[FLUID_IX(i, j)];

        y = j - dt0 * dv[FLUID_IX(i, j)];

        if (x > _NX + 0.5) x = _NX + 0.5f;

        if (x < 0.5)     x = 0.5f;