Julia Vase Generator

%Filip Jander  8/24/2019
%Exports stl from parametric equation, use a program called Meshmixer to
%clean up the stl file for 3d printing. Recommended functions in meshmixer:
%Reduce, Erase and Fill, Smooth and Make Solid

clear
clc
pointDensity = 100;
outputArray = zeros(pointDensity);
iterations = 100;

maxT = 200;
threeDoutput = NaN(pointDensity^2*maxT,3);
% (t/maxT) + sin((t/maxT)*2*pi)/10)
for t = 0:maxT
    clearvars -except pointDensity iterations maxT threeDoutput t
    juliaPoint = 0.32*(sin((t/maxT)*2*pi)) + .3*(cos((t/maxT)*2*pi)-1)*1i;
    juliaPointPlot(t+1) = juliaPoint;
    for m = 1:pointDensity
        for n = 1:pointDensity
            k = 0;
            x = m/((pointDensity-1)/4)-2-4/pointDensity;
            y = n/((pointDensity-1)/4)-2-4/pointDensity;
            point =  x + y*1i;
            zn = point;

            while k < iterations
                if sqrt(real(zn)^2+ imag(zn)^2) > 2
                    k = iterations + 1;
                else
                    k = k + 1;
                    zn = zn^2 + juliaPoint;
                end
            end
            if k == iterations + 1
                outputArray(m,n) = 1;
            end
        end
    end
    clear plotArray count
    count = 1;
    for n = 1:size(outputArray,1)
        for m = 1:size(outputArray,2)
            if outputArray(n,m) == 0
                %if t ~= 1 && t ~= maxT

                if 1==1
                    if outputArray(n+1,m) == 1 || outputArray(n-1,m) == 1 || outputArray(n,m+1) == 1 || outputArray(n,m-1) == 1
                        plotArray(count,1) = m;
                        plotArray(count,2) = n;
                        count = count + 1;
                    end
                else
                    plotArray(count,1) = m;
                    plotArray(count,2) = n;
                    count = count + 1;
                end
            end
        end
    end
    scatter(plotArray(:,1),plotArray(:,2),10,"filled")
    axis equal
    pause(0.005)

    threeDoutput((t+1)*pointDensity^2:(t+1)*pointDensity^2-1+size(plotArray,1),1:2) = plotArray;
    threeDoutput((t+1)*pointDensity^2:(t+1)*pointDensity^2-1+size(plotArray,1),end) = (t+1);
end

count = 0;
for n = 1:size(threeDoutput)
    if ~isnan(threeDoutput(n,1))
        count = count + 1;
    end
end
threeDoutput2 = zeros(count,3);
count = 1;
for n = 1:size(threeDoutput)
    if ~isnan(threeDoutput(n,1))
        threeDoutput2(count,:) = threeDoutput(n,:);
        count = count +1;
    end
end

threeDoutput2(:,3) = threeDoutput2(:,3)/2.5;
plot(juliaPointPlot);
scatter3(threeDoutput2(:,1),threeDoutput2(:,2),threeDoutput2(:,3),1,"filled")
axis equal
%
% out = [threeDoutput2(:,1),threeDoutput2(:,2),threeDoutput2(:,3)];
% save xyzfile.xyz out -ascii

% k = boundary(threeDoutput2,.8);
% hold on
% trisurf(k,threeDoutput2(:,1),threeDoutput2(:,2),threeDoutput2(:,3),'Facecolor','red','FaceAlpha',0.1)
% axis equal

% xSize = max(threeDoutput2(:,1))+1;
% ySize = max(threeDoutput2(:,2))+1;
% zSize = max(threeDoutput2(:,3))+1;
% threeDoutputLogical = false(xSize,ySize,zSize);
% for n = 1:size(threeDoutput2,1)
%     threeDoutputLogical(round(threeDoutput2(n,1))+1,round(threeDoutput2(n,2))+1,round(threeDoutput2(n,3))+1) = 1;
%     %round(threeDoutput2(n+1,1)),round(threeDoutput2(n+1,2)),round(threeDoutput2(n+1,3))
% end
% fv = isosurface(~threeDoutputLogical, 0.5); % Make patch w. faces "out"
% nfv = reducepatch(fv,.25);
% stlwrite('test5.stl',nfv) % Save to binary .stl