Circumcircle and incircle code + test from article

circles.lua

--
-- Permission is hereby granted, free of charge, to any person obtaining
-- a copy of this software and associated documentation files (the
-- "Software"), to deal in the Software without restriction, including
-- without limitation the rights to use, copy, modify, merge, publish,
-- distribute, sublicense, and/or sell copies of the Software, and to
-- permit persons to whom the Software is furnished to do so, subject to
-- the following conditions:
--
-- The above copyright notice and this permission notice shall be
-- included in all copies or substantial portions of the Software.
--
-- THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
-- EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
-- MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
-- IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
-- CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
-- TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
-- SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
--
-- [ MIT license: http://www.opensource.org/licenses/mit-license.php ]
--

local sqrt = math.sqrt

local function Circumcircle (px, py, qx, qy, rx, ry)
	-- Squared side lengths.
	local A2 = (qx - px)^2 + (qy - py)^2
	local B2 = (rx - qx)^2 + (ry - qy)^2
	local C2 = (px - rx)^2 + (py - ry)^2

	-- Center lines.
	local D, E = px - rx, py - ry
	local F = qx * D + qy * E + (A2 - B2) / 2
	local G, H = qx - rx, qy - ry
	local I = px * G + py * H + (A2 - C2) / 2

	-- Solve for the center.
	local cx, cy = (E * I - H * F) / (E * G - H * D)

	if 1 + E^2 == 1 then
		cy = (I - G * cx) / H
	else
		cy = (F - D * cx) / E
	end

	-- Solve for the radius.
	return cx, cy, sqrt((px - cx)^2 + (py - cy)^2)
end

local function Incircle (px, py, qx, qy, rx, ry)
	-- Side lengths.
	local A = sqrt((qx - px)^2 + (qy - py)^2)
	local B = sqrt((rx - qx)^2 + (ry - qy)^2)
	local C = sqrt((px - rx)^2 + (py - ry)^2)

	-- Partial side lengths.
	local Ap = (A + B - C) / 2
	local Bp = (B + C - A) / 2

	-- Find the triangle points on the sides.
	local a1, b1 = Ap / A, Bp / B
	local a2, b2 = 1 - a1, 1 - b1
	local x1, y1 = a1 * px + a2 * qx, a1 * py + a2 * qy
	local x2, y2 = b1 * qx + b2 * rx, b1 * qy + b2 * ry

	-- Find the midpoint from one of those points to the next. Find the direction from one of the
	-- original points toward the center. Get the leg lengths of the right triangle described by
	-- the side point, the midpoint, and the original point.
	local mx, my = (x1 + x2) / 2, (y1 + y2) / 2
	local hx, hy = mx - x1, my - y1
	local vx, vy = mx - qx, my - qy
	local H, V = sqrt(hx^2 + hy^2), sqrt(vx^2 + vy^2)

	-- Find the radius. Use it to scale the direction found in the previous step to reach the
	-- center. Optionally, compare to the circumcircle built from the side points.
	local seca = Ap / V -- sec(alpha)
	local r = H * seca -- r = A' * tan(alpha)
	local scale = seca^2 -- q-to-center : q-to-midpoint = r * sec(a) / V
	local cx, cy = qx + scale * vx, qy + scale * vy

	if true then
		local Cp = (A + C - B) / 2
		local c1 = Cp / C
		local c2 = 1 - c1
		local x3, y3 = c1 * rx + c2 * px, c1 * ry + c2 * py

		print("Compare", cx, cy, r, Circumcircle(x1, y1, x2, y2, x3, y3))
	end

	return cx, cy, r
end

local function RandomPoint ()
	local cx, cy = display.contentCenterX, display.contentCenterY

	return math.random(cx - 150, cx + 150), math.random(cy - 150, cy + 150)
end

local Dist, Mode = 10, "incircle"

local PX, PY = RandomPoint()
local QX, QY, RX, RY

repeat
	QX, QY = RandomPoint()
until (QX - PX)^2 + (QY - PY)^2 > Dist^2

repeat
	RX, RY = RandomPoint()
until (RX - PX)^2 + (RY - PY)^2 > Dist^2 and (RX - QX)^2 + (RY - QY)^2 > Dist^2

local func = Mode == "incircle" and Incircle or Circumcircle
local cx, cy, r = func(PX, PY, QX, QY, RX, RY)
local circ = display.newCircle(cx, cy, r)

local function Line (x1, y1, x2, y2)
	local line = display.newLine(x1, y1, x2, y2)

	line:setStrokeColor(math.random(), math.random(), math.random())

	line.strokeWidth = 2
end

Line(PX, PY, QX, QY)
Line(QX, QY, RX, RY)
Line(RX, RY, PX, PY)

display.newCircle(PX, PY, 5):setFillColor(1, 0, 0)
display.newCircle(QX, QY, 5):setFillColor(0, 1, 0)
display.newCircle(RX, RY, 5):setFillColor(0, 0, 1)

circ:setFillColor(0, 0)
circ:setStrokeColor(.2, .7, .3)

circ.strokeWidth = 3