Image Processing with gluas

We need a coordinate system to describe an image, the coordinate system used to place elements in relation to each other is called user space, since this is the coordinates the user uses to define elements and position them in relation to each other.

Figure 1.1. Coordinate system.

The coordinate system used for all examples in this document has the origin in the upper left, with the x axis extending to the right and y axis extending downwards.

It would have been nice to make a smiling face, instead of the dissatisfied face on the left, by using a bezier curve, or the segment of a circle this could be achieved, this being a text focusing mainly on raster graphics though, that would probably be too complex.

A simple image of a face can be declared as follows:

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.

The preceding description of an image can be seen as a “cooking recipe” for how to draw the image, it contains geometrical primitives like lines, curves and cirles describing color as well as relative size, position and shape of elements. When preparing the image for display is has to be translated into a bitmap image, this process is called rasterization.

A vector image is resolution independent, this means that you can enlarge or shrink the image without affecting the output quality. Vector images are the preferred way to represent Fonts, Logos and many illustrations.

When measuring the value for a pixel, one takes the average color of an area around the location of the pixel. A simplistic model is sampling a square, this is called a box filter, a more physically accurate measurement is to calculate a weighted Gaussian average (giving the value exactly at the pixel coordinates a high weight, and lower weight to the area around it). When perceiving a bitmap image the human eye should blend the pixel values together, recreating an illusion of the continuous image it represents.

Resolution

Resolution is a measurement of sampling density, resolution of bitmap images give a relationship between pixel dimensions and physical dimensions. The most often used measurement is ppi, pixels per inch ^[2].

Figure 1.4. Sampling grid

A rasterized form of the letter 'a' magnified 16 times, where each pixel is represented as a circle instead of a square.

Megapixels

Megapixels refer to the total number of pixels in the captured image, an easier metric is raster dimensions which represent the number of horizontal and vertical samples in the sampling grid. An image with a 4:3 aspect ratio with dimension 2048x1536 pixels, contain a total of 2048x1535=3,145,728 pixels; approximately 3 million, thus it is a 3 megapixel image.

Table 1.1. Common/encountered raster dimensions

Dimensions	Megapixels	Name	Comment
640x480	0.3	VGA	VGA
720x576	0.4	CCIR 601 DV PAL	Dimensions used for PAL DV, and PAL DVDs
768x576	0.4	CCIR 601 PAL full	PAL with square sampling grid ratio
800x600	0.4	SVGA
1024x768	0.8	XGA	The currently (2004) most common computer screen dimensions.
1280x960	1.2
1600x1200	2.1	UXGA
1920x1080	2.1	1080i HDTV	interlaced, high resolution digital TV format.
2048x1536	3.1	2K	Typically used for digital effects in feature films.
3008x1960	5.3
3088x2056	6.3
4064x2704	11.1

Scaling / Resampling

When we need to create an image with different dimensions from what we have we scale the image. A different name for scaling is resampling, when resampling algorithms try to reconstruct the original continous image and create a new sample grid.

Scaling image down

The process of reducing the raster dimensions is called decimation, this can be done by averaging the values of source pixels contributing to each output pixel.

Scaling image up

When we increase the image size we actually want to create sample points between the original sample points in the original raster, this is done by interpolation the values in the sample grid, effectivly guessing the values of the unknown pixels^[3].

Sample depth

The values of the pixels need to be stored in the computers memory, this means that in the end the data ultimately need to end up in a binary representation, the spatial continuity of the image is approximated by the spacing of the samples in the sample grid. The values we can represent for each pixel is determined by the sample format chosen.

Figure 1.5. Sample depth

The same image width varying sample depths, note that high frequency areas (detailed areas) have an OK look earlier than low frequency areas.

8bit

A common sample format is 8bit integers, 8bit integers can only represent 256 discrete values (2^8 = 256), thus brightness levels are quantized into these levels.

12bit

For high dynamic range images (images with detail both in shadows and highlights) 8bits 256 discrete values does not provide enough precision to store an accurate image. Some digital cameras operate with more than 8bit samples internally, higher end cameras (mostly SLRs) also provide RAW images that often are 12bit (2^12bit = 4096).

16bit

The PNG and TIF image formats supports 16bit samples, many image processing and manipulation programs perform their operations in 16bit when working on 8bit images to avoid quality loss in processing.

Floating point

Some image formats used in research and by the movie industry store floating point values. Both "normal" 32bit floating point values and a special format called half which uses 16bits/sample. Floating point is useful as a working format because quantization and computational errors are kept to a minimum until the final render.

Floating point representations often include HDR, High Dynamic Range. High Dynamic Range images are images that include sampling values that are whiter than white (higher values than 255 for a normal 8bit image). HDR allows representing the light in a scene with a greater degree of precision than LDR, Low Dynamic Range images.

Colors

The most common way to model color in Computer Graphics is the RGB color model, this corresponds to the way both CRT monitors and LCD screens/projectors reproduce color. Each pixel is represented by three values, the amount of red, green and blue. Thus an RGB color image will use three times as much memory as a gray-scle image of the same pixel dimensions.

Figure 1.6. RGB bands

Color image built up of bands of red, green and blue color (This is image illustrates how a laptop display is constructed, note that this image is preferred to be viewed on a computer screen, and not in print.)

One of the most common pixel formats used is 8bit rgb where the red, green and blue values are stored interleaved in memory. This memory layout is often referred to as chunky, storing the components in seperate buffers is called planar, and is not as common.

Palette / Indexed images

It was earlier common to store images in a palletized mode, this works similar to a paint by numbers strategy. We store just the number of the palette entry used for each pixel. And for each palette entry we store the amount of red, green and blue light.

Figure 1.7. Indexed image

On the left, image using just 16 colors, on the right the palette used for this image. The way an indexed/paletted image works is similar to how paint by numbers work.

Bitmap images take up a lot of memory, image compression reduces the amount of memory needed to store an image. For instance a 2.1 megapixel, 8bit RGB image (1600x1200) occupies 1600x1200x3 bytes = 5760000 bytes = 5.5 megabytes, this is the uncompressed size of the image.

Compression ratio is the ratio between the compressed image and the uncompressed image, if the example image mentioned above was stored as a 512kb jpeg file the compression ratio would be 0.5mb : 5.5mb = 1:11.

Lossless Image Compression

When an image is losslessly compressed, repetition and predictability is used to represent all the information using less memory. The original image can be restored. One of the simplest lossless image compression methods is run-length encoding. Run-length encoding encodes consecutive similar values as one token in a data stream.

Figure 1.8. Run-length encoding

70,
5, 25,
5, 27,
4, 26,
4, 25,
6, 24,
6, 23,
3, 2, 3, 22,
3, 2, 3, 21,
3, 5, 2, 20,
3, 5, 2, 19,
3, 7, 2, 18,
3, 7, 2, 17,
14, 16,
14, 15,
3, 11, 2, 14,
3, 11, 2, 13,
3, 13, 2, 12,
3, 13, 2, 11,
3, 15, 2, 10,
3, 15, 2, 8,
6, 12, 6, 6,
6, 12, 6, 64
           

In Figure 1.8, “Run-length encoding” a black and white image of a house has been compressed with run length encoding, the bitmap is considered as one long string of black/or white pixels, the encoding is how many bytes of the same color occur after each other. We'll further reduce the amount of bytes taken up by these 72 numerical values by having a maximum span length of 15, and encoding longer spans by using multiple spans separated by zero length spans of the other color.

70,                     15,  0, 15,  0, 15,  0, 10,
 5, 25,                  5, 15,  0, 10,
 5, 27,                  6, 15,  0, 12,
 4, 26,                  4, 15,  0, 11,
 4, 25,                  4, 15,  0, 10,
 6, 24,                  6, 15,  0,  9,
 6, 23,                  6, 15,  0,  8,
 3,  2, 3, 22,           3,  2,  3, 15,  0,  7,
 3,  2, 3, 21,           3,  2,  3, 15,  0,  6,
 3,  5, 2, 20,           3,  5,  2, 15,  0,  5,
 3,  5, 2, 19,           3,  5,  2, 15,  0,  4,
 3,  7, 2, 18,           3,  7,  2, 15,  0,  3,
 3,  7, 2, 17,           3,  7,  2, 15,  0,  2
14, 16,                 14, 15,  0,  1
14, 15,                 14, 15,
 3, 11, 2, 14,           3, 11,  2, 14,
 3, 11, 2, 13,           3, 11,  2, 13,
 3, 13, 2, 12,           3, 13,  2, 12,
 3, 13, 2, 11,           3, 13,  2, 11,
 3, 15, 2, 10,           3, 15,  2, 10,
 3, 15, 2,  8,           3, 15,  2,  8,
 6, 12, 6,  6,           6, 12,  6,  6,
 6, 12, 6, 64            6, 12,  6, 15, 0, 15, 0, 15, 0, 15, 0, 4
  

The new encoding is 113 nibbles long, a nibble i 4bit and can represent the value 0--4, thus we need 57 bytes to store all our values, which is less than the 93 bytes we would have needed to store the image as a 1bit image, and much less than the 750 bytes needed if we used a byte for each pixel. Run length encoding algorithms used in file formats would probably use additional means to compress the RLE stream achieved here.

Lossy Image Compression

Lossy image compression takes advantage of the human eyes ability to hide imperfection and the fact that some types of information are more important than others. Changes in luminance are for instance seen as more significant by a human observer than change in hue.

JPEG is a file format implementing compression based on the Discrete Cosine Transform DCT, together with lossless algorithms this provides good compression ratios. The way JPEG works is best suited for images with continuous tonal ranges like photographs, logos, scanned text and other images with lot's of sharp contours / lines will get more compression artifacts than photographs.

Loss through Generations

Lossy compression algorithms should not be used as a working format, only final copies should be saved as jpeg since loss accumulates over generations.

Figure 1.9. JPEG generation loss

An image specially constructed to show the deficiencies in the JPEG compression algorithm, saved, reopened and saved again 9 times.

JPEG is most suited for photographics content where the adverse effect of the compression algorithm is not so evident.

JPEG is not suited as an intermediate format, only use JPEG for final distribution where filesize actually matters.

A programming tutorial for lua is out of scope for this document, the samples given elsewhere in this document should provide a feel for the language, some peculiarities of lua is explained here, look at The lua website, for a short overview when reading sample code Lua Short Reference might come in handy.

There are functions to set and get individual pixels, they are of the form:

See the section called “Gluas reference” for more functions, and details.

All values used in gluas are floating point values, unless mandated by the color model, the normal range is from 0.0 to 1.0. This range is used instead of an more arbitary range like 0-255, 0-42, 0-100 or 0-65535.

When pixels are set, they don't immediatly overwrite the original pixel values. This is to facilitate area operations. When multipass algorithms are implemented, like a two pass blur, the changes can be comitted to the source buffer by issuing

Refer to Images Processing using gluas (which this gluas introduciton is a part of) for examples of how the gluas environment can be used.

Gluas uses the core lua language for processing, but extends it with some functions specific to interaction with The GIMP.

In order to know the dimensions, and some other aspects of the image, some numeric constants are also imported into the lua enviroment you program in.

As can be observed in the code examples given in this chapter, there is direct mapping between input pixel values and output pixel values for point process operations.

When dealing with 8bit images, the number of input and output values are greatly reduced, if the processing taking place within the transform is expensive, having a precalculated table for all input values is a common optimizations, this is called a LookUpTable (LUT).

When the sample values stored are 16bit the advantage of a lookup table is still present since the number of entries is 65536, but when moving to 32bit floating point, the number of precalculated values needed is exceedingly large.

Write a gluas program that exchanges the red and green component of all pixels in an image.

Is this a point operation?

Create a program performing the transformation in Figure 4.10, “The GIMP's curves tool”.

Create a program that draws a 8x8 checker board on top of the existing image. The checker board should be black and white, square, cover as large area as possible and be centered on the image.

Create a program that has a linearly increasing contrast applied to it. So that the leftmost column of the raster has a contrast of 0.0 applied to it, and the rightmost column of the raster has a contrast of 2.0 applied to it.

Is this a point operation?

When sampling a local area around a pixel, we need to decide what to do when part of the sampling window is outside the source image. The common way to provide a "fake" value for data outside the image are:

In gluas the global variable edge_duplicate can be set to 1, this makes gluas return the nearest valid pixel instead of transparent/black.

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.

A convolution matrix with all weights set to 1.0 is a blur, blurring is an operation that is often performed and can be implemented programmatically faster than the general convolution. A blur can be seen as the average of an area, or a resampling which we'll return to in a later chapter.

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)

A better looking blur is achieved by setting the weights to a gaussian distribution around the active pixel.

A gaussian blur works by weighting the input pixels near the center of ther sampling window higher than the input pixels further away.

One way to approximate a small gaussian blur, is to set the center pixel in a 3×3 matrix to 2, while the surrounding are set to 1, and set the weight to 10.

A blur is a time consuming operation especially if the radius is large. A way to optimize a blur is to observe that first blurring horizontally and then vertically gives the same result as doing both at the same time. This has been done to the blur_fast in Figure 5.4, “blur” and can also be done for gaussian blurs.

Edge preserving smoothing filters are useful for automated feature detection or even just for artistic effects.

Kuwahara

The kuwahara filter is an edge preserving blur filter. It works by calculating the mean and variance for four subquadrants, and chooses the mean value for the region with the smallest variance.

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)

Symmetric Nearest Neighbour

An improvement over the kuwahara filter is the symmetric nearest neighbour filter that compares symmetric pairs/quadruples of pixels with the center pixel and only takes into account the one from this set which is closest in value to the center pixel.

Figure 5.6. snn

-- Symmetric Nearest Neighbour as described by Hardwood et al in Pattern
-- recognition letters 6 (1987) pp155-162.

function snn(radius)
local x,y
for y=0, height-1 do
	for x=0, width-1 do
	    local sum=0.0
	    local count=0.0
	    local u,v
	    
	    local valueC = get_value (x,y)
	    
	    for v=-radius,radius do
	      for u=-radius,radius do
	         local valueA = get_value (x+u,y+v)
	         local valueB = get_value (x-u,y-v)
	         
	         if (math.abs(valueC-valueA) <
	             math.abs(valueC-valueB)) then
               sum = sum + valueA
	         else
               sum = sum + valueB	         
	         end
	         count = count + 1	         
	      end
	    end
		set_value (x,y, sum/count)
	end
	progress (y/height)
end
flush()
end

-- snn(6)

-- and a color version operating in CIE Lab 

function deltaE(l1,a1,b1,l2,a2,b2)
  return math.sqrt( (l1-l2)*(l1-l2) + (a1-a2)*(a1-a2) + (b1-b2)*(b1-b2))
end

function snn_color(radius)
local x,y
for y=0, height-1 do
	for x=0, width-1 do
	    local sumL=0.0
	    local suma=0.0
	    local sumb=0.0
	    local count=0.0
	    local u,v
	    
	    local lc,ac,bc = get_lab (x,y)
	    
	    for v=-radius,radius do
	      for u=-radius,radius do
			 local l1,a1,b1 = get_lab (x+u,y+v)	         
	         local l2,a2,b2 = get_lab (x-u,y-v)
	         
	         if (deltaE(lc,ac,bc,l1,a1,b1) <
	             deltaE(lc,ac,bc,l2,a2,b2)) then	             
               sumL = sumL + l1
               suma = suma + a1
               sumb = sumb + b1
	         else
	           sumL = sumL + l2
               suma = suma + a2
               sumb = sumb + b2	         
	         end
	         count = count + 1	         
	      end
	    end
		set_lab (x,y, sumL/count, suma/count, sumb/count)
	end
	progress (y/height)
end
flush()
end

snn_color(6)

The median filter sorts the sample values in theneighbourhood, and then picks the middle value. The effect is that noise is removed while detail is kept. See also kuwahara. The median filter in The GIMP is called Despeckle ^[5].

Figure 5.7. 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)

Pick an a random pixel within the sampling window.

Figure 5.8. 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)

Extend the convolve example, making it possible for it to accept a 7×7 matrix for it's values, add automatic summing to find the divisor and experiment with various sampling shapes.

Create a filter that returns the original pixel value if the variance of the entire sample window is larger than a provided value. How does this operation compare to kuwahara/snn?

Create a filter that renders a white circle in the middle of the image, the background of the circle should be black.

Translation can be described as moving the image, e.g. move all the contents of the image 50pixels to the right and 50pixels down.

Instead of traversing the original image, and placing the pixels in their new location we calculate which point in the source image ends up at the coordinate calculated.

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)

Rotation is implemented by rotating the unit vectors of the coordinate system, and from the unit vectors calculate the new positions of a given input coordinate, to make the transformation happen in reverse it is enough to negate the angle being passed to our function.

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))

Scaling an image is a matter of scaling both the x and y axes according to a given scale factor when retrieving sample values.

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)

Perspective transformations are useful for texture mapping, texture mapping is outside the scope of this text. Affine transformations are a subset of perspective transformations.

To be able to do texturemapping for use in realtime 3d graphics, the algorithm is usually combined with a routine that only paints values within the polygon defined, such a routine is called a polyfiller.

When stitching multiple images to create a panorama/ fake fish-eye image, correction for the lens is needed, the only difference compared to the earlier mentioned affine transformations are that the parameters used to calculate the distortion are based on properties of the lens.

Bilinear sampling makes the assumption that values in the image changes continously.

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)

With OpenGL and other realtime 3d engines it is common to use mipmaps to increase processing speed and the quality of the resampled texture.

A mipmap is a collection of images of prescaled sizes. The prescaled sizes are 1.0×image dimensions, 0.5×image dimensions, 0.25×image dimensions, 0.125×image dimensions etc.

When retrieving a sample value a bilinear interpolation would be performed at both the version of the image larger than the scaling ratio as well as the image smaller than the scaling ratio. After that a linear interpolation is performed on the two resulting values based on which of the two versions most closely matches the target resolution.

Using The GIMP open a test image, create a duplicate of the original image. (Image/Duplicate). Scale the duplicate image to 25% size (Image/Scale). Scale the duplicate image up 400%. Copy the newly scaled image as a layer on top of the original image. (Ctrl+C, change image Ctrl+V), and press the new layer button when the original image is active. Set the mode of the top layer to difference, Flatten the image (Image/Flatten Image), and automatically stretch the contrast of the image. (Layer/Colors/Auto/Stretch Contrast).

Do this for nearest neighbour, bilienar and bicubic interpolation and compare the results.

Read Appendix A, Assignment, spring 2007 and perhaps make initial attempts at demosaicing the Bayer pattern.

Assuming that we have a good distribution of colors in our scene, the average reflected color should be the color of the light. If the light source is assumed to be white, we know how much the whitepoint should be moved in the color cube.

The compensation calculated by the gray world assumption is an approximation of the measurement taken by digital still and video cameras by capturing an evenly lit white sheet of paper or similar.

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()

Component stretching makes the assumption that either direct reflections, or glossy reflections from surfaces can be found in the image, and that it is amongst the brightest colors in the image. By stretching each R,G and B component to it's full range, as in the section called “Contrast stretching” often leads to a better result than grayworld assumption.

If the image is overexposed, this technique does not work, not having a very reflective object will also bias the results in a non desireable way.

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()

ACE and Retinex are algorithms that performs color balancing and color enhancements baded on statistical measures of the spatial neighbourhood of each pixel; doing a kind of local contrast stretching. These algorithms aim to model some of the automatic adaptation that happens in the human visual system.

Figure 7.4. lce

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 lce_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 lce(radius)

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

lce(32)

	Note
	This implementation of an local contrast enhancement doesn't really show the potential of this type of algorithms but the implementation is kept simple to be readable and understandable.

The artifacts appearing in Figure 7.4, “lce” are due to the sampling of the pixels used to find the maximum and minimum values to do contrast stretching with. How can the filter be improved to give a smoother appearance?

Line doubling is the simplest form of deinterlacing, it replaces the odd lines with the even lines. This gets rid of the interlacing effect, but reduces the vertical resolution in half.

"Linear Blend operates by taking a line, then averaging the pixel values in it with the line below, effectively blurring the frame. This almost completely eliminates the effects of interlacing - at times you may notice slight ghosting (instant transition from one thing to another), the image appears to persist for a fraction of a second. This is absolutely not a serious issue, as far as I am concerned - in fact it's really only noticeable if you are me, and a perfectionist." [quote from http://xaw-deinterlace.sourceforge.net/index.php?sec=about]

By reconstructing motion vectors for the video frames, it is possible to predict the position of object between frames, thus providing a volume video object where any position can be sampled. Such an advanced approach can be used for both deinterlacing and high quality time stretching.