Texture coordinate formula generation

DEFSNG A-Z

CONST last = 2
WIDTH 80,50:CLS ' Clear screen

TYPE vec3
  x AS SINGLE
  y AS SINGLE
  z AS SINGLE
END TYPE
TYPE vec2
  x AS SINGLE
  y AS SINGLE
END TYPE

' Actual sample data from face 35 in leaf 1039 in sp_a1_intro4.bsp
' Result we are looking for:
'   u = x*-1.95841 + y*7.48778 + z*-2.02429 + -3826.37
'   v = x*-7.73966 + y*-2.02429 + z*0 + -11783.1
'     X,       Y,       Z,     U,   V
'DATA -1808,    142.735, 128,   524.0944383000005,    1921.268246849998
'DATA -1791.87, 144.751, 128,   507.6006494799994,    1792.346562409999
'DATA -1775.75, 83.2621, 128,   15.61572463800075,    1792.054608590999
'DATA -1790.86, 77.2148, 128,   -0.07355225599985715, 1921.242360107999

'# U formula: u = x*-7.43992 + y*-2.09195 + z*-14.0097 + 14611
'# V formula: u = x*-13.4867 + y*-3.79221 + z*7.72843 + 27058.5
'     X,       Y,       Z,     U,   V
'DATA 2126.87, -400    , -16  , -151 ,-232
'DATA 2094.61, -409.067, -16  ,  107 , 236
'DATA 2094.61, -409.07 , -176 , 2348 ,-999
'DATA 2126.87, -400    ,-176  , 2089 ,-1469

'     X,       Y,       Z,     U,   V
'DATA 1802227, -1980499, 3734283,   -899,249
'DATA 2621440, -1966080, 4587520 , -1530,923
'DATA -2621440, -1966080, 4587520 , 2564,872
'DATA -1802190, -1980478, 3734276 , 1916,214

DATA  4587520 ,   -368640,    -819200 , -352, 551
DATA  4915200  ,  -655360 ,   -778240 , -364, 1890
DATA  4953447 ,   -239506,    -794673 ,  689, 1167

' Read the coordinate data into p() and u().
DIM p(last) AS vec3, uv(last) AS vec2
FOR n = 0 TO last:  READ p(n).x, p(n).y, p(n).z, uv(n).x, uv(n).y: NEXT

' Goal: Determine a,b,c,d such that
'        x*a + y*b + z*c + d  = u

DIM pa AS vec3, pb AS vec3, pc AS vec3
DIM SHARED re0 AS vec3, re1 AS vec3, re2 AS vec3

' Convert these three 3D points into 2D points (pa,pb,pc).
' We also get a matrix in re0,re1,re2.
Flatten p(0),p(1),p(2), pa,pb,pc

' Invert the matrix into qe0,qe1,qe2.
DIM qe0 AS vec3, qe1 AS vec3, qe2 AS vec3
matinv3 re0,re1,re2, qe0,qe1,qe2

' Debug output: Print both matrices.
PRINT "re0:";re0.x;re0.y;re0.z
PRINT "re1:";re1.x;re1.y;re1.z
PRINT "re2:";re2.x;re2.y;re2.z

PRINT "qe0:";qe0.x;qe0.y;qe0.z
PRINT "qe1:";qe1.x;qe1.y;qe1.z
PRINT "qe2:";qe2.x;qe2.y;qe2.z

' Debug output: Print the 2D coordinates.
PRINT "---------"
PRINT "pa:";pa.x;pa.y;"(";pa.z;")":PRINT "pb:";pb.x;pb.y;"(";pb.z;")": PRINT "pc:";pc.x;pc.y;"(";pc.z;")":PRINT "---"

DIM tmp AS vec3
PRINT "***** FOR U: *******"

' Now that we have 2D points,
' Solve a,b,c in the following equation:  x*a + y*b + c = u
' Using the three 2D points.
Solve3vars pa.x,pa.y, pb.x,pb.y, pc.x,pc.y, uv(0).x,uv(1).x,uv(2).x, a,b,c

' Debug output: Test whether the conversion worked.
PRINT "abc:";a;b;c
res = pa.x*a + pa.y*b + c: PRINT "u0: "; res; "expect:";uv(0).x;"error:";ABS(uv(0).x-res)
res = pb.x*a + pb.y*b + c: PRINT "u1: "; res; "expect:";uv(1).x;"error:";ABS(uv(1).x-res)
res = pc.x*a + pc.y*b + c: PRINT "u2: "; res; "expect:";uv(2).x;"error:";ABS(uv(2).x-res)
' And it did work.

' Now convert the 2D factors (a,b) into 3D by using the inverse matrix.
' The value of c generated above won't actually be needed for anything.
' We want x*a + y*b + z*c + d = u. These are a new set of a,b,c,d so we call them aa,bb,cc,dd.
'
tmp.x = a : tmp.y = b: tmp.z = 0 ' z can be anything.
aa = vdot3(tmp, qe0)
bb = vdot3(tmp, qe1)
cc = vdot3(tmp, qe2)
' We got a,b,c. Solve d by substituting real coordinate data into the equation.
' It doesn't matter which coordinate we choose for the substitution.
dd = uv(0).x - (p(0).x*aa + p(0).y*bb + p(0).z*cc)
PRINT "aabbccdd:";aa;bb;cc;dd
FOR n = 0 TO last
  res = p(n).x * aa + p(n).y * bb + p(n).z * cc + dd
  PRINT "Got:";res;" expect:";uv(n).x;"error:";ABS(uv(n).x-res)
NEXT

PRINT "***** FOR V: *******"

Solve3vars pa.x,pa.y, pb.x,pb.y, pc.x,pc.y, uv(0).y,uv(1).y,uv(2).y, a,b,c
PRINT "abc:";a;b;c
res = pa.x*a + pa.y*b + c: PRINT "u0: "; res; "expect:";uv(0).y;"error:";ABS(uv(0).y-res)
res = pb.x*a + pb.y*b + c: PRINT "u1: "; res; "expect:";uv(1).y;"error:";ABS(uv(1).y-res)
res = pc.x*a + pc.y*b + c: PRINT "u2: "; res; "expect:";uv(2).y;"error:";ABS(uv(2).y-res)

tmp.x = a : tmp.y = b: tmp.z = 0 ' z can be anything.
aa = vdot3(tmp, qe0)
bb = vdot3(tmp, qe1)
cc = vdot3(tmp, qe2)
dd = uv(0).y - (p(0).x*aa + p(0).y*bb + p(0).z*cc)
PRINT "aabbccdd:";aa;bb;cc;dd
FOR n = 0 TO last
  res = p(n).x * aa + p(n).y * bb + p(n).z * cc + dd
  PRINT "Got:";res;" expect:";uv(n).y;"error:";ABS(uv(n).y-res)
NEXT

END



SUB vadd2(a AS vec2, b AS vec2, result AS vec2)
  result.x = a.x + b.x
  result.y = a.y + b.y
END SUB
SUB vsub2(a AS vec2, b AS vec2, result AS vec2)
  result.x = a.x - b.x
  result.y = a.y - b.y
END SUB
SUB vmul2(a AS vec2, b AS vec2, result AS vec2)
  result.x = a.x * b.x
  result.y = a.y * b.y
END SUB
SUB vmul2s(a AS vec2, b, result AS vec2)
  result.x = a.x * b
  result.y = a.y * b
END SUB
SUB vdiv2(a AS vec2, b AS vec2, result AS vec2)
  result.x = a.x / b.x
  result.y = a.y / b.y
END SUB
FUNCTION vdot2(a AS vec2, b AS vec2)
  vdot2 = a.x * b.x + a.y * b.y
END SUB
FUNCTION vlen2(a AS vec2)
  vlen2 = SQR(vdot2(a,a))
END SUB
SUB vnorm2(a AS vec2, result AS vec2)
  l = vlen2(a)
  result.x = a.x / l
  result.y = a.y / l
END SUB
FUNCTION vxs2(a AS vec2, b AS vec2)
  vxs2 = a.x*b.y - a.y*b.x
END FUNCTION

SUB vadd3(a AS vec3, b AS vec3, result AS vec3)
  result.x = a.x + b.x
  result.y = a.y + b.y
  result.z = a.z + b.z
END SUB
SUB vsub3(a AS vec3, b AS vec3, result AS vec3)
  result.x = a.x - b.x
  result.y = a.y - b.y
  result.z = a.z - b.z
END SUB
SUB vmul3(a AS vec3, b AS vec3, result AS vec3)
  result.x = a.x * b.x
  result.y = a.y * b.y
  result.z = a.z * b.z
END SUB
SUB vmul3s(a AS vec3, b, result AS vec3)
  result.x = a.x * b
  result.y = a.y * b
  result.z = a.z * b
END SUB
SUB vdiv3(a AS vec3, b AS vec3, result AS vec3)
  result.x = a.x / b.x
  result.y = a.y / b.y
  result.z = a.z / b.z
END SUB
FUNCTION vdot3(a AS vec3, b AS vec3)
  vdot3 = a.x * b.x + a.y * b.y + a.z * b.z
END SUB
FUNCTION vlen3(a AS vec3)
  vlen3 = SQR(vdot3(a,a))
END SUB
SUB vnorm3(a AS vec3, result AS vec3)
  l = vlen3(a)
  result.x = a.x / l
  result.y = a.y / l
  result.z = a.z / l
END SUB
SUB vxs3(a AS vec3, b AS vec3, result AS vec3)
  DIM temp AS vec2
  temp.x = a.y*b.z - a.z*b.y
  temp.y = a.z*b.x - a.x*b.z
  result.z = a.x*b.y - a.y*b.x
  result.x = temp.x
  result.y = temp.y
END SUB
SUB matmul3(m0 AS vec3,m1 AS vec3,m2 AS vec3, n0 AS vec3,n1 AS vec3,n2 AS vec3, r0 AS vec3,r1 AS vec3,r2 AS vec3)
  r0.x = m0.z*n2.x + m0.y*n1.x + m0.x*n0.x
  r0.y = m0.z*n2.y + m0.y*n1.y + m0.x*n0.y
  r0.z = m0.z*n2.z + m0.y*n1.z + m0.x*n0.z
  r1.x = m1.z*n2.x + m1.y*n1.x + m1.x*n0.x
  r1.y = m1.z*n2.y + m1.y*n1.y + m1.x*n0.y
  r1.z = m1.z*n2.z + m1.y*n1.z + m1.x*n0.z
  r2.x = m2.z*n2.x + m2.y*n1.x + m2.x*n0.x
  r2.y = m2.z*n2.y + m2.y*n1.y + m2.x*n0.y
  r2.z = m2.z*n2.z + m2.y*n1.z + m2.x*n0.z
END SUB
SUB matinv3(m0 AS vec3,m1 AS vec3,m2 AS vec3, r0 AS vec3,r1 AS vec3,r2 AS vec3)
  DIM tmp AS vec3
  vxs3 m1,m2, tmp: invdet = 1! / vdot3(m0, tmp)
  PRINT "matinv3: det=";det
  vxs3 m1,m2, tmp: r0.x = tmp.x * invdet: r1.x = tmp.y * invdet: r2.x = tmp.z * invdet
  vxs3 m2,m0, tmp: r0.y = tmp.x * invdet: r1.y = tmp.y * invdet: r2.y = tmp.z * invdet
  vxs3 m0,m1, tmp: r0.z = tmp.x * invdet: r1.z = tmp.y * invdet: r2.z = tmp.z * invdet
END SUB

SUB Flatten(a AS vec3, b AS vec3, c AS vec3, out1 AS vec3,out2 AS vec3,out3 AS vec3)
  DIM ba AS vec3, ca AS vec3, normal AS vec3
  DIM vX1 AS vec3, vX2 AS vec3
  DIM vZ1 AS vec3, vZ2 AS vec3, vY AS vec3
  DIM m0 AS vec3, m1 AS vec3, m2 AS vec3
  vsub3  b,a, ba
  vsub3  c,a, ca
  vxs3   ba, ca, normal
  vnorm3 normal, vX2
  vX1.x = 0: vX1.y = 0: vX1.z = 1
  vxs3 vX1, vX2, vY
  'PRINT "vX2 (normal):"; vX2.x;vX2.y;vX2.z
  'PRINT "vY:"; vY.x;vY.y;vY.z
  IF vdot3(vY,vY) > 0 THEN
    vnorm3 vY, vY
    vxs3 vX1, vY, vZ1 : vnorm3 vZ1, vZ1
    vxs3 vX2, vY, vZ2 : vnorm3 vZ2, vZ2
    'PRINT "vZ1:"; vZ1.x;vZ1.y;vZ1.z
    'PRINT "vZ2:"; vZ2.x;vZ2.y;vZ2.z
    m0.x = vX1.x: m0.y = vY.x: m0.z = vZ1.x
    m1.x = vX1.y: m1.y = vY.y: m1.z = vZ1.y
    m2.x = vX1.z: m2.y = vY.z: m2.z = vZ1.z
    matmul3 m0,m1,m2, vX2,vY,vZ2, re0,re1,re2
  ELSE
    re0.x = 1 : re0.y = 0: re0.z = 0
    re1.x = 0 : re1.y = 1: re1.z = 0
    re2.x = 0 : re2.y = 0: re2.z = 1
  END IF
  out1.x = vdot3(a,re0)
  out1.y = vdot3(a,re1)
  out1.z = vdot3(a,re2)
  out2.x = vdot3(b,re0)
  out2.y = vdot3(b,re1)
  out2.z = vdot3(b,re2)
  out3.x = vdot3(c,re0)
  out3.y = vdot3(c,re1)
  out3.z = vdot3(c,re2)
END SUB

SUB Solve3vars(x1,y1, x2,y2, x3,y3, u1,u2,u3, a,b,c)
  det = (y1*(x3-x2) + y2*(x1-x3) + y3*(x2-x1))
  invdet = 1! / det
  PRINT "Solve3: det=";det
  a   = (y1*(u3-u2) + y2*(u1-u3) + y3*(u2-u1)) * invdet
  b   = (x1*(u3-u2) + x2*(u1-u3) + x3*(u2-u1)) * -invdet
  c   = (y1*(u2*x3-u3*x2) + y2*(u3*x1-u1*x3) + y3*(u1*x2-u2*x1)) * invdet
END SUB

' Solve4vars is unused.
SUB Solve4vars(x1,y1,z1, x2,y2,z2, x3,y3,z3, x4,y4,z4, u1,u2,u3,u4, a,b,c,d)
  ' ratsimp(solve( [x1*a+y1*b+z1*c+d=u1,x2*a+y2*b+z2*c+d=u2,x3*a+y3*b+z3*c+d=u3,x4*a+y4*b+z4*c+d=u4], [a,b,c,d] ))
  det =(z4*((x2-x1)*y3+(x1-x3)*y2+(x3-x2)*y1) _
       +z3*((x1-x2)*y4+(x2-x4)*y1+(x4-x1)*y2) _
       +z2*((x4-x3)*y1+(x3-x1)*y4+(x1-x4)*y3) _
       +z1*((x3-x4)*y2+(x4-x2)*y3+(x2-x3)*y4))
  a = (z4*(y3*(u2-u1)+y2*(u1-u3)+y1*(u3-u2)) _
      +z3*(y4*(u1-u2)+y2*(u4-u1)+y1*(u2-u4)) _
      +z2*(y4*(u3-u1)+y3*(u1-u4)+y1*(u4-u3)) _
      +z1*(y4*(u2-u3)+y3*(u4-u2)+y2*(u3-u4))) / det
  b = (z4*(x3*(u2-u1)+x2*(u1-u3)+x1*(u3-u2)) _
      +z3*(x4*(u1-u2)+x2*(u4-u1)+x1*(u2-u4)) _
      +z2*(x4*(u3-u1)+x3*(u1-u4)+x1*(u4-u3)) _
      +z1*(x4*(u2-u3)+x3*(u4-u2)+x2*(u3-u4))) / -det
  c = (y4*(x3*(u2-u1)+x2*(u1-u3)+x1*(u3-u2)) _
      +y3*(x4*(u1-u2)+x2*(u4-u1)+x1*(u2-u4)) _
      +y2*(x4*(u3-u1)+x3*(u1-u4)+x1*(u4-u3)) _
      +y1*(x4*(u2-u3)+x3*(u4-u2)+x2*(u3-u4))) / det
  d = (z4*(y3*(u1*x2-u2*x1)+y2*(u3*x1-u1*x3)+y1*(u2*x3-u3*x2)) _
      +z3*(y4*(u2*x1-u1*x2)+y2*(u1*x4-u4*x1)+y1*(u4*x2-u2*x4)) _
      +z2*(y4*(u1*x3-u3*x1)+y3*(u4*x1-u1*x4)+y1*(u3*x4-u4*x3)) _
      +z1*(y4*(u3*x2-u2*x3)+y3*(u2*x4-u4*x2)+y2*(u4*x3-u3*x4))) / det
END SUB