% TRACTION uses Fourier Transform Traction Microscopy to compute the
% tractions that a cell exerts on its substrate. The algorithm is described
% in Butler al, AJP-Cell, 2002. A typical experiment involves one or more
% images of the fluorescent beads embeded in the gels, the corresponding 
% one of more phase contrast images of the cell, and a reference image obtained 
% after trypsinizing the cell.
%
% Variables:
%   SETTINGS contains the user-defined information needed by the algoritm
%   FILE contains information about the each of the files opened
%   RESULTS contains the results of the computation.
%
% QUICK GUIDE
% 1) The program promts the user for the first fluorescent image of the
% sequence, the first phase contrast image of the sequence and the
% reference image.
%
% 2) The user is asked to draw the cell boundary. This boundary is used to
% define the ROI over which the calculation is done but the contour of the
% cell is only used to compute constrained tractions not unconstrained
% tractions
%
% 3) The region of interest is output in the folder 'croppeddata'. The displacement
% fields are output as in the folder 'displacements', the traction fields are output 
% in the folder 'tractions'. 

% by Xavier Trepat, Iva Tolic-Norrelykke, and James Butler, 2000-2008.
% xtrepat@hsph.harvard.edu or xtrepat@gmail.com
% jbutler@hsph.harvard.edu
%modified slightly by Robert Steward plot intercellular stresses as well
% note a few additions:
%pixelation of traction images can be reduced by adjusting Filter.CutFreq


%clear all;
%close all;
%%
%Set variables
obj = 5;
mag = 1; % don't really need for our calculation. our microscope magnification in 1x, some microscopes has different magnification themselves.
Settings.PixelSize = (6.5/mag)/obj; % size of one pixel after magnification in MICRONS (typically physical size of the pixel on the CCD * Magnification)
Settings.Imgfmt='tif'; % image format
Settings.Young = 1.2*1e3; % Young's modulus of the gel. NOT shear modulus! changed from e9 to get proper scale measurements
Settings.Poisson = 0.5; % Poisson ratio
Settings.Filter.Type = 'LP'; % 'LP' for a low pass filter or 'wiener' for an adaptive wiener filter.

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Open the images and store the names in structure FILE. Create folders to
% store results. For a stack of N measurements, image names must be numbered from 1 to N
%only 2 trypsin images needes. first tyrpsin image "trypsin1" and digit at end of last trypsin image must be total number of files. ex. - for a
%stack of 10 images name trypsin file "trypsin10.tif". 
%trypsin file needed - RS

%File = OpenFiles(Settings);

%% Crop images

%Image parameter for cropping
%I = y, J = x
% ROI list x and y limit for cropping
% See excel spreadsheet for calculating cropping limits

% run image_shape code first to find the exact ROI /miku/

% then rename the images

ROI(1).i1=150;
ROI(1).i2=1749; 
ROI(1).j1=150;
ROI(1).j2=1749;

 %j is width. i is height
%% This ROI will effect the block size in disp on block size, there we need to adjust size of this image if not works

%File = cell_cropper_Practice(Settings,File,ROI);

%%
% If I need any new cropping other than specified ROI, than use this
 % specifiy [xmin ymin width height] of cropping
xmin = 475; ymin = 475; % specify the xmin and ymin with crop_size
crop_size = 1048;

%File = New_Cropping(Settings,File,xmin, ymin, crop_size, imagenum);

%%
%Compute Gel Displacement 

%tic

Settings.resolution = 32; % Size of each window of the PIV analysis
Settings.overlap = 0.75; % Overlap between windows for PIV analysis. 0 means no overlap, 1 means full overlap (and of course, crash of the program)

%File = displacement_finder(Settings,File);

%toc
%%
data_Set= 8;
imnum = 100;
Settings.resolution = 32; % Size of each window of the PIV analysis
Settings.overlap = 0.75; % Overlap between windows for PIV analysis. 0 means no overlap, 1 means full overlap (and of course, crash of the program)

%File = displacement_finder_batch(Settings,File,data_Set,imnum);

%%
%Draw boundary around ROI
%Note only need if doing constrained tractions

%File = draw_boundary(Settings,File);
     
%%
%compute stage drift

Settings.disp_max = 4;% sets color bar limit for displacement 

%File = dedrift(Settings,File);%use for bead defrift
%%

%File = dedrift_batch(Settings,File,imnum,data_Set);

%%
%subtract rotation from island(code from Dhananjay)

Settings.DedriftImage = 'yes';% was image dedrifted

%File = subtract_rotation(Settings,File);%subtract rotation from
%displacements using Dhananjay's code

%File = rotation_translation(Settings,File);%subtract rotation from
%displacements using Jae's code


%%
% Compute the traction field
%displacements come out slightly different after this routine and overwrite
%original displacements calculated from displacement_finder. - this part is
%commented out
% results for all files written into Results.dat - rms traction (Pa), netmoment (pj),
% and strain energy (pj)
%code modified to crop out extra expanded matrix data
%I have code that does this and have cut and paste code from chan that does
%the same thing
%if doing edge or interior image matrix must be flipped on all optical edges to
%reduce noise at boundary
%modified for island analysis and for dedrift or disp with rotation
%subtracted - RS 7/2013
Settings.Filter.CutFreq = 1; % cutoff frequency in units of the nyquist frequency originally %original frequency 1 (changed to 0.65), use 0 for rama code
Settings.FlipMatrix = 'no';%flip matrix for edge or interior assay only!
Settings.GelHeight = 100;% Gel Height input if using traction_finite
Settings.CalcRot = 'no';%calculate monolayer rotation using Jae's (traction) code (for island only), not need using disp rot subtraction using dhananjay or jae's code
Settings.height = 'finite';%'infinite','finite','constrained' , tells what traction code to use
Settings.DedriftImage = 'yes';%was image dedrifted above?, set to no if image not dedrtifted or image dedrifted and rotation subtracted
Settings.SubRotation = 'no';% set to 'yes' if rotation was subtracted using jae's or DT's disp subtraction code
%using DT code and 'no' if not, variable no longer used

Settings.TracBoundCond = 'unconstrained';%traction boundary condition constrained or unconstrained

%constrained tractions orientation calculaiton offset by 41, F.Y.I.

%tic

%File = traction_finder(Settings,File,ROI);

%toc
%%

%File = traction_finder_batch(Settings,File,imnum,data_Set);

%%
% Plot the traction results

Settings.reshapeMatrix = 'yes';%reshape image to reduce traction pixelation in plot_traction.m (if thats even a word!)
Settings.PlotDisp = 'yes';%plot displacements yes or no
Settings.PlotTracX = 'yes';%plot tractions in x
Settings.PlotTracY = 'yes';%plot tractions in y
Settings.PlotRmsTrac = 'yes';%plot rms tractions
Settings.cropFig = 'no';

%Traction and Stress colorbar limits for plots
Settings.disp_max = 4;
Settings.trac_max = 60;
Settings.tracx_max = 30;
Settings.tracx_min = 0;

%plot_traction(Settings,File,ROI);
%%

%plot_traction_batch(Settings,File,imnum,data_Set);

%%

%%plot displacements as vectors i.e. velocity
imagenum = 100;
start_time = 5;
interval_period =5;
ending_time = 500;
%plot_velocity(Settings,File,imagenum,start_time,interval_period,ending_time);
%%plot velocity code is changed here, check again and change in lab
%plot_str_moms(start_time,interval_period,ending_time);
%plot_scatter(imagenum,obj); %3d plot of entire data set 

%%
%plot_str_moms_batch(start_time,interval_period,ending_time,data_Set);

%%
%Crop any size and get corresponsind traction and displacements values
in_ROI = 18; %define initial ROI we used for original big inmage
crop_size = 384; %384 pixel equals to ~500um

%section_cropper(Settings,File,in_ROI,imagenum,crop_size,start_time,interval_period,ending_time);

%%
%Make domain file for monolayer stress microscopy
%if domain drawn using imagej both variables must be set to yes
%if domain previously drawn for any other reason only set DomDrawn to yes
%if none of the above done set both variables to no
Settings.DomDrawn = 'no';%was domain for monolayer boundary drawn?
Settings.DomImagej = 'no';%was domain drawn using imagej?

%mark_circular_domain(File,Settings,ROI);

%this practice code takes only one data set point and copies that for all,

%mark_circular_domain_practice(File,Settings,ROI);

%%
%mark_circular_domain_practice_batch(File,Settings,imnum,data_Set);

%%
%Function to run the stress code
% Md. Mydul Islam Miku 2016

%tic
%assign number of images to analyze

%Run_StressCode(imagenum); %running code

%toc
%%
%Run_StressCode_batch(imnum,data_Set);

%%
% check if FEM/tractions/time is empty or not. if empty, you have to copy it from tractions/time to here
%otherwise code will not run. /.miku./ 

%Plot Intercellular Stresses
Settings.reshapeMatrix = 'yes';% set to yes to reduce pixelation
Settings.max_shear_max = 300;
Settings.avg_norm_min = 0;
Settings.avg_norm_max = 300;
Settings.domSize = [196 196];%matrix domain size for stress plots. domsize = CROPPED image matrix size/resolution.
%example if cropped image size is 1152 x 1152 (must be square) with 32 bit
%resolution domSize = 1152/32 = 36
%value for nXPts can be used as well

%plot FEM might not show figures is pixel to micron in model file doesn't
%match, so check if it changed according to objective i.e 5x, 10x
%set both variables to yes to plot intercellular stresses 
%Only set "EvenStresses" variable to yes if you want to calculate stresses,
%but not plot stresses
Settings.PlotStresses = 'yes';%plot stress distribution
Settings.EvenStresses = 'yes';%even stress matrix to match traction matrix

%tic        
%plot_FEM_results(File,Settings);
%toc
%%

%plot_FEM_results_batch(File,Settings,imnum,data_Set);

%%
%raw_data_analysis(imagenum,start_time,interval_period,ending_time);
% make sure time point is consistent with number of images
% miku january 9,2017
%%
%don't need raw_data_analysis anymore
%raw_data_analysis_batch(start_time,interval_period,ending_time,imnum,data_Set);

%Collect_FEM_data_batch(data_Set);

%%
%go to make_movie and make movies from results images
%make_movie(imagenum);

%%
%max_min_stress_plot(imagenum,Settings);
%%
%dont forget to see these
% manipulate movie making
%disp_fixing
%surface_plot_practice
%crop_disp_surf
%image_merging_practice
%phase_velocity_main_routine %% to get phase velocities from crop_pcs 
% prac %%to plot phase velocities
% phase_video_maker