create eigen image

eigen_image.rb

# encoding: utf-8
#

require "narray"
require "png"

def read_pgm(filepath)
  open(filepath) do |io|
    raise "#{filepath}: only P5(Graymap rawbits) is supported" unless io.gets == "P5\n"
    width, height = io.gets.split(/\s+/).map{|i| i.to_i }
    raise "#{filepath}: only 255 depth is supported" unless io.gets == "255\n"
    buf = io.read(width*height).unpack("C*")
    [width, height, buf]
  end
end

# read pgm files
graymap_vectors = ARGV.map do |f|
  read_pgm(f)
end

unless graymap_vectors.all?{|v|
         v[0] == graymap_vectors[0][0] and
         v[1] == graymap_vectors[0][1] }
  raise "graymap image data size mismatch"
end
width = graymap_vectors[0][0]
height = graymap_vectors[0][0]
dimention = width * height
graymap_vectors.map!{|v| v[2]}

# normalize vectors (0..255) => (0.0 .. 1.0)
graymap_vectors.map! do |vec|
  vec.map do |i|
    i/255.0
  end
end

# avarage vector
avarage = Array.new(dimention, 0)
graymap_vectors.each do |vec|
  dimention.times do |i|
    avarage[i] += vec[i]
  end
end
avarage.map! do |val|
  val / graymap_vectors.size
end

# subtract avarage
graymap_vectors.each do |vec|
  dimention.times do |i|
    vec[i] -= avarage[i]
  end
end

# obtain covariance matrix
cmatrix = NMatrix.float(dimention, dimention)
graymap_vectors.each do |vec|
  vertical = NMatrix.float(1, dimention)
  horizontal = NMatrix.float(dimention, 1)
  dimention.times do |i|
    vertical[i] = vec[i]
    horizontal[i] = vec[i]
  end
  cmatrix += vertical * horizontal
end
cmatrix /= graymap_vectors.size

# obtain set of eigen value & eigen vector by Power iteration

def power_iteration(matrix)
  eigen_vec = NVector.float(matrix.shape.first).fill(1.0)
  eigen_value = 0.0
  # iterate until eigen value converge
  100.times do |idx|
    vec = matrix * eigen_vec
    norm = 0.0
    vec.each do |val| norm += val ** 2 end
    eigen_prev = eigen_value
    eigen_value = Math.sqrt(norm)
    eigen_vec = vec / eigen_value
    puts "#{idx+1} iterate...: #{eigen_value}(#{eigen_value-eigen_prev})"
    if (eigen_value-eigen_prev).abs < 0.001
      break
    end
  end
  [eigen_value, eigen_vec]
end

eigens = []
matrix = cmatrix
100.times do |i|
  puts "#{i+1} eigen..."
  eigen_value, eigen_vector = power_iteration(matrix)
  eigens.push([eigen_value, eigen_vector])

  if eigen_value < 1.0
    break
  end

  vertical = NMatrix.float(1, dimention)
  horizontal = NMatrix.float(dimention, 1)
  dimention.times do |i|
    vertical[i] = eigen_vector[i]
    horizontal[i] = eigen_vector[i]
  end
  matrix -= vertical * horizontal * eigen_value
end

eigens.each_with_index do |eigen_pair, index|
  eigen_value, eigen_vector = eigen_pair
  puts "eigen value[#{index}] = #{eigen_value}"
  canvas = PNG::Canvas.new(width, height)
  height.times do |iy|
    width.times do |ix|
      v = eigen_vector[iy*width+ix]
      v = ((v*eigen_value+1.0)*0.5*255).round
      if v < 0
        v = 0
      elsif 255 < v
        v = 255
      end
      canvas[ix, height-1-iy] = PNG::Color.new(v, v, v)
    end
  end
  PNG.new(canvas).save("eigen_image#{index+1}.png")
end

# normalize eigen vectors
eigens.map! do |eigen_val, eigen_vec|
  norm = 0.0
  eigen_vec.each do |v| norm += v ** 2 end
  eigen_vec / norm
end

open("eigen_vectors.dat", "w") do |port|
  Marshal.dump(eigens.map{|eigen_vec| eigen_vec.to_a }, port)
end

# And more...

factor_analyze.rb

# encoding: utf-8
#

require "narray"

def read_pgm(filepath)
  open(filepath) do |io|
    raise "#{filepath}: only P5(Graymap rawbits) is supported" unless io.gets == "P5\n"
    width, height = io.gets.split(/\s+/).map{|i| i.to_i }
    raise "#{filepath}: only 255 depth is supported" unless io.gets == "255\n"
    buf = io.read(width*height).unpack("C*")
    [width, height, buf]
  end
end

def analyze_factor(pgmfile, eigen_vectors)
  width, height, buf = read_pgm(pgmfile)
  # normalize (0..255) => (-1.0..1.0)
  buf.map! do |val|
    (val / 255.0) * 2.0 - 1.0
  end
  vec = NVector.to_narray(buf)
  eigen_vectors.map do |ev|
    vec * ev
  end
end

# load eigen vectors
eigen_vectors = open("eigen_vectors.dat", "rb") do |f|
  Marshal.load(f).map do |v|
    NVector.to_narray(v)
  end
end

# read pgm files
ARGV.each do |f|
  factors = analyze_factor(f, eigen_vectors)
  puts "#{f}: #{factors.inspect}"
end