arrays - 将矢量化的 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 中的 meshgrid 行。

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

x 在 matlab 中使用测试数据的大小为 517,517python 中的 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 中工作得很好我只是把它放在这里所以我可以对自己简洁。

数学实验:

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

python :

    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

我认为你应该放慢速度，阅读一些有关 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) 以进行计算居中。继续该示例，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, :]