function [I2,I] = raytrace_with_planes2(spheres, planes, lights, cop, dov, xres, yres,z,sig)
%Does planes in addition to the things that raytrace does. (Not well documented).
%
%Returns an image of the spheres with the light sources
%from the given center of projection and direction of view.
%The spheres list is n by 14 (center, radius, material
%properties).
%the true color image.
I = zeros(z*xres,z*yres,3);
jI = zeros(z*xres,z*yres,2);
%First, set up the Pwfi matrix.
focal_length = .050;
[Pwfi, foc_v] = Pwfi(cop,dov,focal_length,xres,yres);
%Now the flow of the code will be as such:
%for each ray
%   for each sphere
%       detect intersections
%       if smaller t less than current-small t
%           update closest sphere and current-small t
%   if there was a sphere
%       for each light
%           compute specular and diffuse terms to calculate 
%           lighting at that intersection point to light pixel.
%
%And hopefully that should do it.
q = size(spheres);
n_spheres = q(1);
q = size(lights);
n_lights = q(1);
q = size(planes);
n_planes = q(1);
%A ray is given by cop + t*Pwfi.[u v 1]' where u and v are steps
%in xres and yres horizontal and vertical pixels.
e = cop; %just the eye.  Makes for less typing below.
for i = 1:xres
    for j = 1:yres
    for sami = 1:z
    for samj = 1:z
        
            
        min_t = 1e+90; 
        closest_s = 0; closest_p = 0;  s_or_p = 0;
        %remember t=1 is exactly the image plane and so should 
        %be the minimum allowable intersection. s_or_p keeps track
        %of what type of object was closest, -1 for sphere, 1 for plane.
        
        %Just for fun, jitter the sample a little.
        ii = i + (-.5 + rand);
        jj = j + (-.5 + rand);
        w = z*xres - (z*(j-1)+samj-1);           %this is for the wierd matlab output.
        w2 = z*(i-1) + sami;
       %w = z*(i-1) + sami; w2 = z*(j-1) + samj;
        jI(z*(i-1) + sami,z*(j-1) + samj,:) = [ii jj];
        v = Pwfi*[ii jj 1]';
        for s = 1:n_spheres
            c = spheres(s,1:3);  %current center.
            r = spheres(s,4);    %radius
            %now cop + t*v is the ray and we need to solve it 
            %simultaneously with the sphere.  The following turns
            %out to be the solution to the quadratic equation.
            t = intersect(e,v,c,r);  %intersects the sphere with the ray, returning t.
            if abs(imag(t(1))) < 1e-20
                if (real(t(1)) < min_t) & (real(t(1)) > 1)
                    min_t = real(t(1));
                    closest_s = s;
                    s_or_p = -1;
                end
            end
            if abs(imag(t(2))) < 1e-20
                if (real(t(2)) < min_t) & (real(t(2)) > 1)
                    min_t = real(t(2));
                    closest_s = s;
                    s_or_p = -1;
                end
            end
        end
        %Now check for closer planes.
        for p = 1:n_planes
            p_on_plane = planes(p,1:3);
            p_normal = planes(p,4:6);
            t = p_intersect(e,v',p_on_plane,p_normal);
            if t > 1 & t < min_t
                min_t = t;
                closest_p = p;
                s_or_p = 1;
            end
        end
        %Now, if there was a closest sphere or plane, let's do some light computation.
        if s_or_p == -1 | s_or_p == 1
        %if closest_s > 0
            %Pick up the material properties.
            sum_l = [0 0 0]; %for the light properties.
            if s_or_p == -1
                am = spheres(closest_s, 5:7);
                dif = spheres(closest_s, 8:10);
                spe = spheres(closest_s, 11:13);
                ex = spheres(closest_s, 14);
            else 
                if s_or_p == 1
                    am = planes(closest_p,7:9);
                    dif = planes(closest_p,10:12);
                    spe = planes(closest_p,13:15);
                    ex = planes(closest_p,16);
                else
                    ['error1'];
                end
            end
            %And the normal of the sphere at the point of intersection.
            %Remember am, dif, and spe are 1x3 vectors, for r,g,b.
            p = cop + min_t*v';  %the point that the ray intersects the sphere or plane.
                                 %I'm reusing variables for different purposes: bad.
            if s_or_p == -1
                n = p - spheres(closest_s,1:3);
            else 
                if s_or_p == 1
                    n = planes(closest_p,4:6);
                end
            end    
            n = n/norm(n);      %the unit normal on the sphere.
            view = cop - p;
            view = view/norm(view); %the unit view direction vector.

            %I(w,i,1:3) = am;
            for k = 1:n_lights
                temp2 = zeros(1,1,3);
                l = lights(k,1:3);      %current light.
                l(2) = l(2) + .1*(-.5 + rand);
                l(3) = l(3) + .1*(-.5 + rand);
                sum_l = sum_l + lights(k,4:6);   %the light's color.
                l_dif = lights(k,7:9);  %no idea what a diffuse filter means.
                l_spe = lights(k,10:12);%seems like a softening of the light. 
                l = l - p;
                l = l/norm(l);  %unit light vector from the point.
                %Now we have to do shadow calculation, which is a while loop
                %through all of the other spheres, doing lighting only when 
                %we can.
                s = 1; shad = 0;
                while s <= n_spheres & shad == 0
                    if s ~= closest_s
                        t = intersect(p,l,spheres(s,1:3),spheres(s,4));
                        if (abs(imag(t(1))) < 1e-20 & (real(t(1)) > 0)) | (abs(imag(t(2))) < 1e-20 & (real(t(2)) > 0))
                            shad = 1;
                        end
                    end
                    s = s + 1;
                end
                q = 1;  %go through the planes for shadows.
                if closest_p ~= 1
                    x = [0];
                end
                while q <= n_planes & shad == 0
                    if q ~= closest_p
                        t = p_intersect(p,l,planes(q,1:3),planes(q,4:6));
                        if t > 0 & t < norm(lights(k,1:3) - p)
                            shad = 1;
                        end
                    end
                    q = q + 1;
                end
                if shad == 0
                    do = dot(l,n);
                    if do < 0
                        do = 0;
                    else
                        ref = 2*do*n - l;
                        do2 = dot(view,ref);
                        if do2 < 0
                            do2 = 0;
                        end
                        temp2 = do*dif.*l_dif + (do2^ex)*spe.*l_spe;
                    end
                    %size((I(i,j,1:3)))
                    %I(i,j,1:3) = (I(i,j,1:3)) + do*dif + ((dot(view,ref))^ex)*spe;
                    I(w,w2,1) = (I(w,w2,1)) + temp2(1);
                    I(w,w2,2) = (I(w,w2,2)) + temp2(2);
                    I(w,w2,3) = (I(w,w2,3)) + temp2(3);
                end
            end
            %avg_l = avg_l/n_lights;
            if sum_l(1) > 1
                sum_l(1) = 1;
            end
            if sum_l(2) > 1
                sum_l(2) = 1;
            end
            if sum_l(3) > 1
                sum_l(3) = 1;
            end
            I(w,w2,1) = (I(w,w2,1)) + am(1)*sum_l(1);
            I(w,w2,2) = (I(w,w2,2)) + am(2)*sum_l(2);
            I(w,w2,3) = (I(w,w2,3)) + am(3)*sum_l(3);
            if I(w,w2,1) > 1
                I(w,w2,1) = 1;
            end
            if I(w,w2,2) > 1
                I(w,w2,2) = 1;
            end
            if I(w,w2,3) > 1
                I(w,w2,3) = 1;
            end
        end
    end
    end
    end
end

%At this point I is an image and jI is the corresponding jittered pixel points from which
%the I values were computed.  So now we need to perform a filtering operation.
sum = [0 0 0 0]; g = [0 0 0 0];

I2 = zeros(xres,yres,3);
for i = 1:2
    for j = 1:2
        for r = 1:z
            for s = 1:z
                x1 = z*(i-1)+r; y1 = z*(j-1)+s; x2 = z*(xres-i+1)-r+1; y2 = z*(yres-j+1)-s+1;
                g(1) = g2d(jI(x1,y1,1),jI(x1,y1,2),[1 1],sig);
                g(2) = g2d(jI(x2,y2,1),jI(x2,y2,2),[xres yres],sig);
                g(3) = g2d(jI(x1,y2,1),jI(x1,y2,2),[1 yres],sig);
                g(4) = g2d(jI(x2,y1,1),jI(x2,y1,2),[xres 1],sig);
                sum(1) = sum(1) + g(1);
                sum(2) = sum(2) + g(2);
                sum(3) = sum(3) + g(3);
                sum(4) = sum(4) + g(4);
                I2(1,1,1)= I2(1,1,1) + g(1)*I(x1,y1,1);
                I2(1,1,2)= I2(1,1,2) + g(1)*I(x1,y1,2);
                I2(1,1,3)= I2(1,1,3) + g(1)*I(x1,y1,3);
                I2(xres,yres,1) = I2(xres,yres,1) + g(2)*I(x2,y2,1);
                I2(xres,yres,2) = I2(xres,yres,2) + g(2)*I(x2,y2,2);
                I2(xres,yres,3) = I2(xres,yres,3) + g(2)*I(x2,y2,3);
                I2(1,yres,1) = I2(1,yres,1) + g(3)*I(x1,y2,1);
                I2(1,yres,2) = I2(1,yres,2) + g(3)*I(x1,y2,2);
                I2(1,yres,3) = I2(1,yres,3) + g(3)*I(x1,y2,3);
                I2(xres,1,1) = I2(xres,1,1) + g(4)*I(x2,y1,1);
                I2(xres,1,2) = I2(xres,1,2) + g(4)*I(x2,y1,2);
                I2(xres,1,3) = I2(xres,1,3) + g(4)*I(x2,y1,3);
            end
        end
    end
end
I2(1,1,:) = I2(1,1,:)/sum(1); %I2(1,1,2) = I2(1,1,2)/sum(1); I2(1,1,3) = I2(1,1,3)/sum(1);
I2(xres,yres,:) = I2(xres,yres,:)/sum(2); %I2(xres,yres,2) = I2(xres,yres,2)/sum(2); I2(xres,yres,3) = I2(xres,yres,3)/sum(2); 
I2(1,yres,:) = I2(1,yres,:)/sum(3); %I2(1,yres,2) = I2(1,yres,2)/sum(3); I2(1,yres,3) = I2(1,yres,3)/sum(3); 
I2(xres,1,:) = I2(xres,1,:)/sum(4); %I2(xres,1,2) = I2(xres,1,2)/sum(4); I2(xres,1,3) = I2(xres,1,3)/sum(4); 

%Now I'm assuming a square image.
for i = 2:xres-1
    sum(:) = [0 0 0 0];
    for bi = -1:1:1
        for bj = 0:1
            for r = 1:z
                for s = 1:z
                    x1 = z*(i+bi-1)+r; y1 = z*bj+s; x2 = z*(xres-bj-1)+s; y2 = z*(yres-bj-1)+s;
                    g(1) = g2d(jI(x1,y1,1),jI(x1,y1,2),[i 1],sig);
                    g(2) = g2d(jI(y1,x1,1),jI(y1,x1,2),[1 i],sig);
                    g(3) = g2d(jI(x1,y2,1),jI(x1,y2,2),[i yres],sig);
                    g(4) = g2d(jI(x2,x1,1),jI(x2,x1,2),[xres i],sig);
                    %sum(1) = sum(1) + g(1);
                    %sum(2) = sum(2) + g(2);
                    %sum(3) = sum(3) + g(3);
                    %sum(4) = sum(4) + g(4);
                    sum = sum+g;
                    I2(i,1,1) = I2(i,1,1) + g(1)*I(x1,y1,1);
                    I2(i,1,2) = I2(i,1,2) + g(1)*I(x1,y1,2);
                    I2(i,1,3) = I2(i,1,3) + g(1)*I(x1,y1,3);
                    I2(1,i,1) = I2(1,i,1) + g(2)*I(y1,x1,1);
                    I2(1,i,2) = I2(1,i,2) + g(2)*I(y1,x1,2);
                    I2(1,i,3) = I2(1,i,3) + g(2)*I(y1,x1,3);
                    I2(i,yres,1) = I2(i,yres,1) + g(3)*I(x1,y2,1); 
                    I2(i,yres,2) = I2(i,yres,2) + g(3)*I(x1,y2,2); 
                    I2(i,yres,3) = I2(i,yres,3) + g(3)*I(x1,y2,3); 
                    I2(xres,i,1) = I2(xres,i,1) + g(4)*I(x2,x1,1);
                    I2(xres,i,2) = I2(xres,i,2) + g(4)*I(x2,x1,2);
                    I2(xres,i,3) = I2(xres,i,3) + g(4)*I(x2,x1,3);
                end
            end
        end
    end
    I2(i,1,:) = I2(i,1,:)/sum(1);
    I2(1,i,:) = I2(1,i,:)/sum(2);
    I2(i,yres,:) = I2(i,yres,:)/sum(3);
    I2(xres,i,:) = I2(xres,i,:)/sum(4);
end
%Now the middle.
for i = 2:xres-1
    for j = 2:yres-1
        sum = 0;
        for bi = -1:1:1
            for bj = -1:1:1
                for r = 1:z
                    for s = 1:z
                        x = z*(i+bi-1)+r; y = z*(j+bj-1)+s;
                        g = g2d(jI(x,y,1),jI(x,y,2),[i j],sig);
                        sum = sum + g;
                        I2(i,j,1) = I2(i,j,1) + g*I(x,y,1);
                        I2(i,j,2) = I2(i,j,2) + g*I(x,y,2);
                        I2(i,j,3) = I2(i,j,3) + g*I(x,y,3);
                    end
                end
            end
        end
        I2(i,j,:) = I2(i,j,:)/sum;
    end
end