Image Processing with gluas

Figure 1.1. Coordinate system.

Figure 1.2. Vector image

draw circle
     center        0.5, 0.5
     radius        0.4
     fill-color    yellow
     stroke-color  black
     stroke-width  0.05
draw circle
     center        0.35, 0.4
     radius        0.05
     fill-color    black
draw circle
     center        0.65, 0.4
     radius        0.05
     fill-color    black
draw line
     start         0.3, 0.6
     end           0.7, 0.6
     stroke-color  black
     stroke-width  0.1
         

A vector image of a face, and the instructions used to create the image.

Figure 5.2. Over the edge

The 3×3 blue area in the raster on the left, indicates which pixels contribute to the single red pixel in the destination image.

Figure 5.4. blur

-- a box blur filter implemented in lua

edge_duplicate = 1;

function blur_slow_get_rgba(sample_x, sample_y, radius)
   local x,y
   local rsum,gsum,bsum,asum=0,0,0,0
   local count=0

   for x=sample_x-radius, sample_x+radius do
     for y=sample_y-radius, sample_y+radius do
     
       local r,g,b,a= get_rgba(x,y)
       
       rsum, gsum, bsum, asum = rsum+r, gsum+g, bsum+b, asum+a
       
       count=count+1
     end
   end
   return rsum/count, gsum/count, bsum/count, asum/count
end

function blur_slow(radius)
    local x,y


    for y = 0, height-1 do
        for x = 0, width-1 do
           set_rgba (x,y, blur_slow_get_rgba (x,y,radius))
        end
        progress (y/height)
    end
    flush ()
end


-- seperation into two passes, greatly increases the speed of this
-- filter

function blur_h_get_rgba(sample_x, sample_y, radius)
   local x,y
   local rsum,gsum,bsum,asum=0,0,0,0
   local count=0
   
   y=sample_y
   for x=sample_x-radius, sample_x+radius do
       local r,g,b,a= get_rgba(x,y)
       
       rsum, gsum, bsum, asum = rsum+r, gsum+g, bsum+b, asum+a
       
       count=count+1
   end
   return rsum/count, gsum/count, bsum/count, asum/count
end

function blur_v(radius)
    for y = 0, height-1 do
        for x = 0, width-1 do
           set_rgba (x,y, blur_v_get_rgba(x,y,radius))
        end
        progress(y/height)
    end
    flush()
end

function blur_v_get_rgba(sample_x, sample_y, radius)
   local x,y
   local rsum,gsum,bsum,asum=0,0,0,0
   local count=0
   
   x=sample_x
   for y=sample_y-radius, sample_y+radius do
       local r,g,b,a= get_rgba(x,y)
       
       rsum, gsum, bsum, asum = rsum+r, gsum+g, bsum+b, asum+a
       
       count=count+1
   end
   return rsum/count, gsum/count, bsum/count, asum/count
end

function blur_h(radius)
    for y = 0, height-1 do
        for x = 0, width-1 do
           set_rgba (x,y, blur_h_get_rgba(x,y,radius))
        end
        progress(y/height)
    end
    flush()
end

function blur_fast(radius)
  blur_v(radius)
  blur_h(radius)
end

blur_fast(4)

Figure 5.5. kuwahara

-- Kuwahara filter
--
-- 
-- Performs the Kuwahara Filter. This filter is an edge-preserving filter.
--
--
-- ( a  a  ab   b  b)
-- ( a  a  ab   b  b)
-- (ac ac abcd bd bd)
-- ( c  c  cd   d  d)
-- ( c  c  cd   d  d)
--
-- In each of the four regions (a, b, c, d), the mean brightness and the variance are calculated. The
-- output value of the center pixel (abcd) in the window is the mean value of that region that has the
-- smallest variance.
--
-- description copied from http://www.incx.nec.co.jp/imap-vision/library/wouter/kuwahara.html
--
-- implemented by Øyvind Kolås <oeyvindk@hig.no> 2004


-- the sampling window is:  width=size*2+1  height=size*2+1

size = 4

edge_duplicate = 1;

-- local function to get the mean value, and
-- variance from the rectangular area specified
function mean_and_variance (x0,y0,x1,y1)
 local variance
 local mean

 local min         = 1.0
 local max         = 0.0
 local accumulated = 0
 local count       = 0
   
 local x, y

 for y=y0,y1 do
   for x=x0,x1 do
     local v = get_value(x,y)

     accumulated   = accumulated + v
     count         = count + 1
      
     if v<min then min = v end
     if v>max then max = v end

    end
  end

  variance = max-min
  mean     = accumulated /count
   
  return mean, variance
end

-- local function to get the mean value, and
-- variance from the rectangular area specified
function rgb_mean_and_variance (x0,y0,x1,y1)
 local variance
 local mean
 local r_mean
 local g_mean
 local b_mean 
 
 local min           = 1.0
 local max           = 0.0
 local accumulated_r = 0
 local accumulated_g = 0
 local accumulated_b = 0
 local count         = 0
   
 local x, y

 for y=y0,y1 do
   for x=x0,x1 do
     local v     = get_value(x,y)
     local r,g,b = get_rgb (x,y)

     accumulated_r   = accumulated_r + r
     accumulated_g   = accumulated_g + g
     accumulated_b   = accumulated_b + b
     count         = count + 1
      
     if v<min then min = v end
     if v>max then max = v end

    end
  end

  variance = max-min
  mean_r   = accumulated_r /count
  mean_g   = accumulated_g /count
  mean_b   = accumulated_b /count
    
  return mean_r, mean_g, mean_b, variance
end

-- return the kuwahara computed value
function kuwahara(x, y, size)
  local best_mean     = 1.0
  local best_variance = 1.0

  local mean, variance
	
  mean, variance = mean_and_variance (x-size, y-size, x, y)	

  if variance < best_variance then
    best_mean     = mean
    best_variance = variance
  end
	
  mean, variance = mean_and_variance (x, y-size, x+size,y)

  if variance < best_variance then
    best_mean     = mean
    best_variance = variance
  end	
	
  mean, variance = mean_and_variance (x, y, x+size, y+size)

  if variance < best_variance then
    best_mean     = mean
    best_variance = variance
  end	
	
  mean, variance = mean_and_variance (x-size, y, x,y+size)

  if variance < best_variance then
    best_mean     = mean
    best_variance = variance
  end	

  return best_mean
end


-- return the kuwahara computed value
function rgb_kuwahara(x, y, size)
  local best_r, best_g, best_b
  local best_variance = 1.0

  local r,g,b, variance
	
  r,g,b, variance = rgb_mean_and_variance (x-size, y-size, x, y)	

  if variance < best_variance then
    best_r, best_g, best_b = r, g, b
    
    best_variance = variance
  end
	
  r,g,b, variance = rgb_mean_and_variance (x, y-size, x+size,y)

  if variance < best_variance then
    best_r, best_g, best_b = r, g, b
    best_variance = variance
  end	
	
  r,g,b, variance = rgb_mean_and_variance (x, y, x+size, y+size)

  if variance < best_variance then
    best_r, best_g, best_b = r, g, b
    best_variance = variance
  end	
	
  r,g,b, variance = rgb_mean_and_variance (x-size, y, x,y+size)

  if variance < best_variance then
    best_r, best_g, best_b = r, g, b
    best_variance = variance
  end	

  return best_r, best_g, best_b
end

function kuwahara(radius)
  for y=0, height-1 do
    for x=0, width-1 do
      r,g,b = rgb_kuwahara (x,y, radius)
      set_rgb (x,y, r,g,b)
    end
    progress (y/height)
  end
end

kuwahara(2)

Figure 5.6. median

-- a median rank filter

edge_duplicate = 1;

function median_value(sample_x, sample_y, radius)
   local x,y
   local sample = {}
   local samples = 0

   for x=sample_x-radius, sample_x+radius do
     for y=sample_y-radius, sample_y+radius do
     
       local value = get_value(x,y)

       sample[samples] = value
       samples = samples + 1
     end
   end

   table.sort (sample)
   local mid = math.floor(samples/2)
   if math.mod(samples,2) == 1 then
       median = sample[mid+1]
   else
       median = (sample[mid]+sample[mid+1])/2
   end

   return median
end

function median(radius)
    local x,y

    for y = 0, height-1 do
        for x = 0, width-1 do
           set_value (x,y, median_value (x,y,radius))
        end
        progress (y/height)
    end
    flush ()
end

median(3)

Figure 5.7. jitter

-- a jittering filter implemented in lua

edge_duplicate=1 -- set edge wrap behavior to use nearest valid sample

function jitter(amount)
  local x,y
  var2=var1
  for y = 0, height-1 do
    for x = 0, width-1 do
      local nx,ny
      local r,g,b,a
      nx=x+(math.random()-0.5)*amount
      ny=y+(math.random()-0.5)*amount
      r,g,b,a=get_rgba(nx,ny)
      set_rgba(x,y, r,g,b,a)
    end
    progress(y/height)
  end
end

jitter(10)

Figure 6.1. translation

function translate(dx,dy)
  for y=0, height-1 do
    for x=0, width-1 do
      -- calculate the source coordinates (u,v)
      u = x-dx
      v = y-dy
       
      r,g,b=get_rgba(u,v)
      set_rgb(x,y,r,g,b)
    end
    progress (y/height)
  end
  flush()
end

translate(width*0.25,height*0.25)

Figure 6.2. rotation

function rotate(angle)
  for y=0, height-1 do
    for x=0, width-1 do
      -- calculate the source coordinates (u,v)
      u = x * math.cos(-angle) + y * math.sin(-angle)
      v = y * math.cos(-angle) - x * math.sin(-angle)
       
      r,g,b = get_rgba(u,v)
      set_rgb(x,y,r,g,b)
    end
    progress (y/height)
  end
  flush()
end

rotate(math.rad(15))

Figure 6.3. scaling

function scale(ratio)
  for y=0, height-1 do
    for x=0, width-1 do
      -- calculate the source coordinates (u,v)
      u = x * (1.0/ratio)
      v = y * (1.0/ratio)
       
      r,g,b = get_rgba(u,v)
      set_rgb(x,y,r,g,b)
    end
    progress (y/height)
  end
  flush()
end

scale(3)

Figure 6.4. bilinear

function get_rgb_bilinear (x,y)
  local x0,y0 -- integer coordinates
  local dx,dy -- offset from coordinates
  local r0,g0,b0, r1,g1,b1, r2,g2,b2, r3,g3,b3
  local r,g,b

  x0=math.floor(x)
  y0=math.floor(y)
  dx=x-x0
  dy=y-y0
  
  r0,g0,b0 = get_rgb (x0  ,y0  )
  r1,g1,b1 = get_rgb (x0+1,y0  )
  r2,g2,b2 = get_rgb (x0+1,y0+1)
  r3,g3,b3 = get_rgb (x0  ,y0+1)

  r = lerp (lerp (r0, r1, dx), lerp (r3, r2, dx), dy)
  g = lerp (lerp (g0, g1, dx), lerp (g3, g2, dx), dy)
  b = lerp (lerp (b0, b1, dx), lerp (b3, b2, dx), dy)
 
  return r,g,b
end

function scale(ratio)
  for y=0, height-1 do
    for x=0, width-1 do
      -- calculate the source coordinates (u,v)
      u = x * (1.0/ratio)
      v = y * (1.0/ratio)
      r,g,b=get_rgb_bilinear(u,v)
      set_rgb(x,y,r,g,b)
    end
    progress (y/height)
  end
  flush()
end

function lerp(v1,v2,ratio)
   return v1*(1-ratio)+v2*ratio;
end
-- LERP
-- /lerp/, vi.,n.
--
-- Quasi-acronym for Linear Interpolation, used as a verb or noun for
-- the operation. "Bresenham's algorithm lerps incrementally between the
-- two endpoints of the line." (From Jargon File (4.4.4, 14 Aug 2003)
scale(3)

Figure 6.5. bicubic

function get_rgb_cubic_row(x,y,offset)
  local r0,g0,b0, r1,g1,b1, r2,g2,b2, r3,g3,b3
  r0,g0,b0 = get_rgb(x,y)
  r1,g1,b1 = get_rgb(x+1,y)
  r2,g2,b2 = get_rgb(x+2,y)
  r3,g3,b3 = get_rgb(x+3,y)
  return cubic(offset,r0,r1,r2,r3), cubic(offset,g0,g1,g2,g3), cubic(offset,b0,b1,b2,b3)
end

function get_rgb_bicubic (x,y)
  local xi,yi -- integer coordinates
  local dx,dy -- offset from coordinates
  local r,g,b

  xi=math.floor(x)
  yi=math.floor(y)
  dx=x-xi
  dy=y-yi

  r0,g0,b0 = get_rgb_cubic_row(xi-1,yi-1,dx)
  r1,g1,b1 = get_rgb_cubic_row(xi-1,yi,  dx)
  r2,g2,b2 = get_rgb_cubic_row(xi-1,yi+1,dx)
  r3,g3,b3 = get_rgb_cubic_row(xi-1,yi+2,dx)

  return cubic(dy,r0,r1,r2,r3),
         cubic(dy,g0,g1,g2,g3),
         cubic(dy,b0,b1,b2,b3)
end

function scale(ratio)
  for y=0, height-1 do
    for x=0, width-1 do
      -- calculate the source coordinates (u,v)
      u = x * (1.0/ratio)
      v = y * (1.0/ratio)
      r,g,b=get_rgb_bicubic(u,v)
      set_rgb(x,y,r,g,b)
    end
    progress (y/height)
  end
  flush()
end

function cubic(offset,v0,v1,v2,v3)
  -- offset is the offset of the sampled value between v1 and v2
   return   (((( -7 * v0 + 21 * v1 - 21 * v2 + 7 * v3 ) * offset +
               ( 15 * v0 - 36 * v1 + 27 * v2 - 6 * v3 ) ) * offset +
               ( -9 * v0 + 9 * v2 ) ) * offset + (v0 + 16 * v1 + v2) ) / 18.0;
end

scale(3)

Figure 7.2. grayworld assumption

function get_avg_a_b ()
    sum_a=0
    sum_b=0

    -- first find average color in CIE Lab space
    for y=0, height-1 do
      for x=0, width-1 do
        l,a,b = get_lab(x,y)
        sum_a, sum_b = sum_a+a, sum_b+b
      end
      progress(y/height)
    end

    avg_a=sum_a/(width*height)
    avg_b=sum_b/(width*height)

    return avg_a,avg_b
end
 
function shift_a_b(a_shift, b_shift)
    for y=0, height do
      for x=0, width do
        l,a,b = get_lab(x,y)
        
        -- scale the chroma distance shifted according to amount of
        -- luminance. The 1.1 overshoot is because we cannot be sure
        -- to have gotten the data in the first place.
        a_delta = a_shift * (l/100) * 1.1
        b_delta = b_shift * (l/100) * 1.1
        a,b = a+a_delta, b+b_delta

        set_lab(x,y,l,a,b)
      end
      progress(y/height)
    end
    flush()
end

function grayworld_assumption()
  avg_a, avg_b = get_avg_a_b()
  shift_a_b(-avg_a, -avg_b)
end

grayworld_assumption()

Figure 7.3. component stretching

function get_min_max_r ()
    min, max = 1000, -1000          
    for y=0, height-1 do
      for x=0, width-1 do
        value, temp, temp = get_rgb (x,y)
        if value<min then
          min = value
        end
        if value>max then
          max = value
        end
      end  
    end
    return min,max
end

function get_min_max_g ()
    min, max = 1000, -1000          
    for y=0, height-1 do
      for x=0, width-1 do
        temp, value, temp = get_rgb (x,y)
        if value<min then
          min = value
        end
        if value>max then
          max = value
        end
      end  
    end
    return min,max
end

function get_min_max_b ()
    min, max = 1000, -1000          
    for y=0, height-1 do
      for x=0, width-1 do
        temp, temp, value = get_rgb (x,y)
        if value<min then
          min = value
        end
        if value>max then
          max = value
        end
      end  
    end
    return min,max
end

function remap(v, min, max)
    return (v-min) * 1.0/(max-min)
end

function cs_get_rgb(x,y,min_r,max_r,min_g,max_g,min_b,max_b)
    r,g,b = get_rgb(x,y)
    r = remap(r, min_r, max_r)
    g = remap(g, min_g, max_g)
    b = remap(b, min_b, max_b)
    return r,g,b
end

function component_stretch()
    min_r, max_r = get_min_max_r ()
    min_g, max_g = get_min_max_g ()
    min_b, max_b = get_min_max_b ()

    for y=0, height do
      for x=0, width do
        set_rgb(x,y, cs_get_rgb(x,y,min_r,max_r, min_g, max_g, min_b, max_b))
      end
    end
    flush ()
end

component_stretch()

Figure 7.4. ace

edge_duplicate = 1;

function get_min_max (x0,y0,x1,y1)
    min_r, max_r = 1000, -1000          
    min_g, max_g = 1000, -1000          
    min_b, max_b = 1000, -1000          
    for y=y0,y1 do
      for x=x0,x1 do
        r, g, b= get_rgb (x,y)
        if r<min_r then min_r = r end if r>max_r then max_r = r end
        if g<min_g then min_g = g end if g>max_g then max_g = g end
        if b<min_b then min_b = b end if b>max_b then max_b = b end
      end  
    end
    return min_r,max_r,min_g,max_g,min_b,max_b
end

function remap(v, min, max)
    return (v-min) * 1.0/(max-min)
end

function ace_get_rgb(x,y,radius)
    min_r, max_r, min_g, max_g, min_b, max_b =
          get_min_max (x-radius,y-radius,x+radius,y+radius)
    r,g,b = get_rgb(x,y)
    r = remap(r, min_r, max_r)
    g = remap(g, min_g, max_g)
    b = remap(b, min_b, max_b)
    return r,g,b
end

function ace(radius)

    for y=0, height do
      for x=0, width do
        set_rgb(x,y, ace_get_rgb(x,y,radius))
      end
      progress (y/height)
      print (y/height)
    end
    flush ()
end

ace(32)