function a = plot3d(x, y, z, color)
%  a = plot3d(x, y, z, color)
%
% plot x, y, z data in three dimensions, with sliders attached so that you can adjust
% elevation and azimuth.
% Returns a handle to the figure so created.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% First: if x is numeric, then the person calling this routine wants us to create the figure
if ~ischar(x),
  if nargin<4,
    color = repmat('g',1,length(x));
  end
    
  % Create a UserData structure containing the specified numeric x, y, and z
  UserData = struct('x',x,'y',y,'z',z,'color',color);
  
  % Create the figure
  %
  a = figure('Units','normalized', ...
      'Color',[0.2 0.45 0.6], ...
      'Colormap',repmat([0:255]'/255, [1 3]), ...
      'Name','plot3d', ...
      'Pointer','crosshair', ...
      'Position',[0.3625 0.318359 0.592187 0.576172], ...
      'Renderer','zbuffer', ...
      'RendererMode','manual', ...
      'Tag','Fig1', ...
      'UserData',UserData, ...
      'WindowButtonDownFcn','LodImage zoomdown Image');

  % Sliders to adjust azimuth and elevation
  ImageView = [ 322.5 30 ];
  b = uicontrol('Parent',a, ...
      'Units','normalized', ...
      'BackgroundColor',[0.701961 0.701961 0.701961], ...
      'Callback','plot3d azimuth', ...
      'Max',360, ...
      'Min',0, ...
      'Position',[0.2 0.0 0.6 0.05], ...
      'SliderStep',[0.01 0.05], ...
      'Style','slider', ...
      'Tag','AzimuthSlider', ...
      'Value',ImageView(1));
  b = uicontrol('Parent',a, ...
      'Units','normalized', ...
      'BackgroundColor',[0.701961 0.701961 0.701961], ...
      'Position',[0.0 0.0 0.1 0.05 ], ...
      'String','Azimuth:', ...
      'Tag','AzimuthLabel');
  b = uicontrol('Parent',a, ...
      'Units','normalized', ...
      'BackgroundColor',[0.701961 0.701961 0.701961], ...
      'Position',[0.1 0.0 0.1 0.05 ], ...
      'String',sprintf('%3.3g',ImageView(1)), ...
      'Tag','AzimuthText');
  b = uicontrol('Parent',a, ...
      'Units','normalized', ...
      'BackgroundColor',[0.701961 0.701961 0.701961], ...
      'Callback','plot3d elevation', ...
      'Max',90, ...
      'Min',-90, ...
      'Position',[0.95 0.1 0.05 0.8 ], ...
      'SliderStep',[0.01 0.05], ...
      'Style','slider', ...
      'Tag','ElevationSlider', ...
      'Value', ImageView(2));
  b = uicontrol('Parent',a, ...
      'Units','normalized', ...
      'BackgroundColor',[0.701961 0.701961 0.701961], ...
      'Position',[0.9 0.95 0.1 0.05 ], ...
      'String','Elevation:', ...
      'Tag','ElevationLabel');
  b = uicontrol('Parent',a, ...
      'Units','normalized', ...
      'BackgroundColor',[0.701961 0.701961 0.701961], ...
      'Position',[0.9 0.9 0.1 0.05 ], ...
      'String',sprintf('%3.3g',ImageView(2)), ...
      'Tag','ElevationText');

  % Buttons to change the plot type
  b = uicontrol('Parent',a, ...
      'Units','normalized', ...
      'BackgroundColor',[0.701961 0.701961 0.701961], ...
      'Callback','plot3d plot3d',...
      'Position',[0.0 0.73 0.15 0.07 ], ...
      'String','Plot 3D', ...
      'Tag','Plot 3D Button');
  b = uicontrol('Parent',a, ...
      'Units','normalized', ...
      'BackgroundColor',[0.701961 0.701961 0.701961], ...
      'Callback','plot3d slice axials',...
      'Position',[0.0 0.66 0.15 0.07 ], ...
      'String','Plot Axials',...
      'Tag','Plot Axials Button');
  b = uicontrol('Parent',a, ...
      'Units','normalized', ...
      'BackgroundColor',[0.701961 0.701961 0.701961], ...
      'Callback','plot3d slice coronals',...
      'Position',[0.0 0.59 0.15 0.07 ], ...
      'String','Plot Coronals',...
      'Tag','Plot Coronals Button');
  b = uicontrol('Parent',a, ...
      'Units','normalized', ...
      'BackgroundColor',[0.701961 0.701961 0.701961], ...
      'Callback','plot3d slice sagittals',...
      'Position',[0.0 0.52 0.15 0.07 ], ...
      'String','Plot Sagittals',...
      'Tag','Plot Sagittals Button');
  
  % Text box with the name of the file in which we should save this thing
  OutputMarker = uicontrol('Parent',a, ...
      'Units','normalized', ...
      'BackgroundColor',[0.8 0.8 0.83], ...
      'Position',[0.1 0.95 0.25 0.05 ], ...
      'String', 'Save as File:',...
      'Style','text', ...
      'Tag','OutputMarker');
  OutputFile = uicontrol('Parent',a, ...
      'Units','normalized', ...
      'BackgroundColor',[0.8 0.8 0.83], ...
      'Position',[0.35 0.95 0.55 0.05 ], ...
      'String', fullfile(pwd, 'file1'),...
      'Style','edit', ...
      'Tag','OutputFile');

  % Save Menu, for saving this object as different file types
  b = uimenu('Parent',a, ...
      'Label','&Save', ...
      'Tag','SaveMenu');
  c = uimenu('Parent',b, ...
      'Callback','plot3d save otl', ...
      'Label','Save as OTL', ...
      'Tag','SaveOutlineOTLMenu');
  c = uimenu('Parent',b, ...
      'Callback','plot3d save txt', ...
      'Label','Save as TXT', ...
      'Tag','SaveOutlineTXTMenu');
  c = uimenu('Parent',b, ...
      'Callback','plot3d save csv', ...
      'Label','Save as CSV', ...
      'Tag','SaveOutlineCSVMenu');
  c = uimenu('Parent',b, ...
      'Callback','print gcbf ', ...
      'Label','Print', ...
      'Tag','PrintMenu');
  
  % Image Axes
  %
  ImageAxes = axes('Parent',a, ...
      'Box','on', ...
      'CameraUpVector',[0 1 0], ...
      'Color',[1 1 1], ...
      'Layer','top', ...
      'Position',[0.2 0.1 0.6 0.8], ...
      'Tag','ImageAxes', ...
      'XColor',[0 0 0], ...
      'XLim',[ 0.5 65536.5 ], ...
      'XLimMode','manual', ...
      'YColor',[0 0 0], ...
      'YLim',[ 0.5 65536.5 ], ...
      'YLimMode','manual', ...
      'ZColor',[0 0 0], ...
      'ZLim',[ 0.5 65536.5 ], ...
      'ZLimMode','manual');

  % Now plot the data
  plot3d_plot3d(ImageAxes);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Second: if x is a string, interpret it as a command
elseif ischar(x),
  % First, find the image axes
  ImageAxes = findobj(gcbf,'Tag','ImageAxes');
  
  % Now, switch based on the argument
  switch x,
    case 'azimuth',
      % Change the azimuth angle from which the object is viewed
      plot3d_changeview(ImageAxes, 1, 'AzimuthSlider', 'AzimuthText');
    case 'elevation',
      % Change the elevation angle from which the object is viewed
      plot3d_changeview(ImageAxes, 2, 'ElevationSlider', 'ElevationText');
    
    case 'save',
      % Save the object in the format specified in the 'y' argument
      plot3d_save(ImageAxes, y);
    case 'plot3d',
      % Re-plot the object in 3d plot
      plot3d_plot3d(ImageAxes);
    case 'slice',
      % Plot slices, in orientation specified
      if nargin<2, y = ''; end
      plot3d_slice(ImageAxes, y);
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plot3d_changeview(ImageAxes, Viewnum, Tag, TextTag)
% Change the view from which we display an object
%
ImageView = get(ImageAxes, 'View');
Object = findobj(gcbf,'Tag',Tag);
ImageView(Viewnum) = get(Object,'Value');
set(ImageAxes, 'View', ImageView);
TextWidget = findobj(gcbf,'Tag',TextTag);
set(TextWidget,'String',sprintf('%3.3g',ImageView(Viewnum)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plot3d_save(ImageAxes, fileformat)
% Save the object in a file with the specified file format
%

% List of possible color codes
COLOR_CODES = 'kbgcrmyw';

% Get the figure UserData
FU = get(get(ImageAxes,'Parent'), 'UserData');

% Get the output filename
OutputFile = findobj(get(ImageAxes,'Parent'), 'Tag', 'OutputFile');
if isempty(OutputFile),
  FileName = 'file1';
else
  FileName = get(OutputFile, 'String');
end

% Now switch, depending on the specified output format
switch fileformat,
  case 'otl',
    % In the case of 'otl' format, save as an outline file, with default header
    otlwrite([ FileName '.otl'], [FU.x; FU.y; FU.z; double(FU.color)]);
    
  case { 'csv' 'txt' },
    % In the case of 'csv' or 'txt' format, save as text
    % With 'csv', comma is the separator; with 'txt', tab is the separator
    separator=sprintf('\t');
    if strmatch(fileformat, {'csv','CSV'},'exact'), separator=','; end
    
    % Open the file
    fid=fopen(sprintf('%s.%s', FileName, fileformat ), 'w');
    if fid<0,
      disp(sprintf('Unable to open file %s.%s for writing', FileName, fileformat));
      return;
    end
    for i=1:length(FU.x),
      fprintf(fid, '%g%c%g%c%g%c%c\n', FU.x(i), separator, FU.y(i), separator, FU.z(i), separator, FU.color(i));
    end
    fclose(fid);
    disp(sprintf('Wrote outline in ASCII (%s format) to file %s.%s',FileName, fileformat)); 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plot3d_plot3d(ImageAxes)
% Re-plot the data in 3D
%

% Get the figure UserData
U = get(get(ImageAxes,'Parent'),'UserData');
x = U.x;
y = U.y;
z = U.z;
color = U.color;

% Plot each color, in turn, on the 3D plot
axes(ImageAxes);
while length(color)>0,
  which = (color==color(1));
  ImageDisplay = plot3(x(which),y(which),z(which),[color(1) '.']);
  set(ImageDisplay, 'Tag', 'ImageDisplay');
  set(get(ImageDisplay,'Parent'), 'Tag', 'ImageAxes')
  hold on;
  
  % Remove plotted data from the list, and get next plotted data
  x = x(which==0);
  y = y(which==0);
  z = z(which==0);
  color = color(which==0);
end

% Get the current Image View, and make sure Azimuth is within [0 360]
ImageView = get(ImageAxes,'View');
ImageView(1) = mod(ImageView(1), 360);

% Now make sure that the sliders have the correct endpoints, values, and callbacks
AzimuthSlider = findobj(get(ImageAxes,'Parent'),'Tag','AzimuthSlider');
set(AzimuthSlider,'Max',360,'Min',0,'Value',ImageView(1),...
    'SliderStep',[0.01 0.05],'Callback','plot3d azimuth');
ElevationSlider = findobj(get(ImageAxes,'Parent'),'Tag','ElevationSlider');
set(ElevationSlider,'Max',90,'Min',-90,'Value',ImageView(2),...
    'SliderStep',[0.01 0.05],'Callback','plot3d elevation');

% Reset the value of ElevationLabel
h = findobj(get(ImageAxes,'Parent'),'Tag','ElevationLabel');
set(h,'String','Elevation:');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plot3d_slice(ImageAxes, orientation)
% Re-plot the data in slices with the specified orientation
%
% Get the ElevationSlider handle
ElevationSlider = findobj(get(ImageAxes,'Parent'),'Tag','ElevationSlider');

% If orientation is anything other then '', set the new orientation
if length(orientation)>0,
  % Get the figure UserData
  FU = get(get(ImageAxes,'Parent'),'UserData');

  % Plot slices: first thing to do is figure out which one specifies level, and which ones to plot
  switch orientation,
    case 'axials',
      dat = [FU.x(:) FU.y(:) FU.z(:)];
    case 'coronals',
      dat = [FU.x(:) FU.z(:) FU.y(:)];
    case 'sagittals',
      dat = [FU.y(:) FU.z(:) FU.x(:)];
  end

  % Next, sort and specify the levels, and collect X, Y, and COLOR
  % For now, we round dat(:,3) to the nearest 0.01 of the possible range
  z = dat(:,3);
  RANGE = max(z) - min(z);
  if RANGE > 0,
    z = 0.01*RANGE * round(z/(0.01*RANGE));
  end
  levels = unique(z);
  X = {};
  Y = {};
  COLOR = {};
  L = {};
  for i=1:length(levels),
    L{i} = levels(i);
    which = find(z==levels(i));
    X{i} = dat(which,1);
    Y{i} = dat(which,2);
    COLOR{i} = FU.color(which);
  end

  % Find axis for these plots
  MM = [ min(dat(:,1:2)); max(dat(:,1:2)) ];
  axis_lims = [ [1.05 -0.05]*MM [-0.05 1.05]*MM ];
  axis_lims = axis_lims([1 3 2 4]);
  
  % Put the cell arrays and level specifications into a UserData structure for the axes
  AU = struct('L',L,'X',X,'Y',Y,'COLOR',COLOR,'axis',axis_lims);
  set(ImageAxes,'UserData',AU);
  
  % Next, set up the elevation slider so that it quickly moves back and forth through images
  set(ElevationSlider,'Max',length(levels),'Min',1,'Value',round(length(levels)/2),...
      'SliderStep',[1 2]/length(levels),'Callback','plot3d slice');
  
  % Clear callbacks from the AzimuthSlider
  h = findobj(get(ImageAxes,'Parent'),'Tag','AzimuthSlider');
  set(h,'Callback','');

  % Change the value of ElevationLabel
  h = findobj(get(ImageAxes,'Parent'),'Tag','ElevationLabel');
  set(h,'String','Coordinate:');  
end

% Get the current setting of ElevationSlider, and round it off
level_ind = round(get(ElevationSlider,'Value'));

% Get the Axes UserData structure; find the color, x, y, and level
AU = get(ImageAxes,'UserData');
color = AU(level_ind).COLOR;
x = AU(level_ind).X;
y = AU(level_ind).Y;
level = AU(level_ind).L;
axis_lims = AU(level_ind).axis;

% Clear any previous children of the image axes
for c=get(ImageAxes,'Children'),
  delete(c);
end

% Plot specified x's and y's with specified colors
axes(ImageAxes);
set(ImageAxes,'View',[0 90]);
hold on;
while length(color)>0,
  which = (color==color(1));
  ImageDisplay = plot(x(which),y(which),[color(1) 'x']);
  set(ImageDisplay, 'Tag', sprintf('ImageDisplay%c', color(1)));
  set(get(ImageDisplay,'Parent'), 'Tag', 'ImageAxes')
  
  % Remove plotted data from the list, and get next plotted data
  x = x(which==0);
  y = y(which==0);
  color = color(which==0);
end
axis(axis_lims);
hold off;

% Set ElevationText to specify the level coordinate
ElevationText = findobj(get(ImageAxes,'Parent'),'Tag','ElevationText');
set(ElevationText,'String',sprintf('%3.3g',level));
