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:
[![enter image description here][1]][1]
[![enter image description here][2]][2]


  [1]: https://i.sstatic.net/mPaWn.jpg
  [2]: https://i.sstatic.net/hiCfC.jpg