【Python】 基于numpy求矩阵中锐角扇环形范围内最大值及其行列号

---

前言

在矩阵中给定锐角扇环形范围内，求最大矩阵值，及其在矩阵中的行列号
注：矩阵中先行号后列号，图形中先x（列表示）后y（行表示）

计算函数

输入椭圆圆心坐标、在内椭圆上的的扇形曲线端点pm1 pn1、在外椭圆上的的扇形曲线端点pm2 pn2、锐角扇环形范围内角平分线与曲线交点坐标pct1 pct2（均为整数）及矩阵M
输出椭圆内最大值及其行列号

注：几何区域由两直线线段与两个同心椭圆曲线段构成，椭圆曲线部分的验证可改为其他曲线，

几何条件

对于锐角扇环形，计算判定是否在内部时，首先需满足该点在内椭圆外（某x方加某y方大于0）、且在外椭圆内（某x方加某y方小于0）；然后使用中心线交点pct判定，计算点需与pct在om线同侧、且与pct在on线同侧。由这些筛选该点是否位于几何图形内部。

代码

求最值的代码如下（示例）：

import numpy as np
def searchin4(po, pm1, pn1, pct1, pm2, pn2, pct2, 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 pm2: int 点M2，外椭圆上一点，包含两个数字（整数），第一个x值，第二个y值
    :param pn2: int 点N2，外椭圆上一点，包含两个数字（整数），第一个x值，第二个y值
    :param pct2: int 中心点2，区域内角平分线与椭圆的交点，包含两个数字（整数），第一个x值，第二个y值
    :param M:  ndarray 矩阵
    :return:   height：扇环形内矩阵最大值； highest_point: 最大值在M的x值y值

    注： Dmaxcoord[0]-行号-y-ymax; Dmaxcoord[1]-列号-x-xmax
        扇环形弧可能非取自圆形
    """

    rowmin = np.min(([pm1[1], pn1[1], pct1[1], pm2[1], pn2[1], pct2[1]]))
    # print(rowmin)
    rowmax = np.max(([pm1[1], pn1[1], pct1[1], pm2[1], pn2[1], pct2[1]]))
    colmin = np.min(([pm1[0], pn1[0], pct1[0], pm2[0], pn2[0], pct2[0]]))
    colmax = np.max(([pm1[0], pn1[0], pct1[0], pm2[0], pn2[0], pct2[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])
        # height = 1
        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 veri3 >= 0 and veri4 < 0 and veri1 > 0 and veri2 > 0:
            # if veri3 < 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 veri3 >= 0 and veri4 < 0 and veri1 > 0 and veri2 > 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
 

---

欢迎在评论区讨论或补充说明。