Basic Wiener application
To have the desired effect of image denoising, the Wiener filter for can only be implemented with a point-spread function (blurring filter) that filters out the noise by averaging in some neighborhood of image pixels. As you do not disclose your program text, one cannot say for sure, but most probably your disappointing result comes from failure to apply the averaging in the vicinity of each pixel.
The typical derivation of the Wiener filter in Fourier domain arrives at the expression
$$
H(u,v) = {\frac {B^*(u,v)S_{ff}(u,v)} {|B(u,v)|^2 S_{ff}(u,v) + S_{nn}(u,v)}}
$$
where $H(u,v)$ is the Wiener filter in agreement with your designation, $B$ is the point-spread function (blurring filter), and I rename your $P$'s into $S$'s in order to reserve the character $P$ for a common PSD. In keeping with your use of $f$ subscripts for the original image and $g$ for this image with WGN added, $S_{ff}(u,v)$ is the (supposedly unknown) power spectrum of the original image before the noise was added.
Using the periodogram technique to estimate the power spectrum from the DFT of the observation $X(u,v)$, one can write $S^{periodogram} = |X(u,v)|^2/N^2$ and then arrive at the expression similar to your $H(u,v) = (P_g - P_n)/P_g$:
$$
H(u,v) = {1 \over B}{\frac {S^{periodogram}(u,v) - S_{nn}}{S^{periodogram}(u,v)}}
$$
with the remarkable appearance of the ${1 \over B}$ factor.
Now, there can arise problems with the inverse of the point-spread function $B(u,v)$. The solution is to use matrix pseudoinverse.
Back to your question: you failed to reconstruct a denoised image because you omitted the averaging in the pixel neighborhoods. Experimentation with MATLAB function wiener2(I,[m n]) for the pixel-wise adaptive low-pass Wiener filtering with different sizes of the pixel neighborhood would help you understand the working of the Wiener filter on noisy images. You would see how the choice of optimal m,n values depends on the local image variance.