function [grid_image,corners_x,corners_y,bestangle]=make_grid_v4(gridimage_filename,format)

%Fred Cross June 29 2016. The Rockefeller University

%make a grid with specified format given a digital image of a grid
%alignment plate

global ninespots;
bestangle=0;


grid_image=imread(gridimage_filename);


figure(101);imagesc(grid_image);

ninespots=get_ninespots(grid_image);

[corners_x,corners_y]=combine_four_grids(format);

[corners_x,corners_y,bestangle]=adjust_corners(corners_x,corners_y,grid_image,format);




%%%%%%%%%%%%%%%%%
function ninespots=get_ninespots(input_grid_image)
%given new image with 3 dots at top, mid, bottom of micro plate, find their
%x,y coords ; save in ninespots 3x3 cell array, each entry [x,y]
%%
%grid_image=imread('PC141732.JPG');



%new version for red spots

grid_image_p=255-input_grid_image;
grid_image=double(grid_image_p(:,:,3)-grid_image_p(:,:,1));
grid_image(grid_image<0)=0;
[n,xout]=hist(grid_image(:),0:255);
cutoff=xout(min(find(cumsum(n)/sum(n)>=0.995))); % based on approx rep of red spots in image
grid_image(grid_image<cutoff)=0;
grid_image=grid_image/max(grid_image(:));

% figure(1);imagesc(grid_image);
% set(gcf,'Position',[ 324          83        1474         892]);
% 
% 
% title(sprintf('Click on middle spot','FontSize',10));
% 
% [x,y,b]=ginput(1);
% 
% 
% 
% grid_image=input_grid_image(:,:,1);
% grid_image=(255-double(grid_image))/255;
grid_image_t=grid_image;  %imtophat(grid_image,strel('disk',5));
% figure(1001);
% imagesc(grid_image_t);
% try
% grid_image_t(round(y),round(x))
% grid_image_t(round(y/2),round(x/2))
% grid_image_t(round(3*y/2),round(3*x/2))
% grid_image_t(round(y/2),round(3*x/2))
% grid_image_t(round(3*y/2),round(x/2))
% impixelinfo
% pause
% catch
% end
% end


%clear im;
figure;imagesc(grid_image_t);shg;impixelinfo
%[n,xout]=hist(grid_image_t(:),10000);
%thresh=xout(min(find(cumsum(n)/sum(n)>0.997)));
found_thresh=0;

%%%  automatic threshold finder. A bit clunky but seems to work well and
%%%  quickly enough
low_thresh=0;
high_thresh=1;
ntries=0;
while 1 %~found_thresh
    ntries=ntries+1;
    if ntries>100
        gi_temp=grid_image_t;
        gi_temp(gi_temp<0.2)=0;
        break
    end
    
    thresh=(low_thresh+high_thresh)/2;
    %%
    %thresh=xout(min(find(cumsum(n)/sum(n)>cutoff)));

    gi_temp=grid_image_t;
    %dump dim stuff
    gi_temp(gi_temp<thresh)=0;
    %%%%%
     BW=gi_temp>0;
     %dump small stuff
     BWO=bwareaopen(BW,1000);
    
    L=bwlabel(BWO);
    stats=regionprops(BWO,'Centroid','Area','Eccentricity');
    %dump elongated stuff
    for k=1:length(stats)
        if stats(k).Eccentricity>0.9
            L(L==k)=0;
        end
    end
    gi_temp=L>0;
    stats=regionprops(gi_temp,'Centroid','Area','Eccentricity');
   % imagesc(gi_temp);shg
  %binary search for threshold giving just 9 objects
    if length(stats)<9 %cut out too much
        high_thresh=(thresh+high_thresh)/2;
      %  thresh=thresh/(1+rand/10);
    elseif length(stats)>9 %too many odd bits
        low_thresh=(thresh+low_thresh)/2;
        %thresh=thresh*(1+rand/10);
    
    else 
        break
    end
    
end


   %%%%%
   
   while ~found_thresh
    
    figure;imagesc(gi_temp);shg;impixelinfo
    q=menu('Accept threshold?','y','n');
    if q==1
        found_thresh=1;
        grid_image_t=gi_temp;
    else
        gi_temp=grid_image_t;
        answer=inputdlg({'threshold'},'',1,{num2str(thresh)});
        thresh=str2num(answer{1});
    
    end
end
%%
found_pts=0;
Area_crit=1000;
Ecc_crit=0.9;
while ~found_pts
    L=bwlabel(grid_image_t);
   % imagesc(L);shg
    stats=regionprops(L>0,'Centroid','Area','Eccentricity');
    
    for k=1:length(stats)
        if stats(k).Area<Area_crit | stats(k).Eccentricity>Ecc_crit
            L(L==k)=0;
        end
    end
    
    figure;imagesc(L>0);
    found_pts=0;
%     q=menu('Accept anchor points?','y','n');
%     if q==1
        found_pts=1;
        grid_image_cleaned=double(L>0).*grid_image_t;
        L=bwlabel(grid_image_cleaned);
        stats=regionprops(L>0,'Centroid','Area','Eccentricity');
        clear grid_image L;
        
%     else
%         answer=inputdlg({'Min Area','Max Ecc'},'',1,{num2str(Area_crit),num2str(Ecc_crit)});
%         Area_crit=str2num(answer{1});
%         Ecc_crit=str2num(answer{2});
%         
%     end
    
    
end
L=bwlabel(grid_image_cleaned);



[r,c]=size(L);
while ~found_pts
    
    imagesc(L>0);
    title('click on object to eliminate; right-click to accept');
    [x,y,b]=ginput(1);
    
    killval=L(round(y),round(x));
    L(L==killval)=0;
    if b==3 | x< 1 | x > c | y<1 | y>r
        found_pts=1;
        stats=regionprops(L>0,'Centroid','Area','Eccentricity');
        
    end
    
end
%%

%now find the centroids in order
close all;
tocheck=1:9;ninrow=3;
sorted=[];
while ~isempty(tocheck)
    %sort them by row order
    c=[];
     for k=1:length(tocheck)
         c(k,1)=tocheck(k);
         c(k,2)=stats(tocheck(k)).Centroid(1);
         c(k,3)=stats(tocheck(k)).Centroid(2);
     end
    [~,idx1]=sort(c(:,3));
    
    
    rowcontents_index=idx1(1:ninrow);
    rowc=c(rowcontents_index,:);
    [~,idx2]=sort(rowc(:,2));
    sorted=[sorted,rowc(idx2(1),1)];
    origval=sorted(end);
    tocheck(tocheck==origval)=[];
    ninrow=ninrow-1;
    if ninrow==0
        ninrow=3;
    end
end
stats=stats(sorted);
%%
%figure;imagesc(grid_image);hold on;
row=1;col=1;
for k=1:length(stats)
    ninespots{row,col}=stats(k).Centroid;
 %   plot(stats(k).Centroid(1),stats(k).Centroid(2),'ro');
    col=col+1;
    if col==4
        col=1;
        row=row+1;
        if row==4
            break
        end
    end
end
   