2017-06-16 166 views
1

我需要計算數百萬個粒子的歸一化互相關的最大值。 normxcorr2的兩個參數的大小是56 * 56。我無法並行計算。是否有任何建議加快代碼,特別是我不需要所有的結果,但只有每個互相關的最大值(知道位移)?該算法加速normxcorr2的最大值的計算

%The choice of 170 particles is because in each time 
%the code detects 170 particles, so over 10000 images it's 1 700 000 particles 
particle_1=rand(54,54,170); 
particle_2=rand(56,56,170); 
for i=1:170 
    C=normxcorr2(particle_1(:,:,i),particle_2(:,:,i)); 
    L(i)=max(C(:)); 
end 
+2

NORMXCORR2歸一化的二維互相關。 C = NORMXCORR2(模板,A)計算矩陣TEMPLATE和A的歸一化互相關。矩陣A *必須大於矩陣 模板使歸一化有意義。 – Jed

+0

在內存預分配的循環前加'L =零(1,170);'會加快一點。但是,正如Jed所說,目前還存在一個概念性問題 –

+0

@Jed是的,你說得對,我在大小上犯了一個錯誤,但問題仍然是一樣的 – ransa

回答

1

我沒有MATLAB的

例子,所以我跑在這個網站下面的代碼:https://www.tutorialspoint.com/execute_matlab_online.php這實際上是八度。所以,我實現了「天真」歸一化互相關性,的確是這些小圖像尺寸的幼稚表現更好:

經過時間是2.62645秒 - 爲normxcorr2
經過時間是0.199034秒 - 我naive_normxcorr2

代碼是基於文章http://scribblethink.org/Work/nvisionInterface/nip.pdf,它描述瞭如何使用integral image以有效的方式計算標準化所需的標準偏差,這是box_corr函數。

此外,MATLAB的normxcorr2返回一個填充圖像,所以我採取了無襯墊部分的最大值。

pkg load image 

function [N] = naive_corr(pat,img) 

[n,m] = size(img); 
[np,mp] = size(pat); 

N = zeros(n-np+1,m-mp+1); 

for i = 1:n-np+1 
    for j = 1:m-mp+1 
     N(i,j) = sum(dot(pat,img(i:i+np-1,j:j+mp-1))); 
    end 
end 

end 


%w_arr the array of coefficients for the boxes 
%box_arr of size [k,4] where k is the number boxes, each box represented by 
%4 something ... 
function [C] = box_corr2(img,box_arr,w_arr,n_p,m_p) 

% construct integral image + zeros pad (for boundary problems) 
I = cumsum(cumsum(img,2),1); 
I = [zeros(1,size(I,2)+2); [zeros(size(I,1),1) I zeros(size(I,1),1)]; zeros(1,size(I,2)+2)]; 

% initialize result matrix 
[n,m] = size(img); 
C = zeros(n-n_p+1,m-m_p+1); 
%C = zeros(n,m); 

jump_x = 1; 
jump_y = 1; 

x_start = ceil(n_p/2); 
x_end = n-x_start+mod(n_p,2); 
x_span = x_start:jump_x:x_end; 

y_start = ceil(m_p/2); 
y_end = m-y_start+mod(m_p,2); 
y_span = y_start:jump_y:y_end; 

arr_a = box_arr(:,1) - x_start; 
arr_b = box_arr(:,2) - x_start+1; 
arr_c = box_arr(:,3) - y_start; 
arr_d = box_arr(:,4) - y_start+1; 

% cumulate box responses 
k = size(box_arr,1); % == numel(w_arr) 
for i = 1:k 
    a = arr_a(i); 
    b = arr_b(i); 
    c = arr_c(i); 
    d = arr_d(i); 

    C = C ... 
     + w_arr(i) * (I(x_span+b,y_span+d) ... 
         - I(x_span+b,y_span+c) ... 
         - I(x_span+a,y_span+d) ... 
         + I(x_span+a,y_span+c)); 
end 

end 



function [NCC] = naive_normxcorr2(temp,img) 

    [n_p,m_p]=size(temp); 

    M = n_p*m_p; 

    % compute template mean & std 
    temp_mean = mean(temp(:)); 
    temp = temp - temp_mean; 

    temp_std = sqrt(sum(temp(:).^2)/M); 

    % compute windows' mean & std 
    wins_mean = box_corr2(img,[1,n_p,1,m_p],1/M, n_p,m_p); 
    wins_mean2 = box_corr2(img.^2,[1,n_p,1,m_p],1/M,n_p,m_p); 


    wins_std = real(sqrt(wins_mean2 - wins_mean.^2)); 
    NCC_naive = naive_corr(temp,img); 

    NCC = NCC_naive ./ (M .* temp_std .* wins_std); 
end 

n = 170; 

particle_1=rand(54,54,n); 
particle_2=rand(56,56,n); 

[n_p1,m_p1,c_p1]=size(particle_1); 
[n_p2,m_p2,c_p2]=size(particle_2); 

L1 = zeros(n,1); 
L2 = zeros (n,1); 

tic 
for i=1:n 
    C1=normxcorr2(particle_1(:,:,i),particle_2(:,:,i)); 

    C1_unpadded = C1(n_p1:n_p2 , m_p1:m_p2); 
    L1(i)=max(C1_unpadded(:)); 

end 

toc 

tic 
for i=1:n 

    C2=naive_normxcorr2(particle_1(:,:,i),particle_2(:,:,i)); 
    L2(i)=max(C2(:)); 
end 

toc 
+0

另一件事,因爲你只需要最大值而不是整個NCC,通過天真的實現,你可以決定最大值的閾值。在天真實施的迭代過程中,一旦得到滿足閾值的值,就可以停止。如果閾值是好的,那麼你可以假設你找到的第一個位置是最大值。它可能並不總是成功,但大部分時間。 – dafnahaktana

+0

謝謝你的回答。最後,我也只是通過計算互相關來計算互相關,這也是我感興趣的部分。它加快了速度(對於5x5的區域來說,係數爲10) – ransa