向量化多个高斯计算

Vectorizing multiple gaussian calculations

我正在尝试矢量化我的一些代码,这些代码在图像上添加了许多高斯分布的强度。我目前为每个高斯循环遍历函数 'gaussIt2D',它针对单个 2D 高斯进行矢量化:

windowSize=10;
imSize=[512,512];
%pointsR is an nx2 array of coordinates [x1,y1;x2,y2;...;xn,yn]
pointsR=rand(100,2)*511+1;
%sigmaR is the standard deviation of the gaussian being created
sigmaR = 1;

outputImage=zeros(imSize);

for n=1:size(pointsR,1)
    rangeX = floor(pointsR(n,1)-windowSize):ceil(pointsR(n,1)+ windowSize);
    rangeX = rangeX(rangeX > 0 & rangeX <= imSize(1));
    rangeY = floor(pointsR(n,2)-windowSize):ceil(pointsR(n,2)+windowSize);
    rangeY = rangeY(rangeY > 0 & rangeY <= imSize(2));
    outputImage(rangeX,rangeY) = outputImage(rangeX,rangeY)+gaussIt2D(rangeX(1),rangeX(end),rangeY(1),rangeY(end),[sigmaR,pointsR(n,1),pointsR(n,2)]);
end

function [result] = gaussIt2D(xInit,xFinal,yInit,yFinal,sigma,xCenter,yCenter)
    %Returns gaussian intenisty values for the region defined by [xInit:xFinal,yInit:yFinal] using the gaussian properties sigma,centerX,centerY

    [gridX,gridY]=ndgrid(xInit:xFinal,yInit:yFinal);

    result=exp( -( (gridX-xCenter).^2 + (gridY-yCenter).^2 ) ./ (2*sigma.^2) );

end

我正在尝试通过允许 gaussIt2D 函数接受 x 和 y 值的向量以及 x 和 y 中心的向量并将所有这些放在一起来进一步向量化此过程。到目前为止,我的思考过程一直是尝试堆叠网格并复制中心并进行逐元素高斯计算。对于(简化的)示例,如果:

xInits = [1,2,3];
xFinals = [2,3,4];
xCenters = [1.2,2.8,3.1];
yInits = [1,2,3];
yFinals = [2,3,4];
yCenters = [1.5,2.4,3.6];

然后我想按照以下形式创建网格和中心:

gridX = [1,2
         1,2
         2,3
         2,3
         3,4
         3,4]

xCenters = [1.2,1.2
            1.2,1.2
            2.8,2.8
            2.8,2.8
            3.1,3.1
            3.1,3.1]

然后可以将其用于原始函数中使用的相同高斯方程。但是,生成这些数组让我感到困惑。我现在拥有的是:

function [result]=gaussIt2DVectorized(xInits,xFinals,yInits,yFinals,sigmas,xCenters,yCenters)
    %Incomplete
    %Returns gaussian intenisty values for the region defined by 
    %[xInit:xFinal,yInit:yFinal] using the values array:[sigma,centerX,centerY]
    [gridX,gridY]=arrayfun('ndgrid',xInits:xFinals,yInits:yFinals);
    xCenters = repelem(xCenters,numel(xInits(1):xFinals(1)), numel(yInits(1):yFinals(1)));
    yCenters = repelem(yCenters,numel(xInits(1):xFinals(1)), numel(yInits(1):yFinals(1)));

    result=exp( -( (gridX-xCenters).^2 + (gridY-yCenters).^2 ) ./ (2*sigmas^2) );

end

虽然这实际上行不通,但我也预计难以解释不同长度的范围(即 xInit:xFinal)。

我们将提供任何帮助、提示或替代方法。

谢谢。

由于您无法确定所有网格的大小都相同,因此最好将它们存储在元胞数组中,而不是将它们堆叠在矩阵中。使用元胞数组,您仍然可以 运行 您的计算而无需使用 cellfun.

循环

例如:

function [result] = gaussIt2D_better(xInits,xFinals,yInits,yFinals,sigmas,xCenters,yCenters)
    [gridsX, gridsY] = arrayfun(@(x) ndgrid(xInits(x):xFinals(x), yInits(x):yFinals(x)),1:length(xInits),'UniformOutput',0);
    f=@(gridX, gridY, xCenter, yCenter, sigma) exp( -( (gridX-xCenter).^2 + (gridY-yCenter).^2 ) ./ (2*sigma.^2) );
    result=cellfun(f, gridsX, gridsY, num2cell(xCenters), num2cell(yCenters), num2cell(sigmas), 'UniformOutput',0);
end

请注意,在此示例中,返回值是一个与输入向量长度相同的元胞数组,每个向量一个结果。