This answer only covers Affine Texture Mapping. Consider this a partial solution.
I'll be assuming you can process your obj to a list of triangles. Each triangle having three vertex. And each vertex having their 3D space coordinates (x, y, z) and their texture coordinates (u, v).
For each triangle, sort the points by one texture coordinate - let us say -u. Once sorted, let us call these points A, B, and C.
You are going to define - at most - two ranges: The range from A.u to B.u, and the range from B.u to C.u. And I say at most, because perhaps some of these coordinates are equal. If they are, you don't need that range.
The start and end of each range define two lines, which constraint the texture space on the u axis. On The v axis, you will constraint the texture space by two line equations of the form v = mu + n.
For the first range, the lines equations correspond to the lines l_AB and l_AC. For the second range, the lines are l_AC, l_BC.
You can figure out the line equations with a little algebra. The slope is the good old formula:
m = (B.v - A.v)/(B.u - A.u)
And for the intercept, we can replace one of the points:
A.v = m(A.u) + n
=>
n = A.v - m(A.u)
Remember that we are doing the lines l_AB and l_AC for the first range, which cross at the point A. And the lines l_AC and l_BC for the second, which cross at the point C. Thus, we can use the v coordinates of the other points to determine which line has lower values in each range.
Now, you can iterate over the texture. For example, for the first range:
if (B.v < C.v)
{
l_1 = l_AB
l_2 = l_AC
}
else
{
l_1 = l_AC
l_2 = l_AB
}
for (int u = (int)floor(A.u); u <= (int)ceil(B.u); u++)
{
for (int v = (int)floor(l_1(u)); v <= (int)ceil(l_2(u)); v++)
{
// Something with u and v
}
}
I don't know how you want to handle texels that are partially on a triangle. However, I'm suggesting you include them. Which is the reason behind the rounding (floor and ceil) above.
We need to convert these coordinates to 3D, of course. For that we will have to express these u and v coordinates as a linear combination of the lines of the triangle.
We can use a single conversion for the whole triangle - being the points of a triangle always coplanar. Let us base it on the point A. So we will work assuming A = 0.
We will define vectors in texture space: dB = B - A, dC = C - A and d = uv - A. So we could say:
p * dB + q * dC = d
Which is actually two equations:
p * dB.u + q * dC.u = d.u
p * dB.v + q * dC.v = d.v
Since we have two linear equations and two unknown (p, q) we can solve the system. Let us do elimination, Gauss-Jordan style:
p q
┌ ┐
│ dB.u dC.u │ d.u │
│ │ │
│ dB.v dC.v │ d.v │
└ ┘
Divide the first equation by dB.u =>
p q
┌ ┐
│ dB.u dC.u │ d.u │
│ ------ ------ │ ----- │
│ dB.u dB.u │ dB.u │
│ │ │
│ dB.v dC.v │ d.v │
└ ┘
=>
p q
┌ ┐
│ dC.u │ d.u │
│ 1 ------ │ ------ │
│ dB.u │ dB.u │
│ │ │
│ dB.v dC.v │ d.v │
└ ┘
Subtract from the second equation the first equation dB.v times =>
p q
┌ ┐
│ dC.u │ d.u │
│ 1 ------ │ ------ │
│ dB.u │ dB.u │
│ │ │
│ dC.u │ d.u │
│ dB.v -(dB.v)1 dC.v -(dB.v)------ │ d.v -(dB.v)------ │
│ dB.u │ dB.u │
└ ┘
=>
p q
┌ ┐
│ dC.u │ d.u │
│ 1 ------ │ ------ │
│ dB.u │ dB.u │
│ │ │
│ dB.v * dC.u │ dB.v * d.u │
│ 0 dC.v - ------------- │ d.v - ------------ │
│ dB.u │ dB.u │
└ ┘
Divide the second equation by dC.v - (dB.v * dC.u)/dB.u =>
p q
┌ ┐
│ dC.u │ d.u │
│ 1 ------ │ ------ │
│ dB.u │ dB.u │
│ │ │
│ dB.v * dC.u │ dB.v * d.u │
│ dC.v - ------------- │ d.v - ------------ │
│ dB.u │ dB.u │
│ 0 ---------------------- │ -------------------- │
│ dB.v * dC.u │ dB.v * dC.u │
│ dC.v - ------------- │ dC.v - ------------- │
│ dB.u │ dB.u │
└ ┘
=>
p q
┌ ┐
│ dC.u │ d.u │
│ 1 ------ │ ------ │
│ dB.u │ dB.u │
│ │ │
│ │ dB.v * d.u │
│ │ d.v - ------------ │
│ │ dB.u │
│ 0 1 │ -------------------- │
│ │ dB.v * dC.u │
│ │ dC.v - ------------- │
│ │ dB.u │
└ ┘
=>
p q
┌ ┐
│ dC.u │ d.u │
│ 1 ------ │ ------ │
│ dB.u │ dB.u │
│ │ │
│ │ dB.u * d.v dB.v * d.u │
│ │ ------------ - ------------ │
│ │ dB.u dB.u │
│ 0 1 │ ------------------------------- │
│ │ dB.u * dC.v dB.v * dC.u │
│ │ ------------- - ------------- │
│ │ dB.u dB.u │
└ ┘
=>
p q
┌ ┐
│ dC.u │ d.u │
│ 1 ------ │ ------ │
│ dB.u │ dB.u │
│ │ │
│ │ dB.u * d.v - dB.v * d.u │
│ │ ------------------------- │
│ │ dB.u │
│ 0 1 │ ----------------------------- │
│ │ dB.u * dC.v - dB.v * dC.u │
│ │ --------------------------- │
│ │ dB.u │
└ ┘
=>
p q
┌ ┐
│ dC.u │ d.u │
│ 1 ------ │ ------ │
│ dB.u │ dB.u │
│ │ │
│ │ dB.u * d.v - dB.v * d.u │
│ 0 1 │ --------------------------- │
│ │ dB.u * dC.v - dB.v * dC.u │
└ ┘
Subtract from the first equation the second equation dC.u/dB.u times =>
p q
┌ ┐
│ dC.u dC.u │ d.u dB.u * d.v - dB.v * d.u dC.u │
│ 1 ------ - ------ │ ------ - --------------------------- * ------ │
│ dB.u dB.u │ dB.u dB.u * dC.v - dB.v * dC.u dB.u │
│ │ │
│ │ dB.u * d.v - dB.v * d.u │
│ 0 1 │ --------------------------- │
│ │ dB.u * dC.v - dB.v * dC.u │
└ ┘
=>
p q
┌ ┐
│ │ d.u dB.u * d.v - dB.v * d.u dC.u │
│ 1 0 │ ------ - --------------------------- * ------ │
│ │ dB.u dB.u * dC.v - dB.v * dC.u dB.u │
│ │ │
│ │ dB.u * d.v - dB.v * d.u │
│ 0 1 │ --------------------------- │
│ │ dB.u * dC.v - dB.v * dC.u │
└ ┘
=>
p q
┌ ┐
│ │ dC.v * d.u - dC.u * d.v │
│ 1 0 │ --------------------------- │
│ │ dB.u * dC.v - dB.v * dC.u │
│ │ │
│ │ dB.u * d.v - dB.v * d.u │
│ 0 1 │ --------------------------- │
│ │ dB.u * dC.v - dB.v * dC.u │
└ ┘
This last simplification was computer assisted. Thanks Wolfram|Alpha.
So we have:
k = dB.u * dC.v - dB.v * dC.u
p = (dC.v * d.u - dC.u * d.v) / k
q = (dB.u * d.v - dB.v * d.u) / k
We can incorporate that to our code, and use those components to find the 3D coordinates.
VectorUV dB = B.uv - A.uv
VectorUV dC = C.uv - A.uv
if (B.v < C.v)
{
l_1 = l_AB
l_2 = l_AC
}
else
{
l_1 = l_AC
l_2 = l_AB
}
k = dB.u * dC.v - dB.v * dC.u
for (int u = (int)floor(A.u); u <= (int)ceil(B.u); u++)
{
for (int v = (int)floor(l_1(u)); v <= (int)ceil(l_2(u)); v++)
{
VectorUV d = new VectorUV(u - A.u, v - A.v)
p = (dC.v * d.u - dC.u * d.v) / k
q = (dB.u * d.v - dB.v * d.u) / k
coords = A.xyz + (B.xyz - A.xyz) * p + (C.xyz - A.xyz) * q
// …
}
}
I remind you that the code above only handles the first range of a triangle. For the other range you need another for loop. And, of course, you would be doing that for all triangles.
You would have an array where you can store those coordinates. Make it an array of sets, because you may find that the same texel is mapped to one position in 3D, multiple, or none. And you might find a texel has a position multiple times, in particular on the edges.
By the way. We have a division by zero if k = dB.u * dC.v - dB.v * dC.u = 0. I believe this happens if the shape in texture space does not make a triangle. I suggest you consider these cases as pathological.
