Computing the sum of squared differences

Question

About Data: 3D data (100x100x20).

Computation/Formula and Method: I am trying to compute the sum of squared differences along rows, columns and angles for various time differences.

%% Grid and time paramters
%# Grid parameters
nRows=100;
nCol=100;
InitLag_Row=0;
InitLag_Col=0;
InitLag_Ang=0;
LessLags=20;
nLags_Row=nRows-LessLags;
nLags_Col=nCol-LessLags;
nLags_Ang=min(nRows,nCol)-LessLags;
%# Time parameters
T=1:1:20;
nT=numel(T);
%# Load concentrations data file with 100x100x20 dimensions
load('c.mat')
%% Shift values along rows for all columns and all time steps
for hRow=InitLag_Row:nLags_Row
 c_ShiftedRow(:,:,:,hRow+1)=circshift(c(:,:,:),[-hRow 0]);
end
%% Shift values along columns for all rows and all time steps
for hCol=InitLag_Col:nLags_Col
 c_ShiftedCol(:,:,:,hCol+1)=circshift(c(:,:,:),[0 -hCol]);
end
%% Shift values along NW-SE and NE-SW directions for all time steps
for hAng=InitLag_Ang:nLags_Ang
 c_ShiftedNWSE(:,:,:,hAng+1)=circshift(c(:,:,:),[-hAng -hAng]);
 c_ShiftedNESW(:,:,:,hAng+1)=circshift(c(:,:,:),[-hAng hAng]);
end
idel_t=1; % initialize index for zero time lag
for del_t=0:nT-1 % time lag
 for nLags_t=1:nT-del_t; % number of time lags formed with del_t lag value
 for hRow=InitLag_Row:nLags_Row
 variogramTemp_hRow3D=(cumsum(((c(1:end-hRow,:,nLags_t)-...
 c_ShiftedRow(1:end-hRow,:,nLags_t+del_t,hRow+1)).^2),1))./...
 (2*size(c(1:end-hRow,:,nLags_t),1));
 variogram_hRow3D(hRow+1,:,nLags_t,idel_t)=variogramTemp_hRow3D(end,:)./...
 mean(var(c(1:end-hRow,:,[nLags_t,nLags_t+del_t]),0,1),3);
 end
 %% Variogram along columns for all images
 for hCol=InitLag_Col:nLags_Col 
 variogramTemp_hCol3D=(cumsum(((c(:,1:end-hCol,nLags_t)-...
 c_ShiftedCol(:,1:end-hCol,nLags_t+del_t,hCol+1)).^2),2))./...
 (2*size(c(:,1:end-hCol,nLags_t),2));
 %# normalized variogram values at all times for different lags at all 'nRows'
 variogram_hCol3D(:,hCol+1,nLags_t,idel_t)=variogramTemp_hCol3D(:,end)./...
 mean(var(c(:,1:end-hCol,[nLags_t,nLags_t+del_t]),0,2),3);
 end
 %% Variogram along NW-SE and NE-SW directions for all images
 for hAng=InitLag_Ang:nLags_Ang 
 variogramTemp_h3DNWSE=(cumsum(((c(1:end-hAng,1:end-hAng,nLags_t)-...
 c_ShiftedNWSE(1:end-hAng,1:end-hAng,nLags_t+del_t,hAng+1)).^2),1))./...
 (2*size(c(1:end-hAng,1:end-hAng,nLags_t),1));
 variogramTemp_h3DNESW=(cumsum(((c(1:end-hAng,end:-1:1+hAng,nLags_t)-...
 c_ShiftedNESW(1:end-hAng,end:-1:1+hAng,nLags_t+del_t,hAng+1)).^2),1))./...
 (2*size(c(1:end-hAng,end:-1:1+hAng,nLags_t),1));
 %# normalized variogram values/cumulative sum at all times for different lags along NW-SE and NE-SW directions
 variogram_h3DNWSE(hAng+1,1:nCol-hAng,nLags_t,idel_t)=variogramTemp_h3DNWSE(end,:)./...
 mean(var(c(1:end-hAng,1:end-hAng,[nLags_t,nLags_t+del_t]),0,1),3);
 variogram_h3DNESW(hAng+1,1:nCol-hAng,nLags_t,idel_t)=variogramTemp_h3DNESW(end,:)./...
 mean(var(c(1:end-hAng,end:-1:1+hAng,[nLags_t,nLags_t+del_t]),0,1),3);
 end
 end
 %# Change the zero matrices after nT-del_t to NaN
 variogram_hRow3D(:,:,nT-del_t+1:end,idel_t)=nan;
 variogram_hCol3D(:,:,nT-del_t+1:end,idel_t)=nan;
 variogram_h3DNWSE(:,:,nT-del_t+1:end,idel_t)=nan;
 variogram_h3DNESW(:,:,nT-del_t+1:end,idel_t)=nan;
 %# change to next time lag index
 idel_t=idel_t+1;
end

I would appreciate any help on ways to make this code faster and more efficient.

• tmpearce: As a start, check out pdist - you can probably use that to do the bulk of it.
• user2947: tmpearce: Thanks for the function. I was also thinking to use circshift(), even though it's just a small step for the problem.
• tmpearce: Also, if you could come up with a small example (3 by 3 by 2 perhaps) and the expected results, people may be more interested in playing around with the problem and coming up with a good solution.
• user2947: I've added some explanation and the rest I'd try to do tomorrow. Please let me know if there's something in the question which is not clear.
• user2947: @tmpearce: I've added some illustration to explain my specific problem and how it is done.

Answer 1

Here's a way to do it using pdist to generate both the differences in value, plus logical indices that you can use to select which distances you wish to look at/use further.

D1 = [1 2 3;4 5 6;7 8 9];
D2 = [9 8 7;6 5 4;3 2 1];
D = cat(3,D1,D2);
D(:,:,1) =
 1 2 3
 4 5 6
 7 8 9
D(:,:,2) =
 9 8 7
 6 5 4
 3 2 1
Y = repmat([1 2 3]', [1 3 2]); %# (' in comment to fix SO highlighting)
X = repmat([1 2 3],[3 1 2]);
T = repmat(cat(3,1,2),[3 3 1]);
D_diffsqrd = pdist(D(:),'euclidean').^2;
X_dist = pdist(X(:),'euclidean');
Y_dist = pdist(Y(:),'euclidean');
T_dist = pdist(T(:),'euclidean');
Angle_dist = pdist([X(:) Y(:)],'euclidean');
D_diffsqrd(X_dist==2 & T_dist==0 & Y_dist==0)
 4 4 4 4 4 4
D_diffsqrd(X_dist==0 & T_dist==0 & Y_dist==2)
 36 36 36 36 36 36
D_diffsqrd(X_dist==0 & T_dist==1 & Y_dist==0)
 64 4 16 36 0 36 16 4 64

• user2947: tmpearce: I did the ones for Row_Lag and Col_Lag and I am getting some results, though have to make sure they are correct. I did for Angle_Lab on a single image (2D) but have not tried yet for the 3D. I used circshift() to shift the numbers to take differences. Will post the updated question later. Thanks for your efforts!!