/******************************************************************************************/
/*
/* Generates quartic curvature approximations and uv scale factors on an editable mesh
/* Bakes those values into channels 5 and 6 for use with the curved relief mapping shader.
/*
/* Based on the brilliant work of Fabio Policarpo 
/* Author: Drew Skillman
/* drew@drewskillman.com
/* 1/10/2006
/*
/* Don't blame me for any alien abductions that follow the use of this or any parts of this 
/* script. They were tested in a building that had very low ceilings, so there's almost 
/* certainly no way an alien spaceship could get in there.  If you are in a warehouse, or
/* perhaps outdoors on a laptop, I suggest you wrap foil or a wet towel around your head.
/******************************************************************************************/

/*********************************************************************/
/* Return the texture vertex index of the passed in face and vert
/*********************************************************************/
fn getTVfromV obj f v =
	(
	faceVerts = getface obj f
	
	fi = 0
	for i = 1 to 3 do
		(
		if faceVerts[i] == v then fi = i
		)
				
	if fi == 0 then 
		(
		print "Could not find texture vert index"
		return FALSE
		)w
		
	mapVertIndex = (meshop.getMapFace obj 1 f)[fi]		
	
	debug = false
	if (debug) then
		(
		if v != mapvertindex then print "vert changed"
		)
	
	return (mapVertIndex as integer)

	)

/******************************************************************************************/
/* WE HAVE TO MAKE OUR OWN NXM MATRIX MULTIPLICATION FUNCTIONS ... THANKS DISCREET 
/******************************************************************************************/

fn matrix3toBigMatrix m3 =
(
outB = bigMatrix 4 3
for i = 1 to 4 do
	(
	for j = 1 to 3 do
		(
		outB[i][j] = m3[i][j]
		)
	)
outB
)

fn mMultiply inA inB = 
(

if (inA.columns != inB.rows) then
	(
	print "Wrong Matrix Size...Aborted"
	return FALSE
	)

C = bigMatrix inA.rows inB.columns

for i = 1 to inA.rows do
	(
	for j = 1 to inB.columns do
		(
		sum = 0
		for k = 1 to inA.columns do
			(
			sum += inA[i][k]*inB[k][j]
			)
		C[i][j] = sum	
		)	
	)
C
)

/*********************************************************************/
/* Initialize all our variables and vertex data channels
/*********************************************************************/

(
node = $ --operate on the current selection (MUST BE EDITABLE_MESH)
vCount = meshop.getNumVerts node
meshop.setNumMaps node 1
meshop.setNumMaps node (1 + 6)
meshop.freeMapChannel node 5
meshop.setMapSupport node 5 TRUE
meshop.setNumMapVerts node 5 vCount		
meshop.freeMapChannel node 6
meshop.setMapSupport node 6 TRUE
meshop.setNumMapVerts node 6 vCount	


global tanuArray = #()
global tanvArray = #()
global normalArray = #()
global scaleArray = #()
global curvatureArray = #()
smallDetCount = 0

for v = 1 to vCount do
	(
	append tanuArray [0,0,0]
	append tanvArray [0,0,0]
	append normalArray (getNormal node v)
	append scaleArray [0,0,0]
	)

/*********************************************************************/
/* Each triangle of the mesh, get the 3 tri verts and
/* accumulate the tangent space vectors in the above arrays
/*********************************************************************/

faceCount = meshop.getNumFaces node

for f = 1 to faceCount do
	(
	triverts = (meshop.getVertsUsingFace node f) as array
	
	if triverts.count != 3 then 
		(
		print ("warning. Tricount of " + triverts.count as string)
		continue
		)
	
	index0 = triverts[1]
	index1 = triverts[2]
	index2 = triverts[3]
		
	--work some magic to find out which texture uvs correspond to these verts	
	uvindex0 = getTVfromV node f index0
	uvindex1 = getTVfromV node f index1
	uvindex2 = getTVfromV node f index2	
	
	v0 = meshop.getVert node index0 
	v1 = meshop.getVert node index1
	v2 = meshop.getVert node index2 
	
	v0 = v0* (inverse node.transform)
	v1 = v1* (inverse node.transform)
	v2 = v2* (inverse node.transform)
	
	v1 = v1-v0
	v2 = v2-v0
	 	
	uv0 = meshop.getMapVert node 1 uvindex0 
	uv1 = meshop.getMapVert node 1 uvindex1 
	uv2 = meshop.getMapVert node 1 uvindex2 
	
	uv0[3] = 0.0 	
	uv1[3] = 0.0 
	uv2[3] = 0.0 
		
	uv1 = uv1 - uv0
	uv2 = uv2 - uv0
	
	det = uv1.x*uv2.y-uv2.x*uv1.y

	if (abs det) > 0.000001 then 
	(
	u=0
	v=0
	u = u-uv0.x
	v = v-uv0.y
	p0 = v0 + v1*((u*uv2.y-uv2.x*v)/det)+v2*((uv1.x*v-u*uv1.y)/det);
	
	u=1
	v=0
	u = u-uv0.x
	v = v-uv0.y
	p1 = v0 + v1*((u*uv2.y-uv2.x*v)/det)+v2*((uv1.x*v-u*uv1.y)/det);
	
	u=0
	v=1
	u = u-uv0.x
	v = v-uv0.y
	p2 = v0 + v1*((u*uv2.y-uv2.x*v)/det)+v2*((uv1.x*v-u*uv1.y)/det);
	
	d1 = p1-p0
	d2 = p2-p0	
	l1 = length d1
	l2 = length d2		
	d1*=1.0/l1
	d2*=1.0/l2
		
	scaleArray[index0] += [l1,l2,1.0]
	scaleArray[index1] += [l1,l2,1.0]
	scaleArray[index2] += [l1,l2,1.0]
		
	tanuArray[index0] += d1
	tanvArray[index0] += d2
	tanuArray[index1] += d1
	tanvArray[index1] += d2
	tanuArray[index2] += d1
	tanvArray[index2] += d2
	)
	else (smallDetCount += 1)
)

print "SMALL DETERMINANT COUNT:"
print smallDetCount 

for v = 1 to vCount do
	(

	--normalArray[v] = cross tanuArray[v] tanvArray[v]
	
	tanuArray[v] = normalize tanuArray[v]
	l = length tanvArray[v]
	if l < .01 then
		(
		tanvArray[v] = cross tanuArray[v] normalArray[v]
		)
	tanvArray[v] = normalize tanvArray[v]	
	
	if scaleArray[v][3]==0 then print ("0 scale array at vert " + v as string)

	scaleArray[v] = scaleArray[v]/scaleArray[v][3]
	

	) 

	for vk = 1 to vCount do
	(
			
		--find all the neighbor verts
		sharedFaces = (meshop.getFacesUsingVert node vk) as array
		sharedVerts = (meshop.getVertsUsingFace node sharedFaces) as array
		
		global A = BigMatrix (sharedVerts.count-1) 2	
		global B = BigMatrix (sharedVerts.count-1) 1
		n  = 0
		for v in sharedVerts do
		(
			if v != vk then
			(
				n += 1
				vkPos = meshop.getVert node vk
				vPos = meshop.getVert node v
				vkPos = vkPos * (inverse node.transform)
				vPos = vPos * (inverse node.transform)
				
				vOffset = [0,0,0]
				t = vPos-vkPos
				
				vOffset.x = dot tanuArray[vk] t
				vOffset.y = dot tanvArray[vk] t
				vOffset.z = dot normalArray[vk] t
				
				for i = 1 to 2 do (vOffset[i]  = vOffset[i]*vOffset[i]) --square the x and y components only
				
				A[n][1] = vOffset[1]
				A[n][2] = vOffset[2]
				B[n][1] = -vOffset[3]
				
				
			)
		)
			
		--matrix least squares solution
		a1 = mmultiply (transpose A) A
		inv = invert a1
		a2 = mmultiply inv (transpose A)
		final = mmultiply a2 B
				
		append curvatureArray final
	
	
		--if final is 0, then make it very small instead.
		--this is for planar surfaces
		if (abs final[1][1] < .00001) then final[1][1] = .0000001
		if (abs final[2][1] < .00001) then final[2][1] = .0000001
		
		--print (([final[1][1], final[2][1], 0.0]) as string)
		
		/*********************************************************************/
		/* After all that hard work, let's color some verts
		/*********************************************************************/
		
		meshop.setMapVert node 5 vk ([final[1][1], final[2][1], 0.0])
		meshop.setMapVert node 6 vk ([scaleArray[vk][1],scaleArray[vk][2], 0.0])	
		
	)
		
)