Gradients
0 - Intro
This little text presents some useful formulae about texture mapping. "Texture mapping" is to be taken in its general sense here: it actually concerns anything you want to be linearly interpolated on a 3D-polygon and projected on a 2D-screen. There's some maths involved, but it's worth it. It provides a unified vision of screen gradients as mere vectors, subject to special transform rules when 3D space is transformed.
1 - The problem
Suppose you have a quantity f defined at the n vertex of a n-gon, and linearly interpolated on every other points of this polygon. This quantity ('field') can be anything you want, provided it remains linear in space: UV mapping coordinates, Gouraud illumination, normal, transparency, etc...
I've already talked about projecting linear fields on screen in the article About Texture Mapping, so I will only recall the conclusion: if f(x,y,z) is linear in space, then w.f(x.w,y.w,w) will be linear on screen (where w is the quantity 1/z). Therefore the gradients of the function w.f() will be constants all over the polygon. Up to us, now, to compute them.
2 - Basic formulae
Let's take a triangle, even if the following is also valid for a generic n-gon:
P1 (=X1,Y1,F1)
+
/ \
u / \
/ \
/ \
/ \ v
P0 + \
(=X0,Y0,F0) ---- \
------ \
------- \
-----+ P2 (=X2,Y2,F2)
According to the definition, the gradient along X will be the variation of the quantity F (=w.f(x,y,z)) when we jump from point X to point X+1. Same goes for Y when going from point Y to Y+1.
Because these gradients are constants, let's compute them with this "raw" method: we start from points P0 and P1, where the values of F are respectively F0 and F1, and jump till point P2. We then have:
F2 = F0 + (X2-X0).DxF + (Y2-Y0).DyF
F2 = F1 + (X2-X1).DxF + (Y2-Y1).DyF
where DxF and DyF are our unknown gradients along X and Y.
Using the notations dXij = Xi-Xj, dYij = Yi-Yj, dFij=Fi-Fj (where i,j,k take the values 0,1 or 2), we can re-arrange this in:
dX20.DxF + dY20.DyF = dF20
dX21.DxF + dY21.DyF = dF21
(or even, in a general form: dXij.DxF + dYij.DyF = dFij)
We thus have 2 equations with 2 unknown variables (DxF and DyF). Inverting the whole gives the desired result:
dF20.dY12 - dF12.dY20
DxF = -----------------------
dX20.dY12 - dX12.dY20
(1)
dF12.dX20 - dF20.dX12
DyF = -----------------------
dY12.dX20 - dY20.dX12
3 - Some remarks
3-1. You will have noticed the great symmetry in this result. You can generate a lot of equivalent expressions by using the fact that: dXik = dXij + dXjk (or the same with Y or F).
You can even write something like:
dYij.dFjk - dYjk.dFij = dF/dX . ( dYij.dXjk - dYjk.dXij )
dXij.dFjk - dXjk.dFij = dF/dY . ( dYij.dXjk - dYjk.dXij )
where (i,j,k) is any permutation of indices (0,1,2).
3-2. The denominator is the same for DxF and DyF. You then only have to compute it once for all F you want to interpolate along both X and Y. Moreover, it has the useful property of being proportional the area (in pixel, oriented) of the projected polygon. Indeed, the oriented area of a parallelogram defined by 2 vectors u and v is equal to the half of the norm of the cross-product u^v. And by computing the norm (only one component is non-zero, actually) of the cross-product [dX10, dY10, 0] ^ [dX21, dY21, 0], we get the denominator, up to minus sign. Then we'll use the notation
A = dX20.dY12 - dX12.dY20
just to recall it's the area we're talking about.
Note that if you're consistent with polygons' orientation during modeling, you can do the backface cull by eliminating polygons whose projected area is negative (provided they don't cross the z=0 plane...)
3-3. These expressions for gradients are rather "resistant" to numerical round-off. If your triangle is very small, or has axis-aligned edges, or is very thin, etc., you won't have any particular problems using them. You just have to first check that A is non-zero. If it's the case, it just means that your three points P0, P1, P2 are aligned. You needn't draw the triangle, then...
3-4. Computing A first can be of great help to know better what sort of polygon we're about to draw. First it gives a crude approximation of the surface and allow one to adapt the shading we'll need to use. Or even the mipmap level. For instance, let's write the equivalent area expression in U-V mapping coordinates space (texels space). It's:
At = dU20.dV12 - dU12.dV12
This quantity is a constant if mapping coordinate associated with the polygone remains fixed (not for env-mapping, then) and only needs to be computed once for all at the scene's setup. Then, the ratio R=At/A will give you the average number of texel being mapped for 1 pixel. If R=1, you use 1 texel for 1 pixel, and your texture is well used. If R<1, you don't have enough texels to fit 1 pixel, and bilinear filtering is then required. The other way round, if R>1, your texture is too big for the polygon, and you can switch to mip-mapping. Which mipmap level? Well, each time you go down by one mipmap level, the At quantity is divided by 4. Then, you just have to go down amongst mipmap levels until the ration is around the value 1. Note that it's rather the quantity sqrt(At/A) that should be used here, since an area is quadratic quantity. But it's not really worth the sqrt()...
Eventually, you can try zo extend the above to the anisotropic cases, when the polygone is terribly stretched along one of the axis... Left as an exercise...
3-5. There's another way of deriving the formulae, somehow more obfuscated. Let's consider the 3D-space made of points Qi whose coordinates are (Xi,Yi,Fi). In such a space, our flat triangle on screen is a _plane_. A plane whose 'hyperplane equation' is just
Nx.x + Ny.y + Nz.z + h = 0
where (Nx,Ny,Nz) is the 'normal' to this plane and h is given by: h = -N.Q, provided the point Q belongs to the plane (it can for instance be one of our 'extended' vertex (Xi,Yi,Fi)).
Well, let's compute this normal. By taking the vectorial product Q0Q1 ^ Q1Q2 for instance. Explicitly, it gives, after taking the point Q0 as origin:
| dX10 | dX21 | dY10.dF21 - dY21.dF10
N = Q0Q1^Q1Q2 = | dY10 ^ | dY21 = | dX21.dF10 - dX10.dF21
| dF10 | dF21 | dX10.dY21 - dX21.dY10 (= -A)
Taking Q0 as origin amounts to setting h=0. Then, checking that a point Q belongs to the hyperplane is just a matter of verifying that: N.(Q-Q0) = 0.
For drawing our polygon on screen, we're going to start from a known point (Q0 for instance) and change by a quantity (dX,dY,dF) while staying on this hyperplane. We know how much the change is in X and Y on screen, and we'd like to know what the required value for dF is in order to remain stuck on the hyperplane. But as long as we don't know the gradients, we can't compute dF.
Let's then differentiate the equation N.(Q-Q0)=0. We have N.dQ = 0, where dQ = (dX,dY,dF). But according to linearity (or if you prefer, the constancy of gradients), we can choose any dQ we want. Let's take dQ = (1,0,dF). The restriction N.dQ=0 will then give us the unknown dF required to have a proper mapping:
Nx.1 + Ny.0 + Nz.dF = 0
Hence: dF = -Nx/Nz. And as we've stepped from X to X+1, this dF is actually the DxF. Now, if you take the explicit values above, for Nx and Nz, you'll find back the gradients' formulae.
For DyF, it's just the same, choosing dQ = (0,1,dF).
4 - Back to the 3rd dimension
You now know all what's required to draw your mapped polygon on screen, *provided* you have the projected vertex for computing the gradients. That's what's bad about the method.
For one reason, mainly: clipping. Indeed, when you project your (x,y,z) point on screen with: X=w.x, Y=w.y (these are 'normalized' projection formulae, with w=1/z), we must insure that w is strictly positive. This is absolutely not an obligation. Actually, except for the pathologic case z=0, there's no real problem in dealing with negative values for w, but this is rather tricky. You can even do the 3D-clipping with the quantities (w.x, w.y, w) only, but that's not really satisfying. Especially when you have a very long polygon starting on the z>0 side and ending on the z<0 side.
Moreover, when you clip your projected polygon on the screen side, it can happen that one of its vertex is off-screen by a very small amount. And during the clipping, this pathological vertex will be split in two very near vertex that will make the gradients formulae become rather instable. For instance, imagine you have a vertex with a slightly negative X coordinate. It will then be clipped on the left of the screen, hence producing a zero value for one of the above dXij quantity, and very small value for the corresponding dYij. This case *can* make your evaluation of DxF or DyF rather poor.
We'll then endeavour using only x,y,z in the expression (1) instead of X and Y. Let's re-introduce x,y,z, and f in (1) in replacement for X,Y and F.
With X=w.x, Y=w.y and F=w.f, this gives:
A.DxF = dF20.dY12 - dF12.dY20
= (w2.f2-w0.f0).(w1.y1-w2.y2) - (w1.f1-w2.x2).(x2.f2-x0.y0)
= ...
= f0y1w0w1-f0y2w0w2 + f1y2w1w2-f1y0w0w1 + f2y0w0w2-f2y1w1w2
= w0.w1.w2.[ f0.(y1.z2-y2.z1) + f1.(z0.y2-y0.z2) + f2.(y0.z1-z0.y1) ]
You will have noticed the use of z instead of w in the sum, so we can now factorize by w0.w1.w2.
Letting fi=xi will give us the corresponding expression for A, and lead us to the desired new set of expression for gradients, after simplification by the common factor w0.w1.w2:
f0.(y1.z2-y2.z1) + f1.(y2.z0-y0.z2) + f2.(y0.z1-y1.z0)
DxF = --------------------------------------------------------
A'
f0.(x1.z2-x2.z1) + f1.(x2.z0-x0.z2) + f2.(x0.z1-x1.z0)
(2) DyF = --------------------------------------------------------
A'
with: A' = x0.(y1.z2-y2.z1) + x1.(y2.z0-y0.z2) + x2.(y0.z1-y1.z0)
or A' = y0.(x1.z2-x2.z1) + y1.(x2.z0-x0.z2) + x2.(x0.z1-x1.z0)
or even A' = z0.(x1.y2-x2.y1) + z1.(x2.y0-x0.y2) + z2.(x0.y1-x1.y0)
[I've noted the denominator A' instead of A, just to recall that there's a missing w0w1w2 factor between the two: A=w0w1w2.A']
These expressions have a lot of common factors, so we could try simplifying the notations. Let's define the following vectors: Q0 = P1^P2, Q1=P2^P0 and Q2=P0^P1, where ^ is the cross product, as usual, and P0=(x0,y0,z0), P1=(x1,y1,z1), P2=(x2,y2,z2) are our vertex.
Explicitly writting the components of the Qi allows us to compound the expressions (2) into the compact, vectorial, form:
| DxF
(3) | DyF = 1/A' . ( f0*Q0 + f1*Q1 + f2.Q2 )
| DzF
DzF is of no use, not even defined, for now, but will be needed afterwards. Let's give here the explicit expression for it:
A'.DzF = f0.(x1.y2-x2.y1) + f1.(x2.y0-x0.y2) + f2.(x0.y1-x1.y0)
Well, the above condensed vectorial notation is rather interesting and needs further investigations. Let's use the transform matrix Mo to the basis (Q0,Q1,Q2). It simply is the matrix made with the vectors Q0,Q1,Q2 put in columns: Mo = [Q0,Q1,Q2] = [P1^P2,P2^P0,P0^P1]. Then, the parenthesied expression in (3) is simply a matter of Matrix-vector product:
(4) A'.dF = Mo.f
with f being the vector (f0,f1,f2) and dF being our gradient vector: dF = (DxF, DyF, DzF).
5 - Algebraic intermission
At this point, if you're used to algebrae, you will most probably have caught the hidden goal of all these (obfuscated?) manipulations. This Mo matrix should have given you a hint about the whats-and-whys of all these apparent complications. Anyway, I will here recall some basic, but useful results, concerning matrix manipulation, just to make the notations clear. You can skip this paragraph if you want...
To put it simply, matrices are just a way of specifying *linear* transformations we perform of our 3D space of vectors. Writing P'=M.P is simply a matter of gathering all the required coefficients into a 3x3 array, whilst defining some rules about addition and multiplication of such arrays. These rules must match the geometrical properties of such space transformations. That's why matrix addition and product are defined as:
(A+B)[i][j] = A[i][j] + B[i][j]
and (A.B)[i][j] = sum_over_k( A[i][k] * B[k][j] )
where [i][j] means 'row i, column j'. As can be seen from above, the addition does commute, whereas multiplication doesn't: (A+B) = (B+A), but (A.B) != (B.A).
Transpose and inverse of a matrix M will now be noted tM and iM respectively. The transpose of (A)[i][j] is simply (A)[j][i] (symmetry about the diagonal). The inverse of A is the matrix verifying: A.iA = Id, where Id is the unit matrix having 1 on diagonal, 0 elsewhere.
They must obey the rules: t(A+B) = tA+tB, t(A.B) = tB.tA and i(A.B) = iB.iA in order to conserve group's properties. Note that i(A+B) can't be easily expressed. Hence, we can deduce from the definition A.iA=I, and from the fact that Id does commute with any matrix, the additional properties:
iA.A = A.iA = I -> iA and A are commuting
=> t(A.iA) = t(iA.A) = tI = I
=> t(iA).tA = tA.t(iA) = I -> the inverse of tA is just t(iA)
tA.i(tA) = I = tA.tiA -> but also i(tA). So: i(tA) = t(iA)
etc...
I know it seems trivial properties, but it actually isn't. One has to be careful about 'trivialities' when dealing with non-commuting groups. What follows will provide examples for this warning...
Transpose is useful for computing scalar products. When you write: U.V = Ux.Vx + Uy.Vy + Uz.Vz, it can be seen as a product of 3 x 1 'row' matrix with a 1 x 3 'column' matrix. Matrix are to be understood as general n x m array instead of just square n x n ones. Hence, if we identify 3d vectors with the 1 x n column matrix whose coefficients are just (x,y,z), one can write the scalar product as:
U.V = t(U).(V)
where (U) and (V) are the column matrix associated with the vectors U and V.
The power of this notation shows up when we transform the vectors with a n x n matrix M. This transformation can be anything linear (rotation, scaling, shear,...). One can ask how will the scalar product transform upon this transformation. Since U'=M.U and V'=M.U are our new transformed vectors, we then have:
U'.V' = t(U').(V') = t(M.U).(M.V) = tU.tM.M.V = tU.(tM.M).V
Hence, U'.V' has no good reason of being equal to U.V, except if tM.M = Id (the transformation M is then called an isometry).
Let's now have a look at the cross-product. When you *define* a vector as being: N=u^v, you can't assume that transforming both u and v with a matrix M will have the same effect on N. This is actually not the case, just because N is not an ordinary vector you can choose how ever you want, but must indeed obey the definition N=u^v. Hence you don't have, a priori, N' = M.N if you transform u and v according to u'=M.u and v'=M.v. The exact value for N' is just: N' = u'^v'. Anyway, because of linearity, you can assume that N' has something to do with the original N. Maybe not the simple relation N'=M.N, but something resembling. Let's then write: N' = A.N, where A is an unknown matrix. What is the relation between A and the transform M? Let's find it out by recalling that if N'=u'^v', we then must have the necessary conditions: N'.u' = 0 and N'.v' = 0. Using matrix notations, it amounts to having:
N'.u' = 0
= tN'.U' = t(A.N).(M.U) = tN.(tA.M).U = t(tM.A.N).U
and the same for v': t(tM.A.N).V = 0
Hence the vector (tM.A.N) is orthogonal to *both* u and v. It must then be proportional u^v, since u^v is the only vector having this property too (up to a scaling). So: tM.A.N = s.N, where 's' is a scaling factor. This must be true for any vector N and any transform M, because of unicity of A.
Hence: tM.A = s(M).Id, and conservation of norm and orientation will force s being the determinant of M, noted detM. We thus have: tM.A = detM.Id => A=detM.itM (unitary matrix), leading to:
( M.U ^ M.V ) = detM . itM.( U ^ V )
and certainly *not* M.( U ^ V ). How bad. :)
This is often a confusing point when rotating normals, for instance. Since rotations are isometries, we have iM = tM, and writing N'= M.N is a valid relation, but shouldn't be taken as granted whenever iM != tM (scaling, for instance). In addition, we often need normal up to a scaling, and consequently often discard the 'detM' factor. A good way of checking you've transformed vectors right in your engine is to try to apply negative scaling to your meshes. Results are often... hmm... visually surprising.
Now that we're off with scalar product and cross-product, let's finish this inspection with mixed product. It is defined from three vectors A, B and C by: m = A.(B^C), and represents the volume enclosed between this three vector. It has the symmetry property that: m = A.(B^C) = B.(C^A) = C.(A^B). How does it transform upon M? Just write it with matrix:
m' = A'.(B'^C') = t(MA).(MB^MC) = detM tA.tM.itM.(B^C) = detM tA.(B^C)
= detM . m
It's then 'almost' a constant upon linear transforms and shows that detM measures the space dilatation involved by the transform.
We are now ready to get back to the matrix Mo of the previous paragraph. It was defined from three vectors P0,P1,P2 as: Mo = [P1^P2,P2^P0,P0^P1]. Just for the fun of it, let's compute the product of Mo with the transpose of the matrix [P0,P1,P2]:
[ P0--- ] [ P1^P2, P2^P0, P0^P1 ]
t[P0,P1,P2].Mo = [ P1--- ] .[ | | | ]
[ P2--- ] [ | | | ]
[ P0.(P1^P2) P0.(P2^P0) P0.(P0^P1) ]
= [ P1.(P1^P2) P1.(P2^P0) P1.(P0^P1) ]
[ P2.(P1^P2) P2.(P2^P0) P2.(P0^P1) ]
Most of the mixed products are null, the remaining ones being equal to m = P0.(P1^P2). In fact:
[ m 0 0 ]
t[P0,P1,P2] . Mo = [ 0 m 0 ] = m.Id
[ 0 0 m ]
So, up to a mixed product m=P0.(P1^P2), Mo and t[P0,P1,P2] are just reciprocal matrix. Put otherwise, we have:
i[P0,P1,P2] = 1/m .t[ P1^P2, P2^P0, P0^P1 ]
This provides us a means of finding the inverse of a given [P0,P1,P2] matrix. And if we focus on the quantity m instead, then:
(5) m = Mo . t[P0,P1,P2]
In addition, one has: det[P0,P1,P2] = P0.(P1^P2) = m. Consequently, det(Mo) = m^2.
We are now almost finished with matrix. There's only one last point to take care of: We must adapt the above rules to the case of points instead of vector. Actually, when dealing with points, the general linear transform must account for translations and not just matrix. A transformed *point* will then be P'= M.(P + T), instead of just V' = M.V. We adopt here the convention of performing the translation T first, and then the vectorial transform M. It's just a matter of having slightly simpler expressions.
This additional translation vector T will add some terms in the transform rules. E.g.:
P' ^ Q' = M.(P+T)^M(Q+T) = detM itM.[ (P^Q) + T^(Q-P) ]
A'.(B'^C') = ... = detM A.(B^C) + detM T.( B^C + C^A + A^B )
You will note that this last vector B^C + C^A + A^B is just proportional to the normal of the triangle defined by the points (A,B,C). Just because:
N = AB^BC = (B-A)^(C-B) = B^C + 0 - A^C + A^B = B^C + C^A + A^B
Then, in compact form, we have: m' = detM.( m + T.N ) for transforming the mixed product.
6 - Back to gradients
Now that we're done with the basics, let's get back to expression (4) [A'.dF = Mo.f] and investigate it further. If you look at the exact expression of A', you can recognize a mixed product in it. It's just: A' = P0.(P1^P2). So (4) only contains quantities whose behaviour during space transform is known. Let's write this golden rule:
(6) dF = 1/m Mo.f = it[P0,P1,P2].f
with m = P0.(P1^P2) and Mo = [ P1^P2, P2^P0, P0^P1 ].
The first interesting consequence of (6) appears as we compute the scalar product of dF with any of the vertex Pi defining our triangle:
Pi.dF = 1/A' tPi.[ P1^P2, P2^P0, P0^P1 ].f
= Fi
This is really great: retrieving the original values of our field F to be interpolated, once we have the gradient vector dF, is only a matter of using this scalar product. And by linearity, we can even write, with proper combinations:
(7) F(P) = P.dF
giving the interpolated value of F at *any* point P. This is not really surprising we end with such a result, and we will actually take (7) as a *definition* of the gradient. Something like being a linear form.
The thing to remember is that, once you have explicitly computed (with (6)) the gradient dF for a given F, you can throw out the original values Fi. They aren't needed anymore. If you really need them back, just use Pi.dF. You can even compute the F value on any point of the screen by re-ordering (7) in order to use the real projection rules. Namely, if you project a 3D point on screen using:
X = Xc + Focal * x/z
Y = Yc - Focal * y/z (beware of the minus sign!)
where (Xc,Yc) is the center of your screen, dividing both (7) by z and re-ordering gives:
F(X,Y) = Screen_Vector.dF
with Screen_Vector being an usual 3D-vector whose components are:
| -Xc + X
Screen_Vector = | Yc - Y
| Focal
This is particularly useful when you want to find the values of F at point X,Y on screen without having to find the real 3d-coordinates first. For instance, in a span-drawing function where you only know the starting pixel of the span, not the 3D position of this span on the polygon...
This is good, but only valid for the original triangle, not the transformed one, as far as we know. We need to examine how does the gradient dF transforms as space is linearly transformed with something like P' = M.(P + T). That's straightforward since we have all the formulae in paragraph 5. It gives:
m' dF' = Mo'.f = detM itM.[ dF + ( N^dF )^T ]
where N is the aforementioned vector B^C + C^A + A^B, the unormalized normal. Using the transform rule for the mixed product m', this eventually gives:
(8) dF' = itM.( m.dF+(N^dF)^T ) / ( m + T.N )
You will have noticed that, from our polygon-based point of view, the transform rule for the mixed product (m'=detM*(m+T.N)) used above stands as the transform rule for our polygon's area. Hence, doing backface cull for the transformed polygon is just a matter of testing the sign of m'! If m' is negative or zero, the polygon won't be seen in its new position...
The above expression (8) is still not a ready-to-use expression, because we'd have to introduce the unormalized normal N. Can we find a relation between this normal and some gradients? Yes, indeed: take the field f as being an uniform one, say 1. This only means that we are searching the gradient of w on screen, for some z-buffer drawing... With F now equal to [1,1,1], the w-gradient on screen now reads dw = 1/m. Mo.[1,1,1], from (4). But, according to the exact writing of Mo, the Mo.[1,1,1] vector is just the N vector. Hence: dw = 1/m.N. And we recover the well-known fact that the w-gradient of a polygon is simply given by its normal. Replacing N by m.dw in (8) will give the final result:
dF + T^(dF^dw)
(9-1) dF' = itM . ----------------
1 + T.dw
Although complex, relation (9-1) respects the group properties of space transforms. That is, if you perform a transform (M1,T1) followed by another transform (M2,T2), you can apply successively rule (9-1) with (M1,T1) and then with (M2,T2); or directly use the composition of the transform (M2,T2)o(M1,T1). As you can check directly by writing it down, performing (M1,T1) and then (M2,T2) simply amounts to directly performing the transform (M2.M1, T1+iM2.T2 ).
Now, one can re-arrange this a bit using A^(B^C)=(A.C)B-(A.B)C:
(9-2) dF' = itM.( dF - c*dw )
where 'c' is the scalar:
T.dF
(9-3) c = ---------
1 + T.dw
showing that, up to a correction 'c' along the dw direction due to the translation T, the gradient dF 'almost' transforms like axial vectors: dF' ~= itM.dF. In fact, the only vector that transforms according to X'= itM.X for a given transform (M,T) is: X = T^(dF^dw) / ( T.dw ).
As a special case, note that when dF is actually the gradient dw, relation (9) reads:
itM.dw
(9-4) dw' = ----------
1 + T.dw
We now have have all the tools to deal with pre-stored gradients. This is indeed a new approach we end up with: As strange as it seems, what we are actually doing is pre-projecting the gradients at scene set-up as if the camera was standing at the origin O, looking in the z-direction; and then we transform this stored pre-projection along with meshes, instead of transforming meshes and then recomputing the gradients. But, there's still a potential bug we need to look into further details...
7 - The m=0 "bug"
According to all the above derivations, the idea of pre-projecting gradients in the scene as viewed from the origin is good ... as soon as this operation is possible. And actually there's a degenerate case that can potentially ruin all this: if you look back till expression (2), you'll see that we have implicitly assume that the Mo matrix was invertible. And it actually isn't the case if the mixed product 'm' is zero, or, put otherwise, when our poly's vertex P0,P1,and P2 define a plane containing the origin. This is not surprising: How can you define and store the gradients of a polygon which projects with zero-area of screen, and thus which shouldn't be drawn? We are indeed in big trouble since our major hypothesis (A' is non-zero in expression (2)) is simply wrong.
The cure could be to tell your graphist he should not make polygons pass through the origin O while modelling his scene. Well... Let's better search for workarounds:
First possible workaround is to use an initial, 'hidden', translation that will send spurious polygons far away from this annoying origin. Not really practical...
On another hand, the Right Thing to do, in order to carefully examine this singularity, is to translate the buggy polygon of a very small quantity epsilon, making 'm' being non-zero, express all the formulae with epsilon, and let epsilon tend toward zero to see how they behave. In fact, this approach shows that information about mapping is lost during pre-projection for such polygons, and can't be retrieved later as we use brute force.
As an illustration, let's search for the possibility of having a zero-transformed gradient. Solving the necessary condition dF-cdw = 0 gives: T.N = m + T.N, which is never met, except when m=0, where it is *always* true. Here occurs the loss of information. Anyhow, only gradients proportional to dw will suffer this loss...
Let's have a closer look at this 'dw' vector. We already noticed that it's proportional to the normal and represents the gradient for the uniform field f=1. Anyway, according to the definition (7), we must always have P.dw = 1 for any P belonging to the poly's plane. This constancy involves two things:
First, if we take for P the orthogonal projection on the plane of the origin O, we deduce that, up to a minus sign:
P.dw = ||P||*||dw||
(||.|| being the norm of vectors)
and hence the norm of ||dw|| is 1/||P||.
Aieeee, what if P is of norm 0?!? You've guessed it: if the plane contains the origin, the norm of dw is infinite. In fact, it's just because the visible area of a poly passing through the origin O is zero, and we then have a null denominator in expression (2).
On a second hand, P.dw is always positive. No matter which P we choose. This mean that dw *always* points toward the half-space (whose boundary is our poly's plane) *not* containing the origin. As a consequence, dw has the same orientation as the "true" normal N only if this normal does not point toward the 0-containing half-space.
As a summary, when you translate your poly, its plane *may* cross the origin O, thus reversing the orientation of dw, and even giving it an infinite norm when it exactly contains the origin. There's a so-called singularity, and we cross this singularity as soon as 1 + T.dw is zero in expression (9).
Eventually, this singularity seems an artificial one (remember we've factorized by a quantity w0.w1.w2 while deriving expressions (2)?).
I won't go deeper in this subject, especially as it may seem of really petty concern. Anyway, it actually isn't, and trying to remove this m=0 bug is rewarded by a cool surprise: a new set of expressions for gradients that are *always* resistant to spurious cases such a m=0. It's a good exercise to try to find them, starting back from scratch. It involves splitting the 'uniform' part of the gradient (proportional to dw) and the rest. It actually is a decomposition into dw-orthogonal and dw-parallel components, this later being the only direction holding the critical part of the behaviour...
Eventually, it leads to the following final result:
dw = P0^P1 + P1^P2 + P2^P0 (the unormalized normal)
dF = dw^Df / (dw.dw)
where Df is the vector:
Df = f0*(P2-P1) + f1*(P0-P2) + f2.(P1-P0)
Transform rule (8) must then be used for dF, even if m=0...
It is rather a simple result I could have written first, but I thought it was more interesting to respect erring and questioning of research... sorry. :)
8 - Conclusion
You now are in condition to transform any basic gradient dF any time you either rotate the object or move the camera, without having to care about clipping and projection at first. There's no particular trap once you're ok with all the formulae...
As an example of all the above, you may have a look at the proggy 'grad.exe'. Gradient code is gathered in the file 'grad.c', whereas 'matrix.h' and 'video.h' are just boring funcs you should replace by your own...
The chapter isn't closed yet, and there's still a lot to find about gradients: How can we simply rotate or translate a texture's UV-coordinates? What's the impact on the stored gradients if we dynamically change the input values (f0,f1,f2) for env-mapping? What shall be changed if definition (7) for gradients is replaced by something like: F(P) = P.dF + tP.S.P + ..., S being the symmetric matrix of a quadratic form (this has something to do with quadratic mapping instead of just linear one)? Etc...
Anyway, one shall always prefer faster maths to a faster CPU... :)
Skal / BomB!
Pascal.Massimino@ens.fr
Pascal.Massimino@ptm.u-cergy.fr
http://www.eleves.ens.fr:8080/home/massimin/
http://bomb.citeweb.net/
Thu Feb 18 08:15:06 GMT 1999