function [output]=Localizer(bright_file_name, red_file_name, first_pic, lastpic, interactive_mode, ellipse_spacing)

%Set Filenames
for i=first_pic:lastpic
        I=int2str(i);
    picture(i).ch1_filename=strcat(red_file_name, I, '.tif');
    picture(i).ch2_filename=strcat(bright_file_name, I, '.tif');
    %textname(i, :)=strcat(picname, I);
    if ~exist(picture(i).ch1_filename) & ~exist(picture(i).ch2_filename)
        disp(strcat('File #', int2str(i), ' missing! Try AGAIN'));
        return;
    end
end
picname=red_file_name;


%Now to open each pic and:
% 1) Locate Cells
% 2) For each image define a background region for background subtraction
% 2) Open a Save File if it has been previously analyzed
% 3) Do the analysis which includes:
%  3a) Input the number of cells in each image
%  3b) Crop each cell in each image
%  3c) Determine Cell Boundary with user clicks
%  3d) Determine nuclear membrane with user clicks
%  3e) Locate NPs using pkfnd and cntrd (from Eric Weeks)
% 4) Average results over multiple cells

for i=first_pic:lastpic
    backgroundF(i)= definebackground(picture(i).ch1_filename, picture(i).ch2_filename);
    
    
    if exist(strcat(picname, int2str(i), 'data.txt'), 'file')
        %load file and read how many cells were defined in image
        fileID=fopen(strcat(picname, int2str(i), 'data.txt'), 'r');
        formatspec = '%f';
        num_cells(i)=fscanf(fileID, formatspec);
        fclose(fileID);
        
        for j=1:num_cells(i)
            
            %Read file with cornerpoints and crop main image
            
            picture(i).cell(j).cornerfile=fopen(strcat(picname, int2str(i), '_boundarypoints', int2str(j), '.txt'), 'r');
            picture(i).cell(j).corners=fscanf(picture(i).cell(j).cornerfile, '%f ');
            fclose(picture(i).cell(j).cornerfile);
            [picture(i).cell(j).red, picture(i).cell(j).bright]=cropimageFP_w_points(picture(i).ch1_filename, picture(i).ch2_filename, ...
                picture(i).cell(j).corners(1), picture(i).cell(j).corners(2), picture(i).cell(j).corners(3), picture(i).cell(j).corners(4));
            
            %Read file with boundary
            picture(i).cell(j).boundaryfile=fopen(strcat(picname, int2str(i),'_', int2str(j), 'boundary', '.txt'), 'r');
            picture(i).cell(j).boundary_xy=fscanf(picture(i).cell(j).boundaryfile,'%f ');
            fclose(picture(i).cell(j).boundaryfile);
            
            picture(i).cell(j).boundary_x=transpose(picture(i).cell(j).boundary_xy(1:2:end));
            picture(i).cell(j).boundary_y=transpose(picture(i).cell(j).boundary_xy(2:2:end));
            
            [picture(i).cell(j).tot_NPs_in_region, picture(i).cell(j).tot_F_in_region, picture(i).cell(j).tot_pix_in_region, picture(i).cell(j).final_image]=...
                conc5_w_saved_info(picture(i).cell(j).red, picture(i).cell(j).bright, picture(i).cell(j).boundary_x, picture(i).cell(j).boundary_y, backgroundF(i), ellipse_spacing, interactive_mode, red_file_name, j);
    
        end   
            
     else
        %File Does Not Exist
    
        num_cells(i)=usercellcount(picture(i).ch2_filename);
        
   
        for j=1:num_cells(i)
            [picture(i).cell(j).red, picture(i).cell(j).bright, picture(i).cell(j).corners]=crop_imageFP(picture(i).ch1_filename, picture(i).ch2_filename, j);
        end
    
        for j=1:num_cells(i)
            [picture(i).cell(j).tot_NPs_in_region, picture(i).cell(j).tot_F_in_region, picture(i).cell(j).tot_pix_in_region, picture(i).cell(j).boundary_x, picture(i).cell(j).boundary_y, picture(i).cell(j).final_image]...
                =conc5(picture(i).cell(j).red, picture(i).cell(j).bright, backgroundF(i), interactive_mode, ellipse_spacing, [red_file_name, int2str(i)], j);
            
            %Write File with Corners of cell
            picture(i).cell(j).cornerfile=fopen(strcat(picname, int2str(i), '_boundarypoints', int2str(j), '.txt'), 'w');
            fprintf(picture(i).cell(j).cornerfile, '%f ', picture(i).cell(j).corners);
            fclose(picture(i).cell(j).cornerfile);
            
            %Write Boundary file
            picture(i).cell(j).boundaryfile=fopen(strcat(picname, int2str(i),'_', int2str(j), 'boundary', '.txt'), 'w');
            fprintf(picture(i).cell(j).boundaryfile, '%f ', [picture(i).cell(j).boundary_x; picture(i).cell(j).boundary_y]);
            fclose(picture(i).cell(j).boundaryfile);
        end
        writesavefile=fopen(strcat(picname, int2str(i), 'data.txt'), 'w');
        fprintf(writesavefile, '%f', num_cells(i));
        fclose(writesavefile);    
    end

end


%Everything Below Here is about averaging the results from multiple cells
%and outputting it both the command window as a structure named output as
%well as a text file for easy copy-pasting

%Find the biggest cell
max_num_regions=1;
k=0;
for i=first_pic:lastpic
    for j=1:num_cells(i)
         num_regions_cell(i, j)=size(picture(i).cell(j).tot_NPs_in_region, 2);
         if num_regions_cell(i, j) > max_num_regions
             max_num_regions=num_regions_cell(i, j);
         end
    end
end

all_NPs_in_region=zeros(max_num_regions, 1);
all_F_in_region=zeros(max_num_regions, 1);
all_pix_in_region=zeros(max_num_regions, 1);
num_regions_counted=zeros(max_num_regions, 1);
h=0;

%Add up all the NPs, fluorescence and pixels in each region and how many cells contain that
%region
for i=first_pic:lastpic
    for j=1:num_cells(i)
        for k=1:num_regions_cell(i, j)
            all_NPs_in_region(k)=all_NPs_in_region(k)+picture(i).cell(j).tot_NPs_in_region(k);
            
            all_F_in_region(k)=all_F_in_region(k)+(picture(i).cell(j).tot_F_in_region(k));
            
            all_pix_in_region(k)=all_pix_in_region(k)+picture(i).cell(j).tot_pix_in_region(k);
            
            num_regions_counted(k)=num_regions_counted(k)+1;
        end
    end
end



%Take average

average_NPs_region=all_NPs_in_region./num_regions_counted;
average_F_region=all_F_in_region./num_regions_counted;
average_pix_region=all_pix_in_region./num_regions_counted;
average_F_perpix_region=all_F_in_region./all_pix_in_region;
average_NP_perpix_region=all_NPs_in_region./all_pix_in_region;


%Find Standard Deviation

STD_NPs_in_region=zeros(max_num_regions, 1);
STD_F_in_region=zeros(max_num_regions, 1);
STD_pix_in_region=zeros(max_num_regions, 1);


for i=first_pic:lastpic
    for j=1:num_cells(i)
        %num_regions_cell=size(picture(i).cell(j).tot_NPs_in_region, 2);
        for k=1:num_regions_cell(i, j)
            STD_NPs_in_region(k)=STD_NPs_in_region(k)+(average_NPs_region(k)-picture(i).cell(j).tot_NPs_in_region(k))^2;
            
            STD_F_in_region(k)=STD_F_in_region(k)+(average_F_region(k)-(picture(i).cell(j).tot_F_in_region(k)))^2;
            
            STD_pix_in_region(k)=STD_pix_in_region(k)+(average_pix_region(k)-picture(i).cell(j).tot_pix_in_region(k))^2;
            
        end
    end
end

STD_NPs_in_region=sqrt(STD_NPs_in_region./(num_regions_counted.*(num_regions_counted-1)));
STD_F_in_region=sqrt(STD_F_in_region./(num_regions_counted.*(num_regions_counted-1)));
STD_pix_in_region=sqrt(STD_pix_in_region./(num_regions_counted.*(num_regions_counted-1)));


%Find Error for NP or Fluorescence per pixel 
ERROR_NPs=STD_NPs_in_region;
ERROR_F_in_region=STD_F_in_region;
ERROR_F_perpix_region=average_F_perpix_region.*sqrt((STD_F_in_region./average_F_region).^2+(STD_pix_in_region./average_pix_region).^2);
ERROR_NP_perpix_region=average_NP_perpix_region.*sqrt((STD_NPs_in_region./average_NPs_region).^2+(STD_pix_in_region./average_pix_region).^2);



output.average_NPs_region=transpose(average_NPs_region);
output.ERROR_NPs=(ERROR_NPs);
output.average_F_region=transpose(average_F_region);
output.ERROR_F_in_region=(ERROR_F_in_region);
output.average_F_perpix_region=transpose(average_F_perpix_region);
output.ERROR_F_perpix_region=transpose(ERROR_F_perpix_region);
output.average_NP_perpix_region=1000*transpose(average_NP_perpix_region);
output.ERROR_NP_perpix_region=1000*transpose(ERROR_NP_perpix_region);
tot_cells=sum(num_cells(:));
%Save Data in a Text File
save_file=fopen([red_file_name, '_results.txt'], 'w');
fprintf(save_file, 'Number of Cells Measured: %f\n\n', tot_cells);
fprintf(save_file, 'Average Number of NPs per Region: \n');
fprintf(save_file, '%f\n', output.average_NPs_region);
fprintf(save_file, '\nError:\n');
fprintf(save_file, '%f\n', output.ERROR_NPs);
fprintf(save_file, '\nAverage Fluorescence Intensity Per Region:\n');
fprintf(save_file, '%f\n', output.average_F_region);
fprintf(save_file, '\nError:\n');
fprintf(save_file, '%f\n', output.ERROR_F_in_region);
fprintf(save_file, '\nAverage NPs per region per Pixel:\n');
fprintf(save_file, '%f\n', output.average_NP_perpix_region);
fprintf(save_file, '\nError:\n');
fprintf(save_file, '%f\n', output.ERROR_NP_perpix_region);
fprintf(save_file, '\nAverage Fluorescence Intensity Per Region Per Pixel:\n');
fprintf(save_file, '%f\n', output.average_F_perpix_region);
fprintf(save_file, '\nError:\n');
fprintf(save_file, '%f\n', output.ERROR_F_perpix_region);
fclose(save_file);


end

% This function prompts the user to click on opposing corners of a blank
% region in the fluorescent channel. This allows the program to calculate a
% background value of fluorescence and subtract it from fluorescence
% measurements later.

function average_F_per_pix=definebackground(image1_filename, image2_filename)
    
    image1=imread(image1_filename);
    image2=imread(image2_filename);
    
    subplot(1, 2, 1), subimage(image1);
    subplot(1, 2, 2), subimage(image2);
    disp('Click on corners of background area for normalization')
    p=ginput(2);
    
    sp(1)=min(floor(p(1)), floor(p(2)));
    sp(2)=min(floor(p(3)), floor(p(4)));
    sp(3)=max(ceil(p(1)), ceil(p(2)));
    sp(4)=max(ceil(p(3)), ceil(p(4)));
    
    ch1_region=image1(sp(2):sp(4), sp(1):sp(3));
    
    total_F=sum(sum(ch1_region));
    
    num_pix=abs((sp(3)-sp(1))*(sp(4)-sp(2)));
    
    average_F_per_pix=total_F/num_pix;
    close all
end

%This function opens an image and prompts the user to identify
%how many cells are in the image


function no_cells=usercellcount(filename)
a=imread(filename);
imagesc(a);
colormap('gray');

question1=strcat('How many cells in this image?');
no_cells=input(question1);

end


% This function takes in two image filenames (ch1_filename and ch2_filename) as well as
% the current cell number (cell_num) and asks the usert to click the
% corners of cell for cropping purposes. It outputs the cropped arrays.
% It also outputs the corner of the images for back-up purposes.

function [ch1_region, ch2_region, b]=crop_imageFP(ch1_filename, ch2_filename, cell_num)
    ch1_image=imread(ch1_filename);
    ch2_image=imread(ch2_filename);
    [mx_row, mx_col]=size(ch2_image);
    figure;
    imagesc(ch2_image); axis image
    disp(['Click on Corners of Cell #', int2str(cell_num)])
    p=ginput(2);
    
    sp(1)=min(floor(p(1)), floor(p(2)));
    sp(2)=min(floor(p(3)), floor(p(4)));
    sp(3)=max(ceil(p(1)), ceil(p(2)));
    sp(4)=max(ceil(p(3)), ceil(p(4)));
    
    sp(2)=max(1, sp(2));
    sp(1)=max(1, sp(1));
    sp(3)=min(mx_col, sp(3));
    sp(4)=min(mx_row, sp(4));
    
    ch1_region=ch1_image(sp(2):sp(4), sp(1):sp(3));
    ch2_region=ch2_image(sp(2):sp(4), sp(1):sp(3));
    b=[sp(2), sp(4), sp(1), sp(3)];
    
end
  
% When a file is being analyzed again this version of crop image is used
% which just re-uses the old corners

function [ch1_region, ch2_region]=cropimageFP_w_points(ch1_filename, ch2_filename, sp2, sp4, sp1, sp3)

    ch1=imread(ch1_filename);
    ch2=imread(ch2_filename);

    ch1_region=ch1(sp2:sp4, sp1:sp3);
    ch2_region=ch2(sp2:sp4, sp1:sp3);
  
end






% This program Calculates the concentration of particles and fluorescence 
% intesnity as a function of distance from the nucleus center. 
% The inputs are: 
% redregion - Image array of fluorescent channel
% brightregion - Image array of  bright channel 
% (for cell and nuclear boundary determination)
% backgroundF - A constant representing the background fluorescence value
% Ellipse Spacing - The size of the regions of the cell that will be
% defined
% The outputs are:
% tot_NPs_in_region - Total Number of particles in each region
% tot_F_in_region - Total Fluorescence in each reigon
% tot_pix_in_region - Total Pixels in each region
% boundary_x, boundary_y - x,y values of cell boundary (will be written to
% a save file)
% final_image - An image showing the results of the analysis
% This program requires heavy user input, the user will define the cell
% and nuclear boundary by clicking, and se the parameters for locating
% particles

function [tot_NPs_in_region, tot_F_in_region, tot_pix_in_region, boundary_x, boundary_y]=conc5(redregion, brightregion, backgroundF, interactive_mode, ellipse_spacing, filename, cell_num)
        [height, width]=size(redregion);
        outline_logical='n';
        
        while outline_logical ~= 'y'
            [nuc_bound_x, nuc_bound_y, boundary_x, boundary_y]=boundary_maker(redregion, brightregion, ellipse_spacing);
        
            imshow(brightregion);
            colormap('gray');
            hold on
        
            plot(boundary_x, boundary_y);
            hold on
            plot(nuc_bound_x(:, 1), nuc_bound_y(:, 1));
            hold off
        
            outline_logical=input('Happy with Nucleus and Boundary Outline? (y/n)', 's');
        end
        
        
        %Make the Cell Boundary into a mask
        boundarymask=poly2mask(boundary_x, boundary_y, height, width);
        boundarymaskedgreen_region=(boundarymask).*double(redregion);

        
        
        %First Find all Particles in Image (user will set parameters)
        [pkheight_NP, part_size_NP]=paramfinder2(boundarymaskedgreen_region);
        pk_NP = pkfnd(boundarymaskedgreen_region,pkheight_NP,part_size_NP);
        cntb_NP = cntrd(boundarymaskedgreen_region,pk_NP,part_size_NP);
        
        in_cell_x=cntb_NP(:, 1);
        in_cell_y=cntb_NP(:, 2);
        in_cell_size=cntb_NP(:, 4);

        
        
        num_NP_in_cell=size(in_cell_x);
        
        %Get rid of ones that are too close to the boundary
        
        for i=1:num_NP_in_cell
            distance_to_bound=min_dist(in_cell_x(i), in_cell_y(i), boundary_x, boundary_y);
            if distance_to_bound < 3
                tooclose(i)=0;
            else
                tooclose(i)=1;
            end
        end
        if size(in_cell_x)~=0
            tooclose=logical(tooclose);
            in_cell_x=in_cell_x(tooclose);
            in_cell_y=in_cell_y(tooclose);
        elseif size(x)==0
            x=0;
            y=0;
        end
        
        reg_num=size(nuc_bound_x, 2);
        %Measure NPs and Fluorescence in each region of the cell
     
        for i=1:reg_num-1
            innermask=poly2mask(nuc_bound_x(:, i), nuc_bound_y(:, i), height, width);
          
            outermask=poly2mask(nuc_bound_x(:, i+1), nuc_bound_y(:, i+1), height, width);

            full_regionmask=outermask-innermask;

            regionmasked_picture=(boundarymaskedgreen_region).*full_regionmask;


            region(i).pk_NP = pkfnd(regionmasked_picture,pkheight_NP,part_size_NP);
            if size(region(i).pk_NP, 1) == 0
                tot_NPs_in_region(i)=0;
            else
                region(i).cntb_NP = cntrd(regionmasked_picture,region(i).pk_NP,part_size_NP);
                tot_NPs_in_region(i)=max(0, size(region(i).cntb_NP, 1));
            end
            
            normalized_regionmasked_picture=regionmasked_picture-backgroundF;
            negatives= normalized_regionmasked_picture < 0;
        
            normalized_regionmasked_picture(negatives) = 0;
            
            tot_F_in_region(i)=sum(sum(normalized_regionmasked_picture));
            tot_pix_in_region(i)=sum(sum(full_regionmask.*boundarymask));
            
        end
            

        % Make a pretty picture with results

        if interactive_mode==1
            figure;
            imshow(redregion);
            colormap('gray');
            hold on;
            plot(boundary_x, boundary_y, '-w')
            plot(in_cell_x, in_cell_y, '+r')
            plot(nuc_bound_x, nuc_bound_y);
            hold off;
            pause(5);
            F=getFrame;
            close all;
            [b_mask_rows, b_mask_col]=size(boundarymask);
            final_image(:, :, 1)=uint8(boundarymask).*F.cdata(1:b_mask_rows, 1:b_mask_col, 1);
            final_image(:, :, 2)=uint8(boundarymask).*F.cdata(1:b_mask_rows, 1:b_mask_col, 2);
            final_image(:, :, 3)=uint8(boundarymask).*F.cdata(1:b_mask_rows, 1:b_mask_col, 3);
            frame=F.cdata;
            imshow(final_image, []);
            save_file=getFrame;
            imwrite(save_file.cdata, [filename, '_cell_', int2str(cell_num), '.tif']);
        end
        


end



%This function makes the cell boundary, the nuclear outline and defines
%each region from the cell
function [nuc_bound_x, nuc_bound_y, boundary_x, boundary_y]=boundary_maker(redregion, brightregion, ellipse_spacing)
    


        % Display red channel with particles marked
        subplot(1, 2, 1), imshow(redregion)
        subplot(1, 2, 2), imshow(brightregion)
        pkheight=10;
        part_size=3;
        [height, width]=size(redregion);

    
        % User picks which image to use for boundary determination
        usebright=input('Which image for boundary determination? (b/r) ', 's');
        
        if usebright=='r'
            c = redregion;
        elseif usebright=='b'
            c = brightregion;
        end
        close all;

        
        
        %Optimize boundary_determination image for locating boundaries
        
        imagesc(c);
        colormap('gray');
        axis image;
        %[pkheight2, part_size]=paramfinder3(c);
        pkheight2=10;
        part_size=3;
        pk = pkfnd(double(c),pkheight2,part_size);
        cntb = cntrd(double(c),pk,part_size);
        
        
        
        clickingsides='y';
        num_sides=0;
        
       % Let the user click on sides and indicate when he/she is finsihed
       
        while(clickingsides=='y')
            disp('Click on points of a side and press enter when done')
            num_sides=num_sides+1;
            [side(num_sides).x, side(num_sides).y]=ginput;
            clickingsides=input('More sides to define? (y/n)', 's');
            side(num_sides).num_points=size(side(num_sides).x, 1);
        end
        
        
   
        %Exchange the user clicks for actualy points found on the side of
        %the cell
        firstpoint_index=findclosest(side(1).x(1), side(1).y(1), cntb(:, 1:2));
        side(1).firstpoint_x=cntb(firstpoint_index, 1);
        side(1).firstpoint_y=cntb(firstpoint_index, 2);
        
        
        for i=1:num_sides
            for k=1:side(i).num_points
                    thepoint=findclosest(side(i).x(k), side(i).y(k), cntb(:, 1:2));
                    side(i).xy(k, 1)=cntb(thepoint, 1);
                    side(i).xy(k, 2)=cntb(thepoint, 2);
            end
        end
   
        %Rotate, fit and unrotate each side
    
        for j=1:num_sides
            [side(j).fitx, side(j).fity]=rotate_and_fit(side(j).xy(:, 1), side(j).xy(:, 2));
        end

        %Put all sides into one big array
        
        if r(side(1).firstpoint_x-side(1).fitx(1), side(1).firstpoint_y-side(1).fity(1)) < r(side(1).firstpoint_x-side(1).fitx(end), side(1).firstpoint_y-side(1).fity(end)) 
                boundary_x=side(1).fitx;
                boundary_y=side(1).fity;
        else
              boundary_x=flipdim(side(1).fitx, 2);
              boundary_y=flipdim(side(1).fity, 2);
        end
        old_index=size(side(1).fitx, 2);
        for i=2:num_sides
            distance_last_to_first=r(boundary_x(end)-side(i).fitx(1), boundary_y(end)-side(i).fity(1));
            distance_last_to_last=r(boundary_x(end)-side(i).fitx(end), boundary_y(end)-side(i).fity(end));
            if distance_last_to_first < distance_last_to_last 
                boundary_x(old_index:(old_index+size(side(i).fitx, 2)-1))=side(i).fitx;
                boundary_y(old_index:(old_index+size(side(i).fitx, 2)-1))=side(i).fity;
            elseif distance_last_to_first > distance_last_to_last 
                boundary_x(old_index:(old_index+size(side(i).fitx, 2)-1))=flipdim(side(i).fitx, 2);
                boundary_y(old_index:(old_index+size(side(i).fitx, 2)-1))=flipdim(side(i).fity, 2);
            end
            old_index=old_index+size(side(i).fitx, 2);
        end
        


        boundary_x(end)=boundary_x(1);
        boundary_y(end)=boundary_y(1);

        
        [nuc_bound_x, nuc_bound_y, nuc_cent_x, nuc_cent_y, reg_num, phi, x_axis, y_axis]=nucleusfinder_regionmaker2(brightregion, boundary_x, boundary_y, ellipse_spacing);

end

%This function prompts the user to click around the nuclear membrane and
%then fits the user selected points to an ellipse, next it will generate
%equidistant ellipses from the nuclear membrane

function [nuc_bound_x, nuc_bound_y, x_cent, y_cent, num_ellipse, phi, x_axis, y_axis]=nucleusfinder_regionmaker2(bright_array, boundary_x, boundary_y, spacing)
    
    %How big the regions will be
   
    %Read in Image and prompt user to click around nucleus
    
    
    imagesc(bright_array);
    colormap('gray');
    disp('Click on some points to define the nucleus')
    [x, y]=ginput;
    xy(:, 1)=x;
    xy(:, 2)=y;
    
    a=fit_ellipse(x, y);
    % Output:   ellipse_t - structure that defines the best fit to an ellipse
    %                       a           - sub axis (radius) of the X axis of the non-tilt ellipse
    %                       b           - sub axis (radius) of the Y axis of the non-tilt ellipse
    %                       phi         - orientation in radians of the ellipse (tilt)
    %                       X0          - center at the X axis of the non-tilt ellipse
    %                       Y0          - center at the Y axis of the non-tilt ellipse
    %                       X0_in       - center at the X axis of the tilted ellipse
    %                       Y0_in       - center at the Y axis of the tilted ellipse
    %                       long_axis   - size of the long axis of the ellipse
    %                       short_axis  - size of the short axis of the ellipse
    %                       status      - status of detection of an ellipse
        num_boundary_points=size(boundary_x, 2);
        maxr=1;
        for h=1:num_boundary_points
            radialsize(h)=r(a.X0_in-boundary_x(h), a.Y0_in-boundary_y(h));
            if radialsize(h)> maxr
                 maxr=radialsize(h);
                 max_bound_index=h;
                
            end
        end
       

       
       num_ellipse=floor(maxr/spacing);
    
   
    for k=1:num_ellipse
        major=a.a+(k-1)*spacing;
        minor=a.b+(k-1)*spacing;
    
        % initial points for the non-tilt ellipse
        tempx(1, k)=-(major);
        tempy(1, k)=sqrt(minor^2*(1-(tempx(1, k)^2)/major^2));
        tempx(100, k)=-(major);
        tempy(100, k)=-sqrt(minor^2*(1-(tempx(1)^2)/major^2));
        % Sweeping across, we calculate both y values
        for i=2:101
            tempx(i, k)=tempx(1, k)+(((i-1)/100)*2*major);
            tempx(202-i, k)=tempx(1, k)+(((i-1)/100)*2*major);
            tempy(i, k)=sqrt(minor^2*(1-(tempx(i, k)^2)/major^2));
            tempy(202-i, k)=-sqrt(minor^2*(1-(tempx(i, k)^2)/major^2));
        end
    
  
  
    
        
        for i=1:200
            [nuc_bound_x(i, k), nuc_bound_y(i, k)]=rotate(tempx(i, k), tempy(i, k), -a.phi);
            nuc_bound_x(i, k)=nuc_bound_x(i, k)+a.X0_in;
            nuc_bound_y(i, k)=nuc_bound_y(i, k)+a.Y0_in;
        end
        nuc_bound_x(200, k)=nuc_bound_x(1, k);
        nuc_bound_y(200, k)=nuc_bound_y(1, k);
    end
    
    %Plot it on the image

    x_cent=a.X0_in;
    y_cent=a.Y0_in;
    phi=a.phi;
    x_axis=a.a;
    y_axis=a.b;
end

%This function takes a list of x-y coordinates (input_x, input_y), rotates
%them at different angles, fits the rotated data, determines which angle
%gives the best fit and returns (fit_x, fit_y) which is the fitted data
%after rotation.

function [fit_x, fit_y ] = rotate_and_fit(input_x, input_y )

%Different angles to try
phi(1:8)=pi/4*[0:7];
%Turn polyfit warning off
ws=warning('off', 'all');

%Rotate each all points on the side to the different angles and fit a
%polynomial
for i=1:8
    for j=1:size(input_x, 1)
        [angle(i).rotated(j, 1), angle(i).rotated(j, 2)]=rotate(input_x(j, 1), input_y(j, 1), phi(i));
    end
    
    [angle(i).fit(1:5), angle(i).Error]=polyfit(angle(i).rotated(:, 1), angle(i).rotated(:, 2), 4);
    [angle(i).trash, angle(i).delta]=polyval(angle(i).fit(1:5), angle(i).rotated(:, 1), angle(i).Error);
    
end
%Find which angle gave the best polynomial fit
best_angle=1;
for i=1:8
    if angle(i).delta < angle(best_angle).delta
        best_angle=i;
    end
end


%Get the values of the fit

min_x=min(angle(best_angle).rotated(:, 1));
max_x=max(angle(best_angle).rotated(:, 1));
fine_x=[min(min_x, max_x): abs(min_x-max_x)/100 :max_x];
bestfit_y=polyval(angle(best_angle).fit(1:5), fine_x);

%Rotate them back
for i=1:size(fine_x, 2);
    [fit_x(i), fit_y(i)]=rotate(fine_x(i), bestfit_y(i), -phi(best_angle));
end

warning(ws);
end

%This program reads in an image and then shows the user what particles it would
%identify in the image. The user is allowed to change the peak finding parameters until
%the proper number of particles are found


function [pkheight, partsize, min_radius]=paramfinder2(image_array)

a=image_array;


imagesc(a);
colormap('gray');

maxthresh=max(max(a));
partsize=5;
pkheight=floor(.8*maxthresh);
min_radius=0;
pks=pkfnd(a, pkheight, partsize);
cnt=cntrd(a, pks, partsize);

    too_small= cnt(:, 4) < min_radius;
    cnt(too_small, :)=[];

        
hold on
plot(cnt(:, 1), cnt(:, 2), '+y')
hold off
good='n';
while good~= 'y'
    question1=strcat('Threshold currently at ', int2str(pkheight), ' New Threshold? (max at', int2str(maxthresh), ')');
    pkheight=input(question1);
    if pkheight > maxthresh
        question1=strcat('Threshold too high! ', ' New Threshold? (max at', int2str(maxthresh), ')');
        pkheight=input(question1);
    end
    question2=strcat('Particle size currently at  ', int2str(partsize), ' New Partsize?');
    partsize=input(question2);
    question3=strcat('Minimum Particle Radius is currently at ', int2str(min_radius), ' New Minimum Radius?');
    min_radius=input(question3);
    pks=pkfnd(a, pkheight, partsize);
    cnt=cntrd(a, pks, partsize);
    too_small= cnt(:, 4) < min_radius;
    cnt(too_small, :)=[]; 
    
    imagesc(a)
    hold on
    plot(cnt(:, 1), cnt(:, 2), '+y')
    hold off
    good=input('Good? (y/n)', 's');
end
close all;
end


% Given an x, y coordiante this function determines r

function length =r(x, y)
    length = sqrt((x^2)+(y^2));
end


% This is like conc5 above expcept it uses saved info

function [tot_NPs_in_region, tot_F_in_region, tot_pix_in_region]=conc5_w_saved_info(redregion, brightregion, boundary_x, boundary_y, backgroundF, ellipse_spacing, interactive_mode, filename, cell_num)
    
    % Start with Boundary
    
        a=redregion;
        c=brightregion;
        [height, width]=size(a);
        boundary_x(size(boundary_x, 2))=boundary_x(1);
        boundary_y(size(boundary_y, 2))=boundary_y(1);
        boundarymask=poly2mask(boundary_x, boundary_y, height, width);
        boundarymaskedgreen_region=boundarymask.*double(redregion);
        

        
        
        [nuc_bound_x, nuc_bound_y, nuc_cent_x, nuc_cent_y, reg_num, phi, x_axis, y_axis]=nucleusfinder_regionmaker2(brightregion, boundary_x, boundary_y, ellipse_spacing);
        imagesc(c)
        colormap('gray');
        hold on
        
        plot(boundary_x, boundary_y);
        hold on
        
        plot(nuc_bound_x(:, 1), nuc_bound_y(:, 1))
       
        hold off
        
        better=input('Happy with boundary and Nucleus?', 's');
        
        %Now to find the particles that are inside the cell
        
        %First Find all Particles in Image
        [pkheight_NP, part_size_NP]=paramfinder2(boundarymaskedgreen_region);
        pk_NP = pkfnd(boundarymaskedgreen_region,pkheight_NP,part_size_NP);
        cntb_NP = cntrd(boundarymaskedgreen_region,pk_NP,part_size_NP);
        tot_part_in_cell=size(boundarymaskedgreen_region);
        
        
        in_cell_x=cntb_NP(:, 1);
        in_cell_y=cntb_NP(:, 2);
        in_cell_size=cntb_NP(:, 4);

        
        num_NP_in_cell=size(in_cell_x);
        
        %Get rid of ones that are too close to the boundary
        
        for i=1:num_NP_in_cell
            distance_to_bound=min_dist(in_cell_x(i), in_cell_y(i), boundary_x, boundary_y);
            if distance_to_bound < 3
                tooclose(i)=0;
            else
                tooclose(i)=1;
            end
        end
        if size(in_cell_x)~=0
            tooclose=logical(tooclose);
            size(in_cell_x);
            size(tooclose);
            in_cell_x=in_cell_x(tooclose);
            in_cell_y=in_cell_y(tooclose);
            in_cell_size=in_cell_size(tooclose);
        elseif size(x)==0
            x=0;
            y=0;
        end
        
       
        
        %Find out how many are in each region of nucleus
        max_size=1;
        for i=1:reg_num-1
            innermask=poly2mask(nuc_bound_x(:, i), nuc_bound_y(:, i), height, width);

            outermask=poly2mask(nuc_bound_x(:, i+1), nuc_bound_y(:, i+1), height, width);

            full_regionmask=outermask-innermask;

            regionmasked_picture=(boundarymaskedgreen_region).*full_regionmask;
            

            region(i).pk_NP = pkfnd(regionmasked_picture,pkheight_NP,part_size_NP);
            if size(region(i).pk_NP, 1) == 0
                tot_NPs_in_region(i)=0;
            else
                region(i).cntb_NP = cntrd(regionmasked_picture,region(i).pk_NP,part_size_NP);
                tot_NPs_in_region(i)=max(0, size(region(i).cntb_NP, 1));
            end
            normalized_regionmasked_picture=regionmasked_picture-backgroundF;
            negatives= normalized_regionmasked_picture < 0;
        
            normalized_regionmasked_picture(negatives) = 0;
            
            tot_F_in_region(i)=sum(sum(normalized_regionmasked_picture));
            tot_pix_in_region(i)=sum(sum(full_regionmask.*boundarymask));
            
        end
        close all;
        if interactive_mode==1
            figure;
            imshow(redregion);
            colormap('gray');
            hold on;
            plot(boundary_x, boundary_y, '-w')
            plot(in_cell_x, in_cell_y, '+r')
            plot(nuc_bound_x, nuc_bound_y);
            hold off;
            pause(5);
            F=getFrame;
            [b_mask_rows, b_mask_col]=size(boundarymask);
            final_image(:, :, 1)=uint8(boundarymask).*F.cdata(1:b_mask_rows, 1:b_mask_col, 1);
            final_image(:, :, 2)=uint8(boundarymask).*F.cdata(1:b_mask_rows, 1:b_mask_col, 2);
            final_image(:, :, 3)=uint8(boundarymask).*F.cdata(1:b_mask_rows, 1:b_mask_col, 3);
            close all;
            imshow(final_image, []);
            pause(5);
            save_file=getFrame;
            imwrite(save_file.cdata, [filename, '_cell_', int2str(cell_num), '.tif']);
        end


end

%This function rotates points (x, y) by theta and outputs the (newx, newy)

function [newx, newy]=rotate(x, y, theta)
    newx=x*cos(theta)-y*sin(theta);
    newy=x*sin(theta)+y*cos(theta);
end


%Given an x, y coordinate and a list of other coordinates this program
%determines which coordinates is closes to x, y. (This can be optimized for
%matlab)


function  closest = findclosest(x, y, points)

number_points=size(points);
min=r(x-points(1,1), y-points(1,2));
closest = 1;
for i=1:number_points(1)
    if r(x-points(i,1), y-points(i,2)) < min
        min=r(x-points(i,1), y-points(i,2));
        closest=i;
    end
end
end

%This function finds the shortest distance between a point (x, y) and a
%list of xy points (boundary_x, boundary_y)(This can be optimized for
%matlab)

function shortest=min_dist(x, y, boundary_x, boundary_y)

boundsize=size(boundary_x);
r_min=200;

for i=1:boundsize(2)
    r=((boundary_x(i)-x)^2)+((boundary_y(i)-y)^2);
    
    if r < r_min
        r_min=r;
    end
end
shortest=r_min;
end

