问题总结：给定几个点，绘制从原点到这几个点的垂直平分线，然后几个向量的垂直平分线相交围成的最小面积（三维情况下就是体积）的图形就是需要绘制的图形。如何能绘制出这个图形
背景：打算绘制一个布里渊区示意图；布里渊区图绘制方法：给定两个基矢量b1,b2；绘制所有G=n1*b1+n2*b2，其中n1,n2为整数，比如我们给定n从-10取值到10，一个晶格常数为1的简单格子的布里渊区就是{{-10pi,-10pi},{-10pi,-9pi},...,{-10pi,10pi},{-9pi,-10pi},{-9pi,-9pi},...,{-9pi,10pi},......,{10pi,10pi}}(不严谨但是不影响画图过程)。
第一布里渊区取从{0,0}点到上述格点的最短距离，也就是最近邻格点。绘制这个矢量的垂直平分线，最近邻格点的所有垂直平分线围成的最小面积就是第一布里渊区
第二布里渊区就是次近邻格点重复上述过程；
目前问题：
布里渊区的取值当然比较轻松，
Subscript[b, 1] = (2 \[Pi])/a {1, 0};
Subscript[b, 2] = (2 \[Pi])/a {0, 1};
Subscript[b, 1 t] = Table[n*Subscript[b, 1], {n, -10, 10}];
Subscript[b, 2 t] = Table[n*Subscript[b, 2], {n, -10, 10}];
Subscript[b, 1 t][[1]] + Subscript[b, 2 t][[1]];
input = Flatten[
Table[Subscript[b, 2 t][[j]] + Subscript[b, 1 t][[i]], {i, 1,
21}, {j, 1, 21}], 1];
这样可以得到所有布里渊区的格点的坐标。如果取n=10，则有441个格点。
接下来需要找到最近邻格点，我使用了For循环来找，因为如果直接求Norm，就会返回一个列表的模，我觉得这方面我不是很熟悉，我猜应该有更为简洁省力的方案。但是这个方法至少是可以找出来的。
我可以根据如上的方案找到最近邻、次近邻、次次近邻，.......的格点，绘制出他们的垂直平分线，首先我没有想明白这个垂直平分线该怎么绘制，因为我们需要让最近邻格点的垂直平分线互相相交，然后找到他们围成的最小面积。

尽力了，前两个能画，到第三个就不对了。似乎第三布里渊区不完全按这个规则画？不然应该是一个矩形。部分函数是版本12的，尤其是求中垂线的PerpendicularBisector函数。绘制思路就是采用几何元的堆砌
Clear["`*"]
kMax = 3;
a = 1;
b1 = (2 \[Pi])/a {1, 0};
b2 = (2 \[Pi])/a {0, 1};
base = Flatten[Table[i b1 + j b2, {i, -5, 5}, {j, -5, 5}], 1];
distance = DeleteDuplicates[Norm /@ base] // N;
kthNearest[k_] := With[{dist = distance}, RankedMin[dist, k + 1]]
kthNearestPoints[k_] := Cases[base, v_ /; Norm[v] == kthNearest[k]]
color[k_] := ColorData[35, "ColorList"][[k]]
findpb[vec_?VectorQ] := PerpendicularBisector[{{0, 0}, vec}]
findPoint[v1_?VectorQ, v2_?VectorQ] :=
Block[{x, y},
SolveValues[{x, y} \[Element] findpb[v1] && {x, y} \[Element]
findpb[v2], {x, y}, Reals] // Quiet]
polars[k_] :=
Cases[Flatten[
Table[findPoint[v1, v2], {v1, kthNearestPoints[k]}, {v2,
kthNearestPoints[k]}], 2], {x_?NumericQ, y_?NumericQ}] //
DeleteDuplicates
\[Theta][v_?VectorQ] :=
Piecewise[{{ArcCos[{1, 0} . v/Norm[v]],
v[[2]] >= 0}, {-ArcCos[{1, 0} . v/Norm[v]], v[[2]] < 0}}]
polygon[k_] := {Opacity[0.3], color[k],
Polygon[Sort[polars[k], \[Theta][
#1] > \[Theta][#
2] &]]}
Graphics[{
{Point[base]},
Sequence @@ Table[{color[k], Point[kthNearestPoints[k]]}, {k, kMax}],
Sequence @@ Table[polygon[k], {k, kMax}],
Inset[PointLegend[Table[color[k], {k, kMax}],
Table[ToString[k] <> "th", {k, kMax}]], {45, 1}]
}, ImageSize -> 260]

首先声明LZ在弄的这个我以前没接触过，但是查了一下，LZ这个解说应该是有问题的，当然这个似乎也不能怪LZ，因为诸如
http://www.sciencedirect.com/topics/earth-and-planetary-sciences/brillouin-zone#:~:text=Similarly%2C%20the%20second%20Brillouin%20zone,is%20the%20third%20Brillouin%20zone.
里的解说也差不多是这么套说辞，可是配的插图又不是这么回事。从
http://www.doitpoms.ac.uk/tlplib/brillouin_zones/zone_construction.php
给出的分步图来看，
正确的
绘制方法似乎是要把前几个布里渊区的垂直平分线一并保留，然后找交点，此外还附带相邻和面积相等的限制……显然这个画法也太烦了点，所以实践中大家好像都是用更简明也更严格的那个数学定义来绘制（在这里提了：physics.stackexchange.com/q/264646/18778），参：
mathematica.stackexchange.com/q/255164/1871
demonstrations.wolfram.com/2DBrillouinZones/

……难得有个有点意思的问题，再认真点吧。使用上一楼最后提及的思路，搭配一些高级的区域画图技巧（不过吧里的老人们大概也看腻了），我们可以得到如下代码：
cf =
Compile[{n, dn},
With[{gridlst = Flatten[Table[{i, j}, {i, -n, n}, {j, -n, n}], 1]},
Table[
Sum[If[Total[grid^2]/2 < {px, py} . grid, 1., 0.], {grid, gridlst}], {px, -n,
n, dn}, {py, -n, n, dn}]], RuntimeOptions -> "Speed",
CompilationTarget -> "C"]; // AbsoluteTiming
(* {0.46108, Null} *)
num = 2; dx = 0.01; dat = cf[num, dx]; // AbsoluteTiming
MatrixPlot[dat, DataRange -> {{-1, 1}, {-1, 1}} num, PlotRange -> {1/2, 5}]
(* {0.319835, Null} *)
把dx改小一点可以得到像素风布里渊区（下图为dx=0.1）：
不喜欢上面的配色的话大家可以自行修改ColorRules等选项。