codesamples/Matlab/SphereCam.m

function [in] = SphereCam()
%SphereCam Render & look at a sphere in space
%   In order to create each sphere, takes in the following user input:
%   (vectors should be entered [x;y;z])
%   Viewpoint (the point from which rays are cast - vector)
%   (note: the z-value should be negative; origin - z taken for focal
%   length)
%   Resolution (the image is this many pixels square - scalar)
%   Size (the pixels fit into this distance - scalar)
%   Light direction vector (direction towards the light source at infinity)
%   Center of the sphere (vector)
%   Radius of the sphere (scalar)
%   Albedo of the sphere (how dark or light the surface is - scalar)
%   Gloss value (the proportion of non-Lambertian effect - scalar)
%   Gloss sharpness (closeness to mirror finish - scalar)
%   Ambient light level (how bright the background is - scalar)
%   
%   Then it renders the space specified, including the sphere if it is 
%   in the field of view. The camera can be moved by additional user input
%%%%%%%%%********* check it
%   Track: enter a translation vector [x;y;z]. 
%       x>0 moves view right; x<0 moves view left
%       y>0 moves view up; y<0 moves view down
%       z>0 moves towards sphere; z<0 moves away from sphere
%       [it zooms]
%   Pan: enter an angle in degrees
%       angle > 0 pans right; angle < 0 pans left
%       (pans right == sphere shifts left, and vice versa)
%   Tilt: enter an angle in degrees
%       angle > 0 tilts up; angle < 0 tilts down
%   	(tilts up == sphere shifts down, and vice versa)
%   
%   Once the user is done playing with the camera, they can make a new 
%   sphere or simply exit.

in = input('New sphere? ', 's');
while(in ~= 'n')
    % massive amounts of taking-in of input
    viewpoint = input('Viewpoint: ');
    viewdir = Normalize(viewpoint);
    res = input('View resolution: ');
    size = input('View size: ');
    % focal length is just the distance from the origin on z of view
    focallen = 0 - viewpoint(3);
    retina = MakeRetina(res, size, focallen);
    lightvec = input('Light direction vector: ');
    % light direction is the light vector normalized
    lightdir = Normalize(lightvec);
    center = input('Center of the sphere: ');
    Sr = input('Radius of the sphere: ');
    % I don't know how to make stripes
    albedo = input('Albedo of the sphere (unpatterned): ');
    gloss = input('Gloss value of sphere (0 to 1): ');
    glossSharp = input('Gloss sharpness of sphere: ');
    ambient = input('Ambient light level: ');
    paint = PaintSphere(viewpoint, retina);
    imshow(mat2gray(paint));
    in2 = input('Want to move the camera? ', 's');
    % any input except n will move this forward
    while(in2 ~= 'n')
        mvmt = input('Track (tr), Pan (p), or Tilt (t)? ', 's');
        switch(mvmt)
            case 'tr'
                trans = input('Translation vector: ');
                retina = Track(retina, trans);
                paint = PaintSphere(viewpoint, retina);
                imshow(mat2gray(paint));
                in2 = input('Again? ','s');
            case 'p'
                pan = input('Angle to pan: ');
                retina = Pan(retina, pan);
                paint = PaintSphere(viewpoint, retina);
                imshow(mat2gray(paint));
                in2 = input('Again? ','s');
            case 't'
                tilt = input('Angle to tilt: ');
                retina = Tilt(retina, tilt);
                paint = PaintSphere(viewpoint, retina);
                imshow(mat2gray(paint));
                in2 = input('Again? ','s');
	    case 'n'
		in2 = 'n';
            otherwise
                in2 = input('Sorry, invalid choice. Try again? ', 's');
        end
    end
    in = input('Nother sphere? ', 's');
end

    % the 'retina' is the frame on which the sphere is projected
    function retina = MakeRetina(n, mm, focallen)
        retina = zeros(3,n^2+1);
        k = mm/n; % the distance between pixels
        for l = 1:n^2
	    % because it's basically a 1-D array of 3-vectors, 
	    % instead of a 2-D one, figuring out what row a 
	    % pixel actually belongs to is trickier
            if mod(l, n) == 0
                retina(1, l) = n*k - mm/2; 
            else
                retina(1, l) = mod(l, n) * k - mm/2;
            end
            retina(2, l) = ceil(l/n) * k - mm/2;
            retina(3, l) = focallen;
        end
        retina(1:3, n^2+1) = viewpoint;
        %retina(2, n^2+1) = viewpoint(2);
        %retina(3, n^2+1) = viewpoint(3);
	% the final bit is the viewpoint, to make camera movement work
    end

    % translate the camera in space 
    function newretina = Track(retina, trans)
        n = length(retina);
        newretina = zeros(3, n);
        for l = 1:n
            newretina(1:3,l) = retina(1:3,l) + trans;
            %newretina(2,l) = retina(2,l) + trans(2);
            %newretina(3,l) = retina(3,l) + trans(3);
        end
        viewpoint = newretina(1:3, n);
        %viewpoint(2) = newretina(2, n);
        %viewpoint(3) = newretina(3, n);
    end

    % look from left to right
    function newretina = Pan(retina, angle)
        RotY = [cosd(angle),0,sind(angle);0,1,0;-sind(angle),0,cosd(angle)];
        n = length(retina);
        newretina = zeros(3, n);
        for l = 1:n
            newretina(1:3,l) = retina(1:3,l)-viewpoint;
            newretina(1:3,l) = RotY*newretina(1:3,l);
            newretina(1:3,l) = newretina(1:3,l)+viewpoint;
%             retpt = zeros(3, 1);
%             retpt = retina(1:3,l);
%             %retpt(2) = retina(2,l);
%             %retpt(3) = retina(3,l);
%             newpt = RotY*(retpt - viewpoint) + viewpoint;
%             newretina(1,l) = newpt(1);
%             newretina(2,l) = newpt(2);
%             newretina(3,l) = newpt(3);
        end
        viewpoint = newretina(1:3, n);
        %viewpoint(2) = newretina(2, n);
        %viewpoint(3) = newretina(3, n);
    end

    % look up and down
    function newretina = Tilt(retina, angle)
        RotX = [1,0,0;0,cosd(angle),-sind(angle);0,sind(angle),cosd(angle)];
        n = length(retina);
        newretina = zeros(3, n);
        for l = 1:n
            newretina(1:3,l) = retina(1:3,l)-viewpoint;
            newretina(1:3,l) = RotX*newretina(1:3,l);
            newretina(1:3,l) = newretina(1:3,l)+viewpoint;
%             retpt = zeros(3,1);
%             retpt(1) = retina(1,l);
%             retpt(2) = retina(2,l);
%             retpt(3) = retina(3,l);
%             newpt = RotX*(retpt - viewpoint) + viewpoint;
%             newretina(1,l) = newpt(1);
%             newretina(2,l) = newpt(2);
%             newretina(3,l) = newpt(3);
        end
        viewpoint = newretina(1:3, n);
        %viewpoint(2) = newretina(2, n);
        %viewpoint(3) = newretina(3, n);
    end

    % Make a vector into a unit vector
    function norm = Normalize(point)
        w = sqrt(point(1)^2 + point(2)^2 + point(3)^2);
        norm = [point(1)/w; point(2)/w; point(3)/w];
    end    

    % Paint the sphere onto the retina
    function paintmat = PaintSphere(viewpoint, retina)
        
	% find the surface normal vector of a point
	% normalize the point if the circle were centered on the origin
        function normal = SurfNorm(point)
            normal = Normalize(point - center);
        end
        
	% dot product of two vectors/points
        function dot = MyDot(point1, point2)
            dot = point1(1)*point2(2)+point1(2)*point2(2)+point1(3)*point2(3);
        end
        
	% the brightness at a given point on the sphere
        function bright = Brightness(point, intersect)
            if intersect > 0 % the point is on the sphere
                n = SurfNorm(point);
                if MyDot(n,lightdir) >= 0 % if I'm not self-shadowed
		    % equation courtesy of the assignment materials
                    C = 2*MyDot(n,lightdir)*MyDot(n,viewdir) - MyDot(viewdir,lightdir);
                    bright = albedo*(gloss*C^glossSharp + (1-gloss)*MyDot(lightdir,n)) + ambient/2;
                    if bright < 0 % if somehow the brightness is negative
                        bright = 0+ambient/2;
                    end
                else %self-shadowed brightness det. by ambient light
                    bright = 0 + ambient/2;
                end
            else
                bright = ambient;
            end
        end
        
	% casting the ray of light to see if I hit the sphere
	% returns the point and the intersect value
        function [pointintersect] = Pixel2Point(pixel)
            x1 = viewpoint(1);
            y1 = viewpoint(2);
            z1 = viewpoint(3);
            x2 = pixel(1);
            y2 = pixel(2);
            z2 = pixel(3);
            x3 = center(1);
            y3 = center(2);
            z3 = center(3);
	    % complicated equation from Bourke
	    % to be plugged into au^2 + bu + c
            a = (x2-x1)^2 + (y2-y1)^2 + (z2-z1)^2;
            b = 2*((x2-x1)*(x1-x3)+(y2-y1)*(y1-y3)+(z2-z1)*(z1-z3));
            c = x3^2+y3^2+z3^2+x1^2+y1^2+z1^2-2*(x3*x1+y3*y1+z3*z1) - Sr^2;
            intersect = b^2 - 4*a*c;
            if intersect < 0 % then I don't hit the sphere
                dist = center(3);
                point = pixel;
                point(3) = point(3)+dist;
            else
                if intersect == 0 % I'm on the edge
                    u = -b/(2*a);
                    point = viewpoint + u*(pixel - viewpoint);
                else % intersect at two places
                    u1 = (-b + sqrt(intersect))/(2*a);
                    u2 = (-b - sqrt(intersect))/(2*a);
                    u = max(u1, u2); % pick the larger intersect value
                    point = viewpoint + u*(pixel - viewpoint);
                end
            end
            pointintersect = [point(1);point(2);point(3);intersect];
        end
        
        paintmat = zeros(res, res);
        pixel = zeros(3,1);
        for i = 1:res
            for j = 1:res
		% for each pixel, find its intersect & brightness
                pixel(1) = retina(1,j+(i-1)*res);
                pixel(2) = retina(2,j+(i-1)*res);
                pixel(3) = retina(3,j+(i-1)*res);
                px2pt = Pixel2Point(pixel);
                point = px2pt(1:3);
                intersect = px2pt(4);
                paintmat(i,j) = Brightness(point, intersect);
            end
        end
        
        
    end



end