There are several misunderstandings:
normalize(normal + uniformSphere())
: I don't think this is the correct way to implement cosine-weighted hemisphere sampling. Check this out: Sampling the hemisphere and 13.6.3 Cosine-Weighted Hemisphere Sampling (this is more involved). The method you implemented by normalizing the sum of normal and uniformly sampled vector is neither uniform hemisphere sampling nor cosine-weighted sampling.
- PDF is not
cos(theta)
. PDF is cos(theta) / pi
. If you integrate the PDF over the hemisphere then you will find that the latter one results in 1 (valid). Don't forget to add the term $\sin \theta$ to account for the measure conversion between solid angle and polar coords.
Getting brighter results is caused by not using unbiased sampling method and PDF pairs. Let's break it down (a little bit of math, but I will try to make it easy to understand):
MC integration is aimed at approximating the following integral in an unbiased way:
$$
\int_{\Omega} f(x)d\mu(x)\tag{1}
$$
That is, for the following estimator (let's consider a one-sample case, N = 1):
$$
\hat{I} = \underbrace{f(x)}_{\text{evaluation}} / \underbrace{p(x)}_{\text{some PDF you choose to divide}}\tag{2}
$$
If we take its expectation:
$$
\mathbb{E}(\hat{I}) = \int_{S}\frac{f(x)}{p(x)}p_{act}(x)d\mu(x)\tag{3}
$$
We want $p(x)$ (the PDF you choose as the denominator) and $p_{act}(x)$ (the actual PDF of the samples, determined by the sampling method you use) can cancel each other out, so that the result will be exactly $(1)$. Now that we know this, for the specific case of yours:
If you choose $p(x)$ as $\cos$ or even $\cos / \pi$: since the reflected direction given by your code is not drawn from the cosine-weighted PDF, $p_{act}$ won't cancel itself out with $p$ in the denominator, which results in bias (hence, brighter). Another extreme example is that if you choose $\cos /\pi$ as $p$ and use uniform hemisphere sampling, $p$ and $p_{act}$ won't cancel each other out, either. So, to get the desired output, always choose sampling method ($p_{act}$) and the denominator $p(x)$ wisely.
What you should do here, is to:
- Implement the local reflection direction sampling, according to Cosine-weighted hemisphere.
- Rotate the local vector to the global frame: since you know the normal, and your cosine-weighted sample is drawn w.r.t to vector (0, 0, 1): calculate the rotation from (0, 0, 1) to your normal and apply the rotation to your local vector.
- Devide $f(x)$ by $\cos\theta / \pi$, the $\cos\theta$ term should cancel the foreshortening term in $f(x)$, so you will end up with a simple path throughput: $k_d$.