The tricky part is the surface color and anisotropy to albedo mapping. We simulate scattering for various different albedos and anisotropies, and then do a 2D function fit to invert that.

See section 6 of Path Traced Subsurface Scattering using Anisotropic Phase Functions and Non-Exponential Free Flights.

They do a two step function fit, but I couldn't get that working. For each anisotropy I got randomly different parameters that each fit the data pretty well, but if you try to fit a function to the parameters across all anisotropies it's much too noisy to fit anything useful.

So I fit the entire 2D function in one go, which may not have the best accuracy. It seems to work except for very high anisotropy, I need to look at that more closely. If anyone thinks it's fun to find a better fit feel free to try, the program to generate the data and do the fit is in the patch.

I haven't tried testing performance or cleaning/optimizing the code yet.