function AbbildungenDG

% http://geometrie.hs-nb.de/MATLAB-Programme/AbbildungenDG.m   Bearbeitungsstand: 28.10.2004
% (C) Martin Nitschke, nitschke@hs-nb.de

clc;
global AbbNr
close('all')
deg=pi/180;
a=20; b=75; c=45; d=33;
P=[0  0  b  b   b-d
   0 -a -a  c-a c-a
   0  0  0  0   0]
KS=[0         15  0  0
    (c-a)*0.6  0 15  0 
    0          0  0 15];
N=zeros(size(P));
alph=[45 45 45 30  30]*deg;
P=[P P(:,1)];
dP=diff(P')';
for k=1:length(alph)
  v=dP(:,k)/norm(dP(:,k));
  w=cos(alph(k))*[-v(2);v(1);0]+sin(alph(k))*[0;0;1];
  n=cross(v,w);
  N(:,k)=n/norm(n);
end
First=[EEE(P,N,1,2,5) EEE(P,N,2,4,5),EEE(P,N,2,3,4)]
AbbNr=0;
X=[min(P(1,:)) max(P(1,:))];
Y=[min(P(2,:)) max(P(2,:))];
Z=[0 max(First(3,:))];

AbbP( [0;0;1],[X(1);Y(1);Z(1)-45],[0;0;1],P,First,N,KS);
AbbP(-[0;1;0],[X(1);Y(1)+70;Z(1)],-[0;1;0],P,First,N,KS);
AbbP(-[1;0;0],[X(1)-50;Y(1);Z(1)],-[1;0;0],P,First,N,KS);
Abb2( [0;0;1],[X(1);Y(1);Z(1)-45],[0;0;1],-[0;1;0],[X(1);Y(1)+70;Z(1)],-[0;1;0],P,First,N,KS);
AbbP([0.5;-0.5;1],[X(1)-100;Y(1)-100;Z(1)-100],[0.5;-0.5;1],P,First,N,KS);
AbbP([0.5;-0.5;1],[X(1)-10;Y(1);Z(1)],[2;1;1],P,First,N,KS);
AbbZ(First(1:3,1)+[80;-50;90],First(1:3,1)+[0;0;-100],[0;0;1],P,First,N,KS);
AbbZ([50;-60;100],[8;-12;-90],[0;0.5;1],P,First,N,KS);
[T,PProj,FProj]=AbbZ([50;-60;100],[8;-12;-90],[-0.2;0.5;1],P,First,N,KS)

function [T,PProj,FProj]=AbbZ(A,P0,n,P,First,N,KS)
global ax AbbNr pass
AbbNr=AbbNr+1;
T=Zmtx(A,P0,n);

Name=sprintf('AbbildungenDG-3d-%d',AbbNr);
ax=[inf -inf inf -inf inf -inf];
figure('numbertitle','off','name',Name);
hold('on')
Original(P,First,KS)
[PProj,FProj,AProj]=Bild(T,P,First,N,KS);
Q=[P First];
QProj=[PProj FProj];
for k=1:size(Q,2)
  drawPolygon(@plot3,[A Q(:,k) QProj(1:3,k)],'--k','erase','none')
end
[V,B]=Rahmen(P0,n,Q,PProj,FProj,AProj);
drawPolygon(@plot3,B(:,[1:end 1]),'k');
hold('off')
axis('off')
axis('equal')
view(-24,32)
axis(ax+sqrt(eps)*[-1 1 -1 1 -1 1])
print('-deps',Name);

Name=sprintf('AbbildungenDG-2d-%d',AbbNr);
ax=[inf -inf inf -inf inf -inf];
figure('numbertitle','off','name',Name);
hold('on')
Bild([[V n] [0;0;0];[0 0 0 1]]'*T,P,First,N,KS);
drawPolygon(@plot3,[V n]'*B(:,[1:end 1]),'k');
hold('off')
axis('off')
axis('equal')
view(2)
axis(ax(1:4))
print('-deps',Name);

function AbbP(r,P0,n,P,First,N,KS)
global ax AbbNr pass
AbbNr=AbbNr+1;
T=Pmtx(r,P0,n);

Name=sprintf('AbbildungenDG-3d-%d',AbbNr);
ax=[inf -inf inf -inf inf -inf];
figure('numbertitle','off','name',Name);
hold('on')
Original(P,First,KS)
[PProj,FProj,AProj]=Bild(T,P,First,N,KS);
Q=[P First];
QProj=[PProj FProj];
for k=1:size(Q,2)
  drawPolygon(@plot3,[Q(:,k) QProj(1:3,k)],'--k','erase','none')
end
[V,B]=Rahmen(P0,n,Q,PProj,FProj,AProj);
drawPolygon(@plot3,B(:,[1:end 1]),'k');
hold('off')
axis('off')
view(-24,32)
axis(ax+sqrt(eps)*[-1 1 -1 1 -1 1])
axis('equal')
print('-deps',Name);

Name=sprintf('AbbildungenDG-2d-%d',AbbNr);
ax=[inf -inf inf -inf inf -inf];
figure('numbertitle','off','name',Name);
hold('on')
Bild([[V n] [0;0;0];[0 0 0 1]]'*T,P,First,N,KS);
drawPolygon(@plot3,[V n]'*B(:,[1:end 1]),'k');
hold('off')
axis('off')
axis('equal')
view(2)
axis(ax(1:4))
print('-deps',Name);

function Abb2(r1,P1,n1,r2,P2,n2,P,First,N,KS)
global ax AbbNr pass
AbbNr=AbbNr+1;
T1=Pmtx(r1,P1,n1);
T2=Pmtx(r2,P2,n2);

Name=sprintf('AbbildungenDG-3d-%d',AbbNr);
ax=[inf -inf inf -inf inf -inf];
figure('numbertitle','off','name',Name);
hold('on')
Original(P,First,KS)
Q=[P First];
[PProj1,FProj1,AProj1]=Bild(T1,P,First,N,KS);
[PProj2,FProj2,AProj2]=Bild(T2,P,First,N,KS);
QProj1=[PProj1 FProj1];
QProj2=[PProj2 FProj2];
for k=1:size(Q,2)
  drawPolygon(@plot3,[QProj1(1:3,k) Q(:,k) QProj2(1:3,k)],'--k','erase','none')
end
[V1,B1]=Rahmen(P1,n1,Q,PProj1,FProj1,AProj1);
[V2,B2]=Rahmen(P2,n2,Q,PProj2,FProj2,AProj2);
drawPolygon(@plot3,B1(:,[1:end 1]),'k');
drawPolygon(@plot3,B2(:,[1:end 1]),'k');
hold('off')
axis('off')
view(-24,32)
axis(ax+sqrt(eps)*[-1 1 -1 1 -1 1])
axis('equal')
print('-deps',Name);

Name=sprintf('AbbildungenDG-2d-%d',AbbNr);
ax=[inf -inf inf -inf inf -inf];
figure('numbertitle','off','name',Name);
hold('on')
U=[[V1 n1] [0;0;0];[0 0 0 1]]';
[PProj1,FProj1,AProj1]=Bild(U*T1,P,First,N,KS);
UB1=U*[B1(:,[1:end 1]);ones(1,5)];
U=[1 0 0 0;0 1 0 max(UB1(2,:));0 0 1 0;0 0 0 1]*[[V2 n2] [0;0;0];[0 0 0 1]]';
[PProj2,FProj2,AProj2]=Bild(U*T2,P,First,N,KS);
UB2=U*[B2(:,[1:end 1]);ones(1,5)];
x=unique([UB1(1,:) UB2(1,:)]);
y=unique([UB1(2,:) UB2(2,:)]);
UB=[min(x) min(x) max(x) max(x) min(x)
    min(y) max(y) max(y) min(y) min(y)
    0      0      0      0      0
    1      1      1      1      1];
drawPolygon(@plot3,UB,'k');
plot([min(x) max(x)],[1 1]*(max(PProj1(2,:))+min(PProj2(2,:)))/2,'k')
hold('off')
axis('off')
axis('equal')
view(2)
axis(ax(1:4))
print('-deps',Name);

function [V,B]=Rahmen(P0,n,Q,PProj,FProj,AProj);
QProj=[PProj FProj];
v1=PProj(:,3)-PProj(:,2);
if norm(v1)<sqrt(eps)
  v1=PProj(:,2)-PProj(:,1);
end
v1=v1/norm(v1);
v2=cross(n,v1);
V=[v1 v2/norm(v2)];
PNull=repmat(P0,1,size(Q,2)+4);
R=V'*([QProj AProj]-PNull);
Rand=0.17*min(max(R(1,:))-min(R(1,:)),max(R(2,:))-min(R(2,:)));
r=[[min(R(1,:)) max(R(1,:))]+Rand*[-1 1]
   [min(R(2,:)) max(R(2,:))]+Rand*[-1 1]];  
B=[P0+r(1,1)*V(:,1)+r(2,1)*V(:,2) ...
   P0+r(1,1)*V(:,1)+r(2,2)*V(:,2) ... 
   P0+r(1,2)*V(:,1)+r(2,2)*V(:,2) ...
   P0+r(1,2)*V(:,1)+r(2,1)*V(:,2)];

function Original(P,First,KS)
deg=pi/180;
for k=2:4
  drawPolygon(@plot3,[KS(:,1) KS(:,1)+KS(:,k)],'k')
  Spitze(KS(:,1)+KS(:,k),KS(:,k),6,30*deg)
end
drawPolygon(@fill3,[P(:,1) First(:,1)     P(:,2)],0.5*[1 1 1],'EdgeColor','none')
drawPolygon(@fill3,[P(:,2) First          P(:,3)],0.9*[1 1 1],'EdgeColor','none')
drawPolygon(@fill3,[P(:,3) First(:,3)     P(:,4)],0.5*[1 1 1],'EdgeColor','none')
drawPolygon(@fill3,[P(:,4) First(:,[3 2]) P(:,5)],0.6*[1 1 1],'EdgeColor','none')
drawPolygon(@fill3,[P(:,5) First(:,[2 1]) P(:,1)],0.7*[1 1 1],'EdgeColor','none')

function [PProj,FProj,AProj]=Bild(T,P,First,N,KS);
deg=pi/180;
AProj=T*[KS(:,1) KS(:,1)+KS(:,2) KS(:,1)+KS(:,3) KS(:,1)+KS(:,4);1 1 1 1];
AProj=[AProj(1,:)./AProj(4,:) 
       AProj(2,:)./AProj(4,:)
       AProj(3,:)./AProj(4,:)];
for k=2:4
  drawPolygon(@plot3,[AProj(:,1) AProj(:,k)],'k')
  v=AProj(:,k)-AProj(:,1);
  L=6;
  phi=30*deg;
  if all(v==0)
    r=2*L*tan(phi);
    [xK,yK,zK]=sphere(60);
    plot3(AProj(1,k)+xK,AProj(2,k)+yK,AProj(3,k)+zK,'k')
  else
    Spitze(AProj(:,k),v,L,phi);
  end
end
PProj=T*[P;ones(1,size(P,2))];
PProj=[PProj(1,:)./PProj(4,:) 
       PProj(2,:)./PProj(4,:)
       PProj(3,:)./PProj(4,:)];
FProj=T*[First;ones(1,size(First,2))];
FProj=[FProj(1,:)./FProj(4,:)
       FProj(2,:)./FProj(4,:)
       FProj(3,:)./FProj(4,:)];
S={[PProj(:,1) FProj(:,1)     PProj(:,2)]
   [PProj(:,2) FProj          PProj(:,3)]
   [PProj(:,3) FProj(:,3)     PProj(:,4)]
   [PProj(:,4) FProj(:,[3 2]) PProj(:,5)]
   [PProj(:,5) FProj(:,[2 1]) PProj(:,1)]};
c=[0.5 0.9 0.5 0.6 0.7 0.9];
for s=1:length(S)
  drawPolygon(@fill3,S{s},c(s)*[1 1 1],'EdgeColor','none')
end

function A=EEE(P,N,i,j,k);
A=[N(:,i)';N(:,j)';N(:,k)']\[N(:,i)'*P(:,i);N(:,j)'*P(:,j);N(:,k)'*P(:,k)];

function T=Projmtx(z,u)
T=z*u'-u'*z*eye(4);

function T=Zmtx(A,P0,n) % Abbildungsmatrix Zentralprojektion mit Augenpunkt A, Bildebene (P-P0)'*n=0
T=Projmtx([A;1],[n;-P0'*n]);

function T=Pmtx(r,P0,n) % Abbildungsmatrix Parallelprojektion mit Blickrichtung r, Bildebene (P-P0)'*n=0
n=n/norm(n);
T=Projmtx([r;0],[n;-P0'*n]);

function drawPolygon(draw,P,varargin);
global ax
feval(draw,P(1,:),P(2,:),P(3,:),varargin{:})
if min(P(1,:))<ax(1)
  ax(1)=min(P(1,:));
end
if max(P(1,:))>ax(2)
  ax(2)=max(P(1,:));
end
if min(P(2,:))<ax(3)
  ax(3)=min(P(2,:));
end
if max(P(2,:))>ax(4)
  ax(4)=max(P(2,:));
end
if min(P(3,:))<ax(5)
  ax(5)=min(P(3,:));
end
if max(P(3,:))>ax(6)
  ax(6)=max(P(3,:));
end

function Spitze(P,v,L,phi,varargin);
deg=pi/180;
if nargin<=2
  L=1;
end
if nargin<=3
  phi=7.5*deg;
end
v=L*v./repmat(sqrt(dot(v,v)),size(P,1),1);
Q=P-v;
t=[0:6:360]*deg;
for k=1:size(P,2)
  N=tan(phi/2)*null(v(:,k)');
  for w=2:length(t)
    plot3(...
      [P(1,k) Q(1,k)+cos(t(w))*L(:,k)*N(1,1)+sin(t(w))*L(:,k)*N(1,2) Q(1,k)+cos(t(w-1))*L(:,k)*N(1,1)+sin(t(w-1))*L(:,k)*N(1,2) P(1,k)],...
      [P(2,k) Q(2,k)+cos(t(w))*L(:,k)*N(2,1)+sin(t(w))*L(:,k)*N(2,2) Q(2,k)+cos(t(w-1))*L(:,k)*N(2,1)+sin(t(w-1))*L(:,k)*N(2,2) P(2,k)],...
      [P(3,k) Q(3,k)+cos(t(w))*L(:,k)*N(3,1)+sin(t(w))*L(:,k)*N(3,2) Q(3,k)+cos(t(w-1))*L(:,k)*N(3,1)+sin(t(w-1))*L(:,k)*N(3,2) P(3,k)],...
      'k',varargin{:});
  end
end