%COMP 236 Programming Assignment 1
%Adrian Ilie, Feb. 2002

function pic=cerealBox(cbp,imSize,cop,angle,focalLength)
s=size(cbp);
numPoints=s(1);
imWidth=imSize(2);
imHeight=imSize(1);
Pwi=[1 0 0;0 1 0;0 0 focalLength];%projection matrix
Piw=inv(rotation(Pwi,angle));%inverse of matrix
cop=rotation(cop',angle);%center of projection in pixels
cbpt=cbp-(cop*ones(1,numPoints))';%subtract cop
coords=(Piw*cbpt')';%project onto image plane
ncoords=[coords(:,1)./coords(:,3) coords(:,2)./coords(:,3) ones(1,numPoints)'];%normalize

x = 1 + ncoords(:,1) + imWidth/2;
clamp = find(x>imWidth);
x(clamp) = imWidth;
clamp = find(x<1);
x(clamp) = 1;

y = 1 + ncoords(:,2) + imHeight/2;
clamp = find(y>imHeight);
y(clamp) = imHeight;
clamp = find(y<1);
y(clamp) = 1;

pic=zeros(imHeight,imWidth);%picture matrix
%this would take care of it if the indices were integers
%ind=sub2ind([imHeight imWidth],y,x);
%pic(ind)=1;

%we do anti-aliasing instead
maxd=1.5*1.5*10;%max value distance can have * 10
g=fspecial('gaussian', [1,2*maxd], 10);
maxf=1/g(round(maxd));

for i=1:numPoints
    px=round(x(i));
    py=round(y(i));
    for fx=px-1:px+1
        for fy=py-1:py+1
            dx=x(i)-fx;
            dy=y(i)-fy;
            dist=sqrt(dx*dx+dy*dy);
            if((fx>2)&(fx<imWidth-1)&(fy>2)&(fy<imHeight-1))
                pic(fy,fx)=pic(fy,fx)+g(round(maxd+10*dist))*maxf;%they call this "additive blending" ;)
            end
        end
    end    
end

i=find(pic>1);
pic(i)=1;%cap the values

%imshow(pic);










% rotation about x-axis pointing straight up
% angle given in degrees
function ret=rotation(m, angle)

rot= eye(4);

cosA=cos(angle(1)*pi/180);
sinA=sin(angle(1)*pi/180);
rotX=[1 0 0 0;0 cosA -sinA 0;0 sinA cosA 0;0 0 0 1];
rot = rot * rotX;

cosA=cos(angle(2)*pi/180);
sinA=sin(angle(2)*pi/180);
rotY=[cosA 0 sinA 0;0 1 0 0;-sinA 0 cosA 0;0 0 0 1];
rot=rot * rotY;

cosA=cos(angle(3)*pi/180);
sinA=sin(angle(3)*pi/180);
rotZ=[cosA -sinA 0 0;sinA cosA 0 0;0 0 1 0;0 0 0 1];
rot=rot * rotZ;

mo=[[m; ones(1,size(m,2))] ones(size(m,1)+1,1)];
ret = rot*mo;
ret=ret(1:(size(ret,1)-1),1:(size(ret,2)-1));

