メタボール（Marching Squaresアルゴリズム）

/**
* metaballs (Marching Squares algorithm)
*
* @author aa_debdeb
* @date 2016/06/25
*/

int num = 10;
ArrayList<Particle> particles;
int cellNum = 50;
float cellSize = 10;
float[][] cells;

void setup(){
  size(500, 500);
  cells = new float[cellNum + 1][cellNum + 1];
  particles = new ArrayList<Particle>();
  for(int i = 0; i < num; i++){
    particles.add(new Particle());
  }
  noFill();
  stroke(255, 255, 240);
  strokeWeight(7);
}

void draw(){
  background(128, 128, 120);
  for(int y = 0; y < cellNum + 1; y++){
    for(int x = 0; x < cellNum + 1; x++){
      PVector c = new PVector(x * cellSize, y * cellSize);
      cells[x][y] = 0.0;
      for(Particle p: particles){
        cells[x][y] += p.radious / PVector.sub(c, p.loc).mag();
      }
    }
  }
  
  for(int y = 0; y < cellNum; y++){
    for(int x = 0; x < cellNum; x++){
      PVector c = new PVector(cellSize * x, cellSize * y);
      int state = int((cells[x][y + 1] >= 1 ? 1: 0)
                    + pow(cells[x + 1][y + 1] >= 1 ? 2: 0, 1)
                    + pow(cells[x + 1][y] >= 1 ? 2: 0, 2)
                    + pow(cells[x][y] >= 1 ? 2: 0, 3));
//      float halfSize = cellSize / 2.0;
//      switch(state){
//        case 0:
//          break;
//        case 1:
//          line(c.x, c.y + halfSize, c.x + halfSize, c.y + cellSize);
//          break;
//        case 2:
//          line(c.x + cellSize, c.y + halfSize, c.x + halfSize, c.y + cellSize);
//          break;
//        case 3:
//          line(c.x, c.y + halfSize, c.x + cellSize, c.y + halfSize);
//          break;
//        case 4:
//          line(c.x + halfSize, c.y, c.x + cellSize, c.y + halfSize);
//          break;
//        case 5:
//          line(c.x + halfSize, c.y, c.x, c.y + halfSize);
//          line(c.x + cellSize, c.y + halfSize, c.x + halfSize, c.y + cellSize);
//          break;
//        case 6:
//          line(c.x + halfSize, c.y, c.x + halfSize, c.y + cellSize);
//          break;
//        case 7:
//          line(c.x + halfSize, c.y, c.x, c.y + halfSize);
//          break;          
//        case 8:
//          line(c.x + halfSize, c.y, c.x, c.y + halfSize);
//          break;     
//        case 9:
//          line(c.x + halfSize, c.y, c.x + halfSize, c.y + cellSize);
//          break;
//        case 10:
//          line(c.x + halfSize, c.y, c.x + cellSize, c.y + halfSize);
//          line(c.x, c.y + halfSize, c.x + halfSize, cellSize);
//          break;
//        case 11:
//          line(c.x + halfSize, c.y, c.x + cellSize, c.y + halfSize);
//          break;
//        case 12:
//          line(c.x, c.y + halfSize, c.x + cellSize, c.y + halfSize);
//          break;
//        case 13:
//          line(c.x + cellSize, c.y + halfSize, c.x + halfSize, c.y + cellSize);
//          break;
//        case 14:
//          line(c.x, c.y + halfSize, c.x + halfSize, c.y + cellSize);
//          break;
//        case 15:
//          break;
//      }
      float x1, y1, x2, y2;
      switch(state){
        case 0:
          break;
        case 1:
        case 14:
          x1 = c.x;
          y1 = c.y + cellSize * ((1.0 - cells[x][y]) / (cells[x][y + 1] - cells[x][y]));
          x2 = c.x + cellSize * ((1.0 - cells[x][y + 1]) / (cells[x + 1][y + 1] - cells[x][y + 1]));
          y2 = c.y + cellSize;
          line(x1, y1, x2, y2);
          break;
        case 2:
        case 13:
          x1 = c.x + cellSize;
          y1 = c.y + cellSize * ((1.0 - cells[x + 1][y]) / (cells[x + 1][y + 1] - cells[x + 1][y]));
          x2 = c.x + cellSize * ((1.0 - cells[x][y + 1]) / (cells[x + 1][y + 1] - cells[x][y + 1]));
          y2 = c.y + cellSize;
          line(x1, y1, x2, y2);  
          break;
        case 3:
        case 12:
          x1 = c.x;
          y1 = c.y + cellSize * ((1.0 - cells[x][y]) / (cells[x][y + 1] - cells[x][y]));
          x2 = c.x + cellSize;
          y2 = c.y + cellSize * ((1.0 - cells[x + 1][y]) / (cells[x + 1][y + 1] - cells[x + 1][y]));
          line(x1, y1, x2, y2);  
          break;
        case 4:
        case 11:
          x1 = c.x + cellSize * ((1.0 - cells[x][y]) / (cells[x + 1][y] - cells[x][y]));
          y1 = c.y;
          x2 = c.x + cellSize;
          y2 = c.y + cellSize * ((1.0 - cells[x + 1][y]) / (cells[x + 1][y + 1] - cells[x + 1][y]));
          line(x1, y1, x2, y2);  
          break;
        case 5:
          x1 = c.x + cellSize * ((1.0 - cells[x][y]) / (cells[x + 1][y] - cells[x][y]));
          y1 = c.y;
          x2 = c.x;
          y2 = c.y + cellSize * ((1.0 - cells[x][y]) / (cells[x][y + 1] - cells[x][y]));
          line(x1, y1, x2, y2);  
          x1 = c.x + cellSize;
          y1 = c.y + cellSize * ((1.0 - cells[x + 1][y]) / (cells[x + 1][y + 1] - cells[x + 1][y]));
          x2 = c.x + cellSize * ((1.0 - cells[x][y + 1]) / (cells[x + 1][y + 1] - cells[x][y + 1]));
          y2 = c.y + cellSize;
          line(x1, y1, x2, y2);  
          break;
        case 6:
        case 9:
          x1 = c.x + cellSize * ((1.0 - cells[x][y]) / (cells[x + 1][y] - cells[x][y]));
          y1 = c.y;
          x2 = c.x + cellSize * ((1.0 - cells[x][y + 1]) / (cells[x + 1][y + 1] - cells[x][y + 1]));
          y2 = c.y + cellSize;
          line(x1, y1, x2, y2);  
          break;
        case 7:
        case 8:
          x1 = c.x + cellSize * ((1.0 - cells[x][y]) / (cells[x + 1][y] - cells[x][y]));
          y1 = c.y;
          x2 = c.x;
          y2 = c.y + cellSize * ((1.0 - cells[x][y]) / (cells[x][y + 1] - cells[x][y]));
          line(x1, y1, x2, y2);  
          break;    
        case 10:
          x1 = c.x + cellSize * ((1.0 - cells[x][y]) / (cells[x + 1][y] - cells[x][y]));
          y1 = c.y;
          x2 = c.x + cellSize;
          y2 = c.y + cellSize * ((1.0 - cells[x + 1][y]) / (cells[x + 1][y + 1] - cells[x + 1][y]));
          line(x1, y1, x2, y2);  
          x1 = c.x;
          y1 = c.y + cellSize * ((1.0 - cells[x][y]) / (cells[x][y + 1] - cells[x][y]));
          x2 = c.x + cellSize * ((1.0 - cells[x][y + 1]) / (cells[x + 1][y + 1] - cells[x][y + 1]));
          y2 = c.y + cellSize;
          line(x1, y1, x2, y2);
          break;
        case 15:
          break;      
      }
    }
  }
  for(Particle p: particles){
    p.update();
  }
}

class Particle{
  
  float radious;
  PVector loc, vel;
  
  Particle(){
    radious = random(10, 20);
    loc = new PVector(random(radious, width - radious), random(radious, height - radious));
    float velSize = random(2, 5);
    float velAng = random(TWO_PI);
    vel = new PVector(velSize * cos(velAng), velSize * sin(velAng));
  }
  
  void update(){
    loc.add(vel);
    if(loc.x < radious){
      vel.x *= -1;
      loc.x += vel.x;
    }
    if(loc.x >= width - radious){
      vel.x *= -1;
      loc.x += vel.x;
    }
    if(loc.y < radious){
      vel.y *= -1;
      loc.y += vel.y;
    }
    if(loc.y >= height - radious){
      vel.y *= -1;
      loc.y += vel.y;
    }
  }
}