function MandelBrot

% T. Hilmer
% 2011.02.18    version 2
%   updated the resolution and scaling
% 2011.02.15    version 1

% Efficiency goal is to work at exactly visual resolution, i.e. screen pixel

% The axes clim corresponds to the iteration count.  I've found 0:200 is a
% good range.  Changes with complex region. Can be improved by normalizing
% the iteration distribution. Pointless to iterate past the clim

% keeping track of iteration escape count would be useful to see if
% iteration is adding detail...

% TO DO: 
% colormap changes are broken

close all

h.f(1) = figure('ResizeFcn',@figureUpdate,'WindowButtonUpFcn',@figureUpdate);

hz = zoom;
set(hz,'ActionPostCallback',@figureUpdate)
% set(hz,'ActionPreCallback',@zoomCallback)
% The ButtonDownFilter callback is (apparently) only useful for
% enabling/disabling the zoom functionality (conditionally)
% ActionPostCallback doesn't work for getting the zoom box position
% (released position of mouse click), because it reports the post-zoom
% position.

colormap(jet(256)) % above 256 causes the map to wrap at upper limit.  strange

function figureUpdate(obj,~)

ha = get(obj,'Children');

ha = findobj(ha,'Tag','MBrot');    

if length(ha) > 1
    error('UPDATING MULTIPLE AXES!!')
end

updateAxes(ha)

function updateAxes(obj)

if isempty(obj) % no axes -> create
    xlim = [-2 1];
    ylim = [-1 1];
    
    obj = axes;
    % for some bizarre reason matlab is creating multiple axes if I call
    % the settings at axes creation
    set(obj,'Tag','MBrot','XLim',xlim,'YLim',ylim,'CLim',[0 300],...
        'DataAspectRatio',[1 1 1],...
        'Position',[0 0 1 1],... %     following "imshow", border tight is position=[0 0 1 1]
        'Visible','off',...
        'Nextplot','add');
    grid on
end

% initiates new Mandelbrot computation if axes has changed
axesLimits(obj);

function [xlim,ylim,res] = axesLimits(obj)

disp('checking axes')

recompute = false; % default no Mandelbrot computation

% THESE 3 MUST BE DEFINED PRIOR!!
xlim = get(obj,'XLim');
ylim = get(obj,'YLim');
ilim = max(get(obj,'CLim')); % clim+ is the iteration limit

% By default, graphs display in a rectangular axes that has the same aspect
% ratio as the figure window. This makes optimum use of space available for
% plotting.
% Setting DataAspectRatio prevents auto-scaling to fill the figure window 
% One would hope the solution is to change the axes limits such that the automatically scaled aspect ratio is [1 1 1],
% but querying an auto dataspectratio returns semi-constant results that
% don't reflect the actual axis behavior.  Fucking Matlab.
% My solution is to query the Figure Position in pixels and work with that.

pos = get(get(obj,'Parent'),'Position');
pix = pos(3:4); % pixel size

% Pixels, Limits, and resolution in Pixels are related.  Since the GUI
% manipulates the window size in pixels (dragging), and Limits via
% the zoom tool, I've chosen to make resolution the dependent variable.

res = [diff(xlim) diff(ylim)]./pix;% complex resolution scaled to pixel resolution
% simple solution to resizing is to use the coarser resolution:
res = max(res);

D = get(obj,'UserData');

if ~isfield(D,'xlim') || any(xlim~=D.xlim) || any(ylim~=D.ylim)
    display(sprintf('New Complex Limits: %E : %E , %E : %E ',xlim,ylim))
    recompute = true;
end
if ~isfield(D,'ilim') || ilim~=D.ilim
    display(sprintf('New Iteration Limit: %E',ilim))
    recompute = true;
end
if ~isfield(D,'res') || any(res~=D.res) ||  ~isfield(D,'pix') || any(pix~=D.pix)
    display(sprintf('New Resolution: %.0f : %.0f px, %E : %E',pix,res))
    recompute = true;
end

D.xlim = xlim;
D.ylim = ylim;
D.ilim = ilim;
D.res = res;
D.pix = pix;
set(obj,'UserData',D);

if recompute
    x = xlim(1):res:xlim(2);
    y = ylim(1):res:ylim(2);
    
    dim = [length(y) length(x)];
    
    C = complex(repmat(x,dim(1),1),repmat(y(:),1,dim(2)));
    
    mandelBrot(C,ilim)
else
    disp('axes unchanged')
end

function mandelBrot(C,ilim)

verbose = true;

disp('start mandelbrot')

kind = 'mandel';
% kind = 'buddha';

dim = size(C);
N = prod(dim);

% resolution could be passed in as a variable, but this should be exactly
% the same
x = real(C(1,:));
y = imag(C(:,1));

res = [median(diff(x)) median(diff(y))];
% simple solution to resizing is to use the coarser resolution:
res = max(res);

Z = zeros(dim);% complex double
j = int64(1:N); % Set indices (retained)
% Iteration count.  This variable retains full size, since it will be the
% image:
% Escape count:
ne = zeros(ilim,1); % # escaped for each iteration.  Possibly interesting.
I = zeros(dim,'uint8'); % format auto-grows with index as necessary
I(:) = uint8(Inf); % initial state is infinite iterations (never escape). Mostly for colormap

% The path of C -> Buddhabrot. The trajectory image.  Also full size.

% t = zeros(2^pw,1);

if verbose
%         fout = []; fprintf('\n')
end

n = uint64(-1);% integer accuracy is better approaching 2^64.  necessary for logic

while true
    
    if verbose
        if ~rem(n,100)
%             fprintf(repmat('\b',1,length(fout)))
%             fout = sprintf('%.0f',n);
%             fwrite(1,fout)
        end
    end
    
    if false % ~rem(n,100) % n == 255 % 
%         fout = []; fprintf('\n')
        fprintf('displaying %d iterations\n',n)
        imageBrot(x,y,I,res,n)
    end
    
    if n == ilim
        break
    end
    
    % auto-grow numerical format
    if n == uint8(Inf)
        %         fout = []; fprintf('\n')
        disp('upconverting to 16-bit')
        I = uint16(I);
        I(I==uint8(Inf)) = uint16(Inf); % move Inf up
    elseif n == uint16(Inf)
        %         fout = []; fprintf('\n')
        disp('upconverting to 32-bit')
        I = uint32(I);
        I(I==uint16(Inf)) = uint32(Inf);
    elseif n == uint32(Inf)
        %         fout = []; fprintf('\n')
        disp('upconverting to 64-bit')
        I = uint64(I);
        I(I==uint32(Inf)) = uint64(Inf);
    elseif n == uint64(Inf)
        %         fout = []; fprintf('\n')
        warning('reached the 2^64 iteration limit')
        save([datestr(now,30) '_mbrot'])
    end
    
    %     tic
    % core.  where the magic happens
    Z = Z.^2 + C;
    k = abs(Z) > 2; % C is outside Mandelbrot set domain.
    ne(n+1) = sum(k(:));
    
    % Mandelbrot section:
    %----------------------
    if strcmp('mandel',kind)
        I(j(k)) = n; % Occured at this iteration.
    end
    %----------------------
    
    % Shared between methods
    % discard outside set:
    C(k) = [];
    Z(k) = [];
    j(k) = [];
    
    % Buddhabrot section:
    %----------------------
    % this might be wrong.  think i'm supposed to be tracking the escapes,
    % not retains
    if strcmp('buddha',kind)
        % out-of-range Z already discarded
        % Z positions rounded to the nearest grid indice:
        l = 1+ round((imag(Z)-min(y))/res(2));% vertical dimension
        m = 1+ round((real(Z)-min(x))/res(1));% "min(x)" is the left limit
        %     jp = sub2ind(dim,l,m);
        % profiler said the sub2ind call was expensive
        % 2.633 vs 0.64 s (20 vs 18 s total).  savings is real
        jp = l + (m-1)*dim(1);
        try
        I(jp) = I(jp) + 1; % increment path-over
        catch
            4
            error('fix indices')
        end
    end
    %----------------------
    
    %     t(n+1) = toc; % timekeeping
    n = n + 1;
end
% fprintf('\n')

imageBrot(x,y,I,res,n)
disp('finished mandelbrot')

function imageBrot(x,y,I,res,iter)
% imagesc + higher resolution on top

hi = imagesc(x,y,I);

% store image related information
if 0
    I(I==0) = [];
    pdf = sort(I(:),'descend');
    D.pdf = pdf(floor(0.05*numel(I))); % pdf -> viewing scale
    % don't have a working pdf normalization method coded.  It's on the
    % wikipedia page
end
D.res = res; % resolution
D.iter = iter; % iteration count

set(hi,'UserData',D)

% highest resolution on top:

ha = get(hi,'Parent'); % parent axes
hc = get(ha,'Children'); % same level children

r = get(hc,'UserData');
if iscell(r)
    r = [r{:}];
end
r = cat(1,r.res);
[~,j] = sort(r);
set(ha,'Children',hc(j))

drawnow