Vectorization help
in [Matlab]
|
Prev: How to? Excel RANK function in Matlab?
Next: How to Find Surface Normal of each triangle in a Surface.
From: Jussi Forma on 31 Mar 2010 05:50 Hi, I have 4d-data, namely certain measurement vector V(t) for each point in x,y,z grid. From this set, I need to find one point x1,y1,z1 which is closest to my given reference vector, similarity metric is sum of squared errors between reference and data vector. Data vectors may be shifted also in t-dimension, all vectors are periodic. So in essence I have one vector, and bunch of candidates to compare it with different lags. Here's my first version, let me know how to optimize this for speed. %%%%%%%%%%%%%%%%%%%%% refvec = reshape(rand(100,1),[1 1 1 100]); %reference vector refblock = repmat(refvec,[5 5 5 1]); %replicate it to correspond the whole spatial %grid floatblock = rand([5 5 5 100]); %example data set maxshift = 10; % maximum shift of vectors in t-domain. (not in xyz) %calculate MSE for each vector pair, and shift. MSEs = zeros(size(refblock,1),size(refblock,2),size(refblock,3),2*maxshift+1); % for each shift, get squared errors for the whole spatial grid. for n = -maxshift : maxshift MSEs(:,:,:,n+maxshift+1) = sum(( refblock - circshift(floatblock,[0 0 0 n]) ).^2,4); end
From: Bruno Luong on 31 Mar 2010 06:35 "Jussi Forma" <juss911(a)luukku.com> wrote in message <hov5sc$5pg$1(a)fred.mathworks.com>... > ... sum(( refblock - circshift(floatblock,[0 0 0 n]) ).^2,4); You can rewrite this as sum(( refblock - circshift(floatblock,[0 0 0 n]) ).^2,4) + sum((refblock).^2,4) + sum((floatblock).^2,4) - 2* sum(refblock.*circshift(floatblock,[0 0 0 n])),4) The two first terms do not depend on n and can be computed once. The third term can be rewritten as circular convolution, which can be computed using fft on the fourth dimension (not sure if it's faster because you only need few elements of the convolution result, I let you play with it). Bruno
From: Jussi Forma on 1 Apr 2010 04:45 "Bruno Luong" <b.luong(a)fogale.findmycountry> wrote in message <hov8gp$cc8$1(a)fred.mathworks.com>... > "Jussi Forma" <juss911(a)luukku.com> wrote in message <hov5sc$5pg$1(a)fred.mathworks.com>... > > > ... sum(( refblock - circshift(floatblock,[0 0 0 n]) ).^2,4); > > You can rewrite this as > > sum(( refblock - circshift(floatblock,[0 0 0 n]) ).^2,4) > + sum((refblock).^2,4) + sum((floatblock).^2,4) > - 2* sum(refblock.*circshift(floatblock,[0 0 0 n])),4) > > The two first terms do not depend on n and can be computed once. The third term can be rewritten as circular convolution, which can be computed using fft on the fourth dimension (not sure if it's faster because you only need few elements of the convolution result, I let you play with it). > > Bruno Thanks a lot Bruno Taking those terms out from the loop helped, and calculating the last term with fft's and picking the right elements, roughly doubled the speed from my first attempt. Here's how I calculated the last term: term = -2*ifft(fft(refblock,[],4).*fft(flipdim(floatblock,4),[],4),[],4); term = flipdim(term(:,:,:,end-(2*maxshift+1)+1:end) ,4)
From: Bruno Luong on 1 Apr 2010 07:56
"Jussi Forma" <juss911(a)luukku.com> wrote in message <hp1mei$5ar$1(a)fred.mathworks.com>... > > term = -2*ifft(fft(refblock,[],4).*fft(flipdim(floatblock,4),[],4),[],4); Because refblock is repmat of refvec, it should be faster to compute quantity fft(refblock,[],4) by: fftref = fft(refvec); fftref = reshape(fftref ,[1 1 1 100]); % if needed fftref= repmat(fftref ,[5 5 5 1]); % Bruno |