%COMP 236 Programming Assignment 2
%Adrian Ilie, Mar. 2002

function m=drawTriangle(img1,t)

iy13=0;%# of lines of an empty triang;e

time=zeros(1,5);
tic;

pi=[t(1) t(4) t(7) t(10) t(13) t(16);
    t(2) t(5) t(8) t(11) t(14) t(17);
    t(3) t(6) t(9) t(12) t(15) t(18)];

%construct an empty image
img=struct('Z',zeros(size(img1.Z)),'R',zeros(size(img1.R)),'G',zeros(size(img1.G)),'B',zeros(size(img1.B)));

ps=sortrows(pi,[2 1]);%sort by Y, break ties by X

p1=ps(1,:);
p2=ps(2,:);
p3=ps(3,:);

%precompute values for distances in Y
y12=p2(2)-p1(2);
y13=p3(2)-p1(2);
y23=p3(2)-p2(2);

%distances between points in screen coordinates
d12=sqrt((p2(1)-p1(1))*(p2(1)-p1(1))+y12*y12);
d13=sqrt((p3(1)-p1(1))*(p3(1)-p1(1))+y13*y13);
d23=sqrt((p3(1)-p2(1))*(p3(1)-p2(1))+y23*y23);

time(1)=time(1)+toc;

%compute area of the triangle with Heron's formula (Kahan version)
if (d13 > d12)
    tmp = d12;
    d12 = d13;
    d13 = tmp;
end
if (d23 > d12)
    tmp = d12;
    d12 = d23;
    d23 = tmp;
end
if (d23 > d13)
    tmp = d13;
    d13 = d23;
    d23 = tmp;
end
ar = sqrt((d12+(d13+d23))*(d23-(d12-d13))*(d23+(d12-d13))*(d12+(d13-d23)))/4;

%don't bother for triangles with one of the sides too close to 0
if((d12>=1)&(d13>=1)&(d23>=1))

tic;
    
%precompute values for interior distances in Y
iy12=floor(p2(2))-ceil(p1(2))+1;
iy13=floor(p3(2))-ceil(p1(2))+1;
iy23=floor(p3(2))-ceil(p2(2))+1;

%generate scanline Y cordinates
vY12=linspace(ceil(p1(2)),floor(p2(2)),iy12);
vY13=linspace(ceil(p1(2)),floor(p3(2)),iy13);
vY23=linspace(ceil(p2(2)),floor(p3(2)),iy23);

%evaluate X values from the line equations
vX12=((p2(1)-p1(1))*(vY12-p1(2)))/y12+p1(1);
vX13=((p3(1)-p1(1))*(vY13-p1(2)))/y13+p1(1);
vX23=((p3(1)-p2(1))*(vY23-p2(2)))/y23+p2(1);

%precompute values for distances in 1/Z
z12=p2(3)-p1(3);
z13=p3(3)-p1(3);
z23=p3(3)-p2(3);

%precompute distances to points on the lines form vertices
dp12=sqrt((vX12-p1(1)).*(vX12-p1(1))+(vY12-p1(2)).*(vY12-p1(2)));
dp13=sqrt((vX13-p1(1)).*(vX13-p1(1))+(vY13-p1(2)).*(vY13-p1(2)));
dp23=sqrt((vX23-p2(1)).*(vX23-p2(1))+(vY23-p2(2)).*(vY23-p2(2)));

%depth interpolation
vZ12=z12*dp12/d12+p1(3);
vZ13=z13*dp13/d13+p1(3);
vZ23=z23*dp23/d23+p2(3);

%precompute differences in values for color
cR12=p2(4)-p1(4);
cG12=p2(5)-p1(5);
cB12=p2(6)-p1(6);
cR13=p3(4)-p1(4);
cG13=p3(5)-p1(5);
cB13=p3(6)-p1(6);
cR23=p3(4)-p2(4);
cG23=p3(5)-p2(5);
cB23=p3(6)-p2(6);

%color interpolation
vR12=cR12*dp12/d12+p1(4);
vG12=cG12*dp12/d12+p1(5);
vB12=cB12*dp12/d12+p1(6);
vR13=cR13*dp13/d13+p1(4);
vG13=cG13*dp13/d13+p1(5);
vB13=cB13*dp13/d13+p1(6);
vR23=cR23*dp23/d23+p2(4);
vG23=cG23*dp23/d23+p2(5);
vB23=cB23*dp23/d23+p2(6);

%common line would get sampled twice?!
t=min(iy23,iy13-iy12);
s=iy13-t;
t=iy23-t;

line1=struct('X',vX12,'Y',vY12,'Z',vZ12,'R',vR12,'G',vG12,'B',vB12);
line2=struct('X',vX13(1:iy12),'Y',vY13(1:iy12),'Z',vZ13(1:iy12),'R',vR13(1:iy12),'G',vG13(1:iy12),'B',vB13(1:iy12));
line3=struct('X',vX13(s+1:iy13),'Y',vY13(s+1:iy13),'Z',vZ13(s+1:iy13),'R',vR13(s+1:iy13),'G',vG13(s+1:iy13),'B',vB13(s+1:iy13));
line4=struct('X',vX23(1+t:iy23),'Y',vY23(1+t:iy23),'Z',vZ23(1+t:iy23),'R',vR23(1+t:iy23),'G',vG23(1+t:iy23),'B',vB23(1+t:iy23));

%"point to left of line" test
d=det([p2(1:2) 1;p1(1:2) 1;p3(1:2) 1]);

if(d>0)%which line is to the left?
    lLeft1=line1;
    lRight1=line2;
    lLeft2=line4;
    lRight2=line3;
else
    lLeft1=line2;
    lRight1=line1;
    lLeft2=line3;
    lRight2=line4;
end    

time(1)=time(1)+toc;%Tsetup

tic;

if((length(lLeft1.X)>0)&(length(lRight1.X)>0))
for y=lLeft1.Y %for each scanline
    ind=y-lLeft1.Y(1)+1;
    x1=lLeft1.X(ind);%compute the span
    x2=lRight1.X(ind);
    sx1=ceil(x1);
    sx2=floor(x2);
    if(sx1<=sx2)%don't bother for lines with size close to zero
        sp=x2-x1;
        if(sp>0)         
            vX=[sx1:sx2];%center pixels
            img.Z(y,vX)=((lRight1.Z(ind)-lLeft1.Z(ind))*(vX-x1))/sp+lLeft1.Z(ind);%interpolate
            img.R(y,vX)=((lRight1.R(ind)-lLeft1.R(ind))*(vX-x1))/sp+lLeft1.R(ind);
            img.G(y,vX)=((lRight1.G(ind)-lLeft1.G(ind))*(vX-x1))/sp+lLeft1.G(ind);
            img.B(y,vX)=((lRight1.B(ind)-lLeft1.B(ind))*(vX-x1))/sp+lLeft1.B(ind);
        end%if sp>0
    end%lines>0
end
end

if((length(lLeft2.X)>0)&(length(lRight2.X)>0))
for y=lLeft2.Y %for each scanline
    ind=y-lLeft2.Y(1)+1;
    x1=lLeft2.X(ind);%compute the span
    x2=lRight2.X(ind);
    sx1=ceil(x1);
    sx2=floor(x2);
    if(sx1<=sx2)%don't bother for lines with size close to zero
        sp=x2-x1;
        if(sp>0)
            vX=[sx1:sx2];%center pixels
            img.Z(y,vX)=((lRight2.Z(ind)-lLeft2.Z(ind))*(vX-x1))/sp+lLeft2.Z(ind);%interpolate
            img.R(y,vX)=((lRight2.R(ind)-lLeft2.R(ind))*(vX-x1))/sp+lLeft2.R(ind);
            img.G(y,vX)=((lRight2.G(ind)-lLeft2.G(ind))*(vX-x1))/sp+lLeft2.G(ind);
            img.B(y,vX)=((lRight2.B(ind)-lLeft2.B(ind))*(vX-x1))/sp+lLeft2.B(ind);
        end%if sp>0
    end%lines>0
end
end

time(2)=toc;%/iy13;%Tline

tic;

%cap the color range
img.R(find(img.R<0))=0;
img.G(find(img.G<0))=0;
img.B(find(img.B<0))=0;
img.R(find(img.R>255))=255;
img.G(find(img.G>255))=255;
img.B(find(img.B>255))=255;

time(3)=toc;%Tcap
tic;

%merge with previous image
ind=find(img.Z>=img1.Z);

time(4)=toc;%Tfind
tic;

%merge with previous image
ind=find(img.Z>=img1.Z);
img1.Z(ind)=img.Z(ind);
img1.R(ind)=img.R(ind)/255;
img1.G(ind)=img.G(ind)/255;
img1.B(ind)=img.B(ind)/255;

time(5)=toc;

else%if size is close to zero
    tic;
    p=mean([p1' p2' p3'],2)'; %just take the mean
    pos=[round(p(2)) round(p(1))];
    if(img1.Z(pos)>p(3))      
        if(ar>1) %triangle area is more than one pixel
            img1.Z(pos)=p(3);
            img1.R(pos)=p(4)/255;
            img1.G(pos)=p(5)/255;
            img1.B(pos)=p(6)/255;
        end
    end
    time(5)=toc;%Tpixel

end%if image size is close to zero


%time(5)=time(5)/ar;%Tpixel

img1.T(1:5)=time;%times
img1.T(6)=iy13;%#of lines
img1.T(7)=ar;%#of pixels (area)

m=img1;
