求平方差之和

Question

关于数据: 3D数据(100x100x20)。

计算/公式和方法:我试图计算不同时差沿行、列和角度的平方差之和。

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

对于如何使代码更快、更高效，我将不胜感激。

Answer 1

下面是一种使用pdist生成值差异的方法，以及可以用来选择希望进一步查看/使用的距离的逻辑索引。

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