1
\$\begingroup\$

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.

\$\endgroup\$

migrated from stackoverflow.com Jul 23 '12 at 21:53

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

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

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
\$\endgroup\$
  • \$\begingroup\$ 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!! \$\endgroup\$ – user2947 Jul 19 '12 at 18:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy