below is a sample solution to hw5, including light buffers and transparency. it assumes that our discussion in class is correct: total internal reflection is impossible for spheres. also, notice that some error-checking was included but is now commented out (see the calls to ERROR). ------------------------------------------------------------------------------- c = [ ... 0 -27 -27 -20 0 -19 5 10 0 20 20 100]; r = [20 15 7 10]; c = c*2; r = r*2; % "higher" resolution ka = [.1 .2 .3 .2]; kd = [ ... .7 0 0 .4 0 .8 0 .4 0 0 1 0]; ks = [.7 .5 .4 .1 ]; n = [3 7 8 10 ]; kt = [0 0 0 .95]; eta = [1 1 1 1.2]; ambient = [1 ; 1 ; 1]; l = [ ... 1 -1 -1 -1 -1 -1 ]; lc = [ ... 1 1 1 1 1 1 1 1 1]; tic; blah=raytrace(c,r,ka,kd,ks,n,kt,eta,ambient,l,lc,5); toc; image(blah); axis equal; axis xy; axis tight ------------------------------------------------------------------------------- function [pixel] = raytrace(c,r,ka,kd,ks,n,kt,eta,ambient,l,lc,limit) % [pixel] = raytrace(c,r,ka,kd,ks,n,kt,eta,ambient,l,lc,limit) % return an image pixel (with integer world and view coordinates) % just large enough to show all the spheres under the specified lights % using ray-tracing to a depth of the specified limit. % (limit=0 means do first-hits, but no recursive calls.) % % 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), % transparency coefficient kt(i), index of refraction eta(i) % (the ambient coefficients for each r,g,b component are ka(i)*kd(1:3,i); % the reflect coefficient is kr=kt; % eta(i) is the index with respect to vacuum.). % % the ambient light has color ambient(1:3) % % the j-th light travels in 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; zinfinity = max(c(3,:)) + max(r) + 100; % z-value above entire scene pixel = zeros(high,wide,3)-1; % use -1 to flag background zbuffer = zeros(high,wide) - inf; % z-buffer [itembuffer,bogus1,bogus2] = makebuffer(c,r); v = [0 ; 0 ; 1]; % view vector (towards eye) % normalize and reverse light vector directions % and create light buffers and associated data structures lightTview = zeros(3,3,size(l,2)); % 3x3 Tview matrices for lights % not 4x4 because no translation needed lightbias = zeros(2,size(l,2)); % like xbias, ybias for each light for j = 1:size(l,2) dir = -l(:,j); % direction towards light % rotate dir to set "x" = 0 T = eye(3); if dir(1) ~= 0 scale = sqrt(dir(1)^2 + dir(3)^2); T = [dir(3) 0 -dir(1) ; 0 scale 0 ; dir(1) 0 dir(3)]/scale; dir = T * dir; end % rotate dir to set "y" = 0 scale = sqrt(dir(2)^2 + dir(3)^2); T = [scale 0 0 ; 0 dir(3) -dir(2) ; 0 dir(2) dir(3)] * T / scale; lightTview(:,:,j) = T; l(:,j) = - l(:,j)/norm(l(:,j)); [buf,lightbias(1,j),lightbias(2,j)] = makebuffer(T * c, r); lightbuffer(1:size(buf,1),1:size(buf,2),j) = buf; end % store constant (during ray-trace) data all in one place: data.c = c; data.ambient = ambient; data.r = r; data.l = l; data.ka = ka; data.lc = lc; data.kd = kd; data.lightbias = lightbias; data.ks = ks; data.lightbuffer = lightbuffer; data.n = n; data.lightTview = lightTview; data.kt = kt; data.v = v; data.eta = eta; for x = 1:size(pixel,2) for y = 1:size(pixel,1) knownhit = itembuffer(y,x); if knownhit % trust that sphere knownhit really is visible at this pixel... origin = [x-xbias ; y-ybias ; zinfinity]; pixel(y,x,:) = raycolor(origin,-v,knownhit,data,limit); end % if knownhit end % for y end % for x % clip color -- rescale components to fit inside range 0..1 m = max(pixel(:)); if m>1, pixel = pixel/m; end pixel(pixel<0) = 1; % pity the printer: draw black background as white %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [rgb] = raycolor(origin,direction,knownhit,data,limit) % return color rgb seen by a ray with given origin and direction % with limit recursive calls remaining. % knownhit>0 means assume the ray first hits sphere knownhit; % knownhit=0 means search for which sphere is hit. % knownhit<0 means search but ignore sphere -knownhit. % data holds constant data -- see raytrace above for fields. c = data.c; r = data.r; l = data.l; lc = data.lc; lightbuffer = data.lightbuffer; % compute list of spheres to test for intersection if knownhit>0 ii = knownhit; else ii = 1:size(c,2); if knownhit<0, ii(-knownhit) = []; end end % find nearest intersect: closest, whichclosest, whichwhere are the current % closest distance, closest sphere, pt of intersection, respectively closest = inf; for i = ii [hit,where,bogus] = raysphere(origin, direction, c(:,i), r(i)); if hit distance = norm(where-origin); if distance=1 & yy>=1 & xx<=size(lightbuffer,2) & yy<=size(lightbuffer,1) visible = lightbuffer(yy,xx,j) == i; end if visible % add in contributions for this light dott = normal' * l(:,j); if dott>0 rgb = rgb + dott * data.kd(:,i) .* lc(:,j); reflect = 2*dott*normal - l(:,j); dott = max(0,reflect' * direction); rgb = rgb + dott^data.n(i) * data.ks(i) * lc(:,j); end end % if visible end % for j dott = normal' * direction; reflect = 2*dott*normal - direction; if limit>0 limit = limit-1; ignore = -i; rgb = rgb + data.ks(i) * raycolor(origin,reflect,ignore,data,limit); if data.kt(i)>0 % do transmitted ray iff transparent % note: total internal reflection is impossible for spheres transmit = refract(origin, normal, -direction, 1, data.eta(i)); % compute point where transmitted ray comes out; % use ideas from ray-sphere intersection algorithm from lecture where = origin + (2*r(i)* transmit'*-normal) * transmit; % if abs(norm(where-c(:,i))-r(i))>1e-6, error('ooops'); end mynormal = (c(:,i)-where)/r(i); % note: point INTO the sphere transmit = refract(where,mynormal,transmit,data.eta(i),1); rgb = rgb + data.kt(i)*raycolor(where,transmit,ignore,data,limit); end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [transmit] = refract(where,normal,incident,eta_i,eta_t) % compute normalized, transmitted ray TRANSMIT from an INCIDENT ray % hitting a surface at a point WHERE with given NORMAL; % assume incoming material has index ETA_I of refraction, % and transmitting material has index ETA_T % NOTE: we assume all vectors are normalized and % verify no total internal refraction % use refraction formula from glassner eta_it = eta_i/eta_t; cos_i = normal' * -incident; cos_t = 1 + eta_it^2 * (cos_i^2-1); % if cos_t<0, error('impossible(in): total internal reflection'); end cos_t = sqrt(cos_t); transmit = eta_it*incident + (eta_it*cos_i - cos_t)*normal; % if abs(norm(transmit)-1)>1e-6, error('(in) not normalized'); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [buffer,xbias,ybias] = makebuffer(c,r) % return a buffer indicating which sphere is visible at each pixel; % (store index for a sphere, 0 for no sphere). c,r are as in raytrace y = c(2,:); y = [y-r y+r]; ybias=ceil(1-min(y)); % offset to make y coords >= 1 x = c(1,:); x = [x-r x+r]; xbias=ceil(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; buffer = zeros(high,wide); zbuffer = zeros(high,wide) - inf; % z-buffer 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 closer; update zbuffer, buffer 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; buffer(zindex) = i; end