Source: SurfaceGeneration.js

/**
 * @class
 * An object representing a particle. Particles are carried by the vector field and
 * their position is tracked to construct the surface.
 * @param {THREE.Vector3} pos The initial position of the particle.
 */
function Particle(pos) {
	/** The current position of the particle.*/
	this.position = pos.clone();
	/** The last position of the particle.*/
	this.lastPosition = null;
	
	/** The current normal. This is used, because normals need some kind of smoothing. */
	this.currentNormal = null;
	
	/** The left neighbor particle. */
	this.left = null;
	/** The right neighbor particle. */
	this.right = null;
	
	/** An counter of iterations this particle has survived. */
	this.lifetime = 0;
	
	/** 'true' if the particle has to be removed until the end of the iteration. */
	this.moribund = false;
};

/** 
 * Simple function checking if all values of the vector are inside the interval [0,1]. 
 * @param {THREE.Vector3} pos The vector to be checked.
 */
function isInVolume(pos) {
	return pos.x > 0.0 && pos.x < 1.0 &&
		pos.y > 0.0 && pos.y < 1.0 &&
		pos.z > 0.0 && pos.z < 1.0;
}

/**
 * @class
 * An object representing a surface in the vector field.
 */
function StreamSurface() {
	this.positions = null;
	this.colors = null;
	this.normals = null;
	this.velocities = null;
	
	this.volume = null;
	this.seedLine = null;
	this.particles = null;
	this.numIt = null;

	this.pushPoint = function(position, color, normal) {
		this.positions.push(position.clone());
		this.colors.push(color.clone());
		this.normals.push(normal.clone());
	}
	
	/** 
	 * Inserts a new particle between to neighbors. 
	 * @param {Particle} left The left neighbor.
	 * @param {Particle} right The right neighbor.
	 */
	this.insert = function(left, right) {
		var n = new Particle(
			new THREE.Vector3(
				(left.position.x + right.position.x)/2,
				(left.position.y + right.position.y)/2,
				(left.position.z + right.position.z)/2
			)
		);
		right.left = n;
		left.right = n;
		n.right = right;
		n.left = left;
		this.particles.push(n);
	}
	
	/**
	 * This method samples particle along the seed line.
	 */
	this.initParticles = function() {
		this.particles = [];

		var seedDirection = new THREE.Vector3();
		seedDirection.subVectors(this.seedLine.end, this.seedLine.start);
		var distanceLeft = seedDirection.length();
		seedDirection.normalize();
		seedDirection.multiplyScalar(this.seedLine.interval);

		var currentPoint = this.seedLine.start.clone();

		var previous = null;

		while (distanceLeft > 0) {
			var p = new Particle(currentPoint);
			if (previous) {
				previous.right = p;
				p.left = previous;
			}
			previous = p;	
			this.particles.push(p);
	
			currentPoint = currentPoint.add(seedDirection);
			distanceLeft = distanceLeft - this.seedLine.interval;
		}	
	}
	
	/**
	 * This method updates the position of all particles according to the vector field.
	 */
	this.updateParticles = function() {
		for (var j = this.particles.length - 1; j >= 0; j--) {
			var p = this.particles[j];
			
			p.lastPosition = p.position.clone();
			if (!isInVolume(p.position)) {
				p.moribund = true;
				continue;
			} 
			
			var direction = this.volume.trilinear(p.position).clone();
			p.position = p.position.add(direction.multiplyScalar(this.seedLine.interval));
			p.lifetime = p.lifetime + 1;
		}
	}
	
	/**
	 * This function calculates the normal of all particles and marks particles without
	 * neighbors as moribund.
	 */
	this.processParticles = function() {
		for (var j = this.particles.length - 1; j >= 0; j--) {
			var p = this.particles[j];
			
			if (p.moribund) {
				continue;
			}
			
			var b = new THREE.Vector3();
			if (p.left && p.right) {
				var b1 = new THREE.Vector3();
				var b2 = new THREE.Vector3();
				b1.subVectors(p.position, p.left.position).normalize();
				b2.subVectors(p.right.position, p.position).normalize();
				b.addVectors(b1, b2).divideScalar(2);
			} else if (p.left) {
				b = b.subVectors(p.position, p.left.position).normalize();
			} else if (p.right) {
				b = b.subVectors(p.right.position, p.position).normalize();
			} else {
				p.moribund = true;
			}
			
			if (!p.moribund) {
				var a = new THREE.Vector3();
				a.subVectors(p.position, p.lastPosition).normalize();
			
				p.currentNormal = new THREE.Vector3();
				p.currentNormal.crossVectors(a,b);
			}
		}
	}
	
	/**
	 * For all particles saves the position, normal and color to the buffers.
	 */
	this.pushParticles = function() {
		for (var j = this.particles.length - 1; j >= 0; j--) {
			var p = this.particles[j];
		
			if (p.moribund) {
				continue;
			}
			
			// Normals need to be smoothed, otherwise the result looks terrible.
			// Smoothing is done by taking the average of the normals of this particle,
			// the neighbors' and their neighbors'.
			var normal = new THREE.Vector3();
			normal.add(p.currentNormal);
			var count = 1;
			if (p.left && !p.left.moribund) {
				normal.add(p.left.currentNormal);
				count = count + 1;
				if (p.left.left && !p.left.left.moribund) {
					normal.add(p.left.left.currentNormal);
					count = count + 1;
				}
			} 
			if (p.right && !p.right.moribund) {
				normal.add(p.right.currentNormal);
				count = count + 1;
				if (p.right.right && !p.right.right.moribund) {
					normal.add(p.right.right.currentNormal);
					count = count + 1;
				}
			}
			normal.divideScalar(count);
		
			this.pushPoint(
				p.position.clone(),
				new THREE.Vector3(p.lifetime/this.numIt, 1.0 - p.lifetime/this.numIt, 0.0),
				normal
			);
		}
	}
	
	/**
	 * For all particles, checks if they need to be deleted, or if new particles need
	 * to be inserted. Deletes, if particles are to close to each other or marked as
	 * moribund. Inserts, if particles are to far from each other.
	 */
	this.checkDensity = function() {
		for (var j = this.particles.length - 1; j >= 0; j--) {
			var p = this.particles[j];
			
			if (p.moribund) {
				if (p.right) {
					if (p.left) {
						p.right.left = p.left;
						p.left.right = p.right;
					} else {
						p.right.left = null;
					}
				} else {
					if (p.left) {
						p.left.right = null;
					}
				}
				this.particles.splice(j, 1);
				continue;
			}
		
			// ignore most left and most right particle
			if (!p.left)
				continue;
			if (!p.right)
				continue;
			
			if (p.position.distanceTo(p.left.position) > 2 * this.seedLine.interval) {
				if (p.left.lastPosition) {
					if (p.position.distanceTo(p.left.position) - p.lastPosition.distanceTo(p.left.lastPosition) < 2 * p.position.distanceTo(p.lastPosition)) {
						this.insert(p.left, p);
					}
				} else {
					p.left.right = null;
					p.left = null;				
				}
			} else if (p.position.distanceTo(p.left.position) < 0.8 * this.seedLine.interval) {
				p.left = p.right;
				p.right = p.left;
				this.particles.splice(j, 1);
			}
		
			if (p.position.distanceTo(p.right.position) > 2 * this.seedLine.interval) {
				if (p.right.lastPosition) {
					if (p.position.distanceTo(p.right.position) - p.lastPosition.distanceTo(p.right.lastPosition) < 2 * p.position.distanceTo(p.lastPosition)) {
						this.insert(p, p.right);
					}
				} else {
					p.right.left = null;
					p.right = null;				
				}
			} else if (p.position.distanceTo(p.right.position) < 0.8 * this.seedLine.interval) {
				p.left = p.right;
				p.right = p.left;
				this.particles.splice(j, 1);
			}
		}
	}
	
	/**
	 * Runs the calculation of the surface.
	 * @param {VolumeData} volume The vector field to be used.
	 * @param {SeedLine} seedLine The seed line for sampling the initial particles.
	 * @param {int} numIt The number of iterations for the surface generation algorithm.
	 */
	this.calculate = function(volume, seedLine, numIt) {
		// reset values
		this.positions = [];
		this.colors = [];
		this.normals = [];
		this.velocities = [];
		
		// set up new params
		this.volume = volume;
		this.seedLine = seedLine;
		this.numIt = numIt;
		
		// generation
		this.initParticles(seedLine);
		for (var it = 0; it < numIt; it++) {		
			this.updateParticles();
			this.processParticles();
			this.pushParticles();
			this.checkDensity();
		}
	}
}