3D Cross Correlation for two matrices

12 views (last 30 days)
Hi all.
I want to do some 3D Matrix cross correlation for a time of flight estimation of a flow imaging technique called Wire Mesh Sensor.
What i get as an input is two 3D matrices that resemble a consecutive 2D distribution of void fraction distributions in a two phase flow. The two matrices are measured at a short distance after each other. The values in the matrix vary between 0 and 1 in each matrix. In the ideal case, we have the exact same content in each matrix, but frame shifted (so, the same patterns can be found in both matrixes at different positions in the z-direction).
I was experimenting with the xcorr, xcorr2 functions and also the convn function. convn does work and gives the right values for the cross correlation (by just calling it with the 2nd matrix flipped, so convn(A,B(end:-1:1,...)).
Unfortunatelly, the correlation has to be weighted in order to find the similarity according to the signal content. If this is not taken into account, erreneous results occur. What i would need is the same option that is available in the simple, 1D xcorr function. The method there is called coeff (Normalizes the sequenceso the autocorrelations at zero lag are identically 1).
Unfortunatelly, this is as far as i know not available in the convn version. So, i tried to program a simple code on my own, yet just checking for the general performance for the unbiased cross correlation.
I tried two versions:
function [ CCM ] = cc_3D( A,B,method )
%%CC_3D 3D Matrix Cross Correlation
% Calclates the 3D Cross Correlation for matrices A and B
%%Inputs
%{
A: 1st Matrix for Cross Correlation
B: 2nd Cross Correlation
method: Method for the cross correlation
(1): 'raw': Nomal Crosscorrelation
%}
%%Outputs
%{
CCM: Cross Correlation Matrix
%}
%%Check and Patch the matrices
as=size(A);
bs=size(B);
if bs>=as %Matrix B larger than A
A = padarray(A,(bs-as)/2); % Pad the arrays to the same size
elseif as>=bs
B = padarray(B,(as-bs)/2); % Pad the arrays to the same size
end
% if the matrices have the same size, nothing has to be done
as=size(A);
CCM=zeros(2*as(1)-1,2*as(2)-1,2*as(3)-1);
%%Cross-Correlation Calculation
% Data Loops for all the values to be calculated, given 3D matrix with
% size x:size(A,1)*2-1; y:size(A,2)*2-1; z:size(A,3)*2-1 [cc].
%
if strcmp(method,'raw');
for i=1:as(1)*2-1;
for j=1:as(2)*2-1;
for k=1:as(3)*2-1;
%%Search for the summation boundary values
%%%%%%%%%%%%%%%% x Values
if i<=as(1)
axi=1; axt=i;
bxi=as(1)-i+1; bxt=as(1);
elseif i>as(1)
axi=i-as(1)+1; axt=as(1);
bxi=1; bxt=2*as(1)-i;
end
%%%%%%%%%%%%%%%% y Values
if j<=as(2)
ayi=1; ayt=j;
byi=as(2)-j+1; byt=as(2);
elseif j>as(2)
ayi=j-as(2)+1; ayt=as(2);
byi=1; byt=2*as(2)-j;
end
%%%%%%%%%%%%%%%% z Values
if k<=as(3)
azi=1; azt=k;
bzi=as(3)-k+1; bzt=as(3);
elseif k>as(3)
azi=k-as(3)+1; azt=as(3);
bzi=1; bzt=2*as(3)-k;
end
%%Cross Correlation Operation
% First, Calculate the vector product of the two matrices
CDummy=...
A(axi:axt,ayi:ayt,azi:azt).*B(bxi:bxt,byi:byt,bzi:bzt);
CCM(i,j,k)=sum(CDummy(:));
end
end
end
end
end
This first version goes through the loop with 3 nested for loops - not optimal for sure.
The second version vectorizes - but probably in non optimal way as well.
function [ CCM ] = cc_3D_LINE( A,B,method )
%%CC_3D 3D Matrix Cross Correlation
% Calclates the 3D Cross Correlation for matrices A and B
%%Inputs
%{
A: 1st Matrix for Cross Correlation
B: 2nd Cross Correlation
method: Method for the cross correlation
(1): 'raw': Nomal Crosscorrelation
%}
%%Outputs
%{
CCM: Cross Correlation Matrix
%}
%%Check and Patch the matrices
as=size(A);
bs=size(B);
if bs>=as %Matrix B larger than A
A = padarray(A,(bs-as)/2); % Pad the arrays to the same size
elseif as>=bs
B = padarray(B,(as-bs)/2); % Pad the arrays to the same size
end
% if the matrices have the same size, nothing has to be done
as=size(A);
CCM=zeros(2*as(1)-1,2*as(2)-1,2*as(3)-1);
%%Cross-Correlation Calculation
% Data Loops for all the values to be calculated, given 3D matrix with
% size x:size(A,1)*2-1; y:size(A,2)*2-1; z:size(A,3)*2-1 [cc].
%
if strcmp(method,'raw');
for ii=1:((as(1)*2-1)*(as(2)*2-1)*(as(3)*2-1));
[i,j,k]=ind2sub(size(CCM),ii);
%%Search for the summation boundary values
%%%%%%%%%%%%%%%% x Values
if i<=as(1)
axi=1; axt=i;
bxi=as(1)-i+1; bxt=as(1);
elseif i>as(1)
axi=i-as(1)+1; axt=as(1);
bxi=1; bxt=2*as(1)-i;
end
%%%%%%%%%%%%%%%% y Values
if j<=as(2)
ayi=1; ayt=j;
byi=as(2)-j+1; byt=as(2);
elseif j>as(2)
ayi=j-as(2)+1; ayt=as(2);
byi=1; byt=2*as(2)-j;
end
%%%%%%%%%%%%%%%% z Values
if k<=as(3)
azi=1; azt=k;
bzi=as(3)-k+1; bzt=as(3);
elseif k>as(3)
azi=k-as(3)+1; azt=as(3);
bzi=1; bzt=2*as(3)-k;
end
%%Cross Correlation Operation
% First, Calculate the vector product of the two matrices
CDummy=...
A(axi:axt,ayi:ayt,azi:azt).*B(bxi:bxt,byi:byt,bzi:bzt);
CCM(i,j,k)=sum(CDummy(:));
end
end
end
Basically, the first thing the code does is padding the matrices to the same size if the input matrices are not. Then, it just starts to look for the correct start and end points for all the matrix sums.
I don't know if this is the way to go or if it can be done easier - if somebody could probably throw a hint to me to improve the code, it would be great.
Thanks in advance!
Best

Answers (0)

Categories

Find more on Environment and Settings in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!