function [output] = modelBuilder3DTool(varargin)
% modelBuilder3DTool M-file for modelBuilder3DTool.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 modelBuilder3DTool. ', ...
                '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(['ModelBuilder3DTool: Parameter ''projectDirectory'' not provided. \n', ...
        'Using current directory (.) as default']));
end

if isempty(ad.aam)
    ad.aam = AAM;
    warning(sprintf(['ModelBuilder3DTool: 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.viewVolumePshBtn, 'callback', {@viewVolumeCallback})
set(ad.guiHandles.editPointSetPshBtn, 'callback', {@editPointSetCallback})
% set(ad.guiHandles.volumesListBox, 'callback', {@editPointSetCallback})

% Setup GUI elements
setupGUIElements(ad);

% Save the application data
setappdata(0, 'modelBuilder3DToolData', 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, 'modelBuilder3DToolData');
output = ad.output;
rmappdata(0, 'modelBuilder3DToolData');

function viewVolumeCallback(cbo, event)

% display('viewVolumeCallback executing.');

ad = getappdata(0, 'modelBuilder3DToolData');

% 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.modelBuilder3DToolFig, 'visible', 'off');
drawnow expose;
viewVolume('path', fullfile(ad.projectDirectory, 'Cropped'));
set(ad.guiHandles.modelBuilder3DToolFig, 'visible', 'on');
setappdata(0, 'modelBuilder3DToolData', ad);

function editPointSetCallback(cbo, event)

display('editPointSetCallback executing.');

ad = getappdata(0, 'modelBuilder3DToolData');
ad = editPointSet(ad);
setappdata(0, 'modelBuilder3DToolData', 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 = editPointSet(ad)

showTemplate = get(ad.guiHandles.showTemplateCheckBox, 'Value');
showClipplanes = get(ad.guiHandles.showClippingPlanesCheckBox, 'Value');

indexedLineSetStructures = [];
volumeStructures = [];

% Get the template data that is needed by the volume of interest

% Reload statistical model
% statsModelFullName = fullfile(ad.projectDirectory, get(ad.aam, 'modelDirec'), ...
%     get(ad.aam, 'modelName'));
% m = load(statsModelFullName, '-mat');
% ad.aam = m.aam;

% Get the template data and load template volume
pmt = get(ad.aam, 'PointModelTemplate');
tVolumeName = get(pmt, 'imageName');
tVolumeDir = get(pmt, 'imageDir');
tName = get(pmt, 'name');
tPtsRealSpace = get(pmt, 'pts');
tLoops = get(pmt, 'loops');
% tPrimaries = get(pmt, 'primaries');
% tLoops = get(pmt, 'loops');

 % 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);
%tPts = tPts / 1000; % ChangeHere for real space

% 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

% Code used before
% voxelSize = getVoxelSizeFromString('string', tVolumeName);
% scaleMat = diag(voxelSize.^-1);
% tPts = scaleMat * tPtsRealSpace ./ 1000;

% End original code
% Begin replacement code

% See function newTemplate in file templateBuilder3DTool for an
% explanation of the scaling and the rationale behind it
tVoxelSize = getVoxelSizeFromString('string', tVolumeName);
tMinVoxelDim = min(tVoxelSize);

% Determine the display voxel size
tDisplayVoxelSize = tVoxelSize ./ tMinVoxelDim ./1000;

% Set the display size for template vertices
tPts = tPtsRealSpace ./ tMinVoxelDim ./ 1000;

% End replacement code
% End change js 081114

% Get PointModelTemplate extensions
templatePath = fullfile(ad.projectDirectory, 'Templates');
[fullTemplatePath, name, extension, versn] = ...
    fileparts(fullfile(templatePath, tName));
templateExtensionsName = [name, '.temp_ex_dat'];
temp = load(fullfile(ad.projectDirectory, 'Templates', templateExtensionsName), '-mat');
pmtExtensions = temp.pmtExtensions;
tLineIndices = pmtExtensions.lineIndices;
tVertexColours = pmtExtensions.vertexColours;
tLineColours = pmtExtensions.lineColours;

if showTemplate == 1


    tVolumeFullPath = fullfile(ad.projectDirectory, tVolumeDir, char(tVolumeName));
    tVolume = loadVolumeFromSlicesDir('fullvolumepath', tVolumeFullPath);

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

    % Initialise volume structure for the template
    v1 = getHTVolume();
    v1.name = 'tVolume';
    v1.parentTransform = 'tVolumeTransform';
    v1.volume = tVolume;
    v1.planes = 256;
    
    % Begin change js 081114
    % Begin addition
    v1.voxelSize = tDisplayVoxelSize;
    % End addition
    % End change js 081114
    
    % Initialise point set structure for the template
    % p1 = getHTPointSet();
    % p1.name = 'tPointSet';
    % p1.parentTransform = 'tVolumeTransform';
    % p1.points = tPts;
    % p1CData = zeros(size(p1.points));
    % p1CData(1, :) = 1;
    % p1.colours = p1CData;
    % p1.interactionMode = 'none';

    ils1 = getHTIndexedLineSet();
    ils1.name = 'tIndexedLineSet';
    ils1.parentTransform = 'tVolumeTransform';
    ils1.vertices = tPts;
    % ils.lineIndices = [-1, 0:length(ils.vertices) - 1, -1];
    ils1.lineIndices = tLineIndices;
    ils1.interactionMode = 'none';
    ils1.vertexColours = tVertexColours;
    % ils1.vertexColours(3, :) = 1;
    ils1.lineColours = tLineColours;
    % ils1.lineColours(1, :) = 1;
    ils1.vertexStiffness = 0.5;
    ils1.lineStiffness = 0.5;

    indexedLineSetStructures = [indexedLineSetStructures, ils1];
    volumeStructures = [volumeStructures, v1];
end
% Load the volume of interest
index = get(ad.guiHandles.volumesListBox, 'Value');
iVolumeNames = get(ad.guiHandles.volumesListBox, 'String');
iVolumeName = iVolumeNames{index};
iVolumeFullPath = fullfile(ad.projectDirectory, 'Cropped', iVolumeName);

% Initialise volume structure for volume of interest
iVolume = loadVolumeFromSlicesDir('fullvolumepath', iVolumeFullPath);

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

v2 = getHTVolume();
v2.name = 'iVolume';
v2.parentTransform = 'iVolumeTransform';
v2.volume = iVolume;
v2.planes = 256;

% Begin change js 081114
% Begin addition js 081114

% See function newTemplate in file templateBuilder3DTool for an
% explanation of the scaling and the rationale behind it
iVoxelSize = getVoxelSizeFromString('string', iVolumeName);
iMinVoxelDim = min(iVoxelSize);

% Determine the display voxel size and set it
iDisplayVoxelSize = iVoxelSize ./ iMinVoxelDim ./1000;
v2.voxelSize = iDisplayVoxelSize;

% End addition
% End change js 081114

% Get voxel size from volume name and scale down as voxel sizes in H3DApi
% are in metres
% v2.voxelSize = getVoxelSizeFromString('string', iVolumeName) / 1000;

% Initialise point set structure for volume of interest
% p2 = getHTPointSet();
% p2.name = 'iPointSet';
% p2.parentTransform = 'iVolumeTransform';
% p2.points = tPts; % Get the real points here!!!!!!!!!!

% Begin commented 071116
% Begin code modified from function UpdateTarget in file ModelBuilder.m
% ce = get_element(ad.aam, index);
% cf = get(ce, 'imagefilename');
% [path, cfname, ext, vers] = fileparts(cf); %Not really necessary
% [path, tNameNoExt, ext, vers] = fileparts(tName);
% pointModelFullName = fullfile(ad.projectDirectory, 'PointModels', ...
%     tNameNoExt, [cfname, '_pm.mat']);
% Begin commented 071116


[path, tNameNoExt, ext, vers] = fileparts(tName);
pointModelFullName = fullfile(ad.projectDirectory, 'PointModels', ...
    tNameNoExt, [iVolumeName, '_pm.mat']);

if exist(pointModelFullName)
    
    pm = load(pointModelFullName);
    vertices = pm.pts;

    if isempty(vertices)
        vertices = ils1.vertices;
    else

        % In case the points got saved as an n x 1 matrix and not
        % a 3 x m matrix. In practice this will be the case as we follow
        % the convention used in AutoSave in file ModelBuilder.m, where
        % the point model is saved to disk as an n x 1 matrix.
        if size(vertices, 2) ~= 3 % 3D points
            vertices = reshape(vertices, 3, length(vertices) / 3);
        end

         % 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
        
        % 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', iVolumeName);
%         scaleMat = diag(voxelSize.^-1);
%         vertices = scaleMat * vertices ./ 1000;

        % End original code
        % Begin replacement code

        % Set the display size for template vertices
        vertices = vertices ./ iMinVoxelDim ./ 1000;

        % End replacement code
        % End change js 081114


    end
else
    vertices = tPts;
end

% End code modified from function UpdateTarget in file ModelBuilder.m

% p2.points = vertices;
% p2CData = zeros(size(p2.points));
% p2CData(1, :) = 1;
% p2.colours = p1CData;
% p2.interactionMode = 'move';

ils2 = getHTIndexedLineSet();
ils2.name = 'iIndexedLineSet';
ils2.parentTransform = 'iVolumeTransform';
ils2.vertices = vertices;
% ils.lineIndices = [-1, 0:length(ils.vertices) - 1, -1];
ils2.lineIndices = tLineIndices;
ils2.interactionMode = 'none';
ils2.vertexColours = tVertexColours;
% ils2.vertexColours(3, :) = 1;
ils2.lineColours = tLineColours;
% ils2.lineColours(1, :) = 1;
ils2.vertexStiffness = 0.5;
ils2.lineStiffness = 0.5;

cpt = getHTTransform();
cpt.name = 'clipPlaneTransform';
cpt.parentTransform = 'iVolumeTransform';
% transforms = [transforms cpt];

% cp1 = getHTClipPlane();
% cp1.name = 'clipPlane1';
% cp1.parentTransform = 'iVolumeTransform';
% cp1.point = [0, 0, 0];
% cp1.normal = [0, 1, 0];
% cp1.interactionMode = 'edit';
% % clipPlanes(1) = cp1;

noOfClipPlanes = 4;
cpXOffset = 0.03;
cpXBasePos = 0 - (noOfClipPlanes * cpXOffset / 2);
cpYBasePos = -size(iVolume, 2) / 1000 / 2;
% cpYBasePos = 0;
cpZBasePos = -size(iVolume, 3) / 1000 / 2;

for i = 1:noOfClipPlanes

    cp = getHTClipPlane();
    cp.name = ['clipPlane', num2str(i)];
    cp.parentTransform = 'iVolumeTransform';
    cp.point = [cpXBasePos + cpXOffset * (i - 1), cpYBasePos, cpZBasePos];
    cp.normal = [0, 1, 0];
    cp.interactionMode = 'edit';
    clipPlaneStructures(i) = cp;
end

if noOfClipPlanes == 0
    clipPlaneStructures = [];
end

% The following is deprecated code
% pointSetStructures = [p1, p2];
% indexedLineSetStructures = [ils1, ils2];
% volumeStructures = [v1, v2];

indexedLineSetStructures = [indexedLineSetStructures, ils2];
volumeStructures = [volumeStructures, v2];

% Initialise global scale transform to scale to screen size
tScale = getHTTransform();
tScale.name = 'globalScaleTransform';

if showTemplate == 1

    volumeTransforms = getHorizontalDistributionTransforms(volumeStructures, ...
        size(iVolume, 2) / 1000 / 20, 'globalScaleTransform'); % 0.03
    transformStructures = [tScale, volumeTransforms];
else
    iVolumeTransform = getHTTransform();
    iVolumeTransform.name = 'iVolumeTransform';
    iVolumeTransform.parentTransform = 'globalScaleTransform';
    transformStructures = [tScale, iVolumeTransform];

end

if showClipplanes == 0;
    clipPlaneStructures = [];
end

% clipPlaneStructures = [cp1];

% 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.modelBuilder3DToolFig, 'visible', 'off');
drawnow expose;

[transformsOut, volumesOut, pointSetsOut, indexedLineSetsOut, clipPlanesOut] =...
    hapticTool(transformStructures, volumeStructures, [],...
    indexedLineSetStructures, ...
    clipPlaneStructures);

% The following does not render line sets and clipping planes
% [transformsOut, volumesOut, pointSetsOut, indexedLineSetsOut, clipPlanesOut] =...
%     hapticTool(transformStructures, volumeStructures, [],...
%     indexedLineSetStructures, ...
%     clipPlaneStructures);

% [transformsOut, volumesOut, pointSetsOut, indexedLineSetsOut, clipPlanesOut] =...
%         hapticTool(transformStructures, volumeStructures, [],...
%         indexedLineSetStructures, ...
%         []);
%

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

% 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).
% adjustedPoints = indexedLineSetsOut(2).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(voxelSize1) * vertices .* 1000

% Used code before
% voxelSize = getVoxelSizeFromString('string', iVolumeName);
% scaleMat = diag(voxelSize);
% adjustedPoints = scaleMat * indexedLineSetsOut(end).vertices .* 1000;

% End original code
% Begin replacement code
iVoxelSize = getVoxelSizeFromString('string', iVolumeName);
iMinVoxelDim = min(iVoxelSize);
adjustedPoints = indexedLineSetsOut(end).vertices .* iMinVoxelDim .* 1000;

% End replacement code
% End change js 081114

% adjustedPoints = pointSetsOut(2).points;

% Begin commented 071116
% ce = set(ce, 'pts', adjustedPoints);
% ad.aam = set_element(ad.aam, ce, index);
%
% % Save statistical model
% statsModelFullName = fullfile(ad.projectDirectory, get(ad.aam, 'modelDirec'), ...
%     get(ad.aam, 'modelName'));
% aam = ad.aam;
% save(statsModelFullName, 'aam');
% End commented 071116

% Save point model
pts = adjustedPoints(:); % As done in AutoSave
save(pointModelFullName, 'pts');

