Hello Greg,
Thank you for your input. I think I figured it out, thanks to "Microfacet Models for Refraction through Rough Surfaces" by Bruce Walter, Stephen R. Marschner, Hongsong Li, and Kenneth E. Torrance.
They state explicitly that the sampling is done from D(W_h)cos(theta_h). That makes sense because this cosine factor is needed to make the PDF integrate to 1. In the book this formula is mentioned in
8.4.2.
And then it goes as I described it above - the results from
13.5.3 are used to convert it to spherical coordinates and then formula
13.12 is applied to arrive with PDF for theta.
Well, math is not my strong side all this could be wrong but sounds convincing enough for me to calm me down:)
Regards,