%
%
%  Filament Angles 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 Angles Statistics" icon="Matlab" tooltip="Plot filament angles statistics.">
%          <Command>Matlab::BPFilamentAnglesStatistics(%i)</Command>
%        </Item>
%       </Submenu>
%      </Menu>
%      <SurpassTab>
%        <SurpassComponent name="bpFilament">
%          <Item name="Angles Statistics" icon="Matlab" tooltip="Plot filament angles statistics.">
%            <Command>Matlab::BPFilamentAnglesStatistics(%i)</Command>
%          </Item>
%        </SurpassComponent>
%      </SurpassTab>
%    </CustomTools>
% 
% 
%  Description:
%   
%   Plot two histograms that display the statistic of the angles of a filament.
%   The angles are divided in two types: Angles between the segments and
%       angles between the branches. 
%
%

function BPFilamentAnglesStatistics(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);


vAngle = [];
vBranchAngles = [];

for vSpots = 1:vNumberOfSpots

    [vEdge, vSide] = find(vFilamentEdges==vSpots-1);

    vPrevious = -1;
    vNext = [];
    for vNeighborIndex = 1:length(vEdge)
        vNeighbor = vFilamentEdges(vEdge(vNeighborIndex), 3-vSide(vNeighborIndex))+1;
        if vNeighbor<vSpots
            vPrevious = vNeighbor;
        else
            vNext = [vNext,vNeighbor];            
        end
    end
    
    % angle between previous and next(s)
    if vPrevious>-1 && length(vNext)>0
        vVector1 = vFilamentXYZ(vSpots,:)-vFilamentXYZ(vPrevious,:);
        vNorm1 = norm(vVector1);
        for vSpot = 1:length(vNext)
            vVector2 = vFilamentXYZ(vNext(vSpot),:)-vFilamentXYZ(vSpots,:);
            vAngle = [vAngle, acos(dot(vVector1,vVector2)/(vNorm1*norm(vVector2)))];
        end
    end

    % angle between next i and next j
    if length(vNext)>1
        for vSpot1 = 1:length(vNext)-1
            vVector1 = vFilamentXYZ(vNext(vSpot1),:)-vFilamentXYZ(vSpots,:);
            vNorm1 = norm(vVector1);
            for vSpot2 = vSpot1+1:length(vNext)
                vVector2 = vFilamentXYZ(vNext(vSpot2),:)-vFilamentXYZ(vSpots,:);
                vBranchAngles = [vBranchAngles, ...
                    acos(dot(vVector1,vVector2)/(vNorm1*norm(vVector2)))];
            end
        end    
    end
    
    waitbar(vSpots/vNumberOfSpots);
end

close(vProgressDisplay);

vAngle = vAngle*180/pi;
vBranchAngles = vBranchAngles*180/pi;

subplot(2,1,1);
vSize = 2;
hist(vAngle,vSize/2+vSize*(0:max(ceil(vAngle/vSize))));
title('Angles between segments');
xlabel('Angle');
ylabel('# of segments');


subplot(2,1,2);
vSize = 5;
hist(vBranchAngles,vSize/2+vSize*(0:max(ceil(vBranchAngles/vSize))));
title('Angles between branches');
xlabel('Angle');
ylabel('# of branches');