index-project.lua

local modname = ...
local index = package.modtable(modname)
local ffi = require('ffi')
local databuffer = require('databuffer')
local filestorage = require('rsbin.filestorage')
local fileopenmode = require('rsbin.fileopenmode')
local storagestream = require('rsbin.storagestream')
local pngreader = require('rsbin.pngreader')
local pngwriter = require('rsbin.pngwriter')
local bitmapformat = require('rsbin.bitmapformat')
local invoke = require('base.invoke')

local function readpng(path, format)
	local storage = filestorage:create(path, fileopenmode.read)
	local stream = storagestream:create(storage, true, false, 0)
	local reader = pngreader:create(stream, bitmapformat[format])
	while not reader:poll() do
		coroutine.yield()
	end
	local bitmap = {
		pixels = reader:getpixels(),
		width = reader:getwidth(),
		height = reader:getheight()}
	reader:release()
	stream:release()
	storage:close()
	storage:release()
	return bitmap
end

local function writepng(path, format, bitmap)
	local storage = filestorage:create(path, fileopenmode.create)
	local stream = storagestream:create(storage, false, true, 0)
	local writer = pngwriter:create(
		stream, bitmapformat[format], bitmap.width, bitmap.height, bitmap.pixels)
	while not writer:poll() do
		coroutine.yield()
	end
	writer:release()
	stream:release()
	storage:close()
	storage:release()
end

local function createbitmap(width, height)
	width = math.ceil(width)
	if width < 1 then
		width = 1
	end
	height = math.ceil(height)
	if height < 1 then
		height = 1
	end
	local bitmap = {
		width = width,
		height = height}
	bitmap.pixels = databuffer:create(4*width*height, 4*width*height, nil)
	return bitmap
end

local function urotator(angle, z0, sscale, tscale)
	local cosphi = math.cos(angle)
	local sinphi = math.sin(angle)
	return function(s, t)
		s = s / sscale
		t = t / tscale
		local uvden = z0 / (t * sinphi + z0 * cosphi)
		local u = uvden * cosphi * s
		local v = uvden * t
		local stden = 1 / (z0 - v * sinphi)
		local dsdu = sscale * stden * z0
		local dsdv = sscale * stden*stden * u * (z0 * sinphi)
		local dtdu = 0
		local dtdv = tscale * stden*stden * (z0*z0 * cosphi)
		return u, v, dsdu, dsdv, dtdu, dtdv
	end
end

local function vrotator(angle, z0, sscale, tscale)
	local cosphi = math.cos(angle)
	local sinphi = math.sin(angle)
	return function(s, t)
		s = s / sscale
		t = t / tscale
		local uvden = z0 / (s * sinphi + z0 * cosphi)
		local u = uvden * s
		local v = uvden * cosphi * t
		local stden = 1 / (z0 - u * sinphi)
		local dsdu = sscale * stden*stden * (z0*z0 * cosphi)
		local dsdv = 0
		local dtdu = tscale * stden*stden * v * (z0 * sinphi)
		local dtdv = tscale * stden * z0
		return u, v, dsdu, dsdv, dtdu, dtdv
	end
end

local function lineartransform(su, sv, tu, tv)
	local ddet = 1 / (su*tv - sv*tu)
	local us = ddet * tv
	local vs = - ddet * tu
	local ut = - ddet * sv
	local vt = ddet * su
	return function(s, t)
		local u = us * s + ut * t
		local v = vs * s + vt * t
		local dsdu = su
		local dsdv = sv
		local dtdu = tu
		local dtdv = tv
		return u, v, dsdu, dsdv, dtdu, dtdv
	end
end

local function lanczos2(x)
	x = 0.25*x*x
	if x > 1 then
		return 0
	end
	local xd = 2*x
	local b5 = -0.3312768593109787
	local b4 = 2.5271939547407243 + xd * b5
	local b3 = -9.308061847706597 + xd * b4 - b5
	local b2 = 21.490018064632473 + xd * b3 - b4
	local b1 = -34.340697422247366 + xd * b2 - b3
	local r = 19.962824109891756 + x * b1 - b2
	return r
end

local function matrixinverse(dsdu, dsdv, dtdu, dtdv)
	local ddet = 1 / (dsdu*dtdv - dsdv*dtdu)
	local duds = ddet * dtdv
	local dudt = - ddet * dtdu
	local dvds = - ddet * dsdv
	local dvdt = ddet * dsdu
	return duds, dudt, dvds, dvdt
end

local function saturatecolor(a)
	return 0.5 * (math.abs(a)-math.abs(a-255)+255)
end

local function basisextent(duds, dudt, dvds, dvdt)
	local dumin = math.min(duds+dudt, duds-dudt, -duds+dudt, -duds-dudt)
	local dvmin = math.min(dvds+dvdt, dvds-dvdt, -dvds+dvdt, -dvds-dvdt)
	local dumax = math.max(duds+dudt, duds-dudt, -duds+dudt, -duds-dudt)
	local dvmax = math.max(dvds+dvdt, dvds-dvdt, -dvds+dvdt, -dvds-dvdt)
	return dumin, dvmin, dumax, dvmax
end

local function sample_decimate(
	bitmap, data,
	u, v,
	dsdu, dsdv, dtdu, dtdv,
	duds, dudt, dvds, dvdt
)
	local dumin, dvmin, dumax, dvmax = basisextent(duds, dudt, dvds, dvdt)
	local minu = math.floor(u + math.max(-10, 2*dumin))
	local minv = math.floor(v + math.max(-10, 2*dvmin))
	local maxu = math.ceil(u + math.min(10,  2*dumax))
	local maxv = math.ceil(v + math.min(10,  2*dvmax))
	local totalr = 0
	local totalg = 0
	local totalb = 0
	local totalrgbweight = 0
	local totala = 0
	local totalaweight = 0
	for pv = minv, maxv do
		for pu = minu, maxu do
			local du = pu - u
			local dv = pv - v
			local ds = dsdu*du + dsdv*dv
			local dt = dtdu*du + dtdv*dv
			local aweight = lanczos2(ds) * lanczos2(dt)
			if
				pu >= 0 and pu < bitmap.width and
				pv >= 0 and pv < bitmap.height
			then
				local pixel = data + 4 * bitmap.width * pv + 4 * pu
				local rgbweight = aweight * pixel[3]
				totalr = totalr + pixel[0] * rgbweight
				totalg = totalg + pixel[1] * rgbweight
				totalb = totalb + pixel[2] * rgbweight
				totala = totala + pixel[3] * aweight
				totalrgbweight = totalrgbweight + rgbweight
			end
			totalaweight = totalaweight + aweight
		end
	end
	if totalrgbweight ~= 0 then
		totalr = totalr / totalrgbweight
		totalg = totalg / totalrgbweight
		totalb = totalb / totalrgbweight
	end
	if totalaweight ~= 0 then
		totala = totala / totalaweight
	end
	return
		saturatecolor(totalr),
		saturatecolor(totalg),
		saturatecolor(totalb),
		saturatecolor(totala)
end

local function sample_interpolate(
	bitmap, data,
	u, v,
	dsdu, dsdv, dtdu, dtdv
)
	local minu = math.floor(u - 2)
	local minv = math.floor(v - 2)
	local maxu = math.ceil(u + 2)
	local maxv = math.ceil(v + 2)
	local totalr = 0
	local totalg = 0
	local totalb = 0
	local totalrgbweight = 0
	local totala = 0
	local totalaweight = 0
	for pv = minv, maxv do
		for pu = minu, maxu do
			local du = pu - u
			local dv = pv - v
			local aweight = lanczos2(du) * lanczos2(dv)
			if
				pu >= 0 and pu < bitmap.width and
				pv >= 0 and pv < bitmap.height
			then
				local pixel = data + 4 * bitmap.width * pv + 4 * pu
				local rgbweight = aweight * pixel[3]
				totalr = totalr + pixel[0] * rgbweight
				totalg = totalg + pixel[1] * rgbweight
				totalb = totalb + pixel[2] * rgbweight
				totala = totala + pixel[3] * aweight
				totalrgbweight = totalrgbweight + rgbweight
			end
			totalaweight = totalaweight + aweight
		end
	end
	if totalrgbweight > 1 then
		totalr = totalr / totalrgbweight
		totalg = totalg / totalrgbweight
		totalb = totalb / totalrgbweight
	end
	if totalaweight ~= 0 then
		totala = totala / totalaweight
	end
	return
		saturatecolor(totalr),
		saturatecolor(totalg),
		saturatecolor(totalb),
		saturatecolor(totala)
end

local function sample(bitmap, data, u, v, dsdu, dsdv, dtdu, dtdv)
	local duds, dudt, dvds, dvdt = matrixinverse(dsdu, dsdv, dtdu, dtdv)
	local dssqr = duds*duds + dvds*dvds
	local dtsqr = dudt*dudt + dvdt*dvdt
	if dssqr > 1 and dtsqr > 1 then
		return sample_decimate(
			bitmap, data,
			u, v,
			dsdu, dsdv, dtdu, dtdv,
			duds, dudt, dvds, dvdt)
	else
		return sample_interpolate(
			bitmap, data,
			u, v,
			dsdu, dsdv, dtdu, dtdv)
	end
end

local function reproject(
	source, transform,
	twidth, theight,
	soffset, toffset, uoffset, voffset,
	tag
)
	local target = createbitmap(twidth, theight)
	local sdata = source.pixels:getdata()
	local tdata = target.pixels:getdata()
	soffset = soffset + 0.5 * (target.width - 1)
	toffset = toffset + 0.5 * (target.height - 1)
	uoffset = uoffset + 0.5 * (source.width - 1)
	voffset = voffset + 0.5 * (source.height - 1)
	for y = 0, target.height-1 do
		print(tag, y, '/', target.height)
		local row = tdata + 4 * target.width * y
		for x = 0, target.width-1 do
			local pixel = row + 4 * x
			local s = x - soffset
			local t = y - toffset
			local u, v, dsdu, dsdv, dtdu, dtdv = transform(s, t)
			local r, g, b, a = sample(
				source, sdata,
				u + uoffset, v + voffset, dsdu, dsdv, dtdu, dtdv)
			pixel[0] = r
			pixel[1] = g
			pixel[2] = b
			pixel[3] = a
		end
	end
	return target
end

local function main()
	local source = readpng('local/projsource.png', 'rgba8')
	for zf = 1, 5 do
		for i = -19, 19 do
			local angle = i/40*math.pi
			local z0 = math.max(source.width, source.height) * zf
			local maxs = math.cos(angle) * z0 * source.width
				/ (2 * z0 + math.sin(angle) * source.width)
			local mins = - math.cos(angle) * z0 * source.width
				/ (2 * z0 - math.sin(angle) * source.width)
			local maxt = z0 * source.height
				/ (2 * z0 - math.abs(math.sin(angle) * source.width))
			local sscale = source.width / (maxs-mins)
			local scaled
			if math.abs(angle) < 0.3*math.pi then
				scaled = reproject(
					source,
					vrotator(angle, z0, 1, 1),
					maxs-mins, 2*maxt,
					0.5*(maxs+mins), 0, 0, 0,
					zf..'|'..i)
			else
				local rotated = reproject(
					source,
					vrotator(angle, z0, 3*sscale, 3),
					sscale*(6+3*(maxs-mins)), 6+6*maxt,
					sscale*1.5*(maxs+mins), 0, 0, 0,
					zf..'|'..i..'r')
				scaled = reproject(
					rotated,
					lineartransform(1/(3*sscale), 0, 0, 1/3),
					maxs-mins, 2*maxt,
					0, 0, 0, 0,
					zf..'|'..i..'s')
			end
			writepng(
				string.format('local/proj%i-%.2i.png', zf, i+20), 'rgba8', scaled)
		end
	end
end

invoke(main)