function [output] = templateBuilder3DTool(varargin)
% templateBuilder3DTool M-file for templateBuilder3DTool.fig
%
% Description:
%
%
% Author: Johann Strasser
% Date: 070312

error(nargchk(0, inf, nargin));

if mod(length(varargin), 2) ~= 0
    % input did not come in parameter/value pairs
    error('Arguments have to come in parameter/value pairs.');
end

% Initialise default application data structure to show all the fields.
ad.guiHandles = [];
ad.projectDirectory = '';
ad.aam = [];
ad.output = [];

% Process input arguments
for i=1:2:length(varargin)
    switch lower(varargin{i})
        %         case 'aam'
        %             ad.aam = varargin{i + 1};
        case 'projectdirectory'
            ad.projectDirectory = varargin{i + 1};
        otherwise
            error(['Unknown parameter name passed to templateBuilder3DTool. ', ...
                'Parameter name: ',  varargin{i}]);
    end
end

% Check if required application data was set, if not set default values
% and issue warning
if isempty(ad.projectDirectory)
    ad.projectDirectory = '.';
    warning(sprintf(['TemplateBuilder3DTool: Parameter ''projectDirectory'' not provided. \n', ...
        'Using current directory (.) as default']));
end

% if isempty(ad.aam)
%    ad.aam = AAM;
%    warning(sprintf(['TemplateBuilder3DTool: Parameter ''aam'' not provided. \n', ...
%        'Using empty AAM object as default']));
% end

% Open figure file
fig = openfig(mfilename, 'new');

% Get handles
handles = guihandles(fig);
ad.guiHandles = handles;

% Set the callback functions for all the gui elements
set(fig, 'closeRequestFcn', {@customCloseRequestFcn});
set(ad.guiHandles.editTemplatePshBtn, 'callback', {@editTemplateCallback});
set(ad.guiHandles.newTemplatePshBtn, 'callback', {@newTemplateCallback});
set(ad.guiHandles.exportTemplateVerticesPshBtn, 'callback', {@exportTemplateVerticesCallback});
set(ad.guiHandles.viewVolumePshBtn, 'callback', {@viewVolumeCallback});

% Setup GUI elements
setupGUIElements(ad);

% Save the application data
setappdata(0, 'templateBuilder3DToolData', ad);

% Move the gui and show it
movegui(fig, 'center');
set(fig, 'visible', 'on');

% Make sure this is a modal dialog
try
    % Wait until resume is called or figure is deleted
    uiwait(fig);
catch
    if ishandle(fig)
        delete(fig)
    end
end

% display('M-file execution resumed.');

% Set return arguments here. The fact that the figure has been destroyed
% is irrelevant as the appdata is associated with the root handle graphics
% object
ad = getappdata(0, 'templateBuilder3DToolData');
output = ad;
rmappdata(0, 'templateBuilder3DToolData');

function viewVolumeCallback(cbo, event)

% display('viewVolumeCallback executing.');

ad = getappdata(0, 'templateBuilder3DToolData');

% Set the figure invisible, note that the event is queued as uiwait was
% called on the figure handle. Hence call drawnow expose to process events
% causing graphics object updates only.
set(ad.guiHandles.templateBuilder3DToolFig, 'visible', 'off');
drawnow expose;
viewVolume('path', fullfile(ad.projectDirectory, 'Cropped'));
set(ad.guiHandles.templateBuilder3DToolFig, 'visible', 'on');
setappdata(0, 'templateBuilder3DToolData', ad);

function newTemplateCallback(cbo, event)

display('newTemplateCallback executing.');

ad = getappdata(0, 'templateBuilder3DToolData');
ad = newTemplate(ad);
setappdata(0, 'templateBuilder3DToolData', ad);

function editTemplateCallback(cbo, event)

display('editTemplateCallback executing.');

ad = getappdata(0, 'templateBuilder3DToolData');
ad = editTemplate(ad);
setappdata(0, 'templateBuilder3DToolData', ad);


function exportTemplateVerticesCallback(cbo, event)

display('exportTemplateVerticesCallback executing.');

ad = getappdata(0, 'templateBuilder3DToolData');
ad = exportTemplateVertices(ad);
setappdata(0, 'templateBuilder3DToolData', ad);


function customCloseRequestFcn(cbo, event)

% display('customCloseRequestFcn executing.');

% The following is the code in the default close request function
% (closereq)
if isempty(gcbf)
    if length(dbstack) == 1
        warning('MATLAB:customCloseRequestFcn', ...
            'Calling customCloseRequestFcn from the command line is now obsolete, use close instead');
    end
    close force
else
    delete(gcbf);
end

function setupGUIElements(ad)

% Populate the list box
% volumeNames = get(ad.aam, 'imageFilenames');
% set(ad.guiHandles.volumesListBox, 'String', volumeNames);

function ad = editTemplate(ad)

imreadSupportedFileTypes = getImreadSupportedTypesFilterSpec();

[fileName, path] = uigetfile({'*.temp_dat'}, ...
    'Choose template', fullfile(ad.projectDirectory, ...
    'Templates', filesep));

if isequal(path, 0) && isequal(fileName, 0)
    disp('User selected cancel.');
    userCancel = 1;
else
    [fullVolumePath, name, extension, versn] = fileparts(fullfile(path, fileName));
    templateExtensionsName = [name, '.temp_ex_dat'];

    % Load template
    temp = load(fullfile(path, fileName), '-mat');
    PMTemplate = temp.PMTemplate;
    
    %load template hack Paul
    %PMTemplate = temp.x5000whgfp3_x17p77y17p77z17p77pop
    %end load template hack Paul

    % The following commented out code was used to convert templates made with a changed
    % version of class PointModelTemplate

    %     % If the PMTemplate from disk is an instance of the PMTemplate class
    %     % with additional fields, then convert.
    %     if ~isa(PMTemplate, 'PointModelTemplate')
    %         template = PointModelTemplate();
    %         template = set(template, 'imageName', PMTemplate.imageName);
    %         template = set(template, 'imageDir', PMTemplate.imageDir);
    %         template = set(template, 'name', PMTemplate.name);
    %         template = set(template, 'direc', PMTemplate.direc);
    %         template = set(template, 'pts', PMTemplate.pts);
    %         template = set(template, 'primaries', PMTemplate.primaries);
    %         template = set(template, 'loops', PMTemplate.loops);
    %
    %         % Initialise extension data
    %         lineIndices = PMTemplate.lineIndices;
    %         vertexColours = PMTemplate.vertexColours;
    %         lineColours = PMTemplate.lineColours;
    %
    %         PMTemplate = template;
    %
    %         % Initialise non-extension data
    %         volumeName = get(PMTemplate, 'imageName');
    %         imageDir = get(PMTemplate, 'imageDir');
    %         templateName = get(PMTemplate, 'name');
    %         directory = get(PMTemplate, 'direc');
    %         vertices = get(PMTemplate, 'pts');
    %         primaries = get(PMTemplate, 'primaries');
    %         loops = get(PMTemplate, 'loops');
    %     else
    % Initialise non-extension data
    volumeName = get(PMTemplate, 'imageName');
    imageDir = get(PMTemplate, 'imageDir');
    templateName = get(PMTemplate, 'name');
    directory = get(PMTemplate, 'direc');
    vertices = get(PMTemplate, 'pts');
    primaries = get(PMTemplate, 'primaries');
    loops = get(PMTemplate, 'loops');

    % The following data are extensions to objects of class
    % PointModelTemplate and are stored separately as changing the class
    % file would make the toolbox incompatible with older projects.
    temp = load(fullfile(path, templateExtensionsName), '-mat');
    pmtExtensions = temp.pmtExtensions;
    lineIndices = pmtExtensions.lineIndices;
    vertexColours = pmtExtensions.vertexColours;
    lineColours = pmtExtensions.lineColours;
    
     % Code to change the vertex colours from blue to something else
%     colour = [0, 0.5, 1];
%     intermediatsColumns = vertexColours(1, :) == 0 & ...
%         vertexColours(2, :) == 0 & ...
%         vertexColours(3, :) == 1;
%     vertexColours(1, intermediatsColumns) = colour(1);
%     vertexColours(2, intermediatsColumns) = colour(2);
%     vertexColours(3, intermediatsColumns) = colour(3);
    %         v = vertexColours;
    %         red = v(1, :);
    %         green = v(2, :);
    %         blue = v(3, :);
    %         redNew = ones(size(red));
    %         greenNew = redNew;
    %         greenNew(green == 1) = 165 / 255;
    %         v(1, :) = redNew;
    %         v(2, :) = greenNew;
    %         v(3, :) = 0;
    %         vertexColours = v;
    %     end


    % Load volume
    fullVolumePath = fullfile(ad.projectDirectory, imageDir, volumeName);
    volume = loadVolumeFromSlicesDir('fullVolumePath', fullVolumePath);

    % Make RGB for hapticToolNative if required
    if size(volume, 4) == 1
        volume = repmat(volume, [1, 1, 1, 3]);
    end

    % Initialise the transform structure for the template
    t = getHTTransform();
    t.name = 'tVolumeTransform';

    % Initialise volume structure for the template
    v = getHTVolume();
    v.name = 'tVolume';
    v.parentTransform = 'tVolumeTransform';
    v.volume = volume;
    v.planes = 256;

    % Begin change js 081114
    % Begin original code

    % Vertices are in unit voxels (1 micron cubed), convert to metres for
    % H3DApi then scale by1000 for sensible display size
    %( vertices / 1000000 * 1000 = vertices % / 1000);
    %     vertices = vertices / 1000; % ChangeHere for real space

    % For project template conversion: !!!!!!!!!!!!!!!!!
    % Put them into real space
    % Vertices are in unit voxels (voxels are 1 micron cubed).
    % Convert to real space
    %     voxelSize = getVoxelSizeFromString('string', volumeName);
    %     scaleMat = diag(voxelSize);
    %     vertices = scaleMat * vertices;

    % Vertices are in real space so convert them to unit voxel space (1
    % micron cubed) then convert to metres for H3DApi and scale by 1000 for
    % sensible display size
    % (diag(voxelSize.^-1) * vertices ./ 1000000 .* 1000 =
    % diag(voxelSize.^-1) * vertices ./ 1000
    
    % Used code before
%     voxelSize = getVoxelSizeFromString('string', volumeName);
%     scaleMat = diag(voxelSize.^-1);
%     vertices = scaleMat * vertices ./ 1000;
    
    % End original code
    % Begin replacement code
    
    % See function newTemplate for an explanation of the scaling and the
    % rationale behind it
    voxelSize = getVoxelSizeFromString('string', volumeName);
    minVoxelDim = min(voxelSize);
    
    % Set the display voxel size
    displayVoxelSize = voxelSize ./ minVoxelDim ./ 1000;
    v.voxelSize = displayVoxelSize;
    
    % Set the display size for vertices
    vertices = vertices ./ minVoxelDim ./ 1000;
    
    % End replacement code
    % End change js 081114

    ils = getHTIndexedLineSet();
    ils.name = 'tIndexedLineSet';
    ils.parentTransform = 'tVolumeTransform';
    ils.vertices = vertices;
    ils.lineIndices = lineIndices;
    ils.interactionMode = 'edit';
    ils.vertexColours = vertexColours;
    ils.lineColours = lineColours;
    ils.vertexStiffness = 0.5;
    ils.lineStiffness = 0.5;

    %     cp = getHTClipPlane();
    %     cp.name = 'clipPlane1';
    %     cp.parentTransform = 'tVolumeTransform';
    %     cp.point = [0, -0.05, 0];
    %     cp.normal = [0, 1, 0];
    %     cp.interactionMode = 'edit';

    % Deprecated
    %     noOfClipPlanes = 4;
    %     cpXOffset = 0.03;
    %     cpXBasePos = 0 - (noOfClipPlanes * cpXOffset / 2);
    %     cpYBasePos = -size(volume, 2) / 1000 / 2;
    %
    %     for i = 1:noOfClipPlanes
    %
    %         cp = getHTClipPlane();
    %         cp.name = ['clipPlane', num2str(i)];
    %         cp.parentTransform = 'tVolumeTransform';
    %         cp.point = [cpXBasePos + cpXOffset * (i - 1), cpYBasePos, 0];
    %         cp.normal = [0, 1, 0];
    %         cp.interactionMode = 'edit';
    %         clipPlaneStructures(i) = cp;
    %     end

    % Initalise clipping planes
    noOfClipPlanes = 8;
    % cpXOffset = 0.03;

    if noOfClipPlanes ~= 1
        cpXOffset = size(volume, 1) / (noOfClipPlanes - 1);
        cpXBasePos = 0 - ((noOfClipPlanes - 1) * cpXOffset / 2);
        cpYBasePos = -size(volume, 2) / 2;
    else
        cpXOffset = -size(volume, 1) / 2;
        cpXBasePos = -size(volume, 2) / 2;
        %     cpYBasePos = -size(currentVolume.data, 2) / 2;
        cpYBasePos = -size(volume, 2) / 2;
    end

    for i = 1:noOfClipPlanes

        cp = getHTClipPlane();
        cp.name = ['clipPlane', num2str(i)];
        cp.parentTransform = 'tVolumeTransform';
        cp.point = [cpXBasePos + cpXOffset * (i - 1), cpYBasePos, 0] / 1000;

        % Apply voxel size to get into real space then apply display scale
        %     cp.point = cp.point .* voxelSize .* displayVoxelScale;
        cp.normal = [0, 1, 0];
        cp.interactionMode = 'edit';
        clipPlaneStructures(i) = cp;
    end
    
    if noOfClipPlanes == 0
        clipPlaneStructures = [];
    end

    % Set the figure invisible, note that the event is queued as uiwait was
    % called on the figure handle. Hence call drawnow expose to process events
    % causing graphics object updates only.
    set(ad.guiHandles.templateBuilder3DToolFig, 'visible', 'off');
    drawnow expose;

    % Run hapticTool
    %     [tOut, vOut, pOut, ilsOut, cpOut] =...
    %             hapticTool(t, v, [], ils, cp);

    [tOut, vOut, pOut, ilsOut, cpOut] =...
        hapticTool(t, v, [], ils, clipPlaneStructures);

    set(ad.guiHandles.templateBuilder3DToolFig, 'visible', 'on');

    % Get the volume directory name
    [starts, ends, extents, matches, tokens, names] = ...
        regexpi(fullVolumePath, ['[^', '\', filesep, ']+']);

    % Begin change js 081114
    % Begin original code

    % Vertices are returned in unit metre voxels scaled up by 1000.
    % Convert to microns, then scale down
    % (vertices * 1000000 / 1000 = vertices * 1000).
    % vertices = ilsOut.vertices * 1000;   % ChangeHere for real space

    % Vertices are in unit metre voxels scaled up by 1000.
    % Convert to real space, then to microns then scale down.
    % (diag(voxelSize) * vertices .* 1000000 ./ 1000 =
    % diag(voxelSize) * vertices .* 1000
    
    % Used code before
%     voxelSize = getVoxelSizeFromString('string', volumeName);
%     scaleMat = diag(voxelSize);
%     vertices = scaleMat * ilsOut.vertices .* 1000;
    
    % End original code
    % Begin replacement code
    
    vertices = ilsOut.vertices .* minVoxelDim .* 1000;
    
    % End replacement code
    % End change js 081114

    primaries = [];
    loops = [];

    % Scan through the line indices field and convert to primaries/loops
    % convention. Conversion rules are as follows (simple by choice of
    % native data structure):
    % Polyline end vertices are primaries
    % Each polyline segment is a loop

    indices = ilsOut.lineIndices;
    polylines = {};
    polylineCount = 0;
    currentPolylineIndices = [];

    for i = 2:length(indices)
        if indices(i) ~= -1
            % Beware that the indices are C-style subscripts, hence increment
            currentPolylineIndices = [currentPolylineIndices; indices(i)];
        else
            % We have reached the end of a polyline
            polylineCount = polylineCount + 1;
            polylines{polylineCount} = currentPolylineIndices;
            currentPolylineIndices = [];
        end
    end

    for i = 1:length(polylines)
        currentPolylineIndices = polylines{i};

        if isempty(find(ismember(primaries, currentPolylineIndices(1) + 1)))
            primaries = [primaries, currentPolylineIndices(1) + 1];
        end

        if isempty(find(ismember(primaries, currentPolylineIndices(end) + 1)))
            primaries = [primaries, currentPolylineIndices(end) + 1];
        end

        for j = 1:length(currentPolylineIndices) - 1
            loops = [loops; ...
                [currentPolylineIndices(j), currentPolylineIndices(j + 1)]];
        end
    end


    % Set member variables of PMTemplate object
    loopsCell = {};
    loopsCell{1, 1} = loops + 1;
    PMTemplate = PointModelTemplate;
    PMTemplate = set(PMTemplate, 'imageName', volumeName);
    PMTemplate = set(PMTemplate, 'imageDir', 'Templates');
    PMTemplate = set(PMTemplate, 'name', templateName);
    PMTemplate = set(PMTemplate, 'direc', 'Templates');
    PMTemplate = set(PMTemplate, 'pts', vertices);
    PMTemplate = set(PMTemplate, 'primaries', primaries);
    PMTemplate = set(PMTemplate, 'loops', loopsCell);

    % Store the PointModelTemplate extensions
    pmtExtensions.lineIndices = ilsOut.lineIndices;
    pmtExtensions.vertexColours = ilsOut.vertexColours;
    pmtExtensions.lineColours = ilsOut.lineColours;
    pmtExtensions.volumeSize = size(v.volume);
    
       % Code for fixing the line colour inconsistencies
%     lineColours = ilsOut.vertexColours;
%     lineColours(:, :) = 0;
%     lineColours(1, :) = 1;
%     pmtExtensions.lineColours = lineColours;

    % Name and path to be requested from and/or confirmed by the user
    defaultTemplateName = templateName;

    % Save the template if desired
    button = questdlg('Save template?', 'Save', 'Yes', 'Save as', 'No', 'Yes');

    if isequal(button, 'Save as')
        [templateName, templatePath] = uiputfile('.temp_dat', ...
            'Save template', fullfile(ad.projectDirectory, ...
            'Templates', defaultTemplateName));

        if isequal(templatePath, 0) && isequal(templateName, 0)
            disp('User selected cancel.');
            userCancel = 1;
        else
            % Save template
            save(fullfile(templatePath, templateName), 'PMTemplate');

            % Save template extensions
            [temPath, tName, extension, versn] = ...
                fileparts(fullfile(templatePath, templateName));
            templateExtensionsName = [tName, '.temp_ex_dat'];
            save(fullfile(templatePath, templateExtensionsName), 'pmtExtensions');
        end
    elseif isequal(button, 'Yes')
        templatePath = fullfile(ad.projectDirectory, 'Templates');

        % Save template
        save(fullfile(templatePath, templateName), 'PMTemplate');

        % Save template extensions
        save(fullfile(templatePath, templateExtensionsName), 'pmtExtensions');
    end
end


function ad = exportTemplateVertices(ad)

imreadSupportedFileTypes = getImreadSupportedTypesFilterSpec();

[fileName, path] = uigetfile({'*.temp_dat'}, ...
    'Choose template', fullfile(ad.projectDirectory, ...
    'Templates', filesep));

if isequal(path, 0) && isequal(fileName, 0)
    disp('User selected cancel.');
    userCancel = 1;
else
    [fullVolumePath, name, extension, versn] = fileparts(fullfile(path, fileName));
    templateExtensionsName = [name, '.temp_ex_dat'];

    % Load template
    temp = load(fullfile(path, fileName), '-mat');
    PMTemplate = temp.PMTemplate;
   
    % Initialise non-extension data
    volumeName = get(PMTemplate, 'imageName');
    imageDir = get(PMTemplate, 'imageDir');
    templateName = get(PMTemplate, 'name');
    directory = get(PMTemplate, 'direc');
    vertices = get(PMTemplate, 'pts');
    primaries = get(PMTemplate, 'primaries');
    loops = get(PMTemplate, 'loops');

    % The following data are extensions to objects of class
    % PointModelTemplate and are stored separately as changing the class
    % file would make the toolbox incompatible with older projects.
    temp = load(fullfile(path, templateExtensionsName), '-mat');
    pmtExtensions = temp.pmtExtensions;
    lineIndices = pmtExtensions.lineIndices;
    vertexColours = pmtExtensions.vertexColours;
    lineColours = pmtExtensions.lineColours;

    % Save the vertices in a mat and a text file
    [tName, remain] = strtok(templateName, '.');
    verticesExtension = '.vertices';
    
    [matFileName, filePath] = uiputfile(verticesExtension, ...
        'Save template vertices', fullfile(ad.projectDirectory, ...
        'Templates', [tName, verticesExtension]));

    if isequal(filePath, 0) && isequal(matFileName, 0)
        disp('User selected cancel.');
        userCancel = 1;
    else
        vertices = vertices';
        save(fullfile(filePath, matFileName), 'vertices');
               
        [vFileName, remain] = strtok(matFileName, '.');
        dlmwrite(fullfile(filePath, [vFileName, '.csv']), vertices, ',');       
    end
end


function ad = newTemplate(ad)

imreadSupportedFileTypes = getImreadSupportedTypesFilterSpec();

[fileName, path] = uigetfile(imreadSupportedFileTypes, ...
    'Choose template volume', fullfile(ad.projectDirectory, ...
    'Templates', filesep));

if isequal(path, 0) && isequal(fileName, 0)
    disp('User selected cancel.');
    userCancel = 1;
else
    [fullVolumePath, name, extension, versn] = fileparts(fullfile(path, fileName));
    [p, volumeDirectory, e, v] = fileparts(fullVolumePath);

    % Load volume
    volume = loadVolumeFromSlicesDirTyped('fullVolumePath', fullVolumePath, ...
        'extension', extension);

    % Make RGB for hapticToolNative if required
    if size(volume, 4) == 1
        volume = repmat(volume, [1, 1, 1, 3]);
    end

    % Initialise the transform structure for the template
    t = getHTTransform();
    t.name = 'tVolumeTransform';

    % Initialise volume structure for the template
    v = getHTVolume();
    v.name = 'tVolume';
    v.parentTransform = 'tVolumeTransform';
    v.volume = volume;
    v.planes = 256;
   
    % Begin change js 081114
    % Begin original code
    
     % Get voxel size from volume name and scale down as voxel sizes in H3DApi
    % are in metres
    %     v.voxelSize = getVoxelSizeFromString('string', volumeDirectory) / 1000000;
%     voxelSize = getVoxelSizeFromString('string', volumeDirectory);
    
    % End original code
    % Begin replacement code
    
    % Set the volume's voxel size to reflect non-cubic voxel sizes
    % But we don't want to set the volume's voxel size to the raw voxel
    % size data, as this would mean that the display size would be the
    % real-space size of the volume which entails massively differing
    % display sizes of volumes depending on the voxel sizes. Instead we aim
    % to display the volume at image space size to get the same display
    % size for all volumes. Obviously for non-cubic voxel data that is not
    % possible, so we normalise by the smallest voxel dimension. So the
    % smallest voxel has a size of 1 and the other ones are factors of it.
    % Of course, this means we have to be careful when scaling the point
    % sets to real space, as the point sets we get are not in image space
    % anymore, but in the smallest-voxel-normalised space. We normalise by
    % the smallest voxel size as we expect data from confocal microscopy,
    % where the z size is always bigger than the x and y size due to the
    % lower resolution in the z axis. Normalising by the smallest ensures
    % we maintain the high resolution in the xy plane. If we normalised by
    % the largest value we would unnecessarily scale down the high
    % resolution in the xy plane. 
    % Note that H3DApi expects voxel sizes in metres, hence we divide by
    % 1000000 and because we want a reasonable display size we scale up by
    % 1000, resulting in a division by 1000 below
    voxelSize = getVoxelSizeFromString('string', volumeDirectory);
    minVoxelDim = min(voxelSize);
    displayVoxelSize = voxelSize ./ minVoxelDim ./ 1000;
    v.voxelSize = displayVoxelSize;
    
    % End replacement code
    % End change js 081114
    
    ils = getHTIndexedLineSet();
    ils.name = 'tIndexedLineSet';
    ils.parentTransform = 'tVolumeTransform';
    ils.vertices = [];
    ils.lineIndices = [];
    ils.interactionMode = 'edit';
    ils.vertexColours = []';
    ils.vertexStiffness = 0.5;
    ils.lineStiffness = 0.5;

    %     cp = getHTClipPlane();
    %     cp.name = 'clipPlane1';
    %     cp.parentTransform = 'tVolumeTransform';
    %     cp.point = [0, -0.05, 0];
    %     cp.normal = [0, 1, 0];
    %     cp.interactionMode = 'edit';

 % Initalise clipping planes
    noOfClipPlanes = 8;
    % cpXOffset = 0.03;

    if noOfClipPlanes ~= 1
        cpXOffset = size(volume, 1) / (noOfClipPlanes - 1);
        cpXBasePos = 0 - ((noOfClipPlanes - 1) * cpXOffset / 2);
        cpYBasePos = -size(volume, 2) / 2;
    else
        cpXOffset = -size(volume, 1) / 2;
        cpXBasePos = -size(volume, 2) / 2;
        %     cpYBasePos = -size(currentVolume.data, 2) / 2;
        cpYBasePos = -size(volume, 2) / 2;
    end

    for i = 1:noOfClipPlanes

        cp = getHTClipPlane();
        cp.name = ['clipPlane', num2str(i)];
        cp.parentTransform = 'tVolumeTransform';
        cp.point = [cpXBasePos + cpXOffset * (i - 1), cpYBasePos, 0] / 1000;

        % Apply voxel size to get into real space then apply display scale
        %     cp.point = cp.point .* voxelSize .* displayVoxelScale;
        cp.normal = [0, 1, 0];
        cp.interactionMode = 'edit';
        clipPlaneStructures(i) = cp;
    end
    
    if noOfClipPlanes == 0
        clipPlaneStructures = [];
    end
% Set the figure invisible, note that the event is queued as uiwait was
% called on the figure handle. Hence call drawnow expose to process events
    % causing graphics object updates only.
    set(ad.guiHandles.templateBuilder3DToolFig, 'visible', 'off');
    drawnow expose;

    % Run hapticTool
    %     [tOut, vOut, pOut, ilsOut, cpOut] =...
    %             hapticTool(t, v, [], ils, cp);

    [tOut, vOut, pOut, ilsOut, cpOut] =...
        hapticTool(t, v, [], ils, clipPlaneStructures);

    set(ad.guiHandles.templateBuilder3DToolFig, 'visible', 'on');

    % Get the volume directory name
    [starts, ends, extents, matches, tokens, names] = ...
        regexpi(fullVolumePath, ['[^', '\', filesep, ']+']);

    volumeName = char(matches(length(matches)));

    % Name and path to be requested from and/or confirmed by the user
    defaultTemplateName = [volumeName, '.temp_dat'];

    % Save the template if desired
    button = questdlg('Save new template?', 'Save', 'Yes', 'No', 'Yes');

    if isequal(button, 'Yes')
        [templateName, templatePath] = uiputfile('.temp_dat', ...
            'Save template', fullfile(ad.projectDirectory, ...
            'Templates', defaultTemplateName));

        if isequal(templatePath, 0) && isequal(templateName, 0)
            disp('User selected cancel.');
            userCancel = 1;
        else
            % Begin change js 081114
            % Begin original code
        
            % Vertices are returned in unit metre voxels scaled up by 1000.
            % Convert to microns, then scale down
            % (vertices * 1000000 / 1000 = vertices * 1000).
            % vertices = ilsOut.vertices * 1000;  % ChangeHere for real space

            % Convert from unit voxel space to real space
            % voxelSize = getVoxelSizeFromString('string', volumeDirectory);
            % scaleMat = diag(voxelSize);
            % vertices = scaleMat * vertices;

            % Vertices are in unit metre voxels scaled up by 1000.
            % Convert to real space, then to microns then scale down.
            % (diag(voxelSize.^1) * vertices .* 1000000 ./ 1000 =
            % diag(voxelSize.^1) * vertices .* 1000
%             voxelSize = getVoxelSizeFromString('string', volumeDirectory);
            
%             scaleMat = diag(voxelSize);
%             vertices = scaleMat * ilsOut.vertices .* 1000;

            % End original code
            % Begin replacement code
            
            % Notes: testing with whirlpool_stack_x0p116y0p116z1p0_cropped
            % Image space volume dimensions: 259 * 268 * 44
            % Real space volume dimensions (x * y * z): 30.044 * 31.088 * 44
            
            vertices = ilsOut.vertices .* minVoxelDim .* 1000;
            
            % End replacement code
            % End change js 081114
            
            primaries = [];
            loops = [];

            % Scan through the line indices field and convert to primaries/loops
            % convention. Conversion rules are as follows (simple by choice of
            % native data structure):
            % Polyline end vertices are primaries
            % Each polyline segment is a loop
            indices = ilsOut.lineIndices;
            polylines = {};
            polylineCount = 0;
            currentPolylineIndices = [];

            for i = 2:length(indices)
                if indices(i) ~= -1
                    % Beware that the indices are C-style subscripts, hence increment
                    currentPolylineIndices = [currentPolylineIndices; indices(i)];
                else
                    % We have reached the end of a polyline
                    polylineCount = polylineCount + 1;
                    polylines{polylineCount} = currentPolylineIndices;
                    currentPolylineIndices = [];
                end
            end

            for i = 1:length(polylines)
                currentPolylineIndices = polylines{i};

                if isempty(find(ismember(primaries, currentPolylineIndices(1))))
                    primaries = [primaries, currentPolylineIndices(1) + 1];
                end

                if isempty(find(ismember(primaries, currentPolylineIndices(end))))
                    primaries = [primaries, currentPolylineIndices(end) + 1];
                end

                for j = 1:length(currentPolylineIndices) - 1
                    loops = [loops; ...
                        [currentPolylineIndices(j), currentPolylineIndices(j + 1)]];
                end
            end

            % Old code:
            %             for i = 1:length(ilsOut.lineIndices)
            %                 if ilsOut.lineIndices(i) == -1
            %                     lastPrimary = [];
            %                     nextPrimary = [];
            %                     if i ~= length(ilsOut.lineIndices)
            %                         nextPrimary = ilsOut.lineIndices(i + 1);
            %                     end
            %
            %                     if i ~= 1
            %                         lastPrimary = ilsOut.lineIndices(i - 1);
            %                     end
            %
            %                     % Avoid double recording of a primary point if
            %                     % polylines are joint
            %                     if nextPrimary == lastPrimary
            %                         primaries = [primaries, nextPrimary];
            %                     else
            %                         if ~isempty(nextPrimary)
            %                             primaries = [primaries, nextPrimary];
            %                         end
            %
            %                         if ~isempty(lastPrimary)
            %                             primaries = [primaries, lastPrimary];
            %                         end
            %                     end
            %                 else
            %                     if ilsOut.lineIndices(i + 1) ~= -1
            %                         loops = [loops; ...
            %                             [ilsOut.lineIndices(i), ilsOut.lineIndices(i + 1)]];
            %                     end
            %                 end
            %             end

            loopsCell = {};
            loopsCell{1, 1} = loops + 1;

            % Set member variables of PMTemplate object
            PMTemplate = PointModelTemplate;
            PMTemplate = set(PMTemplate, 'imageName', volumeName);
            PMTemplate = set(PMTemplate, 'imageDir', 'Templates');
            PMTemplate = set(PMTemplate, 'name', templateName);
            PMTemplate = set(PMTemplate, 'direc', 'Templates');
            PMTemplate = set(PMTemplate, 'pts', vertices);
            PMTemplate = set(PMTemplate, 'primaries', primaries);
            PMTemplate = set(PMTemplate, 'loops', loopsCell);

            % Save the template
            save(fullfile(templatePath, templateName), 'PMTemplate');

            % Store the PointModelTemplate extensions
            pmtExtensions.lineIndices = ilsOut.lineIndices;
            pmtExtensions.vertexColours = ilsOut.vertexColours;
            pmtExtensions.lineColours = ilsOut.lineColours;
            pmtExtensions.volumeSize = size(v.volume);


            % Save the template extensions
            [templatePath, tName, extension, versn] = ...
                fileparts(fullfile(templatePath, templateName));
            templateExtensionsName = [tName, '.temp_ex_dat'];
            save(fullfile(templatePath, templateExtensionsName), 'pmtExtensions');
        end
    end
end
