Using your idea, define [itex]u=(e^{ik(x-x_0)+il(y-y0)}+e^{-ik(x-x_0)-il(y-y0)})/2[/itex] and [itex]v=(e^{-ik(x-x_0)-il(y-y0)}+e^{-ik(x-x_0)+il(y-y0)})/2[/itex], then the intergration region A looks simpler: [itex] A: u+v>b[/itex] , but then [itex]\cos(k(x+y))\cos(ly)[/itex] is not possible to write as a function of u and v.

Is there other transform help to have a better shape of region A?
The main difficulty to directly calculate the double integration is that when I break the double integration into two 1-d integration, first x then y, then the integration for x is over a y dependent range, and this in turn makes the integration for y very complicated.