Wiener filter

The Wiener filter is the optimal linear filter for denoising and deblurring.  The filter requires prior knowledge about the autospectrum of the signal, the point spread function of the blur, and the autospectrum of the noise.  In this example white noise is assumed.  The result of the Gaussian blur example is loaded and the Wiener filter is used to estimate the input signal.

See also

• http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/VELDHUIZEN/node15.html

#!/usr/bin/ruby
require 'hornetseye'
include Hornetseye
class Sequence_
  def radius
    w, h = shape[0], shape[1]
    x = Sequence.int( w ).indgen! -w / 2
    y = Sequence.int( h ).indgen! -h / 2
    MultiArray.tensor( 2 ) { |i,j| Math.sqrt( x[i] ** 2 + y[j] ** 2 ) }
  end
  def spectrum( alpha )
    w, h = shape[0], shape[1]
    ( 255.0 / ( 1.0 + ( radius / alpha ) ** 2 ) ).shift -w / 2, -h / 2
  end
  def blur( sigma )
    w, h = shape[0], shape[1]
    retval = MultiArray.ubyte( w, h ).fill!
    retval[ w / 2, h / 2 ] = 1
    retval.gauss_blur( sigma ).shift -w / 2, -h / 2
  end
  def wiener( alpha, noise, sigma = nil )
    if typecode < RGB_
      result = self.class.float.new
      result.r, result.g, result.b = [ r, g, b ].collect do |c|
        c.wiener alpha, noise, sigma
      end
      result
    else
      si = spectrum( alpha ) ** 2
      if sigma
        bf = blur( sigma ).fft
        ( bf.conj * si / ( bf.abs ** 2 * si + noise ** 2 ) * fft ).ifft.real
      else
        ( si / ( si + noise ** 2 ) * fft ).ifft.real
      end
    end
  end
end
if ARGV.size != 2
  raise 'Syntax: ./wiener.rb <input image> <output image>'
end
img = MultiArray.load_ubytergb ARGV[0]
result = img.wiener( img.shape[0] * 0.1, 2.0, 3.0 ).clip
result.save_ubytergb ARGV[1]