please help me simplify the code

I have wirtten a program about measure the distance between two gray image. However,the code rather cost time,and I have no idea how to simplify it. Could anyone give me some advice on this?

The code attatch as follow:

clc;
clear;
theta=1;
X=imread('C:\Users\lwh\Desktop\aerial2.pgm');
Y=imread('C:\Users\lwh\Desktop\astro1.pgm');
[M,N]=size(X);
[M,N]=size(Y);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%input image%%%%%%%%%%%%%%%%%%%%%%%
P=reshape(X',1,M*N);
S=zeros(3,M*N);
    for j=1:M*N
       S(1,j)=P(j);
    end
r=1;  
for K1=0:M-1
    for L1=0:N-1
        %X(t)=K1*N+L1;
        t=2; 
             S(t,r)=K1;
             S(t+1,r)=L1;
        r=r+1;
    end 
end
%%%%%%%%%%%%%%%%%%%%%%%%%quare image%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q=reshape(Y',1,M*N);
T=zeros(3,M*N);
  for j=1:M*N
       T(1,j)=Q(j);
  end
r=1;  
for K2=0:M-1
    for L2=0:N-1
        %X(t)=K1*N+L1;
        t=2; 
             T(t,r)=K2;
             T(t+1,r)=L2;
        r=r+1;
    end 
end                       %%%%line 43
V=zeros(1,M*N);L=zeros(1,M*N);GT=0;
VV=zeros(1,M*N);LL=zeros(1,M*N);U=zeros(1,M*N);
for t=1:M*N      %%%line 73
   tic
        V(1,:)=S(2,t);
        L(1,:)=T(2,:) ; 
        VV(1,:)=S(3,t);
        LL(1,:)=T(3,:);
        U(1,:)=(S(1,t)-T(1,t))*(S(1,:)-T(1,:));

          VL=V-L;VVLL=VV-LL;
          I=(VL.^2+VVLL.^2).^(0.5); 
      gij=exp((-I.^2)/2)/2*pi; G=gij*U';
      GT=G+GT;
           toc
  end            %%%line 88

The code is very time costly from line73~88. because the max of the t is M*N ,so if the image size is very large, it hard to get result.I just want to start simplifing it from the line 73 .

Please help, many thanks!

Answer 1

[Jan]

Avoid all repeated computations in the loop:

% No pre-allocation for: V=zeros(1,M*N);L=zeros(1,M*N);
% No pre-allocation for: VV=zeros(1,M*N);LL=zeros(1,M*N);U=zeros(1,M*N);

GT    = 0;
S1mT1 = S(1, :) - T(1, :);
c     = 1 / 2 * pi; 
for t = 1:M*N      %%%line 73
  U    = (S(1,t) - T(1,t)) * S1mT1;
  VL   = S(2,t) - T(2, :);
  VVLL = S(3,t) - T(3, :);
  I    = sqrt(VL .* VL + VVLL .* VVLL); 
  gij  = exp((I .* I) * (-0.5)) * c;
  GT   = gij * U' + GT;
end            %%%line 88