%  @Copyright Reserved 
%  For non-commercial use, please cite the following papers:
%  LQ Wan, K Ronaldson, M Park, G Taylor, Y Zhang, JM Gimble, and G Vunjak-Novakovic, 2011. Micropatterned Mammalian Cells Exhibit Phenotype-Specific Left-Right Asymmetry, Proceedings of National Academy of Sciences, 108(30), pp. 12295-300.
%  H Zhang, K Ronaldson-Bouchard, G Vunjak-Novakovic, and LQ Wan, A Micropatterning Assay For Measuring Cell Chirality.
%  For any questions, please contact Dr. Leo Wan at wanq@rpi.edu. Thank you!

close all;
clear all;

%%%%%%%%reading files
dirx='/User/Desktop/LAB/Data/1'; % change directory
files=dir([dirx '/*.tif']);


%%%%%%%%%%%%% generating gradient filter:
half_size=6;  % filter parameter related to filter size
sigma=3;      % related to image intensity standard deviation
for i=-half_size:1:half_size
    for j=-half_size:1:half_size
       hx(i+half_size+1,j+half_size+1)=2*i/sigma^2*exp(-(i^2+j^2)/sigma^2);
       hy(i+half_size+1,j+half_size+1)=2*j/sigma^2*exp(-(i^2+j^2)/sigma^2);
    end
end



%%%%%%%%%%%%%%% processing each image
for fno=1:1:length(files)
    
    %%%%%load image and ring center
    filename_image=[dirx '/' files(fno).name];
    filename_position=[filename_image(1:(findstr(filename_image, '.')-1)), '.mat'];       
    image=imread(filename_image);
    load(filename_position);

    
    % image processing
    %%%%%%% applying graidnet filter:
    img_grdx=double(filter2(hx,image))+10^-3;
    img_grdy=double(filter2(hy,image));
    grad_mag=sqrt(img_grdx.^2+img_grdy.^2);
    orient=atan(img_grdy./img_grdx);
    [sizex,sizey]=size(grad_mag);
    
    
    %%%%%%%% determine local orientation
    size_sub=32;
    theta=0:1:179;
    for i=1:1: floor(sizex/size_sub)
        for j=1:1:floor(sizey/size_sub)
            [tempi, tempj]=ndgrid((i-1)*size_sub+1:i*size_sub, (j-1)*size_sub+1:j*size_sub);
            temp_idx=sub2ind([sizex, sizey],tempi(:),tempj(:));
            temp_orient = orient(temp_idx);
            temp_grad = grad_mag(temp_idx);
            AW=sum(exp(2*cos(2*(ones(length(temp_orient),1)*theta./180*pi-temp_orient*ones(1,length(theta)))))/exp(2).*(temp_grad*ones(1,length(theta))),1);
            [c,index]=max(AW);
            if index>91
                angle(i,j)=index-180;
            else
                angle(i,j)=index-1;
            end
        end
    end
% Note: here angle (i,j) is the gradient direction, not cell alignment direction.
% These two directions are perpendicular to each other



%%%%%%%%%%%%%Generating plots for alignments in the iamge
%     figure, imagesc(image,[0 255]); colormap(gray);
%     hold on;
% % note: angle transformation from image to x-y coordinates. 
% % here x & y is the coordinates on the figures
     [y,x]=ndgrid(size_sub/2:size_sub:sizex-size_sub/2, size_sub/2:size_sub:sizey-size_sub/2); 
%     quiver(x,y, cos((angle)*pi/180),-sin((angle)*pi/180),0.5); axis tight  
%     set(gca,'ytick',0:size_sub:sizex);
%     set(gca,'xtick',0:size_sub:sizey);
%     grid on;
%     hold off;


%Change Radius below if image size changes 

%%%%%%% caculate the deviation from the circumferential direciton
    rout=400;  % region of interest - outer radius
    rin=80;  % region of interest - inner radius 
    cent_ring_x=round(position(1)+position(3)/2);
    cent_ring_y=round(position(2)+position(4)/2);
    vect_x=x(:)-cent_ring_x;
    vect_y=y(:)-cent_ring_y;
    length_r=sqrt(vect_x.^2+vect_y.^2);
    % vector from center of window to ring center
    angle_r=atan(vect_y./vect_x)./pi*180; 
    % angles between the circumferential direction and cell aligned direction
    angle_chrality=(-angle(:))-(90+angle_r);
    % angles in the regions of interest between circumferential direction 
    % and cell aligned direction
    angle_list=angle_chrality(find(length_r>rin & length_r<rout)); %note changes here
    % wrap the angles smaller than 90 degree
    index_90=find(angle_list<-90);
    angle_list(index_90)=angle_list(index_90)+180;
%%%%%%%%% linear histgram plot
%     figure,hist(angle_list, -80:10:80);
%%%%%%%% circular histogram plot
%     figure, rose(angle_list./180*pi, 36)



%%%%%%%%%%% orientation variation with radius
%     radius = rin:size_sub:rout;
%     for i=1:1:(size(radius,2)-1)
%        angle_avg_r(i)=circ_mean(angle(find(length_r>radius(i)&length_r<radius(i+1)))./180*2*pi);
%     end
%     figure,plot(radius(2:end)-size_sub/2, angle_avg_r/2*180/pi);
    
%%%%%%%%% statistics about the cell alignment
    [p(fno) mu(fno) ul(fno) ll(fno)] = circ_mtest(double(angle_list/180*pi*2),0); % note using the circular statistics
    [s(fno) s0(fno)] = circ_std(double(angle_list/180*pi*2));
    %image_name(fno) = files(fno).name;
    fno          
end



%%%%%saving data in the matlab mat files
file_output=[dirx '/excelsummary.mat'];
save(file_output,'files','p','mu','ul','ll','s','s0');



%%%%%%%%%%%%%%%%exporting data into .txt file for veiw with microsoft excel
file_order=[dirx '/datatoexcel.txt'];
fid2=fopen(file_order,'a+');
fprintf(fid2,'%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t\n', 'File Name','nonzero?', 'Mean Anlge','Upper Limit','Lower Limit','standard deviation','standard deviation (log)','Chirality');
for i=1:1:length(files)
    fprintf(fid2,'%s\t %d\t %f6\t %f6\t %f6\t %f6\t %f6\t %f\t \n',files(i).name, p(i), mu(i)/2, ul(i)/2, ll(i)/2, s(i)/2, s0(i)/2, sign(p(i)*mu(i)));
end
% summary on the number of chirality in all images
numb_total=length(files);
numb_chral=sum(p);
non_chral=numb_total-numb_chral;
numb_anti=sum(p(find(p.*mu>0))); %+ represents anticlockwise (ie, couter clockwise)
numb_clk=numb_chral-numb_anti; % - represents clockwise
fprintf(fid2,'total:\t %d\t clockwise:\t %d\t nonchiral:\t %d\t anticlockwise:\t %d\t \n', numb_total, numb_clk, non_chral, numb_anti); 
fclose(fid2);

