function I=reconstrHalf(D,M,N)
%
% This reconstruction algorithm use the full positive d of the dipole
% function
%
xmiddle=floor(M/2);
ymiddle=floor(N/2);
[FM FN L L]=size(D);
ddx=(FM-1)/2;
ddy=(FN-1)/2;
s=zeros(M,N);
I=zeros(M,N);

for m=xmiddle-1:M-1
   for n=ymiddle-1:N-1
      for a=1:L
         for b=1:L
            s(m+1,n+1)=s(m+1,n+1)+D(ddx+M-m,ddy+N-n,a,b)*(a-1);
         end;
      end;
   end;
end;
   
for m=xmiddle:M-1
   for n=ymiddle:N-1
      I(m+1,n+1)=s(m+1,n+1)-s(m,n+1)-s(m+1,n)+s(m,n);
   end;
 end;
 
 s=zeros(M,N);

for m=0:xmiddle
   for n=0:ymiddle
      for a=1:L
	 for b=1:L
            s(m+1,n+1)=s(m+1,n+1)+D(-m+ddx+1,-n+ddy+1,a,b)*(a-1);
	 end;
      end;
   end;
end;

for m=0:xmiddle-1
  for n=0:ymiddle-1
     I(m+1,n+1)=s(m+1,n+1)-s(m+2,n+1)-s(m+1,n+2)+s(m+2,n+2);
  end;
end;
 
s=zeros(M,N);

for m=0:xmiddle
   for n=ymiddle-1:N-1
      for a=1:L
	 for b=1:L
	   s(m+1,n+1)=s(m+1,n+1)+D(-m+ddx+1,N+ddy-n,a,b)*(a-1);
	 end;
      end;
   end;
end;
   
for m=0:xmiddle-1
   for n=ymiddle:N-1
      I(m+1,n+1)=s(m+1,n+1)-s(m+1,n)-s(m+2,n+1)+s(m+2,n);
   end;
end;
   
s=zeros(M,N);

for m=xmiddle-1:M-1
   for n=0:ymiddle
      for a=1:L
	 for b=1:L
            s(m+1,n+1)=s(m+1,n+1)+D(ddx+M-m,-n+ddy+1,a,b)*(a-1);
         end;
      end;
   end;
end;
   
for m=xmiddle:M-1
   for n=0:ymiddle-1
      I(m+1,n+1)=s(m+1,n+1)-s(m,n+1)-s(m+1,n+2)+s(m,n+2);
   end;
end;












