% ---------------------------------------------------------------------- %
% Mandelbrot Set w/Activity Detection                                    %
% Jeremy Cowles                                                          %
% EPS 109, Fall 2009                                                     %
% Based on Lab 4                                                         %
% ---------------------------------------------------------------------- %
% Automatically detects interesting points in the Mandelbrot Set by      %
% keeping an activity matrix (ZH). This matrix is a collection of        %
% samples indicating how many points under each sector have diverged.    %
% A magic range (2, 20) is deemed interesting. On every frame a new      %
% center point is calculated based on the highest activity sector that   %
% is closest to the currently selected center point. After this          %
% the view port is shifted toward the current center point using a       %
% smoothing factor to avoid jerky movement.                              %
%                                                                        %
% The activity range (2, 20) was chose because the most interesting      % 
% areas in the Mandelbrot Set are those that lie on borders between      % 
% complete divergence and varying levels of convergence. Since the       % 
% resolution of the activity matrix is 5x5, the maximum value of any     % 
% sector in the matrix is 25, so this is an attempt to find a middle     %
% value in the range. The detection might be improved by adding a        %
% preference for values closest to 12 (not implemented).                 %
% ---------------------------------------------------------------------- %

clf
clear

show_act = 1;
show_center = 0;
speed = 1/30;

%x=[-1.980, -0.930] y=[-0.488, 0.563]
%x=[-4.0, 3.0] y=[-4.0, 3.0]

%xmin = -1.980;
%xmax = -0.930;
%ymin = -0.488;
%ymax = 0.563;

xmin = -4;
xmax = 3;
ymin = -4;
ymax = 3;

%xmin = 0.291;
%xmax=0.448;
%ymin=0.284;
%ymax=0.441;

make_movie = 1;

center = [xmax - xmin, ymax - ymin];

%number of points = resolution in x and y direction
nPoints = 300;
frames  = 200;

%upper limit for |z|, consider series has diverged if exceeded
zmax = 1e2;

% activity matrix sampling ratio
hash_size = 5;

ZZ = zeros(nPoints, nPoints);
center_is_boring = false;
it_max = 100;
x = zeros(1, nPoints);
y = zeros(1, nPoints);

if make_movie
    frame_loc = get(figure(1),'Position');
    frame_loc(1) = 0;
    frame_loc(2) = 0;
end

for ii=1:frames
    ZH = zeros(round(nPoints/hash_size)+1, round(nPoints/hash_size)+1);
    fprintf('ii=%d x=[%.3f, %.3f] y=[%.3f, %.3f]\n', ii, xmin,xmax,ymin,ymax)
    for i=1:nPoints
        % the 'j' index walks along the 'y' direction
        for j=1:nPoints 

            %%%% map the index 'i' onto the interval [xmin,xmax]
            %%%% i=1       must yield 'xmin'
            %%%% i=nPoints must yield 'xmax'
            %x(i) = xmin + (((xmax - xmin) / (nPoints)) * (i));
            x(i) = xmin + ((i-1)/(nPoints-1))*(xmax-xmin);
            %%%% map the index 'j' onto the interval [ymin,ymax]
            %%%% j=1       must yield 'ymin'
            %%%% j=nPoints must yield 'ymax'
            y(j) = ymin + ((j-1)/(nPoints-1))*(ymax-ymin);

            c = complex(x(i),y(j)); % this converts (x,y) to a complex number 'c'
            z = complex(0.0,0.0);   % iteration starts with complex number = 0

            for n = 1:it_max
                %%%% add the iteration formula for the Mandelbrot set
                z = z^2 + c;
                ZZ(j,i)=n;
                
                % get the coords in the activity matrix
                hj = round(j/hash_size)+1;
                hi = round(i/hash_size)+1;

                % Everytime the sample point blows up, record it in the
                % activity matrix. Later look for ranges of activity that
                % might indicate an interesting area of the fractal.
                if (abs(z) > zmax)
                    ZH(hj, hi) = ZH(hj, hi) + 1;
                    break
                end
            end
        end
    end
    
    % is the center of view still interesting?
    if (center_is_boring)
        % dcenter is the dist from center of the new interesting point
        dcenter = [99999,99999];
        
        % loop over the activity matrix, ignoring the outter row and column
        % because their data is wrong due to edge effects
        for hi = 2:size(ZH,1)-1
            for hj = 2:size(ZH,2)-1
                
                % the magic interesting activity range is (2,20)
                if (ZH(hj, hi) > 2 && ZH(hj, hi) < 20)

                    i = hi*hash_size; 
                    j = hj*hash_size; 
                    
                    % found an interesting point, but try to stay as close
                    % to the old focal point as possible, using dcenter
                    ndc = [center(1) - x(i), center(1) - y(j)];
                    
                    if (norm(ndc) < norm(dcenter))
                        % found a new interesting point, update
                        newcenter = [(x(i)+x(i-1))/2, (y(j)+y(j-1))/2];
                        dcenter = ndc;
                        
                        % set in case we want to break out early
                        center_is_boring = false;
                    end
                end
            end
        end
        
        center = newcenter;
    end
    
    axis tight
    axis off
    set(gca,'YDir','normal');
    
    % fractal image plot
    if show_act
        subplot(1,2,1);
    end
    imagesc(x,y,ZZ);
    
    if show_center
        % current focal point
        hold on
        plot(center(1), center(2), 'g+');
        hold off
    end
    if show_act
        % activity matrix plot
        subplot(1,2,2);
        imagesc(1:size(ZH,1),1:size(ZH,2),ZH);
    end
    
    if make_movie
        f = getframe(figure(1), frame_loc); 
        if (ii == 1)
            [im,map] = rgb2ind(f.cdata,256);
            im(1,1,1,frames) = 0;
        end
    
        size(f.cdata);
        size(map);
        im(:,:,1,ii) = rgb2ind(f.cdata, map,'nodither');
    end
    
    % revise the center point on every frame
    % some sort of tollerence could be used instead, 
    % but this seems to work well    
    center_is_boring = true;
    
    % the first frames are unstable, so don't move until they have passed
    if (ii > 4) 
        w = (xmax - xmin)/2;
        h = (ymax - ymin)/2;
        
        % smooth movement toward the center
        smooth = 3;
        xmin = (smooth*xmin + (center(1) - w)) / (smooth+1);
        xmax = (smooth*xmax + (center(1) + w)) / (smooth+1);
        ymin = (smooth*ymin + (center(2) - h)) / (smooth+1);
        ymax = (smooth*ymax + (center(2) + h)) / (smooth+1);
    end
    
    it_max = it_max + 2;

    % determine zoom factor
    dx = (xmax - xmin) * speed;
    dy = (ymax - ymin) * speed;
    
    % Zoom in by zoom factor
    xmin = xmin + dx;
    xmax = xmax - dx;
    ymin = ymin + dy;
    ymax = ymax - dy;
    
    % yield CPU
    pause(0.01)
end

if make_movie
    imwrite(im,map,'DancingPeaks.gif','DelayTime',0,'LoopCount',inf)
end


   