function find_tiles_video ()
% options : file input '/home/wrightc/find/options2'
%       see READMEoptions2 for input format
% outputs : matlab file format data. 
%   xave, yave, N, x, y, theta, rotim, xmax, ymax, hairpin, connections
%
% this function is similar to find_tiles, only it is meant to process
% multiple subsequent images. (like those that could be composed into a
% video). It accepts AFM files rather than image files. 
%
% processes images according to options
% on the first image, a user clicks to locate tiles on the image. This
% information is used in subsequent processing
% process the image to find tile centers
% allow a user to modify these tile centers
% find connections between tiles and allow user to modify these tile
% conenctions
%
% function calls to other programs : [showafm], [flatten], distance, scale,
% process, modify, modify_conn
% [] = other directories
%
% for help on how to use this code go to
% /home/wrightc/find/READMEfind_tiles_video and READMEoptions2

addpath /home/wrightc/SURF05/program/AFMutils

% read data from file
optfile = '/home/wrightc/SURF05/find/options2';
textform = '%q %u %u %q %u %u %u %u %u %u %u';
[inputpath,first,last,outputpath, a,b,c,d,e,f,g] = textread( optfile, textform );
global OPTIONS;
OPTIONS = [a,b,c,d,e,f,g];

figure(1)
% loop through frames of movie from first to last 
for n = first:last
    str = form(n);
    outpath = [outputpath{1}, '.', str];
    
    
    if ( ~isempty(find(OPTIONS == 10)) )
        % load matloab file
        load( inputpath{1}, '-mat');
    else
        % load image from AFM
        image = showafm([inputpath{1},'.', str], 0) ;
        image = flatten(image, 1);
    end
    image = scale(image);
    
    image = preprocess( image, (n == first) ); 
    image = scale(image);
    
    if ( ~isempty(find(OPTIONS == 11)) )
        theta = 0;
        foundtiles = [];
    else
        % information gathering
        % only performed for first image.
        if (n == first)
            [xave, yave, theta,foundtiles] = firstimage(image);
        end
        save( outpath, 'xave', 'yave', '-mat');
    end
    
    process(image, foundtiles, theta, outpath);
        
    
    % find and modify and store connections
    distance( outpath);
    
    
    % allow user to modify centers and connections again
    again = input( 'modify centers and connections again? [y/n]','s');
    while again == 'y'
        modify( outpath );
        modify_conn( outpath );
        if ( ~isempty(find(OPTIONS == 9)) )
            hairpin(outpath);
        end
        again = input( 'modify centers and connections again? [y/n]','s');
    end
    
    
end

close all;
clear;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function str = form ( in )
% str = form ( in )
% in : numerical input such as 52 or 4
% str character string output of length three such as '052' or '004'
%
% convert 'in' to a 3 digit string
% ex:   5 --> '005' ,   13 --> '013' , 324 --> '324'

if (in < 10)
    str = ['00',int2str(in)];
elseif (in < 100)
    str = ['0',int2str(in)];
else
    str = int2str(in);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = preprocess (image, first)
% out = preprocess (image, first)
% image : grayscale image to be processed
% first : boolean if this image is the first in the video
% out : image after processing
% 
% performs various processing methods according to OPTIONS (global)

global OPTIONS;

% option 6 : contrast adjustment
if ( ~isempty(find(OPTIONS == 6)) )
    image = imadjust(image);
end

% option 1 : gaussian filter
if ( ~isempty(find(OPTIONS == 1)) )
    h = fspecial('gaussian', 3);
    image = filter2(h, image); 
end

% option 2 : median filter
if ( ~isempty(find(OPTIONS == 2)) )
    image = medfilt2(image, [3,3]);
end

% option 3 : masking image
if ( ~isempty(find(OPTIONS == 3)) )
    global SELECT;
    if (first == true)
        disp( 'mark area to mask' )

        SELECT = roipoly( image );
    end
    image = SELECT .* double(image);
    
    [rr,cc] = size (image);
    
    for r = rr:-1:1
        if isempty( find(SELECT(r,:)) )
            image(r, :) = [];
        end
    end
    for c = cc:-1:1
        if isempty( find(SELECT(:,c)) )
            image(:, c) = [];
        end
    end    
end
% option 4 : masking image for each frame
if ( ~isempty(find(OPTIONS == 4)) )
    disp( 'mark area to mask' )
    select = roipoly( image );
    image = select .* double(image);
    
    [rr,cc] = size (image);
    
    for r = rr:-1:1
        if isempty( find(select(r,:)) )
            image(r, :) = [];
        end
    end
    for c = cc:-1:1
        if isempty( find(select(:,c)) )
            image(:, c) = [];
        end
    end  
end

out = image;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [xave, yave, theta,foundtiles] = firstimage ( image )
% [xave, yave, theta, mask,foundtiles] = firstimage ( image )
% image : grayscale input image to use
% xave, yave : average distance between adjacent tiles
% theta : rotation angle of crystal (ccw, degrees)
% mask : image of mask to use for convolution
% foundtiles : tiles located by user

% counts number of values currently averaged into each variable
xnum = 0;   ynum = 0;   tnum = 0;   mnum = 0;
% current average value
xave = 0;   yave = 0;   theta = 0;   mask = [];
% array to store tiles that user has already located for us 
foundtiles = [];

[ymax, xmax] = size(image);

figure(1)
    imagesc(image)
    colormap gray
    title 'image file -- locate 5 adjacent tile'
    
    
 disp( 'locate 5 adjacent tiles on image')
    disp ( ' 1   2 ' )
    disp ( '   3   ' )
    disp ( ' 4   5 ' ) 
again ='y';
while again =='y'
    % user enters 5 tiles 
    %  1   2
    %    3
    %  4   5
    % all tiles conncted to tile 3
    figure(1)
        hold on
        x = zeros(1,5);
        y = zeros(1,5);
        for i = 1:5
            [x(i),y(i)] = ginput(1);
            x(i) = round(x(i));
            y(i) = round(y(i));
            plot( x(i),y(i), 'r+')
        end
        hold off
                
    good = input( 'accept chosen tiles? [y/n] ', 's');
    if (good == 'n')
        hold on
            plot( x,y, 'k+')
            hold off
    elseif (good == 'y')
        % store locations of tiles
        foundtiles = [foundtiles; [x(:),y(:)]];
        
        % find first x distance -------------------------
        temp = ( x(1)-x(2) )^2 + ( y(1)-y(2) )^2;
        temp = sqrt(temp);
        xave = ( xave * xnum + temp ) / (xnum + 1);
        xnum = xnum + 1;
        % find second x distance
        temp = ( x(4)-x(5) )^2 + ( y(4)-y(5) )^2;
        temp = sqrt(temp);
        xave = ( xave * xnum + temp ) / (xnum + 1);
        xnum = xnum + 1;

        % find first y distance -------------------------
        temp = ( x(1)-x(4) )^2 + ( y(1)-y(4) )^2;
        % need to divide by two because this distance is twice the actual dist.
        temp = sqrt(temp) / 2;
        yave = ( yave * ynum + temp ) / (ynum + 1);
        ynum = ynum + 1;
        % find second y distance
        temp = ( x(2)-x(5) )^2 + ( y(2)-y(5) )^2;
        temp = sqrt(temp) / 2;
        yave = ( yave * ynum + temp ) / (ynum + 1);
        ynum = ynum + 1; 

        %find first theta value -------------------------
        temp = ( y(2)-y(1) ) / ( x(2)-x(1) );
        temp = atand(temp);
        theta = (theta * tnum + temp) / (tnum +1);
        tnum = tnum + 1;
        % 2nd theta
        temp = ( x(2)-x(5) ) / ( y(5)-y(2) );
        temp = atand(temp);
        theta = (theta * tnum + temp) / (tnum +1);
        tnum = tnum + 1;    
        % 3rd theta
        temp = ( y(5)-y(4) ) / ( x(5)-x(4) );
        temp = atand(temp);
        theta = (theta * tnum + temp) / (tnum +1);
        tnum = tnum + 1;
        % 4th theta
        temp = ( x(1)-x(4) ) / ( y(4)-y(1) );
        temp = atand(temp);
        theta = (theta * tnum + temp) / (tnum +1);
        tnum = tnum + 1; 
        
        % add known tiles to mask------------------------
        % create new mask if none exists
        if( isempty(mask) )
            msizex = oddround(xave);
            msizey = oddround(yave);
            yrange = floor(msizey/2);
            xrange = floor(msizex/2);
            mask = zeros( msizey, msizex );
        end
        
        % rotate only a portion of the image.
        a = 3;
        lowy = y(3)-a*msizey;
        highy = y(3)+a*msizey;
        lowx = x(3)-a*msizex;
        highx = x(3)+a*msizex;
        if lowy < 1
            lowy = 1;
        end
        if lowx < 1
            lowx = 1;
        end
        if highy > ymax
             highy = ymax;
        end
        if highx > xmax
            highx = xmax;
        end
        timage = image(lowy:highy,lowx:highx);
        timage = imrotate(timage, theta, 'bicubic');
        
        temp = zeros(size(image));
        for i=1:5
            temp(y(i),x(i)) = 1;
        end
        temp = temp(lowy:highy,lowx:highx);
        temp = bwmorph( temp, 'dilate');
        temp = imrotate(temp, theta);
        temp = bwmorph(temp > 0, 'shrink', 2);
        [ty,tx] = find(temp);
        
        temp = zeros(msizey,msizex);
        for i = 1:5
            temp = temp + double(timage( (ty(i)-yrange):(ty(i)+yrange), (tx(i)-xrange):(tx(i)+xrange) ));
        end
        mask = (mask * mnum + temp) / (mnum+5);
        mnum = mnum + 5;
                 
        figure(2)
            imagesc(mask)
            colormap gray
            title 'current tile signature'
    end
    
    again = input('go again?[y/n] ', 's');
end

if xnum == 0
    [xave, yave, theta,foundtiles] = firstimage ( image );
else
    % adjust mask & display
    % average about expected axes of symmetry
    mask = mask + fliplr(mask) + flipud(mask) + rot90(mask,2);
    mask = mask / 4;
    mask = scale(mask .* mask);
        figure(2)
        imagesc(mask)
        colormap gray
        title 'current tile signature'

        global MASK;
        MASK = mask;
end
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = oddround(in)
% out = oddround(in)
% in : number to be reounded
% out : odd number nearest to in
%
% rounds to the nearest odd number

out = ceil(in)/2;
if( floor(out) == ceil(out) )
    out = floor(in);
else
    out = ceil(in);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = evenround(in)
% out = oddround(in)
% in : number to be reounded
% out : even number nearest to in
%
% rounds to the nearest even number

out = ceil(in)/2;
if( floor(out) == ceil(out) )
    out = ceil(in);
else
    out = floor(in);
end