Parallel Mandelbrot fractal implementation in Octave.

mandelbrot.m

% Returns a matrix that represents an image of the Mandelbrot with the
% given parameters. More info: http://en.wikipedia.org/wiki/Mandelbrot_set
%
% Depends on: parallel
%
% Usage:
%
%   % Processes the fractal in parallel using 4 workers
%   M = mandelbrot (-0.75+0.1i, 200, 512, 600, 600, 4); % Seahorse valley
%
% To display the fractal in a image window:
%
%   colormap (bone (512));
%   image (M);
%
% To save the resulting image to a .jpg file:
%
%   imwrite (M, bone (512), 'mandelbrot.jpg', 'Quality', 100);

function M = mandelbrot (point, zoom, iters, width, height, nworkers=1)

  scale = 1 / zoom;

  % Keeps the count of mandelbrot iterations for each pixel
  M = z = zeros (height, width);

  % Center at (0,0) with radius 1.5 at zoom=1

  floatWidth  = 3 * scale;
  floatHeight = floatWidth * height / width;

  x = (real (point) - floatWidth/2) + (0:width-1) / width * floatWidth;
  y = (imag (point) - floatHeight/2) + (0:height-1) / height * floatHeight;

  % Flips the y coordinate to match the y axis orientation
  [X, Y] = meshgrid (x, fliplr (y));

  % Computes `c` for each pixel in the output image
  c = X + i*Y;

  tic;
  f = @ (c, M, z) mandeliters (iters, c, M, z);
  M = parcellfun (nworkers, f, num2cell (c, 2), num2cell (M, 2), num2cell (z, 2));
  toc;
endfunction

% Computes the set inclusion intensity for every pixel in the image
function M = mandeliters (iters, c, M, z)
  for zn = 1:iters
    mask = abs (z) < 2;

    M(mask) = M(mask) + 1;
    z(mask) = z(mask).^2 + c(mask);
  endfor

  M(mask) = iters;
endfunction