fig5_19.m

these are the Mtlab functions to plot the figures in the book Communication system (by Haykin)

%plot typical motion fields due to four different types of motion models.

%to plot a figure; first use the matlab script generate the motion field
%then do "copy figure" from the figure window->edit.
%past onto corel draw using "paste special, enhanced metafile.
%then within corel draw, draw the deformed triangle based on MVs.



x=zeros(4,1);y=zeros(4,1);

dx=zeros(4,1);dy=zeros(4,1);
dxx=zeros(4,1);dyy=zeros(4,1);

dx(1)=1;dy(1)=1;
dx(2)=2;dy(2)=-2;
dx(3)=1;dy(3)=1;
dx(4)=-1;dy(4)=-0.5;

dx=3*dx;
dy=3*dy;

x(1)=-5;y(1)=-5; x(2)=5,y(2)=-5;x(3)=5;y(3)=5;x(4)=-5;y(4)=5;
[x,y,dx,dy]

%translation
u=6*ones(11,11);
v=-6*ones(11,11);
figure;
quiver(u,v,1);axis ij;
axis image; axis([-5,16,-5,16]); 
set(gca,'xtick',[]);
set(gca,'ytick',[]);


% affine transform
A=zeros(4,3);
for (i=1:4)
   A(i,1:3)=[1,x(i),y(i)];
end;

as=zeros(3,1);
bs=zeros(3,1);

%least square fit
as=pinv(A)*dx;
bs=pinv(A)*dy;

[as,bs]

for (i=1:4)
   dxx(i)=as(1)+as(2)*x(i)+as(3)*y(i);
   dyy(i)=bs(1)+bs(2)*x(i)+bs(3)*y(i);
end;

[dx,dy,dxx,dyy]

A=zeros(3,3);
for (i=1:3)
   A(i,1:3)=[1,x(i),y(i)];
end;

%fit into first 3 corner

as=zeros(3,1);
bs=zeros(3,1);
as=inv(A(1:3,1:3))*dx(1:3);
bs=inv(A(1:3,1:3))*dy(1:3);

for (i=1:4)
   dxx(i)=as(1)+as(2)*x(i)+as(3)*y(i);
   dyy(i)=bs(1)+bs(2)*x(i)+bs(3)*y(i);
end;

[dx,dy,dxx,dyy]


for (i=1:11) for (j=1:11) 
      xi=j-6;yj=i-6;
      u(i,j)=as(1)+as(2)*xi+as(3)*yj;
      v(i,j)=bs(1)+bs(2)*xi+bs(3)*yj;
      %u(i,j)=a1*i+a2*j+a0; v(i,j)=b1*i+b2*j+b0;
   end;end;
figure;
quiver(u,v,2);axis ij;
axis image; axis([-5,16,-5,16]); 
set(gca,'xtick',[]);
set(gca,'ytick',[]);

%bilinear transform

A=zeros(4,4);
for (i=1:4)
   A(i,1:4)=[1,x(i),y(i),x(i)*y(i)];
end;

as=zeros(4,1);
bs=zeros(4,1);
as=inv(A)*dx;
bs=inv(A)*dy;

for (i=1:4)
   dxx(i)=as(1)+as(2)*x(i)+as(3)*y(i)+as(4)*x(i)*y(i);
   dyy(i)=bs(1)+bs(2)*x(i)+bs(3)*y(i)+bs(4)*x(i)*y(i);
end;

[dx,dy,dxx,dyy]


for (i=1:11) for (j=1:11) 
      xi=j-6;yj=i-6;      
      u(i,j)=as(1)+as(2)*xi+as(3)*yj+as(4)*xi*yj;
      v(i,j)=bs(1)+bs(2)*xi+bs(3)*yj+bs(4)*xi*yj;
      %u(i,j)=a1*i+a2*j+a0; v(i,j)=b1*i+b2*j+b0;
   end;end;
figure;
quiver(u,v,2),axis ij;
axis image; axis([-5,16,-5,16]); 
set(gca,'xtick',[]);
set(gca,'ytick',[]);

% perspective transform

A=zeros(4,8);B=zeros(4,8);
for (i=1:4)
   A(i,1:8)=[1,x(i),y(i),0,0,0, -x(i)*(x(i)+dx(i)), -y(i)*(x(i)+dx(i))];
   B(i,1:8)=[0,0,0,1,x(i),y(i), -x(i)*(y(i)+dy(i)), -y(i)*(y(i)+dy(i))];
end;

as=zeros(8,1);
as=inv([A;B])*[x+dx;y+dy];

for (i=1:4)
   dxx(i)=(as(1)+as(2)*x(i)+as(3)*y(i))/(1+as(7)*x(i)+as(8)*y(i))-x(i);
   dyy(i)=(as(4)+as(5)*x(i)+as(6)*y(i))/(1+as(7)*x(i)+as(8)*y(i))-y(i);
end;

[dx,dy,dxx,dyy]


for (i=1:11) for (j=1:11) 
      xi=j-6;yj=i-6;
      u(i,j)=(as(1)+as(2)*xi+as(3)*yj)/(1+as(7)*xi+as(8)*yj)-xi;;
      v(i,j)=(as(4)+as(5)*xi+as(6)*yj)/(1+as(7)*xi+as(8)*yj)-yj;    
   end;end;
figure;
quiver(u,v,2);axis ij;
axis image; axis([-5,16,-5,16]); 
set(gca,'xtick',[]);
set(gca,'ytick',[]);