A kludgey implementation that struggles at some of the “points” of the zones:

Define a function that takes an equation defining the boundaries and turns it into an inequality satisfied if a point is on the same side of the boundary as the origin is.
PerpBisIneq(pt1_, pt2_) := Block({eq = PerpBis(pt1, pt2)}, If(eq === True, True, If((First(eq) /. {x > 0, y > 0}) > 0, Greater @@ eq, Less @@ eq)) )

Define a function that returns a list of inequalities as functions of the point
{x,y}
.ineqs({x_, y_}) = Flatten(Table(PerpBisIneq({0, 0}, {i, j}), {i, 2, 2}, {j, 2, 2}))

Create a
RegionPlot
defined by the region of the plan in which precisely $n1$ elements ofineqs
areFalse
.brillouinZone(n_) := RegionPlot( Count(ineqs({x, y}), False) == n  1, {x, 5, 5}, {y, 5, 5}, PlotStyle > ColorData(109)(n), PlotPoints > 100)

Plot the Brillouin zones.
Show( Table(ContourPlot(Evaluate@PerpBis({0, 0}, {i, j}), {x, 5, 5}, {y, 5, 5}), {i, 2, 2}, {j, 2, 2}), Table(brillouinZone(n), {n, 1, 5}) )
A few notes on how this could be refined:

The equation for
i = j = 0
always evaluates asTrue
, requiring the first layer of casechecking inPerpBisIneq
. The second layer of casechecking is not necessary ifgamma
is always negative, but I’m not sure if this is the case. 
The setting for
PlotPoints
in the definition ofBrillouinZone
is necessary to get Mathematica to return a “nice” plot. Lower settings cause it to miss some zones entirely. 
There is probably a nicer way of getting Mathematica to logically reduce the condition “precisely $n$ of these inequalities are
False
” than the way that I’ve done it, in which Mathematica is (I think) effectively sampling the plane pointwise. I suspect there’s a way to do it withReduce
or something similar, but I’m not immediately sure what it is. 
The color scheme can be changed by substituting another number for
109
in the definition ofBrillouinZone
.
EDIT: I was able to use Reduce
to find a set of logical predicates for the nth Brillouin zone. But it’s ugly as sin and slow as molasses, so I can’t say I recommend it. It might, however, be helpful if you ever need to rigorously define a Brillouin zone (for, say, some kind of integration), so I’m including it for reference.
Brillouinpred({x_, y_}, n_) :=
Reduce(Or @@ Apply(And, Subsets(ineqs({x,y}),{Length(ineqs({x, y}))  (n + 1)}), {1}) &&
! (Or @@ Apply(And, Subsets(ineqs({x,y}), {Length(ineqs({x, y}))  n}), {1})))
This returns a logical condition on x
and y
which is True
if {x,y}
is in the nth Brillouin zone and False
otherwise. This is done by noting that if we have a set of $N$ predicates $P = {P_1, P_2, …, P_N}$, then the statement that “no more than $m$ of the predicates are true” is equivalent to
$$
mathcal{P}_m = bigvee_{Q subset P,, Q = Nm} left( bigwedge_{P_i in Q} P_i right).
$$
The outer “or” runs over all subsets of cardinality $Nm$, and the inner “and” runs over the elements of each such subset. One can then construct the predicate “exactly $m$ predicates are true” by demanding that $mathcal{P}_m wedge neg mathcal{P}_{m1}$. The code above explicitly constructs this statement and then applies Reduce
to it for simplification purposes. However, note that this creates a set of ${N choose m}$ logical predicates, each containing $Nm$ subpredicates, as an intermediate step; so Mathematica takes some time to solve this.