Introduction
Recently, I was implementing skin shading in my engine. The pre-integrated skin shading technique is chosen because it has a low performance cost and does not require another extra pass. The idea is to pre-bake the scattering effect over a ring into a texture for different curvature to look up during run-time. More information can be found in the GPU Pro 2, the SIGGRAPH slides, and also in the presentation of the game "The Oder:1886". Here is the result implemented in my engine (all screen shots are rendered with filmic tone mapping):

The head on the left is lit with Oren-Nayar shading
The head on the right is lit with pre-integrated skin

Curve Approximation for Direct Lighting
In my engine, iOS is one of my target platform, which only have 8 texture unit to use for the OpenGL ES 2.0 API. This is not enough for the pre-integrated skin look up texture because my engine have already used some slots for light map, shadow map, IBL... So I need to find an approximation to the look up texture.

Unfortunately I don't have a Mathematica license at home. So I think may be I can fit the curve manually by inspecting the shape of the curve. So, I started by plotting the graph of the equation:

Here is the shape of the red channel diffusion curve, plotting with N.L(normalized to [0, 1])and r (from 2 to 16):

My idea for approximating the curve is by finding some simple curves first and then interpolate between them like this:

For the light blue line in the above figure, a single straight line can get a close enough approximation:

curve1= saturate(1.95 * NdotL -0.96)

To approximate the dark blue line, I divide it into 2 parts: linear part and quadratic part:

curve0_linear= saturate(1.75* NdotL -0.76)

curve0_quadratic= 0.65*(NdotL^ 2) + 0.045

blending both linear and quadratic curve will get a curve that is similar to the original function:

curve0= lerp(curve0_quadratic, curve0_linear, NdotL^2)

Now we have 2 curves that is similar to our original function at both end. By mixing them together, we can get something similar to the original function like this:

curve= lerp(curve0, curve1, 1 - (1 - curvature)^4)

By repeating the above steps for the blue and green channels, we can shade the pre-integrated skin without the look up texture. Here is the result:

Lit with a single directional light
From left to right: r = 2, 4, 8, 16

Upper row: shaded with look up texture

Lower row: shaded with approximated function

This is how it looks like when applying to a human head model:

From left to right: shaded with look up texture, approximated function, lambert shader

Upper row: shaded with albedo texture applied

Lower row: showing only lighting result

For your reference, here is the approximated function I used for the RGB channels:

Curve Approximation for Indirect Lighting
In my engine, the indirect lighting is stored in spherical harmonics up to order 2. So the pre-integrated skin BRDF need to be projected into spherical harmonics coefficients and can be store in look up texture just like the presentation of "The Oder:1886" described. One thing to note about the integral range in equation (19) from the paper should be up to π instead of π/2 because to project the coefficients, we need to integrate over the whole sphere domain and the value of D(θ, r) in the lower hemi-sphere may be non zero for small r due to sub-surface scattering, which does not like the clamped cos(θ). So I compute the spherical harmonics coefficient using this:

To make the indirect lighting work on my target platform, an approximate function to this indirect lighting look up texture also need to be found. By using similar methods described above and with some trail and error. Here is my result:

Lit with both directional light and BRDF projected into SH
From left to right: r = 2, 4, 8, 16

Upper row: shaded with look up texture

Lower row: shaded with approximated function

And applying it to the human head model, but this time the approximation is not close enough and lose a bit of red color:

From left to right: shaded with look up texture, approximated function, lambert shader

Upper row: shaded with albedo texture applied

Lower row: showing only lighting result

And some code for your reference, where the ZH is the zonal harmonics coefficient:

Conclusion
In this post, I describe a way to find an approximated function for the pre-integrated skin diffusion profile, which gives a similar result for the direct lighting function while losing a bit of red color for the indirect lighting. The down side of fitting the curve manually is when the function is changed a bit, say changing the function input from radial distance to curvature(i.e. from r to 1/r), all the approximate functions need be re-do again (or the conversion need to be done during run-time just like my code snippet above...). Also the shadow scattering described in the original paper is not implemented, so some artifact may be seen at the direct shadow boundary. Overall, the skin shading result is improved compare to shade with Lambert or Oren-Nayar under environment with a strong directional light source.