Added OTTv2 code for STL/OBJ shapes

+ott/+shapes/AxisymLerp.m

@@ -62,8 +62,9 @@
       % Calculate the maximum particle radius
       r = max(sqrt(shape.rho.^2 + shape.z.^2));
     end
-    
+
     function v = get_volume(shape)
+      % Calculate the volume of the particle
       error('Not yet implemented');
     end
 

+ott/+shapes/Contents.m

@@ -1,4 +1,4 @@
-% ott.shapes objects describing basic shapes
+% ott.shapes Descriptions of geometric shapes
 %
 % Files
 %   AxisymLerp     - AxisymLerp a axisymmetric particle with lerping between points
@@ -10,3 +10,10 @@
 %   Sphere         - Sphere a simple sphere shape
 %   StarShape      - StarShape abstract class for star shaped particles
 %   Superellipsoid - Superellipsoid a simple superellipsoid shape
+%   StlLoader      - StlLoader load a shape from a STL file
+%   TriangularMesh - TriangularMesh base class for triangular mesh objects (such as file loaders)
+%   WavefrontObj   - WavefrontObj load a shape from a Wavefront OBJ file
+%
+% Copyright 2018 Isaac Lenton
+% This file is part of OTT, see LICENSE.md for information about
+% using/distributing this file.

+ott/+shapes/Cube.m

@@ -166,5 +166,26 @@
             shape, varargin{:});
       end
     end
+
+    function writeWavefrontObj(shape, filename)
+      % Write representation of shape to Wavefront OBJ file
+      %
+      % writeWavefrontObj(filename) writes the shape to the given file.
+
+      % Generate array of vertices
+      verts = [ 1, 1, 1; 1, 1, -1; ...
+                1, -1, 1; 1, -1, -1; ...
+                -1, 1, 1; -1, 1, -1; ...
+                -1, -1, 1; -1, -1, -1].' .* shape.width./2;
+
+      % Generate array of faces
+      % Order vertices so normals face outwards
+      faces = [ 1, 5, 7, 3; 2, 4, 8, 6; ...
+                1, 3, 4, 2; 3, 7, 8, 4; ...
+                7, 5, 6, 8; 5, 1, 2, 6 ].';
+
+      % Write the file
+      shape.writeWavefrontObj_helper(filename, verts, faces);
+    end
   end
 end

+ott/+shapes/Shape.m

@@ -1,8 +1,21 @@
 classdef Shape
 %Shape abstract class for optical tweezers toolbox shapes
 %
+% Properties
+%   maxRadius     maximum distance from shape origin
+%   volume        volume of shape
+%
 % Methods (abstract):
-%   inside(shape, ...) determine if point is inside shape
+%   inside(shape, ...) determine if spherical point is inside shape
+%
+% Methods:
+%   writeWavefrontObj(shape, ...) write shape to Wavefront OBJ file
+%       only implemented if shape supports this action.
+%   insideXyz(shape, ...) determine if Cartesian point is inside shape
+%       requires inside(shape, ...) to be implemented.
+%   simple(...) simplified constructor for shape-like objects.
+%
+% See also simple, ott.shapes.Cube, ott.shapes.TriangularMesh.
 %
 % This file is part of the optical tweezers toolbox.
 % See LICENSE.md for information about using/distributing this file.
@@ -24,6 +37,8 @@
       %   'cone-tipped-cylinder'      [ radius height cone_height ]
       %   'cube'            Cube [ width ]
       %   'axisym'          Axis-symetric particle [ rho(:) z(:) ]
+      %   'stl'             load STL file [filename]
+      %   'obj'             load Wavefront OBJ file [filename]
 
       switch lower(name)
         case 'sphere'
@@ -49,16 +64,168 @@
           shape = ott.shapes.Cube(parameters(:));
         case 'axisym'
           shape = ott.shapes.AxisymLerp(parameters(:));
+        case 'obj'
+          shape = ott.shapes.WavefrontObj(parameters(:));
+        case 'stl'
+          shape = ott.shapes.StlLoader(parameters(:));
         otherwise
           error('Unsupported simple particle shape');
       end
     end
   end
 
+  properties (Dependent)
+    maxRadius       % Maximum particle radius (useful for Nmax calculation)
+    volume          % Volume of the particle
+  end
+
   methods (Abstract)
     inside(shape, varargin)   % Determine if point is inside shape
+
+    get_maxRadius(shape, varargin)    % Get max distance from origin
+    get_volume(shape, varargin);      % Get shape volume
+  end
+
+  methods (Access=protected)
+
+    function writeWavefrontObj_helper(shape, filename, verts, faces)
+      % Helper to write faces and vertices to a file
+      %
+      % helper(filename, verts, faces)
+      %   verts 3xN array of floats
+      %   faces mxN array of integers (starting at 1)
+
+      fp = fopen(filename, 'w');
+
+      % Write documentation
+      fprintf(fp, '# Shape description generated by OTT\n');
+
+      % Write verts
+      for ii = 1:size(verts, 2)
+        fprintf(fp, 'v %f %f %f\n', verts(:, ii));
+      end
+
+      % Write faces
+      for ii = 1:size(faces, 2)
+        fprintf(fp, 'f ');
+        fprintf(fp, '%d ', faces(:, ii));
+        fprintf(fp, '\n');
+      end
+
+      % Close file
+      fclose(fp);
+
+    end
   end
 
   methods
+
+    function r = get.maxRadius(shape)
+      % Get the particle max radius
+      r = shape.get_maxRadius();
+    end
+
+    function r = get.volume(shape)
+      % Get the particle volume
+      r = shape.get_volume();
+    end
+
+    function writeWavefrontObj(shape, filename)
+      % Write representation of shape to Wavefront OBJ file
+      %
+      % This is a placeholder, not supported by all shape types
+      error('Shape does not support writing to WavefrontOBJ file');
+    end
+
+    function surf(shape, varargin)
+      % SURF generate a visualisation of the shape
+      %
+      % SURF() displays a visualisation of the shape in the current figure.
+      %
+      % SURF(..., 'surfoptions', {varargin}) specifies the options to
+      % pass to the surf function.
+      error('Shape does not support surf visualisation');
+    end
+
+    function varargout = voxels(shape, varargin)
+      % Generate an array of xyz coordinates for voxels inside the shape
+      %
+      % voxels(spacing) shows a visualisation of the shape with
+      % circles placed at locations on a Cartesian grid.
+      %
+      % xyz = voxels(spacing) returns the voxel locations.
+      %
+      % Optional named arguments:
+      %     'plotoptions'   Options to pass to the plot3 function
+      %     'visualise'     Show the visualisation (default: nargout == 0)
+
+      p = inputParser;
+      p.addOptional('spacing', shape.maxRadius/10);
+      p.addParameter('plotoptions', []);
+      p.addParameter('visualise', nargout == 0);
+      p.parse(varargin{:});
+
+      plotoptions = p.Results.plotoptions;
+      if isempty(plotoptions)
+        plotoptions = {...
+          'MarkerFaceColor', 'w', ...
+          'MarkerEdgeColor', [.5 .5 .5], ...
+          'MarkerSize', 20*p.Results.spacing/shape.maxRadius};
+      end
+
+      % Calculate range of dipoles
+      numr = ceil(shape.maxRadius / p.Results.spacing);
+      rrange = (-numr:numr)*p.Results.spacing;
+
+      % Generate the voxel grid
+      [xx, yy, zz] = meshgrid(rrange, rrange, rrange);
+
+      % Determine which points are inside
+      mask = shape.insideXyz(xx, yy, zz);
+      xyz = [xx(mask).'; yy(mask).'; zz(mask).'];
+
+      % Visualise the result
+      if p.Results.visualise
+        plot3(xyz(1,:), xyz(2,:), xyz(3,:), 'o', plotoptions{:});
+        axis equal
+        title(['spacing = ' num2str(p.Results.spacing) ...
+            ', N = ' int2str(sum(mask))])
+      end
+
+      % Assign output
+      if nargout ~= 0
+        varargout = {xyz};
+      end
+    end
+
+    function b = insideXyz(shape, x, y, z)
+      % INSIDEXYZ determine if Cartesian point is inside the shape
+      %
+      % b = inside(shape, x, y, z) determine if the Cartesian point
+      % [x, y, z] is inside the star shaped object.
+      %
+      % b = insideXyz(shape, xyz) as above, but using a 3xN matrix of
+      % [x; y; z] positions.
+      %
+      % See also INSIDE.
+
+      % Ensure the sizes match
+      if nargin == 4
+        x = x(:);
+        y = y(:);
+        z = z(:);
+        [x, y, z] = ott.utils.matchsize(x, y, z);
+      else
+        y = x(2, :);
+        z = x(3, :);
+        x = x(1, :);
+      end
+
+      % Convert to spherical coordinates
+      [r, t, p] = ott.utils.xyz2rtp(x, y, z);
+
+      % Call the spherical coordinate version
+      b = shape.inside(r, t, p);
+    end
   end
 end

+ott/+shapes/StarShape.m

@@ -10,21 +10,10 @@
 % This file is part of the optical tweezers toolbox.
 % See LICENSE.md for information about using/distributing this file.
 
-  properties
-  end
-
-  properties (Dependent)
-    maxRadius       % Maximum particle radius (useful for Nmax calculation)
-    volume          % Volume of the particle
-  end
-
   methods (Abstract)
     radii(shape, theta, phi);
     normals(shape, theta, phi);
     axialSymmetry(shape);
-
-    get_maxRadius(shape, varargin);
-    get_volume(shape, varargin);
   end
 
   methods
@@ -50,16 +39,6 @@
       end
     end
 
-    function r = get.maxRadius(shape)
-      % Get the particle max radius
-      r = shape.get_maxRadius();
-    end
-
-    function r = get.volume(shape)
-      % Get the particle volume
-      r = shape.get_volume();
-    end
-
     function varargout = locations(shape, theta, phi)
       % LOCATIONS calculates Cartessian coordinate locations for points
       %