使用 bsxfun 在 matlab 中矢量化嵌套循环，并使此代码通过 matlab 运行到 GPU 中

Question

HAUSDIM 返回对象的豪斯多夫分形维数，表示如下二进制图像。我需要并行 GPU 代码来使用并行计算工具箱从 matlab 在 GPU 上测试它。

算法

1. 用背景像素填充图像，使其尺寸为2的幂。
2. 将框大小"e"设置为图像的大小。
3. 计算 N（e），对应于大小为"e"的盒子数量至少包含一个对象像素。
4. 如果 e> 1，则 e = e/2 并重复步骤 3。
5. 计算点对数（N（e）） x log（1/e） 并使用最小二乘法将线拟合到点的方法。
6. 返回的豪斯多夫分形维数 D 是直线的斜率。

法典

function [ D ] = hausDim( I )
maxDim = max(size(I));
newDimSize = 2^ceil(log2(maxDim));
rowPad = newDimSize - size(I, 1);
colPad = newDimSize - size(I, 2);
I = padarray(I, [rowPad, colPad], 'post');
boxCounts = zeros(1, ceil(log2(maxDim)));
resolutions = zeros(1, ceil(log2(maxDim)));
iSize = size(I, 1);
boxSize = iSize;
boxesPerDim = 1;
idx = 0;
while boxSize >= 1
    boxCount = 0;
    minBox = (1: boxSize: (iSize - boxSize) + 1);
    maxBox = (boxSize: boxSize: iSize);
    for boxRow = 1:boxesPerDim
        for boxCol = 1:boxesPerDim
            objFound = false;
            for row = minBox(boxRow) : maxBox(boxRow)
                for col = minBox(boxCol) : maxBox(boxCol)
                    if I(row, col)
                        boxCount = boxCount + 1;
                        objFound = true; % Break from nested loop.
                        break;
                    end;
                end;
                if objFound
                    break; % Break from nested loop.
                end;
            end;
        end;
    end;
    idx = idx + 1;
    boxCounts(idx) = boxCount;
    resolutions(idx) = 1 / boxSize;
    boxesPerDim = boxesPerDim * 2;
    boxSize = boxSize / 2;
end;
D = polyfit(log(resolutions), log(boxCounts), 1);
D = D(1);
end

Answer 1

这篇文章中列出的是一种矢量化方法，除了最外层的循环之外，所有循环都会被杀死，由于每次这样的迭代所需的不同形状的盒子使其变得复杂，可能不值得付出努力。所以，这是用零填充输入数组后的实现 -

Inz = I~=0;
nrows = size(I,1);
N = log2(nrows)+1;
boxCounts_out = zeros(1,N);
for iter = 1:N
    n = 2.^(iter-1);
    matches = any(any(reshape(Inz,nrows/n,n,nrows/n,n),1),3);
    boxCounts_out(iter) = sum(matches(:));
end
resolutions_out = 1./(2.^(N-1:-1:0));

之后，您可以使用polyfit代码 boxCounts_out 和 resolutions_out .

---

为了使代码可在 GPU 上运行，我绝对建议将Inz和boxCounts_out转换为gpuArrays，然后使用前面列出的代码。因此，像这样执行初始化 -

Inz = gpuArray(I~=0);
boxCounts_out = zeros(1,N,'gpuArray');

由于使用这些输入的操作涉及大量的求和减少和相对较少的内存，我真的会打赌那里有良好的运行时和内存效率。