--- a
+++ b/vol3d.m
@@ -0,0 +1,256 @@
+function [model] = vol3d(varargin)
+
+%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.
+%
+% vol3d                   Provide a demo of functionality.
+%
+% 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.
+%                         XxYxZ array represents scaled colormap indices.
+%                         XxYxZx3 array represents truecolor RGB values for
+%                         each voxel (along the 4th dimension).
+%
+% vol3d(...,'Alpha',alpha) XxYxZ array of alpha values for each voxel, in
+%                          range [0,1]. Default: data (interpreted as
+%                          scaled alphamap indices).
+%
+% vol3d(...,'Parent',axH)  Specify parent axes. Default: gca.
+%
+% vol3d(...,'XData',x)  1x2 x-axis bounds. Default: [0 size(data, 2)].
+% vol3d(...,'YData',y)  1x2 y-axis bounds. Default: [0 size(data, 1)].
+% vol3d(...,'ZData',z)  1x2 z-axis bounds. Default: [0 size(data, 3)].
+%
+% vol3d(...,'texture','2D')  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')  Default. 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 (from the original vol3d FEX submission) 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','3D');
+% view(3);  
+% axis tight;  daspect([1 1 .4])
+% alphamap('rampup');
+% alphamap(.06 .* alphamap);
+%
+% See also alphamap, colormap, opengl, isosurface
+
+% Copyright Joe Conti, 2004
+% Improvements by Oliver Woodford, 2008-2011, with permission of the
+% copyright holder.
+
+if nargin == 0
+    demo_vol3d;
+    return
+end
+
+if isstruct(varargin{1})
+    model = varargin{1};
+    if length(varargin) > 1
+       varargin = {varargin{2:end}};
+    end
+else
+    model = localGetDefaultModel;
+end
+
+if length(varargin)>1
+  for n = 1:2:length(varargin)
+    switch(lower(varargin{n}))
+        case 'cdata'
+            model.cdata = varargin{n+1};
+        case 'parent'
+            model.parent = varargin{n+1};
+        case 'texture'
+            model.texture = varargin{n+1};
+        case 'alpha'
+            model.alpha = varargin{n+1};
+        case 'xdata'
+            model.xdata = varargin{n+1}([1 end]);
+        case 'ydata'
+            model.ydata = varargin{n+1}([1 end]);
+        case 'zdata'
+            model.zdata = varargin{n+1}([1 end]);
+    end
+    
+  end
+end
+
+if isempty(model.parent)
+    model.parent = gca;
+end
+
+[model] = local_draw(model);
+
+%------------------------------------------%
+function [model] = localGetDefaultModel
+
+model.cdata = [];
+model.alpha = [];
+model.xdata = [];
+model.ydata = [];
+model.zdata = [];
+model.parent = [];
+model.handles = [];
+model.texture = '3D';
+tag = tempname;
+model.tag = ['vol3d_' tag(end-11:end)];
+
+%------------------------------------------%
+function [model,ax] = local_draw(model)
+
+cdata = model.cdata; 
+siz = size(cdata);
+
+% Define [x,y,z]data
+if isempty(model.xdata)
+    model.xdata = [0 siz(2)];
+end
+if isempty(model.ydata)
+    model.ydata = [0 siz(1)];
+end
+if isempty(model.zdata)
+    model.zdata = [0 siz(3)];
+end
+
+try
+   delete(model.handles);
+catch
+end
+
+ax = model.parent;
+cam_dir = camtarget(ax) - campos(ax);
+[m,ind] = max(abs(cam_dir));
+
+opts = {'Parent',ax,'cdatamapping',[],'alphadatamapping',[],'facecolor','texturemap','edgealpha',0,'facealpha','texturemap','tag',model.tag};
+
+if ndims(cdata) > 3
+    opts{4} = 'direct';
+else
+    cdata = double(cdata);
+    opts{4} = 'scaled';
+end
+
+if isempty(model.alpha)
+    alpha = cdata;
+    if ndims(model.cdata) > 3
+        alpha = sqrt(sum(double(alpha).^2, 4));
+        alpha = alpha - min(alpha(:));
+        alpha = 1 - alpha / max(alpha(:));
+    end
+    opts{6} = 'scaled';
+else
+    alpha = model.alpha;
+    if ~isequal(siz(1:3), size(alpha))
+        error('Incorrect size of alphamatte');
+    end
+    opts{6} = 'none';
+end
+
+h = findobj(ax,'type','surface','tag',model.tag);
+for n = 1:length(h)
+  try
+     delete(h(n));
+  catch
+  end
+end
+
+is3DTexture = strcmpi(model.texture,'3D');
+handle_ind = 1;
+
+% Create z-slice
+if(ind==3 || is3DTexture )    
+  x = [model.xdata(1), model.xdata(2); model.xdata(1), model.xdata(2)];
+  y = [model.ydata(1), model.ydata(1); model.ydata(2), model.ydata(2)];
+  z = [model.zdata(1), model.zdata(1); model.zdata(1), model.zdata(1)];
+  diff = model.zdata(2)-model.zdata(1);
+  delta = diff/size(cdata,3);
+  for n = 1:size(cdata,3)
+
+   cslice = squeeze(cdata(:,:,n,:));
+   aslice = double(squeeze(alpha(:,:,n)));
+   h(handle_ind) = surface(x,y,z,cslice,'alphadata',aslice,opts{:});
+   z = z + delta;
+   handle_ind = handle_ind + 1;
+  end
+
+end
+
+% Create x-slice
+if (ind==1 || is3DTexture ) 
+  x = [model.xdata(1), model.xdata(1); model.xdata(1), model.xdata(1)];
+  y = [model.ydata(1), model.ydata(1); model.ydata(2), model.ydata(2)];
+  z = [model.zdata(1), model.zdata(2); model.zdata(1), model.zdata(2)];
+  diff = model.xdata(2)-model.xdata(1);
+  delta = diff/size(cdata,2);
+  for n = 1:size(cdata,2)
+
+   cslice = squeeze(cdata(:,n,:,:));
+   aslice = double(squeeze(alpha(:,n,:)));
+   h(handle_ind) = surface(x,y,z,cslice,'alphadata',aslice,opts{:});
+   x = x + delta;
+   handle_ind = handle_ind + 1;
+  end
+end
+
+% Create y-slice
+if (ind==2 || is3DTexture)
+  x = [model.xdata(1), model.xdata(1); model.xdata(2), model.xdata(2)];
+  y = [model.ydata(1), model.ydata(1); model.ydata(1), model.ydata(1)];
+  z = [model.zdata(1), model.zdata(2); model.zdata(1), model.zdata(2)];
+  diff = model.ydata(2)-model.ydata(1);
+  delta = diff/size(cdata,1);
+  for n = 1:size(cdata,1)
+
+   cslice = squeeze(cdata(n,:,:,:));
+   aslice = double(squeeze(alpha(n,:,:)));
+   h(handle_ind) = surface(x,y,z,cslice,'alphadata',aslice,opts{:});
+   y = y + delta;
+   handle_ind = handle_ind + 1;
+  end
+end
+
+model.handles = h;
+
+function demo_vol3d
+figure;
+load mri.mat
+vol3d('cdata', squeeze(D), 'xdata', [0 1], 'ydata', [0 1], 'zdata', [0 0.7]);
+colormap(bone(256));
+alphamap([0 linspace(0.1, 0, 255)]);
+axis equal off
+set(gcf, 'color', 'w');
+view(3);
\ No newline at end of file