%
%
%  Track Plot Angles 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="Tracks Functions">
%        <Item name="Plot Angles of the Track" icon="Matlab" tooltip="Plot the angles of the track.">
%          <Command>Matlab::BPTrackPlotAngles(%i)</Command>
%        </Item>
%       </Submenu>
%      </Menu>
%      <SurpassTab>
%        <SurpassComponent name="bpTrack">
%          <Item name="Plot Angles" icon="Matlab" tooltip="Plot the angles of the track.">
%            <Command>Matlab::BPTrackPlotAngles(%i)</Command>
%          </Item>
%        </SurpassComponent>
%      </SurpassTab>
%    </CustomTools>
% 
%
%  Description:
%   
%   Plot the angles of the track. 
%   
%

function BPTrackPlotAngles(aImarisApplicationID)

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

% the user has to create a scene with some tracks
vSurpassScene = vImarisApplication.mSurpassScene;
if isequal(vSurpassScene, [])
    msgbox('Please create some tracks in the surpass scene!');
    return;
end

% get the track
vTrack = vImarisApplication.mFactory.ToTrack(vImarisApplication.mSurpassSelection);

if ~vImarisApplication.mFactory.IsTrack(vTrack)  
    msgbox('Please select a track!');
    return;
end

% get the spots and edges
%vEdges = vTrack.GetEdges;
vSpots = vTrack.GetSpots;

[vSpotsXYZ,vSpotsTime,vSpotsRadius] = vSpots.Get;
vSpotsXYZ = vSpotsXYZ';

vNumberOfSpots = length(vSpotsRadius);

% a few lines to order eventual unordered data (e.g. vSpotsTime == [1,5,2])
vTimeMin = min(vSpotsTime);
vTimeMax = max(vSpotsTime);

vSpotAtTime = zeros(1,vTimeMax-vTimeMin+1);
vSpotAtTime(vSpotsTime-vTimeMin+1) = 1:vNumberOfSpots;

vTime = sort(vSpotsTime);

vSpotsIndex = vSpotAtTime(vTime-vTimeMin+1);

vStepsSize = vSpotsXYZ(:,vSpotsIndex(2:vNumberOfSpots)) - ...
    vSpotsXYZ(:,vSpotsIndex(1:vNumberOfSpots-1));

% get the angle for each pair of contiguous edges
for vSpot = 1:vNumberOfSpots-2
    vVector1 = vStepsSize(:,vSpot);
    vVector2 = vStepsSize(:,vSpot+1);
    
    vAngleXYZ(vSpot) = GetAngle(vVector1,vVector2);
    vAngleXY(vSpot) = GetAngle(vVector1(1:2),vVector2(1:2));
    vAngleXZ(vSpot) = GetAngle(vVector1([1,3]),vVector2([1,3]));
    vAngleYZ(vSpot) = GetAngle(vVector1(2:3),vVector2(2:3));
end
    
% finally plot the results
vStrings = {'XYZ','XY','XZ','YZ'};
for vPlot = 1:4
    subplot(2,2,vPlot)
    plot(vSpotsTime(2:vNumberOfSpots-1)+1,eval(['vAngle',vStrings{vPlot}])*180/pi,'b-');
    title([vStrings{vPlot}, '-Angles of ', vTrack.mName]);
    xlabel('Time');
    ylabel('Amplitude [degrees]');
    hold off
end


%---------------------------------------------------------%

function aAngle = GetAngle(aVector1, aVector2)
vNormProduct = norm(aVector1) * norm(aVector2);

if vNormProduct == 0
    aAngle = 0;
else
    aAngle = acos(dot(aVector1,aVector2)/vNormProduct);
end