% M. Sukop, 2000. MATLAB V. 5.3
% Program computes, displays, and stores randomized Seirpinski Carpets.
% The homogeneous or heterogeneous (binomial/truncated binomial) algorithms 
% can be specified. 
% b is fractal scaling parameter (3 for standard Seirpinski Carpet).
% p is the probability of a pore at each iteration level (1/9 for standard carpet).
% Output is to file Junk.bmp ('truecolor') or to Junk.dat, 
% which is either ARCINFO-compatible Grid ASCII file 
% or ASCII file of x and y coordinates and 1's (solids) or 0's (pores)
% in the following order:
% row 1, column 1; 
% row 1, column 2; ...

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

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

if typemode < 2

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

end

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

nr(maxit)=0;

if typemode == 1

	hettypemode=input('Select heterogeneous mode: 0 = regular binomial; 1 = truncated binomial ')

elseif typemode == 2
      
   a=input('Specify generator; e.g., [0 0 0, 0 1 0, 0 0 0] for standard Seirpinski carpet: ')
   b=sqrt(length(a))
   p=sum(a)/length(a)
   
end

n=round(b^2-p*b^2)
D=log(n)/log(b)

outmode=input('Select output mode:...0 = GRIDASCII; 1 = X,Y ASCII; 2 = BMP(truecolor) ')

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

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

oldcol=1
while it<maxit
for oldrow=1:b^it;
	disp(100*oldrow*oldcol/b^(2*it))
	disp('% done current iteration')
	for oldcol=1:b^it;
		if d(oldrow,oldcol)==uint8(255)%if pore 
			for newrow=oldrow*b-(b-1):oldrow*b;
				for newcol=oldcol*b-(b-1):oldcol*b;
					e(newrow,newcol)=uint8(d(oldrow, oldcol));
				end
			end
		else				
         
         if typemode==0 %Homogeneous Algorithm
            
				a=randperm(b^2);
				for i = 1:b^2;
					if a(i)<=n 
				 		c(i)=uint8(0);%solid - black 
					else 
						c(i)=uint8(255);%pore - white
					end
			   end
   
			elseif typemode==1 %Heterogeneous Algorithm
   
			   if hettypemode==0 %Regular binomial
      
   				a=rand(b^2,1);
					for i = 1:b^2;
						if a(i)<=1-p 
                     c(i)=uint8(0);
                     nr(it+1)=nr(it+1)+1; %Count solids
						else 
		      			c(i)=uint8(255);
	   	   		end
	      		end
      
   			else %Truncated binomial
   		   	a=2;
		      	while min(a)>1-p
			      	a=rand(b^2,1);
      			end
					for i = 1:b^2;
						if a(i)<=1-p 
	 						c(i)=uint8(0);
						else 
      					c(i)=uint8(255);
      				end
   				end
            end
         else %non-random
            c=uint8(255*a);
			end %end algorithm selection

			newrow=oldrow*b-(b-1);
			newcol=oldcol*b-(b-1);
			for row=1:b;
				for col=1:b;
					e(newrow+row-1,newcol+col-1)=uint8(c(row*b-(b-col)));
				end
			end
		end
	end
end
it=it+1
d=uint8(e);
end
clear e;

if outmode==1
	%generate output arrays
	for i=1:(b^it);
		disp(100*i/b^it)
		disp('% done building x and y vectors')
		for j=1:(b^it);
			x((i-1)*b^it+j)=(2*i-1)/(2*b^it);
			y((i-1)*b^it+j)=1-(2*j-1)/(2*b^it);
		end
	end
	disp('Reshaping and transposing')
   z=double(reshape(d,b^(2*it),1));
   z=z/255;
   z=1-z; %exchange 1 and 0 ('solids' = 1)
	x=x';
	y=y';
	out=[x y z];
	disp('writing file')
   save junk.dat out -ascii;
   
elseif outmode==0
	%generate gridascii output
   z=double(d);
   z=z/255;
   z=1-z; %exchange 1 and 0 ('solids' = 1)

fid=fopen('junk.dat','w');
fprintf(fid,'NCOLS %g',b^it);
fprintf(fid,'\rNROWS %g',b^it);
fprintf(fid,'\rXLLCORNER %g',0);
fprintf(fid,'\rYLLCORNER %g',0);
fprintf(fid,'\rCELLSIZE %g',1/b^it);
fprintf(fid,'\rNODATA_VALUE %g',9999);
for i=1:b^it;
	fprintf(fid,'\r');
		fprintf(fid,'%2.0g',z(i,:));
end
status=fclose(fid);

elseif outmode==2
%bmp1 output
e(:,:,1)=d;
e(:,:,2)=d;
e(:,:,3)=d;   
imwrite(e,'junk.bmp')

viewmode=input('View bmp file?: 0 = No; 1 = Yes ')

	if viewmode==1
		image(imread('junk.bmp'))
   	axis('equal')
	end
end