#!/usr/bin/env python
"""Put things in rectangular gaps in a submesh

The first argument is the submesh we are modifying.
You can then add one or more things to do.
(If zero things to do are given, the script will print the location 
of every vertex in the mesh.)

Arguments ``grid x y file`` will put a curved mesh in the gap based on 
the curvature of the edges. In this case, x any y should be the number 
of vertices across and down that the grid of vertices which we want to 
exist will have, including vertices on the edges and the corners of 
this grid which must already exist, and the file must have x*y lines 
specifying which vertices of this grid already exist and where they are, 
arranged in (Western) reading order, like:

1 2 3 4
5 6 7 8

For vertices which alrady exist, a line should be three floating point 
numbers given its location, which specify an area based on the number of 
decimal points, e.g. 4.95 gives anything from 4.945 to 4.955, 4.950 gives 
anything from 4.9495 to 4.9505. If the script needs to create the vertex, 
any non-empty line without numbers works, such as '-'

Arguments nipple file will put a nipple-like bump in the gap. The gap must 
be either 2x2 vertices (so all four original vertices must already exist) 
or 3x3 vertices (so only the centre one can be missing.) The file should 
first list three locations, and then list the vertex grid as for 'grid' 
arguments. The first location should be a notional point inside the nipple, 
for example, in a 3x3 vertex grid, it might be where the centre vertex would 
otherwise be, the second location is the tip of the nipple, and the third 
is a skin location on the edge of the nipple, conventionally directly above 
its centre.
"""

import sys
import copy
import math

import vb
import vector

def dump_locations(path):
	mesh = vb.Mesh(path)
	for v in mesh.vertices:
		print (' '.join(['%0.6f' % x for x in v.pos()]))
	sys.exit(0)

def read_vertex(mesh, grids, line):
	assert(line.strip())
	bits = line.strip().split()
	if bits[0] in grids:
		x, y = [int(x) for x in bits[1:]]
		across, down, vertices = grids[bits[0]]
		pos = y * across + x
		return vertices[pos]
	elif set('0123456789').intersection(set(line)):
		return mesh.find_by_position(bits)
	else:
		return None

def read_location(line, grids):
	bits = line.strip().split()
	if bits[0] in grids:
		x, y = [int(x) for x in bits[1:]]
		across, down, vertices = grids[bits[0]]
		pos = y * across + x
		return vertices[pos][1].pos()
	return [float(x) for x in bits]

def curvature(v, models):
	corners_a = [models[0].pos(), models[1].pos()]
	corners_b = [models[2].pos(), models[3].pos()]
	upcurved = [models[x].pos() for x in range(4, 8)]
	dirv = vector.cross(vector.minus(corners_a[0], corners_a[1]), 
			vector.minus(corners_b[0], corners_b[1]))
	base_point = vector.average(corners_a + corners_b)
	entended_point = vector.average(upcurved)
	deviation = vector.dot(vector.minus(entended_point, base_point), dirv) * 2
	should_be = deviation + vector.dot(base_point, dirv)
	current = vector.dot(v.pos(), dirv)
	v.setpos(vector.add([v.pos(), 
			vector.divide(
				vector.mult(dirv, should_be - current), 
			vector.vlen2(dirv))]))

def do_grid(mesh, existing, across, down, spec):
	with open(spec, 'r') as fi:
		vertices = [read_vertex(mesh, existing, line) for line in fi if line.strip()]
	if across * down != len(vertices):
		raise Exception(f'Wrong vertex list length, {across}*{down}!={len(vertices)}')
	for y in range(1, down - 1):
		for x in range(1, across - 1):
			pos = y * across + x
			if vertices[pos] is not None:
				continue
			references = [(x, vertices[y * across]), 
					(across - x - 1, vertices[(y + 1) * across - 1]),
					(y, vertices[x]),
					(down - y - 1, vertices[x + (across * (down - 1))])]
			for i in references:
				if i[1] is None:
					raise Exception(f'Missing edge or corner vertex for {x} {y}')
			v = vb.new_vertex([(w[1], 1/j) for j, w in references])
			curvature(v, [vertices[x][1] for x in [0, across * down - 1, across - 1, across * (down - 1), y * across, (y + 1) * across - 1, x, x + (across * (down - 1))]])
			vertices[pos] = (len(mesh.vertices), v)
			mesh.vertices.append(v)
	num_originals = mesh.lv
	mesh.lv = len(mesh.vertices)
	need = []
	for y in range(down - 1):
		for x in range(across - 1):
			corner = y * across + x
			need.append((corner, corner + across, corner + 1))
			need.append((corner + 1, corner + across, corner + across + 1))
	need = [[vertices[i][0] for i in face] for face in need]
	for face in need:
		if max(face) < num_originals:
			if face not in mesh.faces:
				mesh.faces.append(face)
		else:
			mesh.faces.append(face)
	mesh.lf = len(mesh.faces)
	return (across, down, vertices)

def clear_square(mesh, svs):
	mapping = mesh.clear_window([x[0] for x in svs])
	retm = None
	if len(svs) == 9:
		retm = mapping
		if svs[4][0] in mapping:
			raise Exception('Failed to delete centre vertex')
		retm[svs[4][0]] = None
		retm[None] = None
		svs = svs[:4] + svs[5:]
	ret = []
	for i, v in svs:
		if i not in mapping:
			raise Exception('Deleted side vertex!')
		ret.append((mapping[i], v))
	return ret, retm

def do_nip(mesh, grids, spec):
	with open(spec, 'r') as fi:
		basepoint = read_location(fi.readline(), grids)
		tip = read_location(fi.readline(), grids)
		firstjoin = read_location(fi.readline(), grids)
		svs = [read_vertex(mesh, grids, line) for line in fi if line.strip()]
	svs, mapping = clear_square(mesh, svs)
	if len(svs) == 8:
		models = [svs[i][1] for i in (1, 3, 4, 6)]
	else:
		models = [x[1] for x in svs]
		insert = [(mesh.lv, mesh.bisect_line(svs[0][0], svs[1][0]))]
		insert.append((mesh.lv, mesh.bisect_line(svs[0][0], svs[2][0])))
		insert.append((mesh.lv, mesh.bisect_line(svs[1][0], svs[3][0])))
		insert.append((mesh.lv, mesh.bisect_line(svs[2][0], svs[3][0])))
		svs = [svs[0], insert[0], svs[1], insert[1], insert[2], svs[2], insert[3], svs[3]]
	model = vb.new_vertex([(x, 1) for x in models])
	tex_centre = model.tex()
	tex_radius = min([vector.distance(tex_centre, x[1].tex()) for x in [svs[1], svs[3], svs[4], svs[6]]])
	print (f'Textures centred {tex_centre} radius {tex_radius}')
	tr1 = tex_radius * 3.0 / 4
	tr2 = tex_radius / 2.0
	tr3 = tex_radius / 4.0
	tex_vector = vector.unit(vector.minus(svs[1][1].tex(), tex_centre))
	tex_sidevec = [tex_vector[1], -tex_vector[0]]
	outvec = vector.minus(tip, basepoint)
	outmax = vector.vlen(outvec)
	unitout = vector.divide(outvec, outmax)
	yvec = vector.minus(firstjoin, basepoint)
	outmin = vector.dot(unitout, yvec)
	outbig = 0.85 * outmax + 0.15 * outmin
	outsmall = 0.1 * outmax + 0.9 * outmin
	yvec = vector.minus(yvec, vector.mult(unitout, outmin))
	xvec = vector.unit(vector.cross(yvec, outvec))
	xvec = vector.mult(xvec, vector.vlen(yvec))
	new_vertices = [model]
	s = math.sqrt(0.75)
	ydegrees = [1, s, 0.5, 0, -0.5, -s, -1, -s, -0.5, 0, 0.5, s]
	xdegrees = [0, 0.5, s, 1, s, 0.5, 0, -0.5, -s, -1, -s, -0.5]
	for i in range(36):
		new_vertices.append(copy.copy(model))
		new_vertices[-1].values = list(model.values)
	for i in range(12):
		new_vertices[i].setpos(vector.combine([basepoint, outmin, unitout, ydegrees[i], yvec, xdegrees[i], xvec]))
		new_vertices[i].settex(vector.combine([tex_centre, tr1 * ydegrees[i], tex_vector, tr1 * xdegrees[i], tex_sidevec]))
	for i in range(12, 24):
		new_vertices[i].setpos(vector.combine([basepoint, outsmall, unitout, 0.55 * ydegrees[i%12], yvec, 0.55 * xdegrees[i%12], xvec]))
		new_vertices[i].settex(vector.combine([tex_centre, tr2 * ydegrees[i%12], tex_vector, tr2 * xdegrees[i%12], tex_sidevec]))
	for i in range(24, 36):
		new_vertices[i].setpos(vector.combine([basepoint, outbig, unitout, 0.45 * ydegrees[i%12], yvec, 0.45 * xdegrees[i%12], xvec]))
		new_vertices[i].settex(vector.combine([tex_centre, tr3 * ydegrees[i%12], tex_vector, tr3 * xdegrees[i%12], tex_sidevec]))
	new_vertices[36].setpos(tip)
	new_faces = [
			[svs[0][0], svs[3][0], mesh.lv + 9], 
			[svs[0][0], mesh.lv + 9, mesh.lv + 10], 
			[svs[0][0], mesh.lv + 10, mesh.lv + 11], 
			[svs[0][0], mesh.lv + 11, mesh.lv + 0], 
			[svs[0][0], mesh.lv + 0, svs[1][0]], 
			[svs[2][0], svs[1][0], mesh.lv + 0], 
			[svs[2][0], mesh.lv + 0, mesh.lv + 1], 
			[svs[2][0], mesh.lv + 1, mesh.lv + 2], 
			[svs[2][0], mesh.lv + 2, mesh.lv + 3], 
			[svs[2][0], mesh.lv + 3, svs[4][0]], 
			[svs[5][0], mesh.lv + 9, svs[3][0]], 
			[svs[5][0], mesh.lv + 8, mesh.lv + 9], 
			[svs[5][0], mesh.lv + 7, mesh.lv + 8], 
			[svs[5][0], mesh.lv + 6, mesh.lv + 7], 
			[svs[5][0], svs[6][0], mesh.lv + 6], 
			[svs[7][0], svs[4][0], mesh.lv + 3], 
			[svs[7][0], mesh.lv + 3, mesh.lv + 4], 
			[svs[7][0], mesh.lv + 4, mesh.lv + 5], 
			[svs[7][0], mesh.lv + 5, mesh.lv + 6], 
			[svs[7][0], mesh.lv + 6, svs[6][0]], 
	]
	for k in [0, 12]:
		for i in range(11):
			new_faces.extend([
					[mesh.lv + k + i + 12, mesh.lv + k + i + 1, mesh.lv + k + i + 0],
					[mesh.lv + k + i + 12, mesh.lv + k + i + 13, mesh.lv + k + i + 1], 
					])
		new_faces.extend([
				[mesh.lv + k + 23, mesh.lv + k + 0, mesh.lv + k + 11],
				[mesh.lv + k + 23, mesh.lv + k + 12, mesh.lv + k + 0], 
				])
	for i in range(11):
		   new_faces.append([mesh.lv + 36, mesh.lv + i + 25, mesh.lv + i + 24])
	new_faces.append([mesh.lv + 36, mesh.lv + 24, mesh.lv + 35])
	mesh.vertices.extend(new_vertices)
	mesh.lv += 37
	mesh.faces.extend(new_faces)
	mesh.lf += len(new_faces)
	return mapping

if __name__ == '__main__':
	if len(sys.argv) == 2:
		dump_locations(sys.argv[1])
	mesh = vb.Mesh(sys.argv[1])
	grids = {}
	pos = 2
	while len(sys.argv) > pos + 1:
		if sys.argv[pos] == 'grid':
			grids[sys.argv[pos + 3]] = do_grid(mesh, grids, int(sys.argv[pos + 1]), int(sys.argv[pos + 2]), sys.argv[pos + 3])
			pos += 4
		elif sys.argv[pos] == 'nipple':
			mapping = do_nip(mesh, grids, sys.argv[pos + 1])
			pos += 2
			if mapping is not None:
				for k, grid in grids.items():
					grids[k] = (grid[0], grid[1], [(mapping[i], v) for i, v in grid[2]])
		else:
			raise Exception('Unknown operation ' + sys.argv[pos])
	mesh.write(sys.argv[-1])


