You can sample using absolutely any distribution you want, as long as you weight the results by dividing by the pdf of the sampled distribution. It will converge to the right answer (as long as the distribution is nonzero everywhere that you want to integrate). Different distributions will give different amounts of variance though. The trick with importance sampling is to find a distribution that minimizes variance while still being cheap to compute. A distribution similar to the shape of the function you're integrating, or equal to some factor of that function (so that it cancels when you divide), works best. In other words you try to guide the sampling toward areas that are more important to the result, hence, "importance" sampling.
So what the author is saying is the image on the left was sampled with $\omega_o$ (or maybe $\omega_m$, I'm not sure) drawn from a cosine hemisphere distribution, just as an example of a distribution that does not match the BRDF very well. Then they calculate the samples as $f(\omega_i, \omega_o) / p_{\cos}(\omega_o)$ where $p_{\cos}$ is the pdf for the cosine distribution. This is just to show you that it gives the "correct" result, but with a ton of variance. (Although it looks like there are some black areas on the left image that suggest it's not really converging to the correct answer...but let's assume that's a bug.)
Then the image on the right is sampled with (presumably) $\omega_m$ drawn from a distribution matching the $D(\omega_m)$ term from the BRDF. Then when you calculate the BRDF you can leave out the $D$ factor since it's canceled by the $D$ factor in the pdf. (The pdf probably still has some normalization and some conversion of probability from $\omega_m$ to $\omega_o$, so you still have to divide by those factors.)