Plastic BRDF brightness differences in sample and evaluate methods

[BashPrince]

Last year I discovered that the spectrum values returned from the plastic brdf are different in the sample and the evaluate methdos for the same direction. This can be seen when looking at an image rendered with path tracing with next event estimation on and off and flicking between them.

NEE on

NEE off
Here everything grey was rendered with the plastic material.

The differences appear in the diffuse component of the material (you can see the differences when setting the specular reflectance parameter for the material to 0).
I finally got around to looking into the issue and I found the causes, but I'm not sure how to fix it:
In the sample method a vector m is sampled with GGXMDF::sample. This vector is used for calculating fresnel values (for incoming and outgoing directions) as well as in the glossy if-branch for reflecting the outgoing direction to get the incoming direction wi. In the diffuse else-branch however wi is determined by sample_hemisphere_cosine. But the vector w is then still used for fresnel calculation to evaluate the spectrum value.
In the evaluate method the vector m is calculated as normalize(wi + wo) and then also used for calculating fresnel values, which are used to evaluate the specular and diffuse spectrum values. But calculating m as the halfway vector in this way does not make sense for the diffuse case where we use sample_hemisphere_cosine during sampling. The fresnel values and the Spectrum values which depend on m will therefore be different for the same incoming direction in the two methods.

Now as for how to fix this I'm not sure since I don't really know what was intended here. I looked at the cited reference (https://hal.inria.fr/hal-01386157) but I'm not sure how these routines are derived from what is described in that paper. It seems to me like sampling m to be used for spectrum values but then still using sample_hemisphere_cosine (using the same 2D [0, 1) sample again) is probably not the right approach?
Perhaps a solution would be to also calculate m as the halfway vector in the diffuse branch of the sample method and use that for fresnel? That would at least make the approach consistent.

Happy for feedback on this.

[dictoon]

Hi Stephen, thanks for looking into this in details, much appreciated!

Hopefully @est77 can shed some light here as I'm not familiar with this code. I'll ask him if he chime in.