function File = traction_finder(Settings,File,ROI,Results);
% TRACTION_FINDER computes the traction fields.
% Xavier Trepat 2008
%Robert Steward 2011

%change directory
cd([File.pathname,'Displacements']);

mkdir([File.pathname,'Results']);


for k = 1:File.NFiles.Fluorescence,   %par for parallel processing 

    
    disp(['Computing Tractions of ',File.Name(k).CropFluorescence]);
    
    mkdir([File.pathname,'Tractions', '/','Time',num2str(k)]);
    

    %if exist(File.Name(k).CropFluorescence),  
        
        % get the data and reshape into matrix

        if (strcmpi(Settings.DedriftImage,'Yes')) == 1

            tempFile = [File.pathname '/Displacements','/', 'DedriftGelDisp_', num2str(k), '.dat'];

       elseif (strcmpi(Settings.SubRotation,'Yes')) == 1 
        
            tempFile = [File.pathname '/Displacements','/', 'RotSubGelDisp_', num2str(k), '.dat'];
            
        else
            
            tempFile = [File.pathname '/Displacements','/', 'GelDisp_', num2str(k), '.dat'];
     
        end

        TempData=load(tempFile);
        x=reshape(TempData(:,1),sqrt(size(TempData,1)),sqrt(size(TempData,1)));
        y=reshape(TempData(:,2),sqrt(size(TempData,1)),sqrt(size(TempData,1)));
        DispX=reshape(TempData(:,3),sqrt(size(TempData,1)),sqrt(size(TempData,1)));
        DispY=reshape(TempData(:,4),sqrt(size(TempData,1)),sqrt(size(TempData,1)));
        
        
        
    %%     
            if (strcmpi(Settings.FlipMatrix,'Yes')) == 1
                                                            
                
                
                %from chan
                % expand domain size by flipping data. in case of monolayer, flipping helps to reduce noise on the boundaries
                % flip is below
                npts = sqrt(size(TempData,1));
                temp123 = [flipud(x(2:floor(npts/2),:)); x; flipud(x(floor(npts/2)+1:npts-1,:))];
                x = [fliplr(temp123(:,2:floor(npts/2))) temp123 fliplr(temp123(:,floor(npts/2)+1:npts-1))];
                
                temp123 = [flipud(y(2:floor(npts/2),:)); y; flipud(y(floor(npts/2)+1:npts-1,:))];
                y = [fliplr(temp123(:,2:floor(npts/2))) temp123 fliplr(temp123(:,floor(npts/2)+1:npts-1))];
                
                temp123 = [flipud(DispX(2:floor(npts/2),:)); DispX; flipud(DispX(floor(npts/2)+1:npts-1,:))];
                DispX = [fliplr(temp123(:,2:floor(npts/2))) temp123 fliplr(temp123(:,floor(npts/2)+1:npts-1))];
                
                temp123 = [flipud(DispY(2:floor(npts/2),:)); DispY; flipud(DispY(floor(npts/2)+1:npts-1,:))];
                DispY = [fliplr(temp123(:,2:floor(npts/2))) temp123 fliplr(temp123(:,floor(npts/2)+1:npts-1))];
                % flip is above
                
                % Even dimensions
                if(mod(size(DispX,1),2)), DispX=DispX(1:end-1,:); DispY=DispY(1:end-1,:); end;
                if(mod(size(DispX,2),2)), DispX=DispX(:,1:end-1); DispY=DispY(:,1:end-1); end;
                if(mod(size(x,1),2)), x=x(1:end-1,:); y=y(1:end-1,:); end;
                if(mod(size(x,2),2)), x=x(:,1:end-1); y=y(:,1:end-1); end;   
                 %}
                
            else
        
               % Even dimensions
               if(mod(size(DispX,1),2)), DispX=DispX(1:end-1,:); DispY=DispY(1:end-1,:); end;
               if(mod(size(DispX,2),2)), DispX=DispX(:,1:end-1); DispY=DispY(:,1:end-1); end;
               if(mod(size(x,1),2)), x=x(1:end-1,:); y=y(1:end-1,:); end;
               if(mod(size(x,2),2)), x=x(:,1:end-1); y=y(:,1:end-1); end;       
     
            end
     
            
            
           % filter displacments
            if(strcmpi(Settings.Filter.Type,'LP') == 1)%='LP'
            %if(Settings.Filter.Type == 'LP')
                DispX = xfilter(DispX, Settings.Filter.CutFreq);
                DispY = xfilter(DispY, Settings.Filter.CutFreq);
            
            else
                
                DispX = wiener2(DispX, [3 3]);
                DispY = wiener2(DispY, [3 3]);
            end
        
    
        %%
        %decide whether to run finite gel thickness unconstrained
        %tractions, infinite thickness unconstrained tractions, or
        %constrained tractions
        %compare strings using strcmp, also strcmpi can be used for case
        %insensitive
        
      if (strcmpi(Settings.height,'infinite')) == 1%compute infinite gel thickness unconstrained
          
        [netMoment, strainEnergy, rmsTraction, orientation, tx, ty] = traction_boussinesq(  x*Settings.PixelSize*1e-6, ...
                                                                                    y*Settings.PixelSize*1e-6, ...
                                                                                    DispX*Settings.PixelSize*1e-6, ...
                                                                                    DispY*Settings.PixelSize*1e-6, ...
                                                                                    Settings.Poisson, ...
                                                                                    Settings.Young);                                                                           
                                                                                      
                                                                                
      end
                                                                                
      if (strcmpi(Settings.height,'finite')) == 1 %compute finite gel thickness unconstrained  
          
            
        [netMoment, netMomentNew, strainEnergy, rmsTraction, tzx, tzy] = traction_finite_ij(  x*Settings.PixelSize*1e-6, ...
                                                                                    y*Settings.PixelSize*1e-6, ...
                                                                                    DispX*Settings.PixelSize*1e-6, ...
                                                                                    DispY*Settings.PixelSize*1e-6, ...
                                                                                    Settings.GelHeight*1e-6,...
                                                                                    Settings.Poisson, ...
                                                                                    Settings.Young); 
                                                                                
      end
     
       if (strcmpi(Settings.height,'constrained')) == 1 %compute constrained tractions 
           
          cellboundary = load([File.pathname 'croppeddata','/','cellboundary', num2str(k), '.mat']);
        
         % The following form of constrained traction computation is still subjected to debugging
         %x,y,DispX,DispY,xrub,yrub all originally multiplied by *Settings.PixelSize*1e-6
         %Settings.PixelSize * 1e-6
         
             
         
        [rmsTraction, netMoment, strainEnergy, orientation] = traction_boussinesq_constrained2(x, ...
                                                                                                y, ...
                                                                                                DispX, ...
                                                                                                DispY, ...
                                                                                                Settings.Poisson, ...
                                                                                                Settings.Young * 1e-12, ...
                                                                                                cellboundary.xrub, ...
                                                                                                cellboundary.yrub, ...
                                                                                                Settings.PixelSize,...
                                                                                                Settings.CalcRot,...
                                                                                                File.pathname,...
                                                                                                k);
                                                                                                
                                                                                                                                                                                
       end 
                                                                                                
                                                                                                                                                                                                                                          
             %flip matrix for interior or edge assay                                                                   
            if (strcmpi(Settings.FlipMatrix,'Yes')) == 1  
             
                DispX = DispX(floor(npts/2):floor(npts/2)+npts-1,floor(npts/2):floor(npts/2)+npts-1);
                DispY = DispY(floor(npts/2):floor(npts/2)+npts-1,floor(npts/2):floor(npts/2)+npts-1); 
                x = x(floor(npts/2):floor(npts/2)+npts-1,floor(npts/2):floor(npts/2)+npts-1);
                y = y(floor(npts/2):floor(npts/2)+npts-1,floor(npts/2):floor(npts/2)+npts-1);
                   
                 if (strcmpi(Settings.height,'infinite')) == 1    
                    %take gels from infinite thickness 
                    tx = tx(floor(npts/2):floor(npts/2)+npts-1,floor(npts/2):floor(npts/2)+npts-1);
                    ty = ty(floor(npts/2):floor(npts/2)+npts-1,floor(npts/2):floor(npts/2)+npts-1);
                 else
                    %take tractions from finite thickness 
                    tx = tzx(floor(npts/2):floor(npts/2)+npts-1,floor(npts/2):floor(npts/2)+npts-1);
                    ty = tzy(floor(npts/2):floor(npts/2)+npts-1,floor(npts/2):floor(npts/2)+npts-1);
                
                 end
             % flip is above    
                %}                                                   
             end                                                                       
                                                                                                                                          
         % prepare for storage
         
        if (strcmpi(Settings.height,'infinite')) == 1
            
          
                
             tract = rmsTraction;
             netMom = netMoment;
             strnEn = strainEnergy;
             angle = mat2vertcol(orientation); 
             DispX  	= 	mat2vertcol(DispX);
             DispY  	= 	mat2vertcol(DispY);
             tx 	= 	mat2vertcol(tx);
             ty 	= 	mat2vertcol(ty);
             x 	= 	mat2vertcol(x);
             y 	= 	mat2vertcol(y); 
                
         
             
        end
             
             
        if (strcmpi(Settings.height,'finite')) == 1  
              
             tract = rmsTraction;
             netMom = netMoment;
             netMomNew = netMomentNew;
             strnEn = strainEnergy;
             DispX  	= 	mat2vertcol(DispX);
             DispY  	= 	mat2vertcol(DispY);
             tx 	= 	mat2vertcol(tzx);
             ty 	= 	mat2vertcol(tzy);
             x 	= 	mat2vertcol(x);
             y 	= 	mat2vertcol(y); 
             
        end
         
         
          
         if (strcmpi(Settings.height,'constrained')) == 1  
              
             tract = rmsTraction;
             netMom = netMoment;
             angle = orientation;
             strnEn = strainEnergy;
             %DispX  	= 	mat2vertcol(DispX);
             %DispY  	= 	mat2vertcol(DispY);
             %tx 	= 	mat2vertcol(ctx);
             %ty 	= 	mat2vertcol(cty);
             %x 	= 	mat2vertcol(cx);
             %y 	= 	mat2vertcol(cy); 
             
         end
             
       
    
       
       
     
%==========================================================================
       %%store infinite unconstrained or finite unconstrained traction only!
       if (strcmpi(Settings.height,'infinite')) == 1
         % store tractions
         f3 = fopen([File.pathname 'Tractions','/', 'Time',num2str(k),'/', 'traction','.dat'],'w');
         for m = 1:size(x,1),
             fprintf(f3,'%15.5e%15.5e%15.5e%15.5e\n', x(m), y(m), tx(m), ty(m));   
         end;
         fclose(f3);   
         
       end
       
       if (strcmpi(Settings.height,'finite')) == 1
         % store tractions
         f3 = fopen([File.pathname 'Tractions','/', 'Time',num2str(k),'/', 'traction','.dat'],'w');
         for m = 1:size(x,1),
             fprintf(f3,'%15.5e%15.5e%15.5e%15.5e\n', x(m), y(m), tx(m), ty(m));   
         end;
         fclose(f3);   
         
       end
%==========================================================================         
      % store rmstraction, contractile moment, and strain energy    
       if (strcmpi(Settings.height,'infinite')) == 1
               
         
         results = [tract netMom strnEn angle];
         
       end

        if (strcmpi(Settings.height,'finite')) == 1

         results = [tract netMom netMomNew strnEn];
         
        end
        
        if (strcmpi(Settings.height,'constrained')) == 1
               
         
         results = [tract netMom strnEn angle];
         
       end
         
         if k == 1,
             
             finalResults = results;
             
         else
             
             finalResults = [finalResults;results]
        
         end
         
        if (k == File.NFiles.Fluorescence)
         
          
           if (strcmpi(Settings.height,'infinite')) == 1   
                
                f3 = ([File.pathname '/Results','/', 'TractionResults_infinite', '.dat']);
           end
           
        if (strcmpi(Settings.height,'finite')) == 1
                 
                 f3 = ([File.pathname '/Results','/', 'TractionResults_finite', '.dat']);
                
        end
         
        if (strcmpi(Settings.height,'constrained')) == 1   
                
                f3 = ([File.pathname '/Results','/', 'TractionResults_constrained', '.dat']);
           end
             
         
          dlmwrite(f3,finalResults,'precision',6);
         
        end
         %}
          
                                                                                                                                                       
                                                                                
         %end                                                                              
              
                                        
end
    
    