%  MATLAB code to display a black/white dot stereogram
%  given a 2-d depth matrix of an image.

%  Assuming the display will have 72 pixels per inch
DPI = 72;       
%  Eye separation in terms of pixels (avg 2.5 inches)
E = 2.5*DPI;    
%  Depth of field (may determine difficulty for convergence of eyes)
u = 1/3;
%  Recast the depthmap to double so arithmetic procedures may be used on
%  it.  Scale the depth map of the image to be between 0 and 1
depthmap = double(depthmap);
depthmap = depthmap./(max(max(depthmap)));
%  Stereo separation corresponding to the far plane
far = 90;
%  Variables to store the stereogram
image = depthmap*0;
repmask = depthmap*0;
%  Size of depth map of the image to be rendered
[sizey sizex] = size(depthmap);
%  Processing line by line to determine what points will be colored
%  the same 
%  Separation Distance depending on where point lies
for j = 1:sizey
    for i = 1:sizex
        if(i < (sizex/2-DPI))
            D = 12*DPI;
            imsize = 5*DPI;
            inchperpix = imsize/640*i;
            one = (imsize/2 - 1.25*DPI) - inchperpix;
            two = (1-u*depthmap(j,i))*D;
            three = (1-u*depthmap(j,i))/(2-u*depthmap(j,i));
            four = sqrt(three^2 + two^2);
            space(j,i) = 2*four/(sqrt(((2-u*depthmap(j,i))*D)^2 + one^2));
        else
            D = 12*DPI;
            imsize = 5*DPI;
            inchperpix = imsize/640*i;
            one = (imsize/2 + 1.25*DPI) - inchperpix;
            two = (1-u*depthmap(j,i))*D;
            three = (1-u*depthmap(j,i))/(2-u*depthmap(j,i));
            four = sqrt(three^2 + two^2);
            space(j,i) = 2*four/(sqrt(((2-u*depthmap(j,i))*D)^2 + one^2));
        end
    end
end

for y = 1:sizey
    %  Variable arrays to store the color of the pixel and right
    %  reference point that must contain the same color, respectively
    colorpix = zeros(1, sizex);
    rightpix = 1:sizex;
    %  For each pixel within the current surveyed line determine its
    %  color and matching right pixel point
    for x = 1:sizex
        %  Calculate the separation between the two rays
        s = round(space(y,x)*90);
        %  Distance between the pixel that hits the image plane
        left = x - round(s/2);
        right = left + s;
        %  Perform hidden surface removal that occurs when the line
        %  of sight from one point is obscured by another point in
        %  the image
        if (1 <= left) & (right <= sizex)
            t = 1;
            visible = 1;
            %  Z coordinate of the current (x,y) point being evaluated
            zcoord = 0;
            %  While the point is still visible remain checking if
            %  any other point is blocking it
            while (visible & (zcoord < 1))
                zcoord = depthmap(y,x) + 2*(2-u*depthmap(y,x))*t/(u*E);
                visible = (depthmap(y,x-t));
                t = t+1;
            end;
            %  If the (x,y) point is not hidden, record
            if (visible)
                l = rightpix(left);
                %  Ensure that either the left point matches itself
                %  or one right point is matched up
                while ((l ~= left) & (l ~= right))
                    if (l)
                        left = l;
                        l = rightpix(left);
                    else
                        rightpix(left) = right;
                        left = right;
                        l = rightpix(left);
                        right = l;
                    end;
                end;
                rightpix(left) = right;
            end; 
        end;
    end;
    
    %  Assign colors to each point on the image on the basis if 
    %  its a random color or restricted to be the same as its 
    %  right or left point
    for x = fliplr(1:sizex)
        if (rightpix(x) == x) 
            colorpix(x) = (rand(1)>0.5);
        else 
            colorpix(x) = colorpix(rightpix(x));
        end;
        image(y,x) = colorpix(x);
    end;
    repmask(y,:) = rightpix;
end;

%  Convert image back to unsigned int to display
image = uint8(image);
imagesc(image);
colormap(gray);