diff --git a/+ott/+shapes/AxisymLerp.m b/+ott/+shapes/AxisymLerp.m index c588cacf..601934a6 100644 --- a/+ott/+shapes/AxisymLerp.m +++ b/+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 diff --git a/+ott/+shapes/Contents.m b/+ott/+shapes/Contents.m index 1c096e18..f6dce10e 100644 --- a/+ott/+shapes/Contents.m +++ b/+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. diff --git a/+ott/+shapes/Cube.m b/+ott/+shapes/Cube.m index 2eb613ff..07a8b17d 100644 --- a/+ott/+shapes/Cube.m +++ b/+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 diff --git a/+ott/+shapes/Shape.m b/+ott/+shapes/Shape.m index 40c29df0..512fd7c2 100644 --- a/+ott/+shapes/Shape.m +++ b/+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 diff --git a/+ott/+shapes/StarShape.m b/+ott/+shapes/StarShape.m index e0b37723..ab83f0d1 100644 --- a/+ott/+shapes/StarShape.m +++ b/+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 % diff --git a/+ott/+shapes/StlLoader.m b/+ott/+shapes/StlLoader.m new file mode 100644 index 00000000..c8fde262 --- /dev/null +++ b/+ott/+shapes/StlLoader.m @@ -0,0 +1,133 @@ +classdef StlLoader < ott.shapes.TriangularMesh +% StlLoader load a shape from a STL file +% +% This class uses a 3rd party STL file reader: +% https://au.mathworks.com/matlabcentral/fileexchange/22409-stl-file-reader +% +% See tplicenses/stl_EricJohnson.txt for information about licensing. +% +% This file is part of the optical tweezers toolbox. +% See LICENSE.md for information about using/distributing this file. + + properties (SetAccess=protected) + filename % Name of the file this object loaded + end + + methods (Hidden, Static) + function [F, V, N] = stlbinary(M) + % Function to read STL file + % + % This is from STL file reader version 1.2.0.0 (1.6 MB) by Eric Johnson + % See tplicenses/stl_EricJohnson.txt for information about licensing. + + F = []; + V = []; + N = []; + + if length(M) < 84 + error('MATLAB:stlread:incorrectFormat', ... + 'Incomplete header information in binary STL file.'); + end + + % Bytes 81-84 are an unsigned 32-bit integer specifying the + % number of faces that follow. + numFaces = typecast(M(81:84),'uint32'); + %numFaces = double(numFaces); + if numFaces == 0 + warning('MATLAB:stlread:nodata','No data in STL file.'); + return + end + + T = M(85:end); + F = NaN(numFaces,3); + V = NaN(3*numFaces,3); + N = NaN(numFaces,3); + + numRead = 0; + while numRead < numFaces + % Each facet is 50 bytes + % - Three single precision values specifying the face normal vector + % - Three single precision values specifying the first vertex (XYZ) + % - Three single precision values specifying the second vertex (XYZ) + % - Three single precision values specifying the third vertex (XYZ) + % - Two unused bytes + i1 = 50 * numRead + 1; + i2 = i1 + 50 - 1; + facet = T(i1:i2)'; + + n = typecast(facet(1:12),'single'); + v1 = typecast(facet(13:24),'single'); + v2 = typecast(facet(25:36),'single'); + v3 = typecast(facet(37:48),'single'); + + n = double(n); + v = double([v1; v2; v3]); + + % Figure out where to fit these new vertices, and the face, in the + % larger F and V collections. + fInd = numRead + 1; + vInd1 = 3 * (fInd - 1) + 1; + vInd2 = vInd1 + 3 - 1; + + V(vInd1:vInd2,:) = v; + F(fInd,:) = vInd1:vInd2; + N(fInd,:) = n; + + numRead = numRead + 1; + end + + end + end + + methods + function shape = StlLoader(filename) + % Construct a new shape from a STL file + % + % StlLoader(filename) loads the face and vertex information + % contained in the file. + % + % Only supports binary STL files. + % + % This function uses 3rd party code, + % see tplicenses/stl_EricJohnson.txt for licensing information. + + assert(nargin == 1, 'Filename not supplied'); + + % + % Begin third-party code + % + + if ~exist(filename,'file') + error(['File ''%s'' not found. If the file is not on MATLAB''s path' ... + ', be sure to specify the full path to the file.'], file); + end + + fid = fopen(filename,'r'); + if ~isempty(ferror(fid)) + error(lasterror); %#ok + end + + M = fread(fid,inf,'uint8=>uint8'); + fclose(fid); + + [faces,verts,~] = ott.shapes.StlLoader.stlbinary(M); + + % + % End third-party code + % + + % Convert to a mesh with unique vertices + + % Start by finding unique vertices + [uverts,~,new_idx] = unique(verts,'rows'); + + % Replace face indices with new indices + faces = new_idx(faces); + + % Setup object + shape = shape@ott.shapes.TriangularMesh(uverts.', faces.'); + shape.filename = filename; + + end + end +end diff --git a/+ott/+shapes/Superellipsoid.m b/+ott/+shapes/Superellipsoid.m index a8ef9922..d024872a 100644 --- a/+ott/+shapes/Superellipsoid.m +++ b/+ott/+shapes/Superellipsoid.m @@ -53,6 +53,11 @@ r = max([shape.a, shape.b, shape.c]); end + function v = get_volume(shape) + % Calculate the volume of the particle + error('Not yet implemented'); + end + function b = isSphere(shape) % ISSPHERE Returns true if the shape is a sphere b = shape.a == shape.b && shape.a == shape.c ... @@ -74,10 +79,6 @@ r = ( ast.^(2/shape.ns) .* cpsp.^(shape.ew/shape.ns) ... + (act/shape.c).^(2/shape.ns) ).^(-shape.ns/2); end - - function v = get_volume(shape) - error('Not yet implemented'); - end function n = normals(shape, theta, phi) % NORMALS calcualtes the normals for each requested point diff --git a/+ott/+shapes/TriangularMesh.m b/+ott/+shapes/TriangularMesh.m new file mode 100644 index 00000000..d69aed10 --- /dev/null +++ b/+ott/+shapes/TriangularMesh.m @@ -0,0 +1,144 @@ +classdef TriangularMesh < ott.shapes.Shape +% TriangularMesh base class for triangular mesh objects (such as file loaders) +% +% Properties (read-only): +% verts 3xN matrix of vertex locations +% faces 3xN matrix of vertex indices describing faces +% +% Faces vertices should be ordered so normals face outwards for +% volume and inside functions to work correctly. +% +% This file is part of the optical tweezers toolbox. +% See LICENSE.md for information about using/distributing this file. + + properties (SetAccess=protected) + verts % Matrix of vertices in the object + faces % Matrix of faces in the object + end + + methods + function shape = TriangularMesh(verts, faces) + % Construct a new triangular mesh representation + % + % TriangularMesh(verts, faces) + % verts 3xN matrix of vertex locations + % faces 3xN matrix of vertex indices describing faces + % + % Faces vertices should be ordered so normals face outwards for + % volume and inside functions to work correctly. + + shape = shape@ott.shapes.Shape(); + + % Verify size of inputs + assert(size(verts, 1) == 3, 'Verts must be matrix of 3xN'); + assert(size(faces, 1) == 3, 'Faces must be matrix of 3xN'); + + shape.verts = verts; + shape.faces = faces; + + % Verify we have sufficient verts for faces + if max(shape.faces(:)) > numel(shape.verts) + error('faces matrix refers to non-existent vertices'); + end + + end + + function r = get_maxRadius(shape) + % Calculate the maximum distance from the shape origin + r = max(sum(shape.verts.^2, 1).^0.5); + end + + function totalVolume = get_volume(shape) + % Calculate the volume of the shape + % + % This function is based on a surface triangulation volume code + % version 1.0.0.0 (1.43 KB) by Krishnan Suresh + % matlabcentral/fileexchange/26982-volume-of-a-surface-triangulation + % See tplicenses/stl_KrishnanSuresh.txt for information about licensing. + + p = shape.verts; + t = shape.faces; + + % Compute the vectors d13 and d12 + d13= [(p(1,t(2,:))-p(1,t(3,:))); (p(2,t(2,:))-p(2,t(3,:))); ... + (p(3,t(2,:))-p(3,t(3,:)))]; + d12= [(p(1,t(1,:))-p(1,t(2,:))); (p(2,t(1,:))-p(2,t(2,:))); ... + (p(3,t(1,:))-p(3,t(2,:)))]; + + cr = cross(d13,d12,1);%cross-product (vectorized) + + area = 0.5*sqrt(cr(1,:).^2+cr(2,:).^2+cr(3,:).^2);% Area of each triangle + totalArea = sum(area); + + crNorm = sqrt(cr(1,:).^2+cr(2,:).^2+cr(3,:).^2); + zMean = (p(3,t(1,:))+p(3,t(2,:))+p(3,t(3,:)))/3; + nz = -cr(3,:)./crNorm;% z component of normal for each triangle + + volume = area.*zMean.*nz; % contribution of each triangle + totalVolume = sum(volume);%divergence theorem + + end + + function b = inside(shape, radius, theta, phi) + % Determine if spherical point point is inside shape + % + % b = inside(shape, radius, theta, phi) determine if the + % point described by radius, theta (polar), phi (azimuthal) + % is inside the shape. + + [x, y, z] = ott.utils.rtp2xyz(radius(:), theta(:), phi(:)); + b = shape.insideXyz(x, y, z); + 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 + + % Using a third-party function for insidexyz + b = ott.utils.inpolyhedron(shape.faces.', shape.verts.', ... + [x(:),y(:),z(:)]); + 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. + + p = inputParser; + p.addParameter('surfoptions', {}); + p.parse(varargin{:}); + + trisurf(shape.faces.', shape.verts(1, :), ... + shape.verts(2, :), shape.verts(3, :), p.Results.surfoptions{:}); + end + + function writeWavefrontObj(shape, filename) + % Write internal representation of shape to Wavefront OBJ file + % + % writeWavefrontObj(filename) writes the shape to the given file. + shape.writeWavefrontObj_helper(filename, shape.verts, shape.faces); + end + end +end diff --git a/+ott/+shapes/WavefrontObj.m b/+ott/+shapes/WavefrontObj.m new file mode 100644 index 00000000..fe9df8c4 --- /dev/null +++ b/+ott/+shapes/WavefrontObj.m @@ -0,0 +1,84 @@ +classdef WavefrontObj < ott.shapes.TriangularMesh +% WavefrontObj load a shape from a Wavefront OBJ file +% +% The file format is described on the Wikipedia page. +% +% This file is part of the optical tweezers toolbox. +% See LICENSE.md for information about using/distributing this file. + + properties (SetAccess=protected) + filename % Name of the file this object loaded + end + + methods + function shape = WavefrontObj(filename) + % Construct a new shape from a Wavefront OBJ file + % + % WavefrontObj(filename) loads the face and vertex information + % contained in the file. Faces are converted to triangles. + + assert(nargin == 1, 'Filename not supplied'); + + fp = fopen(filename, 'r'); + + verts = zeros(3, 0); + faces = zeros(3, 0); + + % Read lines of file + fline = fgets(fp); + while fline ~= -1 + + switch fline(1:2) + case 'v ' + + % Vert lines should have 3 or more coordinates: x y z [w] + v = sscanf(fline(3:end), '%f'); + + if length(v) < 3 + warning('Ignoring vertex line with insufficient coords'); + else + verts(:, end+1) = v(1:3); + end + + case 'f ' + + % Read face indices, they can have multiple formats: + % f 1 2 3 + % f 1/2 2/2 3/2 + % f 1//2 2//2 3//2 + f1 = sscanf(fline(3:end), '%f'); + f2 = sscanf(fline(3:end), '%f%*[^ ]'); + + % Determine which pattern matched + if length(f1) > length(f2) + f = f1; + else + f = f2; + end + + if length(f) < 3 + warning('Ignoring face with fewer than 3 verts'); + else + for ii = 3:length(f) + faces(:, end+1) = [f(1); f(ii-1); f(ii)]; + end + end + + otherwise + % Nothing to do, ignore line + end + + % Get next line of file + fline = fgets(fp); + end + + % Close file + fclose(fp); + + % Setup object + shape = shape@ott.shapes.TriangularMesh(verts, faces); + shape.filename = filename; + + end + end +end diff --git a/tests/shapes/cube.stl b/tests/shapes/cube.stl new file mode 100644 index 00000000..7d3adf4e Binary files /dev/null and b/tests/shapes/cube.stl differ diff --git a/tests/shapes/testShape.m b/tests/shapes/testShape.m new file mode 100644 index 00000000..24b2078d --- /dev/null +++ b/tests/shapes/testShape.m @@ -0,0 +1,24 @@ +function tests = shape + tests = functiontests(localfunctions); +end + +function setupOnce(testCase) + addpath('../../'); +end + +function testSimple(testCase) + + import matlab.unittest.constraints.IsEqualTo; + import matlab.unittest.constraints.IsOfClass; + import matlab.unittest.constraints.RelativeTolerance; + tol = 1.0e-3; + + r = 1.0; + shape = ott.shapes.Shape.simple('sphere', r); + testCase.verifyThat(shape, IsOfClass(?ott.shapes.Sphere), ... + 'Incorrect shape type'); + testCase.verifyThat(shape.perimiter, IsEqualTo(2*pi*r, ... + 'Within', RelativeTolerance(tol)), ... + 'Incorrect sphere radius'); + +end diff --git a/tests/shapes/testShapeSurface.m b/tests/shapes/testShapeSurface.m new file mode 100644 index 00000000..a152d149 --- /dev/null +++ b/tests/shapes/testShapeSurface.m @@ -0,0 +1,172 @@ +function tests = shapesurface + tests = functiontests(localfunctions); +end + +function setupOnce(testCase) + addpath('../../'); +end + +function testValues(testCase) + + import ott.shapes.Shape; + + import matlab.unittest.constraints.IsEqualTo; + import matlab.unittest.constraints.AbsoluteTolerance; + tol = 1.0e-4; + + % Create a grid of test points + [theta, phi] = ott.utils.angulargrid(3, 3); + sz = numel(theta); + + % Test a sphere with radius 1 + shape = Shape.simple('sphere', 1.0); + r = shape.radii(theta, phi); + n = shape.normals(theta, phi); + [~, ~, rotsym] = shape.axialSymmetry(); + testCase.verifyThat(r, IsEqualTo(repmat(1.0, sz, 1), ... + 'Within', AbsoluteTolerance(tol)), ... + 'Sphere radius should be 1.0'); + testCase.verifyThat(n, IsEqualTo(repmat([1.0, 0, 0], sz, 1), ... + 'Within', AbsoluteTolerance(tol)), ... + 'Sphere normals should be [1 0 0]'); + testCase.verifyThat(rotsym, IsEqualTo(0), ... + 'Sphere rotational symmetry incorrect'); + + % Test a ellipsoid + shape = Shape.simple('ellipsoid', [ 1, 2, 3 ]); + r = shape.radii(theta, phi); + n = shape.normals(theta, phi); + [~, ~, rotsym] = shape.axialSymmetry(); + rTarget = [ 1.7321 1.0000 1.7321 2.2780 1.5119, ... + 2.2780 2.2780 1.5119 2.2780 ].'; + nTarget = [ + 0.6547 0.7559 0 + 1.0000 0.0000 0 + 0.6547 -0.7559 0 + 0.6670 0.4892 0.5620 + 0.8030 0.0000 0.5960 + 0.6670 -0.4892 0.5620 + 0.6670 0.4892 -0.5620 + 0.8030 0.0000 -0.5960 + 0.6670 -0.4892 -0.5620 ]; + testCase.verifyThat(r, IsEqualTo(rTarget, ... + 'Within', AbsoluteTolerance(tol)), ... + 'Ellipsoid radii incorrect'); + testCase.verifyThat(n, IsEqualTo(nTarget, ... + 'Within', AbsoluteTolerance(tol)), ... + 'Ellipsoid normals incorrect'); + testCase.verifyThat(rotsym, IsEqualTo(2), ... + 'Ellipsoid rotational symmetry incorrect'); + + % Test a cylinder + shape = Shape.simple('cylinder', [ 1, 1 ]); + r = shape.radii(theta, phi); + n = shape.normals(theta, phi); + [~, ~, rotsym] = shape.axialSymmetry(); + rTarget = [ 0.5774 1.0000 0.5774 0.5774 1.0000, ... + 0.5774 0.5774 1.0000 0.5774 ].'; + nTarget = [ + 0.8660 -0.5000 0 + 1.0000 0.0000 0 + 0.8660 0.5000 0 + 0.8660 -0.5000 0 + 1.0000 0.0000 0 + 0.8660 0.5000 0 + 0.8660 -0.5000 0 + 1.0000 0.0000 0 + 0.8660 0.5000 0 ]; + testCase.verifyThat(r, IsEqualTo(rTarget, ... + 'Within', AbsoluteTolerance(tol)), ... + 'Cylinder radii incorrect'); + testCase.verifyThat(n, IsEqualTo(nTarget, ... + 'Within', AbsoluteTolerance(tol)), ... + 'Cylinder normals incorrect'); + testCase.verifyThat(rotsym, IsEqualTo(0), ... + 'Cylinder rotational symmetry incorrect'); + + % Test a cylinder (x-axis aligned) + % Doesn't work, not actually supported + %[r, n, rotsym] = shapesurface(theta, phi, -1, [ 1, 1 ]); + %disp(r); + %disp(n); + %disp(rotsym); + + % Test a superellipsoid + shape = Shape.simple('superellipsoid', [ 0.2, 0.5, 1, 1, 1 ]); + r = shape.radii(theta, phi); + n = shape.normals(theta, phi); + [~, ~, rotsym] = shape.axialSymmetry(); + rTarget = [0.3780 0.2000 0.3780 0.5714 0.3288, ... + 0.5714 0.5714 0.3288 0.5714 ].'; + nTarget = [ + 0.5587 0.8294 0 + 1.0000 0.0000 0 + 0.5587 -0.8294 0 + 0.4680 0.5460 0.6949 + 0.7131 0.0000 0.7010 + 0.4680 -0.5460 0.6949 + 0.4680 0.5460 -0.6949 + 0.7131 0.0000 -0.7010 + 0.4680 -0.5460 -0.6949 ]; + testCase.verifyThat(r, IsEqualTo(rTarget, ... + 'Within', AbsoluteTolerance(tol)), ... + 'Superellipsoid radii incorrect'); + testCase.verifyThat(n, IsEqualTo(nTarget, ... + 'Within', AbsoluteTolerance(tol)), ... + 'Superellipsoid normals incorrect'); + testCase.verifyThat(rotsym, IsEqualTo(1), ... + 'Superellipsoid rotational symmetry incorrect'); + + % Test a cone-tipped cylinder + shape = Shape.simple('cone-tipped-cylinder', [ 1, 1, 1 ]); + r = shape.radii(theta, phi); + n = shape.normals(theta, phi); + [~, ~, rotsym] = shape.axialSymmetry(); + rTarget = [1.0981 1.0000 1.0981 1.0981 1.0000, ... + 1.0981 1.0981 1.0000 1.0981 ].'; + nTarget = [ + 0.9659 0.2588 0 + 1.0000 0.0000 0 + 0.9659 -0.2588 0 + 0.9659 0.2588 0 + 1.0000 0.0000 0 + 0.9659 -0.2588 0 + 0.9659 0.2588 0 + 1.0000 0.0000 0 + 0.9659 -0.2588 0 ]; + testCase.verifyThat(r, IsEqualTo(rTarget, ... + 'Within', AbsoluteTolerance(tol)), ... + 'Cone-tipped cylinder radii incorrect'); + testCase.verifyThat(n, IsEqualTo(nTarget, ... + 'Within', AbsoluteTolerance(tol)), ... + 'Cone-tipped cylinder normals incorrect'); + testCase.verifyThat(rotsym, IsEqualTo(0), ... + 'Cone-tipped cylinder rotational symmetry incorrect'); + + % Test a cube + shape = Shape.simple('cube', 1); + r = shape.radii(theta, phi); + n = shape.normals(theta, phi); + [~, ~, rotsym] = shape.axialSymmetry(); + rTarget = [0.5774; 0.5000; 0.5774; 0.5774; 0.5774; ... + 0.5774; 0.5774; 0.5774; 0.5774]; + testCase.verifyThat(r, IsEqualTo(rTarget, ... + 'Within', AbsoluteTolerance(tol)), ... + 'Cube radii incorrect'); + nTarget = [ + 0.8660 -0.5000 0 + 1.0000 0.0000 0 + 0.8660 0.5000 0 + 0.8660 -0.5000 0 + 0.8660 0.0000 -0.5000 + 0.8660 0.5000 0 + 0.8660 -0.5000 0 + 0.8660 0.0000 0.5000 + 0.8660 0.5000 0 ]; + testCase.verifyThat(n, IsEqualTo(nTarget, ... + 'Within', AbsoluteTolerance(tol)), ... + 'Cube normals incorrect'); + testCase.verifyThat(rotsym, IsEqualTo(4), ... + 'Cube rotational symmetry incorrect'); + +end diff --git a/tests/shapes/testStlLoader.m b/tests/shapes/testStlLoader.m new file mode 100644 index 00000000..28c0f1ae --- /dev/null +++ b/tests/shapes/testStlLoader.m @@ -0,0 +1,22 @@ +function tests = stlloader + tests = functiontests(localfunctions); +end + +function setupOnce(testCase) + addpath('../../'); +end + +function testLoad(testCase) + + % Load the shape + shape = ott.shapes.StlLoader('cube.stl'); + + import matlab.unittest.constraints.IsEqualTo; + + testCase.verifyThat(numel(shape.verts), IsEqualTo(3*8), ... + 'Incorrect number of vertices'); + + testCase.verifyThat(numel(shape.faces), IsEqualTo(3*2*6), ... + 'Incorrect number of faces'); + +end diff --git a/tests/shapes/testWavefrontObj.m b/tests/shapes/testWavefrontObj.m new file mode 100644 index 00000000..a9b2e42c --- /dev/null +++ b/tests/shapes/testWavefrontObj.m @@ -0,0 +1,34 @@ +function tests = wavefrontobj + tests = functiontests(localfunctions); +end + +function setupOnce(testCase) + addpath('../../'); +end + +function testLoad(testCase) + + fn = tempname(); + + % Generate a cube and write it to a file + cube = ott.shapes.Cube(1.0); + cube.writeWavefrontObj(fn); + + % Load the shape + shape = ott.shapes.WavefrontObj(fn); + + import matlab.unittest.constraints.IsEqualTo; + import matlab.unittest.constraints.RelativeTolerance; + tol = 1.0e-3; + + testCase.verifyThat(shape.maxRadius, IsEqualTo(cube.maxRadius, ... + 'Within', RelativeTolerance(tol)), ... + 'Incorrect radius after writing/loading OBJ file'); + + testCase.verifyThat(shape.volume, IsEqualTo(cube.volume, ... + 'Within', RelativeTolerance(tol)), ... + 'Incorrect volume after writing/loading OBJ file'); + + delete(fn); + +end diff --git a/thirdparty/README.md b/thirdparty/README.md new file mode 100644 index 00000000..121723e6 --- /dev/null +++ b/thirdparty/README.md @@ -0,0 +1,6 @@ + +This directory contains licenses for third party code that has +been used as part of this toolbox. Functions that contain +third party code have a reference to the license file in this +directory. + diff --git a/thirdparty/stl_EricJohnson.txt b/thirdparty/stl_EricJohnson.txt new file mode 100644 index 00000000..09b5ba8e --- /dev/null +++ b/thirdparty/stl_EricJohnson.txt @@ -0,0 +1,27 @@ +Copyright (c) 2011, Eric Johnson +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in + the documentation and/or other materials provided with the distribution + * Neither the name of the The MathWorks, Inc. nor the names + of its contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE. diff --git a/thirdparty/stl_KrishnanSuresh.txt b/thirdparty/stl_KrishnanSuresh.txt new file mode 100644 index 00000000..73e8d45a --- /dev/null +++ b/thirdparty/stl_KrishnanSuresh.txt @@ -0,0 +1,24 @@ +Copyright (c) 2010, Krishnan Suresh +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in + the documentation and/or other materials provided with the distribution + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE. diff --git a/thirdparty/stl_Sven.txt b/thirdparty/stl_Sven.txt new file mode 100644 index 00000000..4e65529b --- /dev/null +++ b/thirdparty/stl_Sven.txt @@ -0,0 +1,24 @@ +Copyright (c) 2015, Sven +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in + the documentation and/or other materials provided with the distribution + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE.