image processing – Color regions defined by many lines

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

  1. 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))
  2. 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}))
  3. Create a RegionPlot defined by the region of the plan in which precisely $n-1$ elements of ineqs are False.

    brillouinZone(n_) := RegionPlot(
      Count(ineqs({x, y}), False) == n - 1, {x, -5, 5}, {y, -5, 5}, 
      PlotStyle -> ColorData(109)(n), PlotPoints -> 100)
  4. Plot the Brillouin zones.

      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})

enter image description here
A few notes on how this could be refined:

  • The equation for i = j = 0 always evaluates as True, requiring the first layer of case-checking in PerpBisIneq. The second layer of case-checking is not necessary if gamma is always negative, but I’m not sure if this is the case.

  • The setting for PlotPoints in the definition of BrillouinZone 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 point-wise. I suspect there’s a way to do it with Reduce 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 of BrillouinZone.

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| = N-m} left( bigwedge_{P_i in Q} P_i right).

The outer “or” runs over all subsets of cardinality $N-m$, 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}_{m-1}$. 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 $N-m$ subpredicates, as an intermediate step; so Mathematica takes some time to solve this.