【问题标题】:Changing the values in a numpy array depending on the area of a circle根据圆的面积更改 numpy 数组中的值
【发布时间】:2013-03-04 03:51:47
【问题描述】:

上下文

我需要在 python 中测量组合圆的面积。我想出了使用 numpy 数组的方法。首先,我用零填充一个网格(numpy 数组),网格中的每个位置对应于 0.5 厘米的长度。然后我将圆心放到网格上,并在网格中将此值更改为 1。我知道圆的半径,所以我可以计算圆的面积,因为我知道圆的面积,我将网格中的零更改为落在圆的区域内。然后我计算网格中的频率并用它来计算组合圆的面积,因为我知道网格中每个位置的长度我可以计算面积。这是目前一种非常粗粒度的方法,我计划在确定算法后将其更改为更细粒度的方法。

示例

如果您查看我在下面发布的图片,可以更好地描述我的想法。我的网格上有两个圆圈(红线),圆圈的中心标有蓝色方块,圆圈占据的区域为浅橙色。我想将标有橙色的区域更改为一个。我目前可以将橙色方块水平和垂直更改为圆心,但中心的对角线框给我带来了麻烦。

当前代码

class area():

    def make_grid(self):
        '''
        Each square in the grid represents 0.5 cm
        '''
        import numpy as np
        grid = np.zeros((10,10))
        square_length = 0.5
        circles = {'c1':[[4,2],1.5],'c2':[[5,6],2.0]}

        print grid
        for key,val in circles.iteritems():
            grid[val[0][0]][val[0][1]] = 1
            area = int((val[1] - square_length)/0.5)            
            for i in xrange(1,area+1):
                grid[val[0][0]][val[0][1]+i] = 1 # Change column vals in +ve direction
                grid[val[0][0]][val[0][1]-i] = 1 # Chnage column vals in -ve direction
                grid[val[0][0]+i][val[0][1]] = 1 # Chnage row vals in +ve direction 
                grid[val[0][0]-i][val[0][1]] = 1 # Chnage row vals in -ve direction

        print ''
        print grid

在上面的字典中,key是圆的名字,value的第一个元素是圆心坐标,第二个元素是圆的半径。

代码输出:

[[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  1.  0.  0.  0.  1.  0.  0.  0.]
 [ 0.  0.  1.  0.  0.  0.  1.  0.  0.  0.]
 [ 1.  1.  1.  1.  1.  0.  1.  0.  0.  0.]
 [ 0.  0.  1.  1.  1.  1.  1.  1.  1.  1.]
 [ 0.  0.  1.  0.  0.  0.  1.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  1.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  1.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]]

【问题讨论】:

  • 您不应该逐步遍历整个网格,计算每个网格点到每个圆心的距离,然后查看它是否在圆内并将网格点计数器加一吗?这样,一些点将有 2 而不是 1,但至少圆内的所有点(舍入问题除外)都将被计算在内。它比较慢,因为您必须逐步完成所有要点,但它可能更简单。
  • 你考虑几何解“测量面积”吗? mathworld.wolfram.com/Circle-CircleIntersection.html
  • 我已经考虑过几何解决方案,但是如果我有多个带有“缺失区域”的圆圈,如本问题所述:stackoverflow.com/questions/1667310/…

标签: python algorithm numpy computational-geometry


【解决方案1】:

根据@Jaime 的 cmets 更新。

我可能会从这样的事情开始。重点是圆内像素的正确计算。

import numpy as np
import matplotlib.pyplot as plt

grid = np.zeros((10,10), dtype=np.bool)
square_length = 0.5
circles = {'c1':[[4,2],1.5],'c2':[[5,6],2.0]}

# Generate arrays of indices/coordiates so we can do the
# calculations the Numpy way, without resorting to loops
# I always get the index order wrong so double check...
xx = np.arange(grid.shape[0])
yy = np.arange(grid.shape[1])

for val in circles.itervalues():
    radius = val[1]
    # same index caveat here
    # Calling Mr Pythagoras: Find the pixels that lie inside this circle
    inside = (xx[:,None] - val[0][0]) ** 2 + (yy[None, :] - val[0][1]) ** 2 <= (radius ** 2)
    # do grid & inside and initialize grid with ones for intersection instead of union
    grid = grid | inside 

plt.imshow(grid)
plt.show()

如果您在交叉点之后,请注意您的圆心相距sqrt(17) ~= 4.123.. 个单位,两个半径之和为3.5,因此实际上没有重叠。

【讨论】:

  • +1 我认为这是要走的路。我会摆脱np.sqrt 并与radius**2 进行比较。此外,通过使用&amp;,我认为您正在确定由 all 圆圈重叠的点。使用 gridzeros| 代替在 any 圈内获取点。您也可以不使用meshgrid,而是使用xx = np.arange(grid.shape[0])yy = np.arange(grid.shape[1]),然后使用广播作为inside = (xx[:, None] - val[0][0])**2 + (yy - val[0][1])**2 &lt;= radius**2 来减少内存占用。
  • 啊,是的,我理解“组合圆圈”是指交叉点,但拥有联合是有意义的。感谢 cmets - 将更新我的答案。
【解决方案2】:

您只是从中心左右直接搜索。要从中心获得对角线正方形,您可能必须使用Pythagorean theorem。使用该定理,您可以通过使用相对于圆心的水平和垂直偏移作为边长来找到正方形的斜边。然后,您可以将斜边与半径进行比较,如果斜边是两者中较短的,则将平方值加一。

半径的使用也有点奇怪,因为它不是以正方形的中间为中心。这使得半径为 2 的圆的直径为 3.5。

【讨论】:

    【解决方案3】:

    您可以使用蒙特卡罗方法找到重叠区域的面积。这样会更准确。 http://www.chem.unl.edu/zeng/joy/mclab/mcintro.html

    找到正方形的边界,然后用随机数填充,并为每个随机数检查随机数是否在圆内。

    if (np.sqrt((xx - val[0][0]) ** 2 + (yy - val[0][1]) ** 2) <= radius) :inside=inside+1
    

    面积=圆内的像素总数/生成的像素总数*正方形的面积

    【讨论】:

      【解决方案4】:

      您可以使用Floodfill algorithm,稍微修改停止条件:同时考虑颜色和到圆心的距离

      附:小区域的破解方法:绘制实心圆并计算颜色像素...

      【讨论】:

        猜你喜欢
        • 1970-01-01
        • 2019-03-04
        • 1970-01-01
        • 2019-07-21
        • 1970-01-01
        • 2021-12-10
        • 1970-01-01
        • 2022-01-07
        • 2020-02-29
        相关资源
        最近更新 更多