B01_getTraj.m

%%
% HT-LMD - High-Throughput Learning, Memory, and Drug-Testing Platform for Small Aquatic Models
% Copyright (C) 2025 Gokul Rajan
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <https://www.gnu.org/licenses/>.

%%
clear all; clc;


%%%%%%%%%%%%%%%%%%%%%%%%% IMPORTANT
% only use the coords flipping in this script if you used stimulus flipping
% in yuour experimental
% protocol.

% if you use coord fliiping, make sure the variable addStim is set correctly based on your stimulus flipping time in the protocol.
% check the stimulus protocol file from athe day of experiment if there is
% any confusion
%%%%%%%%%%%%%%%%%%%%%%%%%


%% %%%%%%%%%%%%%%%%%%%%%% how to navigate

%B01_getTraj -> B02_calAllParams --> B03_cleanData --> run tests like B04_stats

%%%%%%%%%%%%%%%%%%%%%%%%%
%%

camscale_px_per_mm =8.4;
data_rate_hz=1;

if ismac
addpath(genpath('/Users/gokulrajan/Documents/MATLAB/ZebranalysisSystem/zebraFunctions')); 
rootDir='F:\GR\Dc_LnM_data';
elseif ispc
    addpath(genpath('G:\GR\matlab_GR\usefulFunctions'));
    addpath(genpath('G:\GR\matlab_GR\LnM-main'));
    addpath(genpath('F:\GR\Dc_LnM_data'));
    rootDir='F:\GR\Dc_LnM_data\9-arena-data';
end

%% get the calibration points here

% there will always be a "calib" folder inside the "xxxx-xx-xx" experiment
% folder
% irrespective of the number of imgs or the file name, this peogram will
% pick one img and provide a GUI interface to select 9 points of interest

disp('select the file for calibration');

if ismac
myCal = uigetdir(rootDir);
if ~isdir(myCal)
    uiwait(warndlg('The selected directory does not exist'));
    return;
end
calibfilePattern = fullfile(myCal,'*.png');
mycalibFiles = dir(calibfilePattern); %all files
cd(mycalibFiles(1).folder)

elseif ispc
    myCal = uigetdir(rootDir);
if ~isdir(myCal)
    uiwait(warndlg('The selected directory does not exist'));
    return;
end

calibfilePattern = fullfile(myCal,'*.png');
mycalibFiles = dir(calibfilePattern); %all files
cd(mycalibFiles(1).folder);

end

disp('select the 9 mid points');
imshow(mycalibFiles(1).name);
[xc, yc] = ginput(9);
close all;
disp('mid points selected!');
pause(5);

%% assign xc and yx to the respective arena in the arena strcuture

for aa=1:length(xc)
arena(aa).center=[xc(aa),yc(aa)];
end

%% select the folder to analyse vdo

if ismac
myDir = uigetdir(rootDir);
if ~isdir(myDir)
    uiwait(warndlg('The selected directory does not exist'));
    return;
end
filePattern = fullfile(myDir,'*.avi');
myFiles = dir(filePattern); %all txt files
cd(myFiles(1).folder);

elseif ispc
    myDir = uigetdir(rootDir);
if ~isdir(myDir)
    uiwait(warndlg('The selected directory does not exist'));
    return;
end
filePattern = fullfile(myDir,'*.avi');
myFiles = dir(filePattern); %all avi files
cd(myFiles(1).folder);
end

%%
for ff = 1:length(myFiles) % go through each file
    tic;
    
    dirName = myDir;
    fileName = myFiles(ff).name;
    inputName = fullfile(myDir,fileName);

obj = VideoReader(fileName);

%% read video and display

% read the file
mov1=cell(1,obj.NumFrames);

k = 1;
while hasFrame(obj)
    mov1{k} = rgb2gray(readFrame(obj));
    k = k+1;
end

mov3=cell2mat(reshape(mov1,1,1,[]));

clearvars mov1;
% imshow3D(mov3);

%% for subsequent loops

for aa=1:length(xc)
arena(aa).center=[xc(aa),yc(aa)];
end

% for aa=1:length(arena)
% [xc(aa)]=arena(aa).center(1);
% [yc(aa)]=arena(aa).center(2);
% end

%% preprocessing

% median image
medImage=median(mov3(:,:,1:10:300),3); 
% imshow(medImage);


%% roi

disp('select the 3 LEDs');
% slecte the leds
cs1=roipoly(medImage);
cs2=roipoly(medImage);
csFull=roipoly(medImage);

disp('LEDs selected!');
disp('now select the arenas');
pause(5);

%% select the smaller arenas

imshow(medImage);
% Loop over the number of ROIs
for i = 1:9
    % Let the user draw a rectangle on the image
    title(sprintf('Draw ROI %d and double-click when done.', i))
    roi = drawrectangle;

    % Wait for the user to finish drawing the rectangle
    wait(roi);

    % Get the position of the rectangle
    roiPosition = roi.Position;

    % Create a binary mask that is the same size as the image
    mask = false(size(medImage));

    % Set the pixels inside the rectangle to 1
    mask(round(roiPosition(2)):(round(roiPosition(2))+round(roiPosition(4))), ...
         round(roiPosition(1)):(round(roiPosition(1))+round(roiPosition(3)))) = true;
    
    % Save the mask in the arena structure
    arena(i).roi = mask;
end

disp('checkpoint 1'); %uiwait();

%% extract the stim
addStim= 300;
mov4= medImage-mov3;

maskT1=cs1;hold on;
cs1tmp=(mov3.*uint8(maskT1));


maskT2=cs2;hold on;
cs2tmp=(mov3.*uint8(maskT2));


maskT3=csFull;
%csfulltmp;%(mov3.*uint8(maskT3));

maskT=maskT1+maskT2+maskT3;
baselineImg=(medImage.*uint8(maskT));

cs1tmp_=[0; round(diff(squeeze(mean(mean(cs1tmp-baselineImg,1),2)))/norm(diff(squeeze(mean(mean(cs1tmp-baselineImg,1),2))),inf))];
cs2tmp_=[0; round(diff(squeeze(mean(mean(cs2tmp-baselineImg,1),2)))/norm(diff(squeeze(mean(mean(cs2tmp-baselineImg,1),2))),inf))];
%csfulltmp_=[0; round(diff(squeeze(mean(mean(csfulltmp-baselineImg,1),2)))/norm(diff(squeeze(mean(mean(csfulltmp-baselineImg,1),2))),inf))];

cs1_start_idx=find(cs1tmp_==1);
cs1_end_idx=cs1_start_idx+addStim;

cs2_start_idx=find(cs2tmp_==1);
cs2_end_idx=cs2_start_idx+addStim;

clearvars mov1 mov2 mov3;
%% go through each ROI per frame here
% Pre-allocate the X and Y fields in the arena structure
for whichArena = 1:length(arena)
    arena(whichArena).Xtmp = nan(size(mov4,3),1);
    arena(whichArena).Ytmp = nan(size(mov4,3),1);

    arena(whichArena).X = nan(size(mov4,3),1);
    arena(whichArena).Y = nan(size(mov4,3),1);
end

% Loop over frames
for ikk = 1:size(mov4,3)
    
    % Loop over arenas
    for whichArena = 1:length(arena)

        maskROI = arena(whichArena).roi;
        mask=uint8(maskROI);
        img_real = mov4(:,:,ikk).*mask;
        img = img_real;
        
        blob_img= (imgaussfilt(uint8(img),5));
        blob_img=imbinarize(blob_img,0.1);
        SE = strel('diamond',10);
        blob_img=imdilate(blob_img,SE);
        labeledImage = bwlabel(blob_img, 8);
        blobMeasurements = regionprops(labeledImage, img, 'all');  
        numberOfBlobs = size(blobMeasurements, 1);
        allBlobCentroids = [blobMeasurements.Centroid];
        
        % Check if any blobs were found
        if isempty(blobMeasurements)
            % If no blobs were found, assign NaN
            arena(whichArena).Xtmp(ikk) = NaN;
            arena(whichArena).Ytmp(ikk) = NaN;
        else
            % Find the largest blob
            [~, idx] = max([blobMeasurements.Area]);
            
            % Assign the X and Y coordinates of the largest blob
            arena(whichArena).Xtmp(ikk) = blobMeasurements(idx).Centroid(1);
            arena(whichArena).Ytmp(ikk) = blobMeasurements(idx).Centroid(2);
        end

    end

    clearvars -except addStim whichArena xc yc myDir myFiles fileName filePattern rootDir ff mov4 ikk X Y cs2_start_idx cs1_start_idx cs2_end_idx cs1_end_idx arena
   
    %disp((ikk/size(mov4,3))*100);
end

disp('checkpoint 2'); %uiwait();


    %% mirror the coords for the time period when the stimulus was flipped
for whichArena=1:size(arena,2)
    flip_start1 = cs2_start_idx(1); % The stimulus was flipped
    flip_end1 = flip_start1+addStim; % The stimulus flipped period %addStim

    % Flip the x-coordinates during the times when the stimulus was flipped
    flipped_data = arena(whichArena).Xtmp(:);
    flipped_data(flip_start1:flip_end1) = 2*arena(whichArena).center(1) - flipped_data(flip_start1:flip_end1);

      if length(cs2_start_idx)>1
    flip_start2 = cs2_start_idx(2); % The stimulus was flipped again at 900 seconds
    flip_end2 = flip_start2+addStim; % The stimulus was flipped back at 1200 seconds
    flipped_data(flip_start2:flip_end2) = 2*arena(whichArena).center(1) - flipped_data(flip_start2:flip_end2);
      end

    arena(whichArena).X = flipped_data;
    arena(whichArena).Y= arena(whichArena).Ytmp(:);
    clearvars flipped_data;
end

%%
toc; % time for processing one video
save(strcat(myDir,fileName,'.mat'),'cs1_start_idx','cs2_start_idx', 'arena','-v7.3');
clearvars -except addStim myDir myFiles xc yc filePattern rootDir ff

end