matplotlib 使用晶格点覆盖二维打印区域

pkbketx9  于 2023-03-19  发布在  其他
关注(0)|答案(1)|浏览(157)

我的目标是用晶格点覆盖打印区域。
在这个例子中,我们在2D中工作,我们称集合Λ R^2为格,如果存在一个基B R^2,Λ = Λ(B),集合Λ(B)是基向量的所有整数线性组合的集合,所以Λ(B)= {xb1 + yb2|x,y整数}。
换句话说,我们通过计算给定基向量的整数组合得到一个网格,对于标准基E=[[1,0]^T,[0,1]^T],我们可以把点[8,4]^T写为8 * [1,0]^T + 4 * [0,1]^T,其中8和4都是整数,所有整数组合的集合(因此格Λ)看起来像这样:

如果我们改变基,就会得到不同的格,下面是b1=[2,3]^T,b2=[4,5]^T的例子:

为了生成这些图像,我使用了以下Python代码:

import numpy as np
import matplotlib.pyplot as plt
from typing import List, Tuple

def plotLattice(ax: plt.Axes, basis_vectors: List[np.ndarray], 
                   ldown: np.ndarray, rup: np.ndarray, color: str, 
                   linewidth: float, alpha: float) -> List[np.ndarray]:
    """
    Draws a two-dimensional lattice.

    Args:
        ax: The Matplotlib Axes instance to plot on.
        basis_vectors: A list of two NumPy arrays representing the basis vectors of the lattice.
        ldown: A NumPy array representing the lower left corner of the rectangular area to draw the lattice in.
        rup: A NumPy array representing the upper right corner of the rectangular area to draw the lattice in.
        color: A string representing the color of the lattice points and basis vectors.
        linewidth: A float representing the linewidth of the lattice points.
        alpha: A float representing the alpha value of the lattice points.

    Returns:
        A list of NumPy arrays representing the lattice points.
    """
    
    # get the basis vectors
    b1, b2 = np.array(basis_vectors[0]), np.array(basis_vectors[1])
    # list to hold the lattice points
    points = []
    
    # upper bounds for the for loops
    xmax, ymax = 0, 0
    if b1[0] == 0:
        xmax = np.floor(rup[0] / abs(b2[0]))
    elif b2[0] == 0:
        xmax = np.floor(rup[0] / abs(b1[0]))
    else:
        xmax = np.floor(rup[0] / min(abs(b1[0]),abs(b2[0])))
    
    if b1[1] == 0:
        ymax = np.floor(rup[1] / abs(b2[1]))
    elif b2[1] == 0:
        ymax = np.floor(rup[1] / abs(b1[1]))
    else:
        ymax = np.floor(rup[1] / min(abs(b1[1]),abs(b2[1])))
    
    # increase the bounds by 1
    xmax = int(xmax) + 1
    ymax = int(ymax) + 1
    
    # get the lower bounds for the for loops
    xmin, ymin = -int(xmax), -int(ymax)
    
    for i in range(xmin, int(xmax)):
        for j in range(ymin, int(ymax)):
            # make the linear combination
            p = i * b1 + j * b2
            # if the point is within the plotting area, plot it and add the point to the list
            if ldown[0] <= p[0] <= rup[0] and ldown[1] <= p[1] <= rup[1]:
                ax.scatter(p[0], p[1], color=color, linewidths=linewidth, alpha=alpha)
                points.append(p)

    # plot basis vectors
    ax.quiver(0, 0, b1[0], b1[1], color=color, scale_units='xy', scale=1, alpha=1)
    ax.quiver(0, 0, b2[0], b2[1], color=color, scale_units='xy', scale=1, alpha=1)

    return points

if __name__ == '__main__':
    # pick a basis
    b1 = np.array([2, 3])
    b2 = np.array([-4, 5])
    basis_vectors = [b1, b2]
    
    # define the plotting area
    ldown = np.array([-15, -15])
    rup = np.array([15, 15])

    fig, ax = plt.subplots()
    points = plotLattice(ax, basis_vectors, ldown, rup, 'blue', 3, 0.25)

    # resize the plotting window
    mngr = plt.get_current_fig_manager()
    mngr.resize(960, 1080)
    
    # tune axis
    ax.set_aspect('equal')
    ax.grid(True, which='both')
    ax.set_xlim(ldown[0], rup[0])
    ax.set_ylim(ldown[1], rup[1])
    
    # show the plot
    plt.show()

现在我们来看看问题,对于基向量b1=[1,1],b2=[1,2],代码没有覆盖绘图区域:

我们可以通过选择一些不太正交的向量来增加问题:

因此,每当向量彼此越来越接近时,问题就会出现,因此当点积很大时,现在考虑我之前选择的例子:

我的方法是取绝对坐标值的最小值,然后估计一个轴沿着可以有多少点,这也是代码中的方法,取b1=[1,1]b2=[1,2]的x坐标的最小值,得到1,取y坐标的最小值,得到1。我的绘图区域是由ldown=[-15,-15]rup=[15,15]两个点组成的正方形定义的,因此我知道x轴沿着可以有最多floor(rup[0]/1) = 15个点。以及y轴沿着的最大值floor(rup[1]/1) = 15。包括零点在内,每个轴上有31个点,所以我希望在图上看到31*31 = 961个点。因此,我认为我已经完成并设置了xmax=15, xmin=-15, ymax=15, ymin=-15
但这给了我上面给出的结果,所以计算是错误的,然后我说,“好的,我知道点[15,-15]必须在图中”。因此我可以解出系统Bx = [15,-15]^T。这会得到向量x=[45,-30]。现在我可以设置xmax=45, ymin=-30。对点[-15,15]做同样的操作会得到向量x=[-45,30]。因此,我可以设置xmin=-45, ymin=-30,得到的图为:

这看起来很不错,但是你可以注意到图中缺少了[15,-15][-15,15]点,因此我必须通过设置xmax=46, xmin=-46, ymax=31, ymin=-31将边界扩大1,然后覆盖整个区域。
所以,这个机制的缺点是我做了一点手脚。在这里,我只知道点[15,-15]必须在绘图上。我可以解方程组并确定for循环的边界。偶尔,这个点也是离原点最远的点,所以我知道覆盖它应该自动覆盖整个绘图平面。然而,有一些基向量我无法确定这个点,我们可以从前面的图中选择一个:

在这里,我的方法是,我们可以沿着x轴有min(2,4) = 2点,沿y轴有min(3,5)=3点,但我不能简单地说[14,-9]=[7*2,-3*3]点在图上(因为它不是)。而且,我甚至不能说哪些点[12,-12], [12,-15], [14,-12],[14-15]是情节的一部分,而哪些不是,知道了这个图,我看到[12,-15][14,-12]一定在图中,不知道的话,我甚至不知道要解出Bx=b系统的哪个点。
选择不同的基或不同的(非原点中心)绘图区域使问题对我来说变得异常复杂,-即使我们只是在2D平面上操作。
所以,现在当描述这个问题时,我可以把它表述如下:给定绘图区域的点rupldown,基b1, b2,定义for回路的边界xmax, xmin, ymax, ymin,使得整个绘图区域被网格点覆盖。
我现在甚至不要求代码高效,但是xmax = sys.maxsizexmax = 100 * rup[0]类型的解决方案不算在内。
你会怎么做?

svujldwt

svujldwt1#

下面是一个解决方案的概要。
你的问题是你用来构建网格的两个向量不一定是网格中最近的点。所以我们想找到最近的点。知道这些之后,我们就可以构建网格了,方法是取网格上的每个点,然后加上它最近的点,再加上它们最近的点,依此类推,直到我们覆盖了整个绘图区域。
那么我们如何找到网格中最近的点呢?
我们要做的是不断地向队列中添加点,优先考虑那些靠近原点的点,直到每个点都在一个极小点列表中,或者是两个更近的点的和或差。
为了区分优先级,我们可以使用priority queue
现在让我们以(4, 1)(4, 2)为起始向量,看看它是如何工作的,我们将有一个要查找的点的队列upcoming,为了清楚起见,我将按排序顺序来写它。
我们会有一份检查点的清单。
我们会得到一组非极小点。
我们从以下内容开始:

upcoming = [Point(4, 1), Point(4, 2)]
examined = []
not_minimal = set([])

现在我们取出第一个点,将其加到examined,将其与所有检查点的和或差加到upcomingnot_minimal,这取决于它离我们加的点更近还是更远。
现在我们的代码是这样的(未经测试)

while 0 < len(upcoming):
    point = heapq.heappop(upcoming)
    if point not in not_minimal: # Might have found another way to get it.
        examined.append(point)
        for point2 in examined:
            for point3 in [
                point + point2,
                point - point2,
                -point + point2,
                -point - point2
            ]):
                if max(distance(point), distance(point2)) < distance(point3):
                    not_minimal.add(point3)
                else:
                    heapq.heappush(upcoming, point3)
minimal_points = [point for point in examined if point not in not_minimal]

现在你将得到离原点最近的点。现在从你最终答案中想要的一个点开始,并围绕它建立网格。根据网格的不同,你可能需要离开绘图区域max((distance(point) for point in minimal_points))才能发现网格中的一些角点。但你不应该错过任何一个。
有限点的最初发现需要有限时间,其上限为时间O(m log(m)),其中m是在固定圆/球/任何围绕半径为较大向量两倍的原点的向量中有多少个向量。
发现整个网格需要时间O(n * k),其中k是最小点的数量,n是网格中需要查找的点的数量。
这个算法可以在任意维数上工作。
为了好玩,我把它编了代码。这只是返回你想要的点,而不是画它们。修复它很容易。它也将处理超过二维。画它将需要更多的工作。

import numpy as np
import heapq

def distance_to_box (ldown, lup, point):
    deltas = []
    for i in range(len(ldown)):
        (x_min, x_max) = sorted([ldown[i], lup[i]])
        if point[i] < x_min:
            deltas.append(x_min - point[i])
        elif x_max < point[i]:
            deltas.append(point[i] - x_max)
        else:
            deltas.append(0)
    da = np.array(deltas)
    return np.linalg.norm(da)

def grid (basis_vectors, ldown, lup, epsilon = 0.1**8):
    upcoming = []
    i = 0
    for v in basis_vectors:
        point = np.array(v)
        i += 1
        heapq.heappush(upcoming, (np.dot(point, point), i, point))

    minimal = {}
    not_minimal = {}
    while len(upcoming):
        (d, _, p) = heapq.heappop(upcoming)
        if str(p) not in not_minimal:
            minimal[str(p)] = p
            to_delete = []
            for s2, p2 in minimal.items():
                if s2 in not_minimal:
                    to_delete.append(s2)
                else:
                    d2 = np.linalg.norm(p2)
                    for p3 in [p + p2, p - p2, -p + p2, -p - p2]:
                        d3 = np.linalg.norm(p3)
                        if max(d, d2) + epsilon < d3:
                            not_minimal[str(p3)] = p3
                        elif d3 < epsilon:
                            pass # ignore 0 vector
                        else:
                            i += 1
                            heapq.heappush(upcoming, (d3, i, p3))
            for s2 in to_delete:
                minimal.pop(s2)
    directions = minimal.values()

    start = np.array([0 for _ in ldown])
    improved = True
    while improved:
        improved = False
        for direction in directions:
            step = direction
            # Does walking in this direction help?
            while np.linalg.norm(start + step - ldown) < np.linalg.norm(start - ldown):
                improved = True
                start = start + step
                step = step + step

    todo = [start]
    seen = {}
    max_d = max([np.linalg.norm(d) for d in directions]) + epsilon
    answer = []
    while len(todo):
        p = todo.pop()
        d = distance_to_box(ldown, lup, p)
        if d < max_d:
            if str(p) not in seen:
                for direction in directions:
                    todo.append(p + direction)
                if d < epsilon:
                    answer.append(p)
                seen[str(p)] = p
    return answer


for p in grid([(4, 1), (4, 2)], (-4, -4), (4, 4)):
    print(p)

相关问题