前言
在矩阵中给定钝角扇环形范围内，求最大矩阵值，及其在矩阵中的行列号
注：矩阵中先行号后列号，图形中先x（列表示）后y（行表示）
计算函数
输入椭圆圆心坐标、在内外椭圆上的的扇形曲线端点pm1 pn1 pm2 pn2、钝角扇形范围内角平分线与内外曲线交点坐标pct1 pct2（均为整数）及矩阵M
输出椭圆内最大值及其行列号
注：几何区域由两直线线段与椭圆曲线段构成，椭圆曲线部分的验证可改为其他曲线，
几何条件
对于钝角扇形，计算判定是否在内部时，首先需满足该点在内椭圆外（某x方加某y方大于0）、该点在外椭圆内（某x方加某y方大于0）；然后使用中心线交点pct判定，该计算点只要不是同时满足两个条件（1. 与pct在om线不同侧 2. 与pct在on线不同侧）则均在所选区域内。
由此可判定该点是否位于几何图形内部。
代码
求最值的代码如下（示例）：
import numpy as np
def searchin6(po, pm1, pn1, pct1, pa1, pb1, pc1, pd1, pm2, pn2, pct2, pa2, pb2, pc2, pd2, M):
"""
求钝角扇环形内矩阵最大值，及最大值在M的xy值
这部分示例代码还输入了椭圆上下左右点，但计算中非必需
:param po: int 点O，椭圆心，包含两个数字（整数），第一个x值，第二个y值
:param pm1: int 点M1，内椭圆上一点，包含两个数字（整数），第一个x值，第二个y值
:param pn1: int 点N1，内椭圆上一点，包含两个数字（整数），第一个x值，第二个y值
:param pct1: int 中心点1，区域内角平分线与椭圆的交点，包含两个数字（整数），第一个x值，第二个y值
:param pa1: int 点A1，内椭圆上最右点，包含两个数字（整数），第一个x值，第二个y值
:param pb1: int 点B1，内椭圆上最上点，包含两个数字（整数），第一个x值，第二个y值
:param pc1: int 点C1，内椭圆上最左点，包含两个数字（整数），第一个x值，第二个y值
:param pd1: int 点D1，内椭圆上最下点，包含两个数字（整数），第一个x值，第二个y值
:param pm2: int 点M2，外椭圆上一点，包含两个数字（整数），第一个x值，第二个y值
:param pn2: int 点N2，外椭圆上一点，包含两个数字（整数），第一个x值，第二个y值
:param pct2: int 中心点2，区域内角平分线与椭圆的交点，包含两个数字（整数），第一个x值，第二个y值
:param pa2: int 点A2，外椭圆上最右点，包含两个数字（整数），第一个x值，第二个y值
:param pb2: int 点B2，外椭圆上最上点，包含两个数字（整数），第一个x值，第二个y值
:param pc2: int 点C2，外椭圆上最左点，包含两个数字（整数），第一个x值，第二个y值
:param pd2: int 点D2，外椭圆上最下点，包含两个数字（整数），第一个x值，第二个y值
:param M: ndarray 矩阵
:return: height：扇环形内矩阵最大值； highest_point: 最大值在M的x值y值
注： Dmaxcoord[0]-行号-y-ymax; Dmaxcoord[1]-列号-x-xmax
扇环形弧可能非取自圆形
"""
rowmin = np.min(([pa2[1], pb2[1], pc2[1], pd2[1]]))
# print(rowmin)
rowmax = np.max(([pa2[1], pb2[1], pc2[1], pd2[1]]))
colmin = np.min(([pa2[0], pb2[0], pc2[0], pd2[0]]))
colmax = np.max(([pa2[0], pb2[0], pc2[0], pd2[0]]))
x0=po[0]
y0=po[1]
D = M[rowmin:rowmax, colmin:colmax]
vm1 = (pm1[0]-po[0], pm1[1]-po[1]) # 向量
vn1 = (pn1[0]-po[0], pn1[1]-po[1])
vct1 = (pct1[0]-po[0], pct1[1]-po[1])
vm2 = (pm2[0] - po[0], pm2[1] - po[1]) # 向量
vn2 = (pn2[0] - po[0], pn2[1] - po[1])
vct2 = (pct2[0] - po[0], pct2[1] - po[1])
# 算椭圆1
if abs(pm1[0] - pn1[0]) > 0.001 and abs(pm1[1] - pn1[1]) > 0.001:
a11 = (vn1[0] ** 2) * (vm1[1] ** 2) - (vm1[0] ** 2) * (vn1[1] ** 2)
a21 = vm1[1] ** 2 - vn1[1] ** 2
b21 = vn1[0] ** 2 - vm1[0] ** 2
asquare1 = a11/a21 # ((pn[0] ** 2) * (pm[1] ** 2) - (pm[0] ** 2) * (pn[1] ** 2)) / (pm[1] ** 2 - pn[1] ** 2)
bsquare1 = a11/b21 # ((pn[0] ** 2) * (pm[1] ** 2) - (pm[0] ** 2) * (pn[1] ** 2)) / (pn[0] ** 2 - pm[0] ** 2)
else:
a11 = (vct1[0] ** 2) * (vm1[1] ** 2) - (vm1[0] ** 2) * (vct1[1] ** 2)
asquare1 = a11 / (vm1[1] ** 2 - vct1[1] ** 2)
bsquare1 = a11 / (vct1[0] ** 2 - vm1[0] ** 2)
ra1 = asquare1**0.5
rb1 = bsquare1**0.5
# print(ra,rb,asquare,bsquare)
# 算椭圆2
if abs(pm2[0] - pn2[0]) > 0.001 and abs(pm2[1] - pn2[1]) > 0.001:
a12 = (vn2[0] ** 2) * (vm2[1] ** 2) - (vm2[0] ** 2) * (vn2[1] ** 2)
a22 = vm2[1] ** 2 - vn2[1] ** 2
b22 = vn2[0] ** 2 - vm2[0] ** 2
asquare2 = a12 / a22 # ((pn[0] ** 2) * (pm[1] ** 2) - (pm[0] ** 2) * (pn[1] ** 2)) / (pm[1] ** 2 - pn[1] ** 2)
bsquare2 = a12 / b22 # ((pn[0] ** 2) * (pm[1] ** 2) - (pm[0] ** 2) * (pn[1] ** 2)) / (pn[0] ** 2 - pm[0] ** 2)
else:
a12 = (vct2[0] ** 2) * (vm2[1] ** 2) - (vm2[0] ** 2) * (vct2[1] ** 2)
asquare2 = a12 / (vm2[1] ** 2 - vct2[1] ** 2)
bsquare2 = a12 / (vct2[0] ** 2 - vm2[0] ** 2)
ra2 = asquare2 ** 0.5
rb2 = bsquare2 ** 0.5
count1 = 0
while True:
Dmax = np.max(D)
Dmaxcoord = np.where(D == np.max(D))
#print(Dmaxcoord,np.size(Dmaxcoord[0]))
itern = np.size(Dmaxcoord[0])
if itern > 1:
ymax = Dmaxcoord[0][0] + rowmin
xmax = Dmaxcoord[1][0] + colmin
veri1 = (vm1[0] * vct1[1] - vm1[1] * vct1[0]) * (vm1[0] * (ymax - y0) - vm1[1] * (xmax - x0)) # 算线om，是否与pct在om的同一侧
veri2 = (vn1[0] * vct1[1] - vn1[1] * vct1[0]) * (vn1[0] * (ymax - y0) - vn1[1] * (xmax - x0)) # 算线on
veri3 = bsquare1 * ((xmax - x0) ** 2) + asquare1 * ((ymax - y0) ** 2) - asquare1 * bsquare1 # 椭圆方程，需在外部
veri4 = bsquare2 * ((xmax - x0) ** 2) + asquare2 * ((ymax - y0) ** 2) - asquare2 * bsquare2 # 椭圆方程，需在内部
if veri1 < 0 and veri2 < 0:
veri12 = -1 # 在所求区域外
else:
veri12 = 1 # 在所求区域内
if veri3 >= 0 and veri4 < 0 and veri12 > 0:
height = Dmax
highest_point = (xmax + 1, ymax + 1)
return height, highest_point
else:
Dx = Dmaxcoord[0][0]
Dy = Dmaxcoord[1][0]
D[Dx, Dy] = 10
# print('错误结果 %f\n', xmax, ymax)
count1 = count1 + 1
else:
ymax = Dmaxcoord[0] + rowmin
xmax = Dmaxcoord[1] + colmin
veri1 = (vm1[0] * vct1[1] - vm1[1] * vct1[0]) * (vm1[0] * (ymax - y0) - vm1[1] * (xmax - x0)) # 算线om，是否与pct在om的同一侧
veri2 = (vn1[0] * vct1[1] - vn1[1] * vct1[0]) * (vn1[0] * (ymax - y0) - vn1[1] * (xmax - x0)) # 算线on
veri3 = bsquare1 * ((xmax - x0) ** 2) + asquare1 * ((ymax - y0) ** 2) - asquare1 * bsquare1 # 椭圆方程，需在外部
veri4 = bsquare2 * ((xmax - x0) ** 2) + asquare2 * ((ymax - y0) ** 2) - asquare2 * bsquare2 # 椭圆方程，需在内部
if veri1 < 0 and veri2 < 0:
veri12 = -1 # 在所求区域外
else:
veri12 = 1 # 在所求区域内
if veri3 >= 0 and veri4 < 0 and veri12 > 0:
height = Dmax
highest_point = (xmax + 1, ymax + 1)
return height, highest_point
else:
D[Dmaxcoord[0],Dmaxcoord[1]]=10
# print('错误结果 %f\n',xmax,ymax)
count1 = count1+1
if count1 >= ra2*rb2*4:
height = Dmax
highest_point = po
return height, highest_point
欢迎在评论区讨论或补充说明。