translate Munsell colours to RGB triplets

munsell.rb

require 'matrix'
require 'interpolator'

class Munsell
  HUES = %w{R YR Y GY G BG B PB P RP}.map(&:to_sym)
  INVHUES = %w{RY YG GB BP PR}.map(&:to_sym)
  
  # from http://www.color.org/srgb.pdf
  SRGBMAT = Matrix[[ 3.2406, -1.5372, -0.4986],
                   [-0.9689,  1.8758,  0.0415],
                   [ 0.0557, -0.2040,  1.0570]]
  
  def initialize hue, value, chroma
    @hue = Hue.normalise hue
    @value = value.to_f
    @chroma = chroma.to_f
  end
  
  def to_srgb outside_range=true
    xyY = interpolate
    xyy_to_srgb(xyY).map {|c| (c * 255).round }
  end
  alias to_rgb to_srgb
  
  def make_interpolate_map
    # eg. map[:H][:y] maps H to y
    map = {:H => {:x => {}, :y => {}, :Y => {}},
           :V => {:x => {}, :y => {}, :Y => {}},
           :C => {:x => {}, :y => {}, :Y => {}}}
    XYYTABLE.each do |hvc, xyY|
      h,v,c = hvc
      _x, _y, y = xyY
      dh = munsell_hue_to_degrees h
      c = c.to_f
      
      [[:H, dh], [:V, v], [:C, c]].each do |coord|
        mcoord, mval = coord
        [[:x, _x], [:y, _y], [:Y, y]].each do |coord|
          ccoord, mval = coord
          if mcoord == :H
            # H to s: V => C => H
            map[mcoord][ccoord][v] ||= {}
            map[mcoord][ccoord][v][c] ||= {}
            map[mcoord][ccoord][v][c][dh] = mval
          elsif mcoord == :V
            # V to s: H => C => V
            map[mcoord][ccoord][dh] ||= {}
            map[mcoord][ccoord][dh][c] ||= {}
            map[mcoord][ccoord][dh][c][v] = mval
          elsif mcoord == :C
            # C to s: H => V => C
            map[mcoord][ccoord][dh] ||= {}
            map[mcoord][ccoord][dh][v] ||= {}
            map[mcoord][ccoord][dh][v][c] = mval
          end
        end
      end
    end
    
    map.each do |mcoord, ccoords|
      map[mcoord].each do |ccoord, m_0s|
        map[mcoord][ccoord].each do |m_0, m_1s|
          map[mcoord][ccoord][m_0].each do |m_1, m_2s|
            map[mcoord][ccoord][m_0][m_1] = Interpolator::Table.new(m_2s) {|t| t.style = Interpolator::Table::CATMULL }
          end
          map[mcoord][ccoord][m_0] = Interpolator::Table.new(m_1s) {|t| t.style = Interpolator::Table::CATMULL }
        end
        map[mcoord][ccoord] = Interpolator::Table.new(m_0s) {|t| t.style = Interpolator::Table::CATMULL }
      end
    end
    
    map
  end
  
  def interpolate
    @@interp_map ||= make_interpolate_map
    dh = munsell_hue_to_degrees @hue
    _x = []
    _y = []
    y = []
    
    [[:H, dh], [:V, @value], [:C, @chroma]].each do |coord|
      mcoord, mval = coord
      [[:x, _x], [:y, _y], [:Y, y]].each do |coord|
        ccoord, cval = coord
        map = @@interp_map[mcoord][ccoord]
        begin
          if mcoord == :H
            # H to s: V => C => H
            cval << map.interpolate(@value, @chroma, dh)
          elsif mcoord == :V
            # V to s: H => C => V
            cval << map.interpolate(dh, @chroma, @value)
          elsif mcoord == :C
            # C to s: H => V => C
            cval << map.interpolate(dh, @value, @chroma)
          end
        rescue RuntimeError
        end
      end
    end
    
    [[:x, _x], [:y, _y], [:Y, y]].map do |coords|
      name, coords = coords
      if coords.empty?
        raise InterpolationError, "no value found for #{name}"
      end
      l = coords.length
      coords.inject(:+) / l.to_f
    end
  end
  
  def xyy_to_srgb xyY
    x, y, z = xyy_to_xyz xyY
    # 'LINear ColourS'
    lincs = (SRGBMAT * Matrix[[x], [y], [z]])
    [lincs[0,0], lincs[1,0], lincs[2,0]].map {|c| lin_to_srgb c }
  end
  def xyy_to_xyz xyY
    _x, _y, y = xyY
    x = (y * _x) / _y
    z = (y * (1 - _x - _y)) / _y
    [x, y, z]
  end
  def lin_to_srgb c
    if c <= 0.0031308
      12.92 * c
    else
      (1.055 * (c ** (1 / 2.4))) - 0.055
    end
  end
  
  def munsell_hue_to_degrees hue
    codes = {R: 342, YR: 18, Y: 54, GY: 90, G: 126, BG: 162, B: 198, PB: 234, P: 270, RP: 306}
    ((hue.step * 3.6) + codes[hue.hue]) % 360
  end
  
  def inspect
    "#<Munsell: #{@hue.inspect} #@value/#@chroma>"
  end
  
  class Hue
    attr_reader :step, :hue
    def initialize step, hue
      raise ArgumentError unless step.is_a? Float
      raise ArgumentError unless HUES.include? hue
      @step = step
      @hue = hue
    end
    
    def hash
      @step.hash ^ @hue.hash
    end
    def == b
      @hue == b.hue and @step == b.step
    end
    alias eql? ==
    
    def self.parse str
      _, step, hue = str.match(/\A(\d+(?:\.\d+)?)([A-Z]+)\Z/).to_a
      self.new step.to_f, hue.to_sym
    end
    def self.normalise input
      case input
      when Hue
        input
      when String
        self.parse input
      when Array
        self.new *input
      end
    end
    
    def inspect
      "#@step#@hue"
    end
  end
  
  class Error < Exception; end
  class InterpolationError < Error; end
  
  module NumericMixin
    HUES.each do |hue|
      define_method hue do
        Hue.new self.to_f, hue
      end
    end
    INVHUES.each do |hue|
      define_method hue do
        Hue.new self.to_f, hue.to_s.reverse.to_sym
      end
    end
  end
  
  XYYTABLE = {}
  File.open('all.dat') do |f|
    while line = f.gets
      next if line.include? 'h'
      entry = line.split(/\s+/)
      XYYTABLE[ [Hue.parse(entry[1]), entry[2].to_f, entry[3].to_i] ] = [entry[4].to_f, entry[5].to_f, (entry[6].to_f / 100.0)]
    end
  end
end