function find_tiles ()
% goal: be able to locate tiles using only minimal user input. 

% read data from file
optfile = '/home/wrightc/find/options';
[inputfile,outputfile, a,b,c,d,e] = textread( optfile, '%q %q %u %u %u %u %u' );
options = [a,b,c,d,e]

% image = imread( inputfile{1}, 'tiff' );
    % load image from AFM and flatten
    image = showafm(inputfile{1}, 0) ;
    image = flatten(image, 1);
      image = scale(image);
    
% optional processing
if ( ~isempty(find(options == 6)) )
    image = imadjust(image);    
end
   figure(1)
    imagesc(image) 
    colormap gray
    title 'image file'

if ( ~isempty(find(options == 1)) )
    h = fspecial('gaussian', 3);
    image = filter2(h, image); 
end
if ( ~isempty(find(options == 2)) )
    % option: median filter
    image = medfilt2(image, [3,3]);    
end
if ( ~isempty(find(options == 3)) )
    disp( 'mark area to mask' )
    mask = roipoly( image );
    image = mask .* double(image);
end

figure(1)
    imagesc(image)
    colormap gray
    title 'image file'

xnum = 0;   ynum = 0;   tnum = 0;   mnum = 0;
xave = 0;   yave = 0;   theta = 0;   mask = [];
% array to store tiles that user has already located for us (can later be
% used as reference points?)
foundtiles = [];

again ='y';
while again =='y'
    % user enters 5 tiles 
    %  1   2
    %    3
    %  4   5
    disp( 'locate 5 adjacent tiles on image')
    figure(1)
        imagesc(image)
        colormap gray
        title 'image file'
    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,y, 'r+')
    end
    hold off
                
    good = input( 'accept chosen tiles? [y/n] ', 's');
    
    if (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 -------------------------
        %angle needed to rotate ccw 
        temp = ( y(2)-y(1) ) / ( x(2)-x(1) );
        temp = atand(temp);
        theta = (theta * tnum + temp) / (tnum +1);
        tnum = tnum + 1;
        %find second theta value
        temp = ( x(2)-x(5) ) / ( y(5)-y(2) );
        temp = atand(temp);
        theta = (theta * tnum + temp) / (tnum +1);
        tnum = tnum + 1;    
        %find third theta value
        temp = ( y(5)-y(4) ) / ( x(5)-x(4) );
        temp = atand(temp);
        theta = (theta * tnum + temp) / (tnum +1);
        tnum = tnum + 1;
        %find fourth theta value
        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------------------------
        if( isempty(mask) )
            msizex = oddround(xave);
            msizey = oddround(yave);
            yrange = floor(msizey/2);
            xrange = floor(msizex/2);
            mask = zeros( msizey, msizex );
        end
        timage = imrotate(image, theta, 'bicubic');
        % find tx and ty (x and y of rotated matrix)
        temp = zeros(size(image));
        for i=1:5
            temp(y(i),x(i)) = 1;
        end
        temp = bwmorph( temp, 'dilate');
        temp = imrotate(temp, theta);
        temp = bwmorph(temp > 0, 'shrink', 2);

        [ty,tx] = find(temp);
        % average 5 tiles into mask.
        for i = 1:5
            temp = double(timage( (ty(i)-yrange):(ty(i)+yrange), (tx(i)-xrange):(tx(i)+xrange) ));
            mask = (mask * mnum + temp) / (mnum +1);
            mnum = mnum + 1; 
        end
        figure(2)
            imagesc(mask)
            colormap gray
            title 'current tile signature'
    end
    
    again = input('go again?[y/n] ', 's');
end
    
xave
yave
theta
    
% use inputs for a convolution ---------------------

% adjust mask & display
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'
            
% rotate image by determined angle  ----------------------       
rotim = imrotate(image, theta, 'bicubic');
% rotate found tiles along with image
temp = zeros(size(image));
for i=1:length(foundtiles)
    temp(foundtiles(i,2), foundtiles(i,1)) = 1;
end
temp = bwmorph( temp, 'dilate');
temp = imrotate(temp, theta);
temp = bwmorph(temp , 'shrink', 2);
[y,x]  = find(temp);
foundtiles = [x,y];

%convolve mask with image -----------------------------------
convim = conv2( double(rotim), mask, 'same' ); 
            figure(3)
imagesc(mask)


% flatten black spots
th = graythresh(convim);
high = max(max(convim));                              %% HERE
temp = (convim > high *th);
convim = convim .* temp;

        figure(2)
            imagesc(convim)
            colormap gray
            title 'convoluted image with mask' 
            hold on
            plot(foundtiles(:,1), foundtiles(:,2),'r+')
            hold off
            
% top hat --------------------------------------------------
if ( ~isempty(find(options == 4)) )        
    se = strel( 'rectangle', [floor(msizey/2), floor(msizex/2)]);
    topim = imtophat(convim,se);
    
    % remove high peaks
    th = graythresh (topim);
    high = max(max(topim)); 
    temp = (topim < th*high);
    topim = topim .* temp;
            
    % average results (remove tiny spots that are not really centers)
    h = fspecial('average', [yrange, xrange]);
    topim = filter2(h, topim);

    % average again to help combine centers later
    topim = filter2(h, topim);
    
    % mask out things that werent even in orig image
    temp = ones(size(image));
    temp = imrotate(temp, theta);
    temp = bwmorph(temp, 'erode', xrange);
    topim = temp .* topim;
    
            figure(3)
            imagesc( topim )
            title 'tophat image with mask' 
            hold on
            plot(foundtiles(:,1), foundtiles(:,2),'r+')
            hold off
end

% turn into a b/w image. hopefully with only areas round centers chosen---
if ( ~isempty(find(options == 5)) )
    % thresh
    [i,j] = size(topim);
    thres = graythresh(topim( (i/2):(i/2+10*msixey), (j/2):(j/2+10*msixex)) );
    high = max(max(topim));
    bwim = (topim > thres*high);
else
    % local max
    bwim = imregionalmax( topim );
end
% mask out things not in orig image, and edge effects 
temp = bwmorph(temp, 'erode');
bwim = temp .* bwim;
% remove vertical bars that result from top hat averaging
bwim = bwmorph(bwim, 'shrink');
temp = bwmorph(bwim, 'clean');
bwim = (bwim & ~temp);
            figure(4)    
            imagesc(bwim);
            colormap gray
            truesize
            title 'bw image'
                       
% combine centers
bwim = imdilate( bwim, ones(yrange, floor(msizex/3) ));
bwim = bwmorph( bwim, 'thin', 7);           % more?
bwim = bwmorph( bwim, 'shrink', 10);
temp = bwmorph(bwim, 'clean');
bwim = (bwim & ~temp);
        
            figure (5)
            imagesc(bwim);
            colormap gray
            truesize
            title 'bw image'
         
% allow user to modify centeres if option is turned on
if ( ~isempty(find(options == 7)) )
    bwim = modify( bwim, rotim );
end
          
[tempy,tempx] = find(bwim);

            figure(2)
            hold on
            plot( tempx, tempy, 'm+')
            hold off
            
            figure(1)
            imagesc (rotim)
            colormap gray
            hold on
            plot( tempx, tempy, 'm+')
            hold off

% save image into file ---------------------------
save = input('save?[y/n] ', 's');
if save == 'y'
    imwrite( bwim, outputfile{1}, 'tiff' );
%     filename = input('filename: ', 's');
%     imwrite( scale(convim) , ['/home/wrightc/program/images/',filename],'tiff')
end

    % TOPHAT?
    
pause;
close all;
clear

%%-------------------------------------------------------------------------
function out = oddround(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)
% rounds to the nearest even number

out = ceil(in)/2;
if( floor(out) == ceil(out) )
    out = ceil(in);
else
    out = floor(in);
end
%%-------------------------------------------------------------------------
function output_image = modify( input_image , raw_image )

output_image = input_image;

figure(8) 
    imagesc(raw_image)
    hold on
    [i,j] = find( input_image );
    plot( j,i , 'b+');

    colormap gray
    truesize
    title 'LEFT click = create center, RIGHT click = remove center, CENTER to end'
    
    
    [xmax, ymax] = size(input_image);
    button = 1;
while button ~= 2   
    [x,y,button] = ginput(1);
    x = round(x);
    y = round(y);
    
    if( x<0 | y<0 | x> xmax | y > ymax)
        break;
    end
    
    if button == 1
        hold off
        output_image(y,x) = 1;
        hold on
        plot(x,y,'g+')
    elseif button == 3
        hold off
        [y,x] = findclosest( output_image, y,x );
        output_image(y,x) = 0;
        hold on
        plot(x,y,'ro')
    end
end    
    hold off
    
%%-------------------------------------------------------------------------
function [i2,j2] = findclosest( input_image, i , j )

if (input_image(i,j))
    i2=i;
    j2=j;
    return;
end

%create image with only selected point in it
point = zeros( size(input_image) );
point(i,j) = 1;
% find distance form that point for all points
point = bwdist( point ); 
% multiply by orig bw image so that only points already in image are counted
point = point .* input_image;
% find the point with the minimum distance from selected point
indicies = find(point);
[c,I] = min( point(indicies) );
I= indicies(I);
[i2,j2] = ind2sub(size(point), I);