diff --git a/+ciapkg/+behavior/importDeepLabCutData.m b/+ciapkg/+behavior/importDeepLabCutData.m index e606bdb..22dc779 100644 --- a/+ciapkg/+behavior/importDeepLabCutData.m +++ b/+ciapkg/+behavior/importDeepLabCutData.m @@ -1,4 +1,4 @@ -function [mainTbl] = importDeepLabCutData(inputPath,varargin) +function [mainTbl, mainTensor, mainCell] = importDeepLabCutData(inputPath,varargin) % [outputTable] = importDeepLabCutData(inputPath,varargin) % % Loads DeepLabCut processed CSV files and loads into a table or structure for use by other functions. @@ -10,7 +10,9 @@ % inputPath - Str: path to DLC file. % % Outputs - % mainTbl - Struct with fields named after body parts. Each field is [likelihood x y] matrix with rows equal to frames of the movie. + % mainTbl - Struct: with fields named after body parts. Each field is [likelihood x y] matrix with rows equal to frames of the movie. + % mainTensor - Tensor: format [nFeatures nDlcOutputs nFrames]. + % mainCell - Cell: format {1 nFeatures} with each cell containing a [likelihood x y] matrix with rows equal to frames of the movie. % % Options (input as Name-Value with Name = options.(Name)) % % DESCRIPTION @@ -18,6 +20,7 @@ % Changelog % 2022.03.14 [01:47:04] - Added nested and local functions to the example function. + % 2024.02.19 [09:25:58] - Added conversion to tensor and cell to the import function from other functions. % TODO % @@ -47,18 +50,22 @@ % dlcTablePath = fullfile(rootAnalysisPath,dlcFileCSV{vidNo}); t1 = readtable(inputPath); t2 = readlines(inputPath); - bodyPartListMain = strsplit(t2{2},','); - bodyPartList = unique(bodyPartListMain); - nBodyParts = length(bodyPartList)-1; + featuresListMain = strsplit(t2{2},','); + featuresList = unique(featuresListMain); + nfeaturess = length(featuresList)-1; if options.dispPlots==1 figure; set(gcf,'Color','k') end - for partX = 1:nBodyParts - partName = bodyPartListMain{(partX-1)*3+4}; - if partX==nBodyParts + mainCell = {}; + nFrames = length(t1{:,(1-1)*3+2}); + mainTensor = NaN([nfeaturess 3 nFrames]); + + for partX = 1:nfeaturess + partName = featuresListMain{(partX-1)*3+4}; + if partX==nfeaturess fprintf('%s.\n',partName) else fprintf('%s | ',partName) @@ -66,7 +73,14 @@ conf1 = t1{:,(partX-1)*3+4}; dlcX = t1{:,(partX-1)*3+2}; dlcY = t1{:,(partX-1)*3+3}; - mainTbl.(partName) = [conf1(:) dlcX(:) dlcY(:)]; + + mainMatrix = [conf1(:) dlcX(:) dlcY(:)]; + + mainTbl.(partName) = mainMatrix; + + mainCell{partX} = mainMatrix; + + mainTensor(partX,:,:) = mainMatrix'; if options.dispPlots==1 subplot(3,3,partX) @@ -92,7 +106,25 @@ % outputs = ; end end -function [outputs] = localfxn_exampleFxn(arg) - % Always start local functions with "localfxn_" prefix. - % outputs = ; -end \ No newline at end of file +% function [outputs] = localfxn_exampleFxn(arg) +% opts.nFrames = size(inputMovie,3); +% featurePts = cell([1 opts.nFrames]); +% opts.nFields = fieldnames(dlcData); +% featurePtsTensor = NaN([length(opts.nFields) 3 opts.nFrames]); + +% for i = 1:opts.nFrames +% conf1 = []; +% dlcX = []; +% dlcY = []; +% for partX = 1:length(opts.nFields) +% thisField = opts.nFields{partX}; +% dlcDataNew.(thisField) = vertcat(dlcData.(thisField)); +% conf1(partX) = dlcDataNew.(thisField)(i,1); +% dlcX(partX) = dlcDataNew.(thisField)(i,2); +% dlcY(partX) = dlcDataNew.(thisField)(i,3); +% end +% featurePts{i} = [conf1(:) dlcX(:) dlcY(:)]; +% featurePtsTensor(:,:,i) = featurePts{i}; +% end +% disp('Done!') +% end \ No newline at end of file diff --git a/+ciapkg/+classification/matchObjBtwnTrials.m b/+ciapkg/+classification/matchObjBtwnTrials.m index 6b36318..2f12591 100644 --- a/+ciapkg/+classification/matchObjBtwnTrials.m +++ b/+ciapkg/+classification/matchObjBtwnTrials.m @@ -26,6 +26,7 @@ % 2019.10.29 [13:07:18] - Added support for sparse (ndSparse) array inputs. % 2019.12.05 [10:37:17] - Fix for cases in which cells are lost when cropping all movies to the same size before starting cross-session alignment. % 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package. + % 2023.05.09 [12:33:23] - Made ability to switch between different motion correction and image alignment methods user visible. % notes % the cell array of traces allows us to have arbitrary numbers of trials to align automatically, % TODO @@ -70,6 +71,14 @@ options.turboregZeroThres = 1e-6; % Binary: 1 = check correlation matches visually options.checkImageCorr = 0; + % Str: Register images using 'imtransform' or 'imwarp' (Matlab) or 'transfturboreg' (C) + options.registrationFxn = 'transfturboreg'; + % Str: motion correction algorithm. + % 'turboreg' - TurboReg as developed by Philippe Thevenaz. + % 'normcorre' - NoRMCorre as developed by several authors at Flatiron Institute. + % 'imregdemons' - Displacement field alignment based on imregdemons. + options.mcMethod = 'turboreg'; + options = getOptions(options,varargin); %======================== % obtain trial stats @@ -124,6 +133,10 @@ % ======= RegisTypeFinal = options.RegisTypeFinal; % turboreg options. + ioptions.mcMethod = options.mcMethod; + % ioptions.registrationFxn = 'imtransform'; + ioptions.registrationFxn = options.registrationFxn; + ioptions.meanSubtract = 0; ioptions.complementMatrix = 0; ioptions.normalizeType = []; @@ -141,8 +154,6 @@ ioptions.cropCoords = []; ioptions.closeMatlabPool = 0; ioptions.removeEdges = 0; - % ioptions.registrationFxn = 'imtransform'; - ioptions.registrationFxn = 'transfturboreg'; % ======= % turboreg all object maps to particular trial's map if issparse(objectMap{options.trialToAlign}) @@ -211,10 +222,13 @@ ioptions.closeMatlabPool = 0; ioptions.meanSubtract = 0; ioptions.normalizeType = 'divideByLowpass'; - ioptions.registrationFxn = 'transfturboreg'; + % ioptions.registrationFxn = 'transfturboreg'; ioptions.removeEdges = 0; ioptions.cropCoords = []; + ioptions.mcMethod = options.mcMethod; + ioptions.registrationFxn = options.registrationFxn; + for switchNo = 1:length(switchStrArray) switch switchStrArray{switchNo} case 'optional' @@ -352,6 +366,7 @@ OutStruct.registrationCoords = registrationCoords; OutStruct.coords = coords; OutStruct.inputImages = inputImages; + % OutStruct.inputOptions = options; if ~isempty(options.inputSignals) OutStruct.inputSignals = options.inputSignals; end diff --git a/+ciapkg/+demo/cross_session_motion_correction.m b/+ciapkg/+demo/cross_session_motion_correction.m new file mode 100644 index 0000000..2886b7a --- /dev/null +++ b/+ciapkg/+demo/cross_session_motion_correction.m @@ -0,0 +1,482 @@ +% Cross-session motion correction method (CS-MCM) +% Run each cell in succession to process data. +% Assumes that all the imaging data has near equivalent µm/pixel. +% https://doi.org/10.1101/2023.05.22.541477 +% Biafra Ahanonu +% started: 2021.04.16 [14:16:27] + +%% ========================================== +% Load parameters +% Str: path to output data +outputMovieFolder = 'PATH'; +% Str: Path to root folder containing all sub-folders with imaging data. +rootPathHere = 'PATH'; +% [OPTIONAL] Str: if all imaging sessions are in the same folder +rootPathLoad = ''; + +% Str: name of the animal subject whose data is being analyzed +subjectName = 'mouse01'; +% Str: regular expression used to identify folders +folderMatchStr = 'day'; +% Str: regular expression used to find files in each imaging session folder. +fileRegExpHere = 'epi01.*.tif'; + +% Int: Reference frame for within session correction +refFrame = 2; +% Int: Session to use for reference for cross-session alignment +refSession = 102; +% Int vector: Frames to use in the movies, leave blank to load all frames +frameRange = 100:102; + +% 1 = match the size of the images by upscaling to largest image size +sizeMatchFlag = 1; +% Which file number in the folder to load +fileNoLoadForce = 1; +% Str: '1P' or '2P' +imgType = '1P'; +% Int vector: List of sessions to flip dimension #1 (rows). +flipDim1 = []; +% Int vector: List of sessions to flip dimension #2 (columns). +flipDim2 = []; +% Int vector: List of sessions to rotate 90 degrees. +rot90Dim = []; +% Int vector: List of sessions to skip +skipSessions = []; + +% Get the current date-time string for file saving +currentDateTimeStr = datestr(now,'yyyy_mm_dd_HHMMSS','local'); +outputMovieNameStart = [currentDateTimeStr '_' subjectName '_crossDay_day']; + +%% ========================================== +% Start parallel pool +ciapkg.io.manageParallelWorkers(20,'forceParpoolStart',1); + +%% ========================================== +% Load all the session information where each folder contains a different imaging session +try + folderList = cellfun(@(x) ciapkg.io.getFileList(x,folderMatchStr),rootPathLoad,'UniformOutput',0); + folderList = [folderList{:}]; +catch + folderList = ciapkg.io.getFileList(rootPathLoad,folderMatchStr); +end +nFolders = length(folderList); +dayList = []; +folderNoList = []; +runList = []; +runListNew = []; +folderListNew = {}; +for folderNo = 1:nFolders + folderPath = folderList{folderNo}; + movieLoadPath = ciapkg.io.getFileList(folderPath,fileRegExpHere); + % Get the "day" number, e.g. day01 or day04 will yield 1 and 4 for later identification. + dayNum = str2num(cell2mat(strrep(regexp(folderPath,[folderMatchStr '\d+'],'match'),folderMatchStr,''))); + if ~isempty(movieLoadPath)&&any(skipSessions==dayNum)==0 + dayList(end+1) = dayNum; + folderListNew{end+1} = folderPath; + fprintf('ADDING Day %d | \n',dayList(end)) + runList(folderNo) = 1; + runListNew(end+1) = 1; + else + fprintf('SKIPPING Day %d | \n',dayNum) + runList(folderNo) = 0; + end +end +dayList = dayList'; +fprintf('\n'); + +%% ========================================== +% Load all the session data +nFolders = length(folderListNew); +gMovieRaw = cell([nFolders 1]); +errorLog = zeros([1 nFolders]); +cNum = 1; +WaitMessage = parfor_wait(nFolders, 'Waitbar', true); +parfor folderNo = 1:nFolders + try + WaitMessage.Send; + disp('==========================================') + fprintf('Loading day %d/%d (%d/%d).\n',folderNo,dayList(end),folderNo,nFolders); + if runListNew(folderNo)==0 + disp('No data for this day, skipping...') + continue; + end + folderPath = folderListNew{folderNo}; + + if ~isempty(rootPathHere) + dayCheck = dayList(folderNo); + if dayCheck<10 + dayCheck = ['0' num2str(dayCheck)]; + else + dayCheck = num2str(dayCheck); + end + movieLoadPath = ciapkg.io.getFileList(rootPathHere,[folderMatchStr dayCheck]); + else + movieLoadPath = []; + end + + % Priority to manually selected TIF else use from raw folder + if isempty(movieLoadPath) + movieLoadPath = ciapkg.io.getFileList(folderPath,fileRegExpHere); + frameRangeTmp = frameRange; + else + frameRangeTmp = 1:length(frameRange); + end + + % Load imaging data if avaliable in folder + if ~isempty(movieLoadPath) + warning off + fileNoLoad = 1; + if ~isempty(fileNoLoadForce) + fileNoLoad = fileNoLoadForce; + end + disp(fileNoLoad) + if length(frameRangeTmp)>1 + gMovieRaw{folderNo} = ciapkg.io.loadMovieList(movieLoadPath{fileNoLoad},'frameList',frameRangeTmp); + else + gMovieRaw{folderNo} = ciapkg.io.readFrame(movieLoadPath{fileNoLoad},frameRangeTmp); + end + warning on + end + catch err + disp(repmat('@',1,7)) + disp(getReport(err,'extended','hyperlinks','on')); + disp(repmat('@',1,7)) + try + errorLog(folderNo) = 1; + catch + end + end +end + +WaitMessage.Destroy; + +% Save a backup of original data in case need to go back to it +gMovieRawOriginal = gMovieRaw; +nSessions = length(gMovieRaw); +disp('Done loading movies') + +%% ========================================== +% Plot to make sure loaded data is usable +localfun_plotInputMovies(gMovieRaw,nSessions,refSession,refFrame,dayList,6,0,1,0); + +%% ========================================== +% Match size of all images across all movies +gMovieRaw = ciapkg.image.matchImageDimensions(gMovieRaw,'modType','pad'); +disp('Done!') + +%% ========================================== +% Prepare for run, copy raw data to a new matrix +nSessions = length(gMovieRaw); +gMovie = gMovieRaw; +disp('Done gMovie = gMovieRaw') + +%% ========================================== +% Plot a frame from each to confirm data loaded properly +cmapRange = [100 200]; % [min max], pixel intensity range to plot data. +cmapRange = 1; % Plot the 15th and 99th percentile for min and max. +localfun_plotInputMovies(gMovie,nSessions,refSession,refFrame,dayList,10,0,cmapRange,0); + +%% ========================================== +% [OPTIONAL] For certain datasets, rotate all images +rot90Dim = 1:length(dayList); +disp('Done assigning files to rotate') + +%% ========================================== +% [OPTIONAL] For certain datasets, flip all images in dimension #1 (rows) +flipDim1 = 1:length(dayList); +disp('Done assigning files to rotate') + +%% ========================================== +% [OPTIONAL] For certain datasets, flip all images in dimension #2 (columns) +flipDim2 = 1:length(dayList); +disp('Done assigning files to rotate') + +%% ========================================== +% [OPTIONAL] Flip movies in x or y as needed due to camera or other changes +% Index of session to flip in each dimension +for i = 1:nSessions + fprintf('Checking %d/%d\n',i,nSessions); + if any(flipDim1==i) + fprintf('Flip dim1 %d/%d\n',i,nSessions); + gMovie{i} = flip(gMovie{i},1); + end + if any(flipDim2==i) + fprintf('Flip dim2 %d/%d\n',i,nSessions); + gMovie{i} = flip(gMovie{i},2); + end + if any(rot90Dim==i) + fprintf('Rotate -90 deg %d/%d\n',i,nSessions); + gMovie{i} = rot90(gMovie{i},1); + end +end + +%% ========================================== +% Convert all input data to single-precision (32-bit) arrays. +gMovie = cellfun(@single,gMovie,'UniformOutput',0); + +%% ========================================== +% [OPTIONAL] Manually correct motion across movies to a reference frame +% If there is significant motion (e.g. translation or rotation) or the field of view is flipped, can improve automated registration. +[gMovieTmp, outStructTmp] = ciapkg.motion_correction.computeManualMotionCorrection(gMovie,... + 'registerUseOutlines',0,'cellCombineType','mean','gammaCorrection',1.6,'refFrame',refSession,'translationAmt',1,'includeImgsOutputStruct',0); +disp('Done!') + +%% [OPTIONAL] Save output +outputMoviePath = fullfile(outputMovieFolder,[outputMovieNameStart num2str(dayList(1)) '-' num2str(dayList(nSessions)) '_postManual.mat']); +disp(['Saving: ' outputMoviePath]) +save(outputMoviePath,'gMovieTmp',"outStructTmp"); +disp('DONE!') + +%% Transfer over, safer in case manual exists early. +gMovie = gMovieTmp; +clear gMovieTmp; +clear outStructTmp +disp('DONE!') + +%% ========================================== +% Register each movie to itself to reduce complexity of cross-session motion correction +[~,cropCoordsInSession] = ciapkg.image.cropMatrix(gMovie{refSession},'cropOrNaN','crop','inputCoords',[]); +parfor i = 1:nSessions + disp('=====================') + fprintf('Session %d/%d:\n',i,nSessions); + gMovie{i} = ciapkg.motion_correction.turboregMovie(gMovie{i},'refFrame',refFrame,... + 'RegisType',1,'normalizeType','matlabDisk','cropCoords',cropCoordsInSession,... + 'SmoothX',10,'SmoothY',10,'meanSubtract',1,'showFigs',0); +end +disp('Done within-session motion correction') + +%% ========================================== +% Automatic handling of large shifts in the movies across days, to be followed by motion correction to handle smaller shifts. +% First take mean image for each session to conduct initial motion correction. + +try close(9); catch; end + +% ======= +% TurboReg settings +ioptions.SmoothX = 30;%10 +ioptions.SmoothY = 30;%10 +ioptions.minGain = 0.0; +ioptions.Levels = 6; +ioptions.Lastlevels = 1; +ioptions.Epsilon = 1.192092896E-07; +ioptions.zapMean = 0; +% +ioptions.turboregRotation = 1; +ioptions.parallel = 1; +ioptions.closeMatlabPool = 0; +ioptions.meanSubtract = 1; +% ioptions.normalizeType = 'divideByLowpass'; +ioptions.normalizeType = 'matlabDisk'; +ioptions.matlabdiskR1 = 20; +ioptions.matlabdiskR2 = 10; +ioptions.registrationFxn = 'transfturboreg'; +ioptions.removeEdges = 0; +ioptions.complementMatrix = 1; % Switch to 0 for 2P data +% ======= + +% Manually get area to use for registration. +[~,inputCoords] = ciapkg.image.cropMatrix(gMovie{refSession},'cropOrNaN','crop','inputCoords',[]); +ioptions.cropCoords = inputCoords; + +% Store the motion correction coordinates. +registrationCoords = {}; + +% Reference frame is the average of the entire reference session, less sensitive to peculiarities of a single frame. +refMovieFrame = cast(mean(gMovie{refSession},3),'double'); + +inputImagesTranslated = {}; +inputImagesTranslated{refSession} = gMovie{refSession}; +gMovieTmp = gMovie; +for i = 2:nSessions + disp('=====================') + sessionStr = sprintf('Session %d/%d:\n',i,nSessions); + disp(sessionStr); + title(sessionStr); + for roundNo = 1:2 + fprintf('Registration round %d/%d:\n',roundNo,2); + if roundNo==2 + ioptions.RegisType = 2; % 2 - affine, rotation, no skew + else + ioptions.RegisType = 1; % 1 - affine, no rotation, no skew + end + ioptions.altMovieRegister = gMovieTmp{i}; + thisSessionFrame = cast(mean(gMovieTmp{i},3),'double'); + objMapsToTurboreg = cat(3,refMovieFrame,thisSessionFrame); + [gMovieTmp{i}, registrationCoords{i}] = ciapkg.motion_correction.turboregMovie(objMapsToTurboreg,'options',ioptions); + drawnow + end +end +inputImagesTranslated = gMovieTmp; +clear gMovieTmp; +disp('Done correcting large shifts motion correction') + +%% Save output +savePostStr = '_postCrossSessionMotionCorrect.mat'; +outputMoviePath = fullfile(outputMovieFolder,[outputMovieNameStart num2str(dayList(1)) '-' num2str(dayList(nSessions)) savePostStr]); +disp(['Saving...' outputMoviePath]) +save(outputMoviePath,'gMovie','inputImagesTranslated','dayList','registrationCoords'); +disp('DONE!') + +%% ========================================== +% [OPTIONAL] If need to go back to gMovie before cross-session motion correction +inputImagesTranslated = gMovie; + +%% ========================================== +% Plot frames from the motion corrected movie, check that there are not large errors in motion correction +localfun_plotInputMovies(inputImagesTranslated,nSessions,refSession,refFrame,dayList,22,0,1,0); +colormap gray + +%% ========================================== +% [OPTIONAL] Backup motion corrected data +inputImagesTranslatedTmp = inputImagesTranslated; + +%% ========================================== +% [OPTIONAL] Further manual correction to quickly correct any remaining shifts +[inputImagesTranslated, outputStruct] = ciapkg.motion_correction.computeManualMotionCorrection(inputImagesTranslated,... + 'registerUseOutlines',0,'cellCombineType','mean','gammaCorrection',1.6,'refFrame',refSession); + +%% ========================================== +% [OPTIONAL] Normalize the movie to reduce issues with cross-session displays due to background intensity variations +disp('Start normalizing movies!') +inputImagesTranslated = cellfun(@(x) ciapkg.signal_processing.normalizeVector(x,'normRange','zeroToOneSoft','prctileMax',99.5,'prctileMin',0.1),inputImagesTranslated,'UniformOutput',false); +disp('Done normalizing movies!') + +%% ========================================== +% Conduct fine and final motion correction now that large shifts have been made + +% Get location to use for rigid motion correction of the movie +[~,inputCoords] = ciapkg.image.cropMatrix(inputImagesTranslated{:},'cropOrNaN','crop','inputCoords',[]); + +% Perform the motion correction with all movies combined +gStackCorrected = ciapkg.motion_correction.turboregMovie(cat(3,inputImagesTranslated{:}),... + 'refFrame',(refSession-1)*length(frameRange)+1,'RegisType',1,'normalizeType','matlabDisk','cropCoords',inputCoords,... + 'SmoothX',50,'SmoothY',50,'meanSubtract',1); +disp('Done registering.'); + +%% ========================================== +% Play motion corrected movie +ciapkg.view.playMovie(gStackCorrected); + +%% Save output as a MAT file +savePostStr = '_gStackCorrected.mat'; +outputMoviePath = fullfile(outputMovieFolder,[outputMovieNameStart num2str(dayList(1)) '-' num2str(dayList(nSessions)) savePostStr]); +disp(['Saving...' outputMoviePath]) +save(outputMoviePath,'gStackCorrected','dayList','folderList'); +disp('DONE!') + +%% Save output as a TIFF file +savePostStr = '_gStackCorrected_.tif'; +outputMoviePath = fullfile(outputMovieFolder,[outputMovieNameStart num2str(dayList(1)) '-' num2str(dayList(nSessions)) savePostStr]); +disp(['Saving...' outputMoviePath]) +ciapkg.io.saveMatrixToFile(gStackCorrected,outputMoviePath); +disp('DONE!') + +%% ========================================== +function localfun_plotInputMovies(inputImagesTranslated,nSessions,refSession,refFrameHere,dayList,yplot,run2ndPlot,sameLookupTable,cropMovieFlag) + % Plots all input sessions as a movie montage with different methods of displaying intensity. + if cropMovieFlag==1 + [~,cropCoordsInSession] = ciapkg.image.cropMatrix(inputImagesTranslated{refSession},'cropOrNaN','crop','inputCoords',[]); + for k = 1:length(inputImagesTranslated) + inputImagesTranslated{k} = ciapkg.image.cropMatrix(inputImagesTranslated{k},'cropOrNaN','crop','inputCoords',cropCoordsInSession); + end + else + end + + subplotTmp = @(x,y,z) subaxis(x,y,z, 'Spacing', 0.005, 'Padding', 0.00, 'PaddingTop', 0.02, 'MarginTop', 0.02,'MarginBottom', 0.02,'MarginLeft', 0.01,'MarginRight', 0.01); + + xplot = ceil(nSessions/yplot); + figure; + allCaxis = []; + if length(sameLookupTable)>1 + allCaxis = sameLookupTable; + elseif sameLookupTable==1 + tmpI = cat(3,inputImagesTranslated{:}); + allCaxis = [prctile(tmpI(:),15) prctile(tmpI(:),99)]; + clear tmpI; + elseif sameLookupTable==2 + allCaxis = [115 275]; + end + linkAxesList = []; + for i = 1:nSessions + thisImg = inputImagesTranslated{i}(:,:,1); + thisImg = imresize(thisImg,0.5); + linkAxesList(i) = subplotTmp(xplot,yplot,i); + imAlpha = thisImg>1e-6|thisImg~=0; + imagesc(thisImg,'AlphaData',imAlpha); + title([num2str(dayList(i)) ' (' num2str(i) ')']) + axis image + linkAxesList(i) = gca; + box off; + axis off; + hold on; + plot(round(size(thisImg,1)/2),round(size(thisImg,2)/2),'r+') + tmpFrame = thisImg; + tmpFrame = tmpFrame(tmpFrame>1e-6); + if sameLookupTable~=0 + caxis(allCaxis); + else + caxis([prctile(tmpFrame(:),15) prctile(tmpFrame(:),99)]); + end + if mod(i,yplot)==0 + drawnow + end + drawnow + end + if sameLookupTable~=0 + disp('Linking axes') + linkaxes(linkAxesList); + hlinkHere = linkprop(linkAxesList,'CLim'); + KEY_SAVE = 'graphics_linkCLim'; + for zz = 1:length(linkAxesList) + setappdata(linkAxesList(zz),KEY_SAVE,hlinkHere); + end + end + tmpPos = get(gca,'Position'); + h=colorbar('southoutside'); + tmpPosH = get(h,'Position'); + tmpPosH = [tmpPosH(1) tmpPosH(2)-0.1 tmpPosH(3) tmpPosH(4)]; + tmpPosH(2) = 0+tmpPosH(4)++tmpPosH(4)*0.2; + disp(tmpPosH) + set(h,'Position',tmpPosH); + set(gca,'Position',tmpPos); + + set(gcf,'color',[0 0 0]); + ciapkg.view.changeFont(8,'fontColor','w') + colormap([ciapkg.view.customColormap()]); + + if run2ndPlot==0 + return; + end + % Plot ref (green) and comparison (purple) with shift + figure + refFrameHere = inputImagesTranslated{refSession}(:,:,1); + refFrameHere = ciapkg.signal_processing.normalizeVector(single(refFrameHere),'normRange','zeroToOne'); + gammaCorrectionRef = 1.5; + rgbImage = zeros([size(refFrameHere,1) size(refFrameHere,2) 3]); + + refFrameHere = imadjust(refFrameHere,[],[],gammaCorrectionRef); + for i = 1:nSessions + + mainFrameHere = inputImagesTranslated{i}(:,:,1); + mainFrameHere = ciapkg.signal_processing.normalizeVector(single(mainFrameHere),'normRange','zeroToOne'); + + mainFrameHere = imadjust(mainFrameHere,[],[],gammaCorrectionRef); + + rgbImage(:,:,1) = mainFrameHere; %red + rgbImage(:,:,2) = refFrameHere; %green + rgbImage(:,:,3) = mainFrameHere; %blue + + subplotTmp(xplot,yplot,i) + imagesc(rgbImage) + title([num2str(dayList(i)) ' (' num2str(i) ')']) + axis image + axis off; + box off; + if mod(i,yplot)==0 + drawnow + end + end + set(gcf,'color',[0 0 0]); + ciapkg.view.changeFont(10,'fontColor','w') + suptitle('Green = reference | ') +end \ No newline at end of file diff --git a/+ciapkg/+demo/displacementFields.m b/+ciapkg/+demo/displacementFields.m new file mode 100644 index 0000000..905cdc6 --- /dev/null +++ b/+ciapkg/+demo/displacementFields.m @@ -0,0 +1,43 @@ +% Displacement field motion correction example across several images +% +% Biafra Ahanonu +% started: 2022.12.14 [04:02:26] + +%% Get data paths +inputFrames = {... + {'dots_reference.tif','dots_template.tif'}; + {'spinalImaging01_ref.tif','spinalImaging01_template.tif'}; + {'testFace_ref.jpg','testFace_template.jpg'}; +}; +nTests = length(inputFrames); +mainPath = fullfile(ciapkg.getDirPkg('data'),'displacementFields'); + +%% Load data, run displacement field motion correction, display results +for testNo = 1:nTests + disp('==========') + f1 = imread(fullfile(mainPath,inputFrames{testNo}{1})); + f2 = imread(fullfile(mainPath,inputFrames{testNo}{2})); + inputMovie = cat(3,im2gray(f1),im2gray(f2)); + inputMovieTmp = inputMovie; + frameA = 1; + frameB = 2; + + if testNo==2 + inputMovieTmp = ciapkg.movie_processing.downsampleMovie(inputMovie,'downsampleDimension','space','downsampleFactor',2); + end + + [inputMovieTmp2,ResultsOutOriginal] = ciapkg.motion_correction.turboregMovie(inputMovieTmp,... + 'mcMethod','imregdemons',... + 'refFrame',1,... + 'df_AccumulatedFieldSmoothing',0.54,... % 1.3 + 'df_Niter',[500 400 200],... % [500 400 200] + 'df_PyramidLevels',3,... + 'df_DisplayWaitbar',false); + fixedTmp = inputMovieTmp(:,:,1); + movingTmp = inputMovieTmp(:,:,2); + movingReg = inputMovieTmp2(:,:,2); + Dxy = ResultsOutOriginal{2}; + + ciapkg.view.displacementFieldCompare(fixedTmp,movingTmp,movingReg,Dxy); +end +disp('Done!') \ No newline at end of file diff --git a/+ciapkg/+demo/ldmcm.m b/+ciapkg/+demo/ldmcm.m new file mode 100644 index 0000000..3b36b76 --- /dev/null +++ b/+ciapkg/+demo/ldmcm.m @@ -0,0 +1,133 @@ +% LD-MCM motion correction using deep learning and control point registration +% +% Biafra Ahanonu +% started: 2022.10.29 [13:57:45] +% Changelog +% 2024.03.05 [10:15:28] - Updated for CIAtah repo + +%% ================================= +% Settings +opts = struct; +% Int: frame to use for registering images +opts.refFrameNo = 1; +% Binary: 1 = make plots useful for debugging +opts.debugPlotFlag = 1; +% Str: full path to control point tracking (e.g. DeepLabCut) CSV +opts.dlcCsvPath = 'PATH_TO_CSV'; +% Str: full path to input movie +opts.inputMoviePath = 'PATH_TO_MOVIE_FILE'; +% Int: if movie used for tracking if downsampled X amount from movie being registered +opts.dsSpace = 1; +% Int: folder number, keep 1 for now. +opts.folderNo = 1; + +% Automatically create output file name, else directly change below. +[pathT, nameT, extT] = fileparts(opts.inputMoviePath); +opts.outputMoviePath = fullfile(pathT,[nameT '_lcmcm_detrend' extT]); +opts.outputMoviePath_dfof = fullfile(pathT,[nameT '_lcmcm_detrend_dfof' extT]); + +% This will be filled with a spatially filtered movie if user requests it +inputMovieNorm = []; +disp('Done!') + +%% ================================= +% Import DLC CSV data file. Load the dlcTensor to give to LD-MCM function as an input +[dlcStruct,dlcTensor,dlcCell] = ciapkg.behavior.importDeepLabCutData(opts.dlcCsvPath); +disp('Done!') + +%% ================================= +% Load movie +inputMovie = ciapkg.io.loadMovieList(opts.inputMoviePath,'frameList',frameList); +inputMovie = single(inputMovie); +disp('Done!') + +%% ================================= +% Overlay features on the reference frame from the input movie. +thisFrame = squeeze(inputMovie(:,:,opts.refFrameNo)); +folderInfo = ciapkg.io.getFileInfo(opts.inputMoviePath); +figure +localfxn_overlayFeaturesMovie(thisFrame,dlcTensor*opts.dsSpace,opts.refFrameNo,opts.debugPlotFlag,opts.folderNo,folderInfo) +disp('Done!') + +%% ================================= +% Get the crop coordinates for rigid registration +thisFrame = squeeze(inputMovie(:,:,opts.refFrameNo)); +[~,cropCoordsInSession{opts.folderNo}] = ciapkg.image.cropMatrix(thisFrame,'cropOrNaN','crop','inputCoords',[]); +disp('Done!') + +%% ================================= +% [OPTIONAL] Normalize movie now (normalizing after registration can lead to artifacts) +inputMovieNorm = ciapkg.movie_processing.normalizeMovie(... + inputMovie,... + 'normalizationType','lowpassFFTDivisive',... + 'freqLow',0,'freqHigh',4,... + 'bandpassMask','gaussian'); +disp('Done!') + +%% ================================= +% Run LD-MCM control point motion correction followed by rigid motion correction +% maxDist will help determine the max distance to motion correct, it can be limited by determining the maximal rostrocaudal displacement between frames +mt_inputMovie = ciapkg.motion_correction.ldmcm(inputMovie,dlcTensor,... + 'cropCoordsInSession',cropCoordsInSession{opts.folderNo},... + 'maxDist',50,... + 'refFrame',opts.refFrameNo,... + 'rotateImg',0,... + 'dsSpace',1,... + 'runPostReg',1,... + 'inputMovie2',inputMovieNorm,... + 'dimToFix',[]); +disp('Done!') + +%% ================================= +% View motion-corrected movie +ciapkg.view.playMovie(mt_inputMovie,'extraMovie',inputMovie); +disp('Done!') + +%% ================================= +% Create a border around the edges of the movie to handle uneven edges created by movement +mt_inputMovie = ciapkg.image.cropMatrix(mt_inputMovie,'cropType','auto','autoCropType','perSide'); +ciapkg.view.playMovie(mt_inputMovie,'extraMovie',inputMovie); +disp('Done!') + +%% ================================= +%% Detrend movie if needed (3rd order fit often works well) +mt_inputMovie_dt = ciapkg.movie_processing.detrendMovie(mt_inputMovie,'detrendDegree',3); + +% Make plot of mean frame intensity before and after detrending +figure; plot(squeeze(mean(mt_inputMovie,[1 2],'omitnan'))); +hold on; plot(squeeze(mean(mt_inputMovie_dt,[1 2],'omitnan'))); ylabel('Intensity'); xlabel('Frame') +box off; legend({'Raw','Detrended'}); legend boxoff; ciapkg.view.changeFont(20) +disp('Done!') + +%% ================================= +%% Compute dF/F movie +mt_inputMovie_dt_dfof = ciapkg.movie_processing.dfofMovie(mt_inputMovie_dt,'dfofType','dfof'); +disp('Done!') + +%% ================================= +% View processed and dF/F movie +ciapkg.view.playMovie(mt_inputMovie_dt_dfof,'extraMovie',mt_inputMovie_dt); +disp('Done!') + +%% ================================= +% Save motion corrected movie +ciapkg.io.saveMatrixToFile(mt_inputMovie_dt,opts.outputMoviePath); +% Save dF/F movie +ciapkg.io.saveMatrixToFile(mt_inputMovie_dt_dfof,opts.outputMoviePath_dfof); +disp('Done!') + +%% Local functions +function localfxn_overlayFeaturesMovie(thisFrame,dlcTensor,refFrameNo,debugPlotFlag,folderNo,folderInfo) + % Overlay the features for a given frame from the movie. + if debugPlotFlag==1 + imagesc(thisFrame) + axis image + axis tight + box off + colormap gray + hold on; + title(['Movie #' num2str(folderNo) ' | ref frame #' num2str(refFrameNo) 10 strrep([folderInfo.date '_' folderInfo.protocol '_' folderInfo.subjectStr '_' folderInfo.trial],'_','\_')]) + plot(dlcTensor(:,2,refFrameNo),dlcTensor(:,3,refFrameNo),'y.','MarkerSize',15); + drawnow + end +end \ No newline at end of file diff --git a/+ciapkg/+download/downloadGithubRepositories.m b/+ciapkg/+download/downloadGithubRepositories.m index a2673f4..29102a2 100644 --- a/+ciapkg/+download/downloadGithubRepositories.m +++ b/+ciapkg/+download/downloadGithubRepositories.m @@ -10,6 +10,7 @@ % 2021.02.01 [15:19:40] - Update `_external_programs` to call ciapkg.getDirExternalPrograms() to standardize call across all functions. % 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package. % 2022.01.27 [08:49:09] - If multiple repos requested but one already downloaded, previously exited before downloading the other repos. + % 2023.05.29 [16:15:18] - Automatically convert string inputs to character arrays. import ciapkg.api.* % import CIAtah functions in ciapkg package API. @@ -40,6 +41,13 @@ gitRepos = options.gitRepos; outputDir = options.outputDir; gitName = options.gitName; + + % Convert to cell character array if user only inputs a string + if ~iscell(gitNameDisp); gitNameDisp = {gitNameDisp}; end + if ~iscell(gitRepos); gitRepos = {gitRepos}; end + if ~iscell(outputDir); outputDir = {outputDir}; end + if ~iscell(gitName); gitName = {gitName}; end + nRepos = length(outputDir); % If forcing an update, make sure to remove all external programs from the path diff --git a/+ciapkg/+download/example_downloadTestData.m b/+ciapkg/+download/example_downloadTestData.m index 9bd4a4f..5bc2040 100644 --- a/+ciapkg/+download/example_downloadTestData.m +++ b/+ciapkg/+download/example_downloadTestData.m @@ -13,6 +13,7 @@ % 2021.02.02 [11:25:34] - Function now calls data directory via standardized ciapkg.getDirPkg('data') to avoid placing data in incorrect folder. % 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package. % 2022.11.07 [19:09:38] - Update links from Stanford to UCSF box (Stanford migrating away from Box, links are dead). + % 2022.11.14 [17:29:56] - Added additional link. % TODO % Switch Box downloads to Google Drive due to Stanford phase out. @@ -48,6 +49,7 @@ downloadList{end}.fileName = 'concat_recording_20140401_180333.h5'; % downloadList{end}.fileUrl = 'https://stanford.box.com/shared/static/jmld9o9s0oemvn6oionr3lf9lwobqk9l.h5'; downloadList{end}.fileUrl = 'https://ucsf.box.com/shared/static/dtdwpr4xjs1u67dv30qnhu7arf889ysr.h5'; + % 'http://tiny.ucsf.edu/KATTBh' if options.downloadExtraFiles==1 [downloadList] = localfxn_downloadExtra(downloadList); diff --git a/+ciapkg/+download/updatePathCleanup.m b/+ciapkg/+download/updatePathCleanup.m index 3b60132..076795a 100644 --- a/+ciapkg/+download/updatePathCleanup.m +++ b/+ciapkg/+download/updatePathCleanup.m @@ -1,7 +1,7 @@ function [status] = updatePathCleanup(inputDirectory,varargin) % Cleans up the specified directory and sub-directories from MATLAB path along with the Java path. Used to allow updating of directories. % Biafra Ahanonu - % started: 2020.06.28 [‏‎13:35:52] + % started: 2020.06.28 [13:35:52] % inputs % % outputs diff --git a/+ciapkg/+hdf5/appendDataToHdf5.m b/+ciapkg/+hdf5/appendDataToHdf5.m index 0b30f31..16a8584 100644 --- a/+ciapkg/+hdf5/appendDataToHdf5.m +++ b/+ciapkg/+hdf5/appendDataToHdf5.m @@ -17,12 +17,30 @@ function appendDataToHdf5(filename, datasetName, inputData, varargin) % 2014.07.22 % 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package. % 2022.03.13 [23:49:04] - Update comments. + % 2023.09.24 [19:29:43] - Flag to silence command line output. % TODO % + % ======================== + % Binary: 1 = whether to display info on command line. 2 = short output. + options.displayInfo = 1; + % get options + options = ciapkg.io.getOptions(options,varargin); + % disp(options) + % unpack options into current workspace + % fn=fieldnames(options); + % for i=1:length(fn) + % eval([fn{i} '=options.' fn{i} ';']); + % end + % ======================== + import ciapkg.api.* % import CIAtah functions in ciapkg package API. - display(['appending data to: ' filename]) + if options.displayInfo==1 + display(['appending data to: ' filename]) + elseif options.displayInfo==2 + fprintf('+'); + end % open the dataset and append data to the unlimited dimension which in this case is the second dimension as seen from matlab. fid = H5F.open(filename, 'H5F_ACC_RDWR', 'H5P_DEFAULT'); dset_id = H5D.open(fid, datasetName); @@ -70,5 +88,9 @@ function appendDataToHdf5(filename, datasetName, inputData, varargin) H5S.close(space_id); H5D.close(dset_id); H5F.close(fid); - display('done appending.') + if options.displayInfo==1 + display('done appending.') + elseif options.displayInfo==2 + fprintf('1|'); + end end \ No newline at end of file diff --git a/+ciapkg/+hdf5/downsampleHdf5Movie.m b/+ciapkg/+hdf5/downsampleHdf5Movie.m index 2c7dc81..32d5528 100644 --- a/+ciapkg/+hdf5/downsampleHdf5Movie.m +++ b/+ciapkg/+hdf5/downsampleHdf5Movie.m @@ -10,7 +10,7 @@ % note: checked between this implementation and imageJ's scale, they are nearly identical (subtracted histogram mean is 1 vs each image mean of ~1860, so assuming precision error). % % inputs - % inputFilePath - Str: path to movie to downsample into HDF5. + % inputFilePath - Str: path to movie to downsample into HDF5. Supported formats are HDF5, TIFF, and formats loadable by Bio-Formats. % outputs % success - Binary: 1 = completed successfully, 0 = did not complete successfully. % options @@ -22,6 +22,8 @@ % 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package. % 2022.07.10 [19:24:29] - Added Bio-Formats support. % 2022.07.15 [13:34:28] - Switch chunking to automatic by default, generally improved reading from data later than small chunks. + % 2023.09.28 [17:50:26] - Integrated TIFF support from other functions. + % 2023.09.29 [11:01:40] - options.newFilenameTwo is sufficient to initiate secondary movie downsampling. Option to convert data type during downsampling % TODO % Use handles to reduce memory load when doing computations. @@ -72,6 +74,8 @@ options.displayInfo = 1; % Binary: 1 = waitbar/progress bar is shown, 0 = no progress shown. options.waitbarOn = 1; + % Str: convert movie to new type. 'single','uint16','uint8','double' + options.convertType = ''; % get options options = getOptions(options,varargin); % unpack options into current workspace @@ -107,6 +111,8 @@ case 'hdf5' % movie dimensions and subsets to analyze [subsets, dataDim] = getSubsetOfDataToAnalyze(inputFilePath, options, varargin); + case 'tiff' + [subsets dataDim] = getSubsetOfDataToAnalyzeTIFF(inputFilePath, options, varargin); case 'bioformats' disp(['Opening connection to Bio-Formats file:' inputFilePath]) startReadTime = tic; @@ -203,27 +209,36 @@ % get current subset location and size currentSubsetLocation = subsets(currentSubset); lengthSubset = subsetsDiff(currentSubset); - % convert offset to C-style offset for low-level HDF5 functions - offset = [0 0 currentSubsetLocation-1]; - block = [dataDim.x dataDim.y lengthSubset]; display('---') - % display(sprintf(['current location: ' num2str(round(currentSubsetLocation/dataDim.z*100)) '% | ' num2str(currentSubsetLocation) '/' num2str(dataDim.z) '\noffset: ' num2str(offset) '\nblock: ' num2str(block)])); - fprintf('current location: %d%% | %d/%d \noffset: %s \nblock: %s\n',round(currentSubsetLocation/dataDim.z*100),currentSubsetLocation,dataDim.z,mat2str(offset),mat2str(block)); switch options.movieType case 'hdf5' + % convert offset to C-style offset for low-level HDF5 functions + offset = [0 0 currentSubsetLocation-1]; + block = [dataDim.x dataDim.y lengthSubset]; + % display(sprintf(['current location: ' num2str(round(currentSubsetLocation/dataDim.z*100)) '% | ' num2str(currentSubsetLocation) '/' num2str(dataDim.z) '\noffset: ' num2str(offset) '\nblock: ' num2str(block)])); + fprintf('current location: %d%% | %d/%d \noffset: %s \nblock: %s\n',round(currentSubsetLocation/dataDim.z*100),currentSubsetLocation,dataDim.z,mat2str(offset),mat2str(block)); % load subset of HDF5 file into memory inputMovie = readHDF5Subset(inputFilePath,offset,block,'datasetName',options.inputDatasetName); + case 'tiff' + thisFrameList = currentSubsetLocation:(currentSubsetLocation+lengthSubset-1); + disp(['Frames ' num2str(thisFrameList(1)) ' to ' num2str(thisFrameList(end))]) + + inputMovie = ciapkg.io.loadMovieList(inputFilePath,'frameList',thisFrameList); case 'bioformats' thisFrameList = currentSubsetLocation:(currentSubsetLocation+lengthSubset-1); disp(['Frames ' num2str(thisFrameList(1)) ' to ' num2str(thisFrameList(end))]) - [inputMovie] = local_readBioformats(inputFilePath,thisFrameList,fileIdOpen,options); + inputMovie = local_readBioformats(inputFilePath,thisFrameList,fileIdOpen,options); otherwise end + if ~isempty(options.convertType) + disp(['Converting from "' class(inputMovie) '"" to "' options.convertType '"']); + inputMovie = cast(inputMovie,options.convertType); + end % split into second movie if need be - if ~isempty(options.saveFolderTwo) + if ~isempty(options.saveFolderTwo)|~isempty(options.newFilenameTwo) inputMovieTwo = inputMovie; end @@ -249,7 +264,7 @@ toc(loopStartTime); display(sprintf(['downsample dims: ' num2str(size(inputMovie)) '\n-------'])); - if ~isempty(options.saveFolderTwo) + if ~isempty(options.saveFolderTwo)|~isempty(options.newFilenameTwo) display('secondary downsample in progress...') inputMovie = inputMovieTwo; % downsample section of the movie, keep in memory @@ -482,6 +497,22 @@ function downsampleMovieNested(varargin) % get the subsets of the 3D matrix to analyze subsets = floor(linspace(1,dataDim.z,numSubsets)); end + +function [subsets dataDim] = getSubsetOfDataToAnalyzeTIFF(inputFilePath, options, varargin) + import ciapkg.api.* % import CIAtah functions in ciapkg package API. + + dataDim = ciapkg.io.getMovieInfo(inputFilePath); + testFrame = ciapkg.io.readFrame(inputFilePath,1); + + % estimate size of movie in Mbytes + testFrameInfo = whos('testFrame'); + estSizeMovie = (testFrameInfo.bytes/options.bytesToMB)*dataDim.z; + + % get the subsets of the 3D matrix to analyze + numSubsets = ceil(estSizeMovie/options.maxChunkSize)+1; + subsets = floor(linspace(1,dataDim.z,numSubsets)); +end + function hReadInfo = getHdf5Info(hinfo,options) import ciapkg.api.* % import CIAtah functions in ciapkg package API. diff --git a/+ciapkg/+image/cropMatrix.m b/+ciapkg/+image/cropMatrix.m index 4661c8e..9ee2cbc 100644 --- a/+ciapkg/+image/cropMatrix.m +++ b/+ciapkg/+image/cropMatrix.m @@ -1,43 +1,59 @@ function [inputMatrix, coords] = cropMatrix(inputMatrix,varargin) % Crops a matrix either by removing rows or adding NaNs to where data was previously. + % Automatic version to detect boundaries after motion correction and add borders accordingly. % Biafra Ahanonu % 2014.01.23 [16:06:01] % inputs - % inputMatrix - a [m n p] matrix of any class type + % inputMatrix - a [m n p] matrix of any class type. % outputs - % inputMatrix - cropped or NaN'd matrix, same name to reduce memory usage + % inputMatrix - cropped or NaN'd matrix, same name to reduce memory usage. % changelog % 2017.01.14 [20:06:04] - support switched from [nSignals x y] to [x y nSignals]. % 2021.04.18 [14:42:53] - Updated to make imrect the default method of selecting the coordinates. % 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package. % 2021.09.28 [10:37:12] - Updated to allow specifying the size of the rectangle. + % 2023.10.01 [22:32:11] - Update to callback. + % 2024.03.11 [21:04:57] - Added support for automatic detection of crop borders from main pre-processing functions. Remove workspace unpacking of options. % TODO % import ciapkg.api.* % import CIAtah functions in ciapkg package API. %======================== + % Str: 'manual' or 'auto' + options.cropType = 'manual'; + % Str: 'NaN' to add a NaN border, 'crop' to reduce the dimensions of the input matrix. options.cropOrNaN = 'NaN'; - % input coordinates: [left-column top-row right-column bottom-row] + % Int vector: input coordinates as [left-column top-row right-column bottom-row] options.inputCoords = []; - % amount of pixels around the border to crop in primary movie + % Int: amount of pixels around the border to crop in primary movie options.pxToCrop = 0; % Str: title to add. options.title = 'Select a region'; + % Int: figure number to open for GUI cropping. options.figNo = 142568; % Vector: [xmin ymin xmax ymax] options.rectPos = []; - + % Str: 'NaN' or 'zero' for type of value used by registration function to fill space after motion correction. + options.registrationFillValue = 'NaN'; + % Str: 'identical' = all dimensions get cropped to the max movement in any one side, 'perSide' = crop only to max movement on a given side. + options.autoCropType = 'perSide'; + % % get options options = getOptions(options,varargin); % display(options) % unpack options into current workspace - fn=fieldnames(options); - for i=1:length(fn) - eval([fn{i} '=options.' fn{i} ';']); - end + % fn=fieldnames(options); + % for i=1:length(fn) + % eval([fn{i} '=options.' fn{i} ';']); + % end %======================== + if strcmp(options.cropType,'auto') + removeInputMovieEdges(); + return; + end + if options.pxToCrop>0 if size(inputMatrix,2)>=size(inputMatrix,1) coords(1) = options.pxToCrop; %xmin @@ -87,6 +103,61 @@ otherwise return end + function removeInputMovieEdges() + % turboreg outputs 0s where movement goes off the screen + thisMovieMinMask = zeros([size(inputMatrix,1) size(inputMatrix,2)]); + switch options.registrationFillValue + case 'NaN' + reverseStr = ''; + for row=1:size(inputMatrix,1) + % thisMovieMinMask(row,:) = logical(max(isnan(squeeze(inputMatrix(3,:,:))),[],2,'omitnan')); + thisMovieMinMask(row,:) = logical(max(isnan(squeeze(inputMatrix(row,:,:))),[],2,'omitnan')); + reverseStr = cmdWaitbar(row,size(inputMatrix,1),reverseStr,'inputStr','getting crop amount','waitbarOn',1,'displayEvery',5); + end + case 'zero' + reverseStr = ''; + for row=1:size(inputMatrix,1) + thisMovieMinMask(row,:) = logical(min(squeeze(inputMatrix(row,:,:))~=0,[],2,'omitnan')==0); + reverseStr = cmdWaitbar(row,size(inputMatrix,1),reverseStr,'inputStr','getting crop amount','waitbarOn',1,'displayEvery',5); + end + otherwise + % do nothing + end + topVal = sum(thisMovieMinMask(1:floor(end/4),floor(end/2))); + bottomVal = sum(thisMovieMinMask(end-floor(end/4):end,floor(end/2))); + leftVal = sum(thisMovieMinMask(floor(end/2),1:floor(end/4))); + rightVal = sum(thisMovieMinMask(floor(end/2),end-floor(end/4):end)); + tmpPxToCrop = max([topVal bottomVal leftVal rightVal]); + pxToCropPreprocess = tmpPxToCrop; + display(['[topVal bottomVal leftVal rightVal]: ' num2str([topVal bottomVal leftVal rightVal])]) + + % Get the crop regions. + switch options.autoCropType + case 'perSide' + topRowCrop = topVal; % top row + leftColCrop = leftVal; % left column + bottomRowCrop = size(inputMatrix,1)-bottomVal; % bottom row + rightColCrop = size(inputMatrix,2)-rightVal; % right column + case 'identical' + topRowCrop = pxToCropPreprocess; % top row + leftColCrop = pxToCropPreprocess; % left column + bottomRowCrop = size(inputMatrix,1)-pxToCropPreprocess; % bottom row + rightColCrop = size(inputMatrix,2)-pxToCropPreprocess; % right column + otherwise + % do nothing + end + + % rowLen = size(inputMatrix,1); + % colLen = size(inputMatrix,2); + % set leftmost columns to NaN + inputMatrix(1:end,1:leftColCrop,:) = NaN; + % set rightmost columns to NaN + inputMatrix(1:end,rightColCrop:end,:) = NaN; + % set top rows to NaN + inputMatrix(1:topRowCrop,1:end,:) = NaN; + % set bottom rows to NaN + inputMatrix(bottomRowCrop:end,1:end,:) = NaN; + end end function [coords] = getCropCoords(thisFrame,options) import ciapkg.api.* % import CIAtah functions in ciapkg package API. @@ -133,9 +204,9 @@ else % [xmin ymin width height] pos = options.rectPos; - h= imrect(gca,[pos(1) pos(2) pos(3)-pos(1) pos(4)-pos(2)]); + h = imrect(gca,[pos(1) pos(2) pos(3)-pos(1) pos(4)-pos(2)]); end - addNewPositionCallback(h,@(p) title([titleStr 10 mat2str(p,3)])); + addNewPositionCallback(h,@(p) title([titleStr 10 mat2str(p,5) ' | ' num2str([p(3)-p(1) p(4)-p(2)])])); fcn = makeConstrainToRectFcn('imrect',get(gca,'XLim'),get(gca,'YLim')); setPositionConstraintFcn(h,fcn); end \ No newline at end of file diff --git a/+ciapkg/+image/thresholdImages.m b/+ciapkg/+image/thresholdImages.m index 731237b..9b4309b 100644 --- a/+ciapkg/+image/thresholdImages.m +++ b/+ciapkg/+image/thresholdImages.m @@ -25,6 +25,7 @@ % 2021.07.06 [18:57:02] - Fast border calculation using convn for options.fastThresholding==1, less precise than bwboundaries or bwareafilt but works for fast display purposes. % 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package. % 2022.06.27 [19:41:34] - manageParallelWorkers now passed options.waitbarOn value to reduce command line clutter if user request in thresholdImages. + % 2023.02.22 [19:47:27] - fastThresholding respects normalization=0, so threshold functions as absolute value. Also make sure remove unconnected works with fastThresholding at all times. % TODO % @@ -132,10 +133,18 @@ if options.fastThresholding==1 if options_binary==1 - inputImages = inputImages>max(inputImages,[],[1 2])*options_threshold; + if options_normalize==1 + inputImages = inputImages>max(inputImages,[],[1 2])*options_threshold; + else + inputImages = inputImages>options_threshold; + end imgFiltHere = options_imageFilterBinary; else - rmIdx = inputImages0.99),y11(conf11>0.99),'.'); + end + end + + % ======================== + opts.rotateImg = options.rotateImg; + % opts.rotateImg = 1; + opts.flagMakeColorVid = 0; + opts.refIdx = options.refFrame; + + imgRef = inputMovie(:,:,opts.refIdx); + if opts.rotateImg==1 + imgRef = rot90(imgRef,-1); + end + + opts.nFramesHere = size(inputMovie,3); + opts.nFrames = size(inputMovie,3); + if isempty(options.frameList) + opts.frameList = 1:opts.nFrames; + else + opts.frameList = options.frameList; + end + + if opts.rotateImg==1 + mtVideo = zeros([size(inputMovie,2) size(inputMovie,1) length(opts.frameList)],class(inputMovie)); + else + mtVideo = zeros([size(inputMovie,1) size(inputMovie,2) length(opts.frameList)],class(inputMovie)); + end + + if opts.flagMakeColorVid==1 + colorVideo = zeros([size(imgRef,1) size(imgRef,2) 3 opts.nFramesHere],'uint8'); + end + + disp('Done!') + + % ======================== + %% Run control point tracking + opts.confThreshold = options.confThreshold; + opts.correctMotion = options.correctMotion; + opts.dispPlots = options.makePlots; + frameCount = 1; + opts.filtPts = options.filtPts; + opts.maxDist = options.maxDist; + + % Get the reference control points and associated confidence + if opts.rotateImg==1 + pointsRef = ctrlPtTensor(:,[3 2],opts.refIdx); + confRef = ctrlPtTensor(:,1,opts.refIdx); + else + pointsRef = ctrlPtTensor(:,[2 3],opts.refIdx); + confRef = ctrlPtTensor(:,1,opts.refIdx); + end + % pointsRef = ctrlPtTensor(:,[2 3],opts.refIdx); + % confRef = ctrlPtTensor(:,1,opts.refIdx); + % pointsRef = featurePts{opts.refIdx}(:,[2 3]); + % confRef = featurePts{opts.refIdx}(:,1); + + % Run on all frames requested by the user + for fNo = opts.frameList + if mod(fNo,100)==0 + disp(['===' 10]) + end + disp([num2str(fNo) ' | ']) + if opts.rotateImg==1 + pointsReg = ctrlPtTensor(:,[3 2],fNo); + confReg = ctrlPtTensor(:,1,fNo); + else + pointsReg = ctrlPtTensor(:,[2 3],fNo); + confReg = ctrlPtTensor(:,1,fNo); + end + % pointsReg = featurePts{fNo}(:,[2 3]); + % confReg = featurePts{fNo}(:,1); + + % Filter only features with reference and moving likelihood above threshold to avoid bad registration. + confIdx = confReg>=opts.confThreshold & confRef>=opts.confThreshold; + pointsRegFilt = pointsReg(confIdx,:); + pointsRefFilt = pointsRef(confIdx,:); + + % Scale the registration points if the input movie has a different size compared to control point tracking. + if options.dsSpace~=1 + pointsRegFilt = pointsRegFilt*options.dsSpace; + pointsRefFilt = pointsRefFilt*options.dsSpace; + end + + % Use normalized movie + imgReg = inputMovieNorm(:,:,fNo); + + % If not enough registration points, skip motion correction, likely blurred movie + if sum(confIdx)=opts.confThreshold&confRef(filtPts)>=opts.confThreshold; + pointsRegFilt = pointsReg(confIdx,:); + pointsRefFilt = pointsRefTmpTmp(confIdx,:); + + end + % If requested, fix a particular dimension + if ~isempty(options.dimToFix) + pointsRegFilt(:,options.dimToFix) = pointsRefFilt(:,options.dimToFix); + end + + if opts.rotateImg==1 + % imgReg = rot90(imgReg,-1); + imgReg = permute(imgReg,[2 1]); + + end + + if opts.dispPlots==1 + % clf + subplotTmp(2,2,1) + localFxn_showMatchedFeaturesFast(imgRef,imgReg,pointsRefFilt,pointsRegFilt) + box off; + title('Overlap reference and raw') + subplotTmp(2,2,3) + imagesc(imgReg) + box off + axis image + title('Raw image') + colormap(gca,'gray') + end + + if opts.correctMotion==1 + % "rigid", "similarity", "affine", or "projective" + opts.cpRegMethod = options.cpRegMethod; + [imgReg,pointsReg,pointsRefPlot] = localFxn_motionCorrectPtFeatures(imgRef,imgReg,pointsRefFilt,pointsRegFilt,opts.cpRegMethod,options.maxTrials,opts.maxDist); + else + pointsRefPlot = pointsRefFilt; + end + + if opts.dispPlots==1 + subplotTmp(2,2,2) + localFxn_showMatchedFeaturesFast(imgRef,imgReg,pointsRefPlot,pointsReg) + box off; + title('Overlap reference and registered') + %showMatchedFeatures(rot90(imgRef,-1),rot90(imgReg,-1),pointsRefPlot,pointsReg) + %imagesc(imgReg);axis off;box off;axis image;colormap gray; + title(fNo) + + subplotTmp(2,2,4) + imagesc(imgReg) + box off + axis image + title('Registered image') + colormap(gca,'gray') + drawnow + pause(0.01); + end + + mtVideo(:,:,frameCount) = imgReg; + frameCount = frameCount + 1; + + if opts.flagMakeColorVid==1 + g2 = getframe; + g2 = imresize(g2.cdata,[size(imgRef,1) size(imgRef,2)]); + colorVideo(:,:,:,fNo) = g2; + end + end + + if options.runPostReg==1 + % ======================== + %% Get TurboReg registration coordinates + if isempty(options.cropCoordsInSession) + [~,cropCoordsInSession] = ciapkg.image.cropMatrix(single(mtVideo(:,:,opts.refIdx)),'cropOrNaN','crop','inputCoords',[]); + else + cropCoordsInSession = options.cropCoordsInSession; + end + + % ======================== + %% Turboreg to correct affine motion after + % Motion correction + toptions.cropCoords = cropCoordsInSession; + toptions.turboregRotation = 0; + toptions.removeEdges = 1; + toptions.pxToCrop = 10; + toptions.RegisType = 1; + toptions.removeNan = 1; + toptions.refFrame = opts.refIdx; + toptions.refFrameMatrix = mtVideo(:,:,opts.refIdx); + toptions.registrationFxn = 'imwarp'; % transfturboreg has issues with NaNs in displacement field output. + % Pre-motion correction + toptions.complementMatrix = 1; + toptions.meanSubtract = 1; + toptions.meanSubtractNormalize = 1; + toptions.normalizeType = 'matlabDisk'; % matlabDisk + % Spatial filter + toptions.normalizeBeforeRegister = []; %'divideByLowpass' + toptions.freqLow = 0; + toptions.freqHigh = 4; + disp('Done!') + + %% Run TurboReg to correct fine remaining motion + mtVideo = ciapkg.motion_correction.turboregMovie(mtVideo,... + 'options',toptions); + disp('Done!') + end + + catch err + disp(repmat('@',1,7)) + disp(getReport(err,'extended','hyperlinks','on')); + disp(repmat('@',1,7)) + end + + % function [outputs] = nestedfxn_exampleFxn(arg) + % % Always start nested functions with "nestedfxn_" prefix. + % % outputs = ; + % end +end +%% HELPER FUNCTIONS +function [imgBp,pointsBmp,pointsAm] = localFxn_motionCorrectPtFeatures(imgA,imgB,pointsA,pointsB,transformType,MaxNumTrials,MaxDistance) + [tform,inlierIdx] = estgeotform2d(pointsB,pointsA,transformType,'MaxNumTrials',MaxNumTrials,'MaxDistance',MaxDistance); + pointsBm = pointsB(inlierIdx,:); + pointsAm = pointsA(inlierIdx,:); + % [tform,pointsBm,pointsAm] = estimateGeometricTransform(featurePts{2},featurePts{1},'similarity','Confidence',50); + imgBp = imwarp(imgB,tform,'OutputView',imref2d(size(imgB)),'FillValues',NaN('single')); + pointsBmp = transformPointsForward(tform,pointsBm); +end +function [] = localFxn_getMovieShift(imgA,imgB,pointsA,pointsB) + localFxn_showMatchedFeaturesFast(imgA,imgBp,pointsAm,pointsBmp) +end +function localFxn_showMatchedFeaturesFast(image1,image2,matchedPoints1,matchedPoints2) + % Modified version of showMatchedFeatures for fast plotting of comparisons. + + image1 = single(image1); + image2 = single(image2); + rgbImg = zeros([size(image1,1) size(image1,2) 3]); + normImgFun = @(x) (x-min(x(:),[],'omitnan'))/(max(x(:),[],'omitnan')-min(x(:),[],'omitnan')); + rgbImg(:,:,1) = normImgFun(image1); + rgbImg(:,:,2) = normImgFun(image2); + rgbImg(:,:,3) = normImgFun(image2); + + imagesc(rgbImg) + axis image + hold on; + + plot(matchedPoints1(:,1), matchedPoints1(:,2),'o','Color','r'); + plot(matchedPoints2(:,1), matchedPoints2(:,2),'+','Color','g'); + + % Plot by using a single line object with line segments broken by using NaNs. This is more efficient and makes it easier to customize the lines. + lineX = [matchedPoints1(:,1)'; matchedPoints2(:,1)']; + numPts = numel(lineX); + lineX = [lineX; NaN(1,numPts/2)]; + + lineY = [matchedPoints1(:,2)'; matchedPoints2(:,2)']; + lineY = [lineY; NaN(1,numPts/2)]; + + plot(lineX(:), lineY(:), 'y-'); % line + hold off; +end \ No newline at end of file diff --git a/+ciapkg/+motion_correction/turboregMovie.m b/+ciapkg/+motion_correction/turboregMovie.m index 1f5a86b..cb2d727 100644 --- a/+ciapkg/+motion_correction/turboregMovie.m +++ b/+ciapkg/+motion_correction/turboregMovie.m @@ -1,19 +1,19 @@ function [inputMovie, ResultsOutOriginal] = turboregMovie(inputMovie, varargin) % [inputMovie, ResultsOutOriginal] = turboregMovie(inputMovie, varargin) - % - % Motion corrects (using turboreg or NoRMCorre) a movie. - % - Both turboreg (to get 2D translation coordinates) and registering images (transfturboreg, imwarp, imtransform) have been parallelized. - % - Can also turboreg to one set of images and apply the registration to another set (e.g. for cross-day alignment). + % + % Motion corrects (using turboreg, NoRMCorre, displacement fields, etc.) a movie. + % - Both turboreg (to get 2D translation coordinates) and registering images (transfturboreg, imwarp, imtransform) have been parallelized. + % - Can also turboreg to one set of images and apply the registration to another set (e.g. for cross-day alignment). % - Spatial filtering is applied after obtaining registration coordinates but before transformation, this reduced chance that 0s or NaNs at edge after transformation mess with proper spatial filtering. - % + % % Biafra Ahanonu % started 2013.11.09 [11:04:18] - % + % % Parts of code based on that from Jerome Lecoq (2011) and parallel code update by Biafra Ahanonu (2013). - % + % % Input % inputMovie - 3D matrix: [x y frames] matrix containing data to be motion corrected across frames. - % + % % Output % inputMovie - 3D matrix: [x y frames] matrix containing data after motion correction across frames. % ResultsOutOriginal @@ -46,6 +46,14 @@ % 2022.03.09 [17:32:42] - Use custom bahanonu NoRMCorre that is within a package and update to use reference picture instead of template. % 2022.04.23 [18:40:22] - Updated to which('normcorre.normcorre') from which('normcorre') since normcorre now inside a package. % 2022.09.12 [20:46:04] - Additional NoRMCorre support and getNoRMCorreParams checking before running NoRMCorre. + % 2022.09.24 [18:32:04] - Changes to GUI display. + % 2022.10.03 [18:20:34] - Correct user options input to normcorre.apply_shifts. + % 2022.11.15 [14:16:40] - Add displacement field-based motion correction. + % 2022.11.22 [15:34:27] - Add back in ability to remove NaNs from each frame before calculating shifts or displacements. + % 2022.12.05 [18:35:05] - Add support for gradient-based filtering for motion correction input. + % 2023.08.04 [07:34:45] - Support loading of prior coordinates from MAT-file in standard of modelPreprocessMovieFunction. + % 2024.02.18 [18:46:34] - Added additional options for displacement field-based motion correction. + % 2024.03.11 [20:53:19] - Update to row assignment for removeInputMovieEdges. % TO-DO % Add support for "imregtform" based registration. @@ -63,7 +71,6 @@ % ======================== % Check that input is not empty - disp('Starting motion correction (turboreg)...'); if isempty(inputMovie) disp('Empty movie matrix, exiting motion correction...') return; @@ -75,6 +82,7 @@ % Str: motion correction algorithm. % 'turboreg' - TurboReg as developed by Philippe Thevenaz. % 'normcorre' - NoRMCorre as developed by several authors at Flatiron Institute. + % 'imregdemons' - Displacement field alignment based on imregdemons. options.mcMethod = 'turboreg'; % options.mcMethod = 'normcorre'; % Str: dataset name in HDF5 file where data is stored, if inputMovie is a path to a movie. @@ -85,8 +93,12 @@ options.refFrame = 1; % Matrix: same type as the inputMovie, this will be appended to the end of the movie and used as the reference frame. This is for cases in which the reference frame is not contained in the movie. options.refFrameMatrix = []; - % whether to use 'imtransform' or 'imwarp' (Matlab) or 'transfturboreg' (C) + % Str: Register images using 'imtransform' or 'imwarp' (Matlab) or 'transfturboreg' (C) options.registrationFxn = 'transfturboreg'; + % Binary: 1 = remove NaNs before computing transform/displacement matrix, 0 = do not alter NaN state. + options.removeNan = 0; + % Float: any numeric value to replace NaNs with if options.removeNan==1. + options.nanReplaceVal = 0; % display options for verification of input options.displayOptions = 0; % character string, path to save turboreg coordinates @@ -188,7 +200,7 @@ % ======= % 1 = show figures during processing options.showFigs = 1; - % cmd line waitbar on? + % Binary: 1 = cmd line waitbar on. 0 = waitbar off. options.waitbarOn = 1; % ======= % Binary: 1 = return the movie after normalizing, mean subtract, etc. @@ -199,6 +211,19 @@ % NoRMCorre % Struct: NoRMCorre settings options.optsNoRMCorre = ciapkg.motion_correction.getNoRMCorreParams([1 1 1],'guiDisplay',0); + % ======= + % Displacement fields options (imregdemons) + % Int: Indicate which displacement fields to remove. + % [] = include all displacement fields + % 1 = exclude Y displacement fields + % 2 = exclude X displacement fields + options.dfExclude = []; + % Options for displacement fields + options.df_AccumulatedFieldSmoothing = 0.54; + options.df_Niter = [500 400 200]; + options.df_PyramidLevels = 3; + options.df_DisplayWaitbar = false; + % ======= % % get options options = getOptions(options,varargin); @@ -211,6 +236,9 @@ if options.displayOptions==1 fn_structdisp(options) end + % ======================== + + fprintf('Starting motion correction (%s)...\n',options.mcMethod); % ======================== % Verify that turboreg MEX function is in the path. @@ -227,6 +255,9 @@ case 'normcorre' % NoRMCorre is to be run on the entire input matrix. options.cropCoords = []; + case 'imregdemons' + % NoRMCorre is to be run on the entire input matrix. + options.cropCoords = []; otherwise end @@ -266,6 +297,7 @@ options.maxFrame = size(inputMovie,3); movieDim = size(inputMovie); options.Levels=nestFxnCalculatePyramidDepth(min(movieDim(1),movieDim(2))); + % inputMovie(isnan(inputMovie)) = 0; % ======================== % add turboreg options to turboRegOptions structure @@ -279,18 +311,18 @@ turboRegOptions.zapMean = options.zapMean; turboRegOptions.Interp = options.Interp; if any(turboRegOptions.RegisType==[1 2 4]) - TransformationType='affine'; + TransformationType = 'affine'; else - % 3 5 - TransformationType='projective'; + % RegisType = [3 5] + TransformationType = 'projective'; end % ======================== - % register movie and return without using the rest of the function + % Register movie and return without using the rest of the function if ~isempty(options.precomputedRegistrationCooords) disp('Input pre-computed registration coordinates...') ResultsOut = options.precomputedRegistrationCooords; ResultsOutOriginal = ResultsOut; - for resultNo=1:size(inputMovie,3) + for resultNo = 1:size(inputMovie,3) ResultsOutTemp{resultNo} = ResultsOut{options.altMovieRegisterNum}; end ResultsOut = ResultsOutTemp; @@ -310,10 +342,19 @@ ResultsOutOriginal = ResultsOut; return; end + % ======================== - % register movie and return without using the rest of the function + % Register movie and return without using the rest of the function if ~isempty(options.precomputedRegistrationCooordsFullMovie) disp('Input pre-computed registration coordinates for full movie...') + + % Check if MAT-file is input, load into memory + if ischar(options.precomputedRegistrationCooordsFullMovie) + disp(['Loading: ' options.precomputedRegistrationCooordsFullMovie]) + loadTmp = load(options.precomputedRegistrationCooordsFullMovie); + options.precomputedRegistrationCooordsFullMovie = loadTmp.ResultsOutOriginal{1}{1}; + end + ResultsOut = options.precomputedRegistrationCooordsFullMovie; % ResultsOutOriginal = ResultsOut; % for resultNo=1:size(inputMovie,3) @@ -459,9 +500,9 @@ % inputMovie = cat(3,inputMovie{:}); inputMovie = single(inputMovie); - - subfxn_dispMovieFrames(inputMovie,'Registration==1',2); - + if options.showFigs==1 + subfxn_dispMovieFrames(inputMovie,'Registration==1',2); + end % ======================== if options.removeEdges==1 removeInputMovieEdges(); @@ -551,6 +592,8 @@ function turboregMovieParallel() startTurboRegTime = tic; % nFramesToTurboreg = options.maxFrame; + options_removeNan = options.removeNan; + options_nanReplaceVal = options.nanReplaceVal; if options.parallel==1; nWorkers=Inf;else;nWorkers=0;end parfor (frameNo=1:nFramesToTurboreg,nWorkers) % get current frames @@ -561,6 +604,11 @@ function turboregMovieParallel() thisFrameToAlign=single(thisFrame); % thisFrameToAlign=thisFrame; + if options_removeNan==1 + thisFrame(isnan(thisFrame)) = options_nanReplaceVal; + thisFrameToAlign(isnan(thisFrameToAlign)) = options_nanReplaceVal; + end + if ismac % Code to run on Mac platform [ImageOut,ResultsOut{frameNo}] = turboreg(refPic,thisFrameToAlign,mask,imgRegMask,turboRegOptions); @@ -594,6 +642,7 @@ function turboregMovieParallel() bound = 0; refPic = single(squeeze(inputMovieCropped(:,:,options.refFrame))); + % refPic = refPic(bound/2+1:end-bound/2,bound/2+1:end-bound/2,:); inputMovieCroppedTmp = inputMovieCropped(bound/2+1:end-bound/2,bound/2+1:end-bound/2,:); optsNoRMCorre = options.optsNoRMCorre; @@ -605,14 +654,51 @@ function turboregMovieParallel() % Ensure NoRMCorre gets the correct inputs optsNoRMCorre.d1 = size(inputMovieCropped,1); optsNoRMCorre.d2 = size(inputMovieCropped,2); + disp('===') + fn_structdisp(optsNoRMCorre) + optsNoRMCorre.method % fn_structdisp(optsNoRMCorre) - % [optsNoRMCorre] = subfxn_getNoRMCorreParams(inputMovieCropped); + % [optsNoRMCorre2] = subfxn_getNoRMCorreParams(inputMovieCropped); + % disp('===') + % fn_structdisp(optsNoRMCorre2) + % optsNoRMCorre2.method + disp('===') [~,ResultsOut,template2] = normcorre.normcorre_batch(... inputMovieCroppedTmp,... optsNoRMCorre,... refPic); clear inputMovieCroppedTmp - toc(startTurboRegTime); + toc(startTurboRegTime); + case 'imregdemons' + disp('Demon-based displacement field registration...') + nFramesToTurboreg = options.maxFrame; + + fixed = single(squeeze(inputMovieCropped(:,:,options.refFrame))); + if options.parallel==1; nWorkers=Inf;else;nWorkers=0;end + + df_AccumulatedFieldSmoothing = options.df_AccumulatedFieldSmoothing; + df_Niter = options.df_Niter; + df_PyramidLevels = options.df_PyramidLevels; + df_DisplayWaitbar = options.df_DisplayWaitbar; + + parfor (frameNo=1:nFramesToTurboreg,nWorkers) + thisFrame = inputMovieCropped(:,:,frameNo); + moving = single(thisFrame); + % Niter = df_Niter; + % AccumulatedFieldSmoothing = df_AccumulatedFieldSmoothing; + + [ResultsOut{frameNo},movingReg] = imregdemons(moving,fixed,... + df_Niter,... + 'PyramidLevels',df_PyramidLevels,... + 'DisplayWaitbar',df_DisplayWaitbar,... + 'AccumulatedFieldSmoothing',df_AccumulatedFieldSmoothing); + + if ~verLessThan('matlab', '9.2') + send(D, frameNo); % Update + end + end + % [Dxy,movingReg] = imregdemons(moving,fixed,Niter,... + % 'AccumulatedFieldSmoothing',1.3,'DisplayWaitbar',false); otherwise end @@ -621,7 +707,7 @@ function turboregMovieParallel() function [optsNC] = subfxn_getNoRMCorreParams(inputMovieHere) [d1,d2,T] = size(inputMovieHere); bound = 0; - sizeCorFactor = 2; + sizeCorFactor = 1; % sizeCorFactor = 1; optsNC = normcorre.NoRMCorreSetParms(... 'd1', d1-bound,... @@ -649,12 +735,90 @@ function registerMovie() disp(''); switch options.mcMethod + case 'imregdemons' + subsetSize = options.subsetSizeFrames; + nFramesHere = size(inputMovie,3); + % numSubsets = ceil(length(inputMovie)/subsetSize)+1; + numSubsets = ceil(nFramesHere/subsetSize)+1; + % subsetList = round(linspace(1,length(inputMovie),numSubsets)); + subsetList = round(linspace(1,nFramesHere,numSubsets)); + display(['registering sublists: ' num2str(subsetList)]); + if options.turboregRotation==1 + disp('Using rotation in registration') + end + fprintf('Performing %s registration and %s transformation.\n',options.registrationFxn,TransformationType); + % ResultsOut{1}.Rotation + nSubsets = (length(subsetList)-1); + for thisSet=1:nSubsets + subsetStartIdx = subsetList(thisSet); + subsetEndIdx = subsetList(thisSet+1); + if thisSet==nSubsets + movieSubset = subsetStartIdx:subsetEndIdx; + display([num2str(subsetStartIdx) '-' num2str(subsetEndIdx) ' ' num2str(thisSet) '/' num2str(nSubsets)]) + else + movieSubset = subsetStartIdx:(subsetEndIdx-1); + display([num2str(subsetStartIdx) '-' num2str(subsetEndIdx-1) ' ' num2str(thisSet) '/' num2str(nSubsets)]) + end + + % Get a slice of the data. + % movieDataTemp(movieSubset) = inputMovie(movieSubset); + + % movieDataTemp = inputMovie(:,:,movieSubset); + + % loop over and register each frame + if options.parallel==1; nWorkers=Inf;else;nWorkers=0;end + + nMovieSubsets = length(movieSubset); + + parfor (i = movieSubset,nWorkers) + rOutTmp = ResultsOut{i}; + + thisFrameT = inputMovie(:,:,i); + + D_restrict = rOutTmp; + if ~isempty(options.dfExclude) + D_restrict(:,:,options.dfExclude) = 0; + end + + % Warp frame based on estimated displacement field + thisFrameT = imwarp(thisFrameT,D_restrict,'linear','FillValues',NaN); + + inputMovie(:,:,i) = thisFrameT; + end + dispstat('Finished.','keepprev'); + end + case 'normcorre' startTurboRegTime = tic; - [optsNC] = subfxn_getNoRMCorreParams(inputMovie); + % [optsNC] = subfxn_getNoRMCorreParams(inputMovie); bound = 0; - inputMovie = normcorre.apply_shifts(inputMovie,ResultsOut,optsNC,bound/2,bound/2); % apply the shifts to the removed percentile + + optsNoRMCorre = options.optsNoRMCorre; + if isempty(fieldnames(optsNoRMCorre)) + disp('Getting NoRMCorre params...') + optsNoRMCorre = ciapkg.motion_correction.getNoRMCorreParams(size(inputMovie),'guiDisplay',0); + end + + % Ensure NoRMCorre gets the correct inputs + optsNoRMCorre.d1 = size(inputMovie,1); + optsNoRMCorre.d2 = size(inputMovie,2); + disp('===') + fn_structdisp(optsNoRMCorre) + optsNoRMCorre.method + % fn_structdisp(optsNoRMCorre) + % [optsNoRMCorre2] = subfxn_getNoRMCorreParams(inputMovieCropped); + % disp('===') + % fn_structdisp(optsNoRMCorre2) + % optsNoRMCorre2.method + disp('===') + + % apply the shifts to the removed percentile + inputMovie = normcorre.apply_shifts(... + inputMovie,... + ResultsOut,... + optsNoRMCorre,bound/2,bound/2); toc(startTurboRegTime); + case 'turboreg' % Need to register subsets of the movie so parfor won't crash due to serialization errors. % TODO: make this subset based on the size of the movie, e.g. only send 1GB chunks to workers. @@ -797,7 +961,8 @@ function removeInputMovieEdges() case 'imtransform' reverseStr = ''; for row=1:size(inputMovie,1) - thisMovieMinMask(row,:) = logical(max(isnan(squeeze(inputMovie(3,:,:))),[],2,'omitnan')); + % thisMovieMinMask(row,:) = logical(max(isnan(squeeze(inputMovie(3,:,:))),[],2,'omitnan')); + thisMovieMinMask(row,:) = logical(max(isnan(squeeze(inputMovie(row,:,:))),[],2,'omitnan')); reverseStr = cmdWaitbar(row,size(inputMovie,1),reverseStr,'inputStr','getting crop amount','waitbarOn',1,'displayEvery',5); end case 'transfturboreg' @@ -1045,6 +1210,21 @@ function cropAndNormalizeInputMovie() switch options.normalizeType case 'imagejFFT' imagefFftOnInputMovie('inputMovieCropped'); + case 'gradient' + disp('Creating gradient image') + nImages = size(inputMovieCropped,3); + if options.parallel==1; nWorkers=Inf;else;nWorkers=0;end + parfor (imageNo = 1:nImages,nWorkers) + % imageNow = squeeze(inputMovieCropped{imageNo}); + % inputMovieCropped{imageNo} = transform(imageNow); + + imageNow = squeeze(inputMovieCropped(:,:,imageNo)); + inputMovieCropped(:,:,imageNo) = gradient(imageNow); + + % if (mod(imageNo,20)==0|imageNo==nImages) + % reverseStr = cmdWaitbar(imageNo,nImages,reverseStr,'inputStr','fspecial normalizing'); + % end + end case 'matlabDisk' % single(inputMovieCropped) disp('Matlab fspecial disk background removal') @@ -1132,21 +1312,25 @@ function cropAndNormalizeInputMovie() disp('=======') end function subfxn_dispMovieFrames(inputMovieCropped,titleStr,inputMod) - ciapkg.api.openFigure(9019+inputMod, ''); + ciapkg.api.openFigure(9019, ''); colormap gray; - subplot(2,2,1) + rowN = 2; + colN = 3; + spOffset = (inputMod-1)*colN; + subplot(rowN,colN,1+spOffset) imagesc(squeeze(inputMovieCropped(:,:,1))); axis image; box off; title('Images to use to get registration coordinates'); - subplot(2,2,2) + ylabel(titleStr) + subplot(rowN,colN,2+spOffset) imagesc(squeeze(inputMovieCropped(:,:,end))); axis image; box off; title('Last frame in movie') - subplot(2,2,3) + subplot(rowN,colN,3+spOffset) imagesc(squeeze(inputMovieCropped(:,:,1))-squeeze(inputMovieCropped(:,:,end))); axis image; box off; title('Diff image #1 and #2') - ciapkg.overloaded.suptitle(titleStr) + % ciapkg.overloaded.suptitle(titleStr) drawnow end function cropCoords = getCropSelection(thisFrame) diff --git a/+ciapkg/+movie_processing/detrendMovie.m b/+ciapkg/+movie_processing/detrendMovie.m index 7e1ab8c..1a7d8ed 100644 --- a/+ciapkg/+movie_processing/detrendMovie.m +++ b/+ciapkg/+movie_processing/detrendMovie.m @@ -8,13 +8,13 @@ % % changelog - % + % 2023.04.04 [19:36:21] - Fix detrendDegree passing. % TODO % % ======================== % Int: 1 = linear detrend, >1 = nth-degree polynomial detrend - option.detrendDegree = 1; + options.detrendDegree = 1; % Int: maximum frame to normalize. options.maxFrame = size(inputMovie,3); % get options @@ -30,7 +30,7 @@ try inputMovie = ciapkg.movie_processing.normalizeMovie(inputMovie,... 'normalizationType','detrend',... - 'detrendDegree',option.detrendDegree,... + 'detrendDegree',options.detrendDegree,... 'maxFrame',options.maxFrame); catch err disp(repmat('@',1,7)) diff --git a/+ciapkg/+movie_processing/dfofMovie.m b/+ciapkg/+movie_processing/dfofMovie.m index 002684f..80170d5 100644 --- a/+ciapkg/+movie_processing/dfofMovie.m +++ b/+ciapkg/+movie_processing/dfofMovie.m @@ -10,6 +10,7 @@ % 2013.11.22 [17:49:34] % 2021.06.29 [12:11:43] - Add support for minimum and soft minimum dF/F calculation. % 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package. + % 2024.04.3 [12:40:35] - Updated to allow users to adjust movies with negative values. % TODO % @@ -26,6 +27,8 @@ options.minSoftPct = 0.1/100; % Binary: 1 = waitbar on options.waitbarOn = 1; + % Float: Number >1.0 indicating how much to adjust + options.minAdjust = 1.1; % get options options = getOptions(options,varargin); % display(options) @@ -52,6 +55,14 @@ end %======================== + % Adjust for problems with movies that have negative pixel values before dfof. + if isempty(options.minAdjust) + minMovie = min(inputMovie,[],[1 2 3],'omitnan'); + if minMovie<0 + inputMovie = inputMovie + options.minAdjust*abs(minMovie); + end + end + stdList = {'slidingZscore','binnedZscore','dfstd'}; inputMovieClass = class(inputMovie); diff --git a/+ciapkg/+movie_processing/normalizeMovie.m b/+ciapkg/+movie_processing/normalizeMovie.m index fa6c26b..181e63b 100644 --- a/+ciapkg/+movie_processing/normalizeMovie.m +++ b/+ciapkg/+movie_processing/normalizeMovie.m @@ -16,6 +16,12 @@ % 2021.06.02 [20:04:49] - Detrend movie now uses nanmean to get around issues in motion corrected videos. % 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package. % 2022.05.16 [16:39:45] - Switch away from using nanmean to mean('omitnan'). Detrend flight code refactor. + % 2022.11.15 [15:14:14] - matlabdisk no longer converts to cell array, bad memory usage. Improve other memory usage in function. + % 2022.12.05 [14:15:48] - Change how NaN check before filtering is run. Remove on a per-frame basis so that the NaNs can be added back in at the correct pixel locations. + % 2022.12.05 [14:22:03] - Further eliminate nanmean/nanmin/nanmax usage, use dims instead of (:) [waste memory], and general refactoring. + % 2022.12.06 [17:19:16] - Added movie subset support to reduce memory usage with parfor. + % 2023.04.04 [19:36:21] - Fix detrendDegree passing. + % 2023.06.08 [09:17:06] - Added subfxnFft movie subset support to reduce memory usage with parfor. % TODO % @@ -39,6 +45,8 @@ options.maxFrame = size(inputMovie,3); % Binary: 1 = use parallel registration (using matlab pool). options.parallel = 1; + % Int: number of frames to use for subsetting when using parfor, to reduce memory usage. + options.subsetSize = 300; % === % options for fft % Int: for bandpass, low freq to pass @@ -83,7 +91,7 @@ % Binary: 1 = yes, normalize Matlab FFT test movies options.normalizeMatlabFftTestMovies = 0; % Int: 1 = linear detrend, >1 = nth-degree polynomial detrend - option.detrendDegree = 1; + options.detrendDegree = 1; % === % Str: hierarchy name in hdf5 where movie data is located options.inputDatasetName = '/1'; @@ -163,19 +171,19 @@ subfxnDetrend(); case 'meanSubtraction' disp('mean subtracting movie'); - inputMean = nanmean(nanmean(inputMovie,1),2); + inputMean = mean(inputMovie,[1 2],'omitnan'); inputMean = cast(inputMean,class(inputMovie)); inputMovie = bsxfun(@minus,inputMovie,inputMean); case 'meanDivision' disp('mean dividing movie'); - inputMean = nanmean(nanmean(inputMovie,1),2); + inputMean = mean(inputMovie,[1 2],'omitnan'); inputMean = cast(inputMean,class(inputMovie)); inputMovie = bsxfun(@rdivide,inputMovie,inputMean); % inputMean = nansum(nansum(inputMovie,1),2); % inputMean = cast(inputMean,class(inputMovie)) case 'negativeRemoval' disp('Switching movie negative values to positive (abs)'); - inputMin = abs(nanmin(inputMovie(:))); + inputMin = abs(min(inputMovie,[],[1 2 3],'omitnan')); inputMin = cast(inputMin,class(inputMovie)); inputMovie = bsxfun(@plus,inputMovie,inputMin); case 'zeroToOne' @@ -184,9 +192,8 @@ % reverseStr = ''; for frameNo2 = 1:nFrames thisFrame = squeeze(inputMovie(:,:,frameNo2)); - maxVec = nanmax(thisFrame(:)); - minVec = nanmin(thisFrame(:)); - % meanVec = nanmean(thisFrame(:)); + maxVec = max(thisFrame,[],[1 2],'omitnan'); + minVec = min(thisFrame,[],[1 2],'omitnan'); inputMovie(:,:,frameNo2) = (thisFrame-minVec)./(maxVec-minVec); if ~verLessThan('matlab', '9.2') send(D, frameNo2); % Update @@ -233,21 +240,20 @@ function mijiCheck() end function [movieTmp] = addText(movieTmp,inputText,fontSize) nFrames = size(movieTmp,3); - maxVal = nanmax(movieTmp(:)); - minVal = nanmin(movieTmp(:)); + maxVal = max(movieTmp,[],[1 2 3],'omitnan'); + minVal = min(movieTmp,[],[1 2 3],'omitnan'); reverseStr = ''; for frameNo = 1:nFrames - movieTmp(:,:,frameNo) = squeeze(nanmean(... - insertText(movieTmp(:,:,frameNo),[0 0],num2str(inputText(frameNo)),... + tmpFrame = insertText(movieTmp(:,:,frameNo),[0 0],num2str(inputText(frameNo)),... 'BoxColor',[maxVal maxVal maxVal],... 'TextColor',[minVal minVal minVal],... 'AnchorPoint','LeftTop',... 'FontSize',fontSize,... - 'BoxOpacity',1)... - ,3)); + 'BoxOpacity',1); + movieTmp(:,:,frameNo) = squeeze(mean(tmpFrame,3,'omitnan')); reverseStr = cmdWaitbar(frameNo,nFrames,reverseStr,'inputStr','adding text to movie','waitbarOn',1,'displayEvery',10); end - % maxVal = nanmax(movieTmp(:)) + % maxVal = max(movieTmp,[],[1 2 3],'omitnan') % movieTmp(movieTmp==maxVal) = 1; % 'BoxColor','white' end @@ -455,9 +461,6 @@ function subfxnImagejFFT_test() end function subfxnFft() disp('Running Matlab FFT') - if options.runNanCheck==1 - inputMovie(isnan(inputMovie)) = 0; - end % bandpassMatrix = zeros(size(inputMovie)); % get options ioptions.showImages = options.showImages; @@ -512,12 +515,15 @@ function subfxnFft() d = []; end - if isempty(gcp) - opts2.Value = ioptions; - else - opts2 = parallel.pool.Constant(ioptions); - end + %if isempty(gcp) + % opts2.Value = ioptions; + %else + % opts2 = parallel.pool.Constant(ioptions); + %end + opts2 = parallel.pool.Constant(ioptions); % startState = ticBytes(gcp); + options_runNanCheck = options.runNanCheck; + parfor frame = 1:nFramesToNormalize % [percent progress] = parfor_progress;if mod(progress,dispStepSize) == 0;dispstat(sprintf('progress %0.1f %',percent));else;end % thisFrame = squeeze(inputMovie(:,:,frame)); @@ -532,6 +538,11 @@ function subfxnFft() thisFrame = inputMovie(:,:,frame); thisFrame = squeeze(thisFrame); + if options_runNanCheck==1 + nanMask = isnan(thisFrame); + thisFrame(nanMask) = mean(thisFrame,[1 2 3],'omitnan'); + end + if isempty(secondaryNormalizationType) if isempty(coords) thisFrame = fftImage(thisFrame,'options',opts2.Value); @@ -545,13 +556,17 @@ function subfxnFft() tmpFrame = thisFrame; tmpFrame(c:d, a:b) = fftImage(thisFrame(c:d, a:b),'options',opts2.Value); end - tmpFrameMin = nanmin(tmpFrame(:)); + tmpFrameMin = min(tmpFrame,[],[1 2],'omitnan'); if tmpFrameMin<0 tmpFrame = tmpFrame + abs(tmpFrameMin); end thisFrame = thisFrame./tmpFrame; end + if options_runNanCheck==1 + thisFrame(nanMask) = NaN; + end + inputMovie(:,:,frame) = thisFrame; if ~verLessThan('matlab', '9.2') @@ -636,7 +651,7 @@ function subfxnMatlabFFT_test() % FFT each image reverseStr = ''; % options.secondaryNormalizationType - % nanmean(inputMovieFFT(:)) + % mean(inputMovieFFT,[1 2 3],'omitnan') if options.runNanCheck==1 inputMovieFFT(isnan(inputMovieFFT)) = 0; end @@ -657,27 +672,41 @@ function subfxnMatlabFFT_test() % end % inputMovie = cat(3,inputMovie{:}); + options_runNanCheck = options.runNanCheck; + for frame=1:options.maxFrame thisFrame = squeeze(inputMovieFFT(:,:,frame)); - if isempty(options.secondaryNormalizationType) - inputMovieFFT(:,:,frame) = fftImage(thisFrame,'options',ioptions); - else - % tmpFrame = fftImage(thisFrame,'options',ioptions); - % inputMovieFFT(:,:,frame) = thisFrame./tmpFrame; - inputMovieFFT(:,:,frame) = fftImage(thisFrame,'options',ioptions); + + if options_runNanCheck==1 + nanMask = isnan(thisFrame); + thisFrame(nanMask) = mean(thisFrame,[1 2 3],'omitnan'); end + + thisFrame = fftImage(thisFrame,'options',ioptions); + % if isempty(options.secondaryNormalizationType) + % else + % % tmpFrame = fftImage(thisFrame,'options',ioptions); + % % inputMovieFFT(:,:,frame) = thisFrame./tmpFrame; + % thisFrame = fftImage(thisFrame,'options',ioptions); + % end + + if options_runNanCheck==1 + thisFrame(nanMask) = NaN; + end + + inputMovieFFT(:,:,frame) = thisFrame; % bandpassMatrix(:,:,frame) = fftImage(thisFrame,'options',ioptions); % bandpassMatrix(:,:,frame) = imcomplement(bandpassMatrix(:,:,frame)); reverseStr = cmdWaitbar(frame,options.maxFrame,reverseStr,'inputStr','normalizing movie','waitbarOn',options.waitbarOn,'displayEvery',5); % = bsxfun(@ldivide,squeeze(movie20hz(:,:,1)),filteredFrame end % inputImageTest(:,:,freqNo) = inputMovieFFT; - inputMovieFFTMin = nanmin(inputMovieFFT(:)); + inputMovieFFTMin = min(inputMovieFFT,[],[1 2 3],'omitnan'); inputImageTest{freqNo} = inputMovieFFT - inputMovieFFTMin+1; inputMovieDuplicate{freqNo} = cast(inputMovie,outputClass); % divide lowpass from image if options.testDuplicateDfof == 1 - minMovie = min(inputMovieDuplicate{freqNo}(:)); + minMovie = min(inputMovieDuplicate{freqNo},[],[1 2 3],'omitnan'); if minMovie<0 inputMovieDuplicate{freqNo} = inputMovieDuplicate{freqNo} + 1.1*abs(minMovie); end @@ -688,8 +717,8 @@ function subfxnMatlabFFT_test() % else % inputMovieDivide{freqNo} = inputImageTest{freqNo}; % end - % nanmean(inputMovieDuplicate{freqNo}(:)) - % nanmean(inputImageTest{freqNo}(:)) + % mean(inputMovieDuplicate{freqNo}(:),[],'omitnan') + % mean(inputImageTest{freqNo}(:),[],'omitnan') % inputMovieDivide{freqNo} = bsxfun(@minus,inputMovieDuplicate{freqNo},inputImageTest{freqNo}); inputMovieDfof{freqNo} = dfofMovie(inputMovieDivide{freqNo}); else @@ -724,10 +753,10 @@ function subfxnDetrend() frameMeanInputMovie = squeeze(mean(inputMovie,[1 2],'omitnan')); - trendVals = frameMeanInputMovie - detrend(frameMeanInputMovie,option.detrendDegree); + trendVals = frameMeanInputMovie - detrend(frameMeanInputMovie,options.detrendDegree); % meanInputMovie = detrend(frameMeanInputMovie,0); - meanInputMovie = nanmean(frameMeanInputMovie); + meanInputMovie = squeeze(mean(frameMeanInputMovie,[1 2 3],'omitnan')); nFramesToNormalize = options.maxFrame; nFrames = nFramesToNormalize; @@ -740,61 +769,106 @@ function subfxnDetrend() % opts2 = parallel.pool.Constant(ioptions); % end - parfor frame = 1:nFramesToNormalize - inputMovie(:,:,frame) = inputMovie(:,:,frame) - trendVals(frame) + meanInputMovie; + % parfor frame = 1:nFramesToNormalize + % inputMovie(:,:,frame) = inputMovie(:,:,frame) - trendVals(frame) + meanInputMovie; - % thisFrame = inputMovie(:,:,frame); - % thisFrame = squeeze(thisFrame); - % thisFrame = thisFrame - trendVals(frame); - % thisFrame = thisFrame + meanInputMovie; + % % thisFrame = inputMovie(:,:,frame); + % % thisFrame = squeeze(thisFrame); + % % thisFrame = thisFrame - trendVals(frame); + % % thisFrame = thisFrame + meanInputMovie; - % inputMovie(:,:,frame) = thisFrame; + % % inputMovie(:,:,frame) = thisFrame; - if ~verLessThan('matlab', '9.2') - send(D, frame); % Update + % if ~verLessThan('matlab', '9.2') + % send(D, frame); % Update + % end + % end + + % Detrend each frame only using subset to reduce memory usage + subsetSize = options.subsetSize; + movieLength = size(inputMovie,3); + numSubsets = ceil(movieLength/subsetSize)+1; + subsetList = round(linspace(1,movieLength,numSubsets)); + disp(['Detrending sublists: ' num2str(subsetList)]); + nSubsets = (length(subsetList)-1); + + for thisSet = 1:nSubsets + subsetStartIdx = subsetList(thisSet); + subsetEndIdx = subsetList(thisSet+1); + display(repmat('$',1,7)) + if thisSet==nSubsets + movieSubset = subsetStartIdx:subsetEndIdx; + else + movieSubset = subsetStartIdx:(subsetEndIdx-1); + end + disp([num2str(movieSubset(1)) '-' num2str(movieSubset(end)) ' ' num2str(thisSet) '/' num2str(nSubsets)]) + disp(repmat('$',1,7)) + + inputMovieTmp = inputMovie(:,:,movieSubset); + trendValsTmp = trendVals(movieSubset); + nFramesTmp = size(inputMovieTmp,3); + % movieSubset2 = 1:nFramesTmp; + + parfor frameNo = 1:nFramesTmp + % if frameNo==startFrame||mod(frameNo,50)==0 + % fprintf('%d | ',frameNo); + % end + inputMovieTmp(:,:,frameNo) = inputMovieTmp(:,:,frameNo) - trendValsTmp(frameNo) + meanInputMovie; end + disp('==='); + inputMovie(:,:,movieSubset) = inputMovieTmp; + fprintf('\n'); end end function subfxnMatlabDisk() disp('Matlab fspecial disk background removal') %Get dimension information about 3D movie matrix - [inputMovieX, inputMovieY, inputMovieZ] = size(inputMovie); + % [inputMovieX, inputMovieY, inputMovieZ] = size(inputMovie); % reshapeValue = size(inputMovie); %Convert array to cell array, allows slicing (not contiguous memory block) - inputMovie = squeeze(mat2cell(inputMovie,inputMovieX,inputMovieY,ones(1,inputMovieZ))); + %inputMovie = squeeze(mat2cell(inputMovie,inputMovieX,inputMovieY,ones(1,inputMovieZ))); - imageNow = squeeze(inputMovie{1}); + %imageNow = squeeze(inputMovie{1}); + imageNow = squeeze(inputMovie(:,:,1)); [rows,cols] = size(imageNow); r1 = min(rows,cols)/10; r2 = 3; hDisk = fspecial('disk', r1); hDisk2 = fspecial('disk', r2); + + C_hDisk = parallel.pool.Constant(hDisk); + C_hDisk2 = parallel.pool.Constant(hDisk2); + transform = @(A) transform_2(A,hDisk,hDisk2); - reverseStr = ''; + %reverseStr = ''; % nImages = size(inputMovieCropped,3); - nImages = length(inputMovie); + nImages = size(inputMovie,3); + nFrames = nImages; % [percent progress] = parfor_progress(nImages); dispStepSize = round(nImages/20); dispstat('','init'); - parfor imageNo = 1:nImages + parfor frame = 1:nImages % [percent progress] = parfor_progress;if mod(progress,dispStepSize) == 0;dispstat(sprintf('progress %0.1f %',percent));else;end - imageNow = squeeze(inputMovie{imageNo}); - inputMovie{imageNo} = transform(imageNow); + % imageNow = squeeze(inputMovie{imageNo}); + thisFrame = inputMovie(:,:,frame); + thisFrame = squeeze(thisFrame); + thisFrame = transform_2(thisFrame,C_hDisk.Value,C_hDisk2.Value); + % imageNow = transform(imageNow); + % inputMovie{imageNo} = transform(imageNow); + inputMovie(:,:,frame) = thisFrame; % if (mod(imageNo,20)==0|imageNo==nImages) % reverseStr = cmdWaitbar(imageNo,nImages,reverseStr,'inputStr','fspecial normalizing'); % end - end - parfor_progress(0);dispstat('Finished.','keepprev'); - - inputMovie = cat(3,inputMovie{:}); - end - function A_tr = transform_2(A, ssm_filter, asm_filter) - - A_tr = A - imfilter(A, ssm_filter, 'replicate'); - A_tr = imfilter(A_tr, asm_filter); + if ~verLessThan('matlab', '9.2') + send(D, frame); % Update + end + end + %parfor_progress(0); + %dispstat('Finished.','keepprev'); - end + % inputMovie = cat(3,inputMovie{:}); + end function subfxnMatlabFFTTemporal() sigma = 1.5; fsize = 10; @@ -846,14 +920,23 @@ function subfxnImfilterSmooth() % reshapeValue = size(inputMovie); %Convert array to cell array, allows slicing (not contiguous memory block) % inputMovie = squeeze(mat2cell(inputMovie,inputMovieX,inputMovieY,ones(1,inputMovieZ))); - runNanCheck = options.runNanCheck; + options_runNanCheck = options.runNanCheck; parfor frame=1:nFrames thisFrame = squeeze(inputMovie(:,:,frame)); % thisFrame = inputMovie{frame}; - if runNanCheck==1 - thisFrame(isnan(thisFrame)) = nanmean(thisFrame(:)); + + if options_runNanCheck==1 + nanMask = isnan(thisFrame); + thisFrame(nanMask) = mean(thisFrame,[1 2 3],'omitnan'); + end + + thisFrame = imfilter(thisFrame, movieFilter,options_boundaryType); + + if options_runNanCheck==1 + thisFrame(nanMask) = NaN; end - inputMovie(:,:,frame) = imfilter(thisFrame, movieFilter,options_boundaryType); + + inputMovie(:,:,frame) = thisFrame % inputMovie{frame} = imfilter(thisFrame, movieFilter,options_boundaryType); % reverseStr = cmdWaitbar(frame,nFrames,reverseStr,'inputStr','imfilter normalizing movie','waitbarOn',options.waitbarOn,'displayEvery',5); if ~verLessThan('matlab', '9.2') @@ -877,13 +960,23 @@ function subfxnImfilter() nFrames = size(inputMovie,3); inputMovieFiltered = zeros(size(inputMovie),class(inputMovie)); reverseStr = ''; - runNanCheck = options.runNanCheck; + options_runNanCheck = options.runNanCheck; for frame=1:nFrames thisFrame = squeeze(inputMovie(:,:,frame)); - if runNanCheck==1 - thisFrame(isnan(thisFrame))=nanmean(thisFrame(:)); + + if options_runNanCheck==1 + nanMask = isnan(thisFrame); + thisFrame(nanMask) = mean(thisFrame,[1 2 3],'omitnan'); + end + + thisFrame = imfilter(thisFrame, movieFilter,options.boundaryType); + + if options_runNanCheck==1 + thisFrame(nanMask) = NaN; end - inputMovieFiltered(:,:,frame) = imfilter(thisFrame, movieFilter,options.boundaryType); + + inputMovieFiltered(:,:,frame) = thisFrame; + reverseStr = cmdWaitbar(frame,nFrames,reverseStr,'inputStr','normalizing movie','waitbarOn',options.waitbarOn,'displayEvery',5); end % divide each frame by the filtered movie to remove 'background' @@ -908,4 +1001,11 @@ function subfxnImfilter() % end % end % end +end +function A_tr = transform_2(A, ssm_filter, asm_filter) + + A_tr = A - imfilter(A, ssm_filter, 'replicate'); + + A_tr = imfilter(A_tr, asm_filter); + end \ No newline at end of file diff --git a/+ciapkg/+overloaded/suptitle.m b/+ciapkg/+overloaded/suptitle.m index a759572..29ed408 100644 --- a/+ciapkg/+overloaded/suptitle.m +++ b/+ciapkg/+overloaded/suptitle.m @@ -16,6 +16,8 @@ % Warning: If the figure or axis units are non-default, this % will break. + % changelog + % 2023.08.08 [05:27:57] - Added font name option %======================== % Amount of the figure window devoted to subplots @@ -28,6 +30,8 @@ options.Color = 'k'; % Str: interpreter, latex or tex options.Interpreter = 'latex'; + % Str: Name of font family, e.g. Consolas. + options.fontName = []; % get options options = ciapkg.io.getOptions(options,varargin); % display(options) @@ -122,7 +126,11 @@ warning on ht = text('position',[.5 titleypos-1],'Interpreter',options.Interpreter,'String',str); set(ht,'horizontalalignment','center','fontsize',fs,'Color',options.Color); - set(ht,'FontName',get(0,'DefaultAxesFontName')); + if ~isempty(options.fontName) + set(ht,'FontName',options.fontName); + else + set(ht,'FontName',get(0,'DefaultAxesFontName')); + end set(gcf,'nextplot',np); % focus stealing occurs here, turn off set(gcf, 'CurrentAxes', haold) diff --git a/+ciapkg/+settings/cnmfeSettings.m b/+ciapkg/+settings/cnmfeSettings.m index 914eb39..b649749 100644 --- a/+ciapkg/+settings/cnmfeSettings.m +++ b/+ciapkg/+settings/cnmfeSettings.m @@ -2,6 +2,8 @@ % 2019.04.04 [10:35:15] % Settings to run CNMF-E, used in conjunction with modelExtractSignalsFromMovie.m or % as input to "computeCnmfeSignalExtraction" using "computeCnmfeSignalExtraction(inputMovie,'options',cnmfeOpts)" after running "cnmfeSettings" +% changelog + % 2024.03.29 [19:02:12] - Updated merge settings to reduce duplicates in output. % ======================== % OVERALL @@ -22,7 +24,7 @@ cnmfeOpts.batch_frames = []; % ===SPATIAL % Int: pixel, gaussian width of a gaussian kernel for filtering the data. 0 means no filtering -cnmfeOpts.gSig = 3; +cnmfeOpts.gSig = 7; % Int: pixel, neuron diameter cnmfeOpts.gSiz = 13; % Int: spatial downsampling factor @@ -55,15 +57,15 @@ cnmfeOpts.bg_ssub = 2; % ===MERGING % Float: 0 to 1, thresholds for merging neurons; [spatial overlap ratio, temporal correlation of calcium traces, spike correlation] -cnmfeOpts.merge_thr = 0.65; +cnmfeOpts.merge_thr = 0.2; % Char: method for computing neuron distances {'mean', 'max'} -cnmfeOpts.method_dist = 'max'; +cnmfeOpts.method_dist = 'mean'; % Int: minimum distances between two neurons. it is used together with merge_thr -cnmfeOpts.dmin = 5; +cnmfeOpts.dmin = 10; % Int: merge neurons if their distances are smaller than dmin_only. -cnmfeOpts.dmin_only = 2; +cnmfeOpts.dmin_only = 10; % Float vector: merge components with highly correlated spatial shapes (corr=0.8) and small temporal correlations (corr=0.1) -cnmfeOpts.merge_thr_spatial = [0.8, 0.4, -inf]; +cnmfeOpts.merge_thr_spatial = [0.3, 0.1, -inf]; % ===INITIALIZATION % Int: maximum number of neurons per patch. when K=[], take as many as possible. cnmfeOpts.K = []; diff --git a/+ciapkg/+signal_extraction/computeSignalsFromImages.m b/+ciapkg/+signal_extraction/computeSignalsFromImages.m index 0828d24..4aa1b87 100644 --- a/+ciapkg/+signal_extraction/computeSignalsFromImages.m +++ b/+ciapkg/+signal_extraction/computeSignalsFromImages.m @@ -15,7 +15,7 @@ % TODO % - import ciapkg.api.* % import CIAtah functions in ciapkg package API. + % import ciapkg.api.* % import CIAtah functions in ciapkg package API. % ======================== % DESCRIPTION @@ -32,7 +32,7 @@ % ======================== try - [outputSignal, inputImages] = applyImagesToMovie(inputImages,inputMovie, 'passArgs',varargin); + [outputSignal, inputImages] = ciapkg.signal_processing.applyImagesToMovie(inputImages,inputMovie, 'passArgs',varargin); catch err disp(repmat('@',1,7)) disp(getReport(err,'extended','hyperlinks','on')); diff --git a/+ciapkg/+tracking/annotateOpenFieldParameters.m b/+ciapkg/+tracking/annotateOpenFieldParameters.m index 30a52e9..9a6bf0b 100644 --- a/+ciapkg/+tracking/annotateOpenFieldParameters.m +++ b/+ciapkg/+tracking/annotateOpenFieldParameters.m @@ -19,6 +19,9 @@ % Changelog % 2022.12.13 [09:41:48] - add FPS override for cases in which camera metadata is incorrect. % 2023.04.29 [18:47:02] - Arena size automatically calculated from reference distance size. Add GUI font size option. + % TODO + % Clean up code to make into more nested and sub functions + % Extend to elevated plus maze % ======================== % Str: regular expression of files to find if inputFilePath is a folder. diff --git a/+ciapkg/+video/createMontageMovie.m b/+ciapkg/+video/createMontageMovie.m index 924fed7..ddb8fb5 100644 --- a/+ciapkg/+video/createMontageMovie.m +++ b/+ciapkg/+video/createMontageMovie.m @@ -130,6 +130,7 @@ inputMovieNo = 1; for xNo = 1:xPlot for yNo = 1:yPlot + disp('---') if inputMovieNo>length(inputMovies) [behaviorMovie{xNo}] = createSideBySide(behaviorMovie{xNo},NaN(size(inputMovies{1})),'pxToCrop',[],'makeTimeEqualUsingNans',1,'normalizeMovies',0,'displayInfo',options.displayInfo,'increaseToLargestMovie',options.increaseToLargestMovie); elseif yNo==1 @@ -154,13 +155,13 @@ % behaviorMovie = cat(behaviorMovie{:},3) % do something catch err - display(repmat('@',1,7)) + disp(repmat('@',1,7)) disp(getReport(err,'extended','hyperlinks','on')); - display(repmat('@',1,7)) + disp(repmat('@',1,7)) end function subfxnDisp(txt) if options.displayInfo==1 - display(txt); + disp(txt); end end end diff --git a/+ciapkg/+video/createSideBySide.m b/+ciapkg/+video/createSideBySide.m index f931b3d..d520454 100644 --- a/+ciapkg/+video/createSideBySide.m +++ b/+ciapkg/+video/createSideBySide.m @@ -235,14 +235,14 @@ close(writerObj); end catch err - display(repmat('@',1,7)) + disp(repmat('@',1,7)) disp(getReport(err,'extended','hyperlinks','on')); - display(repmat('@',1,7)) + disp(repmat('@',1,7)) end % subfxnDisp: function description function subfxnDisp(txt) if options.displayInfo==1 - display(txt) + disp(txt) end end end diff --git a/+ciapkg/+video/scaleMatrix8bit.m b/+ciapkg/+video/scaleMatrix8bit.m new file mode 100644 index 0000000..1dbf712 --- /dev/null +++ b/+ciapkg/+video/scaleMatrix8bit.m @@ -0,0 +1,110 @@ +function [inputMovie] = scaleMatrix8bit(inputMovie,varargin) + % [inputMovie] = scaleMatrix8bit(inputMovie,varargin) + % + % Scales an input matrix to 8-bit (256) range. Can also save video out to an AVI file. + % + % Biafra Ahanonu + % started: 2022.12.02 [10:19:00] + % + % Inputs + % inputMovie - [x y nFrames] tensor with movie data. + % + % Outputs + % inputMovie - [x y nFrames] tensor with corrected movie data. + % + % Options (input as Name-Value with Name = options.(Name)) + % % DESCRIPTION + % options.exampleOption = ''; + + % Changelog + % 2024.02.11 [15:27:08] - Updates to function and integration into CIAtah. + % TODO + % + + % ======================== + % Int: amount to downsample movie by. 1 = no downsampling. + options.dsFactor = 1; + % Binary: 1 = display plots. + options.dispPlots = 0; + % Str: Path to AVI to save output to + options.savePath = ''; + % Binary: 1 = ignore zeros when calculating min/max, in cases of movies with a lot of motion (e.g. borders = 0) or other reasons for excessive zeros. + options.ignoreZeros = 0; + % get options + options = ciapkg.io.getOptions(options,varargin); + % disp(options) + % unpack options into current workspace + % fn=fieldnames(options); + % for i=1:length(fn) + % eval([fn{i} '=options.' fn{i} ';']); + % end + % ======================== + + try + %% + opts.dsFactor = options.dsFactor; + + prctileMax = 99.99; + prctileMin = 1; + multiRatio = 1.2; + adjustMinus = 0.01; + + % Normalize and correct contrast for AVI for improved DLC tracking + inputMovie = single(inputMovie); + + if opts.dsFactor~=1 + inputMovie = ciapkg.movie_processing.downsampleMovie(inputMovie,... + 'downsampleDimension','space','downsampleFactor',opts.dsFactor); + end + + % Obtain percentile min/max, less likely to be thrown off by extreme values + if options.ignoreZeros==1 + movieMax = prctile(inputMovie,prctileMax,[1 2 3]); + movieMin = prctile(inputMovie,prctileMin,[1 2 3]); + else + inputMovieTmp = inputMovie; + inputMovieTmp(inputMovieTmp==0) = NaN; + movieMax = prctile(inputMovieTmp,prctileMax,[1 2 3]); + movieMin = prctile(inputMovieTmp,prctileMin,[1 2 3]); + end + + % Normalize the movie + inputMovie = (inputMovie-movieMin)/(movieMax-movieMin); + inputMovie = (inputMovie*multiRatio-adjustMinus)*255; + inputMovie = uint8(inputMovie); + + if options.dispPlots==1 + figure + % subplot(xPlot,yPlot,i) + imagesc(squeeze(inputMovie(:,:,1))) + axis image + box off + colormap gray + title(i) + drawnow + end + if ~isempty(options.savePath) + % [~,saveFileName,~] = fileparts(thisMoviePath); + % outputSavePath = fullfile(outputPath,[saveFileName '.avi']); + + outputSavePath = options.savePath; + fprintf('Save to: %s\n',outputSavePath) + ciapkg.io.saveMatrixToFile(inputMovie,outputSavePath); + end + + disp('Done!') + catch err + disp(repmat('@',1,7)) + disp(getReport(err,'extended','hyperlinks','on')); + disp(repmat('@',1,7)) + end + + % function [outputs] = nestedfxn_exampleFxn(arg) + % % Always start nested functions with "nestedfxn_" prefix. + % % outputs = ; + % end +end +% function [outputs] = localfxn_exampleFxn(arg) +% % Always start local functions with "localfxn_" prefix. +% % outputs = ; +% end diff --git a/+ciapkg/+view/changeFont.m b/+ciapkg/+view/changeFont.m index 44428ae..230abf0 100644 --- a/+ciapkg/+view/changeFont.m +++ b/+ciapkg/+view/changeFont.m @@ -16,6 +16,8 @@ % 2022.01.14 [05:53:37] - Updated so doesn't change Axes backgroundcolor when changing font color, only Axes text. % 2022.03.14 [04:06:10] - Also check for matlab.ui.control.UIControl when conducting font color changes and ignore to not cause errors. % 2022.12.03 [21:51:44] - Update to handle legend, since use TextColor property instead of Color for text. Users can also put in a blank fontSize to avoid changing that property, can skip UI that would normally appear. + % TODO + % Add support for changing other font aspects, e.g. figure Font family, command window font, etc. import ciapkg.api.* % import CIAtah functions in ciapkg package API. diff --git a/+ciapkg/+view/displacementFieldCompare.m b/+ciapkg/+view/displacementFieldCompare.m new file mode 100644 index 0000000..f2bd724 --- /dev/null +++ b/+ciapkg/+view/displacementFieldCompare.m @@ -0,0 +1,123 @@ +function [success] = displacementFieldCompare(fixedTmp,movingTmp,movingReg,Dxy,varargin) + % [output1,output2] = displacementFieldCompare(fixedTmp,movingTmp,movingReg,Dxy,varargin) + % + % Displays different frames along with the displacement fields. + % + % Biafra Ahanonu + % started: 2024.03.01 [14:58:30] + % + % Inputs + % fixedTmp + % movingTmp + % movingReg + % Dxy + % + % Outputs + % success - Binary: 1 = successfully run + % + % Options (input as Name-Value with Name = options.(Name)) + % % DESCRIPTION + % options.exampleOption = ''; + + % Changelog + % 2022.03.14 [01:47:04] - Added nested and local functions to the example function. + % TODO + % + + % ======================== + % DESCRIPTION + options.frameA = 1; + options.frameB = 2; + % get options + options = ciapkg.io.getOptions(options,varargin); + % disp(options) + % unpack options into current workspace + % fn=fieldnames(options); + % for i=1:length(fn) + % eval([fn{i} '=options.' fn{i} ';']); + % end + % ======================== + + try + %% DISPLAY PLOTS + frameA = options.frameA; + frameB = options.frameB; + % + success = 0; + subplotTmp = @(x,y,z) subaxis(x,y,z, 'Spacing', 0.005, 'Padding', 0.02, 'PaddingTop', 0.02, 'MarginTop', 0.03,'MarginBottom', 0.03,'MarginLeft', 0.01,'MarginRight', 0.01); + figure; + colormap gray + xT = 4; + axList = []; + axList(end+1) = subplotTmp(2,xT,1); + imagesc(fixedTmp) + axis image + box off; + title(['Reference image | ' num2str(frameA)]) + axList(end+1) = subplotTmp(2,xT,2); + imagesc(movingTmp) + box off + axis image + title(['Moving image | ' num2str(frameB)]) + axList(end+1) = subplotTmp(2,xT,3); + imshowpair(fixedTmp,movingTmp) + axis image + title('Overlap of reference and moving') + axList(end+1) = subplotTmp(2,xT,4); + imagesc(movingReg) + box off; + title(['Registered moving image | ' num2str(frameB)]) + axis image + axList(end+1) = subplotTmp(2,xT,5); + imshowpair(fixedTmp,movingReg) + axis image + title('Overlap of reference and registered moving') + axList(end+1) = subplotTmp(2,xT,6); + imagesc(Dxy(:,:,1)) + box off; + colorbar + axis image + title('X-axis displacement field') + axList(end+1) = subplotTmp(2,xT,7); + imagesc(Dxy(:,:,2)) + box off; + colorbar + axis image + title('Y-axis displacement field') + axList(end+1) = subplotTmp(2,xT,8); + imagesc(movingTmp) + box off + title(['Moving image w/ displacement fields | ' num2str(frameB)]) + hold on; + sub1 = 5; + Dxy_ds = ciapkg.movie_processing.downsampleMovie(Dxy,'downsampleDimension','space','downsampleFactor',sub1); + %[Xdim,Ydim] = meshgrid(1:size(Dxy_ds,2),1:size(Dxy_ds,1)); + [Xdim,Ydim] = meshgrid(round(linspace(1,size(Dxy,2),size(Dxy,2)/sub1)),round(linspace(1,size(Dxy,1),size(Dxy,1)/sub1))); + Udim = squeeze(Dxy_ds(:,:,1)); + Vdim = squeeze(Dxy_ds(:,:,2)); + %sub1 = 10; + %[Xdim,Ydim] = meshgrid(1:sub1:size(Dxy,2),1:size(Dxy,1)); + % Plot the vectors, invert both U and V for proper display + hq = quiver(Xdim,Ydim,-Udim,-Vdim,1,'-r','LineWidth',1); + axis equal + xlim([0 size(movingTmp,2)]) + ylim([0 size(movingTmp,1)]) + linkaxes(axList) + ciapkg.view.changeFont(18) + disp('Done!') + success = 0; + catch err + disp(repmat('@',1,7)) + disp(getReport(err,'extended','hyperlinks','on')); + disp(repmat('@',1,7)) + end + + % function [outputs] = nestedfxn_exampleFxn(arg) + % % Always start nested functions with "nestedfxn_" prefix. + % % outputs = ; + % end +end +% function [outputs] = localfxn_exampleFxn(arg) +% % Always start local functions with "localfxn_" prefix. +% % outputs = ; +% end \ No newline at end of file diff --git a/+ciapkg/+view/playMovie.m b/+ciapkg/+view/playMovie.m index 89875b4..8a6156d 100644 --- a/+ciapkg/+view/playMovie.m +++ b/+ciapkg/+view/playMovie.m @@ -12,6 +12,7 @@ % RGB: [x y C t] matrix of x,y height/width, t frames, and C color channel (3 as RGB) % Path: full file path to a movie containing a [x y t] or [x y C t] matrix. % Path (MAT-file): full path to a MAT-file with a variable containing a [x y t] or [x y C t] matrix. + % Cell array: cell array of [x y t] matrix or path to movies movies. % options % fps - frame rate to display movie. % extraMovie - extra movie to play, [X Y Z] matrix of X,Y height/width and Z frames @@ -52,6 +53,12 @@ % 2023.02.22 [18:24:46] - Added downsample (both type) support for extra movie, mirrors what happens in the 1st movie. Also improved downsampling when just subsampling by removing unnecessary calls to imresize. Improved tracking point handling to ensure in the correct axes and can turn off prior frame tracking. Better handling of quiver plots with tracking. % 2023.04.06 [19:29:42] - Convert all nanmin/nanmax to 'omitnan' and (:) to [1 2 3 4] to ensure compatible with 3- and 4-d tensors. % 2023.10.23 [17:46:25] - Additional Bio-Formats support. + % 2024.02.19 [11:07:30] - Users can now input a cell array of movies to have the extra movie play instead of the 'extraMovie' option. + % 2024.03.29 [09:45:20] - Added support for additional NWB series. + % TODO + % Add support for user clicking the movie to add marker points for reference (e.g. to help with visualizing motion correction accuracy). + % Generalize extra movie so users can input 3, 4, etc. movies. + % Add contrast bar for extra movie. import ciapkg.api.* % import CIAtah functions in ciapkg package API. @@ -130,7 +137,7 @@ % Str: hierarchy name in hdf5 where movie data is located. options.inputDatasetName = '/1'; % Str: default NWB hierarchy names in HDF5 file where movie data is located, will look in the order indicates - options.defaultNwbDatasetName = {'/acquisition/TwoPhotonSeries/data'}; + options.defaultNwbDatasetName = {'/acquisition/TwoPhotonSeries/data','/acquisition/OnePhotonSeries/data'}; % Int: [] = do nothing, 1-3 indicates R,G,B channels to take from multicolor RGB AVI options.rgbChannel = []; % Int vector: list of specific frames to load. @@ -173,6 +180,12 @@ rgbMovieFlag = 0; + % If a cell array is input, convert to format playMovie will use. + if iscell(inputMovie)==1 + options.extraMovie = inputMovie{2}; + inputMovie = inputMovie{1}; + end + % Obtain movie information and connect to file if user gives a path to a movie. if ischar(inputMovie)==1 inputMovieIsChar = 1; @@ -388,7 +401,7 @@ % axHandle.Toolbar.Visible = 'off'; box off; if ~isempty(options.extraLinePlot) - axis equal tight + axis equal tight; else axis equal tight end diff --git a/+ciapkg/VERSION b/+ciapkg/VERSION index 60a2ace..03e6beb 100644 --- a/+ciapkg/VERSION +++ b/+ciapkg/VERSION @@ -1,2 +1,2 @@ -v4.6.1 -2023.12.17 [15:08:16] \ No newline at end of file +v4.11.0 +2024.07.23 [13:45:12] \ No newline at end of file diff --git a/@ciatah/computeManualSortSignals.m b/@ciatah/computeManualSortSignals.m index 9d0441f..03aa41a 100644 --- a/@ciatah/computeManualSortSignals.m +++ b/@ciatah/computeManualSortSignals.m @@ -16,6 +16,7 @@ % 2021.06.18 [21:41:07] - added modelVarsFromFilesCheck() to check and load signals if user hasn't already. % 2021.06.21 [21:03:25] - Fix check to make sure variables are loaded. % 2021.08.10 [09:57:36] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package. + % 2024.03.27 [15:41:31] - Update display of information. % ADDED % ADD PERSONS NAME TO THE FILE - DONE. % TODO @@ -501,18 +502,19 @@ uiTextHandles = {}; uiXIncrement = 0.03; uiYOffset = 0.90; + uiXOffset = 0.02; uiTxtSize = 0.3; uiBoxSize = 0.65; [figHandle figNo] = openFigure(1337, ''); clf - uicontrol('Style','Text','String',inputTitleStr,'Units','normalized','Position',[0.0 uiYOffset-uiXIncrement*(0) 0.3 0.05],'BackgroundColor','white','HorizontalAlignment','Left'); + uicontrol('Style','Text','String',inputTitleStr,'Units','normalized','Position',[uiXOffset uiYOffset-uiXIncrement*(0) 0.3 0.05],'BackgroundColor','white','HorizontalAlignment','Left'); for propertyNo = 1:nPropertiesToChange property = char(propertyList(propertyNo)); - uiTextHandles{propertyNo} = uicontrol('Style','text','String',[regSettingTitles.(property) '' 10],'Units','normalized','Position',[0.0 uiYOffset-uiXIncrement*propertyNo+0.03 uiTxtSize 0.02],'BackgroundColor',[0.9 0.9 0.9],'ForegroundColor','black','HorizontalAlignment','Left'); + uiTextHandles{propertyNo} = uicontrol('Style','text','String',[regSettingTitles.(property) '' 10],'Units','normalized','Position',[uiXOffset uiYOffset-uiXIncrement*propertyNo+0.03 uiTxtSize 0.02],'BackgroundColor',[0.9 0.9 0.9],'ForegroundColor','black','HorizontalAlignment','Left'); % uiTextHandles{propertyNo}.Enable = 'Inactive'; - uiListHandles{propertyNo} = uicontrol('Style', 'popup','String', propertySettingsStr.(property),'Units','normalized','Position', [uiTxtSize uiYOffset-uiXIncrement*propertyNo uiBoxSize 0.05],'Callback',@subfxnSettingsCallback,'Tag',property); + uiListHandles{propertyNo} = uicontrol('Style', 'popup','String', propertySettingsStr.(property),'Units','normalized','Position', [uiTxtSize+uiXOffset uiYOffset-uiXIncrement*propertyNo uiBoxSize 0.05],'Callback',@subfxnSettingsCallback,'Tag',property); end - uicontrol('Style','Text','String','press enter to continue','Units','normalized','Position',[0.0 uiYOffset-uiXIncrement*(nPropertiesToChange+2) 0.3 0.05],'BackgroundColor','white','HorizontalAlignment','Left'); + uicontrol('Style','Text','String','press enter to continue','Units','normalized','Position',[uiXOffset uiYOffset-uiXIncrement*(nPropertiesToChange+2) 0.3 0.05],'BackgroundColor','white','HorizontalAlignment','Left'); % uicontrol('Style','Text','String',inputTitleStr,'Units','normalized','Position',[0.0 uiYOffset 0.15 0.05],'BackgroundColor','white','HorizontalAlignment','Left'); pause diff --git a/@ciatah/computeMatchObjBtwnTrials.m b/@ciatah/computeMatchObjBtwnTrials.m index 34d398e..a0ff423 100644 --- a/@ciatah/computeMatchObjBtwnTrials.m +++ b/@ciatah/computeMatchObjBtwnTrials.m @@ -12,6 +12,8 @@ % 2021.06.18 [21:41:07] - added modelVarsFromFilesCheck() to check and load signals if user hasn't already. % 2021.08.10 [09:57:36] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package. % 2021.11.05 [14:13:36] - Added manual alignment option before automated for cases in which there are large shifts or rotations along with changes in bulk cells identified that might be hard for automated to align. + % 2023.05.08 [18:55:36] - Added the option to only keep cells within a specific region of the FOV, useful for when want to register cells in a single region and don't want other cells affecting alignment. + % 2023.05.09 [11:52:23] - Add support for normcorre and imregdemons for registration. % TODO % @@ -56,7 +58,55 @@ % usrIdxChoice = userDefaults; % end scnsize = get(0,'ScreenSize'); - userDefaults = {'1','5','0.4','','2','1','corr2','0.6','0.3','1e-6','0','1','filtered','0'}; + + % Check for prior user inputs + if isfield(obj.functionSettings,'computeMatchObjBtwnTrials') + uOpts = obj.functionSettings.computeMatchObjBtwnTrials; + userDefaults = { + uOpts.nCorrections; + uOpts.maxDistance; + uOpts.imageThreshold; + uOpts.trialToAlignUserOption; + uOpts.RegisTypeFinal; + uOpts.runImageCorr; + uOpts.imageCorrType; + uOpts.imageCorr; + uOpts.imageCorrThres; + uOpts.checkImageCorr; + uOpts.turboregZeroThres; + uOpts.runViewMatchObjBtwnSessions; + uOpts.imagesType; + uOpts.runManualAlign; + uOpts.runRemoveRegion; + uOpts.sizeThresMax; + uOpts.mcMethod; + uOpts.registrationFxn; + }; + userDefaults = cellfun(@num2str,userDefaults,'UniformOutput',false); + else + % userDefaults = {'1','5','0.4','','2','1','corr2','0.6','0.3','1e-6','0','1','filtered','0','0','','turboreg','transfturboreg'}; + userDefaults = { + '1'; + '5'; + '0.4'; + ''; + '2'; + '1'; + 'corr2'; + '0.6'; + '0.3'; + '1e-6'; + '0'; + '1'; + 'filtered'; + '0'; + '0'; + ''; + 'turboreg'; + 'transfturboreg'; + }; + end + usrIdxChoice = inputdlg({... 'Number of rounds to register images (integer)',... 'Distance threshold to match cells cross-session (in pixels)',... @@ -72,28 +122,44 @@ 'View full results after [viewMatchObjBtwnSessions] (1 = yes, 0 = no)',... 'Type of image to use for cross-session alignment? (filtered, raw)',... 'Manual cross-session alignment? (1 = yes, 0 = no)'... + 'Only keep cells in user defined region of FOV? (1 = yes, 0 = no)',... + 'Int: Maximum size of cell? (leave blank to skip)'... + 'Registration method? "turboreg", normcorre", or "imregdemons"'... + 'Registration function? "imtransform", "imwarp", "transfturboreg"'... },'Cross-session cell alignment options',1,... userDefaults); % options.frameList = [1:500]; s1 = 1; uOpts = struct; - nCorrections = str2num(usrIdxChoice{s1});s1=s1+1; - maxDistance = str2num(usrIdxChoice{s1});s1=s1+1; - imageThreshold = str2num(usrIdxChoice{s1});s1=s1+1; - trialToAlignUserOption = str2num(usrIdxChoice{s1});s1=s1+1; - RegisTypeFinal = str2num(usrIdxChoice{s1});s1=s1+1; - runImageCorr = str2num(usrIdxChoice{s1});s1=s1+1; - imageCorrType = usrIdxChoice{s1};s1=s1+1; - imageCorr = str2num(usrIdxChoice{s1});s1=s1+1; - imageCorrThres = str2num(usrIdxChoice{s1});s1=s1+1; - checkImageCorr = str2num(usrIdxChoice{s1});s1=s1+1; - turboregZeroThres = str2num(usrIdxChoice{s1});s1=s1+1; - runViewMatchObjBtwnSessions = str2num(usrIdxChoice{s1});s1=s1+1; - imagesType = usrIdxChoice{s1};s1=s1+1; + uOpts.nCorrections = str2num(usrIdxChoice{s1});s1=s1+1; + uOpts.maxDistance = str2num(usrIdxChoice{s1});s1=s1+1; + uOpts.imageThreshold = str2num(usrIdxChoice{s1});s1=s1+1; + uOpts.trialToAlignUserOption = str2num(usrIdxChoice{s1});s1=s1+1; + uOpts.RegisTypeFinal = str2num(usrIdxChoice{s1});s1=s1+1; + uOpts.runImageCorr = str2num(usrIdxChoice{s1});s1=s1+1; + uOpts.imageCorrType = usrIdxChoice{s1};s1=s1+1; + uOpts.imageCorr = str2num(usrIdxChoice{s1});s1=s1+1; + uOpts.imageCorrThres = str2num(usrIdxChoice{s1});s1=s1+1; + uOpts.checkImageCorr = str2num(usrIdxChoice{s1});s1=s1+1; + uOpts.turboregZeroThres = str2num(usrIdxChoice{s1});s1=s1+1; + uOpts.runViewMatchObjBtwnSessions = str2num(usrIdxChoice{s1});s1=s1+1; + uOpts.imagesType = usrIdxChoice{s1};s1=s1+1; uOpts.runManualAlign = str2num(usrIdxChoice{s1});s1=s1+1; + uOpts.runRemoveRegion = str2num(usrIdxChoice{s1});s1=s1+1; + uOpts.sizeThresMax = str2num(usrIdxChoice{s1});s1=s1+1; + uOpts.mcMethod = usrIdxChoice{s1};s1=s1+1; + uOpts.registrationFxn = usrIdxChoice{s1};s1=s1+1; + + % Save user settings to object + if isfield(obj.functionSettings,'computeMatchObjBtwnTrials') + % Do nothing + else + obj.functionSettings.computeMatchObjBtwnTrials = struct; + end + obj.functionSettings.computeMatchObjBtwnTrials = uOpts; - for thisSubjectStr=subjectList + for thisSubjectStr = subjectList try display(repmat('=',1,21)) thisSubjectStr = thisSubjectStr{1}; @@ -117,7 +183,7 @@ % obj.folderBaseSaveStr{obj.fileNum} % [rawSignalsTmp rawImagesTmp signalPeaks signalPeaksArray] = modelGetSignalsImages(obj,'returnType','raw'); - [rawSignalsTmp, rawImagesTmp, signalPeaks, signalPeaksArray] = modelGetSignalsImages(obj,'returnType',imagesType); + [rawSignalsTmp, rawImagesTmp, signalPeaks, signalPeaksArray] = modelGetSignalsImages(obj,'returnType',uOpts.imagesType); if ~isempty(rawSignalsTmp) display('adding to alignment...') rawSignals{end+1} = rawSignalsTmp; @@ -166,10 +232,36 @@ display(['rawSignals: ' num2str(size(rawSignals))]) display(['rawImages: ' num2str(size(rawImages))]) - if isempty(trialToAlignUserOption) + if isempty(uOpts.trialToAlignUserOption) trialToAlign = floor(quantile(1:length(validFoldersIdx),0.5)); else - trialToAlign = trialToAlignUserOption; + trialToAlign = uOpts.trialToAlignUserOption; + end + + % Automatically remove + if ~isempty(uOpts.sizeThresMax) + disp('Filtering cell inputs') + filterImageOptions = struct(... + 'minNumPixels', 10,... + 'maxNumPixels', 100,... + 'SNRthreshold', 1.45,... + 'minPerimeter', 5,... + 'maxPerimeter', 50,... + 'minSolidity', 0.8,... + 'maxSolidity', 1.0,... + 'minEquivDiameter', 3,... + 'maxEquivDiameter', 30,... + 'slopeRatioThreshold', 0.04... + ); + nSessions = length(rawImages); + for sessionNo = 1:nSessions + disp('===') + disp(sessionNo) + [~, ~, validAuto, imageSizes, imgFeatures] = ciapkg.image.filterImages(rawImages{sessionNo}, [],'featureList',{'Eccentricity','EquivDiameter','Area','Orientation','Perimeter','Solidity'},'options',filterImageOptions); + signalsToKeep = imageSizes=round(nSessions*fractionShow(matchingNumbers)); + % else + % fractionShow = pctSessionMatch; + % globalMatchFilter = nMatchGlobalIDs>=round(nSessions*fractionShow); + % end + globalMatchFilter = nMatchGlobalIDs>=round(nSessions*fractionShow(matchingNumbers)); + nGlobalIDsMatch = sum(globalMatchFilter); % globalIdxMax = max(globalIDs(globalMatchFilter,:),[],'all'); globalIdxMax = max(find(globalMatchFilter)); + + outputthisCellmap = []; for sessionNo = 1:nSessions try obj.fileNum = validFoldersIdx(sessionNo); @@ -514,6 +528,7 @@ % get % figure;plot(nMatchGlobalIDs==nSessions); keepIDIdx = globalIDs(globalMatchFilter,sessionNo); + maxNumGlobalIDs = size(globalIDs,1); % if matchingNumbers==1 % keepIDIdx = globalIDs(nMatchGlobalIDs==nSessions,sessionNo); % else @@ -531,14 +546,21 @@ % globalToSessionIDs{sessionNo}(keepIDIdx) globalToSessionIDsTmp(setdiff(1:length(globalToSessionIDsTmp),keepIDIdx)) = 1; globalToSessionIDsTmp(keepIDIdx) = globalToSessionIDsTmp(keepIDIdx)+10; - [groupedImagesRates] = groupImagesByColor(inputImages{sessionNo},globalToSessionIDsTmp); + [groupedImagesRates] = groupImagesByColor(inputImages{sessionNo},globalToSessionIDsTmp,'threshold',0.3); % [groupedImagesRates] = groupImagesByColor(inputImages{sessionNo}(keepIDIdx,:,:),globalToSessionIDs{sessionNo}(keepIDIdx)); % [groupedImagesRates] = groupImagesByColor(inputImages{sessionNo},globalToSessionIDs{sessionNo}); % size(inputImages{sessionNo}(,:,:)) % size(globalToSessionIDs{sessionNo}(nMatchGlobalIDs==nSessions)) thisCellmap = createObjMap(groupedImagesRates); thisCellmap(1,1) = 1; - thisCellmap(1,2) = globalIdxMax; + % thisCellmap(1,2) = globalIdxMax; + thisCellmap(1,2) = maxNumGlobalIDs+12; + + if isempty(outputthisCellmap) + outputthisCellmap = thisCellmap; + else + outputthisCellmap = cat(3,outputthisCellmap,thisCellmap); + end setCmapHere = @(nGlobalIDs) colormap([0 0 0; [1 1 1]*0.3; hsv(nGlobalIDs)]); [~, ~] = openFigure(sessionNo, ''); clf @@ -571,6 +593,12 @@ end end + saveDirPath = [obj.videoSaveDir filesep 'matchObjsGlobal\']; + mkdir(saveDirPath); + currentDateTimeStr = datestr(now,'yyyy-mm-dd-HHMMSS','local'); + % [output] = writeHDF5Data(single(outputthisCellmap),[saveDirPath thisSubjectStr '_cellmap_turboreg_main.tif']); + [output] = ciapkg.io.saveMatrixToFile(single(outputthisCellmap),[saveDirPath thisSubjectStr '_cellmap_turboreg_main_' currentDateTimeStr '.tif']); + saveVideoFile = [obj.picsSavePath filesep folderSaveName{matchingNumbers} 'Session' filesep thisSubjectStr '_matchedCells.avi']; display(['save video: ' saveVideoFile]) writerObj = VideoWriter(saveVideoFile,'Motion JPEG AVI'); diff --git a/@ciatah/viewMovie.m b/@ciatah/viewMovie.m index dc56474..53d4664 100644 --- a/@ciatah/viewMovie.m +++ b/@ciatah/viewMovie.m @@ -16,6 +16,7 @@ % 2021.06.20 [00:14:59] - Added support for simple and advanced settings. % 2021.08.10 [09:57:36] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package. % 2021.12.31 [18:59:24] - Updated suptitle to ciapkg.overloaded.suptitle. + % 2023.06.10 [15:48:01] - Fix readHDF5Subset function unable to be found in certain cases. % TODO % @@ -897,7 +898,7 @@ function equalizeMovieHistograms() xDim = hReadInfo.Dims(1); yDim = hReadInfo.Dims(2); % select the first frame from the dataset - thisFrame = readHDF5Subset(inputFilePath,[0 0 options.refCropFrame],[xDim yDim 1],'datasetName',options.datasetName); + thisFrame = ciapkg.hdf5.readHDF5Subset(inputFilePath,[0 0 options.refCropFrame],[xDim yDim 1],'datasetName',options.datasetName); elseif strcmp(ext,'.tif')|strcmp(ext,'.tiff') TifLink = Tiff(inputFilePath, 'r'); %Create the Tiff object thisFrame = TifLink.read();%Read in one picture to get the image size and data type diff --git a/@ciatah/viewObjmaps.m b/@ciatah/viewObjmaps.m index d1c467b..118d3cb 100644 --- a/@ciatah/viewObjmaps.m +++ b/@ciatah/viewObjmaps.m @@ -12,6 +12,7 @@ % 2020.05.07 [14:54:27] - Fix to deal with empty valid folders. % 2021.06.18 [21:41:07] - added modelVarsFromFilesCheck() to check and load signals if user hasn't already. % 2021.08.10 [09:57:36] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package. + % 2024.03.30 [14:13:35] - Update display of cell overlap % TODO % @@ -355,7 +356,7 @@ imagesc(thisCellmap) axis equal tight zoom on - colormap([0 0 0;customColormap([])]) + colormap([0 0 0;0 1 0;customColormap([],'nPoints',max(thisCellmap(:),[],'omitnan')-1)]); cbh = colorbar; ylabel(cbh,'# cells at that pixel location','FontSize',10); titleStr = sprintf('Cell overlap | %d/%d: %s | %s | %s\n %s | %d cells, %d total | Zoom enabled',thisFileNumIdx,nFilesToAnalyze,obj.folderBaseDisplayStr{obj.fileNum},strrep(foldername,'_','\_'),validType,obj.signalExtractionMethod,sum(valid==1),length(valid)); diff --git a/README.md b/README.md index 6bfa6e6..5c4b6ce 100644 --- a/README.md +++ b/README.md @@ -23,7 +23,7 @@ ## Full documentation at https://bahanonu.github.io/ciatah/. -- Documentation for spinal cord motion correction can be found at https://bahanonu.github.io/ciatah/pipeline_detailed_spinal. + Below are recordings and additional documents for users who want to learn more about calcium imaging analysis/experiments and the CIAtah pipeline. @@ -38,8 +38,11 @@ Below are recordings and additional documents for users who want to learn more a __GRINjector__ — A surgical device to help with implanting gradient-refractive index (GRIN) lens probes into the brain or other regions: https://github.com/bahanonu/GRINjector. -__Upcoming motion correction methods__ — Methods for motion correction of spinal imaging data using feature identification (e.g. with DeepLabCut), control point registration, and other methods. Additional updates on integration into CIAtah in the future. +### Spinal cord imaging +__New motion correction methods__ — Methods for motion correction of spinal imaging data using feature identification (e.g. with DeepLabCut), control point registration, and other methods. Additional updates on integration into CIAtah in the future. - Preprint: Ahanonu and Crowther, _et al_. (2023). _Long-term optical imaging of the spinal cord in awake, behaving animals_. bioRxiv (https://www.biorxiv.org/content/10.1101/2023.05.22.541477v1.full). +- Documentation for spinal cord motion correction can be found at https://bahanonu.github.io/ciatah/pipeline_detailed_spinal. +- Other documentation for spinal cord imaging at https://github.com/basbaumlab/spinal_cord_imaging. @@ -55,7 +58,7 @@ Below are recordings and additional documents for users who want to learn more a - [CIAtah features](#ciatah-features) - [CIAtah example features](#ciatah-example-features) - [Quick start guide](#quick-start-guide) -- [Quick start (command-line)](#quick-start-command-line) +- [Quick start (command-line)](#quick-start-command-line-or-gui-less-batch-analysis) - [`CIAtah` main GUI notes](#ciatah-main-gui-notes) - [Acknowledgments](#acknowledgments) - [References](#references) @@ -215,13 +218,14 @@ obj % then hit enter, no semicolon! ### Visualize movies quickly using read from disk +#### Method #1: CIAtah `playMovie` Users can quickly visualize movies in any of the supported formats (HDF5, NWB, AVI, TIFF, ISXD, etc.) using the `playMovie` function. This will read directly from disk, allowing users to scroll through frames to visually check movies before or after processing. Users can run via the command-line: ```MATLAB % Use the absolute path to the movie file or a valid relative path. -ciapkg.api.playMovie('ABSOLUTE\PATH\TO\MOVIE'); +ciapkg.api.playMovie('ABSOLUTE_PATH_TO_MOVIE\MOVIE_NAME.EXTENSION'); ``` When using HDF5 files, check the dataset name containing movie with `h5disp` then input the full dataset name (e.g. below is for a standard NWB-formatted HDF5 file): @@ -229,7 +233,8 @@ When using HDF5 files, check the dataset name containing movie with `h5disp` the ciapkg.api.playMovie('ABSOLUTE\PATH\TO\MOVIE','inputDatasetName','/acquisition/TwoPhotonSeries/data'); ``` -Alternatively, using the `ciatah` GUI class, users can select loaded folders and change the regular expression to match the name of the files in the movie, both for the raw data and for any processed movies in the folder. See below: +#### Method #2: CIAtah GUI +Using the `ciatah` GUI class, users can select loaded folders and change the regular expression to match the name of the files in the movie, both for the raw data and for any processed movies in the folder. See below:

@@ -238,6 +243,13 @@ Alternatively, using the `ciatah` GUI class, users can select loaded folders and

+#### Method #3: Fiji + N5 +Install N5 library to enable N5 viewer in Fiji: go to `Help->Update...` then under `Manage Update Sites` scroll down to select `N5` (https://sites.imagej.net/N5/) then `Apply and Close`). View a large HDF5 movie by running `File->Import->"HDF5/N5/Zarr/OME-NGFF"`. In the resulting dialog box, select `Browse` along with `Open as virtual` on the bottom left. After a moment different a tree should appear listing different datasets + +- https://imagej.net/libs/n5 +- https://github.com/saalfeldlab/n5 +- https://github.com/saalfeldlab/n5-ij +- https://github.com/saalfeldlab/n5-viewer ## Quick start (command line or GUI-less batch analysis) @@ -252,6 +264,10 @@ import ciapkg.api.* % import CIAtah functions in ciapkg package API. ``` +## `CIAtah` main GUI notes + +See more about the CIAtah GUI at https://bahanonu.github.io/ciatah/install/#ciatah-main-gui-notes. + *** ## Acknowledgments @@ -327,7 +343,7 @@ Users with versions of MATLAB earlier than R2019b can download `CIAtah` version ## License -Copyright (C) 2013-2023 Biafra Ahanonu +Copyright (C) 2013-2024 Biafra Ahanonu This project is licensed under the terms of the MIT license. See LICENSE file for details. diff --git a/docs/docs/help_spatial_filtering.md b/docs/docs/help_spatial_filtering.md index 801baa2..90dd642 100644 --- a/docs/docs/help_spatial_filtering.md +++ b/docs/docs/help_spatial_filtering.md @@ -4,14 +4,25 @@ This page documents different features and functions in the {{ site.name }} repo ![2014_04_01_p203_m19_check01_spatialFilter](https://user-images.githubusercontent.com/5241605/111711340-af792480-8808-11eb-85d1-8c76479995bb.gif) +## Why conduct spatial filtering? + +Spatial filtering can have a large impact on the resulting cell activity traces extracted from the movies and can lead to erroneous conclusions if not properly applied during pre-processing. + +For example, below are the correlations between all cell-extraction outputs from PCA-ICA, ROI back-application of ICA filters, and CNMF-e on a miniature microscope one-photon movie. As can be seen, especially in the case of ROI analysis, the correlation between the activity traces is rendered artificially high due to the correlated background noise. This is greatly reduced in many instances after proper spatial filtering. + +![image](https://user-images.githubusercontent.com/5241605/111710928-e864c980-8807-11eb-9d29-81341290d108.png) + + + + -## Filtering movies with `normalizeMovie` +## Filtering movies with `ciapkg.movie_processing.normalizeMovie` -Users can quickly filter movies using the `normalizeMovie` function. See below for usage. +Users can quickly filter movies using the `ciapkg.movie_processing.normalizeMovie` function. See below for usage. ```Matlab % Load the movie (be in TIF, HDF5, AVI, etc.). Change HDF5 input dataset name as needed. -inputMovie = loadMovieList('pathToMovie','inputDatasetName','/1'); +inputMovie = ciapkg.io.loadMovieList('pathToMovie','inputDatasetName','/1'); % Set options for dividing movie by lowpass version of the movie options.freqLow = 1; @@ -31,57 +42,77 @@ options.bandpassMask = 'gaussian'; options.showImages = 0; % Set to 1 if you want to view % Run analysis -inputMovie = normalizeMovie(single(inputMovie),'options',options); +inputMovie = ciapkg.movie_processing.normalizeMovie(single(inputMovie),'options',options); ``` -If users set `options.showImages = 0;`, then `normalizeMovie` will update a figure containing both real and frequency space before and after the filter has been applied along with an example of the filter in frequency space. This allows users to get a sense of what their filter is doing. +If users set `options.showImages = 0;`, then `normalizeMovie` will update a figure containing both real and frequency space before and after the filter has been applied along with an example of the filter in frequency space. This allows users to get a sense of what their filter is doing. See below for examples. -![image](https://user-images.githubusercontent.com/5241605/111677991-049f4100-87dd-11eb-9bb6-1ba46894ea70.png) + -## Why conduct spatial filtering? +### FFT bandpass filtering +Bandpass filtering where only `red` frequencies in `filter` image (FFT of bottom left input image) are kept producing an image as in `fft image`. +![image](https://github.com/bahanonu/ciatah/assets/5241605/78fda2d7-a041-4622-a9b6-3b667f24bc49) -Spatial filtering can have a large impact on the resulting cell activity traces extracted from the movies. +### Divide by lowpass filtering +Another method is `lowpassFFTDivisive`, which involves dividing the image by a lowpass version of itself. In the below example, the `filter` image shows that only low frequencies will be kept. This will produce an image as in `fft image` that when divided or subtracted from the `input image` will produce `difference` image. +![image](https://github.com/bahanonu/ciatah/assets/5241605/cde8556b-e009-4e92-b28f-b277f331a9a5) -For example, below are the correlations between all cell-extraction outputs from PCA-ICA, ROI back-application of ICA filters, and CNMF-e on a miniature microscope one-photon movie. As can be seen, especially in the case of ROI analysis, the correlation between the activity traces is rendered artificially high due to the correlated background noise. This is greatly reduced in many instances after proper spatial filtering. -![image](https://user-images.githubusercontent.com/5241605/111710928-e864c980-8807-11eb-9d29-81341290d108.png) - - - - ## Images from unit test + ### Main filtering functions. Below is a screen grab from a random frame using all the filtering functions. A nice way to quickly see the many differences between each functions filtering. -![image](https://user-images.githubusercontent.com/5241605/32477562-18b7d9d4-c334-11e7-988f-accdf99a22f2.png) + + + +![image](https://github.com/bahanonu/ciatah/assets/5241605/51356dea-deb6-4bbe-904b-2a79577ffcdf) + ### Test function filtering -This function will take a movie and This is currently only for the Matlab fft, but I'll see about expanding to others. +This function will take a movie and conduct multiple spatial filtering operations on it then display for the user. + ```Matlab -unitNormalizeMovie; +ciapkg.unit.unitNormalizeMovie; ``` -![image](https://user-images.githubusercontent.com/5241605/32477571-26620546-c334-11e7-8ce0-aa5269fcb5f3.png) +After running that function, below is an example movie from a prefrontal cortex animal (miniature microscope, GCaMP) showing the difference in results with different spatial filtering. + +![image](https://github.com/bahanonu/ciatah/assets/5241605/a09e959e-da75-49c3-8c89-7ee1e42ac004) ### Matlab test function -* I've also added the ability to test the parameter space of the Matlab fft, use the below command. + +I've also added the ability to test the parameter space of the Matlab fft, use the below command. + ```Matlab -testMovieFFT = normalizeMovie(testMovie,'normalizationType','matlabFFT_test','secondaryNormalizationType','lowpassFFTDivisive','bandpassMask','gaussian','bandpassType','lowpass'); +testMovieFFT = ciapkg.movie_processing.normalizeMovie(testMovie,'normalizationType','matlabFFT_test','secondaryNormalizationType','lowpassFFTDivisive','bandpassMask','gaussian','bandpassType','lowpass'); ``` -* Should get a movie output similar to the below, where there is the original movie, the FFT movie, the original/FFT movie, and the dfof of original/FFT movie. + +Should get a movie output similar to the below, where there is the original movie, the FFT movie, the original/FFT movie, and the dfof of original/FFT movie. + ![image](https://cloud.githubusercontent.com/assets/5241605/11490967/559152e2-9792-11e5-839b-a93811df70ce.png) +This can also be expanded to look at the effects of different spatial frequency filters on the resulting output, as indicated below. + +![image](https://user-images.githubusercontent.com/5241605/32477571-26620546-c334-11e7-8ce0-aa5269fcb5f3.png) + ### Matlab test function movie output Similar to above, showing results when using `lowpassFFTDivisive` normalization (using the `matlab divide by lowpass before registering` option in `modelPreprocessMovie` and `viewMovieRegistrationTest` functions) with `freqLow = 0` and `freqHigh` set to `1`, `4`, and `20`. This corresponds to removing increasingly smaller features from the movie. -![2014_04_01_p203_m19_check01_fft_example-3](https://user-images.githubusercontent.com/5241605/71422606-aec30400-2640-11ea-8ffb-41cdeea771c1.gif) + + + + + +![image](https://github.com/bahanonu/ciatah/assets/5241605/2534544b-7519-4d17-b4f5-088ee5582a87) + ### ImageJ test function To test the ImageJ FFT and determine the best parameters for a given set of movies, run the following function on a test movie matrix: ```Matlab -inputMovieTest = normalizeMovie(inputMovie,'normalizationType','imagejFFT_test'); +inputMovieTest = ciapkg.movie_processing.normalizeMovie(inputMovie,'normalizationType','imagejFFT_test'); ``` The output should look like the below: @@ -98,8 +129,11 @@ If the spatial filter is not properly configured then dark halos will appear aro ![image](https://cloud.githubusercontent.com/assets/5241605/11329062/1232a886-914b-11e5-9cca-85ec6162319b.png) -https://github.com/schnitzer-lab/miniscope_analysis/pull/30 -* FYI, for 4x downsampled movies, `highFreq` parameter of 4 (which corresponds to a `fspecial` gaussian with std of 4) produces the closest results to ImageJ `Process->FFT->Bandpass Filter...` with inputs of `filter_large=10000 filter_small=80 suppress=None tolerance=5` (the current default in `normalizeMovie`). + + +* FYI, for 4x downsampled movies, `highFreq` parameter of 4 (which corresponds to a `fspecial` gaussian with std of 4) produces the closest results to ImageJ `Process->FFT->Bandpass Filter...` with inputs of `filter_large=10000 filter_small=80 suppress=None tolerance=5` (the current default in `ciapkg.movie_processing.normalizeMovie`). + +### Comparison of MATLAB and ImageJ FFT-based spatial filtering * Example frame from ImageJ and Matlab FFTs. ![image](https://cloud.githubusercontent.com/assets/5241605/11519196/fbcd561e-984c-11e5-95b9-3a23085c2e44.png) @@ -110,11 +144,14 @@ https://github.com/schnitzer-lab/miniscope_analysis/pull/30 * This matches the filter that ImageJ says it uses, which is fairly close to the Matlab filter. ![image](https://cloud.githubusercontent.com/assets/5241605/11519329/22b3af8e-984e-11e5-9379-5b5d4092cee0.png) -__Example video: 2015_11_25_p384_m610_openfield01__ +__Example video: Basolateral amygdala miniature microscope imaging in open field__ + + * Below is an example comparison using the following Matlab commands to produce the filtered inputs: + ```Matlab -testMovieFFT = normalizeMovie(testMovie,'normalizationType','lowpassFFTDivisive','freqHigh',7); -testMovieFFTImageJ = normalizeMovie(testMovie,'normalizationType','imagejFFT'); +testMovieFFT = ciapkg.movie_processing.normalizeMovie(testMovie,'normalizationType','lowpassFFTDivisive','freqHigh',7); +testMovieFFTImageJ = ciapkg.movie_processing.normalizeMovie(testMovie,'normalizationType','imagejFFT'); diffMovie = testMovieFFT-testMovieFFTImageJ ; ``` diff --git a/docs/docs/pipeline_detailed_signal_extraction_validation.md b/docs/docs/pipeline_detailed_signal_extraction_validation.md index 4d69977..8206c48 100644 --- a/docs/docs/pipeline_detailed_signal_extraction_validation.md +++ b/docs/docs/pipeline_detailed_signal_extraction_validation.md @@ -10,4 +10,7 @@ Below is an example, with black outlines indicating location of cell extraction ![2014_04_01_p203_m19_check01_raw_viewCellExtractionOnMovie_ezgif-4-57913bcfdf3f_2](https://user-images.githubusercontent.com/5241605/59560798-50033a80-8fcc-11e9-8228-f9a3d83ca591.gif) +![2014_04_01_p203_m19_check01_overlay2-2](https://github.com/bahanonu/ciatah/assets/5241605/ee5c104a-7b4c-4fe5-b2da-63ddac8955fb) + + \ No newline at end of file diff --git a/docs/docs/pipeline_detailed_signal_sorting_manual.md b/docs/docs/pipeline_detailed_signal_sorting_manual.md index 539e506..6c21ac6 100644 --- a/docs/docs/pipeline_detailed_signal_sorting_manual.md +++ b/docs/docs/pipeline_detailed_signal_sorting_manual.md @@ -9,17 +9,46 @@ title: Sorting cell extraction outputs.

- ciapkgMovie + ciapkgMovie

Outputs from most common cell-extraction algorithms like PCA-ICA, CNMF, etc. contain signal sources that are not cells and thus must be manually removed from the output. The repository contains a GUI for sorting cells from not cells. GUI also contains a shortcut menu that users can access by right-clicking or selecting the top-left menu. -Below users can see a list of options that are given before running the code, those highlighted in green +## Resources on manual identification + +![image](img/Ahanonu_Kitch_manualSort01.png) + +The above figure gives an overview of the CIAtah manual sorting GUI along with examples of candidate cells that are accepted or rejected based on a variety of criteria from several cell extraction algorithms (CELLMax, PCA-ICA, and CNMF). We have discussed manual sorting previously, see the below resources: + +- `3.15.1 Manual Neuron Identification` in our miniature microscope book chapter contains a guide on manual sorting: https://link.springer.com/protocol/10.1007/978-1-0716-2039-7_13#Sec20. +- `Fig. 7: Calcium imaging analysis of nociceptive ensemble.` contains example accepted and rejected cells: https://link.springer.com/protocol/10.1007/978-1-0716-2039-7_13/figures/7. + +Below are several potential criteria to use for accepting or rejecting candidate cells output by a cell extraction algorithm: + +- Filter shape—e.g., cell-like depending on if using one- or two-photon imaging). +- The event triggered movie activity—e.g., whether it conformed to prior expectation of one-photon neuron morphology and fluorescent indicator activity. __Note__ This criteria is critical, as some methods output candidate cells whose cell shape and activity trace look like a cell, but when the movie is checked can see that it is not a cell. +- Location within the imaging field of view—e.g., not within a blood vessel. +- The shape of the transient having characteristic fluorescent indicator dynamics, this will depending on the indicator being used, e.g. GCaMP will have a different expected waveform than other indicators. +- Whether cell is a duplicate cell, e.g. some algorithms will "split" a cell into multiple candidate cells. This can be handled by re-running the algorithm with improved parameters, rejected the lower SNR (or otherwise poorer quality) cell, or accepting both cells then conducting a merging operation later (and re-running the cell trace extraction portion of the algorithm if that feature is available). + +## CIAtah manual sorting GUI + +Below users can see a list of options that are given before running the code. Options highlighted in green are those that are changed from the default settings. ![image](https://user-images.githubusercontent.com/5241605/49845107-43322f80-fd7a-11e8-96b9-3f870d4b9009.png) +### Loading in prior manually sorted data + +Decisions during manual sorting are stored in the `private/tmp` folder within the root CIAtah directory (find with `ciapkg.getDir`). Alternatively, previously manually sorted outputs can be re-sorted if new selection criteria are desired. When loading the `computeManualSortSignals` GUI, select one of the two options below in the `Use CIAtah auto classifications?` setting. + +- `Start with TEMP manually chosen classifications (e.g. backups)` - this option will open up a GUI into `private/tmp` and request users select a MAT-file containing the most recent decisions that were being manually sorted. +- `Start with FINISHED manually chosen classifications` - will automatically load already saved manual decisions located in the same folder as the cell extraction outputs. + +![image](img/manualSort_reload01.png) + + ## GUI usage on large imaging datasets - To manually sort on large movies that will not fit into RAM, select the below options (highlighted in green). This will load only chunks of the movie asynchronously into the GUI as you sort cell extraction outputs.