From ff428e9943d4d33f8db1e8bdcee19b2516379979 Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Mon, 25 Jul 2022 14:45:53 +0100 Subject: [PATCH 01/43] funciton name correction --- Browsing Functions/AtlasTransformBrowser.m | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/Browsing Functions/AtlasTransformBrowser.m b/Browsing Functions/AtlasTransformBrowser.m index 5cf501c..a93b803 100644 --- a/Browsing Functions/AtlasTransformBrowser.m +++ b/Browsing Functions/AtlasTransformBrowser.m @@ -1,4 +1,4 @@ -function f = allenAtlasBrowser(f, templateVolume, annotationVolume, structureTree, slice_figure, save_location, save_suffix, plane) +function f = AtlasTransformBrowser(f, templateVolume, annotationVolume, structureTree, slice_figure, save_location, save_suffix, plane) % ------------------------------------------------ % Browser for the allen atlas ccf data in matlab. % ------------------------------------------------ @@ -359,7 +359,7 @@ function hotkeyFcn(f, slice_figure, keydata, allData, save_location, save_suffix template_point = 1; template_points_shown = 0; updateBoundaries(f,ud, allData); ud = get(f, 'UserData'); end - else; disp('transform point mode OFF'); + else; disp('transform point mode OFF'); end % a -- toggle viewing of annotation boundaries case 'a' @@ -432,7 +432,7 @@ function hotkeyFcn(f, slice_figure, keydata, allData, save_location, save_suffix ud.slice_at_shift_start = ud_slice.slice_num; ud.slice_shift = 0; end - disp('switch scroll mode -- scroll along slice images') + disp('switch scroll mode -- scroll along slice images') % h -- toggle viewing of histology / histology overlay case 'h' disp(' '); @@ -732,7 +732,7 @@ function updateSlice(f, evt, allData, slice_figure, save_location) % scroll through slices (left arrow pressed) elseif ud.scrollMode == 3 - set(ud.pointHands_for_transform(:), 'Visible', 'off'); + set(ud.pointHands_for_transform(:), 'Visible', 'off'); ud.showOverlay = 0; delete(ud.overlayAx); ud.overlayAx = []; ud_slice = get(slice_figure, 'UserData'); @@ -776,7 +776,7 @@ function updateSlice(f, evt, allData, slice_figure, save_location) ud.current_pointList_for_transform = transform_data.transform_points{1}; ud_slice.pointList = transform_data.transform_points{2}; else - ud_slice.pointList = []; + ud_slice.pointList = []; end set(slice_figure, 'UserData', ud_slice); From 5d6e40916bcf4238f4385c60369faaf662c49d75 Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Mon, 25 Jul 2022 14:49:32 +0100 Subject: [PATCH 02/43] load and show saved slice points --- Browsing Functions/sliceBrowser.m | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/Browsing Functions/sliceBrowser.m b/Browsing Functions/sliceBrowser.m index e190b6e..9711a15 100644 --- a/Browsing Functions/sliceBrowser.m +++ b/Browsing Functions/sliceBrowser.m @@ -45,6 +45,9 @@ function sliceClickCallback(im, keydata) f = get(get(im, 'Parent'), 'Parent'); ud = get(f, 'UserData'); +if isempty(ud.pointHands) + ud.pointHands = gobjects(0); +end if ud.getPoint clickX = round(keydata.IntersectionPoint(1)); @@ -127,7 +130,13 @@ function SliceAtlasHotkeyFcn(fig, keydata, f) file_transformations = fullfile(ud.processed_images_folder, 'transformations\\' ,... [processed_image_name(1:end-4) '_transform_data.mat']); - set(ud.pointHands(:), 'Visible', 'off'); + try + delete(findobj(ud.im.Parent, 'Type','line','color',[0 0.5 0],'Marker','o')); % delete all the circles + end + + ud.pointHands = gobjects(0); + % set(ud.pointHands(:), 'Visible', 'off'); %TODO WHY hiding? Maybe you + % want to delete/initialize them? ud is from fig.UserData ud.pointList = []; if exist(file_transformations,'file') @@ -136,6 +145,9 @@ function SliceAtlasHotkeyFcn(fig, keydata, f) transform_data = transform_data.save_transform; if ~isempty(transform_data.transform_points{2}) ud.pointList = transform_data.transform_points{2}; + for i = 1:size(ud.pointList,1) + ud.pointHands(end+1) = plot(ud.sliceAx, ud.pointList(i,1), ud.ref_size(1) - ud.pointList(i,2), 'ro', 'color', [0 .5 0],'linewidth',2,'markers',4); + end title_ending = ' (transform points loaded)'; end end From 0f801272106f69aa4e88877d70022c5cb2f68d51 Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Mon, 25 Jul 2022 14:56:15 +0100 Subject: [PATCH 03/43] show scroll mode as XLabel --- Browsing Functions/AtlasTransformBrowser.m | 37 +++++++++++++++++++--- 1 file changed, 33 insertions(+), 4 deletions(-) diff --git a/Browsing Functions/AtlasTransformBrowser.m b/Browsing Functions/AtlasTransformBrowser.m index 5cf501c..ad41e47 100644 --- a/Browsing Functions/AtlasTransformBrowser.m +++ b/Browsing Functions/AtlasTransformBrowser.m @@ -79,6 +79,15 @@ set(f, 'WindowScrollWheelFcn', @(src,evt)updateSlice(f, evt, allData, slice_figure, save_location)) set(f, 'WindowButtonMotionFcn',@(f,k)fh_wbmfcn(f, allData, slice_figure, save_location)); % Set the motion detector. +ud.atlasAx.XLabel.Visible = 'on'; +if strcmp(ud.plane,'coronal') + xlabel(ud.atlasAx, 'Scroll slices along A/P axis') +elseif strcmp(ud.plane,'sagittal') + xlabel(ud.atlasAx, 'Scroll slices along M/L axis') +elseif strcmp(ud.plane,'transverse') + xlabel(ud.atlasAx, 'Scroll slices along D/V axis') +end + % display user controls in the console function display_controls(plane) fprintf(1, '\n Controls: \n'); @@ -359,7 +368,7 @@ function hotkeyFcn(f, slice_figure, keydata, allData, save_location, save_suffix template_point = 1; template_points_shown = 0; updateBoundaries(f,ud, allData); ud = get(f, 'UserData'); end - else; disp('transform point mode OFF'); + else; disp('transform point mode OFF'); end % a -- toggle viewing of annotation boundaries case 'a' @@ -402,29 +411,47 @@ function hotkeyFcn(f, slice_figure, keydata, allData, save_location, save_suffix ud.scrollMode = 1; if strcmp(ud.plane,'coronal') disp('switch scroll mode -- tilt D/V') + xlabel(ud.atlasAx, 'Scroll to tilt D/V') + elseif strcmp(ud.plane,'sagittal') disp('switch scroll mode -- tilt D/V') + xlabel(ud.atlasAx, 'Scroll to tilt D/V') + elseif strcmp(ud.plane,'transverse') disp('switch scroll mode -- tilt M/L') + xlabel(ud.atlasAx, 'Scroll to tilt M/L') + end case 'rightarrow' % scroll angles along A/P axis ud.scrollMode = 2; if strcmp(ud.plane,'coronal') disp('switch scroll mode -- tilt M/L') + xlabel(ud.atlasAx, 'Scroll to tilt M/L') + elseif strcmp(ud.plane,'sagittal') disp('switch scroll mode -- tilt A/P') + xlabel(ud.atlasAx, 'Scroll to tilt A/P') + elseif strcmp(ud.plane,'transverse') disp('switch scroll mode -- tilt A/P') + xlabel(ud.atlasAx, 'Scroll to tilt A/P') + end case 'downarrow' % scroll along A/P axis ud.scrollMode = 0; if strcmp(ud.plane,'coronal') disp('switch scroll mode -- scroll along A/P axis') + xlabel(ud.atlasAx, 'Scroll slices along A/P axis') + elseif strcmp(ud.plane,'sagittal') disp('switch scroll mode -- scroll along M/L axis') + xlabel(ud.atlasAx, 'Scroll slices along M/L axis') + elseif strcmp(ud.plane,'transverse') disp('switch scroll mode -- scroll along D/V axis') + xlabel(ud.atlasAx, 'Scroll slices along D/V axis') + end case 'leftarrow' % scroll along A/P axis ud.scrollMode = 3; @@ -432,7 +459,9 @@ function hotkeyFcn(f, slice_figure, keydata, allData, save_location, save_suffix ud.slice_at_shift_start = ud_slice.slice_num; ud.slice_shift = 0; end - disp('switch scroll mode -- scroll along slice images') + disp('switch scroll mode -- scroll along slice images') + xlabel(ud.atlasAx, 'Scroll along slice images') + % h -- toggle viewing of histology / histology overlay case 'h' disp(' '); @@ -732,7 +761,7 @@ function updateSlice(f, evt, allData, slice_figure, save_location) % scroll through slices (left arrow pressed) elseif ud.scrollMode == 3 - set(ud.pointHands_for_transform(:), 'Visible', 'off'); + set(ud.pointHands_for_transform(:), 'Visible', 'off'); ud.showOverlay = 0; delete(ud.overlayAx); ud.overlayAx = []; ud_slice = get(slice_figure, 'UserData'); @@ -776,7 +805,7 @@ function updateSlice(f, evt, allData, slice_figure, save_location) ud.current_pointList_for_transform = transform_data.transform_points{1}; ud_slice.pointList = transform_data.transform_points{2}; else - ud_slice.pointList = []; + ud_slice.pointList = []; end set(slice_figure, 'UserData', ud_slice); From 575599bd12b75ad3aaa0b20241b7fef6a41bcf2c Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Mon, 25 Jul 2022 15:01:59 +0100 Subject: [PATCH 04/43] use graphic objects --- Browsing Functions/AtlasTransformBrowser.m | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Browsing Functions/AtlasTransformBrowser.m b/Browsing Functions/AtlasTransformBrowser.m index 5cf501c..eca9bc4 100644 --- a/Browsing Functions/AtlasTransformBrowser.m +++ b/Browsing Functions/AtlasTransformBrowser.m @@ -36,7 +36,8 @@ ud.probe_view_mode = false; ud.currentProbe = 0; ud.ProbeColors = [1 1 1; 1 .75 0; .3 1 1; .4 .6 .2; 1 .35 .65; .7 .7 1; .65 .4 .25; .7 .95 .3; .7 0 0; .5 0 .6; 1 .6 0]; ud.ProbeColor = {'white','gold','turquoise','fern','bubble gum','overcast sky','rawhide', 'green apple','red','purple','orange'}; -ud.getPoint_for_transform =false; ud.pointList_for_transform = zeros(0,2); ud.pointHands_for_transform = []; +ud.getPoint_for_transform =false; ud.pointList_for_transform = zeros(0,2); +ud.pointHands_for_transform = gobjects(0); ud.current_pointList_for_transform = zeros(0,2); ud.curr_slice_num = 1; ud.clicked = false; ud.showAtlas = false; From e9e1f10476cd90521f4ed4680dd13dd1c6dc9a2d Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Mon, 25 Jul 2022 15:04:26 +0100 Subject: [PATCH 05/43] load and show saved atlas points --- Browsing Functions/AtlasTransformBrowser.m | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/Browsing Functions/AtlasTransformBrowser.m b/Browsing Functions/AtlasTransformBrowser.m index eca9bc4..22e8050 100644 --- a/Browsing Functions/AtlasTransformBrowser.m +++ b/Browsing Functions/AtlasTransformBrowser.m @@ -775,6 +775,16 @@ function updateSlice(f, evt, allData, slice_figure, save_location) if ~isempty(transform_data.transform_points{1}) && ~isempty(transform_data.transform_points{2}) ud.current_pointList_for_transform = transform_data.transform_points{1}; + delete(ud.pointHands_for_transform);% delete them %TODO not working + ud.pointHands_for_transform = gobjects(0);% initialize + for i = 1:size(ud.current_pointList_for_transform,1) + ud.pointHands_for_transform(i) = plot(ud.atlasAx, ... + ud.current_pointList_for_transform(i,1), ... + ud.current_pointList_for_transform(i,2), ... + 'ro', 'color', [.7 .3 .3],'LineWidth',2,'markers',4); + end + set(ud.pointHands_for_transform(end),'color', [0 .9 0]); + ud_slice.pointList = transform_data.transform_points{2}; else ud_slice.pointList = []; From 2d8497dde2d2f20aced4bf72c7e7e6dcd8b10d6c Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Mon, 25 Jul 2022 15:19:42 +0100 Subject: [PATCH 06/43] chabge ud.pointHands_for_transform to grahic objects error handling --- Browsing Functions/AtlasTransformBrowser.m | 44 +++++++++++++++++----- 1 file changed, 35 insertions(+), 9 deletions(-) diff --git a/Browsing Functions/AtlasTransformBrowser.m b/Browsing Functions/AtlasTransformBrowser.m index 22e8050..efd1992 100644 --- a/Browsing Functions/AtlasTransformBrowser.m +++ b/Browsing Functions/AtlasTransformBrowser.m @@ -355,12 +355,27 @@ function hotkeyFcn(f, slice_figure, keydata, allData, save_location, save_suffix if ~size(ud.current_pointList_for_transform,1)% ) (ud.curr_slice_num ~= (ud.slice_at_shift_start+ud.slice_shift) || && ~ud.loaded ud.curr_slice_num = ud.slice_at_shift_start+ud.slice_shift; %ud_slice.slice_num; ud.current_pointList_for_transform = zeros(0,2); + try set(ud.pointHands_for_transform(:), 'Visible', 'off'); + catch + end num_hist_points = size(ud_slice.pointList,1); template_point = 1; template_points_shown = 0; updateBoundaries(f,ud, allData); ud = get(f, 'UserData'); + else + try + set(ud.pointHands_for_transform(:), 'Visible', 'on'); + catch + end + end + else; disp('transform point mode OFF'); + % hide transform points + + try + set(ud.pointHands_for_transform(:), 'Visible', 'off'); + catch end - else; disp('transform point mode OFF'); + end % a -- toggle viewing of annotation boundaries case 'a' @@ -433,7 +448,7 @@ function hotkeyFcn(f, slice_figure, keydata, allData, save_location, save_suffix ud.slice_at_shift_start = ud_slice.slice_num; ud.slice_shift = 0; end - disp('switch scroll mode -- scroll along slice images') + disp('switch scroll mode -- scroll along slice images') % h -- toggle viewing of histology / histology overlay case 'h' disp(' '); @@ -631,13 +646,13 @@ function hotkeyFcn(f, slice_figure, keydata, allData, save_location, save_suffix % Try to delete only the most recent point ud.current_pointList_for_transform = ud.current_pointList_for_transform(1:end-1,:); - ud_slice.pointList = ud_slice.pointList(1:end-1,:); + ud_slice.pointList = ud_slice.pointList(1:end-1,:); %TODO this should be accompanied by deleting circles on Slive Viewer set(slice_figure, 'UserData', ud_slice); - if ud.pointHands_for_transform + if ~isempty(ud.pointHands_for_transform) % remove circle for most revent point - set(ud.pointHands_for_transform(end), 'Visible', 'off'); + set(ud.pointHands_for_transform(end), 'Visible', 'off');%TODO why hiding rather than deleting??? ud.pointHands_for_transform = ud.pointHands_for_transform(1:end-1); - if ud.pointHands_for_transform + if ~isempty(ud.pointHands_for_transform) % recolor points set(ud.pointHands_for_transform(end), 'color', [0 .9 0]); end @@ -733,7 +748,10 @@ function updateSlice(f, evt, allData, slice_figure, save_location) % scroll through slices (left arrow pressed) elseif ud.scrollMode == 3 - set(ud.pointHands_for_transform(:), 'Visible', 'off'); + try + set(ud.pointHands_for_transform(:), 'Visible', 'off'); + catch + end ud.showOverlay = 0; delete(ud.overlayAx); ud.overlayAx = []; ud_slice = get(slice_figure, 'UserData'); @@ -787,7 +805,7 @@ function updateSlice(f, evt, allData, slice_figure, save_location) ud_slice.pointList = transform_data.transform_points{2}; else - ud_slice.pointList = []; + ud_slice.pointList = []; end set(slice_figure, 'UserData', ud_slice); @@ -811,6 +829,8 @@ function updateSlice(f, evt, allData, slice_figure, save_location) % if no transform, just show reference ud.histology_overlay = 0; ud.current_pointList_for_transform = zeros(0,2); + delete(ud.pointHands_for_transform);% delete them %TODO not working + ud.pointHands_for_transform = gobjects(0);% initialize set(ud.im, 'CData', ud.ref); ud.curr_im = ud.ref; set(f, 'UserData', ud); set(ud.text,'Visible','off'); @@ -845,7 +865,10 @@ function updateSlice(f, evt, allData, slice_figure, save_location) updateOverlay(f, allData, ann, slice_figure, save_location); end ud.ref = uint8(reference_slice); - set(ud.pointHands_for_transform(:), 'Visible', 'off'); +% try +% set(ud.pointHands_for_transform(:), 'Visible', 'off'); +% catch +% end ud.offset_map = zeros(ud.ref_size); set(f, 'UserData', ud); @@ -1071,7 +1094,10 @@ function atlasClickCallback(im, keydata, slice_figure, save_location) end ud.pointList_for_transform(end+1, :) = [clickX, clickY]; ud.current_pointList_for_transform(end+1, :) = [clickX, clickY]; + try set(ud.pointHands_for_transform(:), 'color', [.7 .3 .3]); + catch + end ud.pointHands_for_transform(end+1) = plot(ud.atlasAx, clickX, clickY, 'ro', 'color', [0 .9 0],'LineWidth',2,'markers',4); ud.slice_at_shift_start = ud_slice.slice_num; From 5ca2bc2f96ce1aaf266eeb1d02c3dd9ece092149 Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Wed, 19 Oct 2022 15:31:08 +0100 Subject: [PATCH 07/43] show/hide point counter for slice --- Browsing Functions/sliceBrowser.m | 41 ++++++++++++++++++++++++------- 1 file changed, 32 insertions(+), 9 deletions(-) diff --git a/Browsing Functions/sliceBrowser.m b/Browsing Functions/sliceBrowser.m index 9711a15..5df5fab 100644 --- a/Browsing Functions/sliceBrowser.m +++ b/Browsing Functions/sliceBrowser.m @@ -11,8 +11,8 @@ function sliceBrowser(slice_figure, processed_images_folder, f, reference_size) ud_slice.break = 0; ud_slice.slice_num = 1; ud_slice.key = 1; -ud_slice.pointList = []; -ud_slice.pointHands = []; +ud_slice.pointList = []; % point coordinates +ud_slice.pointHands = []; % graphic objects ud_slice.getPoint = 0; ud_slice.ref_size = reference_size(2:3); @@ -23,6 +23,10 @@ function sliceBrowser(slice_figure, processed_images_folder, f, reference_size) set(ud_slice.sliceAx, 'HitTest', 'off'); ud_slice.im = plotTVslice(zeros(ud_slice.ref_size, 'uint8')); +ud_slice.pointsText = annotation('textbox', [0.88 0.03 0.1 0.05], ... + 'String', '0 point', 'EdgeColor', 'none', 'Color', 'k', 'HorizontalAlignment', 'right'); +ud_slice.pointsText.Visible = 'off'; + % create functions needed to interact with the figure set(ud_slice.im, 'ButtonDownFcn', @(slice_figure,k)sliceClickCallback(slice_figure, k)); set(slice_figure, 'KeyPressFcn', @(slice_figure, keydata)SliceAtlasHotkeyFcn(slice_figure, keydata, f)); @@ -38,6 +42,9 @@ function sliceBrowser(slice_figure, processed_images_folder, f, reference_size) % set up first slice image ud_slice = updateSliceImage(ud_slice); +set(slice_figure, 'UserData', ud_slice); + + % ------------------------------------------------ % Clicking function to register transform points % ------------------------------------------------ @@ -60,6 +67,8 @@ function sliceClickCallback(im, keydata) ud.pointList = []; set(ud.pointHands(:), 'Visible', 'off'); end + + ud.pointsText.String = sprintf('%d point(s)', length(ud.pointHands)); end set(f, 'UserData', ud); @@ -91,14 +100,28 @@ function SliceAtlasHotkeyFcn(fig, keydata, f) % ud.pointList = []; % Try to delete only most recent point - set(ud.pointHands(end), 'Visible', 'off'); - ud.pointHands = ud.pointHands(1:end-1); - ud.pointList = ud.pointList(1:end-1,:); - disp('transform point deleted') + if isempty(ud.pointHands) + disp('There is no transform point to delete.') + else + set(ud.pointHands(end), 'Visible', 'off'); + ud.pointHands = ud.pointHands(1:end-1); + ud.pointList = ud.pointList(1:end-1,:); + disp('transform point deleted') + ud.pointsText.String = sprintf('%d point(s)', length(ud.pointHands)); + + end % t -- transform point mode elseif strcmp(keydata.Key,'t') ud.getPoint = ~ud.getPoint; - if ud.getPoint; disp('transform point mode!'); end + if ud.getPoint + disp('transform point mode!'); + set(ud.pointHands(:), 'Visible', 'on') + ud.pointsText.Visible = 'on'; + ud.pointsText.String = sprintf('%d point(s)', length(ud.pointHands)); + else + ud.pointsText.Visible = 'off'; + set(ud.pointHands(:), 'Visible', 'off') + end else % otherwise -- call function to atlas browser figure(f); @@ -135,8 +158,6 @@ function SliceAtlasHotkeyFcn(fig, keydata, f) end ud.pointHands = gobjects(0); - % set(ud.pointHands(:), 'Visible', 'off'); %TODO WHY hiding? Maybe you - % want to delete/initialize them? ud is from fig.UserData ud.pointList = []; if exist(file_transformations,'file') @@ -151,4 +172,6 @@ function SliceAtlasHotkeyFcn(fig, keydata, f) title_ending = ' (transform points loaded)'; end end + set(ud.pointHands(:), 'Visible', 'off'); %Hide because not in 't' mode + title(['Slice Viewer -- Slice ' num2str(ud.slice_num) '/' num2str(ud.total_num_files) title_ending]) From 3cafbbbf528149aa8b17fbcc2be2bcd8dffa43e4 Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Wed, 19 Oct 2022 15:31:30 +0100 Subject: [PATCH 08/43] show/hide point counter for atlas --- Browsing Functions/AtlasTransformBrowser.m | 24 +++++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/Browsing Functions/AtlasTransformBrowser.m b/Browsing Functions/AtlasTransformBrowser.m index 9b34253..ad39268 100644 --- a/Browsing Functions/AtlasTransformBrowser.m +++ b/Browsing Functions/AtlasTransformBrowser.m @@ -70,6 +70,10 @@ ud.angleText = annotation('textbox', [.7 0.95 0.4 0.05], ... 'EdgeColor', 'none', 'Color', 'k'); +ud.pointsText = annotation('textbox', [0.88 0.03 0.1 0.05], ... + 'String', '0 point', 'EdgeColor', 'none', 'Color', 'k', 'HorizontalAlignment', 'right'); +ud.pointsText.Visible = 'off'; + allData.tv = templateVolume; allData.av = annotationVolume; allData.st = structureTree; @@ -358,6 +362,10 @@ function hotkeyFcn(f, slice_figure, keydata, allData, save_location, save_suffix if ud.getPoint_for_transform disp('transform point mode on'); + % set(ud.pointHands_for_transform(:), 'Visible', 'on') %TODO + ud.pointsText.Visible = 'on'; + ud.pointsText.String = sprintf('%d point(s)', length(ud.pointHands_for_transform)); + ud.currentProbe = 0; % launch transform point mode @@ -382,6 +390,7 @@ function hotkeyFcn(f, slice_figure, keydata, allData, save_location, save_suffix try set(ud.pointHands_for_transform(:), 'Visible', 'off'); + ud.pointsText.Visible = 'off'; catch end @@ -535,11 +544,18 @@ function hotkeyFcn(f, slice_figure, keydata, allData, save_location, save_suffix ud.curr_im = ud.curr_slice_trans; end % if wrong number of points clicked - catch + catch mexc1 + if mexc1.identifier == "images:geotrans:requiredNonCollinearPoints" + disp(['Unable to transform -- ', mexc1.message]) + + else + ref_mode = true; disp(['Unable to transform -- ' num2str(size(ud_slice.pointList,1)) ... ' slice points and ' num2str(size(ud.current_pointList_for_transform,1)) ' reference points']); key_letter = 'h'; + + end end end % if not doing transform, just show reference atlas @@ -687,6 +703,8 @@ function hotkeyFcn(f, slice_figure, keydata, allData, save_location, save_suffix end disp('transform point deleted'); + ud.pointsText.String = sprintf('%d point(s)', length(ud.pointHands_for_transform)); + elseif ud.currentProbe ud.pointList{ud.currentProbe,1} = ud.pointList{ud.currentProbe,1}(1:end-1,:); @@ -1133,6 +1151,10 @@ function atlasClickCallback(im, keydata, slice_figure, save_location) ud.curr_slice_num = ud.slice_at_shift_start+ud.slice_shift; ud.loaded = 0; ud.clicked = true; + + ud.pointsText.String = sprintf('%d point(s)', length(ud.pointHands_for_transform)); + + end set(f, 'UserData', ud); From 13d6c1bb6759fccbb577323b1e4875b7e85aacca Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Wed, 19 Oct 2022 16:14:27 +0100 Subject: [PATCH 09/43] point counter behaviour --- Browsing Functions/sliceBrowser.m | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/Browsing Functions/sliceBrowser.m b/Browsing Functions/sliceBrowser.m index 5df5fab..654dc81 100644 --- a/Browsing Functions/sliceBrowser.m +++ b/Browsing Functions/sliceBrowser.m @@ -172,6 +172,11 @@ function SliceAtlasHotkeyFcn(fig, keydata, f) title_ending = ' (transform points loaded)'; end end - set(ud.pointHands(:), 'Visible', 'off'); %Hide because not in 't' mode - + if ud.getPoint + set(ud.pointHands(:), 'Visible', 'on'); + ud.pointsText.String = sprintf('%d point(s)', length(ud.pointHands)); + else + set(ud.pointHands(:), 'Visible', 'off'); %Hide because not in 't' mode + ud.pointsText.String = sprintf('%d point(s)', length(ud.pointHands)); + end title(['Slice Viewer -- Slice ' num2str(ud.slice_num) '/' num2str(ud.total_num_files) title_ending]) From 8b43f40b7380abec41edd529b4ecbeb790ac4169 Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Wed, 19 Oct 2022 16:34:54 +0100 Subject: [PATCH 10/43] point counter behaviour --- Browsing Functions/AtlasTransformBrowser.m | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/Browsing Functions/AtlasTransformBrowser.m b/Browsing Functions/AtlasTransformBrowser.m index ad39268..ff5e18b 100644 --- a/Browsing Functions/AtlasTransformBrowser.m +++ b/Browsing Functions/AtlasTransformBrowser.m @@ -852,6 +852,7 @@ function updateSlice(f, evt, allData, slice_figure, save_location) ud_slice.pointList = transform_data.transform_points{2}; else ud_slice.pointList = []; + ud.pointHands_for_transform = gobjects(0); end set(slice_figure, 'UserData', ud_slice); @@ -868,6 +869,9 @@ function updateSlice(f, evt, allData, slice_figure, save_location) ud.curr_slice_num = ud.slice_at_shift_start+ud.slice_shift; ud.histology_overlay = 1; + + ud.pointsText.String = sprintf('%d point(s)', length(ud.pointHands_for_transform)); + set(ud.text,'Visible','off'); fill([5 5 250 250],[5 50 50 5],[0 0 0]); ud.text(end+1) = text(5,15,['Slice ' num2str(ud.slice_at_shift_start+ud.slice_shift)],'color','white'); @@ -880,7 +884,9 @@ function updateSlice(f, evt, allData, slice_figure, save_location) set(ud.im, 'CData', ud.ref); ud.curr_im = ud.ref; set(f, 'UserData', ud); set(ud.text,'Visible','off'); - fill([5 5 250 250],[5 50 50 5],[0 0 0]); ud.text(end+1) = text(5,15,['Slice ' num2str(ud.slice_at_shift_start+ud.slice_shift) ' - no transform'],'color','white'); + fill([5 5 250 250],[5 50 50 5],[0 0 0]); ud.text(end+1) = text(5,15,['Slice ' num2str(ud.slice_at_shift_start+ud.slice_shift) ' - no transform'],'color','white'); + ud.pointsText.String = sprintf('%d point(s)', length(ud.pointHands_for_transform)); + end end From 98cc667394a4b8c8e164fbab78a62eef3a580a34 Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Wed, 19 Oct 2022 17:44:30 +0100 Subject: [PATCH 11/43] points show status --- Browsing Functions/AtlasTransformBrowser.m | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/Browsing Functions/AtlasTransformBrowser.m b/Browsing Functions/AtlasTransformBrowser.m index ff5e18b..297760c 100644 --- a/Browsing Functions/AtlasTransformBrowser.m +++ b/Browsing Functions/AtlasTransformBrowser.m @@ -855,7 +855,13 @@ function updateSlice(f, evt, allData, slice_figure, save_location) ud.pointHands_for_transform = gobjects(0); end set(slice_figure, 'UserData', ud_slice); - + + if ud.getPoint_for_transform + set(ud.pointHands_for_transform, 'Visible', 'on') + else + set(ud.pointHands_for_transform, 'Visible', 'off') + end + % load allen ref location ud.currentSlice = transform_data.allen_location{1}; ud.currentAngle = transform_data.allen_location{2}; From 9e614f55f7d53eee9ebef79088bb22d4b3fa744f Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Mon, 31 Oct 2022 17:05:54 +0000 Subject: [PATCH 12/43] Better action. But perhaps there is a bug about coronal sections? --- Browsing Functions/AtlasTransformBrowser.m | 5 +++++ Browsing Functions/sliceBrowser.m | 6 +++++- 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/Browsing Functions/AtlasTransformBrowser.m b/Browsing Functions/AtlasTransformBrowser.m index 297760c..f307164 100644 --- a/Browsing Functions/AtlasTransformBrowser.m +++ b/Browsing Functions/AtlasTransformBrowser.m @@ -190,6 +190,11 @@ function hotkeyFcn(f, slice_figure, keydata, allData, save_location, save_suffix disp(['probe point mode -- selecting probe ' num2str(ud.currentProbe) ' (' ud.ProbeColor{ud.currentProbe} ')']); ud.getPoint_for_transform = false; + try + set(ud.pointHands_for_transform(:), 'Visible', 'off'); + ud.pointsText.Visible = 'off'; + catch + end % show Transformed Slice & Probage Viewer, if not already showing if ~ud.slice_at_shift_start; add = 1; else; add = 0; end diff --git a/Browsing Functions/sliceBrowser.m b/Browsing Functions/sliceBrowser.m index 654dc81..df8ac47 100644 --- a/Browsing Functions/sliceBrowser.m +++ b/Browsing Functions/sliceBrowser.m @@ -61,7 +61,8 @@ function sliceClickCallback(im, keydata) clickY = round(keydata.IntersectionPoint(2)); ud.pointList(end+1, :) = [clickX, ud.ref_size(1) - clickY]; - ud.pointHands(end+1) = plot(ud.sliceAx, clickX, clickY, 'ro', 'color', [0 .5 0],'linewidth',2,'markers',4); + set(ud.pointHands,'Color',[.7 .3 .3]) + ud.pointHands(end+1) = plot(ud.sliceAx, clickX, clickY, 'o', 'color', [0 .9 0],'linewidth',2,'markers',4); if clickX < 100 && (ud.ref_size(1) - clickY) < 100 % if click in corner, break ud.pointList = []; @@ -134,6 +135,7 @@ function SliceAtlasHotkeyFcn(fig, keydata, f) function ud = updateSliceImage(ud) + %TODO occasionally pointList has more items than pointHands title_ending = ''; @@ -169,6 +171,8 @@ function SliceAtlasHotkeyFcn(fig, keydata, f) for i = 1:size(ud.pointList,1) ud.pointHands(end+1) = plot(ud.sliceAx, ud.pointList(i,1), ud.ref_size(1) - ud.pointList(i,2), 'ro', 'color', [0 .5 0],'linewidth',2,'markers',4); end + set(ud.pointHands(end),'color', [0 .9 0]); + title_ending = ' (transform points loaded)'; end end From f2eafc61dc204475d340a453dcf9a4fd9b401bed Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Wed, 1 Feb 2023 18:00:38 +0000 Subject: [PATCH 13/43] deletion of points --- Browsing Functions/AtlasTransformBrowser.m | 10 +++++++--- Browsing Functions/sliceBrowser.m | 4 ++-- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/Browsing Functions/AtlasTransformBrowser.m b/Browsing Functions/AtlasTransformBrowser.m index f307164..1af1bb0 100644 --- a/Browsing Functions/AtlasTransformBrowser.m +++ b/Browsing Functions/AtlasTransformBrowser.m @@ -483,7 +483,7 @@ function hotkeyFcn(f, slice_figure, keydata, allData, save_location, save_suffix xlabel(ud.atlasAx, 'Scroll slices along D/V axis') end - case 'leftarrow' % scroll along A/P axis + case 'leftarrow' % scroll along along slice images ud.scrollMode = 3; if ~ud.slice_at_shift_start ud.slice_at_shift_start = ud_slice.slice_num; @@ -694,8 +694,12 @@ function hotkeyFcn(f, slice_figure, keydata, allData, save_location, save_suffix % disp('current transform erased'); % Try to delete only the most recent point + if isempty(ud.pointHands_for_transform) + disp('There is no transform point to delete.') + else + ud.current_pointList_for_transform = ud.current_pointList_for_transform(1:end-1,:); - ud_slice.pointList = ud_slice.pointList(1:end-1,:); %TODO this should be accompanied by deleting circles on Slive Viewer + ud_slice.pointList = ud_slice.pointList(1:end-1,:); %TODO this should be accompanied by deleting circles on Slice Viewer set(slice_figure, 'UserData', ud_slice); if ~isempty(ud.pointHands_for_transform) % remove circle for most revent point @@ -709,7 +713,7 @@ function hotkeyFcn(f, slice_figure, keydata, allData, save_location, save_suffix disp('transform point deleted'); ud.pointsText.String = sprintf('%d point(s)', length(ud.pointHands_for_transform)); - + end elseif ud.currentProbe ud.pointList{ud.currentProbe,1} = ud.pointList{ud.currentProbe,1}(1:end-1,:); diff --git a/Browsing Functions/sliceBrowser.m b/Browsing Functions/sliceBrowser.m index df8ac47..0e5199a 100644 --- a/Browsing Functions/sliceBrowser.m +++ b/Browsing Functions/sliceBrowser.m @@ -104,7 +104,7 @@ function SliceAtlasHotkeyFcn(fig, keydata, f) if isempty(ud.pointHands) disp('There is no transform point to delete.') else - set(ud.pointHands(end), 'Visible', 'off'); + set(ud.pointHands(end), 'Visible', 'off'); %TODO isn't this better to be deleted? But if you do, you have a problem for callbacks ud.pointHands = ud.pointHands(1:end-1); ud.pointList = ud.pointList(1:end-1,:); disp('transform point deleted') @@ -156,7 +156,7 @@ function SliceAtlasHotkeyFcn(fig, keydata, f) [processed_image_name(1:end-4) '_transform_data.mat']); try - delete(findobj(ud.im.Parent, 'Type','line','color',[0 0.5 0],'Marker','o')); % delete all the circles + delete(findobj(ud.im.Parent, 'Type','line','Marker','o')); % delete all the circles end ud.pointHands = gobjects(0); From ef13d6b7b31b9943b0bae2207bb53c6bb9bf9f62 Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Wed, 1 Feb 2023 18:01:01 +0000 Subject: [PATCH 14/43] causing a bug of mismatching numbers of points --- Browsing Functions/AtlasTransformBrowser.m | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Browsing Functions/AtlasTransformBrowser.m b/Browsing Functions/AtlasTransformBrowser.m index 1af1bb0..2618458 100644 --- a/Browsing Functions/AtlasTransformBrowser.m +++ b/Browsing Functions/AtlasTransformBrowser.m @@ -863,7 +863,8 @@ function updateSlice(f, evt, allData, slice_figure, save_location) ud_slice.pointList = []; ud.pointHands_for_transform = gobjects(0); end - set(slice_figure, 'UserData', ud_slice); + % set(slice_figure, 'UserData', ud_slice); + %TODO causing a bug of mismatching numbers of points if ud.getPoint_for_transform set(ud.pointHands_for_transform, 'Visible', 'on') From 5207316ac629821226df80bedc552a2c2c605b73 Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Wed, 1 Feb 2023 18:01:45 +0000 Subject: [PATCH 15/43] Not sure why points need to be hidden when moving along an axis or tilting --- Browsing Functions/AtlasTransformBrowser.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Browsing Functions/AtlasTransformBrowser.m b/Browsing Functions/AtlasTransformBrowser.m index 2618458..197cd6d 100644 --- a/Browsing Functions/AtlasTransformBrowser.m +++ b/Browsing Functions/AtlasTransformBrowser.m @@ -1003,7 +1003,7 @@ function updateSlice(f, evt, allData, slice_figure, save_location) end ud.ref = uint8(angle_slice); - set(ud.pointHands_for_transform(:), 'Visible', 'off'); + % set(ud.pointHands_for_transform(:), 'Visible', 'off'); %TODO is this needed??? end From 8ad7313c738e969163483d45ea4c9dca5da9b367 Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Sat, 3 Jun 2023 09:44:37 +0100 Subject: [PATCH 16/43] color change and beep for warning --- Browsing Functions/AtlasTransformBrowser.m | 29 ++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/Browsing Functions/AtlasTransformBrowser.m b/Browsing Functions/AtlasTransformBrowser.m index 197cd6d..7e2857a 100644 --- a/Browsing Functions/AtlasTransformBrowser.m +++ b/Browsing Functions/AtlasTransformBrowser.m @@ -713,6 +713,21 @@ function hotkeyFcn(f, slice_figure, keydata, allData, save_location, save_suffix disp('transform point deleted'); ud.pointsText.String = sprintf('%d point(s)', length(ud.pointHands_for_transform)); + + if abs(length(ud.pointHands_for_transform) - length(ud_slice.pointHands)) >= 2 + beep + if length(ud.pointHands_for_transform) - length(ud_slice.pointHands) >= 2 + ud.pointsText.Color = 'red'; + ud_slice.pointsText.Color = 'blue'; + else + ud.pointsText.Color = 'blue'; + ud_slice.pointsText.Color = 'red'; + end + else + ud.pointsText.Color = 'black'; + ud_slice.pointsText.Color = 'black'; + end + end elseif ud.currentProbe @@ -1176,6 +1191,20 @@ function atlasClickCallback(im, keydata, slice_figure, save_location) ud.pointsText.String = sprintf('%d point(s)', length(ud.pointHands_for_transform)); + if abs(length(ud.pointHands_for_transform) - length(ud_slice.pointHands)) >= 2 + beep + if length(ud.pointHands_for_transform) - length(ud_slice.pointHands) >= 2 + ud.pointsText.Color = 'red'; + ud_slice.pointsText.Color = 'blue'; + + else + ud.pointsText.Color = 'blue'; + ud_slice.pointsText.Color = 'red'; + end + else + ud.pointsText.Color = 'black'; + ud_slice.pointsText.Color = 'black'; + end end set(f, 'UserData', ud); From e913028507b4082bbbb69be5b1bbdfcb5457fa1f Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Thu, 3 Aug 2023 16:52:39 +0100 Subject: [PATCH 17/43] sliceBrowser.m Browsing Functions AtlasTransformBrowser.m Browsing Functions --- Browsing Functions/AtlasTransformBrowser.m | 2 +- Browsing Functions/sliceBrowser.m | 35 +++++++++++++++++----- 2 files changed, 28 insertions(+), 9 deletions(-) diff --git a/Browsing Functions/AtlasTransformBrowser.m b/Browsing Functions/AtlasTransformBrowser.m index 7e2857a..af86a1d 100644 --- a/Browsing Functions/AtlasTransformBrowser.m +++ b/Browsing Functions/AtlasTransformBrowser.m @@ -715,7 +715,7 @@ function hotkeyFcn(f, slice_figure, keydata, allData, save_location, save_suffix ud.pointsText.String = sprintf('%d point(s)', length(ud.pointHands_for_transform)); if abs(length(ud.pointHands_for_transform) - length(ud_slice.pointHands)) >= 2 - beep + % beep if length(ud.pointHands_for_transform) - length(ud_slice.pointHands) >= 2 ud.pointsText.Color = 'red'; ud_slice.pointsText.Color = 'blue'; diff --git a/Browsing Functions/sliceBrowser.m b/Browsing Functions/sliceBrowser.m index 0e5199a..d8bc290 100644 --- a/Browsing Functions/sliceBrowser.m +++ b/Browsing Functions/sliceBrowser.m @@ -60,16 +60,25 @@ function sliceClickCallback(im, keydata) clickX = round(keydata.IntersectionPoint(1)); clickY = round(keydata.IntersectionPoint(2)); - ud.pointList(end+1, :) = [clickX, ud.ref_size(1) - clickY]; + % ud.pointList(end+1, :) = [clickX, ud.ref_size(1) - clickY]; % often pointList is shorter than pointHands set(ud.pointHands,'Color',[.7 .3 .3]) ud.pointHands(end+1) = plot(ud.sliceAx, clickX, clickY, 'o', 'color', [0 .9 0],'linewidth',2,'markers',4); - - if clickX < 100 && (ud.ref_size(1) - clickY) < 100 % if click in corner, break - ud.pointList = []; - set(ud.pointHands(:), 'Visible', 'off'); - end - ud.pointsText.String = sprintf('%d point(s)', length(ud.pointHands)); + % make sure what you see are what you have + phx = zeros(length(ud.pointHands), 1); + phy = zeros(length(ud.pointHands), 1); + for i = 1:length(ud.pointHands) + phx(i) = ud.pointHands(i).XData; + phy(i) = ud.pointHands(i).YData; + end + ud.pointList = [phx, ud.ref_size(1) - phy]; %TODO + + if clickX < 100 && (ud.ref_size(1) - clickY) < 100 % if click in corner, break + ud.pointList = []; + set(ud.pointHands(:), 'Visible', 'off'); + end + + ud.pointsText.String = sprintf('%d point(s)', length(ud.pointHands)); end set(f, 'UserData', ud); @@ -106,7 +115,17 @@ function SliceAtlasHotkeyFcn(fig, keydata, f) else set(ud.pointHands(end), 'Visible', 'off'); %TODO isn't this better to be deleted? But if you do, you have a problem for callbacks ud.pointHands = ud.pointHands(1:end-1); - ud.pointList = ud.pointList(1:end-1,:); + + % make sure the length of ud.pointHands and the rows of ud.pointList match + % make sure what you see are what you have + phx = zeros(length(ud.pointHands), 1); + phy = zeros(length(ud.pointHands), 1); + for i = 1:length(ud.pointHands) + phx(i) = ud.pointHands(i).XData; + phy(i) = ud.pointHands(i).YData; + end + ud.pointList = [phx, ud.ref_size(1) - phy]; + disp('transform point deleted') ud.pointsText.String = sprintf('%d point(s)', length(ud.pointHands)); From 4de3d02dd2eafc33333b3772ad87b56177cec1fd Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Thu, 3 Aug 2023 16:52:55 +0100 Subject: [PATCH 18/43] depth correction --- Browsing Functions/accf2pxs_mm.m | 20 ++++++++++ Browsing Functions/allenCCF2pxs_mm.m | 18 +++++++++ Browsing Functions/apdvml2info.m | 56 ++++++++++++++++++++++++++++ 3 files changed, 94 insertions(+) create mode 100644 Browsing Functions/accf2pxs_mm.m create mode 100644 Browsing Functions/allenCCF2pxs_mm.m create mode 100644 Browsing Functions/apdvml2info.m diff --git a/Browsing Functions/accf2pxs_mm.m b/Browsing Functions/accf2pxs_mm.m new file mode 100644 index 0000000..ded57e6 --- /dev/null +++ b/Browsing Functions/accf2pxs_mm.m @@ -0,0 +1,20 @@ +function y_mm = accf2pxs_mm(x_mm) +% convert Allen Common Coordinates into standard values (Paxinos and Franklin) +% x_mm Allen Common Coordinates +% if this is scaler, x_mm should be depth +% if this is n x 3 array, x_mm is an array of the same (or similar) size +% y_mm Paxinos and Franklin coordinates + +if size(x_mm,2) ==1 + y_mm = 0.921 * x_mm - 0.834; % global + % y = 0.908 * x - 0.957; % local + +elseif size(x_mm,2) == 3 + % assume that the 2nd column is the DV + y_mm = x_mm; + y_mm(:,2) = 0.921 * x_mm(:,2) - 0.834; % global + % y(:,2) = 0.908 * x(:,2) - 0.957; % local + +else + error('unexpected size') +end \ No newline at end of file diff --git a/Browsing Functions/allenCCF2pxs_mm.m b/Browsing Functions/allenCCF2pxs_mm.m new file mode 100644 index 0000000..999221e --- /dev/null +++ b/Browsing Functions/allenCCF2pxs_mm.m @@ -0,0 +1,18 @@ +function y_mm = accf2pxs_mm(x_mm) +% convert Allen Common Coordinates into standard values (Paxinos and Franklin) +% x_mm Allen Common Coordinates +% y_mm Paxinos and Franklin coordinates + +if size(x_mm,2) ==1 + y_mm = 0.921 * x_mm - 0.834; % global + % y = 0.908 * x - 0.957; % local + +elseif size(x_mm,2) == 3 + % assume that the 2nd column is the DV + y_mm = x_mm; + y_mm(:,2) = 0.921 * x_mm(:,2) - 0.834; % global + % y(:,2) = 0.908 * x(:,2) - 0.957; % local + +else + error('unexpected size') +end \ No newline at end of file diff --git a/Browsing Functions/apdvml2info.m b/Browsing Functions/apdvml2info.m new file mode 100644 index 0000000..8735ce2 --- /dev/null +++ b/Browsing Functions/apdvml2info.m @@ -0,0 +1,56 @@ +function Tapdvml = apdvml2info(apdvml_points, av, st) +% apdvml_points has three columns for ap, dv, and ml +% bregma is 1 x 3 in size +% atlas_resolution in mm +% av 1320 ,800, 1140 +% st +% +% see also +% Analyze_Clicked_Points.m + +% generate needed values +bregma = allenCCFbregma(); % bregma position in reference data space +atlas_resolution = 0.010; % mm + + +ap = -(apdvml_points(:,1)-bregma(1))*atlas_resolution; +dv = (apdvml_points(:,2)-bregma(2))*atlas_resolution; +ml = (apdvml_points(:,3)-bregma(3))*atlas_resolution; + +% roi_location_curr = [ap dv ml]; + +% initialize array of region annotations +roi_annotation_curr = cell(size(apdvml_points,1),3); + +% loop through every point to get ROI locations and region annotations +tf_rows = true(size(apdvml_points,1),1); +for point = 1:size(apdvml_points,1) + + % find the annotation, name, and acronym of the current ROI pixel + depth = ceil(apdvml_points(point,2)); + if depth >= 1 + + ann = av(ceil(apdvml_points(point,1)), ... + depth, ... %NOTE this can take the value of 0 and cause an error (must be >= 1) + ceil(apdvml_points(point,3))); + name = st.safe_name{ann}; + acr = st.acronym{ann}; + + roi_annotation_curr{point,1} = ann; + roi_annotation_curr{point,2} = name; + roi_annotation_curr{point,3} = acr; + + else + % ignore this point, too dorsal + roi_annotation_curr{point,1} = NaN; + roi_annotation_curr{point,2} = ''; + roi_annotation_curr{point,3} = ''; + + end + +end + +Tapdvml = [table(ap, dv, accf2pxs_mm(dv), ml, 'VariableNames',{'ap_mm','dv_mm', 'dv_mm_paxinos', 'ml_mm'}), ... + cell2table(roi_annotation_curr, 'VariableNames', {'annotation', 'name', 'acronym'})]; + +end \ No newline at end of file From f522cf993e2ddf64eb88c26bd5f994e41b174b47 Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Fri, 4 Aug 2023 17:40:17 +0100 Subject: [PATCH 19/43] bug fix https://github.com/cortex-lab/allenCCF/issues/70 --- Browsing Functions/AtlasTransformBrowser.m | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Browsing Functions/AtlasTransformBrowser.m b/Browsing Functions/AtlasTransformBrowser.m index af86a1d..78cd002 100644 --- a/Browsing Functions/AtlasTransformBrowser.m +++ b/Browsing Functions/AtlasTransformBrowser.m @@ -779,8 +779,9 @@ function hotkeyFcn(f, slice_figure, keydata, allData, save_location, save_suffix disp('atlas location saved') % save transformed histology image - current_slice_image = imread(fullfile(save_location, [slice_name '.tif'])); -% current_slice_image = flip(get(ud_slice.im, 'CData'),1); + % current_slice_image = imread(fullfile(save_location, [slice_name '.tif']));%TODO no imreisze() used, size mismatch creates a huge error in transform + current_slice_image = flip(get(ud_slice.im, 'CData'),1); + R = imref2d(size(ud.ref)); curr_slice_trans = imwarp(current_slice_image, ud.transform, 'OutputView',R); imwrite(curr_slice_trans, fullfile(folder_transformations, [slice_name '_transformed.tif'])) From 6f73199a9e57716beb37f824af2816579d87ee14 Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Tue, 15 Aug 2023 14:28:05 +0100 Subject: [PATCH 20/43] C:\Users\phar0528\OneDrive - Nexus365\Private_Dropbox\Projects\allenCCF\Browsing Functions\plot_and_compute_probe_positions_from2.m --- .../plot_and_compute_probe_positions_from2.m | 376 ++++++++++++++++++ 1 file changed, 376 insertions(+) create mode 100644 Browsing Functions/plot_and_compute_probe_positions_from2.m diff --git a/Browsing Functions/plot_and_compute_probe_positions_from2.m b/Browsing Functions/plot_and_compute_probe_positions_from2.m new file mode 100644 index 0000000..80475a7 --- /dev/null +++ b/Browsing Functions/plot_and_compute_probe_positions_from2.m @@ -0,0 +1,376 @@ +function [probes, fwireframe, fwireframe2, fig_probes, T_probes, Tapdvml_contacts, T_borders] ... + = plot_and_compute_probe_positions_from2(av, st, subject_id, plane,... + processed_images_folder, probe_save_name_suffix,probes_to_analyze, probe_lengths,... + active_probe_length, probe_radius) +% Create 3D plot for probe positions, 2D plots for brain structures along +% each probe, and return tables for structure names and coordinates +% +% +% SYNTAX +% [] = plot_and_compute_probe_positions(fwireframe, black_brain, probes, subject_id, plane,... +% use_tip_to_get_reference_probe_length, probe_insertion_direction) +% [D] = func1(A,B) +% [D] = func1(____,'Param',value) +% +% longer description may come here +% +% INPUT ARGUMENTS +% av uint16 +% 3 dimensional array of integers representing the different +% brain regions in the model space +% +% st table +% +% subject_id char +% Animal name +% +% plane 'coronal' | 'sagittal' | 'transverse' +% +% processed_images_folder +% char +% +% probe_save_name_suffix +% char +% Needed to access stored data +% +% probes_to_analyze +% 'all' | a row vector of positive integers specifying probes +% + +% B 0 (default) | non-negative integers +% (Optional) Description about B comes here. +% +% +% OPTIONAL PARAMETER/VALUE PAIRS +% 'C' 'on' (default) | 'off' +% (Optional) Description about 'C' comes here. +% +% +% OUTPUT ARGUMENTS +% D cell array +% Description about D comes here. +% +% Written by Kouichi C. Nakamura Ph.D. +% MRC Brain Network Dynamics Unit +% University of Oxford +% kouichi.c.nakamura@gmail.com +% 09-Aug-2023 11:12:42 +% +% See also +% doc + + + +% load the reference brain annotations +% if ~exist('av','var') || ~exist('st','var') +% disp('loading reference atlas...') +% av = readNPY(annotation_volume_location); +% st = loadStructureTree(structure_tree_location); +% end + +arguments + + av (:,:,:) uint16 + st (:, 22) table + subject_id (1, :) char + plane (1, :) char {mustBeMember(plane,{'coronal','sagittal','transverse'})} + processed_images_folder (1, :) char {mustBeFolder} + probe_save_name_suffix (1, :) char + probes_to_analyze % 'all' | A vector of integer indices for probes + probe_lengths {mustBeVector} + active_probe_length (1,1) + probe_radius (1,1) + +end + + +scaling_factor = false; +show_region_table = true; +black_brain = true; +probe_insertion_direction = 'down'; +show_parent_category = false; +distance_past_tip_to_plot = .5; + + + +% select the plane for the viewer +if strcmp(plane,'coronal') + av_plot = av; +elseif strcmp(plane,'sagittal') + av_plot = permute(av,[3 2 1]); +elseif strcmp(plane,'transverse') + av_plot = permute(av,[2 3 1]); +end + +% load probe points +probePoints = load(fullfile(processed_images_folder, ['probe_points' probe_save_name_suffix])); +ProbeColors = .75*[1.3 1.3 1.3; 1 .75 0; .3 1 1; .4 .6 .2; 1 .35 .65; .7 .7 .9; .65 .4 .25; .7 .95 .3; .7 0 0; .6 0 .7; 1 .6 0]; +% order of colors: {'white','gold','turquoise','fern','bubble gum','overcast sky','rawhide', 'green apple','red','purple','orange'}; +fwireframe = []; + +% % scale active_probe_length appropriately +% active_probe_length = active_probe_length*100; + +% determine which probes to analyze +if strcmp(probes_to_analyze,'all') + probes = 1:size(probePoints.pointList.pointList,1); +else + probes = probes_to_analyze; +end + + + +%% + +% create a new figure with wireframe +fwireframe = plotBrainGrid([], [], fwireframe, black_brain); +hold on; +fwireframe.InvertHardcopy = 'off'; + +fwireframe2 = plotBrainGrid([], [], fwireframe, black_brain); %TODO what is this one for????????? +hold on; +fwireframe2.InvertHardcopy = 'off'; + +fig_probes = gobjects(1,length(probes)); +borders_table = cell(1,length(probes)) %TODO +t_size = [size(probePoints.pointList.pointList,1),1]; +T_probes = table(repmat(string(subject_id), t_size),... + strings(t_size), ... + NaN(t_size),... + strings(t_size), ... + NaN(t_size),... + NaN(t_size),... + NaN(t_size),... + NaN(t_size),... + NaN(t_size),... + NaN(t_size),... + strings(t_size), ... + 'VariableNames',{'subject_id','session_id','probe_id','probe_AB',... + 'depth_reported_um','depth_from_ACCF_um','scaling', ... + 'scaling2',... + 'tip_active_um', 'top_active_um', 'probe_note'}); +T_probes.probe_id(probes') = probes'; +T_probes.probe_lengths(probes) = probe_lengths(probes)' * 1000; + + +c_Tapdvml_contacts = cell(length(probes),1); + +for selected_probe = probes + + % get the probe points for the currently analyzed probe + if strcmp(plane,'coronal') + curr_probePoints = probePoints.pointList.pointList{selected_probe,1}(:, [3 2 1]); + elseif strcmp(plane,'sagittal') + curr_probePoints = probePoints.pointList.pointList{selected_probe,1}(:, [1 2 3]); + elseif strcmp(plane,'transverse') + curr_probePoints = probePoints.pointList.pointList{selected_probe,1}(:, [1 3 2]); + end + + assert(size(curr_probePoints,1) ==2, 'Only accepts 2 points (the tip and surface) per probe') + + % get user-defined probe length from experiment + if length(probe_lengths) > 1 + probe_length = probe_lengths(selected_probe); + else + probe_length = probe_lengths; + end + + % get the scaling-factor method to use + if scaling_factor + use_tip_to_get_reference_probe_length = false; + reference_probe_length = probe_length * scaling_factor; + disp(['probe scaling of ' num2str(scaling_factor) ' determined by user input']); + else + use_tip_to_get_reference_probe_length = true; + disp(['getting probe scaling from histology data...']); + end + + % get line of best fit through points + % m is the mean value of each dimension; p is the eigenvector for largest eigenvalue + [m,p,s] = best_fit_line(curr_probePoints(:,1), curr_probePoints(:,2), curr_probePoints(:,3)); %TODO + if isnan(m(1)) + disp(['no points found for probe ' num2str(selected_probe)]) + continue + end + + % ensure proper orientation: want 0 at the top of the brain and positive distance goes down into the brain + if p(2)<0 + p = -p; + end + + + % determine "origin" at top of brain -- step upwards along tract direction until tip of brain / past cortex + ann = 10; + out_of_brain = false; + while ~(ann==1 && out_of_brain) % && distance_stepped > .5*active_probe_length) + m = m-p; % step 10um, backwards up the track + ann = av(round(m(1)),round(m(2)),round(m(3))); %until hitting the top + if strcmp(st.safe_name(ann), 'root') + % make sure this isn't just a 'root' area within the brain + m_further_up = m - p*20; % is there more brain 200 microns up along the track? + ann_further_up = av(round(max(1,m_further_up(1))),round(max(1,m_further_up(2))),round(max(1,m_further_up(3)))); + if strcmp(st.safe_name(ann_further_up), 'root') + out_of_brain = true; + end + end + end + %NOTE m is now the highest point + + %% focus on wireframe plot + figure(fwireframe); + + % plot probe points + hp = plot3(curr_probePoints(:,1), curr_probePoints(:,3),... + curr_probePoints(:,2), '.','linewidth',2, 'color',[ProbeColors(selected_probe,:) .2],'markers',10); + + % plot brain entry point + plot3(m(1), m(3), m(2), 'r*','linewidth',1) + + % use the deepest clicked point as the tip of the probe, if no scaling provided (scaling_factor = false) + if use_tip_to_get_reference_probe_length + % find length of probe in reference atlas space + if strcmp(probe_insertion_direction, 'down') + [depth, tip_index] = max(curr_probePoints(:,2)); + elseif strcmp(probe_insertion_direction, 'up') + [depth, tip_index] = min(curr_probePoints(:,2)); + end + reference_probe_length_tip = sqrt(sum((curr_probePoints(tip_index,:) - m).^2)); + + % and the corresponding scaling factor + shrinkage_factor = (reference_probe_length_tip / 100) / probe_length; + + % display the scaling + disp(['probe length of ' num2str(reference_probe_length_tip/100) ... + ' mm in reference atlas space compared to a reported ' num2str(probe_length) ' mm']); + disp(['probe scaling of ' num2str(shrinkage_factor)]); disp(' '); + + T_probes.depth_from_ACCF_um(selected_probe) = reference_probe_length_tip/100 * 1000; + T_probes.scaling(selected_probe) = shrinkage_factor; + + % plot line the length of the probe in reference space + probe_length_histo = round(reference_probe_length_tip); + % probe_length_histo is the distance from the tip to the surface + % measured in Allen CCF space. + + % if scaling_factor is user-defined as some number, use it to plot the length of the probe + else + probe_length_histo = round(reference_probe_length * 100); + end + + % find the percent of the probe occupied by electrodes + % percent_of_tract_with_active_sites = min([active_probe_length*100 / (probe_length*100), 1.0]); %TODO here probe_length is used, but this should be ignored + % active_site_start = probe_length_histo*(1-percent_of_tract_with_active_sites); + + %NOTE ignore the probe_length. Use probe_length_histo instead. + percent_of_tract_with_active_sites = min([active_probe_length*100 / probe_length_histo, 1.0]); + active_site_start = probe_length_histo*(1-percent_of_tract_with_active_sites); + active_probe_position = round([active_site_start probe_length_histo]); %TODO use this? + + % plot line the length of the active probe sites in reference space + plot3(m(1)+p(1)*[active_probe_position(1) active_probe_position(2)], ... + m(3)+p(3)*[active_probe_position(1) active_probe_position(2)], ... + m(2)+p(2)*[active_probe_position(1) active_probe_position(2)], ... + 'Color', ProbeColors(selected_probe,:), 'LineWidth', 1); + % plot line the length of the entire probe in reference space + plot3(m(1)+p(1)*[1 probe_length_histo], m(3)+p(3)*[1 probe_length_histo], m(2)+p(2)*[1 probe_length_histo], ... + 'Color', ProbeColors(selected_probe,:), 'LineWidth', 1, 'LineStyle',':'); + + + %% ---------------------------------------------------------------- + % Get and plot brain region labels along the extent of each probe + % ---------------------------------------------------------------- + + % convert error radius into mm + error_length = round(probe_radius / 10); + + % find and regions the probe goes through, confidence in those regions, and plot them + borders_table{selected_probe} = plotDistToNearestToTip(m, p, av_plot, st, probe_length_histo, ... + error_length, active_site_start, distance_past_tip_to_plot, ... + show_parent_category, show_region_table, plane); % plots confidence score based on distance to nearest region along probe + title(['Probe ' num2str(selected_probe)],'color',ProbeColors(selected_probe,:)) + + T_probes.tip_active_um(selected_probe) = probe_length_histo *10; %TODO + T_probes.top_active_um(selected_probe) = active_site_start * 10; %TODO + + fig_probes(selected_probe) = gcf; + + + %%%% coordinates, cf. Analyze_Clicked_Points.m + % + % m is the brain entry point + + Tapdvml_m = apdvml2info(m, av, st, plane); + + Tapdvml_tip = apdvml2info(curr_probePoints(tip_index,:), av, st, plane); + + % measure the distance + + tip2surface_mm = sqrt((Tapdvml_tip.ap_mm - Tapdvml_m.ap_mm)^2 + ... + (Tapdvml_tip.dv_mm - Tapdvml_m.dv_mm)^2 + ... + (Tapdvml_tip.ml_mm - Tapdvml_m.ml_mm)^2); + + tip2surface_mm_paxinos = sqrt((Tapdvml_tip.ap_mm - Tapdvml_m.ap_mm)^2 + ... + (Tapdvml_tip.dv_mm_paxinos - Tapdvml_m.dv_mm_paxinos)^2 + ... + (Tapdvml_tip.ml_mm - Tapdvml_m.ml_mm)^2); + + top_active = (m * active_probe_length + ... + curr_probePoints(tip_index,:) * (tip2surface_mm - active_probe_length))... + /tip2surface_mm; + + % obtained the information for all the 384 channels + + + a = [linspace(curr_probePoints(tip_index,1), top_active(1), 192)', ... + linspace(curr_probePoints(tip_index,2), top_active(2), 192)', ... + linspace(curr_probePoints(tip_index,3), top_active(3), 192)']; + probe_contact_points = zeros(384, 3); + for i = 1:192 + probe_contact_points(2*i-1,:) = a(i,:); + probe_contact_points(2*i,:) = a(i,:); + end + clear a + + c_t_contacts = cell(384,1); + + theta = acos(p(2)); % Assuming p(2) corresponds to the DV direction + + for i = 1:384 + c_t_contacts{i} = apdvml2info(probe_contact_points(i,:), av, st, plane); + c_t_contacts{i}.contact_id = repmat(i,height(c_t_contacts{i})); + c_t_contacts{i}.probe_id = repmat(selected_probe,height(c_t_contacts{i})); + c_t_contacts{i}.depth_mm = tip2surface_mm - 0.020 * floor((i-1)/2); + + % Project the depth along the line + projected_depth_mm = c_t_contacts{i}.depth_mm * cos(theta); + + % Convert using the transformation for the chosen plane + projected_depth_mm_paxinos = accf2pxs_mm(projected_depth_mm, plane, 'distance'); + + c_t_contacts{i}.depth_mm_paxinos = projected_depth_mm_paxinos; + end + + c_Tapdvml_contacts{selected_probe} = vertcat(c_t_contacts{:}); + clear c_t_contacts + + pause(.05) +end + + +Tapdvml_contacts = vertcat(c_Tapdvml_contacts{:}); +Tapdvml_contacts.probe_AB = strings(height(Tapdvml_contacts), 1); +Tapdvml_contacts.session_id = strings(height(Tapdvml_contacts), 1); +Tapdvml_contacts.subject_id = strings(height(Tapdvml_contacts), 1); + + +T_borders = []; +for i = probes + borders_table{i}.probe_id = repmat(i,height(borders_table{i}),1); + T_borders = [T_borders; join(borders_table{i}, T_probes(i,:))]; +end + +i = T_borders.Properties.VariableNames == "upperBorder"; +T_borders.Properties.VariableNames{i} = 'upperBorder_um'; + +i = T_borders.Properties.VariableNames == "lowerBorder"; +T_borders.Properties.VariableNames{i} = 'lowerBorder_um'; + From 5d827e9a30db77b6cf71cb230f213c2b9a0e3371 Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Tue, 15 Aug 2023 14:28:33 +0100 Subject: [PATCH 21/43] getting there --- Browsing Functions/accf2pxs_mm.m | 96 ++++++++++++++++++++++++++------ Browsing Functions/apdvml2info.m | 14 ++++- 2 files changed, 91 insertions(+), 19 deletions(-) diff --git a/Browsing Functions/accf2pxs_mm.m b/Browsing Functions/accf2pxs_mm.m index ded57e6..37b107e 100644 --- a/Browsing Functions/accf2pxs_mm.m +++ b/Browsing Functions/accf2pxs_mm.m @@ -1,20 +1,84 @@ -function y_mm = accf2pxs_mm(x_mm) -% convert Allen Common Coordinates into standard values (Paxinos and Franklin) -% x_mm Allen Common Coordinates -% if this is scaler, x_mm should be depth -% if this is n x 3 array, x_mm is an array of the same (or similar) size -% y_mm Paxinos and Franklin coordinates +function y_mm = accf2pxs_mm(x_mm, plane, numtype) +% convert Allen Common Coordinates into standard values (Paxinos and Franklin, 2012, 4th edition) +% +% SYNTAX +% y_mm = accf2pxs_mm(x_mm, plane) +% +% INPUT ARGUMENTS +% x_mm Allen Common Coordinates +% if this is scaler, x_mm should be depth. +% if this is n x 3 array, x_mm is an array of the same (or similar) size +% +% plane 'coronal' | 'sagittal' | 'transverse' +% +% numtype 'coordinate' (defalut) | 'distance' +% +% OUTPUT ARGUMENTS +% y_mm Paxinos and Franklin coordinates +% +% Written by Kouichi C. Nakamura Ph.D. +% MRC Brain Network Dynamics Unit +% University of Oxford +% kouichi.c.nakamura@gmail.com +% 09-Aug-2023 11:06:59 +% +% See also +% apdvml2info -if size(x_mm,2) ==1 - y_mm = 0.921 * x_mm - 0.834; % global - % y = 0.908 * x - 0.957; % local +arguments + x_mm + plane (1,1) string {mustBeMember(plane, {'coronal','sagittal','transverse'})} + numtype (1,1) string {mustBeMember(numtype, {'coordinate','distance'})} = "coordinate" -elseif size(x_mm,2) == 3 - % assume that the 2nd column is the DV - y_mm = x_mm; - y_mm(:,2) = 0.921 * x_mm(:,2) - 0.834; % global - % y(:,2) = 0.908 * x(:,2) - 0.957; % local +end -else - error('unexpected size') + +switch plane + + case "coronal" + if size(x_mm,2) ==1 + switch numtype + case "coordinate" + y_mm = 0.825 * x_mm - 0.131; + case "distance" + y_mm = 0.825 * x_mm ; + end + + elseif size(x_mm,2) == 3 + + % assume that the 2nd column is the DV + y_mm = x_mm; + switch numtype + case "coordinate" + y_mm(:,2) = 0.825 * x_mm(:,2) - 0.131; + case "distance" + y_mm = 0.825 * x_mm(:,2) ; + end + else + error('unexpected size') + end + case "sagittal" + if size(x_mm,2) ==1 + switch numtype + case "coordinate" + y_mm = 0.921 * x_mm - 0.834; + case "distance" + y_mm = 0.921 * x_mm ; + end + elseif size(x_mm,2) == 3 + % assume that the 2nd column is the DV + y_mm = x_mm; + switch numtype + case "coordinate" + y_mm(:,2) = 0.921 * x_mm(:,2) - 0.834; + case "distance" + y_mm = 0.921 * x_mm(:,2) ; + end + else + error('unexpected size') + end + case "transverse" + error("not implemented yet") + otherwise + error("wrong input for 'plane'") end \ No newline at end of file diff --git a/Browsing Functions/apdvml2info.m b/Browsing Functions/apdvml2info.m index 8735ce2..4b88cc9 100644 --- a/Browsing Functions/apdvml2info.m +++ b/Browsing Functions/apdvml2info.m @@ -1,4 +1,4 @@ -function Tapdvml = apdvml2info(apdvml_points, av, st) +function Tapdvml = apdvml2info(apdvml_points, av, st, plane) % apdvml_points has three columns for ap, dv, and ml % bregma is 1 x 3 in size % atlas_resolution in mm @@ -6,7 +6,15 @@ % st % % see also -% Analyze_Clicked_Points.m +% Analyze_Clicked_Points.m, accf2pxs_mm + +arguments + apdvml_points (:,3) + av + st + plane (1,1) string {mustBeMember(plane, {'coronal','sagittal','transverse'})} +end + % generate needed values bregma = allenCCFbregma(); % bregma position in reference data space @@ -50,7 +58,7 @@ end -Tapdvml = [table(ap, dv, accf2pxs_mm(dv), ml, 'VariableNames',{'ap_mm','dv_mm', 'dv_mm_paxinos', 'ml_mm'}), ... +Tapdvml = [table(ap, dv, accf2pxs_mm(dv, plane), ml, 'VariableNames',{'ap_mm','dv_mm', 'dv_mm_paxinos', 'ml_mm'}), ... cell2table(roi_annotation_curr, 'VariableNames', {'annotation', 'name', 'acronym'})]; end \ No newline at end of file From 9099b37ffb3e36bcba576296312eb7d4b68844e6 Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Thu, 17 Aug 2023 15:24:54 +0100 Subject: [PATCH 22/43] mask_boundaries prep_Tavdvml_contacts --- Browsing Functions/mask_boundaries.m | 66 +++++++++++++++++ Browsing Functions/prep_Tapdvml_contacts.m | 84 ++++++++++++++++++++++ 2 files changed, 150 insertions(+) create mode 100644 Browsing Functions/mask_boundaries.m create mode 100644 Browsing Functions/prep_Tapdvml_contacts.m diff --git a/Browsing Functions/mask_boundaries.m b/Browsing Functions/mask_boundaries.m new file mode 100644 index 0000000..239132e --- /dev/null +++ b/Browsing Functions/mask_boundaries.m @@ -0,0 +1,66 @@ +function boundary_mask = mask_boundaries(av, margin_radius_um) +% mask_boundaries will return a 3D array of logical which defines the mask +% for structural boundaries +% +% SYNTAX +% mask_boundaries(av, margin_radius_um) +% +% INPUT ARGUMENTS +% av the 3D array provided from Allen Brain Institute as Allen CCF +% Grid interval is 10 um (micrometers) +% +% margin_radius_um +% positive number +% Boundary margin as a radius in um (micrometers) +% +% +% OUTPUT ARGUMENTS +% boundary_mask +% logical +% The same size as av. +% Voxels within margin_radius_um from boundaries are labelled +% as true. Otherwise the values are false. +% +% +% EXAMPLE +% boundary_mask = mask_boundaries(av, 20) +% +% Written by Kouichi C. Nakamura Ph.D. +% MRC Brain Network Dynamics Unit +% University of Oxford +% kouichi.c.nakamura@gmail.com +% 16-Aug-2023 11:22:57 +% +% See also +% doc + + +arguments + av (:, :, :) uint16 + margin_radius_um (1,1) {mustBePositive} +end + +% Initialize boundary_mask +boundary_mask = false(size(av)); + +% Convert margin_um to grid units +margin_units = round(margin_radius_um / 10); + +% Iterate through each voxel in av +for x = 1:size(av, 1) + tic + for y = 1:size(av, 2) + for z = 1:size(av, 3) + % Check neighbors in all three dimensions + neighbors = av(max(1, x-margin_units):min(size(av, 1), x+margin_units), ... + max(1, y-margin_units):min(size(av, 2), y+margin_units), ... + max(1, z-margin_units):min(size(av, 3), z+margin_units)); + + % If any neighbor has a different value, mark as boundary + if any(neighbors(:) ~= av(x, y, z)) + boundary_mask(x, y, z) = true; + end + end + end + toc +end \ No newline at end of file diff --git a/Browsing Functions/prep_Tapdvml_contacts.m b/Browsing Functions/prep_Tapdvml_contacts.m new file mode 100644 index 0000000..0f6e574 --- /dev/null +++ b/Browsing Functions/prep_Tapdvml_contacts.m @@ -0,0 +1,84 @@ +function Tapdvml_contacts = prep_Tapdvml_contacts(T_probes, Tapdvml_contacts) +% short description comes here +% +% SYNTAX +% Tapdvml_contacts = prep_Tapdvml_contacts(T_probes) +% +% longer description may come here +% +% INPUT ARGUMENTS +% T_probes table +% Rows for probes (+ optic fibers) +% With the following variables: +% probe_id +% probe_AB +% session_id +% subject_id +% probe_note +% upward_from_tip_um +% upward_from_tim_um +% +% Tapdvml_contacts +% table +% +% +% OUTPUT ARGUMENTS +% Tapdvml_contacts +% table +% +% Written by Kouichi C. Nakamura Ph.D. +% MRC Brain Network Dynamics Unit +% University of Oxford +% kouichi.c.nakamura@gmail.com +% 17-Aug-2023 09:51:42 +% +% See also +% doc + +arguments + T_probes table + Tapdvml_contacts table +end + +for i = 1:height(T_probes) + + tf = Tapdvml_contacts.probe_id == i; + Tapdvml_contacts.session_id(tf) = string(T_probes.session_id(i)); + Tapdvml_contacts.probe_AB(tf) = string(T_probes.probe_AB(i)); + Tapdvml_contacts.probe_note(tf) = string(T_probes.probe_note(i)); + +end + + +% Bregma's coordinates in grid units + +bregma = allenCCFbregma(); +bregma_ml = bregma(1); +bregma_dv = bregma(2); +bregma_ap = bregma(3); + +%load boundary_mask +load('\\ettina\Magill_lab\Kouichi Nakamura\Analysis\allenCCFdata\boundary_mask_20um_radius.mat', 'boundary_mask') % boundary_mask + +% Add a new column 'name_with_margin' to Tapdvml_contacts +Tapdvml_contacts.name_with_margin = Tapdvml_contacts.name; + +% Iterate through each row in Tapdvml_contacts +for i = 1:height(Tapdvml_contacts) + % Convert the coordinates to grid units and adjust based on bregma + ml = round(Tapdvml_contacts.ml_mm(i) * 100 + bregma_ml); + dv = round(Tapdvml_contacts.dv_mm(i) * 100 + bregma_dv); + ap = round(Tapdvml_contacts.ap_mm(i) * 100 + bregma_ap); + + if dv < 0 + dv = 1; + end + + if ap < 0 + ap = 1; + end + % Check the boundary_mask array + if boundary_mask(ml, dv, ap) + Tapdvml_contacts.name_with_margin(i) = {'boundary'}; + end +end From 47d559738ef4433afd95126d3308901c3afaa6c9 Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Fri, 6 Oct 2023 11:05:49 +0100 Subject: [PATCH 23/43] removing fwireframe2 --- .../plot_and_compute_probe_positions_from2.m | 46 ++++++++++++------- Browsing Functions/prep_Tapdvml_contacts.m | 4 +- 2 files changed, 32 insertions(+), 18 deletions(-) diff --git a/Browsing Functions/plot_and_compute_probe_positions_from2.m b/Browsing Functions/plot_and_compute_probe_positions_from2.m index 80475a7..312c928 100644 --- a/Browsing Functions/plot_and_compute_probe_positions_from2.m +++ b/Browsing Functions/plot_and_compute_probe_positions_from2.m @@ -1,4 +1,4 @@ -function [probes, fwireframe, fwireframe2, fig_probes, T_probes, Tapdvml_contacts, T_borders] ... +function [probes, fwireframe, fig_probes, T_probes, Tapdvml_contacts, T_borders] ... = plot_and_compute_probe_positions_from2(av, st, subject_id, plane,... processed_images_folder, probe_save_name_suffix,probes_to_analyze, probe_lengths,... active_probe_length, probe_radius) @@ -35,20 +35,38 @@ % % probes_to_analyze % 'all' | a row vector of positive integers specifying probes -% - -% B 0 (default) | non-negative integers -% (Optional) Description about B comes here. -% % -% OPTIONAL PARAMETER/VALUE PAIRS -% 'C' 'on' (default) | 'off' -% (Optional) Description about 'C' comes here. +% probe_lengths +% vector +% A vector of the depths of the deepest points for all the +% probes. Previously in the original SharpTrack, this +% information was used to calculate the depth, but now this is +% much less important. +% +% active_probe_length +% Should be 3.84 mm for Neuropixels 1.0 probes +% +% probe_radius scalar +% in micrometers % % % OUTPUT ARGUMENTS -% D cell array -% Description about D comes here. +% probes +% +% fwireframe figure +% +% fig_probes figure +% +% +% T_probes table +% Summary for probe penetrations +% +% Tapdvml_contacts +% table +% A table containig whereabouts of all the recording contacts +% +% T_borders table +% Modified from original SharpTrack output % % Written by Kouichi C. Nakamura Ph.D. % MRC Brain Network Dynamics Unit @@ -57,7 +75,7 @@ % 09-Aug-2023 11:12:42 % % See also -% doc +% apdvml2info @@ -127,10 +145,6 @@ hold on; fwireframe.InvertHardcopy = 'off'; -fwireframe2 = plotBrainGrid([], [], fwireframe, black_brain); %TODO what is this one for????????? -hold on; -fwireframe2.InvertHardcopy = 'off'; - fig_probes = gobjects(1,length(probes)); borders_table = cell(1,length(probes)) %TODO t_size = [size(probePoints.pointList.pointList,1),1]; diff --git a/Browsing Functions/prep_Tapdvml_contacts.m b/Browsing Functions/prep_Tapdvml_contacts.m index 0f6e574..17acdb0 100644 --- a/Browsing Functions/prep_Tapdvml_contacts.m +++ b/Browsing Functions/prep_Tapdvml_contacts.m @@ -70,11 +70,11 @@ dv = round(Tapdvml_contacts.dv_mm(i) * 100 + bregma_dv); ap = round(Tapdvml_contacts.ap_mm(i) * 100 + bregma_ap); - if dv < 0 + if dv <= 0 dv = 1; end - if ap < 0 + if ap <= 0 ap = 1; end % Check the boundary_mask array From 9d6ddba7b637e4851adda58f5b5f4720791e2835 Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Fri, 6 Oct 2023 11:06:04 +0100 Subject: [PATCH 24/43] Checking the status --- Browsing Functions/sliceBrowser.m | 49 ++++++++++++++++++++++++-- Histology Functions/HistologyBrowser.m | 44 ++++++++++++++++++++++- Histology Functions/SliceFlipper.m | 45 +++++++++++++++++++++-- 3 files changed, 133 insertions(+), 5 deletions(-) diff --git a/Browsing Functions/sliceBrowser.m b/Browsing Functions/sliceBrowser.m index d8bc290..a875484 100644 --- a/Browsing Functions/sliceBrowser.m +++ b/Browsing Functions/sliceBrowser.m @@ -27,6 +27,8 @@ function sliceBrowser(slice_figure, processed_images_folder, f, reference_size) 'String', '0 point', 'EdgeColor', 'none', 'Color', 'k', 'HorizontalAlignment', 'right'); ud_slice.pointsText.Visible = 'off'; +check_file_status(ud_slice) + % create functions needed to interact with the figure set(ud_slice.im, 'ButtonDownFcn', @(slice_figure,k)sliceClickCallback(slice_figure, k)); set(slice_figure, 'KeyPressFcn', @(slice_figure, keydata)SliceAtlasHotkeyFcn(slice_figure, keydata, f)); @@ -82,6 +84,7 @@ function sliceClickCallback(im, keydata) end set(f, 'UserData', ud); +end % ------------------------ % react to keyboard press @@ -151,7 +154,7 @@ function SliceAtlasHotkeyFcn(fig, keydata, f) set(fig, 'UserData', ud); - +end function ud = updateSliceImage(ud) %TODO occasionally pointList has more items than pointHands @@ -202,4 +205,46 @@ function SliceAtlasHotkeyFcn(fig, keydata, f) set(ud.pointHands(:), 'Visible', 'off'); %Hide because not in 't' mode ud.pointsText.String = sprintf('%d point(s)', length(ud.pointHands)); end - title(['Slice Viewer -- Slice ' num2str(ud.slice_num) '/' num2str(ud.total_num_files) title_ending]) + title(['Slice Viewer -- Slice ' num2str(ud.slice_num) '/' num2str(ud.total_num_files) title_ending]) + +end + +end + + +function check_file_status(ud_slice) + +proc_file_names = [string({ud_slice.processed_images.name})]'; +processed_images_folder = ud_slice.processed_images_folder; + + +T_files = table(proc_file_names, zeros(length(proc_file_names),1),'VariableNames',{'file_names','status'}); + +if isfile(fullfile(processed_images_folder, 'file_status.mat')) + + S = load(fullfile(processed_images_folder, 'file_status.mat')); + for i = 1:height(T_files) + % apply saved values for the existing images + T_files.status(i) = S.T_files{S.T_files.file_names == T_files.file_names(i) , 2}; + end + +else + error('You have to run HistologyBrowser first for the image %s', this_image_name) + +end + +T_files.Properties.VariableDescriptions = ["",... + "0, raw; 1, HistologyBrowser; 2, SliceFlipper"]; + +for i = 1:height(T_files) + switch T_files.status(i) + case 0 + error('You have to run HistologyBrowser and SliceFlipper for the image %s', this_image_name) + case 1 + error('You have to run SliceFlipper for the image %s', this_image_name) + case 2 + % good to go! + end +end + +end diff --git a/Histology Functions/HistologyBrowser.m b/Histology Functions/HistologyBrowser.m index c8058f4..319e8a1 100644 --- a/Histology Functions/HistologyBrowser.m +++ b/Histology Functions/HistologyBrowser.m @@ -41,6 +41,7 @@ function HistologyBrowser(histology_figure, save_folder, image_folder, image_fil ud.adjusted_image = original_image; imwrite(ud.adjusted_image, fullfile(ud.save_folder, [image_file_names{ud.file_num}(1:end-4) ud.file_name_suffix '.tif'])) +update_file_status(ud, image_file_names, ud.file_num, 1); set(histology_figure, 'UserData', ud); @@ -96,11 +97,13 @@ function HistologyHotkeyFcn(fig, keydata, image_file_names, use_already_downsamp case 's' % save image disp(['saving processed image ' num2str(ud.file_num)]); imwrite(ud.adjusted_image, fullfile(ud.save_folder, [image_file_names{ud.file_num}(1:end-4) ud.file_name_suffix '.tif'])) + update_file_status(ud, image_file_names, ud.file_num, 1); imshow(ud.adjusted_image) case 'leftarrow' % save image and move to previous image disp(['saving processed image ' num2str(ud.file_num)]); imwrite(ud.adjusted_image, fullfile(ud.save_folder, [image_file_names{ud.file_num}(1:end-4) ud.file_name_suffix '.tif'])) - + update_file_status(ud, image_file_names, ud.file_num, 1); + if ud.file_num > 1 ud.file_num = ud.file_num - 1; move_on = true; @@ -110,6 +113,7 @@ function HistologyHotkeyFcn(fig, keydata, image_file_names, use_already_downsamp case 'rightarrow' % save image and move to next image disp(['saving processed image ' num2str(ud.file_num)]); imwrite(ud.adjusted_image, fullfile(ud.save_folder, [image_file_names{ud.file_num}(1:end-4) ud.file_name_suffix '.tif'])) + update_file_status(ud, image_file_names, ud.file_num, 1); if ud.file_num < ud.num_files; ud.file_num = ud.file_num + 1; @@ -138,6 +142,7 @@ function HistologyHotkeyFcn(fig, keydata, image_file_names, use_already_downsamp % save immediately imwrite(ud.adjusted_image, fullfile(ud.save_folder, [image_file_names{ud.file_num}(1:end-4) ud.file_name_suffix '.tif'])) + update_file_status(ud, image_file_names, ud.file_num, 1); end else % if pressing commands while adjusting contrast @@ -158,4 +163,41 @@ function HistologyHotkeyFcn(fig, keydata, image_file_names, use_already_downsamp set(fig, 'UserData', ud); +end + +function update_file_status(ud, image_file_names, K, status) + + proc_file_names = strings(length(image_file_names),1); + for i = 1:length(image_file_names) + proc_file_names{i} = [image_file_names{i}(1:end-4) ud.file_name_suffix '.tif']; + end + T_files = table(proc_file_names, zeros(length(image_file_names),1),'VariableNames',{'file_names','status'}); + + if isfile(fullfile(ud.save_folder, 'file_status.mat')) + + S = load(fullfile(ud.save_folder, 'file_status.mat')); + for i = 1:height(T_files) + % apply saved values for the existing images + if nnz(S.T_files.file_names == string(T_files.file_names{i})) == 1 + T_files.status(i) = S.T_files{S.T_files.file_names == string(T_files.file_names{i}) , 2}; + elseif nnz(S.T_files.file_names == string(T_files.file_names{i})) > 1 + warning('multiple hits') + else + % ignore + end + end + + end + T_files.Properties.VariableDescriptions = ["",... + "0, raw; 1, HistologyBrowser; 2, SliceFlipper"]; + + T_files.status(K) = status; + + save(fullfile(ud.save_folder, 'file_status.mat'),'T_files') + +end + + + +end \ No newline at end of file diff --git a/Histology Functions/SliceFlipper.m b/Histology Functions/SliceFlipper.m index c8ba7e8..0f75ea6 100644 --- a/Histology Functions/SliceFlipper.m +++ b/Histology Functions/SliceFlipper.m @@ -52,12 +52,14 @@ function SliceCropHotkeyFcn(keydata, slice_figure, folder_processed_images) switch lower(keydata.Key) case 'leftarrow' % save and previous slice - imwrite(ud.current_slice_image, fullfile(folder_processed_images, ud.processed_image_name)) + imwrite(ud.current_slice_image, fullfile(folder_processed_images, ud.processed_image_name)) + update_file_status(folder_processed_images, ud, 2) %TODO ud.slice_num = ud.slice_num - 1*(ud.slice_num>1); ud = load_next_slice(ud,folder_processed_images); case 'rightarrow' % save and next slice imwrite(ud.current_slice_image, fullfile(folder_processed_images, ud.processed_image_name)) + update_file_status(folder_processed_images, ud, 2) %TODO ud.slice_num = ud.slice_num + 1*(ud.slice_num < length(ud.processed_image_names)); ud = load_next_slice(ud,folder_processed_images); @@ -136,6 +138,8 @@ function SliceCropHotkeyFcn(keydata, slice_figure, folder_processed_images) set(slice_figure, 'UserData', ud); +end + function ud = load_next_slice(ud,folder_processed_images) ud.processed_image_name = ud.processed_image_names{ud.slice_num}; ud.current_slice_image = imread(fullfile(folder_processed_images, ud.processed_image_name)); @@ -157,7 +161,8 @@ function SliceCropHotkeyFcn(keydata, slice_figure, folder_processed_images) ud.rotate_angle = 0; imwrite(ud.current_slice_image, fullfile(folder_processed_images, ud.processed_image_name)) - + update_file_status(folder_processed_images, ud, 2) %TODO +end % function to rotate slice by scrolling function SliceScrollFcn(fig, evt) @@ -173,5 +178,41 @@ function SliceScrollFcn(fig, evt) set(fig, 'UserData', ud); +end + +end + + +function update_file_status(folder_processed_images, ud, status) +proc_file_names = ud.processed_image_names; +this_image_name = ud.processed_image_name; +K = find(proc_file_names == string(this_image_name)); + +T_files = table(proc_file_names', zeros(length(proc_file_names),1),'VariableNames',{'file_names','status'}); + +if isfile(fullfile(folder_processed_images, 'file_status.mat')) + + S = load(fullfile(folder_processed_images, 'file_status.mat')); + for i = 1:height(T_files) + % apply saved values for the existing images + T_files.status(i) = S.T_files{S.T_files.file_names == string(T_files.file_names{i}) , 2}; + end + +end +T_files.Properties.VariableDescriptions = ["",... + "0, raw; 1, HistologyBrowser; 2, SliceFlipper"]; + +switch T_files.status(K) + case 0 + error('You have to run HistologyBrowser first for the image %s', this_image_name) + case 1 + T_files.status(K) = status; + case 2 + T_files.status(K) = status; +end + +save(fullfile(folder_processed_images, 'file_status.mat'),'T_files') + +end From 3e71d98dfe6c0d8f22a35dbcdd656ac17de232f2 Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Mon, 9 Oct 2023 14:45:46 +0100 Subject: [PATCH 25/43] file status checker behaviour updated If all the images have been processed, MATLAB will ask you whether to proceed (Yes) or skip (No). --- Histology Functions/HistologyBrowser.m | 33 ++++++++++++++++++++++- Histology Functions/SliceFlipper.m | 36 +++++++++++++++++++++++--- 2 files changed, 65 insertions(+), 4 deletions(-) diff --git a/Histology Functions/HistologyBrowser.m b/Histology Functions/HistologyBrowser.m index 319e8a1..104f698 100644 --- a/Histology Functions/HistologyBrowser.m +++ b/Histology Functions/HistologyBrowser.m @@ -33,6 +33,14 @@ function HistologyBrowser(histology_figure, save_folder, image_folder, image_fil ud.channel = min( 3, size(image,3)); original_image = image(:,:,1:ud.channel)*gain; +out = update_file_status(ud, image_file_names, 0, 0); +if out == 0 + fprintf("Aborted\n") + close(); + return +end + + imshow(original_image); title(['Adjusting channel ' num2str(ud.channel) ' on image ' num2str(ud.file_num) ' / ' num2str(ud.num_files)],... 'color',[1==ud.channel 2==ud.channel 3==ud.channel]) @@ -165,7 +173,7 @@ function HistologyHotkeyFcn(fig, keydata, image_file_names, use_already_downsamp end -function update_file_status(ud, image_file_names, K, status) +function out = update_file_status(ud, image_file_names, K, status) proc_file_names = strings(length(image_file_names),1); for i = 1:length(image_file_names) @@ -191,10 +199,33 @@ function update_file_status(ud, image_file_names, K, status) T_files.Properties.VariableDescriptions = ["",... "0, raw; 1, HistologyBrowser; 2, SliceFlipper"]; + if K == 0 % check this for the first run only + if ~any(T_files.status == 0) + answer = questdlg("All the files have been processed by HistologyBrowser. Would you like to do it again?",... + "HistologyBrowser",... + "Yes", "No", "No"); + + switch answer + case "Yes" + % proceed + out = 1; + return + case "No" + % quit the function + out = 0; + return + end + end + + out = 1; + return + end + T_files.status(K) = status; save(fullfile(ud.save_folder, 'file_status.mat'),'T_files') + out = 1; end diff --git a/Histology Functions/SliceFlipper.m b/Histology Functions/SliceFlipper.m index 0f75ea6..78891dd 100644 --- a/Histology Functions/SliceFlipper.m +++ b/Histology Functions/SliceFlipper.m @@ -21,7 +21,15 @@ function SliceFlipper(slice_figure, folder_processed_images, reference_size) title(['Slice ' num2str(ud.slice_num) ' / ' num2str(ud.total_num_files)]) set(slice_figure, 'UserData', ud); - + +out = update_file_status(folder_processed_images, ud, 0); %TODO +if out == 0 + fprintf("Aborted\n") + close(); + return +end + + % key function for slice set(slice_figure, 'KeyPressFcn', @(slice_figure,keydata)SliceCropHotkeyFcn(keydata, slice_figure, folder_processed_images)); @@ -183,7 +191,7 @@ function SliceScrollFcn(fig, evt) end -function update_file_status(folder_processed_images, ud, status) +function out = update_file_status(folder_processed_images, ud, status) proc_file_names = ud.processed_image_names; this_image_name = ud.processed_image_name; @@ -204,6 +212,28 @@ function update_file_status(folder_processed_images, ud, status) T_files.Properties.VariableDescriptions = ["",... "0, raw; 1, HistologyBrowser; 2, SliceFlipper"]; +if status == 0 % check this for the first run only + if all(T_files.status == 2) + answer = questdlg("All the files have been processed by SliceFlipper. Would you like to do it again?",... + "HistologyBrowser",... + "Yes", "No", "No"); + + switch answer + case "Yes" + % proceed + out = 1; + return + case "No" + % quit the function + out = 0; + return + end + end + + out = 1; + return +end + switch T_files.status(K) case 0 error('You have to run HistologyBrowser first for the image %s', this_image_name) @@ -214,5 +244,5 @@ function update_file_status(folder_processed_images, ud, status) end save(fullfile(folder_processed_images, 'file_status.mat'),'T_files') - +out = 1; end From 0c9123447ff9aec667e0dc128cf63211c3354f3a Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Tue, 10 Oct 2023 10:58:53 +0100 Subject: [PATCH 26/43] last point is set green when deleting the points --- Browsing Functions/sliceBrowser.m | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Browsing Functions/sliceBrowser.m b/Browsing Functions/sliceBrowser.m index a875484..1f5a94a 100644 --- a/Browsing Functions/sliceBrowser.m +++ b/Browsing Functions/sliceBrowser.m @@ -128,6 +128,8 @@ function SliceAtlasHotkeyFcn(fig, keydata, f) phy(i) = ud.pointHands(i).YData; end ud.pointList = [phx, ud.ref_size(1) - phy]; + set(ud.pointHands(end),'color', [0 .9 0]); + disp('transform point deleted') ud.pointsText.String = sprintf('%d point(s)', length(ud.pointHands)); From bdcd3fb7649c31a1b0c8d7eb56a8a23183be205f Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Tue, 10 Oct 2023 15:25:10 +0100 Subject: [PATCH 27/43] minor bug fix --- Browsing Functions/sliceBrowser.m | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Browsing Functions/sliceBrowser.m b/Browsing Functions/sliceBrowser.m index 1f5a94a..45cf0aa 100644 --- a/Browsing Functions/sliceBrowser.m +++ b/Browsing Functions/sliceBrowser.m @@ -128,8 +128,9 @@ function SliceAtlasHotkeyFcn(fig, keydata, f) phy(i) = ud.pointHands(i).YData; end ud.pointList = [phx, ud.ref_size(1) - phy]; - set(ud.pointHands(end),'color', [0 .9 0]); - + if ~isempty(ud.pointHands) + set(ud.pointHands(end),'color', [0 .9 0]); + end disp('transform point deleted') ud.pointsText.String = sprintf('%d point(s)', length(ud.pointHands)); From eeb4136b900bb821a9c48926f510b46fd04580b8 Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Mon, 30 Oct 2023 16:10:21 +0000 Subject: [PATCH 28/43] renamed the function --- ...ot_and_compute_probe_positions_from_ccf.m} | 31 ++++++++++++------- 1 file changed, 19 insertions(+), 12 deletions(-) rename Browsing Functions/{plot_and_compute_probe_positions_from2.m => plot_and_compute_probe_positions_from_ccf.m} (93%) diff --git a/Browsing Functions/plot_and_compute_probe_positions_from2.m b/Browsing Functions/plot_and_compute_probe_positions_from_ccf.m similarity index 93% rename from Browsing Functions/plot_and_compute_probe_positions_from2.m rename to Browsing Functions/plot_and_compute_probe_positions_from_ccf.m index 312c928..9e1f441 100644 --- a/Browsing Functions/plot_and_compute_probe_positions_from2.m +++ b/Browsing Functions/plot_and_compute_probe_positions_from_ccf.m @@ -1,5 +1,5 @@ -function [probes, fwireframe, fig_probes, T_probes, Tapdvml_contacts, T_borders] ... - = plot_and_compute_probe_positions_from2(av, st, subject_id, plane,... +function [probes, fwireframe, fig_probes, Tapdvml_contacts, T_borders, p] ... + = plot_and_compute_probe_positions_from_ccf(av, st, subject_id, plane,... processed_images_folder, probe_save_name_suffix,probes_to_analyze, probe_lengths,... active_probe_length, probe_radius) % Create 3D plot for probe positions, 2D plots for brain structures along @@ -7,7 +7,7 @@ % % % SYNTAX -% [] = plot_and_compute_probe_positions(fwireframe, black_brain, probes, subject_id, plane,... +% [] = plot_and_compute_probe_positions_from_ccf(fwireframe, black_brain, probes, subject_id, plane,... % use_tip_to_get_reference_probe_length, probe_insertion_direction) % [D] = func1(A,B) % [D] = func1(____,'Param',value) @@ -56,10 +56,6 @@ % fwireframe figure % % fig_probes figure -% -% -% T_probes table -% Summary for probe penetrations % % Tapdvml_contacts % table @@ -68,6 +64,9 @@ % T_borders table % Modified from original SharpTrack output % +% p array with 3 columns +% eigen vectors for the largest eigen values +% % Written by Kouichi C. Nakamura Ph.D. % MRC Brain Network Dynamics Unit % University of Oxford @@ -123,7 +122,11 @@ % load probe points probePoints = load(fullfile(processed_images_folder, ['probe_points' probe_save_name_suffix])); ProbeColors = .75*[1.3 1.3 1.3; 1 .75 0; .3 1 1; .4 .6 .2; 1 .35 .65; .7 .7 .9; .65 .4 .25; .7 .95 .3; .7 0 0; .6 0 .7; 1 .6 0]; -% order of colors: {'white','gold','turquoise','fern','bubble gum','overcast sky','rawhide', 'green apple','red','purple','orange'}; + +ProbeColors255 = uint8(255 * ProbeColors); +order_of_colors = ["white","gold","turquoise","fern","bubble gum","overcast sky","rawhide", "green apple","red","purple","orange"]'; +T_colors = table(order_of_colors, ProbeColors255); % useful when managing multiple probes + fwireframe = []; % % scale active_probe_length appropriately @@ -146,8 +149,9 @@ fwireframe.InvertHardcopy = 'off'; fig_probes = gobjects(1,length(probes)); -borders_table = cell(1,length(probes)) %TODO +borders_table = cell(1,length(probes)); %TODO t_size = [size(probePoints.pointList.pointList,1),1]; + T_probes = table(repmat(string(subject_id), t_size),... strings(t_size), ... NaN(t_size),... @@ -157,12 +161,13 @@ NaN(t_size),... NaN(t_size),... NaN(t_size),... - NaN(t_size),... strings(t_size), ... 'VariableNames',{'subject_id','session_id','probe_id','probe_AB',... 'depth_reported_um','depth_from_ACCF_um','scaling', ... - 'scaling2',... - 'tip_active_um', 'top_active_um', 'probe_note'}); + 'tip_active_um', 'top_active_um', 'probe_note', ... + }); +P = zeros([size(probePoints.pointList.pointList,1),3]); + T_probes.probe_id(probes') = probes'; T_probes.probe_lengths(probes) = probe_lengths(probes)' * 1000; @@ -260,6 +265,7 @@ T_probes.depth_from_ACCF_um(selected_probe) = reference_probe_length_tip/100 * 1000; T_probes.scaling(selected_probe) = shrinkage_factor; + P(selected_probe,:) = p; % eigenvector for largest eigenvalue % plot line the length of the probe in reference space probe_length_histo = round(reference_probe_length_tip); @@ -388,3 +394,4 @@ i = T_borders.Properties.VariableNames == "lowerBorder"; T_borders.Properties.VariableNames{i} = 'lowerBorder_um'; +p = P; %rename \ No newline at end of file From 347b83703cfbfdc79ed8a006778fcbef673b6b9b Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Mon, 6 Nov 2023 15:04:11 +0000 Subject: [PATCH 29/43] p and m --- Browsing Functions/apdvml2info.m | 30 +++++++ Browsing Functions/apdvml_mm2ccf.m | 90 +++++++++++++++++++ ...lot_and_compute_probe_positions_from_ccf.m | 7 +- 3 files changed, 125 insertions(+), 2 deletions(-) create mode 100644 Browsing Functions/apdvml_mm2ccf.m diff --git a/Browsing Functions/apdvml2info.m b/Browsing Functions/apdvml2info.m index 4b88cc9..2ddb266 100644 --- a/Browsing Functions/apdvml2info.m +++ b/Browsing Functions/apdvml2info.m @@ -5,6 +5,36 @@ % av 1320 ,800, 1140 % st % +% SYNTAX +% Tapdvml = apdvml2info(apdvml_points, av, st, plane) +% +% INPUT ARGUMENTS +% apdvml_points +% n by 3 array +% [ap, dv, ml] +% Allen CCF coordinates in 10 micrometers. +% +% av uint16 (1320 x 800 x 1140) +% 3D array for structure delineation +% +% st table +% For structures +% +% plane 'coronal' | 'sagittal' | 'transverse' +% +% OUTPUT ARGUMENTS +% Tapdvml table +% Including output ap_mm, dv_mm, and ml_mm +% +% Written by Kouichi C. Nakamura Ph.D. +% MRC Brain Network Dynamics Unit +% University of Oxford +% kouichi.c.nakamura@gmail.com +% 06-Nov-2023 14:41:14 +% +% See also +% apdvml_mm2ccf + % see also % Analyze_Clicked_Points.m, accf2pxs_mm diff --git a/Browsing Functions/apdvml_mm2ccf.m b/Browsing Functions/apdvml_mm2ccf.m new file mode 100644 index 0000000..427af1f --- /dev/null +++ b/Browsing Functions/apdvml_mm2ccf.m @@ -0,0 +1,90 @@ +function Tccf = apdvml_mm2ccf(apdvml_mm)%, av, st, plane) +% Opposite of apdvml2info. +% Convert coordinates in mm back to Allen CCF coodinates in 10 micrometers +% +% SYNTAX +% v = apdvml_mm2ccf(apdvml_mm) +% +% +% INPUT ARGUMENTS +% apdvml_mm n by 3 array +% [ap_mm, dv_mm, ml_mm] +% a coordinate in mm +% +% % av uint16 (1320 x 800 x 1140) +% % 3D array for structure delineation +% % +% % st table +% % For structures +% % +% % plane 'coronal' | 'sagittal' | 'transverse' +% +% OUTPUT ARGUMENTS +% Tccf table +% including Allen CCF coordinates in 10 micrometers. +% +% Written by Kouichi C. Nakamura Ph.D. +% MRC Brain Network Dynamics Unit +% University of Oxford +% kouichi.c.nakamura@gmail.com +% 06-Nov-2023 14:35:03 +% +% See also +% apdvml2info + + +arguments + apdvml_mm (:,3) + % av + % st + % plane (1,1) string {mustBeMember(plane, {'coronal','sagittal','transverse'})} +end + +% generate needed values +bregma = allenCCFbregma(); % bregma position in reference data space +atlas_resolution = 0.010; % mm + + +%TODO + +ap = (- apdvml_mm(:,1))/atlas_resolution + bregma(1); +dv = (apdvml_mm(:,2))/atlas_resolution + bregma(2); +ml = (apdvml_mm(:,3))/atlas_resolution + bregma(3); + +% roi_location_curr = [ap dv ml]; + +% % initialize array of region annotations +% roi_annotation_curr = cell(size(apdvml_points,1),3); +% +% % loop through every point to get ROI locations and region annotations +% tf_rows = true(size(apdvml_mm,1),1); +% for point = 1:size(apdvml_mm,1) +% +% % find the annotation, name, and acronym of the current ROI pixel +% depth = ceil(apdvml_points(point,2)); +% if depth >= 1 +% +% ann = av(ceil(apdvml_points(point,1)), ... +% depth, ... %NOTE this can take the value of 0 and cause an error (must be >= 1) +% ceil(apdvml_points(point,3))); +% name = st.safe_name{ann}; +% acr = st.acronym{ann}; +% +% roi_annotation_curr{point,1} = ann; +% roi_annotation_curr{point,2} = name; +% roi_annotation_curr{point,3} = acr; +% +% else +% % ignore this point, too dorsal +% roi_annotation_curr{point,1} = NaN; +% roi_annotation_curr{point,2} = ''; +% roi_annotation_curr{point,3} = ''; +% +% end +% +% end + +Tccf = table(ap, dv, ml, 'VariableNames',{'ap_mm','dv_mm', 'ml_mm'});%, ... + %cell2table(roi_annotation_curr, 'VariableNames', {'annotation', 'name', 'acronym'})]; + +end \ No newline at end of file diff --git a/Browsing Functions/plot_and_compute_probe_positions_from_ccf.m b/Browsing Functions/plot_and_compute_probe_positions_from_ccf.m index 9e1f441..16b6996 100644 --- a/Browsing Functions/plot_and_compute_probe_positions_from_ccf.m +++ b/Browsing Functions/plot_and_compute_probe_positions_from_ccf.m @@ -1,4 +1,4 @@ -function [probes, fwireframe, fig_probes, Tapdvml_contacts, T_borders, p] ... +function [probes, fwireframe, fig_probes, Tapdvml_contacts, T_borders, p, m] ... = plot_and_compute_probe_positions_from_ccf(av, st, subject_id, plane,... processed_images_folder, probe_save_name_suffix,probes_to_analyze, probe_lengths,... active_probe_length, probe_radius) @@ -167,6 +167,7 @@ 'tip_active_um', 'top_active_um', 'probe_note', ... }); P = zeros([size(probePoints.pointList.pointList,1),3]); +M = zeros([size(probePoints.pointList.pointList,1),3]); T_probes.probe_id(probes') = probes'; T_probes.probe_lengths(probes) = probe_lengths(probes)' * 1000; @@ -266,6 +267,7 @@ T_probes.depth_from_ACCF_um(selected_probe) = reference_probe_length_tip/100 * 1000; T_probes.scaling(selected_probe) = shrinkage_factor; P(selected_probe,:) = p; % eigenvector for largest eigenvalue + M(selected_probe,:) = m; % highest point % plot line the length of the probe in reference space probe_length_histo = round(reference_probe_length_tip); @@ -394,4 +396,5 @@ i = T_borders.Properties.VariableNames == "lowerBorder"; T_borders.Properties.VariableNames{i} = 'lowerBorder_um'; -p = P; %rename \ No newline at end of file +p = P; %rename +m = M; \ No newline at end of file From cc91f663b48268e6a47b284259483b60541fb211 Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Mon, 6 Nov 2023 15:53:24 +0000 Subject: [PATCH 30/43] bug fix --- Browsing Functions/prep_Tapdvml_contacts.m | 1 + 1 file changed, 1 insertion(+) diff --git a/Browsing Functions/prep_Tapdvml_contacts.m b/Browsing Functions/prep_Tapdvml_contacts.m index 17acdb0..aaf8924 100644 --- a/Browsing Functions/prep_Tapdvml_contacts.m +++ b/Browsing Functions/prep_Tapdvml_contacts.m @@ -44,6 +44,7 @@ tf = Tapdvml_contacts.probe_id == i; Tapdvml_contacts.session_id(tf) = string(T_probes.session_id(i)); + Tapdvml_contacts.subject_id(tf) = string(T_probes.subject_id(i)); Tapdvml_contacts.probe_AB(tf) = string(T_probes.probe_AB(i)); Tapdvml_contacts.probe_note(tf) = string(T_probes.probe_note(i)); From eaa68f355a339e2d536d52537adda83b94816d04 Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Mon, 13 Nov 2023 11:23:48 +0000 Subject: [PATCH 31/43] Create updown_probe_with_slider.m --- Browsing Functions/updown_probe_with_slider.m | 542 ++++++++++++++++++ 1 file changed, 542 insertions(+) create mode 100644 Browsing Functions/updown_probe_with_slider.m diff --git a/Browsing Functions/updown_probe_with_slider.m b/Browsing Functions/updown_probe_with_slider.m new file mode 100644 index 0000000..04468c8 --- /dev/null +++ b/Browsing Functions/updown_probe_with_slider.m @@ -0,0 +1,542 @@ +function updown_probe_with_slider(session_id, probeAB, plane, sessions_dir,... + task_dir_name, imaging_session_dir, image_folder_name, depth_level) + +arguments + session_id (1,1) string + probeAB (1,1) string {mustBeMember(probeAB,["ProbeA","ProbeB"])} + plane (1,1) string {mustBeMember(plane,["coronal","sagittal", "transverse"])} + sessions_dir (1,1) string {mustBeFolder} + task_dir_name (1,1) string + imaging_session_dir (1,1) string {mustBeFolder} + image_folder_name (1,1) string + depth_level (1,1) double {mustBePositive, mustBeInteger} = 6 +end + +sorter_output_dir = fullfile(sessions_dir, task_dir_name, ... + session_id, "processed\kilosort", probeAB, "sorter_output"); + +recording_cell_metrics_path = fullfile(sorter_output_dir, "recording.cell_metrics.cellinfo.mat"); +s = load(recording_cell_metrics_path); +recording_cell_metrics = s.cell_metrics; + +opts = detectImportOptions(Tapdvml_contacts_path); +opts = setvartype(opts, 'session_id', 'string'); + +Tapdvml_contacts = readtable(Tapdvml_contacts_path, opts ); + +% Tborders_table = readtable(borders_table_path); + +T_probes= readtable(Tprobes_path); + +if ~ ismember('p_1',T_probes.Properties.VariableNames) || ... + ~ ismember('p_2',T_probes.Properties.VariableNames) ||... + ~ ismember('p_3',T_probes.Properties.VariableNames) + warning('The eigen vectors p is not found in T_probes. You cannot update the Tapdvml_contacts.') +end + + +%% Figure + +ufig = uifigure(Position=[800, 150, 720, 1080]); + +uax1 = uiaxes(ufig); +uax1.PositionConstraint = 'innerposition'; +uax1.Position=[100, 150, 350, 800]; +ylabel(uax1, 'Dsitance from the tip (mm)') +uax1.XTick =[]; + +uax2 = uiaxes(ufig); +uax2.PositionConstraint = 'innerposition'; +uax2.Position=[448, 150, 80, 800]; +uax2.YTickLabel =[]; +uax2.XTick =[]; + + +uax3 = uiaxes(ufig); +uax3.PositionConstraint = 'innerposition'; +uax3.Position=[525, 150, 120, 800]; +uax3.YTickLabel =[]; +uax3.XTick =[]; + + +tickdir([uax1,uax2, uax3], 'out') +ticklengthcm([uax1,uax2, uax3], 0.2) + +% usl1 = uislider(ufig, Orientation='vertical' , Limits=[-50, 50]); +% usl1.Position = [50, 200, 3, 300]; + +bg1 = uibuttongroup(ufig, Title='Move'); +bg1.Position = [20, 150, 70, 200]; + +bt1 = uibutton(bg1, "Text","+5", "ButtonPushedFcn", @(src, event) move_vertically(uax1, 5)); +bt1.Position = [10, 140, 50, 30]; + +bt2 = uibutton(bg1, "Text","+1", "ButtonPushedFcn", @(src, event) move_vertically(uax1, 1)); +bt2.Position = [10, 110, 50, 30]; + +bt3 = uibutton(bg1, "Text","-1", "ButtonPushedFcn", @(src, event) move_vertically(uax1, -1)); +bt3.Position = [10, 40, 50, 30]; + +bt4 = uibutton(bg1, "Text","-5", "ButtonPushedFcn", @(src, event) move_vertically(uax1, -5)); +bt4.Position = [10, 10, 50, 30]; + +bt6 = uibutton(bg1, "Text","0", "ButtonPushedFcn", @(src, event) move_back_to_original(uax1)); +bt6.Position = [10, 75, 50, 30]; + + +bt7 = uibutton(ufig, "Text","Save", "ButtonPushedFcn", @(src, event) save_updown_to_table(uax1)); +bt7.Position = [150, 30, 300, 50]; + +tx1 = annotation(ufig, 'textbox', [0.2 0.85 0.5 0.1], 'String', session_id +"/" + probeAB,... + HorizontalAlignment='center',FontSize=18,LineStyle='none'); + + +%% Plot cells +n = length(recording_cell_metrics.firingRate); + +x_val = rand(n,1); + +y = recording_cell_metrics.maxWaveformCh * 0.010; % mm + +c = recording_cell_metrics.firingRate; + +sc1 = scatter(uax3, x_val, y, 36, c, LineWidth=1.5); +uax3.Color = 'k'; + +ylim(uax3, [-0.100, 3.940]) + +cb1 =colorbar(uax3,'eastoutside'); +ylabel(cb1,'Mean firing rates (spikes/s)' ) +cb1.Label.FontSize = 14; +cb1.TickDirection = 'out'; +ticklengthcm(cb1, 0.2) + + +%% Moving average of five channels +% +% SDF like plot with Gaussian kernel convolution is better +cla(uax2,'reset') + +movingaverageval = zeros(384,1); + +for i = 1:384 + + % moving average of five channels + chans = i - 1 : i + 3; + ind = ismember(recording_cell_metrics.maxWaveformCh, chans); + movingaverageval(i) = mean(recording_cell_metrics.firingRate(ind)); + +end + +isc1 = imagesc(uax2, 1, (1:384)'*0.01, movingaverageval); + +uax2.YDir = 'normal'; +ylim(uax2, [-0.100, 3.940]) +uax2.XTick =[]; + +uax2.Color = 'k'; + +uax3.CLim = uax2.CLim; % same color scale + + +%% Read the Allen CCF data +% directory of reference atlas files +annotation_volume_location = "\\ettina\Magill_lab\Kouichi Nakamura\Analysis\allenCCFdata\annotation_volume_10um_by_index.npy"; +structure_tree_location = "\\ettina\Magill_lab\Kouichi Nakamura\Analysis\allenCCFdata\structure_tree_safe_2017.csv"; +template_volume_location = "\\ettina\Magill_lab\Kouichi Nakamura\Analysis\allenCCFdata\template_volume_10um.npy"; + +% load the reference brain and region annotations +if ~exist('av','var') || ~exist('st','var') || ~exist('tv','var') + disp('loading reference atlas...') + av = readNPY(annotation_volume_location); + st = loadStructureTree(structure_tree_location); + tv = readNPY(template_volume_location); +end + + +probeAB = "ProbeA"; + +switch probeAB + case "ProbeA" + probeAB_short = "A"; + case "ProbeB" + probeAB_short = "B"; + otherwise + +end + +tf = T_probes.session_id == session_id & T_probes.probe_AB == probeAB_short; + +probe_id = T_probes{tf, "probe_id"}; + + +depth = 6; + +Tapdvml_contacts_this = Tapdvml_contacts(Tapdvml_contacts.probe_id == probe_id, :); + +%% +anc = preallocatestruct(["index", "anc_id", "anc_name", "anc_color_hex", "anc_color_255", "anc_color"], [height(Tapdvml_contacts_this), 1]); +% profile on +st.name = regexprep(st.name, '["'']',''); % remove clutters + +f = waitbar(0,'1','Name','Approximating pi...',... + 'CreateCancelBtn','setappdata(gcbf,''canceling'',1)'); + +for i = 1:height(Tapdvml_contacts_this) % SLOW + waitbar( i/height(Tapdvml_contacts_this), f, ... + fprintf("%d of %d", i, height(Tapdvml_contacts_this))) + name = Tapdvml_contacts_this.name{i}; + [anc(i)] = find_name_for_depth(st, depth_level, name); %TODO stopped working for i = 61, i = 1327 +end +% profile viewer + +delete(f) + +%% +Tanc = struct2table(anc, AsArray=true); + +Tapdvml_contacts_this = [Tapdvml_contacts_this, Tanc]; + + +% Sample data +anc_id = Tapdvml_contacts_this.anc_id; + +anc_name = Tapdvml_contacts_this.anc_name; + +anc_color = Tapdvml_contacts_this.anc_color; + + +anc_id_diff = diff(anc_id); + +start_ind = [1; find(anc_id_diff) + 1]; + +end_ind = [start_ind(2:end) - 1; length(anc_id)]; + +block_names = anc_name(start_ind); + +block_colors = anc_color(start_ind, :); + + +Tblocks = table(block_names, start_ind, end_ind, block_colors, ... + VariableNames={'block_names', 'start_ind','end_ind', 'block_colors'}); + + +%% Plot structures +cla(uax1) + +uax1.Color = 'k'; + +for i = 1:height(Tblocks) + X = [0 1 1 0]; + Y = [Tblocks{i,'start_ind'}-0.5, Tblocks{i,'start_ind'}-0.5,... + Tblocks{i,'end_ind'}+0.5, Tblocks{i,'end_ind'}+0.5]/100; % in mm + patch(uax1, X, Y, Tblocks.block_colors(i,:)); + text(uax1, 0.5, mean([Tblocks{i,'start_ind'}, Tblocks{i,'end_ind'}])/100, ... + Tblocks.block_names(i),... + 'FontSize',14,... + 'HorizontalAlignment','center', Clipping='on'); +end + +ylim(uax1, [-0.100, 3.940]) + +uax1.YLabel.FontSize = 14; + + + +function move_back_to_original(uax1) + +if isfield(uax1.UserData, 'YLim_shift_mm') + ch = uax1.Children; + + for i = 1:length(ch) + step = - uax1.UserData.YLim_shift_mm; + + if isa(ch(i),'matlab.graphics.primitive.Text') + ch(i).Position(2) = ch(i).Position(2) + step; + + elseif isa(ch(i),'matlab.graphics.primitive.Patch') + ch(i).YData = ch(i).YData + step; + + end + + end + + uax1.UserData.YLim_shift_mm = 0; +else + %nothing to do +end + +end + + +function move_vertically(uax1, step) + +% arguments +% uax1(1,1) +% step (1,1) {mustBeInteger} +% end + +ch = uax1.Children; + +for i = 1:length(ch) + + if isa(ch(i),'matlab.graphics.primitive.Text') + ch(i).Position(2) = ch(i).Position(2) + step*0.01; + + elseif isa(ch(i),'matlab.graphics.primitive.Patch') + ch(i).YData = ch(i).YData + step*0.01; + + end + +end + +if isfield(uax1.UserData, 'YLim_shift_mm') + uax1.UserData.YLim_shift_mm = uax1.UserData.YLim_shift_mm + step*0.01; +else + uax1.UserData.YLim_shift_mm = step*0.01; +end + +end + + + +function save_updown_to_table(uax1, ~) + +if isempty(uax1.UserData) + disp('The probe has not moved on GUI yet. No change made.') + return +else + if isfield(uax1.UserData, 'YLim_shift_mm') + %TODO move + else + disp('YLim_shift_mm is not found. No change made.') + return + end +end + + +%TODO save the backup + +assert(endsWith(Tapdvml_contacts_path,'.xlsx')) + +assert(isfile(Tapdvml_contacts_path)) + +f_info = dir(Tapdvml_contacts_path); + +dt = datetime(f_info.date, 'InputFormat', 'dd-MMM-yyyy HH:mm:ss'); + +dt.Format = 'yyyy-MM-dd_HH-mm-ss'; +suffix = string(dt); + +new_file_path = regexprep(Tapdvml_contacts_path, '\.xlsx$' , "_" + suffix + ".xlsx"); + +original_backup_path = regexprep(Tapdvml_contacts_path, '\.xlsx$' , "_" + suffix + "_original.xlsx"); + +if isfile(original_backup_path) + % save a backup file + + copyfile(Tapdvml_contacts_path, new_file_path) + +else + % this is the first time you change the Excel file + % save a backup file + + copyfile(Tapdvml_contacts_path, original_backup_path) + +end + +%% Job +see plot_and_compute_probe_positions_from_ccf.m +% I need the eigen vector, to be saved in T_probes? + +switch probeAB + case "ProbeA" + probeAB_ = "A"; + case "ProbeB" + probeAB_ = "B"; + case "optic fiber" + probeAB_ = "optic fiber"; +end + +% eigen vector +tf = T_probes.session_id == session_id & ... + T_probes.probe_AB == probeAB_; + +assert(nnz(tf) == 1) + +p = T_probes{T_probes.session_id == session_id & ... + T_probes.probe_AB == probeAB_, ["p_1","p_2","p_3"]}; + +m = T_probes{T_probes.session_id == session_id & ... + T_probes.probe_AB == probeAB_, ["m_1","m_2","m_3"]}; + +%TODO what columns to change or add??? + +processed_images_folder = fullfile(image_folder, 'processed'); + +probePoints = load(fullfile(processed_images_folder, ['probe_points' probe_save_name_suffix])); + +% get the probe points for the currently analyzed probe +if strcmp(plane,'coronal') + curr_probePoints = probePoints.pointList.pointList{selected_probe,1}(:, [3 2 1]); +elseif strcmp(plane,'sagittal') + curr_probePoints = probePoints.pointList.pointList{selected_probe,1}(:, [1 2 3]); +elseif strcmp(plane,'transverse') + curr_probePoints = probePoints.pointList.pointList{selected_probe,1}(:, [1 3 2]); +end + +% m is the brain entry point + +%TODO modifiy the m and curr_probePoints(tip_index,:) by YLim_shift_mm +% p is in 10 µm + +if strcmp(probe_insertion_direction, 'down') + [depth, tip_index] = max(curr_probePoints(:,2)); +elseif strcmp(probe_insertion_direction, 'up') + [depth, tip_index] = min(curr_probePoints(:,2)); +end + +% curr_probePoints(tip_index,:) either to get it from the ProbePoints or +% from Tapdvml_contacts.xlsx (convert ap_mm, dv_mm, ml_mm back to allen coordinates) +% opposite of apdvml2info + + +% Tapdvml_contacts.session_id + +% tip_mm = + +YLim_shift_mm = uax1.UserData.YLim_shift_mm; +new_cp = curr_probePoints(tip_index,:) - YLim_shift_mm*100; % minus to move up + + +Tapdvml_m = apdvml2info(m, av, st, plane); + +Tapdvml_tip = apdvml2info(new_cp, av, st, plane); + + + +% measure the distance + +tip2surface_mm = sqrt((Tapdvml_tip.ap_mm - Tapdvml_m.ap_mm)^2 + ... + (Tapdvml_tip.dv_mm - Tapdvml_m.dv_mm)^2 + ... + (Tapdvml_tip.ml_mm - Tapdvml_m.ml_mm)^2); + +tip2surface_mm_paxinos = sqrt((Tapdvml_tip.ap_mm - Tapdvml_m.ap_mm)^2 + ... + (Tapdvml_tip.dv_mm_paxinos - Tapdvml_m.dv_mm_paxinos)^2 + ... + (Tapdvml_tip.ml_mm - Tapdvml_m.ml_mm)^2); + +top_active = (m * active_probe_length + ... + curr_probePoints(tip_index,:) * (tip2surface_mm - active_probe_length))... + /tip2surface_mm; + +% obtained the information for all the 384 channels + + +a = [linspace(curr_probePoints(tip_index,1), top_active(1), 192)', ... + linspace(curr_probePoints(tip_index,2), top_active(2), 192)', ... + linspace(curr_probePoints(tip_index,3), top_active(3), 192)']; +probe_contact_points = zeros(384, 3); +for i = 1:192 + probe_contact_points(2*i-1,:) = a(i,:); + probe_contact_points(2*i,:) = a(i,:); +end +clear a + +c_t_contacts = cell(384,1); + +theta = acos(p(2)); % Assuming p(2) corresponds to the DV direction + +for i = 1:384 + c_t_contacts{i} = apdvml2info(probe_contact_points(i,:), av, st, plane); + c_t_contacts{i}.contact_id = repmat(i,height(c_t_contacts{i})); + c_t_contacts{i}.probe_id = repmat(selected_probe,height(c_t_contacts{i})); + c_t_contacts{i}.depth_mm = tip2surface_mm - 0.020 * floor((i-1)/2); + + % Project the depth along the line + projected_depth_mm = c_t_contacts{i}.depth_mm * cos(theta); + + % Convert using the transformation for the chosen plane + projected_depth_mm_paxinos = accf2pxs_mm(projected_depth_mm, plane, 'distance'); + + c_t_contacts{i}.depth_mm_paxinos = projected_depth_mm_paxinos; +end + +c_Tapdvml_contacts{selected_probe} = vertcat(c_t_contacts{:}); +clear c_t_contacts + +%TODO writetable() + +keyboard() + +end + +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function [anc, this] = find_name_for_depth(st, depth, name) + +arguments + st table + depth (1,1) double {mustBeInteger,mustBeNonnegative} + name (1,1) string +end + +st.name = regexprep(st.name, '["'']',''); + +id = st.id(st.name == name); + + +cstr = arrayfun(@(x) strsplit(x, "/"), string(st.structure_id_path),'UniformOutput',false); + +% remove empty string "" + +% cemp = cellfun(@(x) x == "" , cstr,'UniformOutput',false); + +for j = 1:size(cstr,1) % i = 741 + % fprintf("j = %d\n",j) + % cstr{i}(cemp{i}) = []; + % v = zeros(1,length(cstr{i}),'int32'); + % for j = 1:length(cstr{i}) + % fprintf("j = %d\n",j) + % v(j) = int32(str2double(cstr{i}(j))); + % end + % cstr{i} = v; + + if cstr{j}(1) == "" + cstr{j}(1) = []; + end + if cstr{j}(end) == "" + cstr{j}(end) = []; + end + cstr{j} = arrayfun(@(x) int32(str2double(x)), cstr{j});%SUPER SLOW + +end + + +st.structure_id_path_vec = cstr; + + +anc_id_path = st.structure_id_path_vec{st.id == id}; +if length(anc_id_path) <= depth + anc.anc_id = anc_id_path(end); +else + anc.anc_id = anc_id_path(depth + 1); +end + +anc.anc_name = string(st.name(st.id == anc.anc_id)); +anc.anc_color_hex = string(st.color_hex_triplet(st.id == anc.anc_id,:)); +anc.anc_color_255 = [hex2dec(anc.anc_color_hex{1}(1:2)), hex2dec(anc.anc_color_hex{1}(3:4)), hex2dec(anc.anc_color_hex{1}(5:6))]; +anc.anc_color = anc.anc_color_255/255; +anc.index = st.index(st.id == anc.anc_id,:); + +this.id = id; +this.name = name; +this.color_hex = string(st.color_hex_triplet(st.id == id,:)); +this.color_255 = [hex2dec(this.color_hex{1}(1:2)), hex2dec(this.color_hex{1}(3:4)), hex2dec(this.color_hex{1}(5:6))]; +this.color = this.color_255/255; + + +end + + From d273bb3de86f8cb0dbf3ee991ea97a432a7eacaf Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Mon, 13 Nov 2023 11:55:41 +0000 Subject: [PATCH 32/43] Update updown_probe_with_slider.m bug fixes --- Browsing Functions/updown_probe_with_slider.m | 81 +++++++++---------- 1 file changed, 40 insertions(+), 41 deletions(-) diff --git a/Browsing Functions/updown_probe_with_slider.m b/Browsing Functions/updown_probe_with_slider.m index 04468c8..61db6ae 100644 --- a/Browsing Functions/updown_probe_with_slider.m +++ b/Browsing Functions/updown_probe_with_slider.m @@ -1,5 +1,5 @@ function updown_probe_with_slider(session_id, probeAB, plane, sessions_dir,... - task_dir_name, imaging_session_dir, image_folder_name, depth_level) + task_dir_name, imaging_session_dir, image_folder_name, active_probe_length, depth_level) arguments session_id (1,1) string @@ -9,6 +9,7 @@ function updown_probe_with_slider(session_id, probeAB, plane, sessions_dir,... task_dir_name (1,1) string imaging_session_dir (1,1) string {mustBeFolder} image_folder_name (1,1) string + active_probe_length (1,1) double = 3.84 % in mm depth_level (1,1) double {mustBePositive, mustBeInteger} = 6 end @@ -170,8 +171,6 @@ function updown_probe_with_slider(session_id, probeAB, plane, sessions_dir,... probe_id = T_probes{tf, "probe_id"}; -depth = 6; - Tapdvml_contacts_this = Tapdvml_contacts(Tapdvml_contacts.probe_id == probe_id, :); %% @@ -278,13 +277,13 @@ function move_vertically(uax1, step) ch = uax1.Children; -for i = 1:length(ch) +for j = 1:length(ch) - if isa(ch(i),'matlab.graphics.primitive.Text') - ch(i).Position(2) = ch(i).Position(2) + step*0.01; + if isa(ch(j),'matlab.graphics.primitive.Text') + ch(j).Position(2) = ch(j).Position(2) + step*0.01; - elseif isa(ch(i),'matlab.graphics.primitive.Patch') - ch(i).YData = ch(i).YData + step*0.01; + elseif isa(ch(j),'matlab.graphics.primitive.Patch') + ch(j).YData = ch(j).YData + step*0.01; end @@ -370,8 +369,6 @@ function save_updown_to_table(uax1, ~) m = T_probes{T_probes.session_id == session_id & ... T_probes.probe_AB == probeAB_, ["m_1","m_2","m_3"]}; -%TODO what columns to change or add??? - processed_images_folder = fullfile(image_folder, 'processed'); probePoints = load(fullfile(processed_images_folder, ['probe_points' probe_save_name_suffix])); @@ -421,24 +418,24 @@ function save_updown_to_table(uax1, ~) (Tapdvml_tip.dv_mm - Tapdvml_m.dv_mm)^2 + ... (Tapdvml_tip.ml_mm - Tapdvml_m.ml_mm)^2); -tip2surface_mm_paxinos = sqrt((Tapdvml_tip.ap_mm - Tapdvml_m.ap_mm)^2 + ... - (Tapdvml_tip.dv_mm_paxinos - Tapdvml_m.dv_mm_paxinos)^2 + ... - (Tapdvml_tip.ml_mm - Tapdvml_m.ml_mm)^2); +% tip2surface_mm_paxinos = sqrt((Tapdvml_tip.ap_mm - Tapdvml_m.ap_mm)^2 + ... +% (Tapdvml_tip.dv_mm_paxinos - Tapdvml_m.dv_mm_paxinos)^2 + ... +% (Tapdvml_tip.ml_mm - Tapdvml_m.ml_mm)^2); top_active = (m * active_probe_length + ... - curr_probePoints(tip_index,:) * (tip2surface_mm - active_probe_length))... + new_cp * (tip2surface_mm - active_probe_length))... /tip2surface_mm; % obtained the information for all the 384 channels -a = [linspace(curr_probePoints(tip_index,1), top_active(1), 192)', ... - linspace(curr_probePoints(tip_index,2), top_active(2), 192)', ... - linspace(curr_probePoints(tip_index,3), top_active(3), 192)']; +a = [linspace(new_cp(1), top_active(1), 192)', ... + linspace(new_cp(2), top_active(2), 192)', ... + linspace(new_cp(3), top_active(3), 192)']; probe_contact_points = zeros(384, 3); -for i = 1:192 - probe_contact_points(2*i-1,:) = a(i,:); - probe_contact_points(2*i,:) = a(i,:); +for j = 1:192 + probe_contact_points(2*j-1,:) = a(j,:); + probe_contact_points(2*j,:) = a(j,:); end clear a @@ -446,27 +443,39 @@ function save_updown_to_table(uax1, ~) theta = acos(p(2)); % Assuming p(2) corresponds to the DV direction -for i = 1:384 - c_t_contacts{i} = apdvml2info(probe_contact_points(i,:), av, st, plane); - c_t_contacts{i}.contact_id = repmat(i,height(c_t_contacts{i})); - c_t_contacts{i}.probe_id = repmat(selected_probe,height(c_t_contacts{i})); - c_t_contacts{i}.depth_mm = tip2surface_mm - 0.020 * floor((i-1)/2); +for j = 1:384 + c_t_contacts{j} = apdvml2info(probe_contact_points(j,:), av, st, plane); + c_t_contacts{j}.contact_id = repmat(j,height(c_t_contacts{j})); + c_t_contacts{j}.probe_id = repmat(selected_probe,height(c_t_contacts{j})); + c_t_contacts{j}.depth_mm = tip2surface_mm - 0.020 * floor((j-1)/2); % Project the depth along the line - projected_depth_mm = c_t_contacts{i}.depth_mm * cos(theta); + projected_depth_mm = c_t_contacts{j}.depth_mm * cos(theta); % Convert using the transformation for the chosen plane projected_depth_mm_paxinos = accf2pxs_mm(projected_depth_mm, plane, 'distance'); - c_t_contacts{i}.depth_mm_paxinos = projected_depth_mm_paxinos; + c_t_contacts{j}.depth_mm_paxinos = projected_depth_mm_paxinos; end -c_Tapdvml_contacts{selected_probe} = vertcat(c_t_contacts{:}); +Tapdvml_contacts_new = vertcat(c_t_contacts{:}); clear c_t_contacts -%TODO writetable() +%% update the relevant rows of Tapdvml_contacts + +%% update values +for j = 1:height(Tapdvml_contacts_new) + probe_id = Tapdvml_contacts_new{j,"probe_id"}; + contact_id = Tapdvml_contacts_new{j,"contact_id"}; + + cols = ["ap_mm","dv_mm","dv_mm_paxinos","ml_mm","annotation","name","acronym","contact_id","probe_id","depth_mm","depth_mm_paxinos"]; + Tapdvml_contacts(Tapdvml_contacts.probe_id == probe_id ... + & Tapdvml_contacts.contact_id == contact_id, cols) = ... + Tapdvml_contacts_new(j, cols); +end + +writetable(Tapdvml_contacts, "Tapdvml_contacts.xlsx",'FileType','spreadsheet') -keyboard() end @@ -491,17 +500,7 @@ function save_updown_to_table(uax1, ~) % remove empty string "" -% cemp = cellfun(@(x) x == "" , cstr,'UniformOutput',false); - -for j = 1:size(cstr,1) % i = 741 - % fprintf("j = %d\n",j) - % cstr{i}(cemp{i}) = []; - % v = zeros(1,length(cstr{i}),'int32'); - % for j = 1:length(cstr{i}) - % fprintf("j = %d\n",j) - % v(j) = int32(str2double(cstr{i}(j))); - % end - % cstr{i} = v; +for j = 1:size(cstr,1) if cstr{j}(1) == "" cstr{j}(1) = []; From d52e20ceef54d6faf5963a2af86065c3ce229857 Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Mon, 13 Nov 2023 13:01:42 +0000 Subject: [PATCH 33/43] Update updown_probe_with_slider.m --- Browsing Functions/updown_probe_with_slider.m | 81 +++++++++++++------ 1 file changed, 58 insertions(+), 23 deletions(-) diff --git a/Browsing Functions/updown_probe_with_slider.m b/Browsing Functions/updown_probe_with_slider.m index 61db6ae..4e334dc 100644 --- a/Browsing Functions/updown_probe_with_slider.m +++ b/Browsing Functions/updown_probe_with_slider.m @@ -1,16 +1,39 @@ -function updown_probe_with_slider(session_id, probeAB, plane, sessions_dir,... - task_dir_name, imaging_session_dir, image_folder_name, active_probe_length, depth_level) +function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_dir,... + task_dir_name, imaging_session_dir, image_folder_name, depth_level, ... + active_probe_length, probe_insertion_direction) +% Move up and down a probe using a slider to match the anatomical borders +% to the changes in the firing rates of neurons +% +% EXAMPLE USAGE +% session_id = "kms058-2023-03-25-184034" +% probeAB ="ProbeA"; +% plane = "sagittal"; +% sessions_dir = "\\ettina\Magill_Lab\Julien\Data\head-fixed\by_sessions"; +% task_dir_name = "reaching_go_spout_bar_nov22" +% imaging_session_dir = "\\ettina\Magill_Lab\Kouichi Nakamura\Analysis\Images from Otto\20230406 kms058"; +% image_folder_name = 'RGB_ignore'; +% +% % load the reference brain and region annotations +% if ~exist('av','var') || ~exist('st','var') || ~exist('tv','var') +% disp('loading reference atlas...') +% av = readNPY(annotation_volume_location); +% st = loadStructureTree(structure_tree_location); +% tv = readNPY(template_volume_location); +% end arguments - session_id (1,1) string - probeAB (1,1) string {mustBeMember(probeAB,["ProbeA","ProbeB"])} - plane (1,1) string {mustBeMember(plane,["coronal","sagittal", "transverse"])} - sessions_dir (1,1) string {mustBeFolder} - task_dir_name (1,1) string - imaging_session_dir (1,1) string {mustBeFolder} - image_folder_name (1,1) string - active_probe_length (1,1) double = 3.84 % in mm + av uint16 + st table + session_id (1,1) string % eg. "kms058-2023-03-25-184034" + probeAB (1,1) string {mustBeMember(probeAB,["ProbeA","ProbeB"])} % eg. "ProbeA" + plane (1,1) string {mustBeMember(plane,["coronal","sagittal", "transverse"])} % eg. "sagittal" + sessions_dir (1,1) string {mustBeFolder} % eg. "\\ettina\Magill_Lab\Julien\Data\head-fixed\by_sessions" + task_dir_name (1,1) string % eg. "reaching_go_spout_bar_nov22" + imaging_session_dir (1,1) string {mustBeFolder} % eg. "\\ettina\Magill_Lab\Kouichi Nakamura\Analysis\Images from Otto\20230406 kms058" + image_folder_name (1,1) string = "RGB" depth_level (1,1) double {mustBePositive, mustBeInteger} = 6 + active_probe_length (1,1) double = 3.84 % in mm + probe_insertion_direction (1,1) string {mustBeMember(probe_insertion_direction, ["down","up"])} = "down" end sorter_output_dir = fullfile(sessions_dir, task_dir_name, ... @@ -20,6 +43,17 @@ function updown_probe_with_slider(session_id, probeAB, plane, sessions_dir,... s = load(recording_cell_metrics_path); recording_cell_metrics = s.cell_metrics; +image_folder =fullfile(imaging_session_dir, image_folder_name); %TODO + +probe_save_name_suffix = '_probe'; + +Tapdvml_contacts_path = fullfile(imaging_session_dir, "Tapdvml_contacts.xlsx"); + +% borders_table_path= fullfile(imaging_session_dir, "borders_table.xlsx"); + +Tprobes_path= fullfile(imaging_session_dir, "T_probes.xlsx"); + + opts = detectImportOptions(Tapdvml_contacts_path); opts = setvartype(opts, 'session_id', 'string'); @@ -178,7 +212,7 @@ function updown_probe_with_slider(session_id, probeAB, plane, sessions_dir,... % profile on st.name = regexprep(st.name, '["'']',''); % remove clutters -f = waitbar(0,'1','Name','Approximating pi...',... +f = waitbar(0,'1','Name','Analysing structure hierarchy...',... 'CreateCancelBtn','setappdata(gcbf,''canceling'',1)'); for i = 1:height(Tapdvml_contacts_this) % SLOW @@ -247,14 +281,14 @@ function move_back_to_original(uax1) if isfield(uax1.UserData, 'YLim_shift_mm') ch = uax1.Children; - for i = 1:length(ch) + for j = 1:length(ch) step = - uax1.UserData.YLim_shift_mm; - if isa(ch(i),'matlab.graphics.primitive.Text') - ch(i).Position(2) = ch(i).Position(2) + step; + if isa(ch(j),'matlab.graphics.primitive.Text') + ch(j).Position(2) = ch(j).Position(2) + step; - elseif isa(ch(i),'matlab.graphics.primitive.Patch') - ch(i).YData = ch(i).YData + step; + elseif isa(ch(j),'matlab.graphics.primitive.Patch') + ch(j).YData = ch(j).YData + step; end @@ -345,7 +379,7 @@ function save_updown_to_table(uax1, ~) end %% Job -see plot_and_compute_probe_positions_from_ccf.m +% see plot_and_compute_probe_positions_from_ccf.m % I need the eigen vector, to be saved in T_probes? switch probeAB @@ -374,12 +408,13 @@ function save_updown_to_table(uax1, ~) probePoints = load(fullfile(processed_images_folder, ['probe_points' probe_save_name_suffix])); % get the probe points for the currently analyzed probe + if strcmp(plane,'coronal') - curr_probePoints = probePoints.pointList.pointList{selected_probe,1}(:, [3 2 1]); + curr_probePoints = probePoints.pointList.pointList{probe_id,1}(:, [3 2 1]); elseif strcmp(plane,'sagittal') - curr_probePoints = probePoints.pointList.pointList{selected_probe,1}(:, [1 2 3]); + curr_probePoints = probePoints.pointList.pointList{probe_id,1}(:, [1 2 3]); elseif strcmp(plane,'transverse') - curr_probePoints = probePoints.pointList.pointList{selected_probe,1}(:, [1 3 2]); + curr_probePoints = probePoints.pointList.pointList{probe_id,1}(:, [1 3 2]); end % m is the brain entry point @@ -446,7 +481,7 @@ function save_updown_to_table(uax1, ~) for j = 1:384 c_t_contacts{j} = apdvml2info(probe_contact_points(j,:), av, st, plane); c_t_contacts{j}.contact_id = repmat(j,height(c_t_contacts{j})); - c_t_contacts{j}.probe_id = repmat(selected_probe,height(c_t_contacts{j})); + c_t_contacts{j}.probe_id = repmat(probe_id,height(c_t_contacts{j})); c_t_contacts{j}.depth_mm = tip2surface_mm - 0.020 * floor((j-1)/2); % Project the depth along the line @@ -474,8 +509,8 @@ function save_updown_to_table(uax1, ~) Tapdvml_contacts_new(j, cols); end -writetable(Tapdvml_contacts, "Tapdvml_contacts.xlsx",'FileType','spreadsheet') - +writetable(Tapdvml_contacts, fullfile(imaging_session_dir,"Tapdvml_contacts.xlsx"),'FileType','spreadsheet') +frpint("Tapdvml_contacts.xlsx has been updated.\n") end From beb8d6462cd34ef8a4584285abc97fb20912785d Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Mon, 13 Nov 2023 15:57:03 +0000 Subject: [PATCH 34/43] Update AtlasTransformBrowser.m fprintf to report probe points number --- Browsing Functions/AtlasTransformBrowser.m | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Browsing Functions/AtlasTransformBrowser.m b/Browsing Functions/AtlasTransformBrowser.m index 78cd002..553fb4f 100644 --- a/Browsing Functions/AtlasTransformBrowser.m +++ b/Browsing Functions/AtlasTransformBrowser.m @@ -653,7 +653,8 @@ function hotkeyFcn(f, slice_figure, keydata, allData, save_location, save_suffix disp('transform loaded -- press ''l'' again now to load probe points'); else % load probe points if ~size(ud.pointList{1,1},1) - probe_points = load(fullfile(save_location, ['probe_points' save_suffix])); disp('probe points loaded') + probe_points = load(fullfile(save_location, ['probe_points' save_suffix])); + fprintf('%d probe points loaded\n', length(probe_points.pointList.pointList)) ud.pointList = probe_points.pointList.pointList; ud.pointHands = probe_points.pointList.pointHands; else From 68a06810d94036b6a435b19a2fd6de39bfdea31a Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Mon, 13 Nov 2023 17:10:01 +0000 Subject: [PATCH 35/43] Update updown_probe_with_slider.m --- Browsing Functions/updown_probe_with_slider.m | 68 ++++++++++++++++--- 1 file changed, 60 insertions(+), 8 deletions(-) diff --git a/Browsing Functions/updown_probe_with_slider.m b/Browsing Functions/updown_probe_with_slider.m index 4e334dc..6e986e8 100644 --- a/Browsing Functions/updown_probe_with_slider.m +++ b/Browsing Functions/updown_probe_with_slider.m @@ -36,6 +36,14 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d probe_insertion_direction (1,1) string {mustBeMember(probe_insertion_direction, ["down","up"])} = "down" end +%TODO add button to go next and previous probe +%TODO use probe_id +% related, use animal_id +% make session_id, probeAB + +%TODO what if optic fiber data is chosen? + + sorter_output_dir = fullfile(sessions_dir, task_dir_name, ... session_id, "processed\kilosort", probeAB, "sorter_output"); @@ -49,15 +57,29 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d Tapdvml_contacts_path = fullfile(imaging_session_dir, "Tapdvml_contacts.xlsx"); +original_backup_path = regexprep(Tapdvml_contacts_path, '\.xlsx$' , "_original.xlsx"); + % borders_table_path= fullfile(imaging_session_dir, "borders_table.xlsx"); Tprobes_path= fullfile(imaging_session_dir, "T_probes.xlsx"); +YLim_shift_mm_path = fullfile(imaging_session_dir,"YLim_shift_mm.mat"); opts = detectImportOptions(Tapdvml_contacts_path); opts = setvartype(opts, 'session_id', 'string'); -Tapdvml_contacts = readtable(Tapdvml_contacts_path, opts ); +YLim_shift_mm = []; +if isfile(original_backup_path) + Tapdvml_contacts = readtable(original_backup_path, opts ); + + if isfile(YLim_shift_mm_path) + S = load(YLim_shift_mm_path); % YLim_shift_mm + YLim_shift_mm = S.YLim_shift_mm; + clear S + end +else + Tapdvml_contacts = readtable(Tapdvml_contacts_path, opts ); +end % Tborders_table = readtable(borders_table_path); @@ -80,6 +102,11 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d ylabel(uax1, 'Dsitance from the tip (mm)') uax1.XTick =[]; +if ~isempty(YLim_shift_mm) + assert(isequal(size(YLim_shift_mm), [1 1])) + uax1.UserData.YLim_shift_mm = YLim_shift_mm; +end + uax2 = uiaxes(ufig); uax2.PositionConstraint = 'innerposition'; uax2.Position=[448, 150, 80, 800]; @@ -209,7 +236,7 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d %% anc = preallocatestruct(["index", "anc_id", "anc_name", "anc_color_hex", "anc_color_255", "anc_color"], [height(Tapdvml_contacts_this), 1]); -% profile on + st.name = regexprep(st.name, '["'']',''); % remove clutters f = waitbar(0,'1','Name','Analysing structure hierarchy...',... @@ -217,16 +244,17 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d for i = 1:height(Tapdvml_contacts_this) % SLOW waitbar( i/height(Tapdvml_contacts_this), f, ... - fprintf("%d of %d", i, height(Tapdvml_contacts_this))) + sprintf("%d of %d", i, height(Tapdvml_contacts_this))) %TODO not shown properly name = Tapdvml_contacts_this.name{i}; [anc(i)] = find_name_for_depth(st, depth_level, name); %TODO stopped working for i = 61, i = 1327 end -% profile viewer delete(f) %% Tanc = struct2table(anc, AsArray=true); +%TODO save a cache file?? + Tapdvml_contacts_this = [Tapdvml_contacts_this, Tanc]; @@ -274,6 +302,22 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d uax1.YLabel.FontSize = 14; +% move items after initialiation +if isfield(uax1.UserData, 'YLim_shift_mm') + uax1.UserData.YLim_shift_mm + + ch = uax1.Children; + + for i = 1:length(ch) + if isa(ch(i),'matlab.graphics.primitive.Text') + ch(i).Position(2) = ch(i).Position(2) + YLim_shift_mm; + + elseif isa(ch(i),'matlab.graphics.primitive.Patch') + ch(i).YData = ch(i).YData + YLim_shift_mm; + + end + end +end function move_back_to_original(uax1) @@ -363,10 +407,8 @@ function save_updown_to_table(uax1, ~) new_file_path = regexprep(Tapdvml_contacts_path, '\.xlsx$' , "_" + suffix + ".xlsx"); -original_backup_path = regexprep(Tapdvml_contacts_path, '\.xlsx$' , "_" + suffix + "_original.xlsx"); - if isfile(original_backup_path) - % save a backup file + % save a backup file for the latest copyfile(Tapdvml_contacts_path, new_file_path) @@ -439,6 +481,14 @@ function save_updown_to_table(uax1, ~) YLim_shift_mm = uax1.UserData.YLim_shift_mm; new_cp = curr_probePoints(tip_index,:) - YLim_shift_mm*100; % minus to move up +%TODO in case the probe has already moved prevously, YLim_shift_mm won't be +% accurate +% +% always read both the original and the latest Excel +% YLim_shift_mm should represent the gap from the original + +%GOOD always read the original Excel and saved value of YLim_shift_mm +% YLim_shift_mm should represent the gap from the original Tapdvml_m = apdvml2info(m, av, st, plane); @@ -510,7 +560,9 @@ function save_updown_to_table(uax1, ~) end writetable(Tapdvml_contacts, fullfile(imaging_session_dir,"Tapdvml_contacts.xlsx"),'FileType','spreadsheet') -frpint("Tapdvml_contacts.xlsx has been updated.\n") +fprintf("Tapdvml_contacts.xlsx has been updated with %.3f mm shift from the original.\n", YLim_shift_mm) +% save the value of YLim_shift_mm +save(YLim_shift_mm_path, 'YLim_shift_mm'); end From 226124e4717ba7fe1b945bbb1f847775f820f570 Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Thu, 16 Nov 2023 10:01:03 +0000 Subject: [PATCH 36/43] Update updown_probe_with_slider.m --- Browsing Functions/updown_probe_with_slider.m | 43 +++++++------------ 1 file changed, 16 insertions(+), 27 deletions(-) diff --git a/Browsing Functions/updown_probe_with_slider.m b/Browsing Functions/updown_probe_with_slider.m index 6e986e8..29564a5 100644 --- a/Browsing Functions/updown_probe_with_slider.m +++ b/Browsing Functions/updown_probe_with_slider.m @@ -37,9 +37,6 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d end %TODO add button to go next and previous probe -%TODO use probe_id -% related, use animal_id -% make session_id, probeAB %TODO what if optic fiber data is chosen? @@ -51,7 +48,7 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d s = load(recording_cell_metrics_path); recording_cell_metrics = s.cell_metrics; -image_folder =fullfile(imaging_session_dir, image_folder_name); %TODO +image_folder =fullfile(imaging_session_dir, image_folder_name); probe_save_name_suffix = '_probe'; @@ -149,9 +146,6 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d bt7 = uibutton(ufig, "Text","Save", "ButtonPushedFcn", @(src, event) save_updown_to_table(uax1)); bt7.Position = [150, 30, 300, 50]; -tx1 = annotation(ufig, 'textbox', [0.2 0.85 0.5 0.1], 'String', session_id +"/" + probeAB,... - HorizontalAlignment='center',FontSize=18,LineStyle='none'); - %% Plot cells n = length(recording_cell_metrics.firingRate); @@ -234,6 +228,10 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d Tapdvml_contacts_this = Tapdvml_contacts(Tapdvml_contacts.probe_id == probe_id, :); +tx1 = annotation(ufig, 'textbox', [0.2 0.85 0.8 0.1], 'String', ... + sprintf("%s/%s (%d)", session_id, probeAB, probe_id),... + HorizontalAlignment='center',FontSize=18,LineStyle='none'); + %% anc = preallocatestruct(["index", "anc_id", "anc_name", "anc_color_hex", "anc_color_255", "anc_color"], [height(Tapdvml_contacts_this), 1]); @@ -244,17 +242,15 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d for i = 1:height(Tapdvml_contacts_this) % SLOW waitbar( i/height(Tapdvml_contacts_this), f, ... - sprintf("%d of %d", i, height(Tapdvml_contacts_this))) %TODO not shown properly + sprintf("%d of %d", i, height(Tapdvml_contacts_this))); name = Tapdvml_contacts_this.name{i}; - [anc(i)] = find_name_for_depth(st, depth_level, name); %TODO stopped working for i = 61, i = 1327 + [anc(i)] = find_name_for_depth(st, depth_level, name); end delete(f) %% Tanc = struct2table(anc, AsArray=true); -%TODO save a cache file?? - Tapdvml_contacts_this = [Tapdvml_contacts_this, Tanc]; @@ -384,7 +380,7 @@ function save_updown_to_table(uax1, ~) return else if isfield(uax1.UserData, 'YLim_shift_mm') - %TODO move + % move else disp('YLim_shift_mm is not found. No change made.') return @@ -392,7 +388,7 @@ function save_updown_to_table(uax1, ~) end -%TODO save the backup +%% save the backup assert(endsWith(Tapdvml_contacts_path,'.xlsx')) @@ -461,7 +457,7 @@ function save_updown_to_table(uax1, ~) % m is the brain entry point -%TODO modifiy the m and curr_probePoints(tip_index,:) by YLim_shift_mm +% modifiy the m and curr_probePoints(tip_index,:) by YLim_shift_mm % p is in 10 µm if strcmp(probe_insertion_direction, 'down') @@ -481,23 +477,16 @@ function save_updown_to_table(uax1, ~) YLim_shift_mm = uax1.UserData.YLim_shift_mm; new_cp = curr_probePoints(tip_index,:) - YLim_shift_mm*100; % minus to move up -%TODO in case the probe has already moved prevously, YLim_shift_mm won't be -% accurate -% -% always read both the original and the latest Excel -% YLim_shift_mm should represent the gap from the original - -%GOOD always read the original Excel and saved value of YLim_shift_mm +%NOTE always read the original Excel and saved value of YLim_shift_mm % YLim_shift_mm should represent the gap from the original - Tapdvml_m = apdvml2info(m, av, st, plane); Tapdvml_tip = apdvml2info(new_cp, av, st, plane); -% measure the distance +%% measure the distance tip2surface_mm = sqrt((Tapdvml_tip.ap_mm - Tapdvml_m.ap_mm)^2 + ... (Tapdvml_tip.dv_mm - Tapdvml_m.dv_mm)^2 + ... @@ -511,8 +500,7 @@ function save_updown_to_table(uax1, ~) new_cp * (tip2surface_mm - active_probe_length))... /tip2surface_mm; -% obtained the information for all the 384 channels - +%% obtain the information for all the 384 channels a = [linspace(new_cp(1), top_active(1), 192)', ... linspace(new_cp(2), top_active(2), 192)', ... @@ -553,14 +541,15 @@ function save_updown_to_table(uax1, ~) probe_id = Tapdvml_contacts_new{j,"probe_id"}; contact_id = Tapdvml_contacts_new{j,"contact_id"}; - cols = ["ap_mm","dv_mm","dv_mm_paxinos","ml_mm","annotation","name","acronym","contact_id","probe_id","depth_mm","depth_mm_paxinos"]; + cols = ["ap_mm","dv_mm","dv_mm_paxinos","ml_mm","annotation",... + "name","acronym","contact_id","probe_id","depth_mm","depth_mm_paxinos"]; Tapdvml_contacts(Tapdvml_contacts.probe_id == probe_id ... & Tapdvml_contacts.contact_id == contact_id, cols) = ... Tapdvml_contacts_new(j, cols); end writetable(Tapdvml_contacts, fullfile(imaging_session_dir,"Tapdvml_contacts.xlsx"),'FileType','spreadsheet') -fprintf("Tapdvml_contacts.xlsx has been updated with %.3f mm shift from the original.\n", YLim_shift_mm) +fprintf("Tapdvml_contacts.xlsx has been updated with %.3f mm shift from the original for %s of %s (%d).\n", YLim_shift_mm, probeAB, session_id, probe_id) % save the value of YLim_shift_mm save(YLim_shift_mm_path, 'YLim_shift_mm'); From 277ba723b86a68ac966e9d77572c4da5f3e8e92e Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Thu, 16 Nov 2023 16:44:43 +0000 Subject: [PATCH 37/43] Update updown_probe_with_slider.m --- Browsing Functions/updown_probe_with_slider.m | 512 ++++++++++-------- 1 file changed, 286 insertions(+), 226 deletions(-) diff --git a/Browsing Functions/updown_probe_with_slider.m b/Browsing Functions/updown_probe_with_slider.m index 29564a5..99665fa 100644 --- a/Browsing Functions/updown_probe_with_slider.m +++ b/Browsing Functions/updown_probe_with_slider.m @@ -5,6 +5,7 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d % to the changes in the firing rates of neurons % % EXAMPLE USAGE +% % session_id = "kms058-2023-03-25-184034" % probeAB ="ProbeA"; % plane = "sagittal"; @@ -41,6 +42,8 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d %TODO what if optic fiber data is chosen? + + sorter_output_dir = fullfile(sessions_dir, task_dir_name, ... session_id, "processed\kilosort", probeAB, "sorter_output"); @@ -60,32 +63,90 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d Tprobes_path= fullfile(imaging_session_dir, "T_probes.xlsx"); -YLim_shift_mm_path = fullfile(imaging_session_dir,"YLim_shift_mm.mat"); +%YLim_shift_mm_path = fullfile(imaging_session_dir,"YLim_shift_mm.mat"); %TODO +Y_shift_path = fullfile(imaging_session_dir,"T_Y_shift.mat"); %TODO + + + + + +% Tborders_table = readtable(borders_table_path); + +T_probes= readtable(Tprobes_path); + +if ~ ismember('p_1',T_probes.Properties.VariableNames) || ... + ~ ismember('p_2',T_probes.Properties.VariableNames) ||... + ~ ismember('p_3',T_probes.Properties.VariableNames) + warning('The eigen vectors p is not found in T_probes. You cannot update the Tapdvml_contacts.') +end + +switch probeAB + case "ProbeA" + probeAB_short = "A"; + case "ProbeB" + probeAB_short = "B"; + otherwise + +end +tf = T_probes.session_id == session_id & T_probes.probe_AB == probeAB_short; + +probe_id = T_probes{tf, "probe_id"}; + +%% deal with the cache file opts = detectImportOptions(Tapdvml_contacts_path); opts = setvartype(opts, 'session_id', 'string'); YLim_shift_mm = []; if isfile(original_backup_path) + Tapdvml_contacts = readtable(original_backup_path, opts ); - - if isfile(YLim_shift_mm_path) - S = load(YLim_shift_mm_path); % YLim_shift_mm - YLim_shift_mm = S.YLim_shift_mm; + Tapdvml_contacts_this = Tapdvml_contacts(Tapdvml_contacts.probe_id == probe_id, :); + + if isfile(Y_shift_path) % use the cache + S = load(Y_shift_path); % YLim_shift_mm + + assert(istable(S.T_Y_shift)) + T_Y_shift = S.T_Y_shift; + YLim_shift_mm = S.T_Y_shift.YLim_shift_mm(probe_id); + Tanc = S.T_Y_shift.Tanc{probe_id}; clear S + else + % cache file is not found + Tanc_ = cell(height(T_probes), 1); + for i = 1:size(Tanc_,1) + Tanc_{i} = table(); + end + + T_Y_shift = table(zeros(height(T_probes), 1), Tanc_,... + 'VariableNames', ["YLim_shift_mm","Tanc"]); % YLim_shift_mm and Tanc + clear Tanc_ + Tanc = compute_Tanc(Tapdvml_contacts_this, st, depth_level); + T_Y_shift.Tanc{probe_id} = Tanc; + save(Y_shift_path, 'T_Y_shift'); end + Tapdvml_contacts_this = [Tapdvml_contacts_this, Tanc]; + else + % the very first run + Tapdvml_contacts = readtable(Tapdvml_contacts_path, opts ); -end + Tapdvml_contacts_this = Tapdvml_contacts(Tapdvml_contacts.probe_id == probe_id, :); -% Tborders_table = readtable(borders_table_path); + Tanc_ = cell(height(T_probes), 1); + for i = 1:size(Tanc_,1) + Tanc_{i} = table(); + end -T_probes= readtable(Tprobes_path); + T_Y_shift = table(zeros(height(T_probes), 1), Tanc_,... + 'VariableNames', ["YLim_shift_mm","Tanc"]); % YLim_shift_mm and Tanc + clear Tanc_ + Tanc = compute_Tanc(Tapdvml_contacts_this, st, depth_level); + T_Y_shift.Tanc{probe_id} = Tanc; + save(Y_shift_path, 'T_Y_shift'); + + Tapdvml_contacts_this = [Tapdvml_contacts_this, Tanc]; -if ~ ismember('p_1',T_probes.Properties.VariableNames) || ... - ~ ismember('p_2',T_probes.Properties.VariableNames) ||... - ~ ismember('p_3',T_probes.Properties.VariableNames) - warning('The eigen vectors p is not found in T_probes. You cannot update the Tapdvml_contacts.') end @@ -195,64 +256,14 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d uax3.CLim = uax2.CLim; % same color scale -%% Read the Allen CCF data -% directory of reference atlas files -annotation_volume_location = "\\ettina\Magill_lab\Kouichi Nakamura\Analysis\allenCCFdata\annotation_volume_10um_by_index.npy"; -structure_tree_location = "\\ettina\Magill_lab\Kouichi Nakamura\Analysis\allenCCFdata\structure_tree_safe_2017.csv"; -template_volume_location = "\\ettina\Magill_lab\Kouichi Nakamura\Analysis\allenCCFdata\template_volume_10um.npy"; - -% load the reference brain and region annotations -if ~exist('av','var') || ~exist('st','var') || ~exist('tv','var') - disp('loading reference atlas...') - av = readNPY(annotation_volume_location); - st = loadStructureTree(structure_tree_location); - tv = readNPY(template_volume_location); -end - - -probeAB = "ProbeA"; - -switch probeAB - case "ProbeA" - probeAB_short = "A"; - case "ProbeB" - probeAB_short = "B"; - otherwise - -end - -tf = T_probes.session_id == session_id & T_probes.probe_AB == probeAB_short; - -probe_id = T_probes{tf, "probe_id"}; - -Tapdvml_contacts_this = Tapdvml_contacts(Tapdvml_contacts.probe_id == probe_id, :); +% Tapdvml_contacts_this = Tapdvml_contacts(Tapdvml_contacts.probe_id == probe_id, :); tx1 = annotation(ufig, 'textbox', [0.2 0.85 0.8 0.1], 'String', ... sprintf("%s/%s (%d)", session_id, probeAB, probe_id),... HorizontalAlignment='center',FontSize=18,LineStyle='none'); -%% -anc = preallocatestruct(["index", "anc_id", "anc_name", "anc_color_hex", "anc_color_255", "anc_color"], [height(Tapdvml_contacts_this), 1]); - -st.name = regexprep(st.name, '["'']',''); % remove clutters - -f = waitbar(0,'1','Name','Analysing structure hierarchy...',... - 'CreateCancelBtn','setappdata(gcbf,''canceling'',1)'); - -for i = 1:height(Tapdvml_contacts_this) % SLOW - waitbar( i/height(Tapdvml_contacts_this), f, ... - sprintf("%d of %d", i, height(Tapdvml_contacts_this))); - name = Tapdvml_contacts_this.name{i}; - [anc(i)] = find_name_for_depth(st, depth_level, name); -end - -delete(f) - -%% -Tanc = struct2table(anc, AsArray=true); - -Tapdvml_contacts_this = [Tapdvml_contacts_this, Tanc]; +% Tapdvml_contacts_this = [Tapdvml_contacts_this, Tanc]; %TODO % Sample data @@ -279,6 +290,7 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d %% Plot structures + cla(uax1) uax1.Color = 'k'; @@ -294,13 +306,25 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d 'HorizontalAlignment','center', Clipping='on'); end +horzline = line(uax1, [0 1], [0 0], 'Color','r'); % distance from the tip 0 mm + ylim(uax1, [-0.100, 3.940]) uax1.YLabel.FontSize = 14; +if isempty(uax1.UserData) || ~isfield(uax1.UserData, 'YLim_shift_mm') + tx2 = annotation(ufig, 'textbox', [0.2 0.03 0.3 0.1], 'String', ... + regexprep(sprintf("%+.2f mm", 0), '-', char(8722)),... + HorizontalAlignment='left',FontSize=12, LineStyle='none'); +else + tx2 = annotation(ufig, 'textbox', [0.2 0.03 0.3 0.1], 'String', ... + regexprep(sprintf("%+.2f mm", uax1.UserData.YLim_shift_mm), '-', char(8722)),... + HorizontalAlignment='left',FontSize=12, LineStyle='none'); +end + % move items after initialiation if isfield(uax1.UserData, 'YLim_shift_mm') - uax1.UserData.YLim_shift_mm + % uax1.UserData.YLim_shift_mm ch = uax1.Children; @@ -311,249 +335,263 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d elseif isa(ch(i),'matlab.graphics.primitive.Patch') ch(i).YData = ch(i).YData + YLim_shift_mm; + elseif isa(ch(i),'matlab.graphics.primitive.Line') + ch(i).YData = ch(i).YData + YLim_shift_mm; + end end end -function move_back_to_original(uax1) + function move_back_to_original(uax1) -if isfield(uax1.UserData, 'YLim_shift_mm') - ch = uax1.Children; + if isfield(uax1.UserData, 'YLim_shift_mm') + ch = uax1.Children; + + for j = 1:length(ch) + step = - uax1.UserData.YLim_shift_mm; - for j = 1:length(ch) - step = - uax1.UserData.YLim_shift_mm; + if isa(ch(j),'matlab.graphics.primitive.Text') + ch(j).Position(2) = ch(j).Position(2) + step; - if isa(ch(j),'matlab.graphics.primitive.Text') - ch(j).Position(2) = ch(j).Position(2) + step; + elseif isa(ch(j),'matlab.graphics.primitive.Patch') + ch(j).YData = ch(j).YData + step; - elseif isa(ch(j),'matlab.graphics.primitive.Patch') - ch(j).YData = ch(j).YData + step; + elseif isa(ch(j),'matlab.graphics.primitive.Line') + ch(j).YData = ch(j).YData + step; + + end + end + + uax1.UserData.YLim_shift_mm = 0; + + tx2.String = regexprep(sprintf("%+.2f mm", uax1.UserData.YLim_shift_mm), '-', char(8722)); + else + %nothing to do end end - uax1.UserData.YLim_shift_mm = 0; -else - %nothing to do -end -end + function move_vertically(uax1, step) + % arguments + % uax1(1,1) + % step (1,1) {mustBeInteger} + % end -function move_vertically(uax1, step) + ch = uax1.Children; -% arguments -% uax1(1,1) -% step (1,1) {mustBeInteger} -% end + for j = 1:length(ch) -ch = uax1.Children; + if isa(ch(j),'matlab.graphics.primitive.Text') + ch(j).Position(2) = ch(j).Position(2) + step*0.01; -for j = 1:length(ch) + elseif isa(ch(j),'matlab.graphics.primitive.Patch') + ch(j).YData = ch(j).YData + step*0.01; - if isa(ch(j),'matlab.graphics.primitive.Text') - ch(j).Position(2) = ch(j).Position(2) + step*0.01; + elseif isa(ch(j),'matlab.graphics.primitive.Line') + ch(j).YData = ch(j).YData + step*0.01; - elseif isa(ch(j),'matlab.graphics.primitive.Patch') - ch(j).YData = ch(j).YData + step*0.01; + end - end + end -end + if isfield(uax1.UserData, 'YLim_shift_mm') + uax1.UserData.YLim_shift_mm = uax1.UserData.YLim_shift_mm + step*0.01; + else + uax1.UserData.YLim_shift_mm = step*0.01; + end -if isfield(uax1.UserData, 'YLim_shift_mm') - uax1.UserData.YLim_shift_mm = uax1.UserData.YLim_shift_mm + step*0.01; -else - uax1.UserData.YLim_shift_mm = step*0.01; -end + tx2.String = regexprep(sprintf("%+.2f mm", uax1.UserData.YLim_shift_mm), '-', char(8722)); -end + end + function save_updown_to_table(uax1, ~) -function save_updown_to_table(uax1, ~) + if isempty(uax1.UserData) + disp('The probe has not moved on GUI yet. No change made.') + return + else + if isfield(uax1.UserData, 'YLim_shift_mm') + % move + else + disp('YLim_shift_mm is not found. No change made.') + return + end + end -if isempty(uax1.UserData) - disp('The probe has not moved on GUI yet. No change made.') - return -else - if isfield(uax1.UserData, 'YLim_shift_mm') - % move - else - disp('YLim_shift_mm is not found. No change made.') - return - end -end + %% save the backup -%% save the backup + assert(endsWith(Tapdvml_contacts_path,'.xlsx')) -assert(endsWith(Tapdvml_contacts_path,'.xlsx')) + assert(isfile(Tapdvml_contacts_path)) -assert(isfile(Tapdvml_contacts_path)) + f_info = dir(Tapdvml_contacts_path); -f_info = dir(Tapdvml_contacts_path); + dt = datetime(f_info.date, 'InputFormat', 'dd-MMM-yyyy HH:mm:ss'); -dt = datetime(f_info.date, 'InputFormat', 'dd-MMM-yyyy HH:mm:ss'); + dt.Format = 'yyyy-MM-dd_HH-mm-ss'; + suffix = string(dt); -dt.Format = 'yyyy-MM-dd_HH-mm-ss'; -suffix = string(dt); + new_file_path = regexprep(Tapdvml_contacts_path, '\.xlsx$' , "_" + suffix + ".xlsx"); -new_file_path = regexprep(Tapdvml_contacts_path, '\.xlsx$' , "_" + suffix + ".xlsx"); + if isfile(original_backup_path) + % save a backup file for the latest -if isfile(original_backup_path) - % save a backup file for the latest + copyfile(Tapdvml_contacts_path, new_file_path) - copyfile(Tapdvml_contacts_path, new_file_path) + else + % this is the first time you change the Excel file + % save a backup file -else - % this is the first time you change the Excel file - % save a backup file + copyfile(Tapdvml_contacts_path, original_backup_path) - copyfile(Tapdvml_contacts_path, original_backup_path) + end -end + %% Job + % see plot_and_compute_probe_positions_from_ccf.m + % I need the eigen vector, to be saved in T_probes? + + switch probeAB + case "ProbeA" + probeAB_ = "A"; + case "ProbeB" + probeAB_ = "B"; + case "optic fiber" + probeAB_ = "optic fiber"; + end -%% Job -% see plot_and_compute_probe_positions_from_ccf.m -% I need the eigen vector, to be saved in T_probes? + % eigen vector + tf = T_probes.session_id == session_id & ... + T_probes.probe_AB == probeAB_; -switch probeAB - case "ProbeA" - probeAB_ = "A"; - case "ProbeB" - probeAB_ = "B"; - case "optic fiber" - probeAB_ = "optic fiber"; -end + assert(nnz(tf) == 1) -% eigen vector -tf = T_probes.session_id == session_id & ... - T_probes.probe_AB == probeAB_; + p = T_probes{T_probes.session_id == session_id & ... + T_probes.probe_AB == probeAB_, ["p_1","p_2","p_3"]}; -assert(nnz(tf) == 1) + m = T_probes{T_probes.session_id == session_id & ... + T_probes.probe_AB == probeAB_, ["m_1","m_2","m_3"]}; -p = T_probes{T_probes.session_id == session_id & ... - T_probes.probe_AB == probeAB_, ["p_1","p_2","p_3"]}; + processed_images_folder = fullfile(image_folder, 'processed'); -m = T_probes{T_probes.session_id == session_id & ... - T_probes.probe_AB == probeAB_, ["m_1","m_2","m_3"]}; + probePoints = load(fullfile(processed_images_folder, ['probe_points' probe_save_name_suffix])); -processed_images_folder = fullfile(image_folder, 'processed'); + % get the probe points for the currently analyzed probe -probePoints = load(fullfile(processed_images_folder, ['probe_points' probe_save_name_suffix])); + if strcmp(plane,'coronal') + curr_probePoints = probePoints.pointList.pointList{probe_id,1}(:, [3 2 1]); + elseif strcmp(plane,'sagittal') + curr_probePoints = probePoints.pointList.pointList{probe_id,1}(:, [1 2 3]); + elseif strcmp(plane,'transverse') + curr_probePoints = probePoints.pointList.pointList{probe_id,1}(:, [1 3 2]); + end -% get the probe points for the currently analyzed probe + % m is the brain entry point -if strcmp(plane,'coronal') - curr_probePoints = probePoints.pointList.pointList{probe_id,1}(:, [3 2 1]); -elseif strcmp(plane,'sagittal') - curr_probePoints = probePoints.pointList.pointList{probe_id,1}(:, [1 2 3]); -elseif strcmp(plane,'transverse') - curr_probePoints = probePoints.pointList.pointList{probe_id,1}(:, [1 3 2]); -end + % modifiy the m and curr_probePoints(tip_index,:) by YLim_shift_mm + % p is in 10 µm -% m is the brain entry point + if strcmp(probe_insertion_direction, 'down') + [depth, tip_index] = max(curr_probePoints(:,2)); + elseif strcmp(probe_insertion_direction, 'up') + [depth, tip_index] = min(curr_probePoints(:,2)); + end -% modifiy the m and curr_probePoints(tip_index,:) by YLim_shift_mm -% p is in 10 µm + % curr_probePoints(tip_index,:) either to get it from the ProbePoints or + % from Tapdvml_contacts.xlsx (convert ap_mm, dv_mm, ml_mm back to allen coordinates) + % opposite of apdvml2info -if strcmp(probe_insertion_direction, 'down') - [depth, tip_index] = max(curr_probePoints(:,2)); -elseif strcmp(probe_insertion_direction, 'up') - [depth, tip_index] = min(curr_probePoints(:,2)); -end -% curr_probePoints(tip_index,:) either to get it from the ProbePoints or -% from Tapdvml_contacts.xlsx (convert ap_mm, dv_mm, ml_mm back to allen coordinates) -% opposite of apdvml2info + % Tapdvml_contacts.session_id + % tip_mm = -% Tapdvml_contacts.session_id + YLim_shift_mm = uax1.UserData.YLim_shift_mm; + new_cp = curr_probePoints(tip_index,:) - YLim_shift_mm*100; % minus to move up + %NOTE always read the original Excel and saved value of YLim_shift_mm + % YLim_shift_mm should represent the gap from the original -% tip_mm = + Tapdvml_m = apdvml2info(m, av, st, plane); -YLim_shift_mm = uax1.UserData.YLim_shift_mm; -new_cp = curr_probePoints(tip_index,:) - YLim_shift_mm*100; % minus to move up -%NOTE always read the original Excel and saved value of YLim_shift_mm -% YLim_shift_mm should represent the gap from the original + Tapdvml_tip = apdvml2info(new_cp, av, st, plane); -Tapdvml_m = apdvml2info(m, av, st, plane); -Tapdvml_tip = apdvml2info(new_cp, av, st, plane); + %% measure the distance + tip2surface_mm = sqrt((Tapdvml_tip.ap_mm - Tapdvml_m.ap_mm)^2 + ... + (Tapdvml_tip.dv_mm - Tapdvml_m.dv_mm)^2 + ... + (Tapdvml_tip.ml_mm - Tapdvml_m.ml_mm)^2); -%% measure the distance + % tip2surface_mm_paxinos = sqrt((Tapdvml_tip.ap_mm - Tapdvml_m.ap_mm)^2 + ... + % (Tapdvml_tip.dv_mm_paxinos - Tapdvml_m.dv_mm_paxinos)^2 + ... + % (Tapdvml_tip.ml_mm - Tapdvml_m.ml_mm)^2); -tip2surface_mm = sqrt((Tapdvml_tip.ap_mm - Tapdvml_m.ap_mm)^2 + ... - (Tapdvml_tip.dv_mm - Tapdvml_m.dv_mm)^2 + ... - (Tapdvml_tip.ml_mm - Tapdvml_m.ml_mm)^2); + top_active = (m * active_probe_length + ... + new_cp * (tip2surface_mm - active_probe_length))... + /tip2surface_mm; -% tip2surface_mm_paxinos = sqrt((Tapdvml_tip.ap_mm - Tapdvml_m.ap_mm)^2 + ... -% (Tapdvml_tip.dv_mm_paxinos - Tapdvml_m.dv_mm_paxinos)^2 + ... -% (Tapdvml_tip.ml_mm - Tapdvml_m.ml_mm)^2); + %% obtain the information for all the 384 channels -top_active = (m * active_probe_length + ... - new_cp * (tip2surface_mm - active_probe_length))... - /tip2surface_mm; + a = [linspace(new_cp(1), top_active(1), 192)', ... + linspace(new_cp(2), top_active(2), 192)', ... + linspace(new_cp(3), top_active(3), 192)']; + probe_contact_points = zeros(384, 3); + for j = 1:192 + probe_contact_points(2*j-1,:) = a(j,:); + probe_contact_points(2*j,:) = a(j,:); + end + clear a -%% obtain the information for all the 384 channels + c_t_contacts = cell(384,1); -a = [linspace(new_cp(1), top_active(1), 192)', ... - linspace(new_cp(2), top_active(2), 192)', ... - linspace(new_cp(3), top_active(3), 192)']; -probe_contact_points = zeros(384, 3); -for j = 1:192 - probe_contact_points(2*j-1,:) = a(j,:); - probe_contact_points(2*j,:) = a(j,:); -end -clear a + theta = acos(p(2)); % Assuming p(2) corresponds to the DV direction -c_t_contacts = cell(384,1); + for j = 1:384 + c_t_contacts{j} = apdvml2info(probe_contact_points(j,:), av, st, plane); + c_t_contacts{j}.contact_id = repmat(j,height(c_t_contacts{j})); + c_t_contacts{j}.probe_id = repmat(probe_id,height(c_t_contacts{j})); + c_t_contacts{j}.depth_mm = tip2surface_mm - 0.020 * floor((j-1)/2); -theta = acos(p(2)); % Assuming p(2) corresponds to the DV direction + % Project the depth along the line + projected_depth_mm = c_t_contacts{j}.depth_mm * cos(theta); -for j = 1:384 - c_t_contacts{j} = apdvml2info(probe_contact_points(j,:), av, st, plane); - c_t_contacts{j}.contact_id = repmat(j,height(c_t_contacts{j})); - c_t_contacts{j}.probe_id = repmat(probe_id,height(c_t_contacts{j})); - c_t_contacts{j}.depth_mm = tip2surface_mm - 0.020 * floor((j-1)/2); - - % Project the depth along the line - projected_depth_mm = c_t_contacts{j}.depth_mm * cos(theta); + % Convert using the transformation for the chosen plane + projected_depth_mm_paxinos = accf2pxs_mm(projected_depth_mm, plane, 'distance'); - % Convert using the transformation for the chosen plane - projected_depth_mm_paxinos = accf2pxs_mm(projected_depth_mm, plane, 'distance'); + c_t_contacts{j}.depth_mm_paxinos = projected_depth_mm_paxinos; + end - c_t_contacts{j}.depth_mm_paxinos = projected_depth_mm_paxinos; -end + Tapdvml_contacts_new = vertcat(c_t_contacts{:}); + clear c_t_contacts -Tapdvml_contacts_new = vertcat(c_t_contacts{:}); -clear c_t_contacts + %% update the relevant rows of Tapdvml_contacts -%% update the relevant rows of Tapdvml_contacts + %% update values + for j = 1:height(Tapdvml_contacts_new) + probe_id = Tapdvml_contacts_new{j,"probe_id"}; + contact_id = Tapdvml_contacts_new{j,"contact_id"}; -%% update values -for j = 1:height(Tapdvml_contacts_new) - probe_id = Tapdvml_contacts_new{j,"probe_id"}; - contact_id = Tapdvml_contacts_new{j,"contact_id"}; + cols = ["ap_mm","dv_mm","dv_mm_paxinos","ml_mm","annotation",... + "name","acronym","contact_id","probe_id","depth_mm","depth_mm_paxinos"]; + Tapdvml_contacts(Tapdvml_contacts.probe_id == probe_id ... + & Tapdvml_contacts.contact_id == contact_id, cols) = ... + Tapdvml_contacts_new(j, cols); + end - cols = ["ap_mm","dv_mm","dv_mm_paxinos","ml_mm","annotation",... - "name","acronym","contact_id","probe_id","depth_mm","depth_mm_paxinos"]; - Tapdvml_contacts(Tapdvml_contacts.probe_id == probe_id ... - & Tapdvml_contacts.contact_id == contact_id, cols) = ... - Tapdvml_contacts_new(j, cols); -end + writetable(Tapdvml_contacts, fullfile(imaging_session_dir,"Tapdvml_contacts.xlsx"),'FileType','spreadsheet') + fprintf("Tapdvml_contacts.xlsx has been updated with %.3f mm shift from the original for %s of %s (%d).\n", YLim_shift_mm, probeAB, session_id, probe_id) + % save the value of YLim_shift_mm + T_Y_shift.YLim_shift_mm(probe_id) = YLim_shift_mm; + save(Y_shift_path, 'T_Y_shift'); -writetable(Tapdvml_contacts, fullfile(imaging_session_dir,"Tapdvml_contacts.xlsx"),'FileType','spreadsheet') -fprintf("Tapdvml_contacts.xlsx has been updated with %.3f mm shift from the original for %s of %s (%d).\n", YLim_shift_mm, probeAB, session_id, probe_id) -% save the value of YLim_shift_mm -save(YLim_shift_mm_path, 'YLim_shift_mm'); -end + end end @@ -615,3 +653,25 @@ function save_updown_to_table(uax1, ~) end +function Tanc = compute_Tanc(Tapdvml_contacts_this, st, depth_level) + +anc = preallocatestruct(["index", "anc_id", "anc_name", "anc_color_hex", "anc_color_255", "anc_color"], [height(Tapdvml_contacts_this), 1]); + +st.name = regexprep(st.name, '["'']',''); % remove clutters + +f = waitbar(0,'1','Name','Analysing structure hierarchy...',... + 'CreateCancelBtn','setappdata(gcbf,''canceling'',1)'); + +for i = 1:height(Tapdvml_contacts_this) % SLOW + waitbar( i/height(Tapdvml_contacts_this), f, ... + sprintf("%d of %d", i, height(Tapdvml_contacts_this))); + name = Tapdvml_contacts_this.name{i}; + [anc(i)] = find_name_for_depth(st, depth_level, name); +end + +delete(f) +beep(); + +Tanc = struct2table(anc, AsArray=true); + +end \ No newline at end of file From 8effe0d8ee0783df370b8f3381c89f892f8d926a Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Thu, 16 Nov 2023 17:03:04 +0000 Subject: [PATCH 38/43] Update updown_probe_with_slider.m --- Browsing Functions/updown_probe_with_slider.m | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/Browsing Functions/updown_probe_with_slider.m b/Browsing Functions/updown_probe_with_slider.m index 99665fa..e855aaa 100644 --- a/Browsing Functions/updown_probe_with_slider.m +++ b/Browsing Functions/updown_probe_with_slider.m @@ -110,6 +110,12 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d T_Y_shift = S.T_Y_shift; YLim_shift_mm = S.T_Y_shift.YLim_shift_mm(probe_id); Tanc = S.T_Y_shift.Tanc{probe_id}; + if isempty(Tanc) + Tanc = compute_Tanc(Tapdvml_contacts_this, st, depth_level); + T_Y_shift.Tanc{probe_id} = Tanc; + save(Y_shift_path, 'T_Y_shift'); + end + clear S else % cache file is not found From ef21e009bd34aafc2cf621055cdff7d56a4d2fbd Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Fri, 17 Nov 2023 09:15:36 +0000 Subject: [PATCH 39/43] Update updown_probe_with_slider.m --- Browsing Functions/updown_probe_with_slider.m | 87 +++++++++++++++++-- 1 file changed, 78 insertions(+), 9 deletions(-) diff --git a/Browsing Functions/updown_probe_with_slider.m b/Browsing Functions/updown_probe_with_slider.m index e855aaa..d6c1a56 100644 --- a/Browsing Functions/updown_probe_with_slider.m +++ b/Browsing Functions/updown_probe_with_slider.m @@ -4,6 +4,56 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d % Move up and down a probe using a slider to match the anatomical borders % to the changes in the firing rates of neurons % +% SYNTAX +% updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_dir,... +% task_dir_name, imaging_session_dir, image_folder_name, depth_level, ... +% active_probe_length, probe_insertion_direction) +% +% +% INPUT ARGUMENTS +% av uint16 +% part of Allen CCF data +% +% st table +% part of Allen CCF data +% +% session_id string +% eg. "kms058-2023-03-25-184034" +% +% probeAB "ProbeA" | "ProbeB" +% +% plane "coronal" | "sagittal" | "transverse" +% +% sessions_dir +% folder path +% "\\ettina\Magill_Lab\Julien\Data\head-fixed\by_sessions" +% +% task_dir_name +% folder name +% "reaching_go_spout_bar_nov22" +% +% imaging_session_dir +% folder_path +% "\\ettina\Magill_Lab\Kouichi Nakamura\Analysis\Images from Otto\20230406 kms058" +% +% image_folder_name +% "RGB" (default) +% folder_name +% +% depth_level +% 6 (default) +% The hierarchical level of the anatomical structures to be +% shown. +% +% active_probe_length +% 3.840 (default) +% +% probe_insertion_direction +% "down" (default) | "up" +% +% +% +% % EXAMPLE USAGE % % session_id = "kms058-2023-03-25-184034" @@ -21,6 +71,16 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d % st = loadStructureTree(structure_tree_location); % tv = readNPY(template_volume_location); % end +% +% +% Written by Kouichi C. Nakamura Ph.D. +% MRC Brain Network Dynamics Unit +% University of Oxford +% kouichi.c.nakamura@gmail.com +% 17-Nov-2023 09:07:25 +% +% See also +% plot_and_compute_probe_positions_from_ccf arguments av uint16 @@ -241,17 +301,30 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d cla(uax2,'reset') movingaverageval = zeros(384,1); +celldensity = zeros(384,1); for i = 1:384 % moving average of five channels - chans = i - 1 : i + 3; + chans = i - 2 : i + 2; + ind = ismember(recording_cell_metrics.maxWaveformCh, chans); + celldensity(i) = nnz(ind)/(length(chans)*10); % cells per micrometers movingaverageval(i) = mean(recording_cell_metrics.firingRate(ind)); end -isc1 = imagesc(uax2, 1, (1:384)'*0.01, movingaverageval); +for i = 1:384 + + % ten channels + chans = i - 4 : i + 5; + + ind = ismember(recording_cell_metrics.maxWaveformCh, chans); + celldensity(i) = nnz(ind)/(length(chans)*10); % cells per micrometers + +end + +isc1 = imagesc(uax2, [0 1], (1:384)'*0.01, movingaverageval); uax2.YDir = 'normal'; ylim(uax2, [-0.100, 3.940]) @@ -261,18 +334,14 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d uax3.CLim = uax2.CLim; % same color scale - - -% Tapdvml_contacts_this = Tapdvml_contacts(Tapdvml_contacts.probe_id == probe_id, :); +ln1 = line(uax2, celldensity/max(celldensity), (1:384)'*0.01, 'Color', 'w'); +xlim(uax2,[0 1]) tx1 = annotation(ufig, 'textbox', [0.2 0.85 0.8 0.1], 'String', ... sprintf("%s/%s (%d)", session_id, probeAB, probe_id),... HorizontalAlignment='center',FontSize=18,LineStyle='none'); -% Tapdvml_contacts_this = [Tapdvml_contacts_this, Tanc]; %TODO - - -% Sample data +%% anc_id = Tapdvml_contacts_this.anc_id; anc_name = Tapdvml_contacts_this.anc_name; From 22a538867e6b26cad9aae6d60d46b0df77ca3733 Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Mon, 20 Nov 2023 15:04:33 +0000 Subject: [PATCH 40/43] Update updown_probe_with_slider.m --- Browsing Functions/updown_probe_with_slider.m | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/Browsing Functions/updown_probe_with_slider.m b/Browsing Functions/updown_probe_with_slider.m index d6c1a56..b82c23d 100644 --- a/Browsing Functions/updown_probe_with_slider.m +++ b/Browsing Functions/updown_probe_with_slider.m @@ -91,7 +91,7 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d sessions_dir (1,1) string {mustBeFolder} % eg. "\\ettina\Magill_Lab\Julien\Data\head-fixed\by_sessions" task_dir_name (1,1) string % eg. "reaching_go_spout_bar_nov22" imaging_session_dir (1,1) string {mustBeFolder} % eg. "\\ettina\Magill_Lab\Kouichi Nakamura\Analysis\Images from Otto\20230406 kms058" - image_folder_name (1,1) string = "RGB" + image_folder_name (1,1) string {mustBeFolder} = "RGB" depth_level (1,1) double {mustBePositive, mustBeInteger} = 6 active_probe_length (1,1) double = 3.84 % in mm probe_insertion_direction (1,1) string {mustBeMember(probe_insertion_direction, ["down","up"])} = "down" @@ -546,9 +546,18 @@ function save_updown_to_table(uax1, ~) assert(nnz(tf) == 1) - p = T_probes{T_probes.session_id == session_id & ... - T_probes.probe_AB == probeAB_, ["p_1","p_2","p_3"]}; - + try + p = T_probes{T_probes.session_id == session_id & ... + T_probes.probe_AB == probeAB_, ["p_1","p_2","p_3"]}; + catch ME + if ME.identifier == "MATLAB:table:UnrecognizedVarName" + error("It looks like m you have to run plot_and_compute_probe_positions_from_ccf()" ... + + " and then save T_probes with p and m values. " + ... + "See scr20230727_112451_SharpTrackIgnore_kms058.mlx for example.") + else + throw(ME); + end + end m = T_probes{T_probes.session_id == session_id & ... T_probes.probe_AB == probeAB_, ["m_1","m_2","m_3"]}; From a2f6bd4cde8d0df582eb0812d3e7db2c6ded2529 Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Mon, 20 Nov 2023 16:06:21 +0000 Subject: [PATCH 41/43] Update updown_probe_with_slider.m --- Browsing Functions/updown_probe_with_slider.m | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Browsing Functions/updown_probe_with_slider.m b/Browsing Functions/updown_probe_with_slider.m index b82c23d..310a916 100644 --- a/Browsing Functions/updown_probe_with_slider.m +++ b/Browsing Functions/updown_probe_with_slider.m @@ -91,7 +91,7 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d sessions_dir (1,1) string {mustBeFolder} % eg. "\\ettina\Magill_Lab\Julien\Data\head-fixed\by_sessions" task_dir_name (1,1) string % eg. "reaching_go_spout_bar_nov22" imaging_session_dir (1,1) string {mustBeFolder} % eg. "\\ettina\Magill_Lab\Kouichi Nakamura\Analysis\Images from Otto\20230406 kms058" - image_folder_name (1,1) string {mustBeFolder} = "RGB" + image_folder_name (1,1) string = "RGB" depth_level (1,1) double {mustBePositive, mustBeInteger} = 6 active_probe_length (1,1) double = 3.84 % in mm probe_insertion_direction (1,1) string {mustBeMember(probe_insertion_direction, ["down","up"])} = "down" @@ -113,6 +113,8 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d image_folder =fullfile(imaging_session_dir, image_folder_name); +assert(isfolder(image_folder), "Check the name of image_folder_name") + probe_save_name_suffix = '_probe'; Tapdvml_contacts_path = fullfile(imaging_session_dir, "Tapdvml_contacts.xlsx"); From 95cf48a4aa34aa5e37cd60a63aeb52035d78e73e Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Wed, 22 Nov 2023 16:41:18 +0000 Subject: [PATCH 42/43] Update updown_probe_with_slider.m --- Browsing Functions/updown_probe_with_slider.m | 38 +++++++++++++++++-- 1 file changed, 34 insertions(+), 4 deletions(-) diff --git a/Browsing Functions/updown_probe_with_slider.m b/Browsing Functions/updown_probe_with_slider.m index 310a916..2b7c001 100644 --- a/Browsing Functions/updown_probe_with_slider.m +++ b/Browsing Functions/updown_probe_with_slider.m @@ -19,6 +19,8 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d % % session_id string % eg. "kms058-2023-03-25-184034" +% "JC318L-2022-12-14-164509, JC318L-2022-12-14-174050" +% multiple sessions, separated by commas % % probeAB "ProbeA" | "ProbeB" % @@ -97,11 +99,22 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d probe_insertion_direction (1,1) string {mustBeMember(probe_insertion_direction, ["down","up"])} = "down" end +%TODO in order to move upwards, we need to load strucural data for the +% region below the tip of the probe, eg. 1 mm beyond the tip +% This will need p and m and retrieving the anatomical labels from st. + + %TODO add button to go next and previous probe %TODO what if optic fiber data is chosen? +% accept multiple sessions separated by comma +session_ids = strip(strsplit(session_id,",")); + +if length(session_ids) > 1 + session_id = session_ids(1); +end sorter_output_dir = fullfile(sessions_dir, task_dir_name, ... @@ -339,10 +352,15 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d ln1 = line(uax2, celldensity/max(celldensity), (1:384)'*0.01, 'Color', 'w'); xlim(uax2,[0 1]) -tx1 = annotation(ufig, 'textbox', [0.2 0.85 0.8 0.1], 'String', ... - sprintf("%s/%s (%d)", session_id, probeAB, probe_id),... - HorizontalAlignment='center',FontSize=18,LineStyle='none'); - +if length(session_ids) > 1 + tx1 = annotation(ufig, 'textbox', [0.2 0.85 0.8 0.1], 'String', ... + sprintf("%s/%s (%d)", session_ids, probeAB, probe_id),... + HorizontalAlignment='center',FontSize=18,LineStyle='none'); +else + tx1 = annotation(ufig, 'textbox', [0.2 0.85 0.8 0.1], 'String', ... + sprintf("%s/%s (%d)", session_id, probeAB, probe_id),... + HorizontalAlignment='center',FontSize=18,LineStyle='none'); +end %% anc_id = Tapdvml_contacts_this.anc_id; @@ -361,6 +379,9 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d block_colors = anc_color(start_ind, :); +%TODO show names for regions below the tip as well +% may take 1- 2 hrs + Tblocks = table(block_names, start_ind, end_ind, block_colors, ... VariableNames={'block_names', 'start_ind','end_ind', 'block_colors'}); @@ -381,6 +402,8 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d Tblocks.block_names(i),... 'FontSize',14,... 'HorizontalAlignment','center', Clipping='on'); + + %TODO show names for regions below the tip as well end horzline = line(uax1, [0 1], [0 0], 'Color','r'); % distance from the tip 0 mm @@ -500,6 +523,13 @@ function save_updown_to_table(uax1, ~) end end + if length(session_ids) > 1 + %TODO + warning('Save for multilpe sessions not implemented yet.') + keyboard + + end + %% save the backup From eda3c21e40eae1a5b3c8be36890c6e6148cf412a Mon Sep 17 00:00:00 2001 From: "Kouichi C. Nakamura" Date: Tue, 28 Nov 2023 15:48:32 +0000 Subject: [PATCH 43/43] Update updown_probe_with_slider.m --- Browsing Functions/updown_probe_with_slider.m | 26 ++++++------------- 1 file changed, 8 insertions(+), 18 deletions(-) diff --git a/Browsing Functions/updown_probe_with_slider.m b/Browsing Functions/updown_probe_with_slider.m index 2b7c001..49deb6b 100644 --- a/Browsing Functions/updown_probe_with_slider.m +++ b/Browsing Functions/updown_probe_with_slider.m @@ -161,10 +161,14 @@ function updown_probe_with_slider(av, st, session_id, probeAB, plane, sessions_d case "ProbeB" probeAB_short = "B"; otherwise - + probeAB_short = "optic fiber"; end -tf = T_probes.session_id == session_id & T_probes.probe_AB == probeAB_short; +T_probes.session_ids = (cellfun(@(x) string(strsplit(x, ', ')), ... + T_probes.session_id, UniformOutput=false)); % handling multiple session_ids for one probe + +tf = cellfun(@(x) ismember(session_id, x), T_probes.session_ids) ... + & T_probes.probe_AB == probeAB_short; probe_id = T_probes{tf, "probe_id"}; @@ -563,24 +567,11 @@ function save_updown_to_table(uax1, ~) % see plot_and_compute_probe_positions_from_ccf.m % I need the eigen vector, to be saved in T_probes? - switch probeAB - case "ProbeA" - probeAB_ = "A"; - case "ProbeB" - probeAB_ = "B"; - case "optic fiber" - probeAB_ = "optic fiber"; - end - % eigen vector - tf = T_probes.session_id == session_id & ... - T_probes.probe_AB == probeAB_; - assert(nnz(tf) == 1) try - p = T_probes{T_probes.session_id == session_id & ... - T_probes.probe_AB == probeAB_, ["p_1","p_2","p_3"]}; + p = T_probes{tf, ["p_1","p_2","p_3"]}; catch ME if ME.identifier == "MATLAB:table:UnrecognizedVarName" error("It looks like m you have to run plot_and_compute_probe_positions_from_ccf()" ... @@ -590,8 +581,7 @@ function save_updown_to_table(uax1, ~) throw(ME); end end - m = T_probes{T_probes.session_id == session_id & ... - T_probes.probe_AB == probeAB_, ["m_1","m_2","m_3"]}; + m = T_probes{tf, ["m_1","m_2","m_3"]}; processed_images_folder = fullfile(image_folder, 'processed');