This is basically a call for anyone more mathematically literate than myself to verify my results.
I want my shader system to be designed such that it is easy to mix-and-match various specular models and diffuse models while retaining proper conservation between specular and diffuse.

The way I have set it up is so that the specular function returns a factor by which to multiply the diffuse in order to enforce conservation.
So if you look at tri-Ace’s research papers on physically plausible Blinn-Phong, you see that they include a second Fresnel factor (using L·N) in the diffuse component. Since that factor is designed to ensure energy conservation, I actually don’t want to make it part of the diffuse model, but rather the diffuse multiplier returned by the specular model, because that is basically saying how much is left over for diffuse after specular is determined.

In other words, I don’t want to do this:

void BlinnPhongDir( in DEF_LIGHT_ARGS _dlaArgs, inout vec4 _vDiffuse, inout vec4 _vSpecular ) {
	// Diffuse normalization (division by PI) done once outside of the BRDF's.
	vec4 vDiffuse;
	vec4 vSpecular;

	// *********************
	// SPECULAR
	// *********************
	float fS = 0.0397436 * _dlaArgs.fShininess + 0.0856832;
	
	float fSpec = SchlickFresnel( _dlaArgs.fSpecReflectance, _dlaArgs.fVdotH ) * pow( _dlaArgs.fNdotH, _dlaArgs.fShininess );
	fSpec /= max( _dlaArgs.fLdotN, _dlaArgs.fVdotN );
	vSpecular = max( fSpec * fS, 0.0 ) * _dlaArgs.vLightDiffuse * _dlaArgs.fLdotN;


	// *********************
	// DIFFUSE
	// *********************	
	float fD = saturate( (1.0 - SchlickFresnel( _dlaArgs.fSpecReflectance, _dlaArgs.fLdotN )) );
	vDiffuse = _dlaArgs.fLdotN * fD * _dlaArgs.vLightDiffuse;
	
	_vDiffuse += vDiffuse;
	_vSpecular += vSpecular;
}

I want to keep them separated this way so that they can be interchanged with other specular/diffuse models:

float BlinnPhongSpecular( in DEF_LIGHT_ARGS _dlaArgs, out vec4 _vSpecular ) {
	float fS = 0.0397436 * _dlaArgs.fShininess + 0.0856832;
	
	float fSpec = SchlickFresnel( _dlaArgs.fSpecReflectance, _dlaArgs.fVdotH ) * pow( _dlaArgs.fNdotH, _dlaArgs.fShininess );
	fSpec /= max( _dlaArgs.fLdotN, _dlaArgs.fVdotN );
	_vSpecular = max( fSpec * fS, 0.0 ) * _dlaArgs.vLightDiffuse * _dlaArgs.fLdotN;
	
	float fD = saturate( (1.0 - SchlickFresnel( _dlaArgs.fSpecReflectance, _dlaArgs.fLdotN )) );
	//float fD = 1.0 - _dlaArgs.fSpecReflectance;
	return fD;
}

void BlinnPhongDiffuse( in DEF_LIGHT_ARGS _dlaArgs, out vec4 _vDiffuse ) {
	LambertDiffuse( _dlaArgs, _vDiffuse );
	//BurleyDiffuse( _dlaArgs, _vDiffuse );
}

void BlinnPhongDir( in DEF_LIGHT_ARGS _dlaArgs, inout vec4 _vDiffuse, inout vec4 _vSpecular ) {
	// Diffuse normalization (division by PI) done once outside of the BRDF's.
	vec4 vDiffuse;
	vec4 vSpecular;
	float fDFactor = BlinnPhongSpecular( _dlaArgs, vSpecular );
	BlinnPhongDiffuse( _dlaArgs, vDiffuse );
	vDiffuse *= fDFactor; // DIFFUSE REDUCED BY AMOUNT LOST TO SPECULAR HERE.
	
	_vDiffuse += vDiffuse;
	_vSpecular += vSpecular;
}

This way I can do this:

// *********************
// Oren-Nayar/Blinn-Phong.
// *********************
void OrenNayarBlinnPhongDir( in DEF_LIGHT_ARGS _dlaArgs, inout vec4 _vDiffuse, inout vec4 _vSpecular ) {
	// Diffuse normalization (division by PI) done once outside of the BRDF's.
	vec4 vDiffuse;
	vec4 vSpecular;
	float fDFactor = BlinnPhongSpecular( _dlaArgs, vSpecular );
	OrenNayarDiffuse( _dlaArgs, vDiffuse );
	vDiffuse *= fDFactor; // DIFFUSE REDUCED BY AMOUNT LOST TO SPECULAR HERE.

	_vDiffuse += vDiffuse;
	_vSpecular += vSpecular;
}

I want to apply this also to Ashikhmin-Shirley.
Here is how it would be if you separated the specular and diffuse in the “typical” way:

void AshikhminShirleySpecular( inout DEF_LIGHT_ARGS _dlaArgs, out vec4 _vSpecular ) {
	const vec3 vEpsilon = vec3( 1.0, 0.0, 0.0 );
	vec3 vTangent = normalize( cross( _dlaArgs.vNormal, vEpsilon ) );
    vec3 vBiTangent = normalize( cross( _dlaArgs.vNormal, vTangent ) );
	float HdotX = dot( _dlaArgs.vLightHalfVector, vTangent );
    float HdotY = dot( _dlaArgs.vLightHalfVector, vBiTangent );
	
	// _dlaArgs.vAshikFactors.w = Rs.
	float fF = SchlickFresnel( _dlaArgs.vAshikFactors.w, _dlaArgs.fVdotH );
	// _dlaArgs.vAshikFactors.y = sqrt((Nu+1)*(Nv+1)).
	float fNormS = _dlaArgs.vAshikFactors.y * (1.0 / (8.0 * PI));
	float fN = (_dlaArgs.vAshikFactors.x * Sqr( HdotX ) + _dlaArgs.vAshikFactors.y * Sqr( HdotY )) / (1.0 - Sqr( _dlaArgs.fNdotH ));
	float fDenom = max( (_dlaArgs.fVdotH * max( _dlaArgs.fVdotN, _dlaArgs.fLdotN )), LSE_EPSILON );
	float fFinal = max( fNormS * fF * pow( _dlaArgs.fNdotH, fN ) / fDenom, 0.0 );
	
	_vSpecular = (fFinal * _dlaArgs.fLdotN) * _dlaArgs.vLightDiffuse;
}

void AshikhminShirleyDiffuse( in DEF_LIGHT_ARGS _dlaArgs, out vec4 _vDiffuse ) {
	float fFactorV = 1.0 - Pow5( 1.0 - _dlaArgs.fVdotN * 0.5 );
	float fFactorL = (1.0 - Pow5( 1.0 - _dlaArgs.fLdotN * 0.5 ));
	// (23.0 * PI) changed to just 23.0 because diffuse is divided by PI later.
	_vDiffuse = (1.0 - _dlaArgs.vAshikFactors.w) * (28.0 / 23.0) * (fFactorV * fFactorL) * _dlaArgs.fLdotN * _dlaArgs.vLightDiffuse;
}

This is how the specular and diffuse terms are isolated in the papers.

For my purposes I need to extract the part of the diffuse that is actually meant to conserve energy from the specular, and this is my result:

float AshikhminShirleySpecular( inout DEF_LIGHT_ARGS _dlaArgs, out vec4 _vSpecular ) {
	const vec3 vEpsilon = vec3( 1.0, 0.0, 0.0 );
	vec3 vTangent = normalize( cross( _dlaArgs.vNormal, vEpsilon ) );
    vec3 vBiTangent = normalize( cross( _dlaArgs.vNormal, vTangent ) );
	float HdotX = dot( _dlaArgs.vLightHalfVector, vTangent );
    float HdotY = dot( _dlaArgs.vLightHalfVector, vBiTangent );
	
	// _dlaArgs.vAshikFactors.w = Rs.
	float fF = SchlickFresnel( _dlaArgs.vAshikFactors.w, _dlaArgs.fVdotH );
	// _dlaArgs.vAshikFactors.y = sqrt((Nu+1)*(Nv+1)).
	float fNormS = _dlaArgs.vAshikFactors.y * (1.0 / (8.0 * PI));
	float fN = (_dlaArgs.vAshikFactors.x * Sqr( HdotX ) + _dlaArgs.vAshikFactors.y * Sqr( HdotY )) / (1.0 - Sqr( _dlaArgs.fNdotH ));
	float fDenom = max( (_dlaArgs.fVdotH * max( _dlaArgs.fVdotN, _dlaArgs.fLdotN )), LSE_EPSILON );
	float fFinal = max( fNormS * fF * pow( _dlaArgs.fNdotH, fN ) / fDenom, 0.0 );
	
	_vSpecular = (fFinal * _dlaArgs.fLdotN) * _dlaArgs.vLightDiffuse;
	
	// According to the paper, (28.0 / (23.0 * PI)) is a magic number designed to ensure energy conservation.
	return (1.0 - _dlaArgs.vAshikFactors.w) *		// Constant ratio of specular to diffuse.
		(28.0 / 23.0) *					// (23.0 * PI) changed to just 23.0 because diffuse is divided by PI later.
		(1.0 - Pow5( 1.0 - _dlaArgs.fLdotN * 0.5 ));	// Light lost from diffuse by increased specular at glancing angles.
}

void AshikhminShirleyDiffuse( in DEF_LIGHT_ARGS _dlaArgs, out vec4 _vDiffuse ) {
	float fFactorV = 1.0 - Pow5( 1.0 - _dlaArgs.fVdotN * 0.5 );
	_vDiffuse = (fFactorV * _dlaArgs.fLdotN) * _dlaArgs.vLightDiffuse;
}

void AshikhminShirleyDir( in DEF_LIGHT_ARGS _dlaArgs, inout vec4 _vDiffuse, inout vec4 _vSpecular ) {
	// Diffuse normalization (division by PI) done once outside of the BRDF's.
	vec4 vDiffuse;
	vec4 vSpecular;
	float fDFactor = AshikhminShirleySpecular( _dlaArgs, vSpecular );
	AshikhminShirleyDiffuse( _dlaArgs, vDiffuse );
	vDiffuse *= fDFactor; // DIFFUSE REDUCED BY AMOUNT LOST TO SPECULAR HERE.
	
	_vDiffuse += vDiffuse;
	_vSpecular += vSpecular;
}

Can anyone verify this is correct?
28.0 / 23.0 is a magic number I don’t know how was obtained, but the paper says it is for conservation so I put it in my conservation part.
1.0 - _dlaArgs.vAshikFactors.w (1 - Rs) I know is part of the exchange between specular and diffuse.
1.0 - Pow5( 1.0 - _dlaArgs.fLdotN * 0.5 ) I moved to the conservation area because it appears to be the part that compensates for increasing specular at glancing angles.

The remaining part I left in the diffuse term because having diffuse decrease at glancing angles to the viewer is a property that specifically Ashikhmin-Shirley diffuse wants. It’s not related to conserving energy with specular.

The basic question I would like anyone better at math than myself to answer is, if I want to mix Ashikhmin-Shirley specular with Oren-Nayar diffuse while retaining energy conservation (as I did with Blinn-Phong/Oren-Nayar), is this correct?


L. Spiro

tl;dr. Post the (target) BRDF equations (using math-notation) instead of shader code. Integrate your BRDF and see if Diffuse + Specular <= 1 applies. Or just look at some discrete cases (including N=V) and see if it breaks.