I have created a 3D plot and a contour plot for the same (Wigner) function, but as you see the contour plot is not right. It has some "extra" regions in gray color (ovals). What is going wrong? The 3D plot code is

Plot3D[2/((2 + 2/E^8) \[Pi]) (Exp[-2 (4 - 4 x + x^2 + y^2)] + 
Exp[-2 (4 + 4 x + x^2 + y^2)] + 
2 Exp[-2 (x^2 + y^2)] Cos[8 y]), {x, -4, 4}, {y, -4, 4}, 
PlotRange -> Full, Exclusions -> None, PlotPoints -> 180, 
FaceGrids -> {{0, -1, 0}, {1, 0, 0}, {0, 0, -1}}, 
FaceGridsStyle -> Directive[Blue, Dotted], Boxed -> False, 
BoxRatios -> {1, 1, 1}, 
AxesLabel -> {"x", "y", "\!\(\*SubscriptBox[\(W\), \(c\)]\)(x,y)"}, 
LabelStyle ->   {Bold, Black, FontSize -> 16}]

The contour plot code is

ContourPlot[
2/((2 + 2/E^8) \[Pi]) (Exp[-2 (4 - 4 x + x^2 + y^2)] + 
Exp[-2 (4 + 4 x + x^2 + y^2)] + 
2 Exp[-2 (x^2 + y^2)] Cos[8 y]), {x, -3, 3}, {y, -3, 3}, 
PlotLegends -> Automatic, 
PlotPoints -> 180]

And here are the figures:

The real contour plot is taken from a textbook (Quantum Optics by Girish Agarwal, Cambridge U Press, 2013, page 80):

To emphasize that it is not a problem for a special case, here I have added another code and figure. You can see the "ugly" ovals are definitely incorrect.

ContourPlot[
1/((2 + 2/E^16) \[Pi])
2 (E^(((2 + 2 I) - x - I y) ((-2 + 2 I) + x - I y)) + 
E^(((-2 - 2 I) - x - I y) ((2 - 2 I) + x - I y)) + 
E^((-x - I y) (x - I y)) Cos[4 Im[(2 - 2 I) (x + I y)]]), {x, -4, 
4}, {y, -4, 4}, PlotLegends -> Automatic, 
ColorFunction -> "Rainbow", FrameLabel -> {"x", "y"}, Frame ->  True,
FrameStyle -> Directive[Black, 18], PlotPoints -> 180, 
PlotRange -> All, PlotRangePadding -> 0]