% M. Sukop and E. Perfect, University of Kentucky, 2000
% Program computes, displays, and stores randomized Menger Sponges. 
% The homogeneous or heterogeneous (binomial/truncated binomial) algorithms 
% can be specified. 
% b is fractal scaling parameter (3 for standard Menger Sponge).
% p is the probability of a pore at each iteration level (7/27 for standard sponge).
% Output to bmp files by layer.

clear
rand('state',sum(100*clock))

b=input('b value ? ')
p=input('p value ? ')

n=round(b^3-p*b^3)

D=log(n)/log(b)

maxit=input('Maximum number of iterations ? ')

minpore=input('Do you want a mimimum pore size of b voxels (enter 1) or 1 voxel (enter 0) ? ')

typemode=input('Select fractal type: 0 = homogeneous; 1 = heterogeneous; 2 = non-random ')

if typemode == 1
    
    hettypemode=input('Select heterogeneous mode: 0 = regular binomial; 1 = truncated binomial ')
    
elseif typemode == 2
    
    a=[input('Specify generator;    e.g., [1 1 1, 1 0 1, 1 1 1, 1 0 1, 0 0 0, 1 0 1, 1 1 1, 1 0 1, 1 1 1] for standard Menger sponge')]
    
end


it=1
if typemode==0 %Homogeneous Algorithm
    
    a=randperm(b^3);
    for i = 1:b^3;
        if a(i)<=n 
            c(i)=1;%solid
        else 
            c(i)=0;%pore
        end
    end
elseif typemode==1 %Heterogeneous Algorithm
    
    if hettypemode==0 %Regular binomial
        
        a=rand(b^3,1);
        for i = 1:b^3;
            if a(i)<=1-p 
                c(i)=1;
            else 
                c(i)=0;
            end
        end
        
    else %Truncated binomial
        a=2;
        while min(a)>1-p
            a=rand(b^3,1);
        end
        for i = 1:b^3;
            if a(i)<=1-p 
                c(i)=1;
            else 
                c(i)=0;
            end
        end
    end
    
else %non-random
    c=a;
end %end algorithm selection


for row=1:b;
    for col=1:b;
        for lay=1:b;
            d(row,col,lay)=c(row*b^2-(b^2-b*col)-(b-lay));
        end;
    end;
end;
oldcol=1

while it<maxit
    for oldlay=1:b^it;
        disp(100*oldlay*oldcol/b^(2*it))
        disp('% done current iteration')
        for oldcol=1:b^it;
            for oldrow=1:b^it;
                if d(oldrow,oldcol,oldlay)==0;%if pore 
                    for newlay=oldlay*b-(b-1):oldlay*b;
                        for newcol=oldcol*b-(b-1):oldcol*b;
                            for newrow=oldrow*b-(b-1):oldrow*b;
                                e(newrow,newcol,newlay)=0;%d(oldrow,oldcol,oldlay);
                            end   
                        end
                    end  
                elseif (minpore&(it==maxit-1)) %if minimum pore size desired
                    if d(oldrow,oldcol,oldlay)==1%if solid 
                        for newlay=oldlay*b-(b-1):oldlay*b;
                            for newcol=oldcol*b-(b-1):oldcol*b;
                                for newrow=oldrow*b-(b-1):oldrow*b;
                                    e(newrow,newcol,newlay)=1;%d(oldrow, oldcol,oldlay);
                                end
                            end
                        end
                    end
                else
                    
                    if typemode==0 %Homogeneous Algorithm
                        a=randperm(b^3);
                        for i = 1:b^3;
                            if a(i)<=n 
                                c(i)=1;
                            else 
                                c(i)=0;
                            end
                        end
                    elseif typemode==1 %Heterogeneous Algorithm
                        
                        if hettypemode==0 %Regular binomial
                            
                            a=rand(b^3,1);
                            for i = 1:b^3;
                                if a(i)<=1-p 
                                    c(i)=1;
                                else 
                                    c(i)=0;
                                end
                            end
                            
                        else %Truncated binomial
                            a=2;
                            while min(a)>1-p
                                a=rand(b^3,1);
                            end
                            for i = 1:b^3;
                                if a(i)<=1-p 
                                    c(i)=1;
                                else 
                                    c(i)=0;
                                end
                            end
                        end
                    else %non-random
                        c=a;
                    end %end algorithm selection
                    
                    newrow=oldrow*b-(b-1);
                    newcol=oldcol*b-(b-1);
                    newlay=oldlay*b-(b-1);
                    for row=1:b;
                        for col=1:b;
                            for lay=1:b;
                                e(newrow+row-1,newcol+col-1,newlay+lay-1)=c(row*b^2-(b^2-b*col)-(b-lay));
                            end   
                        end
                    end
                end
            end   
        end
        if it==maxit-1
            for newlay=oldlay*b-(b-1):oldlay*b;
                % Measure D of slice
                Dslice(newlay)=log(sum(sum(e(:,:,newlay))))/log(b^(it+1));
                
                imwrite(e(:,:,newlay),sprintf('%s%04d%s','newlay',newlay,'.bmp'),'bmp');
            end
        end
    end

    it=it+1
    d=e;
end
plot(Dslice)
%close


% figure(1)
% hold on
% fid = fopen('junk.dat','w');
% col=0;
% for row=1:b^maxit;
%     disp(100*row*col/b^(2*it))
%     disp('% done generating plot')
%     for col=1:b^maxit;
%         for lay=1:b^maxit;
%             if d(row,col,lay)==1;
%                 xmin=(row-1)/b^maxit;
%                 xmax=row/b^maxit;
%                 ymin=(col-1)/b^maxit;
%                 ymax=col/b^maxit;
%                 zmin=lay/b^maxit;
%                 zmax=(lay-1)/b^maxit; 
%                 fprintf(fid,'%6.5f %6.5f %6.5f \n',xmin+(xmax-xmin)/2,ymin+(ymax-ymin)/2,zmin+(zmax-zmin)/2); vert=[xmin,ymin,zmin;xmax,ymin,zmin; xmax,ymax,zmin; xmin,ymax,zmin; xmin,ymin,zmax; xmax,ymin,zmax; xmax,ymax,zmax; xmin,ymax,zmax];
%                 fac=[1 2 6 5; 2 3 7 6; 3 4 8 7; 4 1 5 8; 1 2 3 4; 5 6 7 8];
%                 rgbcolors=[.5 .5 .5 ; 0 0 0; .5 .5 .5 ; 0 0 0; 1 1 1; 1 1 1];
%                 %rgbcolors=[1 0 0 ; 1 0 0 ; 1 0 0  ; 1 0 0; 1 0 0 ; 1 0 0];
%                 
%                 patch('faces',fac,'vertices',vert,'FaceVertexCData',rgbcolors,'FaceColor','flat','EdgeColor',[0 0 0],'BackFaceLighting', 'reverselit');
%                 view(3)
%             end;   
%         end;
%     end;
% end;    
%     fclose(fid)
%     light('Position',[-1,0,0]);
%     light('Position',[0,-1,0]);
%     light('Position',[0,0,1]);
%     
%     hold off