landscape_matrix.m

利用MATLAB作的分形地形源代码

function [x1,x2,y] = landscape_matrix(n, contraction, flag)

% [X1,X2,Y] = landscape_matrix(N,contraction,FLAG) recusively computes 
% the co-ordinates of the fractal landscape in 3D.
% This file is used in landscape.M.
%
% n gives the number of iterations; factor is the contraction factor;
% flag indicates whether the height difference is chosen randomly - 
% flag = 0 is random, every other flag gives the deterministic picture.

% lowest case of the recursion
if n == 1
	x1 = [  0 0 1/2 ; 1/2 1/2 1 ; 0   0 1/2 ; 1/2 1/2   0];
	x2 = [1/2 0   0 ; 1/2   0 0 ; 1 1/2 1/2 ;   0 1/2 1/2];
	y  = [  0 0   0 ;   0   0 0 ; 0   0   0 ;   0   0   0];
else
    % recursive method call
    [x1,x2,y] = landscape_matrix(n-1, contraction, flag);
    % compute the x-values
    x1 = [x1/2 ; (x1/2 + 1/2) ;   x1/2     ;       (-x1/2 + 1/2)];
    x2 = [x2/2 ;  x2/2      ;    (x2/2 + 1/2)  ;   (-x2/2 + 1/2)];
    if (flag == 0)
        % random case
        % (n+1)x(n+1) matrix with normally distributed random values
        z = (randn(2^n+1))*contraction^n;
    else
        % deterministic case
        z = zeros(2^n + 1);
        z = z + contraction^n;
    end
    yges = [];
    for k = 1:size(y,1) % size(y,1) is the number of rows of y 
		yteil = y(k,:); % y-values of the three nodes
	   	% interpolating values at the middle of the edges
		yteil(4:5) = 0.5 * (yteil(1:2)+yteil(2:3));
 		yteil(6)   = 0.5 * ( yteil(3) + yteil(1) );
		x1teil  = (x1(4*k-3:4*k, :))'; 
        % For k=1 this gives the x1- and x2-values of the left bottom
        % square.
        x2teil  = (x2(4*k-3:4*k, :))';
		% computes the indices of the random matrix
        indices = 2^(n) * [x1teil(:) x2teil(:)] + 1;
        % reads the correct values out of the random matrix
        for i=1:12
			zteil(i) = z(indices(i,1), indices(i,2));
        end
        % yhelp gives the corresponding y-values to the x-values
        % dependent on k: At the middle points of the edges the 
        % corresponding random value is added.
		yhelp = [yteil(4) + zteil(1) , yteil(2)            , ...
                 yteil(5) + zteil(3); ...
		         yteil(6) + zteil(4) , yteil(5) + zteil(5) , ...
                 yteil(3); ...
                 yteil(1)            , yteil(4) + zteil(8) , ...
                 yteil(6) + zteil(9);
		         yteil(5) + zteil(10), yteil(6) + zteil(11), ...
                 yteil(4) + zteil(12)];                                 
       yges  = [yges ; yhelp];
    end
    y = yges;
end