%
%
%  Filament Spines Statistics Function for Imaris
%
%  Copyright Bitplane AG 2006
%
%
%  Installation:
%
%  - Copy this file into the XTensions folder in the Imaris installation directory.
%  - You will find this function in the Image Processing menu
%
%    <CustomTools>
%      <Menu>
%       <Submenu name="Filament Functions">
%        <Item name="Filament Spines Statistics" icon="Matlab" tooltip="Plot filament spines statistics.">
%          <Command>Matlab::BPFilamentSpinesStatistics(%i)</Command>
%        </Item>
%       </Submenu>
%      </Menu>
%      <SurpassTab>
%        <SurpassComponent name="bpFilament">
%          <Item name="Spines Statistics" icon="Matlab" tooltip="Plot filament spines statistics.">
%            <Command>Matlab::BPFilamentSpinesStatistics(%i)</Command>
%          </Item>
%        </SurpassComponent>
%      </SurpassTab>
%    </CustomTools>
% 
%
%  Description:
%
%   Plot three histograms that display the distribution of the length, 
%       the average width and the volume (log scale) of the spines of a
%       filament object. 
%   In addition, a plot shows the distribution of the average width to 
%       respect to the length of the spine. 
%   For "average width" is intended the mean diameter of a single spine. 
%
%

function BPFilamentSpinesStatistics(aImarisApplicationID)

% connect to Imaris Com interface
if ~isa(aImarisApplicationID, 'COM.Imaris_Application')
  vImarisServer = actxserver('ImarisServer.Server');
  vImarisApplication = vImarisServer.GetObject(aImarisApplicationID);
else
  vImarisApplication = aImarisApplicationID;
end

% get the filament
vFilament = vImarisApplication.mFactory.ToFilament(vImarisApplication.mSurpassSelection);

% search the filament if not previously selected
vSurpassScene = vImarisApplication.mSurpassScene;
if ~vImarisApplication.mFactory.IsFilament(vFilament)        
    for vChildIndex = 1:vSurpassScene.GetNumberOfChildren
        vDataItem = vSurpassScene.GetChild(vChildIndex - 1);
        if vImarisApplication.mFactory.IsFilament(vDataItem)
            vFilament = vDataItem;
            break;
        end
    end
    % did we find the filament?
    if isequal(vFilament, [])
        msgbox('Please create some filament!');
        return;
    end
end

% get the filament coordinates
[vFilamentXYZ,vFilamentRadius,vFilamentEdges] = vFilament.Get;

vProgressDisplay = waitbar(0,'Computing statistics');

vNumberOfSpots = length(vFilamentRadius);

% find the terminals
vNumberOfTerminals = 0;
vTerminals = [];
vIsBranch = zeros(1,vNumberOfSpots);

for vSpots = 1:vNumberOfSpots
    [vEdge, vSide] = find(vFilamentEdges==vSpots-1);
    if length(vEdge)==1
        vNumberOfTerminals = vNumberOfTerminals + 1;
        vTerminals(vNumberOfTerminals) = vSpots;
    elseif length(vEdge)>2
        vIsBranch(vSpots) = true;
    end
end

waitbar(0.3);

vLength = 1:vNumberOfTerminals;
vAverageWidth = vLength;
vNumberOfSpots = vLength;

for vTerminalIndex = 1:vNumberOfTerminals
    vTerminal = vTerminals(vTerminalIndex);

    vLength(vTerminalIndex) = 0;
    vAverageWidth(vTerminalIndex) = vFilamentRadius(vTerminal)*2;
    vNumberOfSpots(vTerminalIndex) = 1;

    
    vFound = true;
    while vFound
        [vEdge, vSide] = find(vFilamentEdges==vTerminal-1);
        
        vFound = false;
        for vNeighborIndex = 1:length(vEdge)
            vNeighbor = vFilamentEdges(vEdge(vNeighborIndex), 3-vSide(vNeighborIndex))+1;
            if vNeighbor<vTerminal
                vNewTerminal = vNeighbor;
                vFound = true;
           end
        end

        if vFound
            vLength(vTerminalIndex) = vLength(vTerminalIndex) + ...
                sqrt(sum((vFilamentXYZ(vTerminal,:)-vFilamentXYZ(vNewTerminal,:)).^2));
            vNumberOfSpots(vTerminalIndex) = vNumberOfSpots(vTerminalIndex) + 1;
            vTerminal = vNewTerminal;
            vAverageWidth(vTerminalIndex) = (vAverageWidth(vTerminalIndex)* ...
                (vNumberOfSpots(vTerminalIndex)-1) + vFilamentRadius(vTerminal)*2) / ...
                vNumberOfSpots(vTerminalIndex);
        end
        
        if vIsBranch(vTerminal)
            vFound = false;
        end
    end
    
    waitbar(0.3+0.7*vTerminalIndex/vNumberOfTerminals);
end

close(vProgressDisplay);

vHistResolution = max([5,ceil(sqrt(vNumberOfTerminals)+vNumberOfTerminals/15)]);

subplot(2,2,3);
plot(vLength, vAverageWidth, 'b.');
xlabel('Length');
ylabel('Average width');

subplot(2,2,4);
% V = h * pi*r^2
vVolume = pi * (vAverageWidth/2).^2 .* vLength;

vRobustLog = log(vVolume(vVolume>0.00001));

hist(vRobustLog,vHistResolution);
%semilogy(sort(vLength./vAverageWidth));
xlabel('log(Volume)');
ylabel('# of spines');

vLength = sort(vLength);
vAverageWidth = sort(vAverageWidth);

% subplot(3,2,1);
% plot(vLength)
% xlabel('Spine');
% ylabel('Length');

subplot(2,2,1);
hist(vLength,vHistResolution);
xlabel('Length');
ylabel('# of spines');

title(sprintf('Spines statistics (%i spines)',vNumberOfTerminals));

% subplot(3,2,3);
% plot(vAverageWidth)
% xlabel('Spine');
% ylabel('Average width');

subplot(2,2,2);
hist(vAverageWidth,vHistResolution)
xlabel('Average width');
ylabel('# of spines');