function tdata=remove_planar(data)

%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%Remove planar trend of an interferogram
%The planar is defined as:
%          Trend = a0 + a1 * X + a2 * Y
%Parameters are regressed using all available phases by the above equation
%and removed from the phases.
%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


index=find(data==0);

num=0; sumx=0; sumy=0; sumxy=0; sumx2=0;sumy2=0;sumsx=0; sumsy=0; sums=0;

[len,wid]=size(data);

k=1;
for i=1:len
    for j=1:wid   
      if(data(i,j)==0  || isnan(data(i,j)))
          continue;
      end
      gridx=j;
      gridy=i;
      phase=data(i,j);
      sumx = sumx + gridx; 
      sumy = sumy + gridy;
      sumx2 = sumx2 + gridx * gridx; 
      sumy2 = sumy2 + gridy * gridy;
	  sumxy = sumxy + gridx * gridy;
	  sumsx = sumsx + phase * gridx; 
      sumsy = sumsy + phase * gridy;
	  sums = sums + phase;
      num=num+1;
    end
end

A(1,1) = sumx2; A(1,2)= sumxy;  A(1,3) = sumx;
A(2,1) = sumxy; A(2,2) = sumy2; A(2,3) = sumy;
A(3,1) = sumx;  A(3,2) = sumy;  A(3,3) = num;

L(1) = sumsx; L(2) = sumsy; L(3) = sums;
L=L';

R=A\L;

tdata=zeros(len,wid);

for i=1:len
    for j=1:wid
      if(data(i,j)==0)
          continue;
      end
      gridx=j;
      gridy=i;
      tdata(i,j)=data(i,j)-R(1)*gridx - R(2) * gridy - R(3);     
    end
end

tdata(index)=0;




