diff --git a/Browsing Functions/AtlasTransformBrowser.m b/Browsing Functions/AtlasTransformBrowser.m index 5cf501c..553fb4f 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. % ------------------------------------------------ @@ -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; @@ -69,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; @@ -79,6 +84,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'); @@ -176,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 @@ -348,18 +367,38 @@ 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 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'); + else; disp('transform point mode OFF'); + % hide transform points + + try + set(ud.pointHands_for_transform(:), 'Visible', 'off'); + ud.pointsText.Visible = 'off'; + catch + end + end % a -- toggle viewing of annotation boundaries case 'a' @@ -402,37 +441,56 @@ 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 + 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; 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(' '); @@ -491,11 +549,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 @@ -588,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 @@ -629,20 +695,41 @@ 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,:); + 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 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 end 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 ud.pointList{ud.currentProbe,1} = ud.pointList{ud.currentProbe,1}(1:end-1,:); @@ -693,8 +780,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'])) @@ -732,7 +820,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'); @@ -774,12 +865,30 @@ 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 = []; + 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') + 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}; @@ -793,6 +902,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'); @@ -800,10 +912,14 @@ 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'); - 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 @@ -834,7 +950,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); @@ -901,7 +1020,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 @@ -1060,7 +1179,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; @@ -1068,6 +1190,24 @@ 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)); + + 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); diff --git a/Browsing Functions/accf2pxs_mm.m b/Browsing Functions/accf2pxs_mm.m new file mode 100644 index 0000000..37b107e --- /dev/null +++ b/Browsing Functions/accf2pxs_mm.m @@ -0,0 +1,84 @@ +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 + +arguments + x_mm + plane (1,1) string {mustBeMember(plane, {'coronal','sagittal','transverse'})} + numtype (1,1) string {mustBeMember(numtype, {'coordinate','distance'})} = "coordinate" + +end + + +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/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..2ddb266 --- /dev/null +++ b/Browsing Functions/apdvml2info.m @@ -0,0 +1,94 @@ +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 +% 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 + +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 +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, 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 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/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/plot_and_compute_probe_positions_from_ccf.m b/Browsing Functions/plot_and_compute_probe_positions_from_ccf.m new file mode 100644 index 0000000..16b6996 --- /dev/null +++ b/Browsing Functions/plot_and_compute_probe_positions_from_ccf.m @@ -0,0 +1,400 @@ +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) +% 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_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) +% +% 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 +% +% 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 +% probes +% +% fwireframe figure +% +% fig_probes figure +% +% Tapdvml_contacts +% table +% A table containig whereabouts of all the recording contacts +% +% 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 +% kouichi.c.nakamura@gmail.com +% 09-Aug-2023 11:12:42 +% +% See also +% apdvml2info + + + +% 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]; + +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 +% 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'; + +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),... + strings(t_size), ... + 'VariableNames',{'subject_id','session_id','probe_id','probe_AB',... + 'depth_reported_um','depth_from_ACCF_um','scaling', ... + '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; + + +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; + 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); + % 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'; + +p = P; %rename +m = M; \ 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..aaf8924 --- /dev/null +++ b/Browsing Functions/prep_Tapdvml_contacts.m @@ -0,0 +1,85 @@ +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.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)); + +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 diff --git a/Browsing Functions/sliceBrowser.m b/Browsing Functions/sliceBrowser.m index e190b6e..45cf0aa 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,12 @@ 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'; + +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)); @@ -38,6 +44,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 % ------------------------------------------------ @@ -45,21 +54,37 @@ 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)); 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); - - if clickX < 100 && (ud.ref_size(1) - clickY) < 100 % if click in corner, break - ud.pointList = []; - set(ud.pointHands(:), 'Visible', 'off'); - end + % 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); + + % 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); +end % ------------------------ % react to keyboard press @@ -88,14 +113,41 @@ 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'); %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); + + % 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]; + 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)); + + 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); @@ -105,9 +157,10 @@ function SliceAtlasHotkeyFcn(fig, keydata, f) set(fig, 'UserData', ud); - +end function ud = updateSliceImage(ud) + %TODO occasionally pointList has more items than pointHands title_ending = ''; @@ -127,7 +180,11 @@ 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','Marker','o')); % delete all the circles + end + + ud.pointHands = gobjects(0); ud.pointList = []; if exist(file_transformations,'file') @@ -136,7 +193,61 @@ 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 + set(ud.pointHands(end),'color', [0 .9 0]); + title_ending = ' (transform points loaded)'; end end - title(['Slice Viewer -- Slice ' num2str(ud.slice_num) '/' num2str(ud.total_num_files) title_ending]) + 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]) + +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/Browsing Functions/updown_probe_with_slider.m b/Browsing Functions/updown_probe_with_slider.m new file mode 100644 index 0000000..49deb6b --- /dev/null +++ b/Browsing Functions/updown_probe_with_slider.m @@ -0,0 +1,783 @@ +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 +% +% 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" +% "JC318L-2022-12-14-164509, JC318L-2022-12-14-174050" +% multiple sessions, separated by commas +% +% 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" +% 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 +% +% +% 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 + 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 + +%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, ... + 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; + +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"); + +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"); %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 + probeAB_short = "optic fiber"; +end + +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"}; + +%% 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 ); + 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}; + 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 + 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 ); + Tapdvml_contacts_this = Tapdvml_contacts(Tapdvml_contacts.probe_id == probe_id, :); + + 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'); + + Tapdvml_contacts_this = [Tapdvml_contacts_this, Tanc]; + +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 =[]; + +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]; +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]; + + +%% 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); +celldensity = zeros(384,1); + +for i = 1:384 + + % moving average of five channels + 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 + +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]) +uax2.XTick =[]; + +uax2.Color = 'k'; + +uax3.CLim = uax2.CLim; % same color scale + +ln1 = line(uax2, celldensity/max(celldensity), (1:384)'*0.01, 'Color', 'w'); +xlim(uax2,[0 1]) + +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; + +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, :); + +%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'}); + + +%% 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'); + + %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 + +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 + + 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; + + 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) + + if isfield(uax1.UserData, 'YLim_shift_mm') + ch = uax1.Children; + + 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; + + 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 + + + function move_vertically(uax1, step) + + % arguments + % uax1(1,1) + % step (1,1) {mustBeInteger} + % end + + ch = uax1.Children; + + for j = 1:length(ch) + + 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.Patch') + ch(j).YData = ch(j).YData + step*0.01; + + elseif isa(ch(j),'matlab.graphics.primitive.Line') + ch(j).YData = ch(j).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 + + tx2.String = regexprep(sprintf("%+.2f mm", uax1.UserData.YLim_shift_mm), '-', char(8722)); + + 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') + % move + else + disp('YLim_shift_mm is not found. No change made.') + return + end + end + + if length(session_ids) > 1 + %TODO + warning('Save for multilpe sessions not implemented yet.') + keyboard + + end + + + %% 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"); + + if isfile(original_backup_path) + % save a backup file for the latest + + 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? + + % eigen vector + assert(nnz(tf) == 1) + + try + 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()" ... + + " 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{tf, ["m_1","m_2","m_3"]}; + + 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{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 + + % m is the brain entry point + + % 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 + %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 + + 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 + ... + new_cp * (tip2surface_mm - active_probe_length))... + /tip2surface_mm; + + %% 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)', ... + 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 + + c_t_contacts = cell(384,1); + + theta = acos(p(2)); % Assuming p(2) corresponds to the DV direction + + 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'); + + c_t_contacts{j}.depth_mm_paxinos = projected_depth_mm_paxinos; + end + + Tapdvml_contacts_new = vertcat(c_t_contacts{:}); + clear c_t_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"}; + + 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'); + + + 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 "" + +for j = 1:size(cstr,1) + + 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 + + +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 diff --git a/Histology Functions/HistologyBrowser.m b/Histology Functions/HistologyBrowser.m index c8058f4..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]) @@ -41,6 +49,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 +105,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 +121,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 +150,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 +171,64 @@ function HistologyHotkeyFcn(fig, keydata, image_file_names, use_already_downsamp set(fig, 'UserData', ud); +end + +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) + 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"]; + + 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 + + + +end \ No newline at end of file diff --git a/Histology Functions/SliceFlipper.m b/Histology Functions/SliceFlipper.m index c8ba7e8..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)); @@ -52,12 +60,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 +146,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 +169,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 +186,63 @@ function SliceScrollFcn(fig, evt) set(fig, 'UserData', ud); +end + +end + +function out = 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"]; + +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) + 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') +out = 1; +end