n阶三元曲面拟合以及代码

最新推荐文章于 2024-06-24 09:30:20 发布

爱敲代码的研究生

最新推荐文章于 2024-06-24 09:30:20 发布

阅读量614

点赞数 2

本文链接：https://blog.csdn.net/qq_41502889/article/details/126706458

版权

机器学习专栏收录该内容

1 篇文章 0 订阅

订阅专栏

import numpy as np
import matplotlib.pyplot as plt #导入库
from mpl_toolkits.mplot3d import Axes3D #导入库
from mpl_toolkits import mplot3d
import numpy as np
import matplotlib.pyplot as plt



class PointsFitting(object):
    # n为曲面拟合阶次
    def __init__(self, x_list, y_list, z_list, n):
        # print('xlist',len(x_list))
        # print('rand',np.random.randn(1,len(x_list)).shape)
        self.x_list = x_list
        self.y_list = y_list
        x_list = np.array(x_list) - (2 * np.random.randn(len(x_list)) - 1) * 0.001
        y_list = np.array(y_list) - (2 * np.random.randn(len(x_list)) - 1) * 0.001
        self.z_list = np.array(z_list) - (2 * np.random.randn(len(x_list)) - 1) * 0.001

        self.n = n
        self.x_n = self.cal_x(x_list)
        self.y_n = self.cal_y(y_list)
        self.q = self.cal_xyn(x_list, y_list)
        self.P = self.GetResults()

    def cal_xyn(self, x_list, y_list):
        all_lists = []
        length = len(x_list)
        for m in range(length):
            list = []
            for j in range(1, self.n):
                for i in range(1, self.n - j):
                    list.append(y_list[m] ** j * x_list[m] ** i)
            all_lists.append(list)
        return all_lists

    def cal_x(self, x_list):
        all_lists = []
        length = len(x_list)
        for m in range(length):
            list = []
            # print('xlist',x_list[m])
            for j in range(0, self.n):
                list.append(x_list[m] ** j)
            all_lists.append(list)
        return all_lists

    def cal_y(self, y_list):
        all_lists = []
        length = len(y_list)
        for m in range(length):
            list = []
            for j in range(1, self.n):
                list.append(y_list[m] ** j)
            all_lists.append(list)
        return all_lists

    def large_matrix(self):

        x_n = np.array(self.x_n)
        y_n = np.array(self.y_n)
        q = np.array(self.q)
        x_n_T = x_n.T
        XTX = np.dot(x_n_T, x_n)
        XTY = np.dot(x_n_T, y_n)
        XTQ = np.dot(x_n_T, q)
        y_n_T = y_n.T
        YTY = np.dot(y_n_T, y_n)
        YTQ = np.dot(y_n_T, q)
        q_n_T = q.T
        QTQ = np.dot(q_n_T, q)
        w11 = XTX
        w12 = XTY
        w13 = XTQ
        w21 = XTY.T
        w22 = YTY
        w23 = YTQ
        w31 = XTQ.T
        w32 = YTQ.T
        w33 = QTQ
        W1 = np.concatenate((w11, w12, w13), axis=1)
        W2 = np.concatenate((w21, w22, w23), axis=1)
        W3 = np.concatenate((w31, w32, w33), axis=1)
        W = np.concatenate((W1, W2, W3), axis=0)
        return W

    # 计算矩阵右边的结果
    def values(self):

        x_n = np.array(self.x_n)
        y_n = np.array(self.y_n)
        q = np.array(self.q)
        z = self.z_list
        rex = np.dot(x_n.T, z)
        rey = np.dot(y_n.T, z)
        req = np.dot(q.T, z)
        RESULTS = np.concatenate((rex, rey, req), axis=0)
        return RESULTS

    # 计算曲面参数向量
    def GetResults(self):
        W = self.large_matrix()
        # print(W.shape)
        RESULTS = self.values()
        # print(RESULTS.shape)
        # # 计算参数向量结果
        # print('w',np.linalg.matrix_rank(W))
        P = np.dot(np.linalg.inv(W), RESULTS)
        # print(P)
        return P

    # input须为二位数组,形状为n*2
    def calculate(self, input):
        input = np.array(input)
        x_n = np.array(self.cal_x(input[:, 0]))
        y_n = np.array(self.cal_y(input[:, 1]))
        xyn = self.cal_xyn(np.squeeze(input[:, 0]), np.squeeze(input[:, 1]))
        cat = np.concatenate((x_n, y_n, xyn), axis=1)
        result = np.dot(cat, self.P)
        # print(result)
        return result

    # 计算切点的平面法向量
    def CalGradient(self, x, y):
        list1 = []
        list2 = []
        sum_k = 0

        for i in range(1, self.n):
            list1.append(self.P[i] * i * x ** (i - 1))
            list2.append(self.P[i + self.n - 1] * i * y ** (i - 1))
        #     print('i',i)
        #     print('i+n',i + self.n-1)
        # print('p',len(self.P))
        # print('n',self.n)
        # print('l1',len(list1))
        # print('l2',len(list2))
        for j in range(1, self.n):
            for i in range(1, self.n - j):
                # print(2 * self.n + sum_k + i)
                list1.append(self.P[2 * self.n + sum_k + i - 2] * i * x ** (i - 1) * y ** j)
                list2.append(self.P[2 * self.n + sum_k + i - 2] * j * y ** (j - 1) * x ** i)
            sum_k = sum_k + (self.n - j - 1)
        gx = np.sum(np.array(list1))
        gy = np.sum(np.array(list2))
        gz = -1
        return gx, gy, gz

if __name__ =='__main__':
    x=-np.random.randn(500)*2+1
    y=-np.random.randn(500)*2+1
    z=x**3+y**2
    x=list(x)
    y=list(y)
    z=list(z)
    PF=PointsFitting(x,y,z,6)

    input=np.concatenate((np.expand_dims(np.array(PF.x_list), axis=1), np.expand_dims(np.array(PF.y_list), axis=1)), axis=1)

    z1=PF.calculate(input)
    # print(z1)
    fig = plt.figure()
    print('point:(9,18)')
    print('阶次',5)
    print('梯度',PF.CalGradient(9,18))
    # 创建绘图区域

    ax = plt.axes(projection='3d')
    ax.scatter3D(x, y, z1)
    ax.set_title('3d Scatter plot')
    plt.show()

主要原理应用的是多项式拟合加上最小二乘法，制作不易，请勿白嫖，欢迎点赞