% This Matlab script quantifies polarotactic responses of tethered flies in
% resonse to a constantly rotating evector. All .TXT files within the Matlab
% folder will be processed. The .TXT files are the output
% files of the tracking process using FIJI (macro with instructions can be downloaded at 
% www.flygen.org/skylight-navigation). Please make sure to define 
% the correct data column (Angle1) to be extracted from the .TXT files.
% For each processed file this script will automatically save a .CSV with 
% processed parameters as well as an image for each file containing plots 
% visializing these polarotactic parameters for better overview. 
% Some individual values (eg. fly number, sex,time of experiment) are 
% extracted from the respective filenames. Those filenames were manually
% entered for each fly when recording the videos.



sad = transpose( dir( '*.txt' ) );

    for file = sad
        disp( file.name )
        filename = (file.name)

delimiter = '\t';

%% Read columns of data as strings:
% For more information, see the TEXTSCAN documentation.
formatSpec = '%*s%s%*s%*s%*s%s%[^\n\r]';

%% Open the text file.
fileID = fopen(filename,'r');

%% Read columns of data according to format string.
% This call is based on the structure of the file used to generate this
% code. If an error occurs for a different file, try regenerating the code
% from the Import Tool.
dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter,  'ReturnOnError', false);

%% Close the text file.
fclose(fileID);

%% Convert the contents of columns containing numeric strings to numbers.
% Replace non-numeric strings with NaN.
raw = repmat({''},length(dataArray{1}),length(dataArray)-1);
for col=1:length(dataArray)-1
    raw(1:length(dataArray{col}),col) = dataArray{col};
end
numericData = NaN(size(dataArray{1},1),size(dataArray,2));

% Converts strings in the input cell array to numbers. Replaced non-numeric
% strings with NaN.
rawData = dataArray{2};
for row=1:size(rawData, 1);
    % Create a regular expression to detect and remove non-numeric prefixes and
    % suffixes.
    regexstr = '(?<prefix>.*?)(?<numbers>([-]*(\d+[\,]*)+[\.]{0,1}\d*[eEdD]{0,1}[-+]*\d*[i]{0,1})|([-]*(\d+[\,]*)*[\.]{1,1}\d+[eEdD]{0,1}[-+]*\d*[i]{0,1}))(?<suffix>.*)';
    try
        result = regexp(rawData{row}, regexstr, 'names');
        numbers = result.numbers;
        
        % Detected commas in non-thousand locations.
        invalidThousandsSeparator = false;
        if any(numbers==',');
            thousandsRegExp = '^\d+?(\,\d{3})*\.{0,1}\d*$';
            if isempty(regexp(numbers, thousandsRegExp, 'once'));
                numbers = NaN;
                invalidThousandsSeparator = true;
            end
        end
        % Convert numeric strings to numbers.
        if ~invalidThousandsSeparator;
            numbers = textscan(strrep(numbers, ',', ''), '%f');
            numericData(row, 2) = numbers{1};
            raw{row, 2} = numbers{1};
        end
    catch me
    end
end


%% Split data into numeric and cell columns.
rawNumericColumns = raw(:, 2);
rawCellColumns = raw(:, 1);


%% Replace non-numeric cells with NaN
R = cellfun(@(x) ~isnumeric(x) && ~islogical(x),rawNumericColumns); % Find non-numeric cells
rawNumericColumns(R) = {NaN}; % Replace non-numeric cells

%% Allocate imported array to column variable names
Label1 = rawCellColumns(:, 1);
Angle1 = cell2mat(rawNumericColumns(:, 1));


%% Clear temporary variables
clearvars filename delimiter formatSpec fileID dataArray ans raw col numericData rawData row regexstr result numbers invalidThousandsSeparator thousandsRegExp me rawNumericColumns rawCellColumns R;



% Define range of angular data to be analyzed
Angle1a=Angle1(1:18100)

% Define e-vector positions over time within the angular reference frame of the 
% recording/tracking. The way we inserted the polarizer into the filter cassette, 
% rotation starts at 90° with approximately -0,1°/frame.
evectorasrow=linspace(180,0,1810)
evector2column=evectorasrow'
evectora=repmat(evector2column,11,1)
evector=evectora(906:19005)


% Angle between fly and evector 
Anglediff=Angle1a-evector
Anglediff(Anglediff<-90)=Anglediff(Anglediff<-90)+180
Anglediff(Anglediff>90)=Anglediff(Anglediff>90)-180

% Preferred Evector 
PreferredEvector2=Anglediff+180
PreferredEvector2(PreferredEvector2>180)=PreferredEvector2(PreferredEvector2>180)-180
FinalPreferredEvector=-PreferredEvector2+180

% E-vector speed (deg/frame)
evectornew=-0.099502487562191

% Calulate angular speed
speedarad=deg2rad(Angle1(1:18001))
speedaraddouble=speedarad*2
speedaradunwrap=unwrap(speedaraddouble)
speedaradunwrapdeg=rad2deg(speedaradunwrap)
speedafinal2=diff(speedaradunwrapdeg)
speedafinal3=speedafinal2/2
speedafinalmovmean=movmean(speedafinal3,20)
speedafinalmovmean2=movmean(speedafinalmovmean,20)
speedafinalmovmean3=movmean(speedafinalmovmean2,20)

speedafinaljumpmean=reshape(speedafinal3,600,[])
yy=sum(speedafinaljumpmean,1)./size(speedafinaljumpmean,1)
speedafinaljumpmeanfinal=reshape(repmat(yy, size(speedafinaljumpmean,1),1),1,[])

% Quantify polarotxis when the difference between angular speeds of the fly 
% and the e-vector is smaller than 0.05 deg/frame
speeddiffdiffdiff=evectornew-speedafinaljumpmeanfinal
numberframespolarotaxis=sum(speeddiffdiffdiff>-0.05 & speeddiffdiffdiff<0.05)
number10secondsofpolarotaxis=numberframespolarotaxis/600
polarotaxisonly=speeddiffdiffdiff>-0.05 & speeddiffdiffdiff<0.05

% Define polarotactic segments
Angle1b=Angle1a*2
evectorb=evector*2
normDeg=mod(Angle1b-evectorb,360)
absDiffDeg=min(360-normDeg, normDeg)
a=absDiffDeg/2

polarotaxisonlydiffeneceangles=a(polarotaxisonly)
polarotaxisonlydiffeneceanglesmean=mean(polarotaxisonlydiffeneceangles)
polarotaxisonly180=polarotaxisonly*180

a3=a(1:18000)
a3(polarotaxisonly<0.5)=NaN

speedafinaljumpmeanpolarotaxisonly2=reshape(a3,600,[])
speedafinaljumpmeanpolarotaxisonlymean2=mean(speedafinaljumpmeanpolarotaxisonly2)
FINALPreferredevectorovertime=transpose(speedafinaljumpmeanpolarotaxisonlymean2)

% Quantify preferred e-vector during polarotaxis only
arad=deg2rad(FinalPreferredEvector(polarotaxisonly))
adoublerad=arad*2
adoubleradpiminpi=shiftAnglesFromMinusPIToPI(adoublerad)
meanadoubleradpiminpi = AverageAnglesRadians(adoubleradpiminpi);
meanadoubleradpiminpideg = rad2deg(meanadoubleradpiminpi);
meanadoubleraddeg0360=meanadoubleradpiminpideg + (meanadoubleradpiminpideg < 0)*360

% Final Preferred E-vector
FinalMeanPreferredEvectorPolarotaxisonly=meanadoubleraddeg0360/2

% Final Preferred E-vector per 10s segment over time
a4=FinalPreferredEvector(1:18000)
a4=a4(1:18000)
a4(polarotaxisonly<0.5)=NaN
a4=reshape(a4,600,[])
yypolarotaxisonly2=sum(a4,1)./size(a4,1)
speedafinaljumpmeanfinalpolarotaxisonly2=reshape(repmat(yypolarotaxisonly2, size(a4,1),1),1,[])
a4rad=deg2rad(a4)
a4doublerad=a4rad*2
a4doubleradpiminpi=shiftAnglesFromMinusPIToPI(a4doublerad)
meana4doubleradpiminpi = AverageAnglesRadians(a4doubleradpiminpi);
meana4doubleradpiminpideg = rad2deg(meana4doubleradpiminpi);
meana4doubleraddeg0360=meana4doubleradpiminpideg + (meana4doubleradpiminpideg < 0)*360
a4speedafinaljumpmeanpolarotaxisonlymean2=meana4doubleraddeg0360/2
FINALPreferredevectorovertime=transpose(a4speedafinaljumpmeanpolarotaxisonlymean2)

% Modify plotting data of the fly's heading and e-vector for continuous
% plotting
Plotanglediff=abs(diff(Angle1a))
Plotangle=Angle1a
Plotangle(Plotanglediff>140)=NaN
Plotevectordiff=abs(diff(evector))
Plotevector=evector
Plotevector(Plotevectordiff>140)=NaN

% Create figures
figure
ax2=subplot(4,3,[1 3])
plot(Plotangle)
ax = gca
ax.XTick = [0 1810 3620 5430 7240 9050 10860 12670 14480 16290 18100]
ax.YTick = [0 45 90 135 180]
ax.FontSize = 5
ax.XAxis.Exponent = 0
ax.XAxis.TickLabelFormat = '%6.0f'
axis tight
ylabel('Heading');
ylim([0 180]);
xlim([0 18000])
grid('on');
hold on 
area(polarotaxisonly180,'FaceAlpha',0.1,'EdgeAlpha',0)
hold off
hold on
plot(Plotevector)


% Extract metadata from filename
filenamestring = regexprep(Label1{3,1},'.avi','.jpg');
filename1 = [filenamestring];

time=regexp(filename1,'(?<=-17_).*(?=.000)','match')
time=time{1}
time(regexp(time,'[-]'))=[]
time=str2double(time)

flyNumber=regexp(filename1,'(?<=fly).*(?=.000)','match')
flyNumber=strtok(flyNumber,'_')
flyNumber=cell2mat(flyNumber)
flyNumber=str2double(flyNumber)

femaleyes = '_m_';
genderfemale=isempty(strfind(filename1, femaleyes))

% Create matrix containing values you would like to save
Matrix={[flyNumber],[genderfemale],[time],[speedafinalmovmean3],[number10secondsofpolarotaxis],[polarotaxisonlydiffeneceanglesmean],[FINALPreferredevectorovertime],[FinalMeanPreferredEvectorPolarotaxisonly]}



% Continue plotting
box={['Number of 10 Seconds of Polarotaxis=' num2str(number10secondsofpolarotaxis) '/30'],['Mean Angle Diff during Polarotaxis=' num2str(polarotaxisonlydiffeneceanglesmean) '°'],['Preferred Evector=' num2str(FinalMeanPreferredEvectorPolarotaxisonly) '°'],['Time= ' num2str(time)]}
text(-180.519566052228, -200.662920634921, 0, box)
ax.FontSize = 5

subplot(5,3,12)
bar(speedafinaljumpmeanpolarotaxisonlymean2)
ax = gca
ylabel('Mean diff per 10s pol');
ylim([0 90]);
xlim([0 30])
ax.FontSize = 5
ax.XTick = [0 30]
ax.YTick = [0 45 90]

subplot(5,3,15)
bar(a4speedafinaljumpmeanpolarotaxisonlymean2)
ax = gca
ylabel('Preferred Evector');
ylim([0 180]);
xlim([0 30])
ax.FontSize = 5
ax.XTick = [0 30]
ax.YTick = [0 45 90 135 180]

subplot(5,3,15)
bar(a4speedafinaljumpmeanpolarotaxisonlymean2)
ax = gca
ylabel('Preferred Evector');
ylim([0 180]);
xlim([0 30])
ax.FontSize = 5
ax.XTick = [0 30]
ax.YTick = [0 45 90 135 180]

% Save plot as .JPG
saveas(gcf, filename1, 'jpg');
filename2 = sprintf('%s.csv',[filenamestring]);


% Save results as .CSV
lens = cellfun('length', Matrix);
rmat = NaN(length(Matrix), max(lens));
for i = 1:length(Matrix)
    rmat(i, 1:lens(i)) = Matrix{i};
end
rmat=transpose(rmat)

csvwrite(filename2,rmat)
    end

