adding Point counters

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
+% [email protected]
+% 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

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

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
+% [email protected]
+% 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

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
+% [email protected]
+% 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

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
+% [email protected]
+% 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