calcDensity.m

%% ObjectFinder - Recognize 3D structures in image stacks
%  Copyright (C) 2016,2017,2018 Luca Della Santina
%
%  This file is part of ObjectFinder
%
%  ObjectFinder is free software: you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation, either version 3 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program.  If not, see <https://www.gnu.org/licenses/>.
%
function Density = calcDensity(Settings, Dots, Skel, showPlot)
Dots = getFilteredObjects(Dots, Dots.Filter); % Work only on validated objects

% Mids are the middle point of each skeleton's segment
AllSegCut = cat(2, Skel.SegStats.Seg(:,2,:), ...
                   Skel.SegStats.Seg(:,1,:), ...
                   Skel.SegStats.Seg(:,3,:));

DPos=round(Dots.Pos);
DPos(:,1:2)=DPos(:,1:2)*Settings.ImInfo.xyum; 
DPos(:,3)=DPos(:,3)*Settings.ImInfo.zum;

%Extract Dend positions
Mids=mean(AllSegCut,3); % segment xyz position calculated as mean of two node positions
Length=sqrt((AllSegCut(:,1,1)-AllSegCut(:,1,2)).^2 ...
          + (AllSegCut(:,2,1)-AllSegCut(:,2,2)).^2 ...
          + (AllSegCut(:,3,1)-AllSegCut(:,3,2)).^2);

Nearest = zeros(size(DPos, 1),1);
for i = 1:size(DPos,1)
    Ndist=dist2(Mids,DPos(i,:)); %find dist from dot to all nodes
    Near=min(Ndist); %find shortest distance
    Nearest(i)=find(Ndist==Near,1); %get node at that distance
end
NN=Mids(Nearest,:); %assign that node to NearestNode list for dots

% Image is shrunck so that number of 
yNumVox=fix(Settings.ImInfo.yNumVox*Settings.ImInfo.xyum)+1; % yNumVox is the rounded um of the image in y direction.
xNumVox=fix(Settings.ImInfo.xNumVox*Settings.ImInfo.xyum)+1; % xNumVox is the rounded um of the image in y direction.
Msize=[yNumVox xNumVox];     % 2048*2048 with 0.103um xy pixel size becomes Settings.ImInfo.yNumVox=Settings.ImInfo.xNumVox=211.

% create Nearest Node list (nearest node for each dot and distance)
Mids(Mids<1)=1;
NN(NN<1)=1;
Mids(Mids(:,1)>yNumVox,1)=yNumVox;
NN(NN(:,1)>yNumVox,1)=yNumVox;
Mids(Mids(:,2)>xNumVox,2)=xNumVox;
NN(NN(:,2)>xNumVox,2)=xNumVox;

clear sMids sNN
sMids=Mids;
sNN=NN;

sMids=round(sMids); %This rounds up all the sMids to 1um step.
sNN=round(sNN); %%This rounds up all the sNN to 1um step.

% Generates DotMap and DendMap at 1um resolution. DotMap and
% DendMap will have pixels of values >1 if the number of dots or the
% length of arbor found in the 1um pixel field is >1. In other words,
% this part is calculating the density of dots (#/um2) or the density of
% arbors (um/um2) in each 1um pixel area.

DotMap=zeros(Msize(1),Msize(2)); %DotMap will be 211*211 if the image is 2048*2048 with 0.103um xy pixel size.
for i=1:size(sNN,1)
    DotMap(sNN(i,1),sNN(i,2))=DotMap(sNN(i,1),sNN(i,2))+1; %if you have >1 dot within the 1um pixel field, you gain more value.
end

DendMap=zeros(Msize(1),Msize(2)); %DotMap will be 211*211 if the image is 2048*2048 with 0.103um xy pixel size.
for i=1:size(sMids,1)
    DendMap(sMids(i,1),sMids(i,2))=DendMap(sMids(i,1),sMids(i,2))+Length(i); %if you have >1 arbor within the 1um pixel field, you gain more value. Also, if the arbor is long, you gain more value.
end

% Generate DotFilt and DendFilt by convolving a disk (10um radius) and
% averaging DotMap and DendMap, respectively. 
% So, instead of each 1um pixel representing the density of that particular
% pixel, it represents the average density of dots within the 10um area

AreaS=10; %filter is 10um RADIUS disk (because one pixel is 1um*1um).
Disk=fspecial('disk',AreaS); %fspecial geneartes averaging filter. It averages the pixel values in the 21*21 circular area.
DotFilt=imfilter(DotMap,Disk,'same');
DendFilt=imfilter(DendMap,Disk,'same');

% Find territory
% This part generates Dendritic territory, which is the convolution of
% 5um radius (smaller than DendFilt and DotFilt!) disk averaging filter
% with DendMap and all the pixels within the territory converted to 1
% and those outside the territory converted to 0, so the territory acts
% like a filter when used in the calculation of DotDist and DendDist.
% DotDist and DendDist is generated by filtering DotFilt and DendFilt
% with Territory. This is the P/A and D/A GC figures in Josh's paper,
% and the element-wise DotDist/DendDist is the P/D GC figure. Also mean
% P/A, mean D/A and mean P/D were calculated by averaging all the pixels
% WITHIN THE TERRITORY.

Disk2=fspecial('disk',5); % For territory, use 5um RADIUS averaging disk filter.
Territory=imfilter(DendMap,Disk2,'same');
DotDist = DotFilt.*Territory;
DendDist = DendFilt.*Territory;

Density.DotMap = DotMap;
Density.DendMap = DendMap;
Density.Territory = Territory;
Density.DotDist = DotDist;
Density.DendDist = DendDist;

if showPlot
    figure;
    cmap=jet(256);
    cmap(1,:)=0;
    colormap(cmap);
    
    % Draw Dot maps
    image(DotMap*200);
    title('Skeletonized dot density map (blue-red = 0-1.28puncta/um2) (1um pixel size)');
    pause(2);    
    % Draw Dend map
    image(DendMap*100);
    title('Skeletonized dendrite density map (blue-red = 0-2.55um/um2) (1um pixel size)');
    pause(2);    
    % Draw Territory with perimeter highlighted
    image(Territory*200);
    title('Territory (dendrite skeleton filtered by 5um radius disk)');
    pause(2);
    image((Territory+bwperim(Territory,8))*100);
    title('Perimeter highlighted');
    pause(2);    
    % Drap heatmap of linear density (P/D)
    image(DotDist./DendDist*256);
    title('Puncta/Dendrite (linear density) (0-1 dots/µm)');
    colorbar('Ytick',[]); %Adam added color bar 3/11/2011
    axis image; % adding color bar changes image x y dimention, so resize the image.
    axis off    % remove the axis
    saveas(gcf,[Settings.TPN 'images' filesep 'HeatmapLinearDensity.tif']);
end
end

function[d]=dist2(A,B)
% Calculate distance between vectors A (n,3,n) and B (1,3)
A2=zeros(size(A,1),3,size(A,3));
B2=zeros(size(B,1),3);
A2(:,1:size(A,2),:)=A;
B2(:,1:size(B,2))=B;
A=A2;
B=B2;
d=sqrt((A(:,1,:)-B(1)).^2 + (A(:,2,:)-B(2)).^2 + (A(:,3,:)-B(3)).^2);
end