! Fraser Millar (FGM) - CS4
! Graphics Practical 2
! 20/2/87

%include "inc:maths.imp"
%include "iff:graphinc.imp"
!!%include "level1:graphinc.imp"

%begin

%string (255) param
%const %byte colon=58,dollar=36,hash=35,comma=44
%const %integer red=1,green=32,blue=32*32
%const %integer ambient=8
%real %array x(1:500),y(1:500),z(1:500)
%real %array normal (1:500,1:3)   {point normals}
%integer %array dx(1:500),dy(1:500),bright(1:500)
%real %array polynorm (1:500,1:3)  {polygon normals}
%real vx,vy,vz,tx,ty,tz,avz1,avz2,d,ztemp,modc,modn
%real %array xrot(0:3,0:3),yrot(0:3,0:3),zrot(0:3,0:3),pers(0:3,0:3)
%real %array a(1:3),b(1:3),c(1:3),no(1:3),sum(1:3)
%integer %array pt(1:500,1:4)
%byte ch
%integer i,npoints,n,p,npolys,dummy,d1,d2,gl,pc,polys,g2,g4,y2,y4,xa,ya
%integer divisions
%half %integer %name cmp
%half %integer %array cm(0:255)

%routine drawpoly(%integer p,color,xoffset,yoffset)
  %integer %array sx(1:4),sy(1:4),inten(1:4)
  %integer di,dj,dk,dx1,dx2,dy1,i1,i2
  %real ex1,ey1,ex2,ey2
  sx(di)=dx(pt(p,di))+xoffset %for di=1,1,4
  sy(di)=dy(pt(p,di))+yoffset %for di=1,1,4
  inten(di)=bright(pt(p,di)) %for di=1,1,4
  {Will try to divide poly into a divisions*divisions area}
  %for di=0,1,divisions-1 %cycle
    ex1=(sx(1)+di*(sx(2)-sx(1))/divisions)
    ey1=(sy(1)+di*(sy(2)-sy(1))/divisions)
    i1=inten(1)+(2*di+1)*(inten(2)-inten(1))//(divisions*2)
    ex2=(sx(4)+di*(sx(3)-sx(4))/divisions)
    ey2=(sy(4)+di*(sy(3)-sy(4))/divisions)
    i2=inten(4)+(di*2+1)*(inten(3)-inten(4))//(divisions*2)
    %for dj=0,1,divisions-1 %cycle
      dx1=int(ex1+dj*(ex2-ex1)/divisions)
      dy1=int(ey1+dj*(ey2-ey1)/divisions)
      dk=i1+(dj*2+1)*(i2-i1)//(divisions*2)
      colour(color*32+dk)
      poly(dx1,dy1)
      poly(dx1+int((sx(2)-sx(1))/divisions),dy1+int((sy(2)-sy(1))/divisions))
      poly(dx1+int((sx(3)-sx(1))/divisions),dy1+int((sy(3)-sy(1))/divisions))
      poly(dx1+int((sx(4)-sx(1))/divisions),dy1+int((sy(4)-sy(1))/divisions))
      close poly
    %repeat
  %repeat
%end

%routine transform (%real %array %name m (0:3,0:3),%real %array %name x (0:3))
! applies the transformation denoted by M (homogeneous) to point x
! Is actually pre-multiplication of vector x by matrix m.         
  %real %array r(0:3)
  %for d1=0,1,3 %cycle
    r(d1)=0
    %for d2=0,1,3 %cycle
      r(d1)=r(d1)+m(d1,d2)*x(d2)
    %repeat
  %repeat
  %for d1=0,1,3 %cycle
    x(d1)=r(d1)
  %repeat
%end

%routine transpt (%real %array %name mat (0:3,0:3),%integer p)
! calls TRANSFORM to transform point p by matrix m 
  %real %array xp(0:3)
  xp(0)=x(p)
  xp(1)=y(p)
  xp(2)=z(p)
  xp(3)=1
  transform (mat,xp)
  z(p)=xp(2)/xp(3)
  y(p)=xp(1)/xp(3)
  x(p)=xp(0)/xp(3)
%end

!Set up color map for primary coplor scales

%for i=0,1,255 %cycle
  gl=i&31
  pc=i//32
  cm(i)=0
  %if pc&1=1 %then cm(i)=cm(i)+gl*red
  %if pc&2=2 %then cm(i)=cm(i)+gl*green
  %if pc&4=4 %then cm(i)=cm(i)+gl*blue
%repeat
param=cli param
p=iff open frame("torus","torus",3)
cmp==cm(0)
update colour map (cmp)

! Actually starts here folks!
! First of all, open the file with the data.

open input (1,param)
print string ("Reading file ".param)
newline
select input (1)

i=1
%cycle
  %cycle          
    readsymbol(ch)
  %repeat %until ch=colon    {skip garbage until a colon found}

  read(x(i))
  read(y(i))
  read(z(i))
  i=i+1
  skipsymbol  {should be CR}
%repeat %until nextsymbol=dollar

npoints=i-1
write(npoints,5)
print string (" points read.")
newline

i=1
%cycle
  %cycle
    readsymbol(ch)
  %repeat %until ch=colon              {skip garbage until colon found}

  read(n)    
  %if n<>4 %then %start
                 printstring ("Cockup in polygon")
                 write(i,4)
                 printstring (" - not a quadrilateral")
                 newline
                 %finish
  %for p=1,1,4 %cycle
    %cycle
      readsymbol(ch)
    %repeat %until ch=84   {capital T}
    read(pt(i,p))
  %repeat
  i=i+1
  skipsymbol     {should be CR again}
%repeat %until nextsymbol=dollar
                 
npolys=i-1
write(npolys,5)  
printstring(" polygons read.")
newline

%cycle
  readsymbol (ch)
%repeat %until ch=hash

read(dummy)              {should be a 1}

read(vx)
read(vy)
read(vz)                 {viewpoint coordinates}
vz=vz*4                  {modify vz to avoid excessive image distortion}

read(tx)
read(ty)
read(tz)                 {view angles}

{**** That's the stuff read in, now the nitty gritty... ****}

! Define transformation matrices...
! Will keep them separate for ease of programming, rather
! than composing them into a single transformation matrix

print string ("Defining Transformations")
newline

ty=ty/DtoR; tx=tx/DtoR; tz=tz/DtoR     {convert angles to radians}

xrot(0,0)=1;        xrot(0,1)=0;        xrot(0,2)=0;        xrot(0,3)=0
xrot(1,0)=0;        xrot(1,1)=cos(tx);  xrot(1,2)=-sin(tx); xrot(1,3)=0
xrot(2,0)=0;        xrot(2,1)=sin(tx);  xrot(2,2)=cos(tx);  xrot(2,3)=0
xrot(3,0)=0;        xrot(3,1)=0;        xrot(3,2)=0;        xrot(3,3)=1

yrot(0,0)=cos(ty);  yrot(0,1)=0;        yrot(0,2)=sin(ty);  yrot(0,3)=0
yrot(1,0)=0;        yrot(1,1)=1;        yrot(1,2)=0;        yrot(1,3)=0
yrot(2,0)=-sin(ty); yrot(2,1)=0;        yrot(2,2)=cos(ty);  yrot(2,3)=0
yrot(3,0)=0;        yrot(3,1)=0;        yrot(3,2)=0;        yrot(3,3)=1

zrot(0,0)=cos(tz);  zrot(0,1)=-sin(tz); zrot(0,2)=0;        zrot(0,3)=0
zrot(1,0)=sin(tz);  zrot(1,1)=cos(tz);  zrot(1,2)=0;        zrot(1,3)=0
zrot(2,0)=0;        zrot(2,1)=0;        zrot(2,2)=1;        zrot(2,3)=0
zrot(3,0)=0;        zrot(3,1)=0;        zrot(3,2)=0;        zrot(3,3)=1

%for d1=0,1,3 %cycle     {perspective transformation is mostly zeros  }
  %for d2=0,1,3 %cycle   {so easier to define this way...             }
    pers(d1,d2)=0
    pers(d2,d2)=1
  %repeat
%repeat
pers(2,2)=0    
pers(3,2)=1/vz

print string ("Transforming")
newline

%for i=1,1,npoints %cycle     {apply transformations in order to each point}
  transpt (xrot,i)
  transpt (yrot,i)
  transpt (zrot,i)
  ztemp=z(i)
  transpt (pers,i)     {but need to retain z for depth sorting or normal calculation}
  z(i)=ztemp           {since z(i)=0 after perspective projection}
%repeat

{map real points onto 2D integer points}
%for i=1,1,npoints %cycle
  dx(i)=int((x(i)*6+vx)/2)         {5 is arbitrary scale factor to make the}
  dy(i)=int((y(i)*6+vy)/2)         {image a reasonable size}
%repeat

clear
print string ("Calculating Surface Normals")
newline
%for p=1,1,npolys %cycle
! Calculate surface normals

  a(1)=x(pt(p,3))-x(pt(p,1))
  a(2)=y(pt(p,3))-y(pt(p,1))
  a(3)=z(pt(p,3))-z(pt(p,1))        {ie diagonal a=P3-P1}
  
  b(1)=x(pt(p,4))-x(pt(p,2))
  b(2)=y(pt(p,4))-y(pt(p,2))
  b(3)=z(pt(p,4))-z(pt(p,2))        {ie diagonal b=P4-P2}
  
  no(1)=a(2)*b(3)-a(3)*b(2)
  no(2)=a(3)*b(1)-a(1)*b(3)
  no(3)=a(1)*b(2)-a(2)*b(1)         {ie normal=A x B}

  polynorm(p,i)=no(i) %for i=1,1,3

%repeat
  
!Now calculate point normals
print string ("Calculating Vertex normals")
newline
%for p=1,1,npoints %cycle
  {need to average intensities of 4 polygons around this point}
  {so find them....}
  polys=0
  sum(i)=0 %for i=1,1,3
  %for i=1,1,npolys %cycle
    %for d1=1,1,4 %cycle
      %if pt(i,d1)=p %then %start
                           polys=polys+1
                           sum(d2)=sum(d2)+polynorm(i,d2) %for d2=1,1,3
                           %finish
    %repeat
  %repeat
  normal(p,i)=sum(i)/polys %for i=1,1,3
  c(1)=75
  c(2)=75
  c(3)=-75
  modc=sqrt(c(1)*c(1)+c(2)*c(2)+c(3)*c(3))
  modn=sqrt(normal(p,1)*normal(p,1)+normal(p,2)*normal(p,2)+normal(p,3)*normal(p,3))
  d=(c(1)*normal(p,1)+c(2)*normal(p,2)+c(3)*normal(p,3))/(modc*modn)
  %if d<0 %then bright(p)=ambient %c
          %else bright(p)=ambient+int((31-ambient)*d)
%repeat

xa=2; ya=1
%cycle
!Now can draw polygons...
  print string ("How many polygon divisions ? (1..16ish) ")
  select input (0)
  read (divisions)
  print string ("Drawing...")
  
  %if divisions<>0 %then %start
      %for p=1,1,npolys %cycle
write(p, 3); newline
        d=polynorm(p,3)*(z(p)+vz)   {project to (0,0,vz)}
        %if d<0 %then %start
!!                      %for xa=0,1,3 %cycle
!!                        %for ya=0,1,1 %cycle
                          drawpoly(p,xa+4*ya,xa*160-80,ya*160+50)
!!                        %repeat
!!                      %repeat
                      %finish
      %repeat
      %finish

  print string ("complete.")
  newline
%repeat %until divisions=0
newline
iff close frame
%endofprogram