Energy band plotting

Is there some way in which I could plot energy bands using Mathematica in a manner that is given in the image below? I have $E$ as a function of $k_x$ and $k_y$. How can the plot given below be made with these specific points on the $(k_x,k_y)$ axis?

My function for the energy $E$ is:

energy[kx_,ky_]:=Sqrt[1-(Cos[kx/2]Cos[ky/2])^2]

and I can plot parts of the graph with something like:

Plot[energy[kx,0],{kx,0,Pi}]
Plot[energy[kx,0],{kx,Pi,2 Pi}]

How do I stitch these together?

Answer

Construct an Interpolation function that translates from the $x$-value to the $\{k_x,k_y\}$ coordinates:

xkData = {{0,{Pi,0}},{1,{0,0}},{2,{Pi,-Pi}},{3,{Pi,0}},{4,{2Pi,0}}};
if = Interpolation[xkData, InterpolationOrder->1];

Then, use the above interpolating function to construct the plot, and use xkData again to create the ticks:

Plot[
    energy @@ if[x],
    {x, 0, 4},
    Frame -> True,
    FrameTicks->{

        {Automatic,None},
        {xkData, None}
    },
    GridLines->{Range[0,4], None}
]