mandelbrot

Result

Code

  # feel free to improve on the code!

  require 'narray'
  require 'png'
  
  class NArray
    # returns all the indices to access the values of the NArray.  if start == 1,
    # then the first dimension (row) values are not returned, if start == 2,
    # then the first two dimensions are skipped.
    #
    # if a block is provided, the indices are yielded one at a time
    # [obviously, this could be made into a fast while loop instead of
    # recursive ... someone have at it]
    def indices(start=0, ar_of_indices=[], final=[], level=shape.size-1, &block)
      if level >= 0
        (1...(shape[level])).each do |s|
          new_indices = ar_of_indices.dup
          new_indices.unshift(s)
          if (new_indices.size == (shape.size - start))
            block.call(new_indices)
            final << new_indices 
          end
          indices(start, new_indices, final, level-1, &block)
        end
      end
      final
    end
  end
  
  class SimpleColormapper
    # basic data from pylab:
    #
    #     for i in arange(256):
    #         print cm.jet(i),
    JET_COLORMAP_256_CHAR_BASE64 = "AAB/AACEAACIAACNAACRAACWAACaAACfAACjAACoAACsAACxAAC2AAC6AAC/\nAADDAADIAADMAADRAADVAADaAADeAADjAADoAADsAADxAAD1AAD6AAD+AAD/\nAAD/AAD/AAD/AAT/AAj/AAz/ABD/ABT/ABj/ABz/ACD/ACT/ACj/ACz/ADD/\nADT/ADj/ADz/AED/AET/AEj/AEz/AFD/AFT/AFj/AFz/AGD/AGT/AGj/AGz/\nAHD/AHT/AHj/AHz/AID/AIT/AIj/AIz/AJD/AJT/AJj/AJz/AKD/AKT/AKj/\nAKz/ALD/ALT/ALj/ALz/AMD/AMT/AMj/AMz/AND/ANT/ANj/ANz+AOD6AOT3\nAuj0BezxCPDtDPTqD/jnEvzkFf/hGP/dHP/aH//XIv/UJf/QKf/NLP/KL//H\nMv/DNv/AOf+9PP+6P/+3Qv+zRv+wSf+tTP+qT/+mU/+jVv+gWf+dXP+aX/+W\nY/+TZv+Qaf+NbP+JcP+Gc/+Ddv+Aef99fP95gP92g/9zhv9wif9sjf9pkP9m\nk/9jlv9fmv9cnf9ZoP9Wo/9Tpv9Pqv9Mrf9JsP9Gs/9Ct/8/uv88vf85wP82\nw/8yx/8vyv8szf8p0P8l1P8i1/8f2v8c3f8Y4P8V5P8S5/8P6v8M7f8I8fwF\n9PgC9/QA+vAA/u0A/+kA/+UA/+IA/94A/9oA/9cA/9MA/88A/8sA/8gA/8QA\n/8AA/70A/7kA/7UA/7EA/64A/6oA/6YA/6MA/58A/5sA/5gA/5QA/5AA/4wA\n/4kA/4UA/4EA/34A/3oA/3YA/3MA/28A/2sA/2cA/2QA/2AA/1wA/1kA/1UA\n/1EA/00A/0oA/0YA/0IA/z8A/zsA/zcA/zQA/zAA/ywA/ygA/yUA/yEA/x0A\n/xoA/xYA/hIA+g8A9QsA8QcA7AMA6AAA4wAA3gAA2gAA1QAA0QAAzAAAyAAA\nwwAAvwAAugAAtgAAsQAArAAAqAAAowAAnwAAmgAAlgAAkQAAjQAAiAAAhAAA\nfwAA\n"
    # simple colormap lending RGB triplet, each from 0 to 255
    JET_COLORMAP = JET_COLORMAP_256_CHAR_BASE64.unpack("m*")[0].unpack("C*").each_slice(3).to_a
  
    attr_accessor :maxval
  
    def initialize(maxval=255)
      @maxval = maxval.to_f
    end
  
    def jet(val)
      index = (val/@maxval*255).floor.to_i
      JET_COLORMAP[index]
    end
  end
  
  def mandelbrot(w, h, zoom=3, maxits=20)
    z = ( NArray.scomplex(w,1).indgen!/w-0.65 )*zoom + 
        ( NArray.scomplex(1,h).indgen!/h-0.5 )*(zoom*1.im)
    c = z.dup
    a = NArray.sint(h,w)
    idx = NArray.int(h,w).indgen!
  
    for i in 1..maxits
      z = z**2+c
      idx_t,idx_f = (z.abs>2).where2
      a[idx[idx_t]] = i
      break if idx_f.size==0
      idx = idx[idx_f]
      z = z[idx_f]
      c = c[idx_f]
    end
    a
  end
  
  height = 400
  width = 400
  
  narr = mandelbrot(width, height, 3.2, 20)
  
  cm = SimpleColormapper.new(narr.max)
  
  canvas = PNG::Canvas.new(width, height)
  
  narr.indices do |n,m|
    canvas[m,n] = PNG::Color.new(*cm.jet(narr[n,m]))
  end
  
  #(0...narr.shape.last).each do |m|
  #  (0...narr.shape.first).each do |n|
  #    canvas[m,n] = PNG::Color.new(*cm.jet(narr[n,m]))
  #  end
  #end
  
  PNG.new(canvas).save "mandelbrot.png"