Computer Graphics                                     Solutions to Homework 2
CS417/CS418 Fall 1998                        Due Wednesday, February 17, 1999

1.

function [pixel] = drawspheres(c,r,ka,kd,ks,n,ambient,l,lc)
% [pixel] = drawspheres(c,r,ka,kd,ks,n,l,lc):
%     return an image pixel (with integer world and view coordinates)
%     just large enough to show all the spheres under the specified lights.
%
%     the i-th sphere has center c(1:3,i), radius r(i),
%     ambient coefficient ka(i), diffuse coefficients kd(1:3,i),
%     specular coefficient ks(i), specular exponent n(i)
%     (the ambient coefficients for each r,g,b component are ka(i)*kd(1:3,i)).
%
%     the ambient light has color ambient(1:3)
%
%     the j-th light comes from direction l(1:3,j), with color lc(1:3,j)

y = c(2,:); y = [y-r y+r]; ybias = 1-min(y); % offset to make y coords >= 1
x = c(1,:); x = [x-r x+r]; xbias = 1-min(x); % offset to make x coords >= 1
high = ceil(max(y))-floor(min(y))+1;
wide = ceil(max(x))-floor(min(x))+1;

pixel = zeros(high,wide,3)-1; % use -1 to flag background
zbuffer = zeros(high,wide) - inf; % z-buffer

v = [0 ; 0 ; 1]; % view vector (towards eye)

% normalize and reverse light vector directions
for j = 1:size(l,2), l(:,j) = - l(:,j)/norm(l(:,j)); end

% render each sphere, doing hidden surface removal using zbuffer
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% straightforward but slow way:
% for i = 1:size(c,2)
%   for x = ceil(c(1,i) - r) : floor(c(1,i) + r)
%     for y = ceil(c(2,i) - r) : floor(c(2,i) + r)
%        z = r(i)^2 - (x-c(1,i))^2 - (y-c(2,i))^2;
%        if z>=0
%            z = sqrt(z) + c(3,i);
%            if z>zbuffer(y+ybias,x+xbias)
%                zbuffer(y+ybias,x+xbias) = z;
%                color = kd(:,i)*ka(i).*ambient;
%                normal = [x ; y ; z] - c(:,i);
%                normal = normal/norm(normal);
%                for j = 1:size(l,2)
%                    dot = normal' * l(:,j);
%                    if dot>=0
%                        color = color + dot * kd(:,i) .* lc(:,j);
%                        reflect = 2*dot*normal - l(:,j);
%                        dot = reflect' * v;
%                        if dot >= 0
%                            color = color + dot^n(i) * ks(i) * lc(:,j);
%                        end
%                    end
%                end % for j
%                pixel(y+ybias,x+xbias,:) = color;
%            end % if z>zbuffer
%        end % if z>=0
%     end % for y
%   end % for x
% end % for i
?
(page 2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fancier, faster way.
% problem: matlab doesn't allow logical indices into sub-arrays.
% therefore, we must keep switching between 1-d and 2-d arrays
% and 1-d indices and 2-d subscripts (yuck!)
for i = 1:size(c,2)
    x = ceil(c(1,i) - r(i)) : floor(c(1,i) + r(i));
    y = ceil(c(2,i) - r(i)) : floor(c(2,i) + r(i));
    x = ones(length(y),1) * x;
    y = y' * ones(1,size(x,2));
    x = x(:); y = y(:); % reshape to 1-d
    z = r(i).^2 - (x-c(1,i)).^2 - (y-c(2,i)).^2;
    % identify pixels where sphere is seen and is closer; update zbuffer
          ok = find(z>=0);
          x = x(ok); y = y(ok); z = z(ok);
          z = sqrt(z) + c(3,i);
          zindex = sub2ind(size(zbuffer),y+ybias,x+xbias);
          ok = find(z>zbuffer(zindex)); % new "closest points"
          zindex = zindex(ok); x=x(ok); y = y(ok); z = z(ok);
          zbuffer(zindex) = z;
    % compute normalized normals and ambient colors
          normal = [x'-c(1,i) ; y'-c(2,i) ; z'-c(3,i)];
          lengths = sqrt(sum(normal.^2));
          normal = normal ./ [lengths ; lengths ; lengths];
          color = kd(:,i)*ka(i).*ambient;
          color = color * ones(1,length(x));
    % get contributions of each light source
          for j = 1:size(l,2)
              dott = normal' * l(:,j);
              ok = find(dott>=0); dott = dott(ok)';
              color(:,ok) = color(:,ok) + (kd(:,i) .* lc(:,j)) * dott;
              reflect = 2*normal(:,ok).*[dott;dott;dott] ...
                      - l(:,j) * ones(1,length(ok));
              dott = reflect' * v; dott = dott';
              ok2 = dott >= 0; ok = ok(ok2);
              color(:,ok) = color(:,ok) + ks(i) * lc(:,j) * dott(ok2).^n(i);
          end % for j
  % store colors
          [y,x] = ind2sub(size(zbuffer),zindex); rgb = x;
          rgb(:)=1; pixel(sub2ind(size(pixel),y,x,rgb)) = color(1,:);
          rgb(:)=2; pixel(sub2ind(size(pixel),y,x,rgb)) = color(2,:);
          rgb(:)=3; pixel(sub2ind(size(pixel),y,x,rgb)) = color(3,:);
end % for i
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% clip color -- rescale components to fit inside range 0..1
m = max(pixel(:)); if m>1, pixel = pixel/m; end

pixel(pixel<0) = 1; % set background to white

###############################################################################
?
(page 3)

2.

function [hit,where,normal] = raysphere(origin, direction, center, radius)

% [HIT,WHERE,NORMAL] = RAYSPHERE(ORIGIN,DIRECTION,CENTER,RADIUS)
%
% compute the intersection of a ray with a sphere
%
% ORIGIN:        start of ray
% DIRECTION:     normalized direction of ray
% CENTER:        center of sphere
% RADIUS:        radius of sphere
%
% HIT:          whether ray hits sphere ( 0 = no, 1 = yes )
% WHERE:        where the intersection is
% NORMAL:       normalized normal of sphere at intersection
%
% NOTE: assume the origin in strictly outside the sphere --
% no guarantees on result if origin is on or inside

straightline = center - origin;
nearestpoint = straightline * direction';
hit = 0; where = 0; normal = 0;
if nearestpoint <= 0, return, end % pt is at or behind ray origin

nearestpoint = nearestpoint * direction + origin;
distance = norm(nearestpoint - center);
if distance > radius, return, end % no intersection

adjust = sqrt(radius^2 - distance^2);
where = nearestpoint - adjust * direction;
normal = where - center;
normal = normal / norm(normal);

hit = 1;