Computing the sum of squared differences 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.
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 &amp; T_dist==0 &amp; Y_dist==0)
    4     4     4     4     4     4

D_diffsqrd(X_dist==0 &amp; T_dist==0 &amp; Y_dist==2)
    36    36    36    36    36    36

D_diffsqrd(X_dist==0 &amp; T_dist==1 &amp; Y_dist==0)
    64     4    16    36     0    36    16     4    64