%code to correct for microscope stage shift
%dedrifts image based on bead displacments assuming you are using island or
%edge assay

%Steps for stage dedrift are as follows:

%1.initially prints color figure of original displacements (figure 1) and
   %a binary image of the displacements (figure 2)

%2. Manually draw around area you want to use to dedrift image on figure 2. Note this
    %will be an area on the figure free of cells with 0 displacement. Also
    %the x and y corridnates for upper and lower limit of croppped area will
    %be used as an input.
    
%3. Cropped area will be printed for verification.

%4. Raw data file of dedrifted displacements will be saved and new displacements
    %will be printed on the screen. 

%5. Program will then loop through all files until File.NFiles is reached
%
%by Robert Steward 7/2013
function File = dedrift(Settings,File)

for k = 1:File.NFiles.Fluorescence
    
    %{
    ButtonName1=questdlg(strcat('Do you want to dedrift image ',num2str(k),' ?'), ...
                       'subtract stage drift','Yes');

   switch ButtonName1,
     case 'Yes',
    %}
    
    %load original displacement file
     tempFile = [File.pathname 'Displacements','/', 'GelDisp_', num2str(k), '.dat'];
     TempData=load(tempFile);
     
     %% 
        %reshape original displacment into square matrix and plot     
        x = reshape(TempData(:,1),sqrt(size(TempData(:,1),1)),sqrt(size(TempData(:,1),1)));
        y = reshape(TempData(:,2),sqrt(size(TempData(:,2),1)),sqrt(size(TempData(:,2),1)));
        DispX = reshape(TempData(:,3),sqrt(size(TempData(:,3),1)),sqrt(size(TempData(:,3),1)));%1e6 and (*Settings.PixelSize) removed for both DISP 
        DispY = reshape(TempData(:,4),sqrt(size(TempData(:,4),1)),sqrt(size(TempData(:,4),1)));
        Disp = sqrt(DispX.^2+DispY.^2);
        
        
       figh2=figure(1); set(figh2,'name','Original Displacements')
        h=pcolor(x,y,Disp); 
        set(gca,'YDir','reverse');colormap jet;
        view(2); colormap jet; shading flat; axis image; axis off; set(gca,'CLim',[0 Settings.disp_max]); colorbar; hold on;
        
        
        %loop below makes sure first image is dedrifted and susbequent
        %images are dedrifted using the same area. This is fine for
        %micropatterns, but may not be able to be used for every image in
        %edge assays.
        
        if (k==1)
        %manually draw roi that will be cropped
        figure(2)
        [BW, xrub, yrub] = roipoly(Disp);
        
      
     %crop area with no cells, make into matrix, and plot matrix
    
     Settings.RowCropTop = ceil (yrub(1,1));%top limit of row (y upper limit)
     Settings.RowCropBot = ceil (yrub(3,1));%row crop bottom  (y bottom limit)
     Settings.ColCropLeft = ceil (xrub(1,1));%begin col crop  (x upper limit)
     Settings.ColCropRight = ceil (xrub(3,1));%end col crop   (x bottom limit)
     DispXCrop = DispX(Settings.RowCropTop:Settings.RowCropBot,Settings.ColCropLeft:Settings.ColCropRight);
     DispYCrop = DispY(Settings.RowCropTop:Settings.RowCropBot,Settings.ColCropLeft:Settings.ColCropRight);
     xCrop = x(Settings.RowCropTop:Settings.RowCropBot,Settings.ColCropLeft:Settings.ColCropRight);
     yCrop = y(Settings.RowCropTop:Settings.RowCropBot,Settings.ColCropLeft:Settings.ColCropRight);
     DispCrop = sqrt(DispXCrop.^2+DispYCrop.^2);
     
     %view cropped area to verify there are no cells  
       figh2=figure(3); set(figh2,'name','Cropped Displacements')
        h=pcolor(xCrop,yCrop,DispCrop); 
        set(gca,'YDir','reverse');
        view(2); colormap jet; shading flat; axis image; axis off; set(gca,'CLim',[0 Settings.disp_max]); colorbar; hold on;
      
        
   %average x and y displacement for subtraction
   
   %DispCrop = mat2vertcol(DispCrop);
   DispXCrop = mat2vertcol(DispXCrop);
   DispYCrop = mat2vertcol(DispYCrop);
   xAvg = mean(DispXCrop);
   yAvg = mean(DispYCrop);
   %DispAvg = mean(DispCrop);
        
   DispX2 = (mat2vertcol(DispX)) - xAvg;
   DispY2 = (mat2vertcol(DispY)) - yAvg;
   
         %save dedrifted displacements first then  plot new displacements   
         f3 = fopen([File.pathname 'Displacements','/', 'DedriftGelDisp_', num2str(k), '.dat'],'w');
         for m = 1:size(TempData,1),
             fprintf(f3,'%15.5e%15.5e%15.5e%15.5e\n', x(m), y(m), DispX2(m), DispY2(m));
         end;
         fclose(f3); 
   
   
        %reshape new displacements into matrix for plotting
        %x = reshape(TempData(:,1),sqrt(size(TempData(:,1),1)),sqrt(size(TempData(:,1),1)));
        %y = reshape(TempData(:,2),sqrt(size(TempData(:,2),1)),sqrt(size(TempData(:,2),1)));
        DispX2 = reshape(DispX2,sqrt(size(TempData(:,3),1)),sqrt(size(TempData(:,3),1)));%1e6 and (*Settings.PixelSize) removed for both DISP 
        DispY2 = reshape(DispY2,sqrt(size(TempData(:,4),1)),sqrt(size(TempData(:,4),1)));
        Disp2 =  sqrt(DispX2.^2+DispY2.^2);
   
  
       %plot 
       figh2=figure(4); set(figh2,'name','New Displacements')
        h=pcolor(x,y,Disp2); 
        set(gca,'YDir','reverse');
        view(2); colormap jet; shading flat; axis image; axis off; set(gca,'CLim',[0 Settings.disp_max]); colorbar; hold on;
       
        pause 
        
        else
 %%    
    %crop area with no cells, make into matrix, and plot matrix
    
    clearvars -except xrub yrub DispX DispY x y File Settings k TempData tempFile
    
 
    
    
     Settings.RowCropTop = ceil (yrub(1,1));%top limit of row (y upper limit)
     Settings.RowCropBot = ceil (yrub(3,1));%row crop bottom  (y bottom limit)
     Settings.ColCropLeft = ceil (xrub(1,1));%begin col crop  (x upper limit)
     Settings.ColCropRight = ceil (xrub(3,1));%end col crop   (x bottom limit)
     DispXCrop = DispX(Settings.RowCropTop:Settings.RowCropBot,Settings.ColCropLeft:Settings.ColCropRight);
     DispYCrop = DispY(Settings.RowCropTop:Settings.RowCropBot,Settings.ColCropLeft:Settings.ColCropRight);
     xCrop = x(Settings.RowCropTop:Settings.RowCropBot,Settings.ColCropLeft:Settings.ColCropRight);
     yCrop = y(Settings.RowCropTop:Settings.RowCropBot,Settings.ColCropLeft:Settings.ColCropRight);
     DispCrop = sqrt(DispXCrop.^2+DispYCrop.^2);
   
   %average
   DispXCrop = mat2vertcol(DispXCrop);
   DispYCrop = mat2vertcol(DispYCrop);
   xAvg = mean(DispXCrop);
   yAvg = mean(DispYCrop);
   
        
   DispX2 = (mat2vertcol(DispX)) - xAvg;
   DispY2 = (mat2vertcol(DispY)) - yAvg;
   
         %save dedrifted displacements first then  plot new displacements   
         f3 = fopen([File.pathname 'Displacements','/', 'DedriftGelDisp_', num2str(k), '.dat'],'w');
         for m = 1:size(TempData,1),
             fprintf(f3,'%15.5e%15.5e%15.5e%15.5e\n', x(m), y(m), DispX2(m), DispY2(m));
         end;
         fclose(f3); 
    
        end    
         
   end
   
   close all
         
end