function figh = juliaSetExplorer(z0)
%juliaSetExplorer An interface for exploring the Mandelbrot Julia Set
%
% juliaSetExplorer() launches a simple interface that allows the Julia
% Set for any location in the Mandelbrot set to be viewed. The location
% can be varied by clicking or dragging in a small overview window in the
% top left. Close the window to exit.
%
% figh = juliaSetExplorer() also returns a handle to the interface window
% so that it can be captured or close programmatically.
% Author: Ben Tordoff
% Copyright 2010-2011 The MathWorks, Inc.
% First make sure we are capable of running this
matlabVersionCheck()
% Create the shared data structure
data = struct();
data.WriteVideo = false;
data.DoAnimation = false;
data.LocationList = makeLocationList();
data.NextLocation = 2;
if nargin>0
data.z0 = z0;
else
data.z0 = data.LocationList(1);
end
x = gpuArray.linspace( -2, 2, 800 );
y = gpuArray.linspace( -2, 2, 600 );
[x,y] = meshgrid(x, y);
data.z = complex( x, y );
clear x y;
% Build the interface
gui = createGUI();
% Are we capturing video
if data.WriteVideo
data.VideoWriter = VideoWriter('juliaSetExplorer.avi');
data.VideoWriter.FrameRate = 20;
data.VideoWriter.Quality = 90;
open( data.VideoWriter );
end
% Set the initial view
drawMandelbrot( gui.MandelImage );
onLimitsChanged();
% Animate a set path
if data.DoAnimation
doAnimation();
end
% Return a handle to the figure if requested
if nargout > 0
figh = gui.Window;
end
function gui = createGUI()
gui = struct();
gui.Window = figure( ...
'Name', 'Julia Set explorer v1.1', ...
'Renderer', 'ZBuffer', ...
'NumberTitle', 'off', ...
'HandleVisibility', 'off', ...
'MenuBar', 'none', ...
'ToolBar', 'figure', ...
'CloseRequestFcn', @onFigureClose );
gui.JuliaAxes = axes( ...
'Parent', gui.Window, ...
'Position', [0 0 1 1], ...
'XLim', [-2 2], 'YLim', [-2 2], ...
'XTick', [], 'YTick', [], ...
'DataAspectRatio', [1 1 1] );
% Add listeners so that we can redraw when the axes are moved
axHandle = handle( gui.JuliaAxes );
gui.Listeners = [
handle.listener( axHandle, findprop( axHandle, 'YLim' ), 'PropertyPostSet', @onLimitsChanged )
]; %#ok<NBRAK>
% Also add a small set of axes for showing the Mandelbrot set
gui.MandelAxes = axes( ...
'Parent', gui.Window, ...
'Position', [0 0.8 0.2 0.2], ...
'XLim', [-2 0.5], 'YLim', [-1.25 1.25], ...
'XTick', [], 'YTick', [], ...
'ButtonDownFcn', @onMandelButtonDown );
gui.JuliaImage = imagesc( 'Parent', gui.JuliaAxes, ...
'CData', nan, ...
'XData', [-2 2], 'YData', [-2 2] );
colormap( gui.JuliaAxes, jet2(1000) );
% Add a line so that zooming works. Strange but true.
line( 'Parent', gui.JuliaAxes, ...
'XData', [-2 2], 'YData', [-2 2], ...
'Visible', 'off', ...
'HitTest', 'off' );
gui.MandelImage = imagesc( 'Parent', gui.MandelAxes, 'CData', nan, ...
'XData', [-2 0.5], 'YData', [-1.25 1.25], 'HitTest', 'off' );
gui.MandelCrosshair = [
line( 'Parent', gui.MandelAxes, 'XData', [-2 0.5], 'YData', [0 0], 'Color', 'w', ...
'HitTest', 'off', 'Tag', 'CrossHairH' );
line( 'Parent', gui.MandelAxes, 'YData', [-1.25 1.25], 'XData', [0 0], 'Color', 'w', ...
'HitTest', 'off', 'Tag', 'CrossHairV' );
];
% Remove some things we don't want from the toolbar and add a
% toggle to the toolbar to hide the mandelbrot view
tb = findall( gui.Window, 'Type', 'uitoolbar' );
delete( findall( tb, 'Tag', 'Standard.FileOpen' ) );
delete( findall( tb, 'Tag', 'Standard.NewFigure' ) );
delete( findall( tb, 'Tag', 'Standard.EditPlot' ) );
delete( findall( tb, 'Tag', 'Exploration.Brushing' ) );
delete( findall( tb, 'Tag', 'Exploration.DataCursor' ) );
delete( findall( tb, 'Tag', 'Exploration.Rotate' ) );
delete( findall( tb, 'Tag', 'DataManager.Linking' ) );
delete( findall( tb, 'Tag', 'Plottools.PlottoolsOn' ) );
delete( findall( tb, 'Tag', 'Plottools.PlottoolsOff' ) );
delete( findall( tb, 'Tag', 'Annotation.InsertLegend' ) );
delete( findall( tb, 'Tag', 'Annotation.InsertColorbar' ) );
gui.ShowMandelToggle = uitoggletool( ...
'Parent', tb, ...
'CData', readIcon( 'icon_mandel.png' ), ...
'TooltipString', 'Show/hide the Mandelbrot view', ...
'State', 'on', ...
'Separator', 'on', ...
'ClickedCallback', @onMandelTogglePressed );
gui.AnimToggle = uitoggletool( ...
'Parent', tb, ...
'CData', readIcon( 'icon_play.png' ), ...
'TooltipString', 'Play an animation', ...
'State', 'off', ...
'ClickedCallback', @onAnimTogglePressed );
% Set the resize function (we can't do this on construction as it
% would fire!
set( gui.Window, 'ResizeFcn', @onFigureResize );
end % createGUI
function cdata = readIcon( filename )
[cdata,~,alpha] = imread( filename );
idx = find( ~alpha );
page = size(cdata,1)*size(cdata,2);
cdata = double( cdata ) / 255;
cdata(idx) = nan;
cdata(idx+page) = nan;
cdata(idx+2*page) = nan;
end % readIcon
function onMandelButtonDown( ~, ~ )
pos = get( gui.MandelAxes, 'CurrentPoint' );
% fprintf( 'Click at (%f,%f)\n', pos(1,1), pos(1,2) );
set( gui.Window, ...
'WindowButtonMotionFcn', @onMandelButtonMotion, ...
'WindowButtonUpFcn', @onMandelButtonUp );
updatePosition( complex( pos(1,1), pos(1,2) ) )
end % onMandelButtonDown
function onMandelButtonMotion( ~, ~ )
pos = get( gui.MandelAxes, 'CurrentPoint' );
updatePosition( complex( pos(1,1), pos(1,2) ) )
end % onMandelButtonMotion
function onMandelButtonUp( ~, ~ )
pos = get( gui.MandelAxes, 'CurrentPoint' );
set( gui.Window, ...
'WindowButtonMotionFcn', [], ...
'WindowButtonUpFcn', [] );
updatePosition( complex( pos(1,1), pos(1,2) ) )
end % onMandelButtonUp
function updatePosition( z0 )
if ~ishandle( gui.Window )
return;
end
data.z0 = z0;
drawMandelbrotCrosshair( gui.MandelCrosshair, data.z0 );
drawJulia( gui.JuliaImage, gui.JuliaAxes, data.z, data.z0, data.DoAnimation );
% Capture!
if data.WriteVideo
currFrame = getframe( gui.Window );
writeVideo( data.VideoWriter, currFrame );
else
drawnow();
end
end % updatePosition
function onLimitsChanged( ~, ~ )
% To work out what to draw and at what resolution we need the axis
% limits and pixel counts.
xlim = get(gui.JuliaAxes,'XLim');
ylim = get(gui.JuliaAxes,'YLim');
pixpos = getpixelposition( gui.JuliaAxes );
x = gpuArray.linspace( xlim(1), xlim(2), pixpos(3) );
y = gpuArray.linspace( ylim(1), ylim(2), pixpos(4) );
[x,y] = meshgrid(x, y);
data.z = complex( x, y );
set( gui.JuliaImage, 'XData', xlim, 'YData', ylim );
% Use "update position to force a redraw"
updatePosition( data.z0 );
end % onLimitsChanged
function onFigureClose( ~, ~ )
% Clear up
if data.WriteVideo
close( data.VideoWriter );
end
delete( gui.Window );
end % onFigureClose
function onFigureResize( ~, ~ )
pos = getpixelposition( gui.JuliaAxes );
aspect = pos(4)/pos(3);
% Make sure the Mandelbrot preview is square and not too big
mandelX = min( 0.2, 200/pos(3) );
mandelY = mandelX/aspect;
set( gui.MandelAxes, 'Position', [0 0 mandelX mandelY] );
% Set the Julia-set axes to exactly fill the figure and the
% resolution to match the axes size
pos = get( gui.Window, 'Position' );
xlim = get( gui.JuliaAxes, 'XLim' );
ylim = get( gui.JuliaAxes, 'YLim' );
delta_ylim = ( diff( xlim )*pos(4)/pos(3) - diff( ylim ) ) / 2;
data.WindowPixelSize = pos(3:4);
% Set the YLim to give the correct aspect. This will trigger a
% redraw
set( gui.JuliaAxes, 'YLim', ylim + delta_ylim*[-1 1] );
end % onFigureResize
function onMandelTogglePressed( ~, ~ )
% Toggle the Mandelbrot view on and off
state = get( gui.ShowMandelToggle, 'State' );
set( gui.MandelAxes, 'Visible', state, 'HitTest', state );
set( gui.MandelImage, 'Visible', state );
set( gui.MandelCrosshair, 'Visible', state );
end
function onAnimTogglePressed( ~, ~ )
% Toggle the Mandelbrot view on and off
data.DoAnimation = strcmpi( get( gui.AnimToggle, 'State' ), 'on' );
if data.DoAnimation
doAnimation();
end
end
function doAnimation()
while data.DoAnimation
if ~ishandle( gui.Window )
return;
end
updatePosition( data.LocationList(data.NextLocation) );
data.NextLocation = mod( data.NextLocation, numel( data.LocationList ) ) + 1;
end
% Do a final redraw with animation off
updatePosition( data.LocationList(data.NextLocation) );
end % doAnimation
end
function drawMandelbrot( imh )
escapeRadius2 = 400; % Square of escape radius
maxIters = 500;
x = gpuArray.linspace( -2, 0.5, 500 );
y = gpuArray.linspace( -1.2, 1.2, 500 );
[x0,y0] = meshgrid(x, y);
logCount = arrayfun( @processMandelbrotElement, x0, y0, escapeRadius2, maxIters );
set( imh, ...
'CData', gather(logCount), ...
'XData', [-2 0.5], 'YData', [-1.2 1.2] );
set( ancestor( imh, 'axes' ), 'XLim', [-2 0.5], 'YLim', [-1.2 1.2] );
end % drawMandelbrot
function drawMandelbrotCrosshair( l, z0 )
hline = l(1);
vline = l(2);
set( vline, 'XData', real(z0)*[1 1] );
set( hline, 'YData', imag(z0)*[1 1] );
end % drawMandelbrotCrosshair
function drawJulia( imh, axh, z, z0, animating )
maxIterations = 200;
escapeRadius2 = 100; % Square of escape radius
if animating && numel(z)>1e6
z = z(1:2:end,1:2:end);
end
logCount = arrayfun( @processJuliaSetElement, z, z0, escapeRadius2, maxIterations );
set( imh, 'CData', gather( logCount ) );
set( axh, 'CLim', [1 log(maxIterations)] );
end % drawJulia
function [logCount,t] = processMandelbrotElement( x0, y0, escapeRadius2, maxIterations )
% Evaluate the Mandelbrot function for a single element
t = 1;
z0 = complex( x0, y0 );
z = z0;
count = 0;
while count <= maxIterations && (z*conj(z) <= escapeRadius2)
z = z*z + z0;
count = count + 1;
end
magZ2 = max(real(z).^2 + imag(z).^2,escapeRadius2);
logCount = log( count + 1 - log( log( magZ2 ) / 2 ) / log(2) );
end % processMandelbrotElement
function logCount = processJuliaSetElement( z, z0, escapeRadius2, maxIterations )
% Evaluate the Julia set function for a single element
count = 0;
magZ2 = real(z)^2 + imag(z)^2;
while count <= maxIterations && (magZ2 <= escapeRadius2)
z = z*z + z0;
count = count + 1;
magZ2 = real(z)^2 + imag(z)^2;
end
% Iterate twice more to help smoothing
z = z.*z + z0; z = z.*z + z0; count = count + 2;
% Now adjust the count using the magnitude to get smoothing
magZ2 = max( real(z)^2 + imag(z)^2, escapeRadius2 );
logCount = log( count + 2 - log( log( magZ2 ) / 2 ) / log(2) );
end % processJuliaSetElement
function matlabVersionCheck()
% R2011b is v7.13
majorMinor = sscanf( version, '%d.%d' );
if (majorMinor(1)<7) || (majorMinor(1)==7 && majorMinor(2)<13)
error( 'mandelbrotViewer:MATLABTooOld', 'mandelbrotViewer requires MATLAB R2011b or above.' );
end
end % matlabVersionCheck
function cmap = jet2(m)
% Jet colormap with added fade to black
% A list of break-point colors
colors = [
0.0 0.0 0.5
0.0 0.0 1.0
0.0 0.5 1.0
0.0 1.0 1.0
0.5 1.0 0.5
1.0 1.0 0.0
1.0 0.5 0.0
1.0 0.0 0.0
0.5 0.0 0.0
0.5 0.0 0.0
1.0 0.0 0.0
1.0 0.5 0.0
1.0 1.0 0.0
0.5 1.0 0.5
0.0 1.0 1.0
0.0 0.5 1.0
0.0 0.0 1.0
0.0 0.0 0.5
0.0 0.0 0.0
];
% Now work out the indices into the map
N = size( colors, 1 );
idxIn = 1:N;
idxOut = linspace( 1, N, m );
cmap = [
interp1( idxIn, colors(:,1), idxOut )
interp1( idxIn, colors(:,2), idxOut )
interp1( idxIn, colors(:,3), idxOut )
]';
end % jet2
function locations = makeLocationList()
breakpoints = [
0.29
0.47 + 0.157143i
0.4442 + 0.3286i
0.354911 + 0.585714i
0.2103 + 0.5551i
-0.0827 + 0.8355i
-0.158482 + 1.042857i
-0.5826 + 0.6143i
-0.7533 + 0.1178i
-1.008216 + 0.35i
-1.25 + 0.2i
-1.25
0.29
];
numSteps = 10000;
N = numel(breakpoints);
re = interp1( 1:N, real( breakpoints ), linspace( 1, N, numSteps ) );
im = interp1( 1:N, imag( breakpoints ), linspace( 1, N, numSteps ) );
locations = complex( re, im );
end
