Charlie Harvey

Seven More Languages: Julia day three

The final day of Julia is a lot more challenging than the first 2 — consisting of a larger example of some image processing and a little bit about macros along with some wrap up material and some challenging exercises.

Throughout the whole thing Julia struck me as a very usable language. Intuitive, nicely structured and with just enough syntax to make it very pleasant to use. The quote from Graydon Hoare is very aposite — he calls Julia a Goldilocks language because the designers have got it just right

Exercises

Easy exercises

Write a macro that runs a block of code backwards

Not quite sure what running code backwards really means here. I just deconstruct an expression, reverse its arguments and map eval over that. I suppose that is what was required.

julia> macro backwards(b) quote map(eval,reverse($b.args)) return end end julia> @backwards :(begin println("one") println("two") end) two one

cameraman image contrasty

Experiment with modifying frequencies and observing the effect on an image. What happens when you set some high frequencies to large values?

Scaling up the higher numbers makes the image background much less detailed, effectively it ends up being white as you can see from the picture. My implementation is below.

function blockdct6(img) pixels = convert(Array{Float32}, img.data) y,x = size(pixels) outx = ifloor(x/8) outy = ifloor(y/8) bx = 1:8:outx*8 by = 1:8:outy*8 mask = zeros(8,8) mask[1:3,1:3] = [1 1 1; 1 1 0; 1 0 0] ## keep stuff marked 1, drop stuff that is 0 freqs = Array(Float32, (outy*8,outx*8)) for i=bx, j=by freqs[j:j+7, i:i+7] = dct(pixels[j:j+7, i:i+7]) freqs[j:j+7, i:i+7] .*= mask end map(scaleupHighNumber, freqs) end function scaleupHighNumber(n) if n > 2 n * 144 else n end end

What happens if you add lots of noise? (Hint: try adding scale * rand(size(freqs))

As you can see, adding noise makes the image look significantly blockier.

text image blocky

Here is how I did that.

function blockdct6(img) pixels = convert(Array{Float32}, img.data) y,x = size(pixels) outx = ifloor(x/8) outy = ifloor(y/8) bx = 1:8:outx*8 by = 1:8:outy*8 mask = zeros(8,8) mask[1:3,1:3] = [1 1 1; 1 1 0; 1 0 0] ## keep stuff marked 1, drop stuff that is 0 freqs = Array(Float32, (outy*8,outx*8)) for i=bx, j=by freqs[j:j+7, i:i+7] = dct(pixels[j:j+7, i:i+7]) freqs[j:j+7, i:i+7] .*= mask end freqs .* rand(size(freqs)) end

Medium exercises

Modify the code to allow masking arbitrarily many coefficient, but always the N most important ones. Instead of calling blockdct6(img), you would call blockdct(img,6)

In implementing this one, I made the assumption that arbitrarily many would be bounded by our block size, so up a maximum of to 8x8 == 64 cells could be masked by calling make_mask(16) as we are working with 8x8 blocks. The algorithim could remain substantially the same, if you wanted to support different sized blocks.

I make a subsidiary function to deal with the mask making, it takes the width of the block being processed and an n for the number of "significant" coefficients.

function make_mask(width,n) row_max = n arr=zeros(width,width) println(arr) for i=1:row_max for j=1:row_max if(j<n-(i+1)) arr[i,j]=1 end end end arr end function blockdct(img, n) pixels = convert(Array{Float32}, img.data) y,x = size(pixels) outx = ifloor(x/8) outy = ifloor(y/8) bx = 1:8:outx*8 by = 1:8:outy*8 mask = make_mask(8,n) freqs = Array(Float32, (outy*8,outx*8)) for i=bx, j=by freqs[j:j+7, i:i+7] = dct(pixels[j:j+7, i:i+7]) freqs[j:j+7, i:i+7] .*= mask end freqs end

Our codec outputs a frequency array as big as its input, even though most frequencies are zero. Instead, output only the non-zero frequencies for each block so that the output is smaller than the input. Modify the decoder to use the smaller input as well.

I spent a fair while trying to implement this myself, then in the process of reading the manual came across sparse and full which basically did what was neded. So in the spirit of not doing difficult shit when there is a library funtion to do it.

function smaller_blockdct(img,n) sparse(blockdct(img,n)) end function smaller_blockidct(freqs,n) blockidct(full(freqs)) end

Experiment with different block sizes to see how block size affects the appearance of coding artefacts. Try a large block size on an image containing lots of text and see what happens.

Everything gets a lot blockier! I parameterized the original functions on blocksize first of all. Now the code read as follows.

function blockdct(img, n, bs) pixels = convert(Array{Float32}, img.data) y,x = size(pixels) outx = ifloor(x/bs) outy = ifloor(y/bs) bx = 1:bs:outx*bs by = 1:bs:outy*bs mask = make_mask(bs,n) freqs = Array(Float32, (outy*bs,outx*bs)) for i=bx, j=by freqs[j:j+bs-1, i:i+bs-1] = dct(pixels[j:j+bs-1, i:i+bs-1]) freqs[j:j+bs-1, i:i+bs-1] .*= mask end freqs end function blockidct(freqs,bs) y,x = size(freqs) bx = 1:bs:x by = 1:bs:y pixels = Array(Float32, size(freqs)) for i=bx, j=by pixels[j:j+bs-1, i:i+bs-1] = idct(freqs[j:j+bs-1, i:i+bs-1]) end grayim(pixels) end

Hard exercises

The code currently only works on grayscale images, but the same technique works on colour too. Modify the code to work on colour images like testimage("mandrill")

This was actually a pretty straightforward refactor.

For RGB images, the Images.jl library stores each colour component as a grey-like array so you go from having a 2d representation of an image to a 3d one with the additional 3rd dimension representing colour. Which means that processing colour images is just running the grey code for each colour dimension.

I was able to hack a working implementation together at breakfast one day. Here’s my final code for the whole codec. THe colour parts are the rgb_ functions.

module Codec using Images, Colors, ImageView function smaller_blockdct(img,n) sparse(blockdct(img,n)) end function smaller_blockidct(freqs,n) blockidct(full(freqs)) end function make_mask(width,n) row_max = n arr=zeros(width,width) for i=1:row_max for j=1:row_max if(j<n-(i+1)) arr[i,j]=1 end end end arr end function rgb_blockdct (img, n, bs) freqs = convert(Array{Float32}, data(separate(img))) y,x,cdim = size(freqs) for i in 1:cdim freqs[1:y,1:x,i] = inner_blockdct(freqs[1:y,1:x,i],n,bs) end freqs end function rgb_blockidct (freqs, bs) y,x,cdim = size(freqs) for i in 1:cdim freqs[1:y,1:x,i] = blockidct(freqs[1:y,1:x,i],bs) end convert(Image{RGB}, freqs) end function blockdct(img, n, bs) pixels = convert(Array{Float32}, img.data) inner_blockdct(pixels, n, bs) end function blockdct(img, n) blockdct(img, n, 8) end function blockdct6(img) blockdct(img, 6) end # this is the actual work part that calls dct (we crack out this part as it is called # per colo[u]r dimension for colo[u]r images function inner_blockdct(pixels, n, bs) y,x = size(pixels) outx = ifloor(x/bs) outy = ifloor(y/bs) bx = 1:bs:outx*bs by = 1:bs:outy*bs mask = make_mask(bs,n) freqs = Array(Float32, (outy*bs,outx*bs)) for i=bx, j=by freqs[j:j+bs-1, i:i+bs-1] = dct(pixels[j:j+bs-1, i:i+bs-1]) freqs[j:j+bs-1, i:i+bs-1] .*= mask end freqs end function blockidct(freqs,bs) y,x = size(freqs) bx = 1:bs:x by = 1:bs:y pixels = Array(Float32, size(freqs)) for i=bx, j=by pixels[j:j+bs-1, i:i+bs-1] = idct(freqs[j:j+bs-1, i:i+bs-1]) end grayim(pixels) end function blockidct(freqs) blockidct(freqs,8) end function scaleupHighNumber(n) if n > 2 n * 144 else n end end end # //module

As you can see I also refactored some of the other functions to simplify the whole codec.

JPEG does prediction of the first coefficient, called the DC offset. The previous block’s DC value is subtracted from the current block’s DC value. This encodes an offset instead of a number with full range, saving valuable bits. Try implementing this in code.

Hmmm. I think this is an exercise to return to. I will need to re-read the wikipedia article and probably dig round some implementations to understand exactly what it is asking. I am going to leave it for another day and move on.

Wrapping up

So, what have we learned?

  • Julia is a really nice language that gets a lot of things right
  • My blog’s syntax highlighting can’t quite deal with its syntax ;-)
  • Image processing can be fun and relatively performant.

Julia is certainly a language to look out for in the future. The syntax is lovely, its macros are nice, the package system is fantastic and it is very pleasant to use.

On the down side it is still relatively new, so there needs to be some effort put into getting the libraries you need to get shit done. That is just a matter of time.

I have really enjoyed the language. It felt comfortable to hack in, like an old pair of boots and I can imagine myself using it for a variety of tasks, especially for jobs where R feels sluggish.


Comments

  • Be respectful. You may want to read the comment guidelines before posting.
  • You can use Markdown syntax to format your comments. You can only use level 5 and 6 headings.
  • You can add class="your language" to code blocks to help highlight.js highlight them correctly.

Privacy note: This form will forward your IP address, user agent and referrer to the Akismet, StopForumSpam and Botscout spam filtering services. I don’t log these details. Those services will. I do log everything you type into the form. Full privacy statement.