# Computing the sum of squared differences of the concentration profile

About Data: I have a 3D data (100x100x20) of space and time (`c`) where first 2 dimensions are concentration values along space while 3rd dimension contains concentration values at different time from 1 till 20.

Computation/Formula and Method: I am trying to compute the sum of squared differences of the concentration profile, also called variogram, 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

%% 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

%% Variogram analysis along space and time

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
%% Variogram along rows for all images
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));

%# normalized variogram values at all times for different lags at all 'nCol'
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.

-
As a start, check out `pdist` - you can probably use that to do the bulk of it. –  tmpearce Jul 18 '12 at 1:40
tmpearce: Thanks for the function. I was also thinking to use `circshift()`, even though it's just a small step for the problem. –  Pupil Jul 18 '12 at 1:59
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. –  tmpearce Jul 18 '12 at 4:17
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. –  Pupil Jul 18 '12 at 5:56
@tmpearce: I've added some illustration to explain my specific problem and how it is done. –  Pupil Jul 18 '12 at 21:10

## migrated from stackoverflow.comJul 23 '12 at 21:53

This question came from our site for professional and enthusiast programmers.

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
``````
-
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!! –  Pupil Jul 19 '12 at 18:00

It's not a complete answer, but because of the generic form of the question I would suggest that you use matlab profiler, and after you get the results you can ask a more specific question.

-