使用 numpy (Python 3) 生成软圆形遮罩

Question

我需要生成代表人眼感受野的圆形蒙版。主要问题是我不确定如何计算部分在圆外的像素覆盖的区域百分比。例如，考虑以下用于生成圆形掩码的简单代码：

def circle(size=19, r=7, x_offset=0, y_offset=0):
    x0 = y0 = size // 2
    x0 += x_offset
    y0 += y_offset
    y, x = np.ogrid[:size, :size]
    y = y[::-1]

    return ((x - x0)**2 + (y-y0)**2)<= r**2

生成这个圈子：

如您所见，输出是二进制的，如果圆圈覆盖了至少 50% 的区域，则其权重设置为 1，如果覆盖率低于 50%，则设置为 0。这是一个叠加层，显示更清楚：

另一种生成看起来更平滑的圆的常用方法是使用 2D 高斯：

def gaussian2d(size=19, sigma = 3, amplitude=1, x_offset=0, y_offset=0):
    x = np.arange(0, size, 1, float)
    y = x[:,np.newaxis][::-1]
    x0 = y0 = size // 2
    x0 += x_offset
    y0 += y_offset

    return amplitude*np.exp(-(((x-x0)**2 /(2*sigma**2)) + ((y-y0)**2 /(2*sigma**2))))

产生如下所示的输出：

这解决了二进制输出的问题（权重设置为 1 或 0），但问题是即使对于 100% 被圆圈覆盖的像素，权重也会开始下降，这是不正确的。

理想情况下，面具应该是这样的（这是手绘的）：

100%被圆圈覆盖的像素设置为1，部分覆盖的像素设置为圆圈覆盖的百分比面积，例如：

效率（在运行时和内存使用方面）并不是那么重要，我只需要一种方法来计算像素区域的百分比，即圆所覆盖的区域。准确度也不是很重要，1% 的误差应该没问题。

Answer 1

一个盒子和一个圆圈可以重叠也可以不重叠。如果四个角在圆圈内，那么这就是您的“面积百分比”目标的 100%。

所以第一步是计算四个角到圆心的距离。

di = sqrt((xi-xc)^2 + (yi-yc)^2)

如果所有di 都大于或等于圆的半径R，则框位于圆外。同样，如果所有角都验证di <= R，则单元格完全在圆圈内。

其余的箱子在圆圈内有一到三个角。
让我们看看这个“外面的一个角落”案例：

点是计算点P,Q的坐标。

因为您知道 A-B 是水平的，所以您可以通过使用 y 坐标获得 xQ 坐标（A 和 B 相同）

xQ= xc + sqrt(R^2-(yA-yc)^2)

其中xc,yc 是圆心的坐标。注意yQ = yA = yB

请注意，您可以获得两个值，具体取决于您使用的 sqrt 的符号。选择xA <= xQ <= xB。

类似，使用xA可以找到xP。

在这一刻你知道所有的坐标。
关注直P,Q。它下面的区域很简单，只是一些三角形。或单元格总面积减去三角形P,A,Q。

直线和圆周之间的面积是一个扇形的面积（您可以浏览互联网找到它）减去一个三角形。
您需要找到该扇区的角度。为此，我建议使用atan2 而不是atan，因为前者给出的角度在 (0, 2pi) 范围内。并且您必须使用归一化为 (0,pi) 范围的减法 atan2(f1)-atan2(f2)。

您可以绘制其他案例（例如三个角）并使用上述方法找到您需要的区域。

例如，您有一个带有两个角的盒子，但不知道圆圈交叉的边。你可以通过这个来测试一个方面：

• 找一面（我们称它为 AB）。如果它是水平的，您知道它的 y= yA=yB 坐标。
• 在圆圈中找到此y 的x 坐标（两种解决方案，+x，-x）
• 测试xA < x xB。如果两种解决方案都失败，则圆圈不会越过这一边。

如果侧面是垂直的，则使用其y 坐标找到x 之一。
如果一侧的两个角都在外面，那么您可以避免对其进行测试。

一个简单的近似

非常类似于Montecarlo method

对于那些穿过圆周的单元格（一些但不是所有角都在外面），您可以将每个单元格细分为n x n 网格，并测试每个子单元格。

测试可以是子单元格的中心（例如其对角线的平均点）到圆心的距离。

如果dist > R，则分配 area=0，否则分配 area=1。
那么单元格的“分配”区域就是平均值；这是n^2个子单元格的所有分配面积的总和，除以单元格的面积（即nxn）。

在最坏的情况下，大约有n 子单元格分配错误的区域。这是n/n^2 = 1/n 错误。因此，如果您想计算误差小于 1%，则需要为单元格中的网格设置 100x100 的网格。

Answer 2

我还没有完全测试过代码，但乍一看，它似乎可以工作。我会把它贴在这里，以防将来有人需要类似的东西。几乎不可能对所有可能的情况进行记录和评论，对此感到抱歉，请按原样使用代码：

def getAngle(x,y):
    if x == 0:
        if y>0:
            return 90
        else:
            return 270
    angle = np.arctan(y/x)*180/np.pi
    if x>0:
        if y<0:
            angle+=360
    else:
        angle+=180
    return angle

def getAngleRad(x,y):
    if x == 0:
        if y>0:
            return 0.5*np.pi
        else:
            return 1.5*np.pi
    angle = np.arctan(y/x)
    if x>0:
        if y<0:
            angle+=2*np.pi
    else:
        angle+=np.pi
    return angle

def circularMask(size=3, r=1, x_offset=0, y_offset=0):
    x0 = y0 = size /2.0
    x0 += x_offset
    y0 -= y_offset
    img = np.zeros((size,size))

    for j in range(img.shape[0]):
        y = y0-j
        for i in range(img.shape[1]):
            x = i-x0
            d1 = np.sqrt(x**2 + y**2)
            d2 = np.sqrt((x+1)**2 + y**2)
            d3 = np.sqrt(x**2 + (y-1)**2)
            d4 = np.sqrt((x+1)**2 + (y-1)**2)
            
            corners_inside = 0
            for d in [d1,d2,d3,d4]:
                if d<r:
                    corners_inside+=1

            if corners_inside==4:
                img[j,i]=1

            if corners_inside==3 or corners_inside==1:
                if x>0 and y>0:
                    if corners_inside==3:
                        xp = x+1
                        yq = y
                    else:
                        xp = x
                        yq = y-1
                    xq = np.sqrt(r**2-yq**2)
                    yp = np.sqrt(r**2-xp**2)

                    traingle_area = 0.5*(xp-xq)*(yq-yp)

                    theta_p = getAngleRad(xp,yp)
                    theta_q = getAngleRad(xq,yq)
                    delta_theta = np.abs(theta_p - theta_q)

                    circle_segment_area = 0.5*(delta_theta-np.sin(delta_theta))*(r**2)

                    if corners_inside==3:
                        img[j,i]=1-traingle_area+circle_segment_area
                    else:
                        img[j,i]=traingle_area+circle_segment_area
                    
                if x<0 and y>0:
                    if corners_inside==3:
                        xp = x
                        yq = y
                    else:
                        xp = x+1
                        yq = y-1
                    xq = -np.sqrt(r**2-yq**2)
                    yp = np.sqrt(r**2-xp**2)

                    traingle_area = 0.5*(xq-xp)*(yq-yp)

                    theta_p = getAngleRad(xp,yp)
                    theta_q = getAngleRad(xq,yq)
                    delta_theta = np.abs(theta_p - theta_q)

                    circle_segment_area = 0.5*(delta_theta-np.sin(delta_theta))*(r**2)

                    
                    if corners_inside==3:
                        img[j,i]=1-traingle_area+circle_segment_area
                    else:
                        img[j,i]=traingle_area+circle_segment_area
                    
                if x<0 and y<0:
                    if corners_inside==3:
                        xp = x
                        yq = y-1
                    else:
                        xp = x+1
                        yq = y
                    xq = -np.sqrt(r**2-yq**2)
                    yp = -np.sqrt(r**2-xp**2)

                    traingle_area = 0.5*(xq-xp)*(yp-yq)

                    theta_p = getAngleRad(xp,yp)
                    theta_q = getAngleRad(xq,yq)
                    delta_theta = np.abs(theta_p - theta_q)

                    circle_segment_area = 0.5*(delta_theta-np.sin(delta_theta))*(r**2)
                    
                    if corners_inside==3:
                        img[j,i]=1-traingle_area+circle_segment_area
                    else:
                        img[j,i]=traingle_area+circle_segment_area
                    
                if x>0 and y<0:
                    if corners_inside==3:
                        xp = x+1
                        yq = y-1
                    else:
                        xp = x
                        yq = y
                    xq = np.sqrt(r**2-yq**2)
                    yp = -np.sqrt(r**2-xp**2)

                    traingle_area = 0.5*(xp-xq)*(yp-yq)

                    theta_p = getAngleRad(xp,yp)
                    theta_q = getAngleRad(xq,yq)
                    delta_theta = np.abs(theta_p - theta_q)

                    circle_segment_area = 0.5*(delta_theta-np.sin(delta_theta))*(r**2)

                    if corners_inside==3:
                        img[j,i]=1-traingle_area+circle_segment_area
                    else:
                        img[j,i]=traingle_area+circle_segment_area

            if corners_inside==2:
                angle = getAngle(x+0.5,y-0.5)

                if angle<22.5 or angle>337.5:
                    yp = y
                    yq = y-1
                    xp = np.sqrt(r**2-yp**2)
                    xq = np.sqrt(r**2-yq**2)
                    
                    trapezoid_area = (0.5*(xp+xq))-x
                    
                    theta_p = np.arctan(yp/xp)
                    theta_q = np.arctan(yq/xq)
                    delta_theta = np.abs(theta_p - theta_q)

                    circle_segment_area = 0.5*(delta_theta-np.sin(delta_theta))*(r**2)

                    img[j,i] = trapezoid_area+circle_segment_area
                    

                if 22.5<angle<67.5:
                    if d1<r:
                        yp = y
                        yq = y-1
                        xp = np.sqrt(r**2-yp**2)
                        xq = np.sqrt(r**2-yq**2)
                        trapezoid_area = (0.5*(xp+xq))-x
                    else:
                        xp = x
                        xq = x+1
                        yp = np.sqrt(r**2-xp**2)
                        yq = np.sqrt(r**2-xq**2)
                        trapezoid_area = (0.5*(yp+yq))-y+1

                    theta_p = getAngleRad(xp,yp)
                    theta_q = getAngleRad(xq,yq)
                    delta_theta = np.abs(theta_p - theta_q)

                    circle_segment_area = 0.5*(delta_theta-np.sin(delta_theta))*(r**2)

                    img[j,i] = trapezoid_area+circle_segment_area

                if 67.5<angle<112.5:
                    xp = x
                    xq = x+1
                    yp = np.sqrt(r**2-xp**2)
                    yq = np.sqrt(r**2-xq**2)
                    
                    trapezoid_area = (0.5*(yp+yq))-y+1

                    theta_p = getAngleRad(xp,yp)
                    theta_q = getAngleRad(xq,yq)
                    delta_theta = np.abs(theta_p - theta_q)

                    circle_segment_area = 0.5*(delta_theta-np.sin(delta_theta))*(r**2)
                    
                    img[j,i] = trapezoid_area+circle_segment_area

                if 112.5<angle<157.5:
                    if d3< r:
                        xp = x+1
                        xq = x
                        yp = np.sqrt(r**2-xp**2)
                        yq = np.sqrt(r**2-xq**2)
                        trapezoid_area = (0.5*(yp+yq))-y+1
                    else:
                        yp = y
                        yq = y-1
                        xp = -np.sqrt(r**2-yp**2)
                        xq = -np.sqrt(r**2-yq**2)
                        trapezoid_area = x+1-(0.5*(xp+xq))

                    theta_p = getAngleRad(xp,yp)
                    theta_q = getAngleRad(xq,yq)
                    delta_theta = np.abs(theta_p - theta_q)

                    circle_segment_area = 0.5*(delta_theta-np.sin(delta_theta))*(r**2)
                    
                    img[j,i] = trapezoid_area+circle_segment_area

                if 157.5<angle<202.5:
                    yp = y
                    yq = y-1
                    xp = -np.sqrt(r**2-yp**2)
                    xq = -np.sqrt(r**2-yq**2)
                    
                    trapezoid_area = x+1-(0.5*(xp+xq))

                    theta_p = getAngleRad(xp,yp)
                    theta_q = getAngleRad(xq,yq)
                    delta_theta = np.abs(theta_p - theta_q)

                    circle_segment_area = 0.5*(delta_theta-np.sin(delta_theta))*(r**2)
                    
                    img[j,i] = trapezoid_area+circle_segment_area

                if 202.5<angle<247.5:
                    if d4<r:
                        yp = y
                        yq = y-1
                        xp = -np.sqrt(r**2-yp**2)
                        xq = -np.sqrt(r**2-yq**2)
                        trapezoid_area = x+1-(0.5*(xp+xq))
                    else:
                        xp = x
                        xq = x+1
                        yp = -np.sqrt(r**2-xp**2)
                        yq = -np.sqrt(r**2-xq**2)
                        trapezoid_area = y-(0.5*(yp+yq))

                    theta_p = getAngleRad(xp,yp)
                    theta_q = getAngleRad(xq,yq)
                    delta_theta = np.abs(theta_p - theta_q)

                    circle_segment_area = 0.5*(delta_theta-np.sin(delta_theta))*(r**2)
                    
                    img[j,i] = trapezoid_area+circle_segment_area

                if 247.5<angle<292.5:
                    xp = x
                    xq = x+1
                    yp = -np.sqrt(r**2-xp**2)
                    yq = -np.sqrt(r**2-xq**2)
                    
                    trapezoid_area = y-(0.5*(yp+yq))

                    theta_p = getAngleRad(xp,yp)
                    theta_q = getAngleRad(xq,yq)
                    delta_theta = np.abs(theta_p - theta_q)

                    circle_segment_area = 0.5*(delta_theta-np.sin(delta_theta))*(r**2)
                    
                    img[j,i] = trapezoid_area+circle_segment_area

                if 292<angle<337.5:
                    if d3<r:
                        yp = y
                        yq = y-1
                        xp = np.sqrt(r**2-yp**2)
                        xq = np.sqrt(r**2-yq**2)
                        trapezoid_area = (0.5*(xp+xq))-x
                    else:
                        xq = x
                        xp = x+1
                        yp = -np.sqrt(r**2-xp**2)
                        yq = -np.sqrt(r**2-xq**2)
                        trapezoid_area = y-(0.5*(yp+yq))

                    theta_p = getAngleRad(xp,yp)
                    theta_q = getAngleRad(xq,yq)
                    delta_theta = np.abs(theta_p - theta_q)

                    circle_segment_area = 0.5*(delta_theta-np.sin(delta_theta))*(r**2)
                    
                    img[j,i] = trapezoid_area+circle_segment_area

    return img

输出与问题中发布的值的预期结果相匹配：

编辑：我发现了一些错误并解决了它们（代码已更新），但请注意，我可能还没有发现更多错误。使用代码需要您自担风险。