%  this Matlab code generates an autosterogram of a Torus
%  instead of using random pattern.
% 'Z' is the depth map of the image, normalized from 0 to 1
% 'pattern' is the bitmap [y,x,3] : the background pattern you want to "cover up"
%  the image depth map
%  only restriction is that the pattern bitmap must be at least 90 pixels
%  wide

% first read in the depth map
Z=imread('torusDepth.jpg');
% next we must first change the type to double, and then normalize 
% to have range from 0 to 1
Z=double(Z);
Z=Z./max(max(Z));
% now we read in the pattern
pattern=imread('RANDOM.jpg');
pattern=double(pattern);

%%%%%%%%%%%% now we are ready to actually start implmenting the stereogram

% first define a couple of variables we will need

DPI = 72;   % this is the number of pixels per inch we will assume our monitor to have
EyePix = 2.5*DPI;  % to find the # of pixels spanned by the eye separation, we multiply DPI with 2.5 in, the assumed eye
              % separation
fractionvd = 1/3;  % this variable represents the fraction of viewing distance = Depth of field

posZ = 90;				    % stereo separation corresponding to position Z, and corresponding to far plane Z=0
                         % separation(0) = 90
[ypixMax xpixMax] = size(Z);  % this represents the ypixMax by xpixMax pixels of the image and object

% check size of bitmap pattern to run properly
for j = 1:ypixMax
    for i = 1:xpixMax
        if(i < (xpixMax/2-DPI))
        E = 2.5*DPI;
        D = 12*DPI;
        imsize = 5*DPI;
        inchperpix = imsize/640*i;
        one = (imsize/2 - 1.25*DPI) - inchperpix;
        two = (1-fractionvd*Z(j,i))*D;
        three = (1-fractionvd*Z(j,i))/(2-fractionvd*Z(j,i)) * one;
        four = sqrt(two^2 + three^2);
        r2(j,i) = 2*four/(sqrt(((2-fractionvd*Z(j,i))*D)^2 + one^2));
        else
        E = 2.5*DPI;
        D = 12*DPI;
        imsize = 5*DPI;
        inchperpix = imsize/640*i;
        one = (imsize/2 + 1.25*DPI) - inchperpix;
        two = (1-fractionvd*Z(j,i))*D;
        three = (1-fractionvd*Z(j,i))/(2-fractionvd*Z(j,i)) * one;
        four = sqrt(three^2+two^2);
        r2(j,i) = 2*four/(sqrt(((2-fractionvd*Z(j,i))*D)^2 + one^2));
        end
    end
end



[patternY,patternX,dep] = size(pattern);
if patternX < posZ
    disp('pattern width has to be greater than separation');
    return;
end


imageout = zeros(ypixMax,xpixMax,dep);  % this will be the output variable of our stereogram generation
repmask = zeros(ypixMax,xpixMax); % this is just to be able to understand the contraints.


% the outer "for" loop takes each horizontal line one by one 
for y = 1:ypixMax

    colorpix = zeros(1,xpixMax,dep);    % color of this pixel

    % pattern copied into right most portion of colorpix array

    reyex = xpixMax;  % X-values corresponding to right eye
    leyex = xpixMax-patternX+1; % X-values corresponding to left eye
    vertline = mod((y-1),patternY)+1; % variable indicating the vertical line of the pattern in this loop

    colorpix(1,leyex:reyex,:) = pattern(vertline,:,:); % copy pattern into right side of colorpix array
       
    % copy a reversed line of the pattern into rest of image
    % a reversed line is copied because this avoids black patches : to avoid
    % accidental alignment

    patternrev = pattern(vertline,linspace(patternX,1,patternX),:);
    finished = 0;
    reyex = leyex-1;
    leyex  = reyex-patternX+1;
    while ~finished
  
        colorpix(1,leyex:reyex,:) = patternrev(1,patternX-(reyex-leyex):patternX,:);
        reyex = leyex-1;
        leyex = reyex-patternX+1;
        if reyex<1
            finished = 1;
        else
            if leyex<1
                leyex = 1;
            end
        end
    end



    record = 1:xpixMax;      % initially, each pixel is linked to itself

    
    for x = 1:xpixMax
        sep = round(r2(y,x)*90);
        left = x - round(sep/2);
        right = left + sep;
        if (1 <= left) & (right <= xpixMax)
            % we first implement the hidden-surface removal
            % we will check the points (x-index,y) and (x+index,y)
            index = 1;      
            vis = 1; 
            zindex = 0;   % Z-coord of ray at these 2 points, initialized at 0
            while (vis & ( zindex < 1 ))   % condition to continue hidden-surface removal
                zindex = Z(y,x) + 2*(2 - fractionvd*Z(y,x))*index/(fractionvd*EyePix);
                vis = (Z(y,x-index));
                index = index+1;
            end

            if (vis)        %record the fact that pixels at
                lrecord = record(left); % left and right are record
                while ((lrecord ~= left) & (lrecord~= right))
                    % move pointers until record[left]=left or
                    % record[left]=right
                    if (lrecord)
                        left = lrecord;
                        lrecord = record(left);
                    else
                        record(left) = right;
                        left = right;
                        lrecord = record(left);
                        right = lrecord;
                    end
                end % while
                record(left) = right;   % this records it
            end % if vis
        end % if (lrecord<= left
    end % for x

    %disp(record);

    for x = fliplr(1:xpixMax)  % set pixels on this scan line
        if (record(x) == x)
            % colorpix(1,x,:) = (rand(1)>0.5);  -- comment out this line, only
            % for Random Dot Stereograms
        else
            colorpix(1,x,:) = colorpix(1,record(x),:);   % constrained choice
        end
        imageout(y,x,:) = colorpix(1,x,:);
    end
    disp(['line' num2str(y)]);   % this shows in the Matlab output the line scanned
    repmask(y,:) = record;
end % for y

imageout=uint8(imageout);
imagesc(imageout)
