将矢量化 Matlab 代码传输到 python、numpy

Question

我正在将我的 matlab 代码传输到 python。我试图在 python 和 numpy 中找到很多东西的替代品

Matlab 代码：

    [m,n]=size(Image);

    canvas=zeros(m,n);

    U_res_sel=squeeze(loading);
    [s1,s2]=size(U_res_sel);

    if mod(s1,2)==0
        dy=s1/2-1;
    else
        dy=(s1-1)/2;
    end
    if mod(s2,2)==0
        dx=s2/2-1;
    else
        dx=(s2-1)/2;
    end

    xmin=dx+1;
    ymin=dy+1;
    xmax=n-(s2-dx-1);
    ymax=m-(s1-dy-1);

    [x,y]=meshgrid(xmin:xmax,ymin:ymax);

    ind=sub2ind([m,n],y(:),x(:));

    nps=repmat(((-dx+(0:s2-1))*m-dy),s1,1)+repmat((0:(s1-1)).',1,s2);

    ind=repmat(ind,1,numel(nps))+repmat(nps(:).',numel(ind),1);

    A=(Image(ind)-repmat(mean(Image(ind),2),1,numel(nps)));

    B=repmat((U_res_sel(:)-mean(U_res_sel(:))).',size(ind,1),1);

    canvas(ymin:ymax,xmin:xmax)=reshape(sum(A.*B,2)./sqrt(sum((A.^2),2).*sum((B.^2),2)),ymax-ymin+1,[]);

    canvas = canvas(y1+1:z1+y1+1,y2+1:z2+y2+1);

我想我会逐行解释出现问题的地方。

我正在使用

将 numpy 导入为 np

对于 numpy 函数

1.

变量工作正常，直到我到达 python 中的网格线。

    [x,y] = np.mgrid[xmin:xmax,ymin:ymax]

matlab 中的 x 使用测试数据的大小为 517,517 python 中的 x 大小为 516 x 516，所以我用

更改了代码的 python 部分

    xmax=n-(s2-dx-1) + 1
    ymax=m-(s1-dy-1) + 1

    [y,x] = np.mgrid[xmin:xmax,ymin:ymax]

使其与 matlab 代码具有相同的数据集。但我不确定索引是否合理。如果我在 matlab 中有与 numpy 数组完全相同的矩阵，它们是否等效？

2.

matlab 的下一行对我来说是一团糟。

    ind=sub2ind([m,n],y(:),x(:));

对于我正在使用的 x(:) 和 y(:)

    x = np.reshape(x,(x.size,1))
    y = np.reshape(y,(y.size,1))
    x = np.int64(x)
    y = np.int64(y)

对于我在 python 中使用的 sub2ind 函数

    ind = np.ravel_multi_index((y,x), dims=(m,n) )

但这就是数字混乱的地方。

在 matlab 中，对于某个数据集，我得到一个范围为 3723 到 278760 的列向量 并且对于相同的数据集 在 python 中，我得到一个 4264 到 279292 的列向量，中间有不同的步进。

不过，它们的大小都与 (267289,1) 相同。

3.

这行我在 matlab 和 python 中运行良好，我只是把它放在这里，这样我就可以对自己简洁了。

matlab:

    nps=repmat(((-dx+(0:s2-1))*m-dy),s1,1)+repmat((0:(s1-1)).',1,s2);

蟒蛇：

    dtx = (-dx + np.arange(0,s2,1))
    dtx_2 = np.arange(0,s1,1)
    dtx_2 = np.reshape(dtx_2,(dtx_2.size,1))

    nps = np.tile(   dtx*m-dy,(s1,1)  ) + np.tile(   dtx_2  ,(1,s2)  )

(4)。

matlab中的行

    ind=repmat(ind,1,numel(nps))+repmat(nps(:).',numel(ind),1);

在python中我正在尝试

    a = np.tile(ind,(1,nps.size))
    b = np.tile( np.transpose(dtind) , (ind.size,1) )
    ind = a + b

我将它分成 a 和 b 以使其更具可读性。

但这并没有真正按预期工作。

(5)。

我不确定如何在 python 中通过索引访问变量。

在 matlab 中我可以只做 Image(ind)，但是我的代码现在在 python 中没用，因为我找不到替代方法吗？

如果您尝试运行 matlab 代码，请注意一点，如果您在大数据集上运行它，它将导致您的计算机和 matlab 在没有警告的情况下冻结和崩溃。我通过将代码放在一个包装器中来解决这种情况，该包装器索引数据的较小部分以获得防止内存溢出的最终大图像。

我希望我把我乱七八糟的代码弄得够清楚了。这段代码在 matlab 中运行良好，但是 matlab 非常糟糕，主要是因为我不能将我的代码提供给其他人。

编辑：

这个函数是一个矢量化程序，基本上可以：（这是伪代码，所以索引可能不是我的想法）

此段中也没有填充。我使用的加载变量是高斯矩阵或拉帕克矩阵，范围从 6x6 到 64x64。它实际上可以是任何尺寸，只要它小于图像即可。

    correlation_coeficcient_surface = function(Image,loading)
    [m,n] = size(image)
    [u1,u2] = size(loading)
    canvas = zeros(size(image))
    for yii = 1:n
         for xii = 1:m
              image_segment = Image(yii-floor(u1/2):yii+ceil(u1/2),xii-floor(u2/2):xii+ceil(u2/2));
              if(size(image_segment) == size(loading)
                  canvas(yii-floor(u1/2):yii+ceil(u1/2),xii-floor(u2/2):xii+ceil(u2/2)) = corr2(Image_segment,loading);
              end
          end
     end


    end

它必须被矢量化，因为在每个元素上进行迭代使得处理大图像需要很长时间。

编辑编辑：

这是我的过滤器的实际作用。

这是一个示例图像

http://i.imgur.com/o9kV3nK.png

这是我用来关联的示例加载形状

http://i.imgur.com/oYW3k2K.png

这是在我的 matlab 过滤器完成后，图像未对齐，我只是裁剪示例以向您展示它在形状方面的作用。

http://i.imgur.com/aa4ljue.png

这是 scipy.signal.convolve2d，它做了我不打算做的事情。

http://i.imgur.com/vJhXMam.png

Answer 1

除了使用我的其他答案中的提示从头开始翻译代码之外，这看起来就像一个卷积，您可能只使用 convolve2d（此处使用零填充）：

scipy.signal.convolve2d(Image, loading, mode='full', fillvalue=0.0)

如果您想做一些比零填充更花哨的事情，请参阅http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.convolve2d.html。

编辑：有关卷积的更多详细信息，请参阅2D Convolution in Python similar to Matlab's conv2

编辑 2：以下代码应计算滑动窗口上的局部相关系数。它应该正确处理边界。

from scipy.signal import correlate2d
import numpy as np

# Generate random image
Image = np.random.rand(100,100)

# make kernel
x = np.arange(-10,11)
loading = np.exp( -(x[:,None]**2 + x[None,:]**2 )/10)


Image = np.tile(loading, (10,10))
m,n = loading.shape

oneskernel = np.ones(loading.shape)

# get number of points within each box of the sliding window
num = correlate2d(np.ones(Image.shape), oneskernel, mode='same')

# get mean over sliding window
Image_mu = correlate2d(Image, oneskernel, mode='same')/num

# Get sig of sliding window
Image_sig = np.sqrt(correlate2d((Image - Image_mu)**2,
                        oneskernel, mode='same')
                        / num)

loading_mu  =loading.mean() # mean of kernel
loading_sig = loading.std() # sig  of kenrel


# calculate sliding corrcoeff
C = correlate2d(( Image-Image_mu) / Image_sig , (loading -loading_mu)/loading_sig, mode='same') / num

Answer 2

关于 x 和 y - 值的顺序是不同的，如果您尝试将它们展平，您会看到：

x(:)'
1 1 ... 2 2 ....

x.flatten()
array([ 1,  2,  ... 10,  1,  2,...])

即 MATLAB 数组的排列方式类似于 'F' numpy 数组，而不是默认的 'C'。

---

对于小尺寸，我可以匹配 octave 和 numpy：

"""
octave:47> [x,y]=meshgrid(1:3,3:4)
x =
   1   2   3
   1   2   3
y =
   3   3   3
   4   4   4
octave:48> ind=sub2ind([5,5],y,x)
    3    8   13
    4    9   14
"""
Y,X=np.mgrid[2:4,0:3]
"""
array([[2, 2, 2],
       [3, 3, 3]])
array([[0, 1, 2],
       [0, 1, 2]])
"""
ind = np.ravel_multi_index((X,Y),(5,5))
# np.ravel_multi_index((Y,X),(5,5),order='F')
"""
array([[ 2,  7, 12],
       [ 3,  8, 13]])

"""

我对@9​​87654325@ 的问题感到困惑。 Image 是一个 [m,n] 数组，对吧？

Answer 3

我认为你应该放慢速度，阅读一些有关 Python 数组基础知识的材料，例如索引和广播。首先，我会在http://www.sam.math.ethz.ch/~raoulb/teaching/PythonTutorial/intro_numpy.html 阅读基础教程。 http://mathesaurus.sourceforge.net/matlab-numpy.html 还包含一个表格，其中包含某些 numpy 和 matlab 操作之间的对应关系。但是，总的来说，我会保持开放的心态，并意识到 matlab 方式通常不是 numpy 中的最佳方式。

我不会直接回答你的所有问题，但这里有以下想法。

1. Python 索引为零索引，.因此 matlab arr(1:100) 与 numpy arr[0:100] 相同。
2. 您可以使用逻辑数组或索引数组来索引 numpy 数组，就像在 matlab 中一样

3. 一般来说，repmat 的功能由智能广播在 numpy 中处理，不需要手动调用tile。例如，以下代码创建一个随机 100x100 数组，计算行均值，然后从行中减去行均值（例如，将数据居中）：

   arr = np.random.rand(100,100)
   mu  = arr.mean(axis=1)
   centered = arr - mu[:,None]
   

   mu[:,None] 数组的大小为 (100,1)，numpy 足够聪明，可以将其“广播”到大小 (100,100) 以计算 centered。继续该示例，mu[:,None,None] 的大小为 (100,1,1)。

4. Matlab 的 size(arr) 与 numpy 的 arr.shape 相同。

祝你好运！

编辑：例如，您可以更简洁地执行 #3：

nps = (-dx+np.arange(s2)*m -du)[None,:] + np.arange(s1)[:,None]

还有#4：

ind=  ind[:, None] + nps.ravel()[None, :]