function [volumeRenderModel] = renderVolume(varargin)
% function [volumeRenderModel] = renderVolume(varargin)
% 
% Based on Vol3d Copyright Joe Conti, 2004, The Mathworks
% Source: Matlab Central File Exchange
% 
% Adapted by Johann Strasser for rgb volume rendering.
%
% Take the maximum intensity as the alpha value, i.e. for each pixel
% alpha = max([r, g, b])
% When AlphaDataMapping is set to direct we use max(slice, [], 3) and set
% an appropriate alphamap for the figure
% Do not set to scaled (default) as each slice's alpha data would then be
% scaled to fit into the limits defined by the 'ALim' property of the axes.
% This results in slices in a volume getting too opaque.
% Note that the point of using alphamaps is to be able to change the alpha
% mapping. Note that even we pass a uint8 array als AlphaData, it is
% evident from the plot editor at runtime that matlab casts it to double.
% CDataMapping is set to direct for similar reasons as setting
% AlphaDataMapping to direct. When a single channel volume is used the data
% in is used to index in the colourmap. If the default setting for
% CDataMapping was used the intensitiy values of each slice would be scaled
% separately to span the range of the colourmap. This leads to
% non-intuitive display of intensity values as the same value may be scaled
% differently in each slice. Note that it is assumed that the colourmap is
% setup correctly (i.e.) it is set to the range of the intensity values.
% Furthermore note that the alphamap is a property of the figure. A
% colourmap on the other hand is bound to an axes!
%
% For original vol3d description see below
%
%H = VOL3D Volume render 3-D data. 
% VOL3D uses the orthogonal plane 2-D texture mapping technique for 
% volume rending 3-D data in OpenGL. Use the 'texture' option to fine 
% tune the texture mapping technique. This function is best used with
% fast OpenGL hardware.
%
% H = vol3d('CData',data) Create volume render object from input 
%                         3-D data. Use interp3 on data to increase volume
%                         rendering resolution. Returns a struct 
%                         encapsulating the pseudo-volume rendering object. 
%
% vol3d(...'Parent',axH) Specify parent axes
%
% vol3d(...,'texture','2D')  Default. Only render texture planes parallel
%                            to nearest orthogonal viewing plane. Requires
%                            doing vol3d(h) to refresh if the view is
%                            rotated (i.e. using cameratoolbar).
%
% vol3d(...,'texture','3D')  Render x,y,z texture planes simultaneously. 
%                            This avoids the need to refresh the view but 
%                            requires faster OpenGL hardware peformance.
%
% vol3d(H)  Refresh view. Updates rendering of texture planes 
%           to reduce visual aliasing when using the 'texture'='2D'
%           option.
%
% NOTES
% Use vol3dtool for editing the colormap and alphamap. 
% Adjusting these maps will allow you to explore your 3-D volume 
% data at various intensity levels. See documentation on 
% alphamap and colormap for more information.
%
% Use interp3 on input date to increase/decrease resolution of data
%
% Examples:
%
% % Visualizing fluid flow
% v = flow(50);
% h = vol3d('cdata',v,'texture','2D');
% view(3); 
% % Update view since 'texture' = '2D'
% vol3d(h);  
% alphamap('rampdown'), alphamap('decrease'), alphamap('decrease')
% 
% % Visualizing MRI data
% load mri.mat
% D = squeeze(D);
% h = vol3d('cdata',D,'texture','2D');
% view(3); 
% % Update view since 'texture' = '2D'
% vol3d(h);  
% axis tight;  daspect([1 1 .4])
% alphamap('rampup');
% alphamap(.06 .* alphamap);
%
% See also vol3dtool, alphamap, colormap, opengl, isosurface

% Copyright Joe Conti, 2004

if isstruct(varargin{1})
    volumeRenderModel = varargin{1};
    if length(varargin) > 1
       varargin = {varargin{2:end}};
    end
else
    volumeRenderModel = getVolumeRenderModel();
end


if length(varargin)>1
    for n = 1:2:length(varargin)
        switch(lower(varargin{n}))
            case 'cdata'
                volumeRenderModel.cdata = varargin{n+1};
            case 'axis'
                volumeRenderModel.axis = varargin{n+1};
            case 'parent'
                volumeRenderModel.parent = varargin{n+1};
            case 'planespacing'
                volumeRenderModel.planeSpacing = varargin{n+1};
            case 'texture'
                volumeRenderModel.texture = varargin{n+1};
        end

    end
end

if isempty(volumeRenderModel.axis)
    volumeRenderModel.parent = gca;
end

if isempty(volumeRenderModel.parent)
    volumeRenderModel.parent = volumeRenderModel.axis;
end

% choose default 3-D view
ax = volumeRenderModel.axis;
% axis(ax,'vis3d');
% axis(ax,'tight');

[volumeRenderModel] = local_draw(volumeRenderModel);


function [volumeRenderModel,ax] = local_draw(volumeRenderModel)

cdata = volumeRenderModel.cdata; 
siz = size(cdata);

% Define [x,y,z]data
if isempty(volumeRenderModel.xdata)
    volumeRenderModel.xdata = [0 siz(2)];
end
if isempty(volumeRenderModel.ydata)
    volumeRenderModel.ydata = [0 siz(1)];
end
if isempty(volumeRenderModel.zdata)
    volumeRenderModel.zdata = [0 siz(3)];
end

% try,
%    delete(volumeRenderModel.handles);
% end

ax = volumeRenderModel.axis;
parent = volumeRenderModel.parent;
cam_dir = camtarget(ax) - campos(ax);
[m,ind] = max(abs(cam_dir));


% h = findobj(ax,'type','surface','tag','vol3d');
% for n = 1:length(h)
%   try,
%      delete(h(n));
%   end
% end

is3DTexture = strcmpi(volumeRenderModel.texture,'3D');
handleInd = 1;

% Create z-slice
if(ind==3 || is3DTexture )    
  x = [volumeRenderModel.xdata(1), volumeRenderModel.xdata(2); volumeRenderModel.xdata(1), volumeRenderModel.xdata(2)];
  y = [volumeRenderModel.ydata(1), volumeRenderModel.ydata(1); volumeRenderModel.ydata(2), volumeRenderModel.ydata(2)];
  z = [volumeRenderModel.zdata(1), volumeRenderModel.zdata(1); volumeRenderModel.zdata(1), volumeRenderModel.zdata(1)];
  diff = volumeRenderModel.zdata(2)-volumeRenderModel.zdata(1);
  delta = diff/size(cdata,3);
  
  for n = 1:volumeRenderModel.planeSpacing(3):size(cdata,3)      
      slice = squeeze(cdata(:, :, n, :));
      
      if size(slice, 3) == 3 % True colour slice
          alphaData = max(slice, [], 3);
      else
          if islogical(slice)
              slice = uint8(slice) * 255;
          end
          
          alphaData = slice; % Grayscale, hence use intensities as alpha indices
          %slice = repmat(slice, [1, 1, 3]);
      end

      if isequal(volumeRenderModel.mode, 'new')
          
          % Note that CDataMapping set to direct will only have an
          % effect if a single channel is specified.
          volumeRenderModel.zPlaneHandles(handleInd) = surface(x,y,z,'Parent',parent);
          set(volumeRenderModel.zPlaneHandles(handleInd),'facecolor','texture','cdata', slice,...
              'FaceAlpha', 'texturemap', 'AlphaData', alphaData, ...
              'AlphaDataMapping', 'direct', 'EdgeColor', 'none', 'CDataMapping', 'direct'); %, 'Tag', 'vol3d', 'LineWidth', 1);

      elseif isequal(volumeRenderModel.mode, 'update')

          set(volumeRenderModel.zPlaneHandles(handleInd), ...
              'CData', slice, 'AlphaData', alphaData);
      end

      z = z + delta * volumeRenderModel.planeSpacing(3);
      handleInd = handleInd + 1;
  end

end

handleInd = 1;

% Create x-slice
if (ind==1 || is3DTexture )
    x = [volumeRenderModel.xdata(1), volumeRenderModel.xdata(1); volumeRenderModel.xdata(1), volumeRenderModel.xdata(1)];
    y = [volumeRenderModel.ydata(1), volumeRenderModel.ydata(1); volumeRenderModel.ydata(2), volumeRenderModel.ydata(2)];
    z = [volumeRenderModel.zdata(1), volumeRenderModel.zdata(2); volumeRenderModel.zdata(1), volumeRenderModel.zdata(2)];
    diff = volumeRenderModel.xdata(2)-volumeRenderModel.xdata(1);
    delta = diff/size(cdata,2);
    for n = 1:volumeRenderModel.planeSpacing(1):size(cdata,2)
        slice = squeeze(cdata(:, n, :, :));
        
        if size(slice, 3) == 3 % True colour slice
            alphaData = max(slice, [], 3);
        else            
            if islogical(slice)
                slice = uint8(slice) * 255;
            end

            alphaData = slice; % Grayscale, hence use intensities as alpha indices
            %slice = repmat(slice, [1, 1, 3]);
        end

        if isequal(volumeRenderModel.mode, 'new')

            volumeRenderModel.xPlaneHandles(handleInd) = surface(x,y,z,'Parent',parent);
            
            % Note that CDataMapping set to direct will only have an
          % effect if a single channel is specified.
            set(volumeRenderModel.xPlaneHandles(handleInd),'facecolor','texture','cdata', slice,...
                'FaceAlpha', 'texturemap', 'AlphaData', alphaData, ...
                'AlphaDataMapping', 'direct', 'EdgeColor', 'none', 'CDataMapping', 'direct'); %, 'Tag', 'vol3d', 'LineWidth', 1);

        elseif isequal(volumeRenderModel.mode, 'update')

            set(volumeRenderModel.xPlaneHandles(handleInd), ...
                'CData', slice, 'AlphaData', alphaData);
        end

        x = x + delta * volumeRenderModel.planeSpacing(1);
        handleInd = handleInd + 1;
    end
end

handleInd = 1;

% Create y-slice
if (ind==2 || is3DTexture)
    x = [volumeRenderModel.xdata(1), volumeRenderModel.xdata(1); volumeRenderModel.xdata(2), volumeRenderModel.xdata(2)];
    y = [volumeRenderModel.ydata(1), volumeRenderModel.ydata(1); volumeRenderModel.ydata(1), volumeRenderModel.ydata(1)];
    z = [volumeRenderModel.zdata(1), volumeRenderModel.zdata(2); volumeRenderModel.zdata(1), volumeRenderModel.zdata(2)];
    diff = volumeRenderModel.ydata(2)-volumeRenderModel.ydata(1);
    delta = diff/size(cdata,1);
    for n = 1:volumeRenderModel.planeSpacing(2):size(cdata,1)
        slice = squeeze(cdata(n, :, :, :));
        
      if size(slice, 3) == 3 % True colour slice
          alphaData = max(slice, [], 3);
      else
          if islogical(slice)
              slice = uint8(slice) * 255;
          end

          alphaData = slice; % Grayscale, hence use intensities as alpha indices
          %slice = repmat(slice, [1, 1, 3]);
      end

        if isequal(volumeRenderModel.mode, 'new')

            volumeRenderModel.yPlaneHandles(handleInd) = surface(x,y,z,'Parent',parent);
            
            % Note that CDataMapping set to direct will only have an
          % effect if a single channel is specified.
            set(volumeRenderModel.yPlaneHandles(handleInd),'facecolor','texture','cdata', slice,...
                'FaceAlpha', 'texturemap', 'AlphaData', alphaData, ...
                'AlphaDataMapping', 'direct', 'EdgeColor', 'none', 'CDataMapping', 'direct'); %, 'Tag', 'vol3d', 'LineWidth', 1);

        elseif isequal(volumeRenderModel.mode, 'update')
            
            set(volumeRenderModel.yPlaneHandles(handleInd), ...
                'CData', slice, 'AlphaData', alphaData);
        end

        y = y + delta * volumeRenderModel.planeSpacing(2);
        handleInd = handleInd + 1;
    end
end

%volumeRenderModel.handles = h;