Hi everybody,
I started to code a Lafortune shader this week, since I’ve got some time.(I’m not sure I’m in the good section). And I’m encountering some various problems 
Before explaining them, I’ll try to give an overview of what is going on with this shading technique (otherwise, it’s hard to understand the problems I’m fighting against
)
Everything is explained in this paper : http://www.graphics.cornell.edu/pubs/1997/LFTG97.pdf
However, if you don’t want to read it, here’s the concept : the Phong specular shader computes the cosine between the view vector (v) and the vector of the incident light (l), mirrored around the surface normal vector (n). The cosine of this angle (a on the picture below) is equal to the dot product between r and v, provided those vectors are normalized, of course.

Then, one elevates this cosine to a chosen power, depending on the wanted roughness of the material to be rendered : the higher this power (called n in the following lines) is, the harsher the highlight will look. This is the hardness parameter of Blender (hardness is an integer between 1 and 512 though, but it’s basically the same). One can also multiply this result by a number between 0 and 1, which could be the albedo of the material (in Blender, the “Spec” factor is between 0 and 200%… :D)
It gives this : f(r,v)=albedocos(r,v)^n = albedo(xrxv + yryv + zr*zv)^n
The mirroring around the normal n can be done using a householder matrix (mentioned in the paper). I do not understand what it is… We even didn’t learn anything about matrices in school, so, I’m not used to handle with them. I only know the very basics 
What’s new, is that, by replacing the householder matrix by a general matrix, some calcul lines later, you’ve got the dot product l.v multiplied by a diagonal matrix 33, whose coefficients can be chosen. It means you can give a different weight to the terms of the dot product, by varying a Cx, a Cy, and a Cz coefficient. The l and v vectors must be expressed in local coordinates, with z axis aligned to the surface normal n, and x/y axis aligned to the principal directions of anisotropy (I won’t be able to code anisotropy I think… It’s a lot of geometry, in space what’s more…
)
To sum up, here is the formula : f(r,v)=albedo(Cx * xlxv + Cy * ylyv + Cz * zl*zv)^n
(as I said, I won’t code anisotropy, so Cx = Cy : the vectors are equally deformed along x and y axis).
This shading method allows to get a more realistic rendering of the highlights, especially at grazing angles, for which the ordinary specular shaders aren’t convincing at all :

(top left, Phong specular shader, top right, raytraced soft reflections, bottom left, Lafortune specular shader… The weakness of the Phong shader is obvious here. Even with a Ward isotropic shader, or a Cook Torrance shader, the shape of the highlights is badly described… 
Thas was the theory, here’s my code now. I commented everything to be commented, for my code to be clear. I just show the shading function : the others lines of code I added do not need to be shown (adding buttons in the interface, that kind of thing)
/* Lafortune isotropic spec */
static float Lafortune_Spec( float *n, float *l, float *v, int hard , int tangent, float Cxy, float Cz )
{
float rslt ;
float x[3] ;
float y[3] ;
//find x axis perpendicular to surface normal : x[0]*n[0]+x[1]*n[1]+x[2]*n[2]=0 //
if (n[0] == 0.0)
{
x[0] = 0.0 ;
x[1] = n[2] ;
x[2] = -n[1] ;
}
else if (n[1] == 0.0)
{
x[0] = n[2] ;
x[1] = 0.0 ;
x[2] = -n[0] ;
}
else
{
x[0] = n[1] ;
x[1] = -n[0] ;
x[2] = 0.0 ;
}
Normalize (x);
//find y axis perpendicular to x axis and surface normal, with cross product between n and x//
y[0] = n[1] * x[2] - n[2] * x[1] ;
y[1] = n[2] * x[0] - n[0] * x[2] ;
y[2] = n[0] * x[1] - n[1] * x[0] ;
Normalize (y) ;
//get view vector in new local coordinates system//
v[0] = v[0]*x[0] + v[1]*x[1] + v[2]*x[2] ; // = v.x = cos (v, x) since x and v are normalized
v[1] = v[0]*y[0] + v[1]*y[1] + v[2]*y[2] ; // = v.y = cos (v, y)
v[2] = v[0]*n[0] + v[1]*n[1] + v[2]*n[2] ; // = v.n = cos (v, n)
Normalize (v) ;
//get light vector in new local coordinates system//
l[0] = x[0]*l[0] + x[1]*l[1] + x[2]*l[2] ; // = l.x = cos (l, x) since x and l are normalized.
l[1] = y[0]*l[0] + y[1]*l[1] + y[2]*l[2] ; // = l.y = cos (l, y)
l[2] = n[0]*l[0] + n[1]*l[1] + n[2]*l[2] ; // = l.n = cos (l, n)
Normalize (l) ;
rslt = Cxy*l[0]*v[0] + Cxy*l[1]*v[1] + Cz*l[2]*v[2] ;
// weighted dot product between light vector and view vector //
if (tangent) rslt = sasqrt( 1.0f - rslt*rslt ) ; // in case tangent shading is turned on
if( rslt > 0.0f ) rslt= spec(rslt, hard); // take the hardness parameter into account
else rslt = 0.0f; // clamp negative values to 0
return rslt;
}
To sum up : I first find the local x axis which has to be perpendicular to the normal n (the code works quite simply : imagine n is : (3, 4, 5), it will give x (-4, 3, 0), so that n.x = 0). Then, thanks to the cross product between x and n, I find the y axis perpendicular to both x and n. I can then project l and v on x and y founded before, and on n. Doing so, I’ve got l and v in local coordinates, just like in the paper. Finally, I compute the weighted dot product between l and v, in which Cxy weights the x and y terms of the dot product, and Cz weights the z term of the dot product. My coordinates transform is quite strange, but by intuition, it should work. What do you guys think of it ? It could and would surely be useful to use matrices for this, in order for the code not to be burdened…
Now, here are my problems (finally
). I would say I’ve got 3 problems, and they are related I think…
1- I wonder how vectors are… handled in blender. For instance : I setup a scene with a perfectly vertical Sun, toward the ground, and a perfectly horizontal plane, directed toward the sky, so, its normal is vertical, toward the sky. The Sun lamp throws rays from an infinite distance, so, they are all parallel ; the light vector is thus constant.
I decide to display xn, yn, zn, xl, yl, and zl in the console, and here’s what I get :
> n(-0.010817, 0.895343, 0.445245)
> l(0.010817, -0.895343, -0.445245)
It’s reassuring and strange at the same time… Reassuring because the vectors are actually opposite, but strange because I can’t figure out how vertical vectors can have such weird xyz components…
2- I wonder if my coordinates change is operational. Logically, yes, but the point is that at the moment, what I have IS NOT a lafortune shader…
It must be related to problem #1…
Example : we are in a local coordinates system, with n = z axis. Which means, by choosing Cxy = -1 and Cz = 1, we should get a Phong shader : setting -1 for Cxy will mean mirroring the l vector around the normal, and this is what Phong does before computating the cosine between the 2 vectors. And, with those precise parameters, I get some nice speculars, but not a Phong shader, without a doubt.
And yet, some renders that I get give nice things, which look like what the Lafortune shader is capable of :

So, I was wondering whether my coordinates change is valid or not.
3- The third problem is completely different from the 2 others : the specular highlights do not add correctly. With images :

As we see, when shadows are activated, it’s better, but still hardly bugged
When an area is shadowed, blender seems to consider that speculars cannot be seen anyway, so it allows to show the highlights resulting from others lamps (picture on the right). But when the shadows are turned off, even an intensity of 0 for one specular highlight is enough to overwhelm the others. I’m not sure that it’s clear, that’s why I showed those images^^. (By the way, here we can see that the shading isn’t properly done… On the left image, it looks like tangent shading is turned on…)
I first thought that it didn’t come from the shading function itself, but rather from a wrong declaration somewhere, that kind of thing. But actually, it comes from the inside of the shading function I show here : I copied and pasted the code for the Phong shader into my own function, and everything worked perfectly… So, I have absolutely no idea of what is happening wrong.
That’s all folks, if you came at the end of this post, congrats
I hope I’ll find a response to my problems, 'cause the few promising results I got give me the will to go on… 
…
and for fun, a useless video with my shading function : http://www.vimeo.com/1441982
oh, btw, I posted exactly the same thing on a french forum, and I was asked why I’m not using PyNodes instead, which fit my needs exactly… The problem is that I don’t know anything about Python 
I already took a look at those pages though… I even wrote the Lambert shader. I understood what I was coding, but I couldn’t invent the functions myself
Because I don’t know anything about Python, neither do I know how python scripts are integrated to blender.