---
layout: post
title: "Curves and Beer: The Hough Transform"
categories: [Hacking]
tags: [Featured, Software, _Phase4]
---

<figure class="align-center">
	<canvas width="392" height="392" id="km_sample5" class="align-center">
	</canvas>
  <figcaption id='htStats'>Lines from Pixels</figcaption>
</figure>
<button id="reloader" onclick="htSample1.load_new_image();" class="btn--primary align-center">Press for a Different Random Photo, or to Restart</button>

<script>
	// Hough Transform:
	// Sample Points that bounce around in their little windows
	function bouncePt(x,y,vx,vy,w,h,x0,y0,decay) {
		this.x = x ? x : 0;
		this.y = y ? y : 0;
		this.vx = vx ? vx : 0;
		this.vy = vy ? vy : 0;
		this.w = w ? w : 128;
		this.h = h ? h : 128;
		this.x0 = x0 ? x0 : 0;
		this.y0 = y0 ? y0 : 0;
		this.decay = (decay===undefined) ? 0.999 : decay;
	}
	bouncePt.prototype.advance = function()
	{
		this.x += this.vx;
		this.y += this.vy;
		if ((this.x<0)||(this.x>=this.w)) {
			this.vx = -this.vx + (Math.random()*3-1.5);
			this.vy += (Math.random()-0.5);
		}
		if ((this.y<0)||(this.y>=this.h)) {
			this.vx += Math.random()-0.5;
			this.vy = -this.vy + (Math.random()*3-1.5);
		}
		this.vx *= this.decay;
		this.vy *= this.decay;
		this.vx += 0.05*(Math.random()-0.5); // just a little wandering...
		this.vy += 0.05*(Math.random()-0.5); // just a little wandering...
		this.x = this.x<0 ? -this.x : this.x;
		this.x = this.x>=this.w ? 2*this.w-this.x-1 : this.x;
		this.y = this.y<0 ? -this.y : this.y;
		this.y = this.y>=this.h ? 2*this.h-this.y-1 : this.y;
	}
	bouncePt.prototype.crosshair = function(ctx) {
		ctx.strokeStyle = '#b0ffb0ff';
		var nx = this.x0+this.x;
		var ny = this.y0+this.y;
		ctx.beginPath();
		ctx.moveTo(nx,this.y0);
		ctx.lineTo(nx,this.y0+this.h-1);
		ctx.moveTo(this.x0,ny);
		ctx.lineTo(this.x0+this.w-1,ny);
		ctx.stroke();
		ctx.strokeStyle = '#a0a0a0ff';
		ctx.beginPath();
		ctx.rect(this.x0,this.y0,this.w,this.h);
		var d = 2;
		ctx.moveTo(nx-d,ny-d);
		ctx.lineTo(nx+d,ny+d);
		ctx.moveTo(nx-d,ny+d);
		ctx.lineTo(nx+d,ny-d);
		ctx.stroke();
		ctx.fillStyle = 'black';
		ctx.fillRect(nx,ny,1,1);  	
	}
	//
	var htSample1 = {
		canvas: null,
		ctx: null,
		canvasB: null,
		ctxB: null,
		canvasC: null,
		ctxC: null,
		canvasD: null,
		ctxD: null,
		canvasE: null,
		ctxE: null,
		hiddenCanvas: null,
		iter: 0,
		pauseCount: 0,
		srcImage: null,
		paused: false,
		// various imageData objects
		origPxls: null,
		grayPxls: null,
		edgePxls: null,
		mixdPxls: null,
		linePxls: null,
		mixdPxls: null,
		overPxls: null,
		aSize: {x:256, y:128},
		aPxls: null,
		cosa: null,
		sina: null,
		maxP: null,
		accumulator: null,
		samplePtB: null,
		samplePtC: null,
		stats: null,
		debug: false,
		placeHolderURLs: ['https://loremflickr.com/g/128/128/vespa,scooter','https://picsum.photos/g/128/128/?random',"https://www.fillmurray.com/128/128"],
		//
		filters: { // based on https://www.html5rocks.com/en/tutorials/canvas/imagefilters/
			filterImage: function(filter, image, var_args) {
			  var args = [this.getPixels(image)];
			  for (var i=2; i<arguments.length; i++) {
			    args.push(arguments[i]);
			  }
			  return filter.apply(null, args);
			},
			grayscale: function(pixels, args) {
			  var d = pixels.data;
			  for (var i=0; i<d.length; i+=4) {
			    var v = (d[i]+d[i+1]+d[i+2])/3;
			    d[i] = d[i+1] = d[i+2] = v;
			  }
			  return pixels;
			},
			grayscale_lum: function(pixels, args) {
			  var d = pixels.data;
			  for (var i=0; i<d.length; i+=4) {
			    var r = d[i];
			    var g = d[i+1];
			    var b = d[i+2];
			    // CIE luminance for the RGB
			    // The human eye is bad at seeing red and blue, so we de-emphasize them.
			    var v = 0.2126*r + 0.7152*g + 0.0722*b;
			    d[i] = d[i+1] = d[i+2] = v;
			  }
			  return pixels;
			},
			convolute: function(pixels, weights, opaque) {
			  var side = Math.round(Math.sqrt(weights.length));
			  var halfSide = Math.floor(side/2);
			  var src = pixels.data;
			  var sw = pixels.width;
			  var sh = pixels.height;
			  // pad output by the convolution matrix
			  var w = sw;
			  var h = sh;
			  var output = Filters.createImageData(w, h);
			  var dst = output.data;
			  // go through the destination image pixels
			  var alphaFac = opaque ? 1 : 0;
			  for (var y=0; y<h; y++) {
			    for (var x=0; x<w; x++) {
			      var sy = y;
			      var sx = x;
			      var dstOff = (y*w+x)*4;
			      // calculate the weighed sum of the source image pixels that
			      // fall under the convolution matrix
			      var r=0, g=0, b=0, a=0;
			      for (var cy=0; cy<side; cy++) {
				for (var cx=0; cx<side; cx++) {
				  var scy = sy + cy - halfSide;
				  var scx = sx + cx - halfSide;
				  if (scy >= 0 && scy < sh && scx >= 0 && scx < sw) {
				    var srcOff = (scy*sw+scx)*4;
				    var wt = weights[cy*side+cx];
				    r += src[srcOff] * wt;
				    g += src[srcOff+1] * wt;
				    b += src[srcOff+2] * wt;
				    a += src[srcOff+3] * wt;
				  }
				}
			      }
			      dst[dstOff] = r;
			      dst[dstOff+1] = g;
			      dst[dstOff+2] = b;
			      dst[dstOff+3] = a + alphaFac*(255-a);
			    }
			  }
			  return output;
			},
			convoluteFloat32: function(pixels, weights, opaque) {
	          var side = Math.round(Math.sqrt(weights.length));
	          var halfSide = Math.floor(side/2);

	          var src = pixels.data;
	          var sw = pixels.width;
	          var sh = pixels.height;

	          var w = sw;
	          var h = sh;
	          var output = {
	            width: w, height: h, data: new Float32Array(w*h*4)
	          };
	          var dst = output.data;

	          var alphaFac = opaque ? 1 : 0;

	          for (var y=0; y<h; y++) {
	            for (var x=0; x<w; x++) {
	              var sy = y;
	              var sx = x;
	              var dstOff = (y*w+x)*4;
	              var r=0, g=0, b=0, a=0;
	              for (var cy=0; cy<side; cy++) {
	                for (var cx=0; cx<side; cx++) {
	                  var scy = Math.min(sh-1, Math.max(0, sy + cy - halfSide));
	                  var scx = Math.min(sw-1, Math.max(0, sx + cx - halfSide));
	                  var srcOff = (scy*sw+scx)*4;
	                  var wt = weights[cy*side+cx];
	                  r += src[srcOff] * wt;
	                  g += src[srcOff+1] * wt;
	                  b += src[srcOff+2] * wt;
	                  a += src[srcOff+3] * wt;
	                }
	              }
	              dst[dstOff] = r;
	              dst[dstOff+1] = g;
	              dst[dstOff+2] = b;
	              dst[dstOff+3] = a + alphaFac*(255-a);
	            }
	          }
	          return output;
	        },
			sobel: function(px,edgePx) {
				// px = Filters.grayscale(px); // let's assume it's already gray
				// original method created edgePx, but we can avoid reallocations
				// Note that ImageData values are clamped between 0 and 255, so we need
				// to use a Float32Array for the gradient values because they
				// range between -255 and 255.
				var vertical = this.convoluteFloat32(px,
				  [ -1, 0, 1,
				    -2, 0, 2,
				    -1, 0, 1 ]);
				var horizontal = this.convoluteFloat32(px,
				  [ -1, -2, -1,
				     0,  0,  0,
				     1,  2,  1 ]);
				//var edgePx = Filters.createImageData(vertical.width, vertical.height);
				for (var i=0; i<edgePx.data.length; i+=4) {
				  // make the vertical gradient red
				  var v = Math.abs(vertical.data[i]);
				  edgePx.data[i] = v;
				  // make the horizontal gradient green
				  var h = Math.abs(horizontal.data[i]);
				  edgePx.data[i+1] = h;
				  // and mix in some blue for aesthetics
				  edgePx.data[i+2] = (v+h)/4;
				  edgePx.data[i+3] = 255; // opaque alpha
				}
				return edgePx;
			},
		},
		accumulate_hough_points: function(y) {
		  // look for potential edge points in row "y" of edgePxls, add to our accumulator
			var t = 0;
			for (var x=0; x<this.edgePxls.width; x+=1) {
				var i = x + y*this.edgePxls.width;
				var j = i*4;
				var hh = this.aSize.y / 2;
				var v = this.edgePxls.data[j] + this.edgePxls.data[j+1];
				if (v > 190) {	// TODO: threshhold value
					for(k=0; k<this.aSize.x; k+=1) { // iterate through possible angles
						var p = x*this.cosa[k] + y*this.sina[k];
						var pn = Math.floor(hh + hh * p / this.maxP);
						if (p<0) {
							t += 1;
						}
						var accIndex = pn*this.aSize.x + k;
						this.accumulator[accIndex] += 1;
					}
				}
			}
		},
		update_hough_pixels: function() {
		  // interpret our floating-point accumulator as an 8-bit image
			this.maxA = Math.max.apply(Math, this.accumulator);
			if (this.maxA == 0.0) {
				// TODO
				return;
			}
			for (var i=0; i<this.accumulator.length; i+=1) {
				var j = i*4;
				var v = Math.floor(255*this.accumulator[i]/this.maxA);
				this.aPxls.data[j] = this.aPxls.data[j+1] = this.aPxls.data[j+2] = v;
			}
		},
		update_hough_line_overlay: function(ctx, pixels, x0, y0) {
			ctx.save();
			// this version draws red lines and a "scanner" line over the grayscale image
			// ctx.fillStyle = 'white';
			// ctx.fillRect(x0,y0,this.edgePxls.width,this.edgePxls.height);
			var w = pixels.width;
			var h = pixels.height;
			ctx.putImageData(pixels, x0, y0);
			ctx.beginPath();
			ctx.rect(x0,y0,w,h);
			ctx.clip();
			if (this.iter< h) {
				ctx.strokeStyle = 'blue';
				ctx.moveTo(x0,   y0+this.iter);
				ctx.lineTo(x0+w, y0+this.iter);
				ctx.stroke();
			}
			// this.maxA = Math.max.apply(Math, this.accumulator); // also calc'd in update_hough_pixels()
			if (this.maxA < 1) {
				ctx.restore();
				return 0;
			}
			ctx.strokeStyle = '#20ff0040';
			var i = 0;
			var n = 0;
			var mx = w-1;
			var hh = this.aSize.y / 2;
			for (var y=0; y<this.aSize.y; y+=1) {
				for (var x=0; x<this.aSize.x; x+=1) {
					if (this.accumulator[i] >= (this.maxA*0.8)) { // TODO: tune this
						var rho = this.maxP*((y-hh)/hh);
						ctx.beginPath();
						if (this.sina[x] != 0.0) {
							var yz = (rho-0*this.cosa[x])/this.sina[x];
							var ym = (rho-mx*this.cosa[x])/this.sina[x];
							ctx.moveTo(x0+0, y0+yz);
							ctx.lineTo(x0+mx,y0+ym);
						} else {
							ctx.moveTo(x0+rho,y0);
							ctx.lineTo(x0+rho,y0+h);
						}
						ctx.stroke();
						n += 1;
					}
					i += 1;
				}
			}
			ctx.restore();
			// console.log(n+' strokes');
			return n;
		},
		update_hough_lines_only: function(x0, y0) {
			// this version draws white lines on black, no "scanner"
			var w = this.grayPxls.width;
			var h = this.grayPxls.height;
			this.ctx.save();
			this.ctx.fillStyle = 'black';
			this.ctx.fillRect(x0,y0,w,h);
			// this.maxA = Math.max.apply(Math, this.accumulator); // also calc'd in update_hough_pixels()
			if (this.maxA < 1) {
				this.ctx.restore();
				return 0;
			}
			this.ctx.beginPath();
			this.ctx.rect(x0,y0,w,h);
			this.ctx.clip();
			this.ctx.strokeStyle = '#ffffffff';
			var i = 0;
			var n = 0;
			var mx = w-1;
			var hh = this.aSize.y / 2;
			for (var y=0; y<this.aSize.y; y+=1) {
				for (var x=0; x<this.aSize.x; x+=1) {
					if (this.accumulator[i] >= (this.maxA*0.8)) { // TODO: tune this
						var p = this.maxP*((y-hh)/hh);
						this.ctx.beginPath();
						if (this.sina[x] != 0.0) {
							var yz = (p-0*this.cosa[x])/this.sina[x];
							var ym = (p-mx*this.cosa[x])/this.sina[x];
							this.ctx.moveTo(x0+0,y0+yz);
							this.ctx.lineTo(x0+mx,y0+ym);
						} else {
							this.ctx.moveTo(x0+p,y0);
							this.ctx.lineTo(x0+p,y0+h);
						}
						this.ctx.stroke();
						n += 1;
					}
					i += 1;
				}
			}
			this.ctx.restore();
			// console.log(n+' strokes');
			return n;
		},
		edge_mix: function() {
			for (var i=0; i<this.mixdPxls.data.length; i+=4) {
				var l = this.linePxls.data[i];
				var e = Math.min(255,this.edgePxls.data[i]+this.edgePxls.data[i+1]);
				var m = l * e / 255;
				this.mixdPxls.data[i] =
				this.mixdPxls.data[i+1] =
				this.mixdPxls.data[i+2] = m;
				this.mixdPxls.data[i+3] = 255;
				//var g = this.grayPxls.data[i];
				var f = m > 128;
				this.overPxls.data[i]   = f ? m : this.origPxls.data[i]  ; // g;
				this.overPxls.data[i+1] = f ? 0 : this.origPxls.data[i+1]; // g;
				this.overPxls.data[i+2] = f ? 0 : this.origPxls.data[i+2]; // g;
				this.overPxls.data[i+3] = 255;
			}
		},
		update_canvas: function() {
			// draw points AND means
			this.ctx.fillStyle = 'white';
			this.ctx.fillRect(0,0,this.canvas.width,this.canvas.height);
			if (this.origPxls == null) { // should be a very rare case
				this.ctx.drawImage(this.srcImage, 0,0);
				this.ctx.strokeStyle = 'red';
				this.ctx.beginPath();
				this.ctx.moveTo(128,0);
				this.ctx.lineTo(255,127);
				this.ctx.lineTo(383,0);
				this.ctx.moveTo(128,127);
				this.ctx.lineTo(255,0);
				this.ctx.lineTo(383,127);
				this.ctx.stroke();
				return;
			}
			this.ctx.putImageData(this.origPxls, 0, 0);
			this.ctx.putImageData(this.grayPxls, 132, 0);
			this.ctx.putImageData(this.edgePxls, 262, 0);
			this.ctx.putImageData(this.aPxls, 0, 132);
			var nLines = this.update_hough_line_overlay(this.ctx,this.grayPxls, 262,132);
			var nLines2 = this.update_hough_lines_only(0,262);
			this.linePxls = this.ctx.getImageData(0,262,this.edgePxls.width,this.edgePxls.height);
			this.edge_mix();
			this.ctx.putImageData(this.mixdPxls, 132, 262);
			this.ctx.putImageData(this.overPxls, 262, 262);
			this.update_stats(nLines);
		},
		update_canvasE: function() {
			if (! this.canvasE) return;
			this.ctxE.putImageData(this.aPxls, 0, 0);
			var nLines = this.update_hough_line_overlay(this.ctxE,this.origPxls, 262,0);
		},
		update_canvasD: function() {
			if (! this.canvasD) return;
			this.ctxD.putImageData(this.origPxls, 0, 0);
			this.ctxD.putImageData(this.grayPxls, 132, 0);
			this.ctxD.putImageData(this.edgePxls, 262, 0);
		},
		update_canvasB: function() {
			// canvasB shows a sample point in x/y and draws a matching curve of 
			//    possible theta/rho values across the range of theta of +/- 90 degrees
			if (! this.canvasB) return;
			var w = this.edgePxls.width;
			var h = this.edgePxls.height;
			var aw = this.aSize.x;
			var ah = this.aSize.y;
			var hh = ah / 2;
			this.ctxB.fillStyle = 'white';
			this.ctxB.fillRect(0,0,w,h);
			this.ctxB.fillStyle = 'black';
			this.ctxB.fillRect(w+4,0,aw,ah);
			for (var i=0; i<this.samplePtB.length; i+=1) {
				this.samplePtB[i].advance(w,h);
				this.samplePtB[i].crosshair(this.ctxB,0,0,w,h);
				this.ctxB.fillStyle = '#ffffff80';
				var x = this.samplePtB[i].x;
				var y = this.samplePtB[i].y;
				for(var k=0; k<aw; k+=1) { // iterate through possible angles
					var p = x*this.cosa[k] + y*this.sina[k];
					var pn = Math.floor(hh + hh * p / this.maxP);
					this.ctxB.fillRect(w+4+k,pn,1,1);
				}
			}
			this.ctxB.fillStyle = 'white';
			this.ctxB.font = "Source Sans Pro 16px";
			this.ctxB.fillText("ρ",w+10,ah/2+6);
			this.ctxB.fillText("θ",w+4+aw/2,ah-4);
		},
		update_canvasC: function() {
			// canvasC shows a simple point in theta/rho and draws a line in x/y
			if (! this.canvasC) return;
			var aw = this.aSize.x;
			var ah = this.aSize.y;
			var w = this.edgePxls.width;
			var h = this.edgePxls.height;
			this.samplePtC.advance(aw,ah);
			var x = this.samplePtC.x;
			var y = this.samplePtC.y;
			this.ctxC.fillStyle = 'white';
			this.ctxC.fillRect(0,0,aw,ah);
			this.samplePtC.crosshair(this.ctxC,0,0,aw,ah);
			//
			var x0 = aw+4;
			var mx = w-1;
			var cx = w/2;
			var cy = h/2;
			var hh = this.aSize.y / 2;
			var k = Math.floor(this.samplePtC.x);
			var ck = this.cosa[k];
			var sk = this.sina[k];
			var rho = 0.35*this.maxP*((y-hh)/hh); // cheat for demo: we scale-down rho
			this.ctxC.fillStyle = 'black';
			this.ctxC.fillRect(x0,0,w,h);
			this.ctxC.strokeStyle = 'white';
			this.ctxC.save();
			this.ctxC.beginPath();
			this.ctxC.rect(x0,0,w,h);
			this.ctxC.clip();
			this.ctxC.beginPath();
			if (sk==0) { // vertical
				var x1 = x0 + rho + cx; // cheat for demo: we center the origin
				this.ctxC.moveTo(x1,0);
				this.ctxC.lineTo(x1,h);
			} else {
				var yz = (rho-(-cx)*ck)/sk;
				var ym = (rho-(mx-cx)*ck)/sk;
				this.ctxC.moveTo(x0+0, cy+yz);
				this.ctxC.lineTo(x0+mx,cy+ym);
			}
			this.ctxC.stroke();
			this.ctxC.restore();
			this.ctxC.font = "Source Sans Pro 16px";
			this.ctxC.fillText("ρ",6,ah/2+6);
			this.ctxC.fillText("θ",aw/2,ah-4);
		},
		update_stats: function(counter) {
			if (!this.stats) {
				this.stats = document.getElementById('htStats');
			}
			text = 'Found '+counter+' strongest line';
			if (counter>1) {
				text = text+'s';
			}
			// if (this.debug) {
			// 	text = text + ', maxA: '+Math.floor(this.maxA);
			// }
			this.stats.innerHTML = text;
		},
		update_all: function() {
			if (this.iter >= this.origPxls.height) {
				return 0;
			}
			this.accumulate_hough_points(this.iter);
			this.update_hough_pixels();
			return 1;
		},
		looper: function(timestamp) {
			// called by requestAnimationFrame() forever
			if (this.paused) {
				window.requestAnimationFrame(this.looper.bind(this));
				return;
			}
			var ch = this.update_all();
			if (this.pauseCount < 1) {
				this.update_canvas();
			}
			this.update_canvasB();
			this.update_canvasC();
			this.update_canvasD();
			this.update_canvasE();
			if (ch > 0) {
				this.iter += 1;
			} else {
				this.pauseCount += 1;
			}
			if (this.pauseCount > 200) {
				window.setTimeout(this.load_new_image.bind(this),30);
			} else {
				window.requestAnimationFrame(this.looper.bind(this));
			}
		},
		startup: function() {
			// also called whenever we re-start animation
			this.iter = this.pauseCount = 0;
			window.requestAnimationFrame(this.looper.bind(this));
		},
		collect_gray_pixels: function(context,x0,y0,w,h) {
			var i;
			this.origPxls = context.getImageData(x0,y0,w,h);
			if (this.grayPxls === null) {
				this.grayPxls = context.createImageData(w, h);
			}
			if (this.edgePxls === null) {
				this.edgePxls = context.createImageData(w, h);
			}
			if (this.mixdPxls === null) {
				this.mixdPxls = context.createImageData(w, h);
			}
			if (this.overPxls === null) {
				this.overPxls = context.createImageData(w, h);
			}
			for (i = 0; i < this.origPxls.data.length; i+=1) {
				this.grayPxls.data[i] = this.origPxls.data[i];
			}
			//this.grayPxls = this.filters.grayscale(this.grayPxls);
			this.filters.grayscale_lum(this.grayPxls);
			this.filters.sobel(this.grayPxls, this.edgePxls);
			// prep accumulation for Hough transform
			this.maxP = Math.sqrt(this.edgePxls.width * this.edgePxls.width +
									this.edgePxls.height * this.edgePxls.height);
			this.maxA = 0.0;
			if (this.accumulator === null) {
				this.accumulator = new Float32Array(this.aSize.x * this.aSize.y);
				this.aPxls = context.createImageData(this.aSize.x, this.aSize.y);
				this.cosa = new Float32Array(this.aSize.x);
				this.sina = new Float32Array(this.aSize.x);
				for (i = 0; i < this.aSize.x; i+=1) {
					var a = Math.PI*(i/(this.aSize.x-1)) - (Math.PI/2);
					this.cosa[i] = Math.cos(a);
					this.sina[i] = Math.sin(a);
				}
			}
			for (i = 0; i < this.accumulator.length; i+=1) {
				this.accumulator[i] = 0.0;
			}
			for (i = 0; i < this.aPxls.data.length; i+=4) {
				this.aPxls.data[i] =
				this.aPxls.data[i+1] =
				this.aPxls.data[i+2] = 0;
				this.aPxls.data[i+3] = 255;
			}
		},
		toggle_pause: function() {
			// user can click to stop/start the animation
			this.paused = ! this.paused;
		},
		main: function(canvID,srcImg) {
			// called once we have a source image
			var w = this.srcImage.width;
			var h = this.srcImage.height;
			this.canvas = document.getElementById(canvID);
			// adjust canvas width if on a narrow screen (e.g., a phone)
			var p = this.canvas.parentElement;
			if (p.offsetWidth < (this.canvas.width-4)) {
				this.canvas.width = p.offsetWidth - 4;
			}
			this.ctx = this.canvas.getContext('2d');
			this.canvas.onclick = this.toggle_pause.bind(this);
			// other canvases might not be visible, if the post is part of an index page
			this.canvasB = document.getElementById(canvID+'B');
			if (this.canvasB) {
				this.ctxB = this.canvasB.getContext('2d');
				if (!this.samplePtB) {
					this.samplePtB = [new bouncePt(64, 34,
											(Math.random()*5)-2.5,
											(Math.random()*5)-2.5,
											w,h,0,0,0.995),
										new bouncePt(64, 64,
											(Math.random()*5)-2.5,
											(Math.random()*5)-2.5,
											w,h,0,0,0.995)];
				}
			}
			this.canvasC = document.getElementById(canvID+'C');
			if (this.canvasC) {
				this.ctxC = this.canvasC.getContext('2d');
				if (!this.samplePtC) {
					this.samplePtC = new bouncePt(80, 64,
									(Math.random()*7)-3.5,
									(Math.random()*5)-2.5,
									this.aSize.x, this.aSize.y);
				}
			}
			this.canvasD = document.getElementById(canvID+'D');
			if (this.canvasD) {
				this.ctxD = this.canvasD.getContext('2d');
			}
			this.canvasE = document.getElementById(canvID+'E');
			if (this.canvasE) {
				this.ctxE = this.canvasE.getContext('2d');
			}
			if (!this.srcImage) {
				this.srcImage = document.getElementById('htSrcImg');
			}
			if (!this.hiddenCanvas) {
				this.hiddenCanvas = document.createElement('canvas');
				this.hiddenCanvas.width = w;
				this.hiddenCanvas.height = h;
			}
			var hiddenCtx = this.hiddenCanvas.getContext('2d');
			hiddenCtx.drawImage(this.srcImage, 0,0);
			this.collect_gray_pixels(hiddenCtx,0,0,w,h);
			this.startup();
		},
		load_new_image: function() {
			if (this.srcImage) {
				//this.srcImage.src = this.debug ? "/img/misc/c-scoot.jpg" : ('https://picsum.photos//128/128/?random&junk='+new Date().getTime());
				this.srcImage.src = ('https://picsum.photos//128/128/?random&junk='+new Date().getTime());
			} else {
				// how could we have gotten here???
				console.log("Warning! button pressed too soon?");
			}
		},
		begin: function() {
			var env = "{{ jekyll.environment }}";
			this.debug = (env == 'development');
			// load an interesting pic before we try anything else
			this.srcImage = new Image();
			this.srcImage.setAttribute('id','htSrcImg');
			this.srcImage.setAttribute('crossOrigin','anonymous');
			this.srcImage.onload = function() {this.main("km_sample5"); }.bind(this);
			this.load_new_image();
		}
	};
	htSample1.begin();
</script>

_Part of a brief series that started [here.]({{ site.baseurl}}{% post_url 2018-12-05-K-Means %})_
{: .notice--info}

The world is full of patterns: art, physics, weather, music. In some camps, it's been proposed that human ability to grasp at patterns is the core of conscious experience. We understand our surroundings, and one another, as a complex carpet of overlapping and expected patterns. The patterns might be hard to fathom at first, then... like the sudden twist in a mystery story, _Aha!_ the experience of understanding. It all makes sense now.

With that possibility for epiphany in mind, here's a paradoxical pattern to consider: _when is a point a line, and when is line a point, especially when that line is.. a curve? And how, may I ask, does this relate to beer?_

May we present: _The Hough Transform..._

<!--more-->
### The Idea

The story of the Hough Transform, like many great ideas, begins with our last question above: _the beer._

#### Bubbles

In 1952 physicist [Donald A. Glaser](https://mcb.berkeley.edu/news-and-events/department-news/glaser-bubble-chamber) invented the ["bubble chamber,"](https://en.wikipedia.org/wiki/Bubble_chamber) a device for tracing the motions of charged particles moving through a liquid. [_Some of the initial experiments used beer inside the chamber._](https://www2.lbl.gov/Publications/Currents/Archive/Jul-21-2006.html#6) Of course they did, this was at a university. But the only reported results from the high-energy beer physics experiments, alas, were some bad smells (seriously). Switching from beer to liquid hydrogen led to Glaser winning the Nobel Prize only eight years later, in 1960.

The year before that prize, 1959, another researcher, [Paul Hough,](https://www.gf.org/fellows/all-fellows/paul-v-c-hough/) published "Machine Analysis of Bubble Chamber Pictures," in the popular _Proceedings of the International Conference on High Energy Accelerators and Instrumentation._ If you missed that, he also filed a [1962 patent.](https://patents.google.com/patent/US3069654A/en)

The problem Hough was trying to get at was this: how to find straight lines in very low-quality, sketchy, incomplete pictures? Including ones where there's a lot of extra noise from spurious bubbles, low-quality signals (this is 1959, using a 1950's-era TV camera), and how to do it _quickly_ (at TV camera speed)?

<figure class="align-center">
	<img alt="Fermilab Bubble Chamber Image" src="{{ site.url}}/img/pix2018/fermilab.jpg">
  <figcaption>Not a Hokusai or da Vinci sketch, but a bubble chamber photo from Fermilab</figcaption>
</figure>

Breaking down our perceptions into patterns is a big part of what has been termed "artificial intelligence" at least as far back as those days in the 1950's. Hough, or you or I, could look at a picture with broken lines and instantly understand how the pieces connect into a single continuous stroke. But how might this perception be achieved automatically?

#### Lines Can Be Points

Hough had the genius idea to think of _Lines as Points._

Pictures like those on a TV camera or a computer display are made from points in X and Y dimensions: pixels, arranged across the width and height of the picture. Just like Mr. Evan's Anaylytic Geometry class back in Junior High. As Mr. Evans probably told you, in algebra any perfect, infinite line on an XY picture can be described by a little equation that you might recall from that classroom: _y&nbsp;=&nbsp;m&sdot;x&nbsp;+&nbsp;b_ where _m_ is the slope (how steep is the line?), and _b_ is the "y-intercept" where the line crosses the x=0 axis (that is: changing _b_ moves the line up and down on the chart).

Since a perfect infinite line doesn't care about the _specific_ x or y, we can just express it as _m_ and _b_ &mdash; which in turn could just be a single point _(m,b)_ on its own 2D graph, where, say, the horizontal axis is slope _m_ and the vertical one is _b._ Which is exactly what Hough did.

A problem with this description is that vertical lines, where the slope is infinite, can't be drawn. This worked out okay for Hough's specific case. In fact it worked _great_ for him. But it's kind of bad if we want to use this to look at most everyday pictures, which have lots of vertical lines relating to buildings, trees, and so forth.

A decade later, two researchers at Stanford, [Duda &amp; Hart,](http://www.ai.sri.com/pubs/files/tn036-duda71.pdf) realized you could just use a different equation for straight lines to fix that problem.

The equation they prefered uses an angle between zero and 180 degrees to define the line's _tilt_ (called _theta_ or &theta;),and the _radial distance_ of the line from the (0,0) origin (called _rho_ or &rho;). The &rho; distance for lines that appear in any real picture can never be farther than the distance between a picture's opposite corners.

The illustration below shows the relationships. Any line that can appear in the dark image box on the right is also represented by the single highlighted point in the white rectangle of (&theta;, &rho;) space on the left. You should be able to see how changing values of tilt and distance affect where the picture line is drawn.

The specific line equation for theta and rho, which you're not at all required to recall later, is: _x&sdot;cos&theta;&nbsp;+&nbsp;y&sdot;sin&theta;&nbsp;=&nbsp;&rho;_

<figure class="align-center">
	<canvas width="392" height="130" id="km_sample5C" class="align-center">
	</canvas>
  <figcaption id='htStatsC'>Similarly, any point on the &theta;/&rho; plane represents line in the image plane.</figcaption>
</figure>

#### Points and Lines to Curves

Okay, we have some interesting ways to write about straight lines, but how do we go from a _picture,_ made up of many single points, to finding meaningful straight lines that join those points? Especially since we don't know which points are "good" important points, and which ones are just other random parts of the picture?

First, Hough knew that the lines he cared about were along _edges_ in the picture. In the bubble chamber, a picture was contrasty black-and-white, like in the photo above (but lower quality). So he just watched for sudden changes in brightness along each tiny part of the TV signal, and only paid attention to places where the values jumped between light and dark or back. In our sample code, we use a pretty similar edge finder called the _Sobel filter,_ likewise looks for changes in neighboring pixel brightnesses and mixes edges found in both vertical and horizontal directions (we can talk about the specifics of Sobel in another post).

<figure class="align-center">
	<canvas width="392" height="130" id="km_sample5D" class="align-center">
	</canvas>
  <figcaption id='htStatsD'>Color image, grayscale, vertical/horizontal edges (marked in color)</figcaption>
</figure>

If a pixel isn't on an edge, we can ignore it -- a big time saver for some pictures.

So now that we've picked some points (pixels) as best candidates to be part of a line, let's consider just one pixel's worth: one isolated point. There's an infinite number of _possible_ lines that could pass through that point, with tilts ranging around a full 180 degrees. Lines tilting up, tilting sideways, anywhere between.

If we look at it in our &theta;/&rho; space, we can plot a graph of all the possible values of the tilts and corresponding distances for the possible lines that pass through that point, for every angle from -90 degrees to + 90 degrees. The resulting graph shape will show a smooth, sinuous curve of different possible distances, one for each angle we consider.

Time for another Hough insight.

#### Curves to Lines

Let's now consider a second point on the same picture. It too has an infinite number of possible lines in the picture that could cross it, but only _one_ of those lines will cross _both_ the first and second points.

If we add to our graph the same sort of curve described above for that new point, the two curves of possibilities will always cross one another, exactly _once,_ somewhere on our graph of &theta;/&rho; space. And where the curves cross is _exactly the point that describes the line between those two original picture points._

_Mind blown._ (ymmv)

<figure class="align-center">
	<canvas width="392" height="130" id="km_sample5B" class="align-center">
	</canvas>
  <figcaption id='htStatsB'>For any pixel on the image plane, there will be infinitely many possible combinations of &theta; (horizontal) and &rho; (vertical) that describe lines that pass through that point. Here we show the possible values between -90 and 90 degrees of &theta;, and where the two possiblities coincide &mdash; which will be at the &theta;/&rho; representation of the straight image-space line connecting the two points.</figcaption>
</figure>

#### Curves to Lines: Taking It to the Voters

You might be asking: _Couldn't we have found a line between two points in a much, much, easier way?_

Sure, for just two points. But what about for _lots_ of points?

Let's add a third point. It might not be along the same connecting line as those first two, but if it is, we have that much more reason to think that our line is a good one. Its curve will intersect in just the same spot. The more points that line up, the stronger our conviction that this particular combination of theta and rho is one we like. We can keep doing this for every candidate point.

So let's let every pixel vote!

To allow voting, we can just accumulate all the curves on one gray-scale graph of pixels, adding up the intersections for all possible tilt angles and for all the edge points in the picture. The brightest spot(s) in the graph will correspond to the strongest lines in the original picture.

In our sample, we've slowed the process down a bit for the benefit of human viewing -- in practice it's very quick, and even quicker if you can exploit a GPU (as is done, say, in the cameras of a self-driving car, but not done here).

<figure class="align-center">
	<canvas width="392" height="130" id="km_sample5E" class="align-center">
	</canvas>
  <figcaption id='htStatsE'>Scanning through each row, we find more and more curves that overlap. The brightest areas are those with the largest number of overlaps. On the original picture, we overlay the lines from the strongest values found-so-far as we scan from top to bottom.</figcaption>
</figure>

There you have it! Hough's brilliant idea.

* Lines can be points
* An infinite collection of possible lines can be a line (of lines)
  * ...or a curve (of lines)
* When drawn in a collection of pixel points, curves made from points can vote for the best lines that might connect those points

Since 1959 Hough's gymnastic transform been extended for many other uses and different sorts of shapes. His paper is one of the most-cited in all of computer vision.

I have no idea if it's been used to trace actual beer bubbles. Maybe someone should try, the effort could at least be entertaining while the test supplies last.

---

### The Example Code

As usual, all code here is contained within the HTML &mdash; no calling-out to external resources like _d3_ etc. Selecting "view source" or looking at [this github page](https://raw.githubusercontent.com/joker-b/botzo/master/_posts/2018-12-17-Hough.md) will reveal all. The single `<script>` tag contains the main demo object as `htSample1`

You may have noticed we skipped past the details of the _Sobel filter_ edge extraction. Extracting edges can be straightforward using _image convolutions._ Convolutions are cool but not the direct topic of this post, so I've borrowed and tweaked the convolution code [here from "html5rocks.com."](https://www.html5rocks.com/en/tutorials/canvas/imagefilters/) To make that clearer all such code is organized in a sub-object `htSample1.filters`

A nice attribute of using a raster image to collect the possibility curves is that it lets us precalculate all the potentially-slow _sin_ and _cos_ functions exactly once as the page loads. A big time saver.

This code loads random images from the internet. Every now and then, the image doesn't actually exist, and the code will halt. Just press the button near the top figure to skip forward to the next picture, and things will carry on again.