examples/region_based_shape_matching.html

<!doctype html>
<html class="no-js" lang="en">
<head>
	<meta charset="utf-8">
	<style>
		body {font-family: Helvetica, sans-serif;}
		table {background-color:#CCDDEE;text-align:left}
	</style>
	<script type="text/x-mathjax-config">
	  MathJax.Hub.Config({
		extensions: ["tex2jax.js"],
		jax: ["input/TeX", "output/HTML-CSS"],
		tex2jax: {
		  inlineMath: [ ['$','$'], ["\\(","\\)"] ],
		  displayMath: [ ['$$','$$'], ["\\[","\\]"] ],
		  processEscapes: true
		},
		"HTML-CSS": { fonts: ["TeX"] }
	  });
	</script>
	<script type="text/javascript" aync src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.4/MathJax.js"></script>
	<title>Region-Based Shape Matching</title>
</head>
<body>
<main>
	<h1 style="text-align:center">Region-Based Shape Matching</h1>
	<table style="align_center;border-radius: 20px;padding: 20px;margin:auto">
		<col width="1100"> 
		<col width="400"> 
		<tr>
			<td>
				<canvas id="simCanvas" width="1024" height="768" style="border:2px solid #000000;border-radius: 20px;background-color:#EEEEEE">Your browser does not support the HTML5 canvas tag.</canvas>
			</td>
			<td>
				<table>		
				<col width="180" style="padding-right:10px"> 
				<col width="100"> 
				<tr>
					<td><label>Current time</label></td>
					<td><span id="time">0.00</span> s</td>
				</tr>
				<tr>
					<td><label>Time per sim. step</label></td>
					<td><span id="timePerStep">0.00</span> ms</td>
				</tr>
				<tr>
					<td><label># particles</label></td>
					<td><span id="numParticles">0</span></td>
				</tr>
				<tr>
					<td><label for="widthInput">Width</label></td>
					<td><input onchange="gui.restart()" id="widthInput" type="number" value="40" step="1"></td>
				</tr>
				<tr>
					<td><label for="heightInput">Height</label></td>
					<td><input onchange="gui.restart()" id="heightInput" type="number" value="30" step="1"></td>
				</tr>
				<tr>
					<td><label for="fixedParticlesInput"># fixed particles</label></td>
					<td><select onchange="gui.restart()" id="fixedParticlesInput">
					  <option value="1">1</option>
					  <option value="2" selected="selected">2</option>
					  <option value="4">4</option>
					</select></td>
				</tr>
				<tr>
					<td><label for="regionHalfWidthInput">Region halfwidth</label></td>
					<td><input onchange="gui.restart()" id="regionHalfWidthInput" type="number" step="1" min="1" value="1"></td>
				</tr>
				<tr>
					<td><label for="timeStepSizeInput">Time step size</label></td>
					<td><input onchange="gui.sim.timeStepSize=parseFloat(value)" id="timeStepSizeInput" type="number" value="0.01" step="0.001"></td>
				</tr>
				<tr>
					<td><label for="stiffnessInput">Stiffness</label></td>
					<td><input onchange="gui.sim.stiffness=parseFloat(value)" id="stiffnessInput" type="number" value="1.0" step="0.1"></td>
				</tr>
				<tr>
					<td><label for="gravityInput">Gravity</label></td>
					<td><input onchange="gui.sim.gravity=parseFloat(value)" id="gravityInput" type="number" value="-9.81" step="0.01"></td>
				</tr>
				<tr>
					<td></td>
					<td><button onclick="gui.restart()" type="button" id="restart">Restart</button></td>
				</tr>
				<tr>
					<td></td>
					<td><button onclick="gui.doPause()" type="button" id="Pause">Pause</button></td>
				</tr>
				</table>
			</td>
		</tr>
		<tr><td>
			<h2>Region-Based Shape Matching algorithm:</h2>
				This example shows the region-based shape matching simulation method for deformable solids [BMM17].
				<ol>
					<li>time integration to predict particle positions</li>
					<li class="nostyle"><b>for each region $r$</b></li>
					<li style="margin-left:40px">compute center of mass of current particle configuration</li>
					<li style="margin-left:40px">compute moment matrix</li>
					<li style="margin-left:40px">extract rotation</li>
					<li style="margin-left:40px">determine goal positions</li>
					<li>update positions and velocities</li>
				</ol>

				<p>Shape matching is a meshless method which is geometrically motivated. So it can only generate visually plausible results. The general idea is to match the initial configuration of the particles to the current deformed configuration to get goal positions. Then the particles are pulled towards these goal positions to simulate elasticity. </p>
				
				<p>However, the original shape matching algorithm allows only for
				small deviations from the initial shape. The region-based method performs shape matching in overlapping regions and averages the resulting goal positions for each particle in order to enable the simulation of large deformations. </p>

			<h3>1. Time integration</h3>
				In the first step the particle positions are advected using a symplectic Euler method to obtain predicted positions $\mathbf x^*$:
				$$\begin{align*}
					\mathbf v^* &= \mathbf v(t) + \Delta t \mathbf a_\text{ext}(t) \\
					\mathbf x^* &= \mathbf x(t) + \Delta t \mathbf v^*,
				\end{align*}$$
				where $\mathbf a_\text{ext}$ are accelerations due to external forces.
				Moreover, we initialize the goal positions as $\mathbf g_i = \mathbf 0$.

			<h3>2. Loop over all regions</h3>
				<p>Shape matching is performed for overlapping regions of particles. The resulting goal positions of a particle (one for each region that contains the particle) are averaged. In contrast to performing shape matching for all particles in one single large region, this approach enables the simulation of large deformations.</p>

				<p>Note that in this example we defined one region per particle.</p>


			<h3>3. Compute center of mass</h3>
			
				<p>In order to determine the goal positions for all particles $i$ we search for a rotation matrix $\mathbf R$ and translation vectors $\mathbf t, \mathbf t_0$ which minimize the following equation: 
				$$\begin{equation*}
				\sum_i m_i (\mathbf R(\mathbf X_i - \mathbf t_0) + \mathbf t - \mathbf x^*_i)^2,
			   \end{equation*}
				$$
				where $m$ is the mass, $\mathbf X$ defines the initial configuration and $\mathbf x^*$ the predicted deformed configuration.
				</p>
				<p>
				The optimal translation vectors for region $i$ are the centers of mass of the initial configuration and the deformed configuration:
				$$
				\begin{eqnarray*}
				  \mathbf t_i^0 &=& \frac{1}{\tilde M_i} \sum_{j \in \Re_i} \tilde m_j \mathbf X_j \\
				  \mathbf t_i &=& \frac{1}{\tilde M_i} \sum_{j \in \Re_i} \tilde m_j \mathbf x^*_j,
			    \end{eqnarray*}
			    $$	
				where $\tilde M_i = \sum_{j \in \Re_i} \tilde m_j$.
				Here we use modified particle masses $\tilde m_i = m_i / | \Re_i|$ to consider  the fact that a particle $i$ is part of $| \Re_i|$ regions, where $\Re_i$ is the set of regions containing particle $i$. 
				</p>

				
			<h3>4. Compute moment matrix</h3>
			
				<p>To compute the optimal rotation we move the both configurations so that their centers of mass lie in the origin:
				$$\begin{align*}
					\mathbf q_i &= \mathbf X_i- \mathbf t^0 \\
					\mathbf p_i &= \mathbf x^*_i- \mathbf t
				\end{align*}
				$$
				Moreover, to simplify the problem we search for the optimal linear transformation $\mathbf A$ instead of the optimal rotation $\mathbf R$. This means we have to minimize:
				$$
				\begin{equation*}
				  \sum_i m_i (\mathbf A \mathbf q_i -\mathbf p_i)^2 
				\end{equation*}
				$$
				Setting the derivatives with respect to the coefficients of A to zero gives us the optimal transformation:
				$$\begin{equation*}
					\mathbf A = \left (\sum_{j \in \Re_i} \tilde m_j \mathbf p_j \mathbf q_j^T \right ) \left (\sum_{j \in \Re_i} \tilde m_j \mathbf q_j \mathbf q_j^T \right )^{-1} = \mathbf A_{pq} \mathbf A_{qq}.
			    \end{equation*}
			    $$
				</p>
				
			<h3>5. Extract rotation</h3>
				
				<p>The term $\mathbf A_{qq}$ is symmetric and therefore contains no rotation. Thus, the optimal rotation is the rotational part of $\mathbf A_{pq}$  which can be determined by a polar decomposition:
				$$\mathbf A_{pq} = \mathbf R \mathbf S,$$
				where $\mathbf R$ is an orthogonal matrix and $\mathbf S$ a symmetric matrix. If the configuration is not inverted, $\mathbf R$ is a rotation matrix. Otherwise $\mathbf R$ contains a reflection.					
				</p>
				
			<h3>6. Determine goal positions</h3>
			
				<p>In this step we compute the goal position for each particle of the current region $r$. Since multiple regions can contain the same particle $i$, we sum up all contributions and divide by the number of regions which contain the particle: 
					$$
					\mathbf g_i := \mathbf g_i + \frac{1}{|\Re_i|} (\mathbf R \mathbf q_i + \mathbf t).
					$$
				</p>

			<h3>7. Update positions and velocities</h3>
			
			Finally, we update the positions and velocity by the following modified time integration scheme: 
			$$
			\begin{eqnarray*}
				\mathbf v_i(t+h) &=& \mathbf v_i(t) + \alpha \frac{\mathbf g_i(t) - \mathbf x_i(t)}{h} + \Delta t \frac{\mathbf f_{ext}(t)}{m} \\
				\mathbf x_i(t+h) &=& \mathbf x_i(t) + \Delta t \mathbf v_i(t+h),
			\end{eqnarray*}
			$$
			where $\alpha \in [0,1]$ is a user-defined stiffness parameter. A body is rigid if $\alpha=1$.

			<h3>References</h3>
			<ul>
				<li>[BMM17] Jan Bender, Matthias Müller, Miles Macklin. A Survey on Position Based Dynamics, 2017. Eurographics Tutorials, 2017</li>
			</ul>
			
		</td></tr>
	</table>	
	
</main>

<script id="simulation_code" type="text/javascript">
	class Particle 
	{
		constructor (x, y) 
		{
			// mass 
			this.m = 1.0;
			// reference coordinates
			this.X = x;
			this.Y = y;
			// current coordinates
			this.x = x;
			this.y = y;
			// old positions
			this.xOld = x;
			this.yOld = y;
			// goal position
			this.goal_x = 0.0;
			this.goal_y = 0.0;
			// velocity
			this.vx = 0.0;
			this.vy = 0.0;
			// Number of regions that this particle belongs to
			this.particleRegionCount = 0;
			// mass normalized with respect to region count
			this.normalizedMass = 1.0;
		}
	}

	class Region
	{
		constructor(indexList)
		{
			// mass of this region
			this.m = 0.0;
			// center of mass of this region (initial configuration)
			this.t0_x = 0.0;
			this.t0_y = 0.0;
			// particles in this region
			this.particleIndices = indexList;
		}
	}
	
	class Simulation
	{
		constructor(width, height, fixed, regionHalfWidth)
		{
			// list of particles
			this.particles = [];
			// list of regions
			this.regions = [];
			// dimension of the model
			this.width = width;
			this.height = height;
			// stiffness coefficient 
			this.stiffness = 1;
			this.timeStepSize = 0.01;
			this.time = 0;
			this.gravity = -9.81;
			// region size
			this.regionHalfWidth = regionHalfWidth;
			// indices of boundary particles (only required for rendering)
			this.boundaryIndices = [];			
			this.numFixedParticles = fixed;

			this.init();
		}
		
		init()
		{			
			// create particles and one region per particle
			let i;
			let j;
			let w = this.width;
			let h = this.height;
			let s = 0.1;
			for (i = 0; i < h; i++) 
			{
				for (j = 0; j < w; j++) 
				{
					this.particles.push(new Particle(s*2*j, -s*2*i));

					// init region
					let particleIndices = [];
					for (let k = Math.max(0, i - this.regionHalfWidth); k < Math.min(h, i + this.regionHalfWidth + 1); ++k)
					{
						for (let l = Math.max(0, j - this.regionHalfWidth); l < Math.min(w, j + this.regionHalfWidth + 1); ++l)
						{
							let index = k * w + l;
							particleIndices.push(index);
						}
					}
					this.regions.push(new Region(particleIndices));
				}
			}

			// set mass of inactive particles to high value so that they
			// almost do not move in the shape matching correction
			for (let i = 0; i < this.particles.length; ++i)
			{
				if (!this.isActive(i))
					this.particles[i].m = 100000.0;
			}
			
			// normalize particle masses and compute center of mass of each region
			this.computeParticleRegionCount();
			this.normalizeParticleMasses();
			this.computeCentersOfMass();

			
			// Determine which indices form the boundary. They will be rendered in order,
			// so we have to carefully pick the indices in the right order too (clockwise in this case).

			// Top boundary
			for (let j = 0; j < w; ++j) 
			{
				this.boundaryIndices.push(j);
			}

			// Right boundary
			for (let i = 0; i < h; ++i) 
			{
				const j = w - 1;
				const index = i * w + j;
				this.boundaryIndices.push(index);
			}

			// Bottom boundary
			for (let j = w - 1; j >= 0; --j) {
				const i = h - 1;
				const index = i * w + j;
				this.boundaryIndices.push(index);
			}

			// Left boundary
			for (let i = h - 1; i >= 0; --i)
			{
				const j = 0;
				const index = i * w + j;
				this.boundaryIndices.push(index);
			}

		}

		// For each particle determine in how many regions it is included.
		// This value is required to compute the normalized region masses.
		computeParticleRegionCount()
		{
			for (let i = 0; i < this.particles.length; ++i)
				this.particles[i].particleRegionCount = 0;

			for (const region of this.regions)
			{
				for (const particleIndex of region.particleIndices)
				{
					this.particles[particleIndex].particleRegionCount += 1;
				}
			}
		}

		// Divide the particle mass by the number of regions in which it is included 
		// to obtain the normalized mass.
		normalizeParticleMasses()
		{
			for (let i = 0; i < this.particles.length; ++i)
			{
				let p = this.particles[i];
				p.normalizedMass = p.m / p.particleRegionCount;
			}
		}

		// compute center of mass of reference configuration for each region
		computeCentersOfMass()
		{
			for (let i = 0; i < this.regions.length; i++) 
			{
				let r = this.regions[i];
				r.t0_x = 0.0;
				r.t0_y = 0.0;
				r.m = 0.0;
				for (let j = 0; j < r.particleIndices.length; j++) 
				{
					let p = this.particles[r.particleIndices[j]];
					r.t0_x += p.normalizedMass * p.X;
					r.t0_y += p.normalizedMass * p.Y;
					r.m += p.normalizedMass;
				}
				r.t0_x /= r.m;
				r.t0_y /= r.m;	
			}
		}

		isActive(i)
		{
			// particle 0 and the selected particle are fixed
			if (this.numFixedParticles == 1)
				return (i != 0) && (i != gui.selectedParticle);
			else if (this.numFixedParticles == 2)
				return (i != 0) && (i != this.width-1) && (i != gui.selectedParticle);
			else 
				return (i != 0) && (i != this.width-1) && (i != (this.height-1)*this.width) && (i != this.height*this.width-1) && (i != gui.selectedParticle);
		}
		
		// perform a 2D polar decomposition: M = R*S
		// where R is an orthogonal matrix and S is a symmetric matrix
		polarDecomposition(M)
		{
			// determinant
			let det = M[0][0]*M[1][1] - M[0][1]*M[1][0];
			
			let R = [[0,0],[0,0]];
			if (det >= 0.00001)
			{
				R[0][0] = M[0][0] + M[1][1];
				R[0][1] = M[0][1] - M[1][0];
				R[1][0] = M[1][0] - M[0][1];
				R[1][1] = M[1][1] + M[0][0];
			}
			else
			{
				R[0][0] = M[0][0] - M[1][1];
				R[0][1] = M[0][1] + M[1][0];
				R[1][0] = M[1][0] + M[0][1];
				R[1][1] = M[1][1] - M[0][0];
			}

			let dl1 = R[1][0]*R[1][0] + R[0][0]*R[0][0];
			let dl2 = R[1][1]*R[1][1] + R[0][1]*R[0][1];

			if ((dl1 < 1.0e-12) || (dl2 < 1.0e-12))
			{
				R = [[1,0], [0,1]];
				return R;
			}

			let l1 = 1.0/Math.sqrt(dl1);
			let l2 = 1.0/Math.sqrt(dl2);

			R[0][0] *= l1;
			R[1][0] *= l1;
			R[0][1] *= l2;
			R[1][1] *= l2;
			return R;
		}

		// Compute the sum A+B
		matrixAdd(A, B)
		{
			return [[A[0][0]+B[0][0], A[0][1]+B[0][1]], 
					[A[1][0]+B[1][0], A[1][1]+B[1][1]]];
		}
		
		// Compute the product A*x
		matrixVecProd(A, x)
		{
			return [A[0][0]*x[0] + A[0][1]*x[1], 
					A[1][0]*x[0] + A[1][1]*x[1]];
		}
				
		// simulation step
		simulationStep()
		{	
			let dt = this.timeStepSize;
						
			// compute preview using symplectic Euler
			for (let i = 0; i < this.particles.length; i++) 
			{
				let p = this.particles[i];
				
				// store old position for velocity update after the position correction
				p.xOld = p.x;	
				p.yOld = p.y;

				// reset goal position
				p.goal_x = 0.0;
				p.goal_y = 0.0;
						
				// integrate position
				if (this.isActive(i))
				{
					p.x = p.x + dt * p.vx;
					p.y = p.y + dt * p.vy + dt*dt*this.gravity;
				}
			}
	
			// region-based shape matching
			for (let ri = 0; ri < this.regions.length; ++ri)
			{
				// compute center of mass of current configuration in region i
				let r = this.regions[ri];
				let t_x = 0.0;
				let t_y = 0.0;
				for (let i = 0; i < r.particleIndices.length; i++) 
				{
					let p = this.particles[r.particleIndices[i]];
					t_x += p.normalizedMass * p.x;
					t_y += p.normalizedMass * p.y;
				}
				t_x /= r.m;
				t_y /= r.m;	

				// compute moment matrix
				let Apq = [[0,0], [0,0]];
				for (let i = 0; i < r.particleIndices.length; i++) 
				{
					let part = this.particles[r.particleIndices[i]];
					let m = part.normalizedMass;
					let q = [part.X - r.t0_x, part.Y - r.t0_y];
					let p = [part.x - t_x, part.y - t_y];
					let dyadic = [[m*p[0]*q[0], m*p[0]*q[1]], [m*p[1]*q[0], m*p[1]*q[1]]];
					Apq = this.matrixAdd(Apq, dyadic);
				}

				// extract rotation
				const R = this.polarDecomposition(Apq);

				// update goal positions
				for (let i = 0; i < r.particleIndices.length; i++) 
				{
					let part = this.particles[r.particleIndices[i]];
					let q = [part.X - r.t0_x, part.Y - r.t0_y];

					// Add contribution from this region to goal position of particle j
					let Rq = this.matrixVecProd(R, q);
					part.goal_x += Rq[0] + t_x;
					part.goal_y += Rq[1] + t_y;			
				}
			}		

			// update positions & velocities
			for (let i = 0; i < this.particles.length; i++) 
			{
				let p = this.particles[i];

				// normalize goal positions 
				p.goal_x /= p.particleRegionCount;
				p.goal_y /= p.particleRegionCount;

				if (this.isActive(i))
				{
					p.vx += this.stiffness/dt * (p.goal_x - p.x);
					p.vy += this.stiffness/dt * (p.goal_y - p.y) + dt*this.gravity;
					p.x = p.xOld + dt * p.vx;	
					p.y = p.yOld + dt * p.vy;		
				}
			}

			// update simulation time
			this.time = this.time + dt;
		}
	}
	
	class GUI
	{
		constructor()
		{
			this.canvas = document.getElementById("simCanvas");
			this.c = this.canvas.getContext("2d");
			this.requestID = -1;
			
			this.timeSum = 0.0;
			this.counter = 0;
			this.pause = false;
			this.origin = { x : this.canvas.width / 2, y : this.canvas.height/2};
			this.zoom = 50;
			this.particleRadius = 0.025;
			this.selectedParticle = -1;
			
			// register mouse event listeners (zoom/selection)
			this.canvas.addEventListener("mousedown", this.mouseDown.bind(this), false);
			this.canvas.addEventListener("mousemove", this.mouseMove.bind(this), false);
			this.canvas.addEventListener("mouseup", this.mouseUp.bind(this), false);	
			this.canvas.addEventListener("wheel", this.wheel.bind(this), false);		
		}
		
		// set simulation parameters from GUI and start mainLoop
		restart()
		{
			window.cancelAnimationFrame(this.requestID);
			let w = parseInt(document.getElementById('widthInput').value);
			let h = parseInt(document.getElementById('heightInput').value);
			let f = parseInt(document.getElementById('fixedParticlesInput').value);
			let hw = parseInt(document.getElementById('regionHalfWidthInput').value);
			delete this.sim;
			this.sim = new Simulation(w, h, f, hw);
			
			this.timeSum = 0.0;
			this.counter = 0;

			this.sim.stiffness = parseFloat(document.getElementById('stiffnessInput').value);
			this.sim.gravity = parseFloat(document.getElementById('gravityInput').value);
			this.sim.timeStepSize = parseFloat(document.getElementById('timeStepSizeInput').value);

			document.getElementById("numParticles").innerHTML = this.sim.particles.length;
			this.mainLoop();
		}
		
		drawCoordinateSystem()
		{
			// draw x-axis
			this.c.strokeStyle = "#FF0000";
			this.c.beginPath();
			this.c.moveTo(this.origin.x, this.origin.y);
			this.c.lineTo(this.origin.x+1*this.zoom, this.origin.y);
			this.c.stroke();	

			// draw y-axis
			this.c.strokeStyle = "#00FF00";	
			this.c.beginPath();
			this.c.moveTo(this.origin.x, this.origin.y);
			this.c.lineTo(this.origin.x, this.origin.y-1*this.zoom);
			this.c.stroke();		
		}
		
		numToHex(n) 
		{ 
			let hex = Number(n).toString(16);
			if (hex.length < 2)
			   hex = "0" + hex;
			return hex;
		}
		
		hsvToRgb(h, s, v)
		{
			let i = Math.floor(h * 6);
			let f = h * 6 - i;
			let p = v * (1 - s);
			let q = v * (1 - f * s);
			let t = v * (1 - (1 - f) * s);

			let 	r = 0;
			let 	g = 0;
			let 	b = 0;
			switch (i % 6)
			{
				case 0: r = v, g = t, b = p; break;
				case 1: r = q, g = v, b = p; break;
				case 2: r = p, g = v, b = t; break;
				case 3: r = p, g = q, b = v; break;
				case 4: r = t, g = p, b = v; break;
				case 5: r = v, g = p, b = q; break;
			}
			let rh = this.numToHex(Math.round(r * 255));
			let gh = this.numToHex(Math.round(g * 255));
			let bh = this.numToHex(Math.round(b * 255));
			//console.log(rh)
			return "#" + rh + gh + bh;
		}
		
		draw()
		{
			this.c.clearRect(0, 0, this.canvas.width, this.canvas.height);
			
			this.drawCoordinateSystem();

			// draw boundary
			if (this.sim.boundaryIndices.length > 0)
			{
				this.c.strokeStyle = "#666666";
				this.c.fillStyle = "#556677";
				this.c.beginPath();
				let p = this.sim.particles[this.sim.boundaryIndices[0]];
				this.c.moveTo(this.origin.x + p.x*this.zoom, this.origin.y - p.y*this.zoom);
				for (let i = 1; i < this.sim.boundaryIndices.length; i++) 
				{
					let p = this.sim.particles[this.sim.boundaryIndices[i]];
					this.c.lineTo(this.origin.x + p.x*this.zoom, this.origin.y - p.y*this.zoom);
				}
				this.c.fill();
			}
					
			// draw particles as circles
			for (let i = 0; i < this.sim.particles.length; i++) 
			{
				let p = this.sim.particles[i];
				let r = this.particleRadius;
				if (i == this.selectedParticle)
				{
					// draw selected particle in red with larger radius
					this.c.fillStyle = "#FF0000";	
					r = 3*r;
				}
				else 
					this.c.fillStyle = "#0000FF";	
				let px = this.origin.x + p.x * this.zoom;
				let py = this.origin.y - p.y * this.zoom;
				
				this.c.beginPath();			
				this.c.arc(px, py, r * this.zoom, 0, Math.PI*2, true); 
				this.c.closePath();
				this.c.fill();
			}
		}
		
		mainLoop()
		{
			// perform multiple sim steps per render step
			for (let i=0; i < 8; i++)
			{
				let t0 = performance.now();
				
				if (!this.pause)
				this.sim.simulationStep();

				let t1 = performance.now();
				this.timeSum += t1 - t0;
				this.counter += 1;
					
				if (this.counter % 50 == 0)				
				{
					this.timeSum /= this.counter;
					document.getElementById("timePerStep").innerHTML = this.timeSum.toFixed(2);		
					this.timeSum = 0.0;
					this.counter = 0;
				}	
				document.getElementById("time").innerHTML = this.sim.time.toFixed(2);	
			}			

			this.draw();
			if (!this.pause)
				this.requestID = window.requestAnimationFrame(this.mainLoop.bind(this));		
		}
		
		doPause()
		{
			this.pause = !this.pause;
			if (!this.pause)
				this.mainLoop();
		}
		
		mouseDown(event)
		{
			// left mouse button down
			if (event.which == 1)
			{
				let mousePos = this.getMousePos(this.canvas, event);
				for (let i = 0; i < this.sim.particles.length; i++) 
				{
					let p = this.sim.particles[i];
					let px = this.origin.x + p.x * this.zoom;
					let py = this.origin.y - p.y * this.zoom;
					let dx = px - mousePos.x
					let dy = py - mousePos.y
					let dist2 = Math.sqrt(dx * dx + dy * dy)
					if (dist2 < 10)
					{
						this.selectedParticle = i;
						break;
					}				
				}			
			}
		}
		
		getMousePos(canvas, event) 
		{
			let rect = canvas.getBoundingClientRect();
			return {
				x: event.clientX - rect.left,
				y: event.clientY - rect.top
			};
		}
	
		mouseMove(event)
		{
			if (this.selectedParticle != -1)
			{
				let mousePos = this.getMousePos(this.canvas, event);		
				this.sim.particles[this.selectedParticle].x = (mousePos.x - this.origin.x) / this.zoom;
				this.sim.particles[this.selectedParticle].y = -(mousePos.y - this.origin.y) / this.zoom;
			}
		}
		
		mouseUp(event) 
		{
			this.selectedParticle = -1;
		}
		
		wheel(event)
		{
			event.preventDefault();
			this.zoom += event.deltaY * -0.05;
			if (this.zoom < 1) 
				this.zoom = 1;
		}
	}
	
	gui = new GUI();
	gui.restart();
</script>

</body>
</html>