function i=reconstr(D)
%
% This reconstruction algorithm use the full positive d of the dipole
% function
%
[row col level level]=size(D);
row=(row+1)/2;
col=(col+1)/2;
s=zeros(row,col);
i=zeros(row,col);
for m=0:row-1
   for n=0:col-1
      for p=1:level
	 for q=1:level
	    s(m+1,n+1)=s(m+1,n+1)+D(2*row-1-m,2*col-1-n,p,q)*(p-1);
         end;
      end;
   end;
end;

i(1,1)=s(1,1);
for n=2:col
   i(1,n)=s(1,n)-s(1,n-1);
 end;
 for m=2:row
   i(m,1)=s(m,1)-s(m-1,1);
 end;
 for m=2:row
   for n=2:col
     i(m,n)=s(m,n)-s(m-1,n)-s(m,n-1)+s(m-1,n-1);
   end;
 end;
 
 
 