当先锋百科网

首页 1 2 3 4 5 6 7

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


前言

在矩阵中给定钝角扇环形范围内,求最大矩阵值,及其在矩阵中的行列号
注:矩阵中先行号后列号,图形中先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
 

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