时间戳

首页 > 解决方案 > 如何在python中将二维插值函数显示为矩阵？

python - 如何在python中将二维插值函数显示为矩阵？

问题描述

我环顾四周，但很难找到答案。基本上，当一个插值 v -> w 时，您通常会使用许多插值函数之一。但我想得到对应的矩阵Av = w。

在我的情况下，w 是一个 200x200 矩阵，其中 v 是 w 的随机子集，点数是 w 的一半。我真的不关心花哨的数学，它可以像用距离平方加权已知点一样简单。我已经尝试用一些 for 循环来实现它，但它只适用于小值。但也许它有助于解释我的问题。

from random import sample

    def testScatter(xbig, ybig):
        NumberOfPoints = int(xbig * ybig / 2) #half as many points as in full Sample

    #choose random coordinates
    Index = sample(range(xbig * ybig),NumberOfPoints)
    IndexYScatter = np.remainder(Index,  xbig)
    IndexXScatter = np.array((Index - IndexYScatter) / xbig, dtype=int)

    InterpolationMatrix = np.zeros((xbig * ybig , NumberOfPoints), dtype=np.float32)
    WeightingSum = np.zeros(xbig * ybig )
    coordsSamplePoints = []
    for i in range(NumberOfPoints): #first set all the given points (no need to interpolate)
        coordsSamplePoints.append(IndexYScatter[i] + xbig * IndexXScatter[i])
        InterpolationMatrix[coordsSamplePoints[i], i] = 1
        WeightingSum[coordsSamplePoints[i]] = 1
    
    for x in range(xbig * ybig): #now comes the interpolation
        if x not in coordsSamplePoints:
            YIndexInterpol = x % xbig      #xcoord in interpolated matrix
            XIndexInterpol = (x - YIndexInterpol) / xbig  #ycoord in interp. matrix
            for y in range(NumberOfPoints):
                XIndexScatter = IndexXScatter[y]
                YIndexScatter = IndexYScatter[y]
                distanceSquared = (np.float32(YIndexInterpol) - np.float32(YIndexScatter))**2+(np.float32(XIndexInterpol) - np.float32(XIndexScatter))**2
                InterpolationMatrix[x,y] = 1/distanceSquared
                WeightingSum[x] += InterpolationMatrix[x,y]

    return InterpolationMatrix/ WeightingSum[:,None] , IndexXScatter, IndexYScatter

标签： pythonmatrix2dinterpolation

解决方案

您需要花一些时间阅读从本页顶部开始的 Numpy 文档，然后逐步进行。在 SO 上研究答案，询问如何在使用 Numpy 数组时对操作进行矢量化处理会对您有所帮助。如果您发现您正在迭代索引并使用 Numpy 数组执行计算，那么可能有更好的方法。

第一次切割...
第一个 for 循环可以替换为：

coordsSamplePoints = IndexYScatter + (xbig * IndexXScatter)
InterpolationMatrix[coordsSamplePoints,np.arange(coordsSamplePoints.shape[0])] = 1
WeightingSum[coordsSamplePoints] = 1

这主要利用元素算术和索引数组- 应该阅读完整的索引教程

您可以通过增强函数并执行 for 循环以及Numpy 方式来测试这一点，然后比较结果。

...
    IM = InterpolationMatrix.copy()
    WS = WeightingSum.copy()
    for i in range(NumberOfPoints): #first set all the given points (no need to interpolate)
        coordsSamplePoints.append(IndexYScatter[i] + xbig * IndexXScatter[i])
        InterpolationMatrix[coordsSamplePoints[i], i] = 1
        WeightingSum[coordsSamplePoints[i]] = 1
    cSS = IndexYScatter + (xbig * IndexXScatter)
    IM[cSS,np.arange(cSS.shape[0])] = 1
    WS[cSS] = 1
    # TEST Validity
    print((cSS == coordsSamplePoints).all(),
          (IM == InterpolationMatrix).all(),
          (WS == WeightingSum).all())
...

外循环：

...
    for x in range(xbig * ybig): #now comes the interpolation
        if x not in coordsSamplePoints:
            YIndexInterpol = x % xbig      #xcoord in interpolated matrix
            XIndexInterpol = (x - YIndexInterpol) / xbig  #ycoord in interp. matrix
            ...

可以替换为：

...
    space = np.arange(xbig * ybig)
    mask = ~(space == cSS[:,None]).any(0)
    iP = space[mask]    # points to interpolate
    yIndices = iP % xbig
    xIndices = (iP - yIndices) / xbig
    ...

完整的解决方案：

import random
import numpy as np
def testScatter(xbig, ybig):
    NumberOfPoints = int(xbig * ybig / 2) #half as many points as in full Sample

    #choose random coordinates
    Index = random.sample(range(xbig * ybig),NumberOfPoints)
    IndexYScatter = np.remainder(Index,  xbig)
    IndexXScatter = np.array((Index - IndexYScatter) / xbig, dtype=int)
    InterpolationMatrix = np.zeros((xbig * ybig , NumberOfPoints), dtype=np.float32)
    WeightingSum = np.zeros(xbig * ybig )
    
    coordsSamplePoints = IndexYScatter + (xbig * IndexXScatter)
    InterpolationMatrix[coordsSamplePoints,np.arange(coordsSamplePoints.shape[0])] = 1
    WeightingSum[coordsSamplePoints] = 1

    IM = InterpolationMatrix
    cSS = coordsSamplePoints
    WS = WeightingSum

    space = np.arange(xbig * ybig)
    mask = ~(space == cSS[:,None]).any(0)
    iP = space[mask]    # points to interpolate
    yIndices = iP % xbig 
    xIndices = (iP - yIndices) / xbig
    dSquared = ((yIndices[:,None] - IndexYScatter) ** 2) + ((xIndices[:,None] - IndexXScatter) ** 2)
    IM[iP,:] = 1/dSquared
    WS[iP] = IM[iP,:].sum(1)

    return IM / WS[:,None], IndexXScatter, IndexYScatter

与使用 (100,100) 的参数相比，我得到了大约 200 倍的改进。可能还有其他一些小的改进，但它们不会显着影响执行时间。

广播是另一个必须的 Numpy 技能。

推荐阅读