function [S,Selected, Interm, Interm2, I,S1,Sel1] = mitchell_raytrace(spheres, planes, lights, cop, dov, xres, yres)
%This program attempts to implement Mitchell's Siggraph '87 paper on antialiased
%images at low (adaptive) sampling densities.  This is as opposed to raytrace_with_
%planes2, which will for example shoot 64 jittered samples per pixel, which is probably
%more than we could get away with for a similar quality image.
%The general method is:
%1. Use some good nonuniform sampling method that averages about 1 sample per pixel.
%2. Compute contrasts over 8 or 9 sample point areas.  If the contrast is high, increase
%   increase the sample density in that area.  For simplicity, high contrast will be 
%   supersampled at the highest possible density.
%3. Reconstruction: We'll see about that.  Mitchell's algortithm seems depressing, just
%   averaging the samples over larger areas.
S = zeros(4*xres,4*yres,3);
%S will be a sparse array of samples, that will be averaged for the final image. The -1's
%are so that we don't include unsampled points in the reconstruction.
Selected = zeros(4*xres,4*yres);    %an image of zeros and ones to tell sample concentrations.
%The D is a diffusion grid that will be used to select sample points.
D = zeros(4*xres+2,4*yres+1);
D(1,:) = rand(1,4*yres+1);%-.5*ones(1,4*yres+1) + 
D(:,1) = rand(4*xres+2,1);%-.5*ones(4*xres+2,1) + 
D(4*xres+2,:) = rand(1,4*yres+1);%-.5*ones(1,4*yres+1) + 
%Get the projection matrix.
[Pwfi, foc_v] = Pwfi(cop,dov,0.05,xres,yres);
%Now get the low density samples.
for j = 2:4*yres+1
    for i = 2:4*xres+1
        T = (4*D(i-1,j)+D(i-1,j-1)+2*D(i,j-1)+D(i+1,j-1))/8 + (1/16 + (-1/64 + rand(1)/32));
        sel = 0;
        if T >= .5
            sel = 1;
            S(i-1,j-1,:) = single_raytrace(spheres,planes,lights,cop,Pwfi,(i-1)/4,(j-1)/4);
            Selected(i-1,j-1) = 1;
        end
        D(i,j) = T - sel;
    end
end  
S1 = S;
Sel1 = Selected;
%Now we have to supersample areas based on contrast of the low density samples within that area.
%For each area, we have to find the min greater than zero in S and max, and determine by thresholding
%(Imax-Imin)/(Imax+Imin), .4 for red, .3 for green, .6 for blue.
for i = 1:12:4*xres-12  %for every area
    for j = 1:12:4*yres-12
        %Test this 3x3ish area for contrast and supersample if high.  We have to find the min and max
        %of each of the three channels, r g b.
        min = [1.1 1.1 1.1]; max = [-1 -1 -1];
        %for every color, for the area.
            for s = 0:1:11
                for t = 0:1:11
                    x = i+s; y = j+t;
                    if Selected(x,y) > 0
                        for c = 1:3
                            if S(x,y,c) < min(c)
                                min(c) = S(x,y,c);
                            end
                            if S(x,y,c) > max(c)
                                max(c) = S(x,y,c);
                            end
                        end
                    end
                end
            end
        %Now we have the min and max intensities for the three colors in this area.
        C = zeros(1,3);
        if min(1) < .0001 & max(1) < .0001
            C(1) = 0;
        else
            C(1) = (max(1)-min(1))/(max(1)+min(1));
        end
        if min(2) < .0001 & max(2) < .0001
            C(2) = 0;
        else
            C(2) = (max(2)-min(2))/(max(2)+min(2));
        end
        if min(3) < .0001 & max(3) < .0001
            C(3) = 0;
        else
            C(1) = (max(3)-min(3))/(max(3)+min(3));
        end
        %C = (max-min)./(max+min);
%         if (C(1) > .4 | C(2) > .3 | C(3) > .6) & (C(1) < 1 & C(2) < 1 & C(3) < 1)
        if (C(1) > .35 & min(1) < 1.1 & max(1) > 0) | (C(2) > .25 & min(2) < 1.1 & max(2) > 0) | (C(3) > .5 & min(3) < 1.1 & max(3) > 0)
            %min
            %max
            %These contrast thresholds were determined 
            %experimentally by Mitchell.
            %The contrast is high so supersample the entire 3x3 area. Pick 72 more out of the 
            %144 possibilities in the area (some are already low density samples of course).
            adsamp = 0; err = 0; %number of additionally chosen samples. 
            while adsamp < 72 & err < 2000
                pi = round(i + 11*rand(1)); pj = round(j + 11*rand(1));   %potential ray casting direction.
                if Selected(pi,pj) == 0
                    Selected(pi,pj) = 1;
                    S(pi,pj,:) = single_raytrace(spheres,planes,lights,cop,Pwfi,pi/4,pj/4);
                    adsamp = adsamp + 1;
                else
                    err = err + 1;
                end
            end
            if err == 2000
                ['error1'];
            end
        end
    end
end
%Now we have a bunch of non-uniform samples that well represent the scene.  We must now reconstruct.
%Mitchell uses successively larger averaging boxes, which sounds really bad, but hey, we'll give it
%a shot.  Remember since the grid is 16 grid points per pixel, the final image will be xres*yres.
Interm = S;
SI = Selected;
wei = [.4 .22 .22 .16];
for i = 1:4*xres-1
    for j = 1:4*yres-1
        sum_weights = 0;
        if Selected(i,j) == 1
            Interm(i,j,:) = wei(1)*Interm(i,j,:);
            sum_weights = sum_weights + wei(1);
        end
        if Selected(i+1,j) == 1
            Interm(i,j,:) = Interm(i,j,:) + wei(2)*Interm(i+1,j,:);
            sum_weights = sum_weights + wei(2);
        end
        if Selected(i,j+1) == 1
            Interm(i,j,:) = Interm(i,j,:) + wei(3)*Interm(i,j+1,:);
            sum_weights = sum_weights + wei(3);
        end
        if Selected(i+1,j+1) == 1
            Interm(i,j,:) = Interm(i,j,:) + wei(4)*Interm(i+1,j+1,:);
            sum_weights = sum_weights + wei(4);
        end
        if sum_weights > 0
            Interm(i,j,:) = Interm(i,j,:)/sum_weights;
            SI(i,j) = 1;
        end
    end
end
Interm2 = Interm;
for i = 1:4*xres-1
    for j = 1:4*yres-1
        sum_weights = 0;
        if SI(i,j) == 1
            Interm2(i,j,:) = wei(1)*Interm2(i,j,:);
            sum_weights = sum_weights + wei(1);
        end
        if SI(i+1,j) == 1
            Interm2(i,j,:) = Interm2(i,j,:) + wei(2)*Interm2(i+1,j,:);
            sum_weights = sum_weights + wei(2);
        end
        if SI(i,j+1) == 1
            Interm2(i,j,:) = Interm2(i,j,:) + wei(3)*Interm2(i,j+1,:);
            sum_weights = sum_weights + wei(3);
        end
        if SI(i+1,j+1) == 1
            Interm2(i,j,:) = Interm2(i,j,:) + wei(4)*Interm2(i+1,j+1,:);
            sum_weights = sum_weights + wei(4);
        end
        if sum_weights > 0
            Interm2(i,j,:) = Interm2(i,j,:)/sum_weights;
            SI(i,j) = 1;
        end
    end
end
I = zeros(xres,yres,3); %This is our final image.  Now To fill it up.
%This last filter is like: [.9 1 1 .9;1 4 4 1;1 4 4 1; .9 1 1 .9]
for i = 1:4:4*xres-3
    for j = 1:4:4*yres-3
        sum_weights = 0;
        x = (i+3)/4; y = (j+3)/4;
        for s = 0:1:3
            if SI(i+s,j) == 1
                if s == 0 | s == 3
                    q = .9;
                else
                    q = 1;
                end
                I(x,y,:) = I(x,y,:) + q*Interm2(i+s,j,:);
                sum_weights = sum_weights + q;
            end
            if SI(i,j+s) == 1 & s > 0
                if s == 3
                    q = .9;
                else
                    q = 1;
                end
                I(x,y,:) = I(x,y,:) + q*Interm2(i,j+s,:);
                sum_weights = sum_weights + q;
            end
            if SI(i+(3-s),j+3) == 1 & s < 3
                if s == 0
                    q = .9;
                else
                    q = 1;
                end
                I(x,y,:) = I(x,y,:) + q*Interm2(i+3-s,j+3,:);
                sum_weights = sum_weights + q;
            end
            if SI(i+3,j+3-s) == 1 & s < 3 & s > 0
                I(x,y,:) = I(x,y,:) + Interm2(i+3,j+3-s,:);
                sum_weights = sum_weights + 1;
            end
        end
        if SI(i+1,j+1) == 1
            I(x,y,:) = I(x,y,:) + 1*Interm2(i+1,j+1,:);
            sum_weights = sum_weights + 1;
        end
        if SI(i+1,j+2) == 1
            I(x,y,:) = I(x,y,:) + 1*Interm2(i+1,j+2,:);
            sum_weights = sum_weights + 1;
        end
        if SI(i+2,j+1) == 1
            I(x,y,:) = I(x,y,:) + 1*Interm2(i+2,j+1,:);
            sum_weights = sum_weights + 1;
        end
        if SI(i+2,j+2) == 1
            I(x,y,:) = I(x,y,:) + 1*Interm2(i+2,j+2,:);
            sum_weights = sum_weights + 1;
        end
        if sum_weights > 0
            I(x,y,:) = I(x,y,:)/sum_weights;
        end
    end
end