%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                                                                %%%
%%%      GENESIS - Surface Generation Prototypes, v0.9 (BETA)      %%%
%%%      ----------------------------------------------------      %%%
%%%      Harris Georgiou - mailto:xgeorgiou@yahoo.com              %%%
%%%                                                                %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


clear all;

disp('Terrain Surface Generator, v0.9B - Harris Georgiou (c) 2000');
disp(' ');

N=100;			% terrain size (NxN)
SLR=0.2;		% Sea-Level-Ratio vs size N (0..1)
HBR=0.85;		% Height-to-Base-Ratio for mountains (0..1)
L=16;			% smoothing mask size (LxL)
L2=L/2;			% =L/2 (constant)

disp('generating random peaks plane...');

%... initialize variables ...

randn('seed',sum(100*clock));

v1=0.5*randn;
v2=0.5*randn+0.25*v1;
v3=0.5*randn+0.25*v2+0.17*v1;
v=0.5*randn+0.25*v3+0.17*v2+0.08*v1;

%... generate peaks plane ...
tern=zeros(N,N);
for i=1:N,
    for j=1:N,
	v1=v2;  v2=v3;  v3=v;
	v=0.5*randn+0.25*v3+0.17*v2+0.08*v1;
	tern(i,j)=v*N;
    end;
end;

%... scale to normal height ...
mH=max(max(tern));
mL=min(min(tern));
for i=1:N,
    for j=1:N,
	tern(i,j)=(tern(i,j)-mL)/(mH-mL)*N;
    end;
end;

disp('generating processed torus plane...');

%... peak smoothing (level-1), torus processing ...
smth=zeros(N,N);
for i=1:N,
    for j=1:N,
	sum=0;
	for x=-L2:L2,
	    for y=-L2:L2,
		ii=i+x;
		if (ii<1) ii=ii+N;
		elseif (ii>N) ii=ii-N;
		end;

		jj=j+y;
		if (jj<1) jj=jj+N;
		elseif (jj>N) jj=jj-N;
		end;
		
		sum=sum+tern(ii,jj);
	    end;
	end;
	smth(i,j)=sum/(L*L);
    end;
end;
    
%... add small noise factors ...
for i=1:N,
    for j=1:N,
	v1=v2;  v2=v3;  v3=v;
	v=0.5*randn+0.25*v3+0.17*v2+0.08*v1;
	smth(i,j)=smth(i,j)+v*L/N;
    end;
end;

%... scale to normal height ...
mH=max(max(smth));
mL=min(min(smth));
for i=1:N,
    for j=1:N,
	smth(i,j)=(smth(i,j)-mL)/(mH-mL)*N;
    end;
end;

disp('generating final terrain surface...');

%... surface smoothing (level-2), torus processing ...
smth2=zeros(N,N);
for j=1:N,
    for i=1:N,
	sum=0;
	for y=-L2:L2,
	    for x=-L2:L2,
		ii=i+x;
		if (ii<1) ii=ii+N;
		elseif (ii>N) ii=ii-N;
		end;

		jj=j+y;
		if (jj<1) jj=jj+N;
		elseif (jj>N) jj=jj-N;
		end;
		
		sum=sum+smth(ii,jj);
	    end;
	end;
	smth2(i,j)=sum/(L*L);
    end;
end;
    
%... scale to normal height ...
mH=max(max(smth2));
mL=min(min(smth2));
for i=1:N,
    for j=1:N,
	smth2(i,j)=(smth2(i,j)-mL)/(mH-mL)*N;
    end;
end;

%... calculate color mapping limits ...
MTOP=max(max(smth2));
MBTM=0;
CTOP=256;

%... process sea level and mountain slope rates ...
smth2=round(smth2*HBR-N*SLR);
for i=1:N,
    for j=1:N,
	if (smth2(i,j)<0) smth2(i,j)=0; end;
    end;
end;

disp('plotting 3D graph...');

%... create custom color mapping ...
cmap=zeros(N,3);
for k=MBTM+1:CTOP,
	cmap(k,:)=[k/CTOP*(0.35+0.35*k/CTOP)+0.30 k/CTOP*(0.70-0.35*k/CTOP)+0.30 0]; end;
cmap(1,:)=[0 0 0.4];	% solid sea color (deep blue)
cmap(CTOP,:)=[1 1 1];	% solid mountain peak color (pure white)
colormap(cmap);

%... plot 3-D terrain surface image ...
figure(1);
colormap(cmap);
surf(smth2);
axis([0 N 0 N 0 N]);
shading interp;
view([-N -N 3*N]);

%... plot 2-D terrain elevation image ...
figure(2);
colormap(cmap);
axis([0 N 0 N 0 N]);
surf(smth2);
shading interp;
view([0 0 N]);
colorbar;

disp('all done!');
disp(' ');