A Digital Differential Analyzer-based algorithm for tracing circles on hexagonal grids

hex_circle.m

function [as, bs, x, y] = hex_circle(r)
	%% a,b represents a complex number in the form a + b * zeta, where zeta is the principal third root of unity.
	%% therefore, abs(a + b * zeta) = a^2 - a*b + b^2, which gives the DDA equations
	a = round(r);
	b = 0;
	as = [a];
	bs = [b];
	%% initialize DDA
	e = a^2 - r;
	da = 2 * a + 1 - b;
	db = 2 * b + 1 - a;

	while 2 * a > b
		e1 = e - da;
		e2 = e + db;
		e3 = e + da + db - 1;
		[~,w] = min(abs([e1, e2, e3]));
		switch w
			case 1
				e = e1;
				da -= 2;
				db += 1;
				a  -= 1;
			case 2
				e = e2;
				db += 2;
				da -= 1;
				b  += 1;
			case 3
				e = e3;
				db += 1;
				da += 1;
				a  += 1;
				b  += 1;
		end
		as(end+1) = a;
		bs(end+1) = b;
	end
	
	%% reflections 
	ar = as(end-1:-1:2);
	br = bs(end-1:-1:2);
	as = [as, - ar + br, -as, ar - br];
	bs = [bs,        br, -bs,    - br];

	%% to cartesian coordinates
	x = as - 0.5 * bs;
	y = bs * sqrt(3) / 2;
end