skeletonAlignmentTrial2.m

% Import Fluorescent traces from ROIs, with information about the plane of origin (nxm: n ROI number, m is total frame number/4, nx1: rows are plane the ROI is in, 4xm: m number of frames in ROI data, 4 planes)


% TimeFrame matrix: ()
% 
% 
% 
% Bonsai 2P file: FrameTime 

% Concatenate files
options.OPhysPath = 'Z:\ibn-vision\DATA\SUBJECTS\M25135\OPhys\20260219\';
allFilesinFolder = dir(options.OPhysPath);
allFilesinFolder = struct2table(allFilesinFolder);

contrastSessions = {'Contrast_'};
firstContrastSession = {'_00001_'};

contrastFilesIdx = find(contains(allFilesinFolder.name,contrastSessions));
contrastFiles = allFilesinFolder.name(contrastFilesIdx);

firstContrastFilesIdx = find(contains(contrastFiles,firstContrastSession));
secondContrastFilesIdx = find(~contains(contrastFiles,firstContrastSession));

firstContrastFiles = contrastFiles(firstContrastFilesIdx);
secondContrastFiles = contrastFiles(secondContrastFilesIdx);

% Name the Tiff files to create, because we will append to them - use the
% 'a' option 



for thisContrastFile = 10:11
    thisTiff = tiffreadVolume(append(options.OPhysPath,firstContrastFiles{thisContrastFile}));
    tPlane1 = thisTiff(:,:,1:4:end);
    tPlane2 = thisTiff(:,:,2:4:end);
    tPlane3 = thisTiff(:,:,3:4:end);
    tPlane4 = thisTiff(:,:,4:4:end);

    for thisSlice = 1:size(tPlane1,3)
        % Create the file to append to
        plane1Tiff = Tiff("Z:\ibn-vision\DATA\SUBJECTS\M25135\OPhys\20260219\ConcatenatedPlanes\M25135_Contrast_20260219_00001_plane12.tif",'a');

        % Set parameters for writing image
        tagstruct.ImageLength = 512;
        tagstruct.ImageWidth = 512;
        tagstruct.Photometric = Tiff.Photometric.MinIsBlack;
        tagstruct.PlanarConfiguration = Tiff.PlanarConfiguration.Chunky;
        tagstruct.BitsPerSample = 16;
        tagstruct.SampleFormat = Tiff.SampleFormat.Int;
        tagstruct.SamplesPerPixel = 1;

        % Append
        setTag(plane1Tiff,tagstruct)
        write(plane1Tiff,tPlane1(:,:,thisSlice))
        clear plane1Tiff tagstruct
    end
    % Clear variables

    clear thisTiff tPlane1 tPlane2 tPlane3 tPlane4 

end

%%
% Load 2P Bonsai info 
twoPLog = readtable("Z:\ibn-vision\DATA\SUBJECTS\M25135\Bonsai\20260219\M25135_Contrast_20260219_00001_2P2026-02-19T12_46_39.csv");

% Load ROI 
ROImatrix = readtable("Z:\ibn-vision\DATA\SUBJECTS\M25135\Processed2P\20260219\M25135_Contrast_20260219_00001_plane1_ROImeasurement_combined.csv");
ROImatrix = removevars(ROImatrix,'Var1');
ROImatrix = table2array(ROImatrix);
ROImatrix = ROImatrix';



% plane Info 
planeMatrix = repmat(1,size(ROImatrix,1),1);

% Split twoPLog TwoPFrametime by plane 
planeFrameTime = twoPLog.TwoPFrameTime(1:4:end,:);

% Interpolate ROI measurements to full 2pFrameTime 
for thisROI = 1:height(ROImatrix)
    newROImatrix(thisROI,:) = interp1(planeFrameTime,ROImatrix(thisROI,:),twoPLog.TwoPFrameTime);
end

% Calculate deltaF
F0 = prctile(newROImatrix,10,2);
denominator = repmat(F0,1,length(newROImatrix));
deltaFOverF = newROImatrix-denominator;
deltaFOverF = deltaFOverF./denominator;

% Remove NaNs that were generated by the interpolation
idxNaN = find(isnan(deltaFOverF(:,:)));
%deltaFOverF(:,idxNaN) = [];


% Load Wheel 
WheelDat = readtable("Z:\ibn-vision\DATA\SUBJECTS\M25135\Bonsai\20260219\M25135_Contrast_20260219_00001_Wheel2026-02-19T12_46_44.csv");
newWheelDat = interp1(WheelDat.ArduinoTime,WheelDat.Wheel,twoPLog.TwoPFrameTime);

% Pupil Data 
EyeTimestamps = readtable("Z:\ibn-vision\DATA\SUBJECTS\M25135\EyeTracking\20260219\M25135_Contrast_20260219_00001_EyeCamTimeStamps2026-02-19T12_46_36.csv");
% Remove duplicates of samples in EyeTimeStamps
[tt,tp] = unique(EyeTimestamps.ArduinoTime);
EyeTimestamps = EyeTimestamps(tp,:);

EyeData = readtable("Z:\ibn-vision\DATA\SUBJECTS\M25135\EyeTracking\20260219\M25135_Contrast_20260219_00001_EyeCamLog_2026-02-19T12_46_32.csv");
EyeTimestamps.newArea = interp1(EyeData.Item1_eyeMillisSinceStartOfDay,EyeData.Item2_Area,EyeTimestamps.EyeCamTime,'linear','extrap');

uSyncEye = unique(EyeTimestamps.LastSyncPulseTime);
uSyncTwoP = unique(twoPLog.LastSyncPulseTime);
EyeTimestamps.newArduinoTime = align2PSyncPulses(uSyncEye,uSyncTwoP,EyeTimestamps.ArduinoTime);
newPupilArea = interp1(EyeTimestamps.newArduinoTime,EyeTimestamps.newArea,twoPLog.TwoPFrameTime,'linear','extrap');


% Load PD 

PD = readtable("Z:\ibn-vision\DATA\SUBJECTS\M25139\Bonsai\20260220\M25139_Contrast_20260220_00001_Photodiode2026-02-20T13_00_38.csv");

% Find up and downphases from PD data 
options.TwoPSampleRate = 30;
% Interpolate to twoPFrameTime 
newPD.PDOutput = interp1(PD.ArduinoTime, PD.PDOutput, twoPLog.TwoPFrameTime);

[upPhases, downPhases] = PDUpDownPhases(newPD.PDOutput,options);

% Load stim info 
stimInfo = readtable("Z:\ibn-vision\DATA\SUBJECTS\M25135\Bonsai\20260219\M25135_Contrast_20260219_00001_stimEvents2026-02-19T12_46_39.csv");
stimInfo = renamevars(stimInfo,"Var5",'ContrastLevel');

% Find trials 
columnIdx = [-8:40];
stimFrame = [];
for thisStim = 1:size(upPhases,1)
    [t,idx] = min(abs(twoPLog.TwoPFrameTime(upPhases(thisStim,1))-twoPLog.TwoPFrameTime(1:end-3)));
    stimFrame(thisStim,1) = idx;
end

lowContrastFrameTime = stimFrame(stimInfo.ContrastLevel<0.5);
highContrastFrameTime = stimFrame(stimInfo.ContrastLevel>0.5);


lowContrastTrials = repmat(lowContrastFrameTime,1,length(columnIdx)) + repmat(columnIdx,length(lowContrastFrameTime),1);
highContrastTrials = repmat(highContrastFrameTime,1,length(columnIdx)) + repmat(columnIdx,length(highContrastFrameTime),1);
tdat = newROImatrix(4,:);

deltaFOverF = calculateDeltaF(newROImatrix,twoPLog.TwoPFrameTime,[-5 0],'deltafOverf');
tdat = mean(deltaFOverF,1);
lowContrastData = tdat(lowContrastTrials);
highContrastData = tdat(highContrastTrials);

subplot(421)
imagesc(columnIdx,1:size(lowContrastData,1),lowContrastData)
subplot(423)
imagesc(columnIdx,1:size(highContrastData,1),highContrastData)
subplot(223)
plot(columnIdx,mean(lowContrastData,1),'b'); hold on 
plot(columnIdx,mean(highContrastData,1),'r')