0
votes

I'm trying to output Voronoi cells to a format which can be used for 3D printing.

MATLAB generates voronoi cells from lists of X and Y coordinates. My script generates such a list of points, but getting to a format I can export is seeming problematic.

My main hopes lie with stlwrite, http://www.mathworks.com/matlabcentral/fileexchange/20922-stlwrite-write-binary-or-ascii-stl-file

This function/script requires a surface to export.

function [lines, LineData, pOut] = makeSurfaceFromVorVerts(Lx, Ly, varargin)
    p = inputParser;
    addRequired(p,'Lx',...
        @(x) (isequal(class(x),'double') && isequal(size(x,1),2)));
    addRequired(p,'Ly',...
        @(x) (isequal(class(x),'double') && isequal(size(x,1),2)));
    defaultResolution = 100;
    addOptional(p,'Resolution',defaultResolution,...
        @(x) (isequal(class(x),'double') && isequal(size(x),[1 1])));
    defaultBoundary = [0 110; 0 110];
    addOptional(p,'Boundaries',defaultBoundary,...
        @(x) (isequal(class(x),'double') && isequal(size(x),[2 2])));
    parse(p,Lx,Ly,varargin{:});
    pOut = p;

    LX = p.Results.Lx;
    LY = p.Results.Ly;
    Bounds = p.Results.Boundaries;

    % Strip high values
    reducedXdat = [];
    reducedYdat = [];
    for i=1:size(LX,2)
        if LX(1,i) > Bounds(1,1) && LX(1,i) < Bounds(1,2) && ... % X value of start of line
           LX(2,i) > Bounds(2,1) && LX(2,i) < Bounds(2,2) && ... % Y value of start of line
           LY(1,i) > Bounds(1,1) && LY(1,i) < Bounds(1,2) && ... % X value of end of line
           LY(2,i) > Bounds(2,1) && LY(2,i) < Bounds(2,2),       % Y value of end of line
            reducedXdat = [reducedXdat, LX(:,i)];
            reducedYdat = [reducedYdat, LY(:,i)];
        end
    end

    % Initialise a grid of points
    %sXnew = (Bounds(1,2) - Bounds(1,1)) * p.Results.Resolution;
    %sYnew = (Bounds(2,2) - Bounds(2,1)) * p.Results.Resolution;
    %Z = zeros(sXnew, sYnew,'uint8');


    %for x=1:size(X,1)
    %    for y=1:size(Y,1)
    %        nX = floor(X(x)*p.Results.Resolution);
    %        nY = floor(Y(y)*p.Results.Resolution);
    %        Z(nX,nY) = 1;
    %    end
    %end
    %surface = Z;
    %coords = [X,Y];
    lines = line(reducedXdat,reducedYdat);
    LineData = [reducedXdat; reducedYdat];
end

My processing script, above, takes the points generated by the command

[Lx, Ly] = voronoi(xValuesOfCellCentres, yValuesOfCellCentres);

along with an optional 'boundary' matrix (there's also a check for resolution, for the commented section) and then outputs lines.

I would like these lines to form the surface. I considered creating a mesh using a binary Z value (1 for points, 0 for everywhere else), but I don't know how I could also include the positions between points, ie those covered by the lines.

I expect that there is some relevant middle step I can take, to create a frame based on extrusion of the lines drawn (either by this script, which has cut away the extra lines to infinity, or by voronoi(X,Y), but I can't work it out.

1

1 Answers

0
votes

Found a way to do this.

Changed the script, new script is pasted at the bottom.

Workflow goes:

  1. Matlab (create line plot, then for i=1:size(lineHandles,1), set(lineHandles(i),'lineWidth',4),end). Change the line width as required.)
  2. Save as .png file.
  3. Gimp (crop file to a nice rectangle)
  4. InkScape (open png file, click image, menu Path -> Trace Bitmap -> Color Quantization (2 colors) -> Drag the path aside (it shows path in the status bar at the bottom when you have this selected). Click image and delete. Place path (with nodes) back onto the page (there's an x/y coordinates setting in a toolbar at the top, set these both to 0.
  5. Save as a .dxf file. To get these to input more nicely, use the extension provided at http://www.thingiverse.com/thing:14221. This should be installed by copying the .py and .inx files to /usr/share/inkscape/extensions (or similar)
  6. Open in OpenSCAD using the command import("/path/to/filename.dxf");
  7. This object can be extruded using the linear_extrude(height = x) where x is a height in mm (other lengths probably configurable)
  8. Render using CGAL (F6 is the shortcut for this in OpenSCAD.)
  9. Export to .stl file (menu Design > Export as STL...)

MATLAB script (edit as needed and take outputs as required, if you want the line handles you'll need to put in almost all of the output arguments (or reorder them):

%voronoiScriptHex generates voronoi cells in a hexagonal tesselation grid
%   [X, Y, Fig, Axis, Lx, Ly, lH, lD] = voronoiScript(Bnd, Spc, D, ...)
%   
%   Output variables
%       X is the x coordinate of the voronoi cell centres
%       Y is the y coordinate of the voronoi cell centres
%       Fig is the handle of the figure generated
%       Axis is the handle of the axes used
%       Lx gives the start and end x coordinates of the voronoi lines
%       Ly gives the start and end y coordinates of the voronoi lines
%       lH is the set of handles for the voronoi lines
%       lD is constructed from [Lx; Ly], it is a [4 by Length] array
%
%   Bnd specifies the boundaries for the region to be covered. It should be
%   either one number, or a pair of numbers in a [1x2] vector, specifying
%   [maxX, maxY]. 0 is taken as the minimum value.
%
%   Spc specifies the average spacing. For a hex grid, it only accepts one
%   value.
%
%   D specifies the variation from a uniform grid. It is multiplied by a 
%   random number between -0.5 and 0.5 and added to the [x,y] coordinates
%   of a point. If size(D) = [1x2], then the values are for [Dx, Dy].
%
%   Optional arguments can be used to place some points exactly on the grid
%   they would lie on with D[x/y] = 0. The first should be 'PartFixed' -
%   this is a boolean and if true, some points are fixed to the grid.
%
%   The second argument is 'FractionFixed'. This is an integer value
%   (double class variables are accepted, but floor(arg) must be equal to
%   (arg)). It specifies inversely how often points should be fixed, eg a
%   value of 1 fixes every point, whilst a value of 5 fixes 1/5 of the
%   points.
%
%   PlotScatter is another boolean value, which sets if a scatter plot of
%   the X,Y values corresponding to the cell centres should be included in
%   the figure.

function [X, Y, Figure, Axis, Lx, Ly, lineHandles, lineData] = ...
    voronoiScriptHex(Boundary, Spacing, Delta, varargin)

    p = inputParser;
    p.FunctionName = 'voronoiScript';
    addRequired(p, 'Boundary', @checkTypes);
    addRequired(p, 'Spacing', @isnumeric);
    addRequired(p, 'Delta', @checkTypes);
    defaultPart = false;
    addOptional(p, 'PartFixed', defaultPart, @islogical);
    defaultFraction = 2;
    addOptional(p, 'FractionFixed', defaultFraction, @isAnInt);
    defaultScatter = false;
    addOptional(p, 'PlotScatter', defaultScatter, @islogical);
    parse(p, Boundary, Spacing, Delta, varargin{:});


    % Get values for boundaries and delta
    % (Can be vectors or scalars)
    if isequal(size(p.Results.Boundary),[1,2])
        % Boundary is a vector [maxX, maxY]
        BoundaryY = p.Results.Boundary(1,2);
    else
        BoundaryY = p.Results.Boundary(1,1);
    end
    if isequal(size(p.Results.Delta),[1,2])
        % Delta is a vector [dX, dY]
        DeltaY = p.Results.Delta(1,2);
    else
        DeltaY = p.Results.Delta(1,1);
    end

    Spacing = p.Results.Spacing;
    BoundaryX = p.Results.Boundary(1,1);
    DeltaX = p.Results.Delta(1,1);
    D1 = [2*Spacing*cosd(30), Spacing];
    numP = [ceil(BoundaryX/D1(1,1)) ceil(BoundaryY/D1(1,2))];

    D2 = D1 ./ 2;

    % Create the values
    counter = 1;
    xList(numP(1,1)*numP(1,2)) = 0;
    yList(numP(1,1)*numP(1,2)) = 0;
    for x=1:numP(1,1)
        for y = 1:numP(1,2)
            xList(counter) = (getPointValue(x, D1(1,1), DeltaX)-D2(1,1));
            xList(counter+1) = getPointValue(x, D1(1,1), DeltaX);
            yList(counter) = (getPointValue(y, D1(1,2), DeltaY)-D2(1,2));
            yList(counter+1) = getPointValue(y, D1(1,2), DeltaY);
            counter = counter + 2;
        end
    end

    % Set some of the points to be without random change
    if (p.Results.PartFixed),
        for counter=1:p.Results.FractionFixed:size(xList,2),
            [x, y] = getXYfromC(counter, numP(1,2));
            xList(counter) = x*Spacing;
            yList(counter) = y*Spacing;
        end
    end

    X = xList;
    Y = yList;

    % Set manual ticks for the figure axes
    ticksX = zeros(1,numP(1,1)+1);
    for i=1:numP(1,1)+1,
        ticksX(i) = i*D1(1,1);
    end
    ticksY = zeros(1,numP(1,2)+1);
    for i=1:numP(1,2)+1,
        ticksY(i) = i*D1(1,2);
    end

    BoundCoeff = 1.08;
    Bounds = [0 BoundCoeff*BoundaryX; 0 BoundCoeff*BoundaryY];

    % Give the figure a handle that is returned, and set axes values
    Figure = figure;
    Axis = axes;
    axis equal;
    minor = 'off';
    gridtoggle = 'off';
    set(Axis,'XTickMode','manual','YTickMode','manual', ...
        'XGrid',gridtoggle,'YGrid',gridtoggle, ...
        'XMinorGrid',minor,'YMinorGrid',minor, ...
        'XTick',ticksX,'YTick',ticksY, ...
        'XMinorTick',minor,'YMinorTick',minor, ...
        'XLim',[0 Bounds(1,2)],'YLim',[0 Bounds(2,2)]);

    %set(Axis,'XLim',[0 Bounds(1,2)],'YLim',[0 Bounds(2,2)]);

    % Create the voronoi cells, returning the line points
    [Lx, Ly] = voronoi(X,Y);

    % Strip high values
    counter = 1;
    reducedXdat = zeros(2,size(Lx,2));
    reducedYdat = zeros(2,size(Lx,2));
    for i=1:size(Lx,2)
        if Lx(1,i) > Bounds(1,1) && Lx(1,i) < Bounds(1,2) && ... % X value of start of line
           Lx(2,i) > Bounds(2,1) && Lx(2,i) < Bounds(2,2) && ... % Y value of start of line
           Ly(1,i) > Bounds(1,1) && Ly(1,i) < Bounds(1,2) && ... % X value of end of line
           Ly(2,i) > Bounds(2,1) && Ly(2,i) < Bounds(2,2),       % Y value of end of line
            reducedXdat(:,counter) = Lx(:,i);
            reducedYdat(:,counter) = Ly(:,i);
            counter = counter + 1;
        end
    end
    Lx = reducedXdat(:,1:counter-1);
    Ly = reducedYdat(:,1:counter-1);

    % Plot the voronoi lines
    lineHandles = line(Lx, Ly);

    % Set colours to black
    if (1)
        for i=1:size(lineHandles,1)
            set(lineHandles(i),...
                'LineWidth',3, ...
                'LineSmoothing','on', ...
                'Color',[0 0 0]);
        end
    end

    lineData = [Lx; Ly];
    if (p.Results.PlotScatter)
        hold on;
        scatter(X,Y);
    end
end

function bool = checkTypes(arg)
    bool = (isequal(class(arg),'double') && ...
        (isequal(size(arg),[1,1]) || isequal(size(arg),[1,2])));
end

function bool = isAnInt(arg)
    bool = (isequal(floor(arg),arg) && ...
        isequal(size(arg),[1,1]) && ...
        isnumeric(arg));
end

function val = getPointValue(intV, spacing, delta)
    val = (((rand(1)-0.5)*delta)+(intV*spacing));
end

function [x,y] = getXYfromC(counter, sizeY)
    x = floor(counter/sizeY)+1;
    y = counter - ((x-1)*sizeY);
end

SCAD script file to place the voronoi cells within a pair of cylinders for 3D printing a grid. Edit as needed with path, or changes to shapes etc:

$fn = 360;
inkFile = "/path/to/my/file";
scl = 1.05; // Scale variable
transVec = [-100, -98, 0]; // Move the imported image pattern so that it's centred.
union(){
    makeRing([0,0,3],inR=96, outR=99, height=6);
    makeRing([0,0,1.5], inR=99, outR=102, height=3);

    intersection() {
        #linear_extrude(height = 6.05, convexity=40, center=false) // Removing the # will get rid of the overlay, which helps see where the grid is. 
            scale([scl, scl, 1])
            translate(transVec)
            import(inkFile);
        makeCylinder([0,0,3], 96, 6, $fn=360);
    }
}
module makeCylinder(centre, radius, height) {
    translate(centre){
        cylinder(h = height, r = radius, center=true);
    }
}
module makeRing(centre,inR, outR, height) {
    difference() {
        makeCylinder(centre, outR, height);
        makeCylinder(centre, inR, height+0.1);
    }
}