使用 numpy 生成软圆形掩码 (Python 3)
Generating a soft circluar mask using numpy (Python 3)
我需要生成代表人眼感受野的圆形蒙版。主要问题是我不确定如何计算部分位于圆圈外的像素所覆盖的面积百分比。例如，考虑这个用于生成圆形遮罩的简单代码：
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% 的误差应该没问题。
一个方框和一个圆可以重叠也可以不重叠。如果四个角在圆圈内，那么这就是您的“面积百分比”目标的 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 的网格作为单元格中的网格。
我还没有完全测试代码，但乍一看，它似乎可以工作。我会 post 它在这里，以防将来有人需要类似的东西。几乎不可能对每种可能的情况进行记录和评论，对此深表歉意，请按原样使用代码：
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
输出与问题中输入的值 post 的预期结果相匹配：
编辑：我发现了一些错误并解决了它们（代码已更新），但请注意，可能还有更多我尚未发现的错误。使用代码需要您自担风险。
我需要生成代表人眼感受野的圆形蒙版。主要问题是我不确定如何计算部分位于圆圈外的像素所覆盖的面积百分比。例如，考虑这个用于生成圆形遮罩的简单代码：
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% 的误差应该没问题。
一个方框和一个圆可以重叠也可以不重叠。如果四个角在圆圈内，那么这就是您的“面积百分比”目标的 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 的网格作为单元格中的网格。
我还没有完全测试代码，但乍一看，它似乎可以工作。我会 post 它在这里，以防将来有人需要类似的东西。几乎不可能对每种可能的情况进行记录和评论，对此深表歉意，请按原样使用代码：
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
输出与问题中输入的值 post 的预期结果相匹配：
编辑：我发现了一些错误并解决了它们（代码已更新），但请注意，可能还有更多我尚未发现的错误。使用代码需要您自担风险。