Biome detection algorithm?

Wow this is very informative! Thank you for all the effort , I’m working right now but will read through in a few hours. :slight_smile:

1 Like

You’re welcome! BTW, while I was working on this function I realized that it is important to note that since I’m using exponential decay, the resulting functions won’t be usable to compute a weighted area that the player’s position lies in. (Integration over an asymptote is a no-no). But, you can find the area immediately around the precise point. You’ll just need to add a check into the algorithm for determining the area of a biome. If the biome contains the player’s position then when you will need to split it somehow… I’d just draw a line through the biome. Or you could just say the position is 100% that biome and be done with it!

It has been a long couple days. Unfortunately I couldn’t find a way to solve those integrals nicely, and the many online calculators fried/broke trying to calculate them. It was difficult to find a Taylor series too, so I’ve included an implementation of Keplar’s rule. You can use a rather low number of segments to approximate the integrals since only the ratio of the integrals matters and the application doesn’t require a high degree of precision. For the most optimal precision/speed ratio I’d increase the number of segments as the radius difference increases (deltaR), but for the sake of speed I’d leave it as it is.

Please note that I’ve unit tested the Kepler’s rule implementation, and the pie arc length. The weight is so simple I felt it didn’t need tests, and the other arc lengths are basic mutations of the pie arc length, so I didn’t test them either. If anyone sees any errors or gets funky results, please let me know!

Here it is!

local acos = math.acos
local sqrt = math.sqrt
local clamp = math.clamp

--[[
	Unfortunately experimentation has shown that some values that should be equilvilant
	arn't quite so. I assumed that it was a precision error from .Magnitude, but it may
	also be that my equation for points a and b are incorrect. My intuition tells me that
	as the radius increases, we linearly move from p to p + offset[A|B], but I have not
	priven this. This slight error is why we use clamp.
	
	Also the integrals of these functions got rather complex. integral-calculator.com, 
	symbolab.com, and wolframalpha.com now hate me. Also, the derivatives are ugly, so
	a Taylor series was out of the question. I've settled for Kepler's rule (or simpson's
	if you prefer to name it after the person who stole the idea and published it).
 
	~IdiomicLanguage
--]]

local function keplersRule(f, n)
	assert(n % 2 == 0 and n > 0, "Kepler's rule requires an even number of segments")
	local sum = 0
	local dX = 1/n
	for i = 1, n - 1, 2 do
		sum = sum + f(dX * i)
	end
	sum = sum * 2
	for i = 2, n - 2, 2 do
		sum = sum + f(dX * i)
	end
	return (f(0) + sum * 2 + f(1)) / (3 * n)
end

--[[ Expectations:
	offsetA and offsetB are equidistant from center
	p, offsetA, and offsetB all exist on the same plane
	weight is a function that accepts any positive non-zero radius
	center, nor its projection, lays inside of p-offsetA-offsetB
	center's projection onto the plane is on the same side of the line
		offsetA-offsetB as p. (this will only be an issue if p isn't
		closer to the center than offsetA or offsetB)
	offsetA and offsetB are offsets from p. Their magnitudes
		are both non-zero
	I guess in general, don't do stupid
]]

local function weightedPieArcLength(weight, center, p, offsetA, offsetB)
	-- Welcome, variables. This cycle we are operating on a crossection of an
	-- infintely large sphere, or at least the projection of it onto our plane.
	local norm = offsetA:Cross(offsetB)
	center = center - norm * norm.Unit:Dot(center - p)
	local minR = (p - center).Magnitude
	local deltaR = (p + offsetA - center).Magnitude - minR
	return function(r)
		local a = offsetA * r + p - center
		local b = offsetB * r + p - center
		r = r * deltaR + minR
		return weight(r) * r * acos(clamp(a:Dot(b) / (r * r), -1, 1))
	end
end

-- same expectations as weightedPieArcLength but pa and pb act like
-- p, and must be equidistant from the center.
local function weightedTubeArcLength(weight, center, pa, pb, offsetA, offsetB)
	local norm = offsetA:Cross(offsetB)
	center = center - norm * norm.Unit:Dot(center - p)
	local minR = (pa - center).Magnitude
	local deltaR = (pa + offsetA - center).Magnitude - minR
	local minTheta = acos(clamp(pa:Dot(pb) / (minR * minR), -1, 1))
	return function(r)
		local a = offsetA * r + pa - center
		local b = offsetB * r + pb - center
		r = r * deltaR + minR
		return weight(r) * r * (minTheta + acos(clamp(a:Dot(b) / (r * r), -1, 1)))
	end
end

-- p must be farther from center than the offsets.
local function weightedInvertedPieArcLength(weight, center, p, offsetA, offsetB)
	local norm = offsetA:Cross(offsetB)
	center = center - norm * norm.Unit:Dot(center - p)
	local maxR = (p - center).Magnitude
	local deltaR = maxR - (p + offsetA - center).Magnitude
	return function(r)
		local a = offsetA * r + p - center
		local b = offsetB * r + p - center
		r = r * deltaR + maxR
		return weight(r) * r * acos(clamp(a:Dot(b) / (r * r), -1, 1))
	end
end

local n = 0.5 -- Weight decay ratio. 0 = no weight (distance doesn't matter)
local function weight(r)
	assert(r > 0, 'radius is always positive, and unless we want to blow up the universe, non-zero')
	return 1 / r ^ n
end

-- This is how to use one of these functions. I've named the arguments for your convenience.
local center = Vector3.new(0, 0, 0)
local piePoint = Vector3.new(10, 10, 0)
local offsetA = Vector3.new(0, 1, 0)
local offsetB = Vector3.new(1, 0, 0)

local f = weightedPieArcLength(weight, center, piePoint, offsetA, offsetB)
print(keplersRule(f, 10))
3 Likes