function [] = alignTracks(calling_button, eventData)
%This function finds similar tracks and connects them
%   This function begins by finding a global minima for connecting all tracks'
%   coordinates for consecutive time points for all time points. The
%   function then averages the connected points to a new point. All the new
%   points are connected to new tracks which contain the coordinates of two
%   channels and their intensities for all time points.

alignWaitBar = msgbox('Wow. Look at me. I am aligning tracks all by myself','Program is running');
%notify the user the program is running

trackAlignmentFigureHandle = figure('Position', [200, 500, 250, 350],...
    'tag','trackAlignmentFigure','color', 'w', 'MenuBar','none',...
    'name','Aligned tracks','numbertitle','off','Visible', 'off');

uitable('parent', trackAlignmentFigureHandle,...
        'unit', 'norm',...
        'position',[0.03 0.15 0.9 0.7],... %[left,bottom,width,height]
        'Visible', 'on',...
        'ColumnName', {'Track','Show'},...
        'ColumnFormat', {'numeric','logical'},...
        'ColumnWidth',{70},...
        'tag', 'tracksTable4',...
        'ColumnEditable', [false true],...
        'Data', {});
uicontrol('Style','text',...
    'parent', trackAlignmentFigureHandle,...
    'String','Choose all aligned tracks',...
    'HorizontalAlignment', 'left',...
    'fontsize', 12,...
    'Visible', 'on',...
    'backgroundColor','w',...
    'foregroundColor','r',...
    'unit', 'norm',...
    'Position',[0.12 0.86 0.9 0.05],... %[left,bottom,width,height]
    'tag','showTrack4Title');
uicontrol('Style','checkbox',...
    'parent', trackAlignmentFigureHandle,...
    'Visible', 'on',...
    'Value', 1,...
    'backgroundColor','w',...
    'foregroundColor','r',...
    'unit', 'norm',...
    'Position',[0.01 0.86 0.1 0.05],... %[left,bottom,width,height]
    'tag','showAllTracks4',...
    'callback', {@showAll, 'tracksTable4'});

uicontrol('style'   ,'checkbox',....
    'parent', trackAlignmentFigureHandle,...
    'unit', 'norm',...
    'string','Enable plot aligned tracks',...
    'fontsize', 12,...
    'foregroundColor','r',...
    'position',[0.05 0.025 0.9 0.075],... %[left,bottom,width,height]
    'Enable', 'on',...
    'tag','toggleAlignedTracks');

%find user defined parametres
maxLinkingDistanceHandle = findobj('tag','maxLinkingDistance');
maxLinkingDistance = get(maxLinkingDistanceHandle, 'value');
if maxLinkingDistance < 0 || mod(maxLinkingDistance,1) ~= 0
    close(waitBar);
    errordlg('Max linking distance has to be an integer greater than or equal to 0', 'False initial parametres');
    error('Max linking distance has to be an integer greater than or equal to 0');
end

pixelSizeHandle = findobj('tag','pixelSize');
pixelSize = get(pixelSizeHandle, 'value');
if (pixelSize > 0) == 0
    close(alignWaitBar);
    errordlg('Pixel size has to be greater than 0', 'False initial parametres');
    error('Pixel size has to be greater than 0');
end
%look up user input for pixel size

stepSizeHandle = findobj('tag','stepSize');
stepSize = get(stepSizeHandle, 'value');
if (stepSize > 0) == 0
    close(alignWaitBar);
    errordlg('Step size has to be greater than 0', 'False initial parametres');
    error('Step size has to be greater than 0');
end
%look up user input for step size

minimumTrackLengthHandle = findobj('tag','minTrackLength');
minimumTrackLength = get(minimumTrackLengthHandle, 'value');
if minimumTrackLength < 0 || mod(minimumTrackLength,1) ~= 0
    close(waitBar);
    errordlg('Minimal track length has to be an integer greater than or equal to 0', 'False initial parametres');
    error('Minimal track length has to be an integer greater than or equal to 0');
end
%look up user input for minimum acceptable track length

useChannel1Handle = findobj('tag','chooseChannel1');
useChannel1 = get(useChannel1Handle, 'value');
useChannel2Handle = findobj('tag','chooseChannel2');
useChannel2 = get(useChannel2Handle, 'value');
useChannel3Handle = findobj('tag','chooseChannel3');
useChannel3 = get(useChannel3Handle, 'value');
%look up input image channles from user

%find which channels the user tracked
tracksVector = cell(2,3);
countChannels = 0;
if useChannel1 == 1
Channel1Handle = findobj('tag','Channel1');
imageMatrixName = get(Channel1Handle,'string');
if any(strcmp(evalin('base','who'),imageMatrixName)) == 0
    close(alignWaitBar);
    errordlg('Channel not found. Please match channel name to an existing file or variable', 'Channel selection error');
    error('Channel not found. Please match channel name to an existing file or variable');
end
imageMatrixName = ['finalTracks_' imageMatrixName];
if any(strcmp(evalin('base','who'),imageMatrixName)) == 0
    close(alignWaitBar);
    errordlg('Channel track not found. Please run program again for appropriate channels', 'Channel selection error');
    error('Channel track not found. Please run program again for appropriate channels')
end
%find data set
data = evalin('base',imageMatrixName);
%save data set
tracksVector{1,1} = data;
tracksVector{2,1} = imageMatrixName;
countChannels = countChannels + 1;
end
if useChannel2 == 1
Channel2Handle = findobj('tag','Channel2');
imageMatrixName = get(Channel2Handle,'string');
if any(strcmp(evalin('base','who'),imageMatrixName)) == 0
    close(alignWaitBar);
    errordlg('Channel not found. Please match channel name to an existing file or variable', 'Channel selection error');
    error('Channel not found. Please match channel name to an existing file or variable');
end
imageMatrixName = ['finalTracks_' imageMatrixName];
if any(strcmp(evalin('base','who'),imageMatrixName)) == 0
    close(alignWaitBar);
    errordlg('Channel track not found. Please run program again for appropriate channels', 'Channel selection error');
    error('Channel track not found. Please run program again for appropriate channels');
end
%find data set
data = evalin('base',imageMatrixName);
%save data set
tracksVector{1,2} = data;
tracksVector{2,2} = imageMatrixName;
countChannels = countChannels + 1;
end
if useChannel3 == 1
Channel3Handle = findobj('tag','Channel3');
imageMatrixName = get(Channel3Handle,'string');
if any(strcmp(evalin('base','who'),imageMatrixName)) == 0
    close(alignWaitBar);
    errordlg('Channel not found. Please match channel name to an existing file or variable', 'Channel selection error');
    error('Channel not found. Please match channel name to an existing file or variable');
end
imageMatrixName = ['finalTracks_' imageMatrixName];
if any(strcmp(evalin('base','who'),imageMatrixName)) == 0
    close(alignWaitBar);
    errordlg('Channel track not found. Please run program again for appropriate channels', 'Channel selection error');
    error('Channel track not found. Please run program again for appropriate channels');
end
%find data set
data = evalin('base',imageMatrixName);
%save data set
tracksVector{1,3} = data;
tracksVector{2,3} = imageMatrixName;
countChannels = countChannels + 1;
end

if countChannels ~= 2
    close(alignWaitBar);
    errordlg('Please select two channels for channel alignment', 'Channel alignment error');
    error('Please select two channels for channel alignment');
end

channels2use = ~cellfun(@isempty,tracksVector);
use = find(channels2use(1,:));
alignedTracksNames = tracksVector(2,use);
assignin('base', 'alignedTracksNames', alignedTracksNames);
%find which channel tracks to use

%two blocks to create targets and sources for creating track links between
%channels
channel1Matrix = cell2mat(tracksVector(1,use(1)));
logicalIndices = channel1Matrix~=0;
elements = numel(find(logicalIndices~=0));
matrix = zeros(elements/4, 4);
startRow = 1;
for i = 1:size(logicalIndices,3)
    rows =  numel(find(logicalIndices(:,1,i) == true));
    matrix(startRow:startRow+rows-1,1:4) = channel1Matrix(1:rows,1:4,i);
    startRow = startRow + rows;
end
channel1Matrix = sortrows(matrix,4);

channel2Matrix = cell2mat(tracksVector(1,use(2)));
logicalIndices = channel2Matrix~=0;
elements = numel(find(logicalIndices~=0));
matrix = zeros(elements/4, 4);
startRow = 1;
for i = 1:size(logicalIndices,3)
    rows =  numel(find(logicalIndices(:,1,i) == true));
    matrix(startRow:startRow+rows-1,1:4) = channel2Matrix(1:rows,1:4,i);
    startRow = startRow + rows;
end
channel2Matrix = sortrows(matrix,4);
%channel#matrix has all points for all tracks of one channel sorted
%according to timepoint

timePoints = max(max(channel1Matrix(:,4)),max(channel2Matrix(:,4)));
%finds number of time points in data set

targets = repmat(-1,size(channel1Matrix,1),1);
row = 1;
for timePoint = 1:timePoints
    matrix1 = channel1Matrix((channel1Matrix(:,4) == timePoint),1:3);
    matrix2 = channel2Matrix((channel2Matrix(:,4) == timePoint),1:3);
    target_indices = hungarianlinker(matrix1, matrix2, maxLinkingDistance, pixelSize, stepSize);
    target_indices = target_indices';
    %find closest point in other channel
    endRow = (size(target_indices,1));
    %This next part generates the actual positions of the targets
for j = 1:size(target_indices,1)
    if target_indices(j)>0
        coordinate = [matrix2(target_indices(j),1:3) timePoint];
        realPosition = find(ismember(channel2Matrix, coordinate, 'rows') ~= 0);
        target_indices(j) = realPosition;
    else
        continue
    end
end
    %keep original target position values
    targets(row:endRow+row-1) = target_indices;
    row = row+endRow;
end

saverMatrix = compMiddle(channel1Matrix, channel2Matrix, targets);
midMatrix = saverMatrix(:,1:4);
place = (1:size(midMatrix,1))';
midMatrix = [place midMatrix];
%creates a matrix with all mid point coordinates and their original
%placement (original place,Y,X,Z,t)
sortedMidMatrix = sortrows(midMatrix,5);
%sort midMatrix according to time point

[frequentValue, frequency] = mode(sortedMidMatrix(:,5));
linksSize = max(frequency(:));
%finds max number of connectible particles in all time points
links = repmat(-1, timePoints, linksSize);
%creates a matrix that saves the links between two particles
coordinateSaver = repmat(-1, [frequency, 3, timePoints]);
%creates a matrix that saves coordinates for each time point

for i = 1:timePoints-1
matrix1Index = find(sortedMidMatrix(:,5) == i);
matrix2Index = find(sortedMidMatrix(:,5) == i+1);
timePoint1Matrix = sortedMidMatrix(matrix1Index, 2:4);
timePoint2Matrix = sortedMidMatrix(matrix2Index, 2:4);

 timePoint1Elements = numel(timePoint1Matrix);
 timePoint2Elements = numel(timePoint2Matrix);
 coordinateSaver(1:(timePoint1Elements/3), :, i) = timePoint1Matrix;
 coordinateSaver(1:(timePoint2Elements/3), :, i+1) = timePoint2Matrix;
 %At every loop iteration the two metrices are cleared. Here they are 
 %saved in the coordinates saving metrix that was declared outside the loop. 
 %The number of coordinates is third the number of elements as each 
 %coordinate has three values.
 
 target_indices = hungarianlinker(timePoint1Matrix, timePoint2Matrix, maxLinkingDistance, pixelSize, stepSize);
 %The hungarian algorithm receives two metrices with connectible
 %coordinates and a max linking distance value. It uses these inputs to
 %find the global minima for particle connections. The output is specified
 %in the algorithm itself.
 
 targetElements = numel(target_indices);
 links(i, 1:targetElements) = target_indices;
 %At every iteration the target_indices metrix is cleared. Here it is saved 
 %in the links saving metrix that was declared outside the loop.
 
end

alignedTracks  = linksToTracks( coordinateSaver, links );

maxGap = 3;
finalAlignedTracks = gapCloser( alignedTracks, maxGap, (maxLinkingDistance+maxGap), pixelSize, stepSize, minimumTrackLength );
originalPositions = findPosition(finalAlignedTracks, sortedMidMatrix);

originalCoordinates = zeros(size(finalAlignedTracks,1), 8, size(finalAlignedTracks,3));
positions = find(originalPositions~=0);

for i = 1:numel(positions);
    index = positions(i);
    originalPosition = originalPositions(index);
    coordinates = saverMatrix(originalPosition, 5:12);
    [row, column, track] = ind2sub(size(originalPositions), index);
    originalCoordinates(row, :, track) = coordinates;
end

finalAlignedTracks = [finalAlignedTracks originalCoordinates];
assignin('base', 'finalAlignedTracks', finalAlignedTracks);
%save final aligned tracks globally

%set tracks to table data to allow user to show or hide tracks
dims = ndims(finalAlignedTracks);
tableHandle = findobj('tag', 'tracksTable4');
if dims == 2
    set (tableHandle, 'data', {1 true});
else
    numberTracks = size(finalAlignedTracks);
    iVector = 1:numberTracks(3);
    iVector = iVector';
    trueVector = true(numberTracks(3),1);
    data = [num2cell(iVector) num2cell(trueVector)];
    set (tableHandle, 'data', data);
end


figureHandle = findobj('tag', 'trackAlignmentFigure');
set(figureHandle, 'Visible', 'on');

boxHandle = findobj('tag', 'alignedAnalysis');
set(boxHandle, 'Enable', 'on');

close(alignWaitBar);
end

function [ midPoints ] = compMiddle(sourceMatrix, targetMatrix, linkerMatrix)
%This function receives a source, target and linker matrix and calculates
%the center of each two linked points. For target_indices each row
%corrosponds to source row amd the number specifies target row (at least 
%I think so).

indices = find(linkerMatrix ~= -1);
rows = size(indices,1);
midPoints = zeros(rows,12);
%midpoint(YXZt), source coordinate (YXZt), target coordinate (YXZt)

for i = 1:rows
    sourceIndex = indices(i);
    source = sourceMatrix(sourceIndex, 1:4);
    targetIndex = linkerMatrix(sourceIndex);
    target = targetMatrix(targetIndex, 1:4);
    midPoints(i,1:4) = (source + target)./2;
    midPoints(i,5:8) = source;
    midPoints(i,9:12) = target;
end


end

function [ positionMatrix ] = findPosition(finalTracks, fullMatrix)
%This function receives a coordinate matrix and finds their original
%position

positionMatrix = zeros((size(finalTracks,1)),1,size(finalTracks,3));

for currentTrack = 1:(size(finalTracks,3))
     coordinates = finalTracks(:, :, currentTrack);
     isMember = ismember(fullMatrix(:,2:5),coordinates,'rows');
     indices = fullMatrix(isMember,1);
     positionMatrix(1:numel(indices),:,currentTrack) = indices;
end

end

function [] = showAll (calling_button, eventData, table)
    %callback
    %either shows all tracks or shows none
    tableHandle = findobj('tag', table);
    tableData = get(tableHandle, 'Data');
    rows = size(tableData);
    rows = rows(1);
    if get(calling_button, 'value') == 1
        tableData(:,2) = num2cell(true(rows, 1));
    else
        tableData(:,2) = num2cell(false(rows, 1));
    end
    set(tableHandle, 'Data', tableData);
end