向量化多個高斯計算

Question

我正在嘗試將我的一些代碼向量化，這些代碼在圖像上添加了許多高斯分佈的強度。我目前循環在函數「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）不同的長度。

任何幫助，提示或替代方法將appriciated。

謝謝。

Answer 1

由於您不能確定網格的大小都是相同的，因此最好將它們存儲在單元陣列中，而不是將它們堆疊在矩陣中。使用單元格數組，您仍然可以運行計算而無需使用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 

注意，在該示例中，所返回的值是一個單元陣列具有相同的長度作爲輸入向量，一個結果爲每個。