Selective Colouring with ImageMagick

Following a recent Photojojo newsletter, I thought I would do some experimentation with black and white. But rather than use Photoshop or GIMP, I decided to use ImageMagick: the command line may seem like more effort at first but it gives more control and you can also put everything in a script for later re-use.

Choosing a good photograph

There are some photographs that come out great as black and white, others that don't. Generally, the best candidates are the ones that look a bit uninteresting in colour, generally because there isn't much colour in them to start with. It's easy to get that sort of shots in London: just go out on a cloudy day and you'll get likely candidates. And so, I shot this picture on Sunday afternoon while on the South Bank:

The Globe Theatre

Although the building is nice, the picture itself is quite dull: there is little colour in it apart from the red flags so it's an ideal candidate for black and white treatment.

Separating the channels

Following the Photojojo article, rather than transform that photograph into grey scale, I started by separating the red, green and blue channels, which is very easy to do with ImageMagick:

$ convert img.jpg -separate img_rgb_%d.jpg

Assuming your original image is called img.jpg, this will produce three images called img_rgb_0.jpg, img_rgb_1.jpg and img_rgb_2.jpg that correspond to the red, green and blue channels respectively:

Red, green and blue channels

The last one, the blue channel, is the best black and white in my opinion as it's got the strongest shadows and lightest sky. Having said that, the red flags in the original photograph were the main spots of bright colour and it would be interesting to bring them back in the black and white shot through selective colouring. So let's do that.

Creating a mask

To selectively colour the flags red in the black and white picture, I am going to combine the original shot with the blue channel image. For this, I need to create a mask that will tell ImageMagick what part of the original shot to use and what part of the blue channel image to use. For the mask, I need a black and white picture that is white where the flags are and black everywhere else: this way, I can tell ImageMagick to use pixels from the original image where the mask is white and pixels from the blue channel image where the mask is black.

Because the flags are red, the idea is to find all the pixels in the original image that are red and no other colour. We can do that by combining the three channel images. Pixels that are a shade of grey, anywhere between black and white will show as the same shade of grey on all three channels but pixels that are of a specific colour, such as red, will be of a light grey in their primary channel and very dark in the other two. Indeed, if you look at the three channel images, the flags show as very light in the red channel but nearly black in the other two. So to create our mask, we can subtract the green and blue channels from the red one.

First, let's subtract the green channel from the red one:

$ convert img_rgb_1.jpg img_rgb_0.jpg -compose minus \
 -composite img_rgb_0-1.jpg

Which produces the following image:

Green channel subtracted from the red one

We can already see the flags pop out of the picture but we can also still see the outline of the theatre. So let's subtract the blue channel from the image we just obtained:

$ convert img_rgb_2.jpg img_rgb_0-1.jpg -compose minus \
 -composite img_rgb_0-1-2.jpg

Green and blue channels subtracted from the red one

That looks good, except that the flags are now light grey rather than really white. It is still a good mask but as the grey pixels will mix both images we use, it will dilute the effect. To ensure that we have a full effect, we need to stretch the contract of the mask to have the background area really black and the flags area really white.

$ convert img_rgb_0-1-2.jpg -level 10%,30% img_mask.jpg

The values given after the -level option specify the percentage of white under which a pixel is considered black first followed by the percentage of white over which a pixel is considered white. There is no real perfect value for those so you need to experiment. In this example, the second value is quite low because we really want most grey pixels to be turned white for a sharp mask.

Stretched mask

You will note that there are some small dots of white at the bottom of the mask. That's slightly unexpected but we'll see in a minute what they are.

Creating the final picture

Now that we have a mask, we need to combine the blue channel image with the original image by applying the mask. This is the simplest use of the -composite operator:

$ convert img_rgb_2.jpg img.jpg img_mask.jpg \
 -composite img_final.jpg

This produces the final image. We can now see that the small dots of white at the bottom of the mask have let out some of the colour of the red flowers on the lamppost, which is a nice extra touch and which would have been forgotten had we created the mask by hand with a drawing tool. Such side effects are typical of such a tool that works on the whole image and can be welcome or not depending on the situation. In any way, it is always possible to tweak the mask manually in a drawing tool before applying it.

Final image

Bootnote

The photo in its original size is now visible on flickr.

Fractals with Octave: Trigonometric and Exponential Functions

The story so far

• Classic Mandelbrot and Julia
• More on the Classic Sets
• Other Polynomials

In this fourth instalment of this series on fractals with Octave, I'll have a look at the Mandelbrot and Julia sets generated by non-polynomial series. Polynomials are fine but they can be a bit boring so what about introducing a sine or exponential in the mix? The original idea for this came from an excellent set of articles and examples by Paul Bourke.

Sine function

The first series we will use is defined thus:

z_n+1=c sin(z_n)

That's simple enough and very easy to code in Octave so as our mandelbrot and julia functions now have the ability to take a generating series as argument, it should be no problem. But before doing that, we need to know what value to give to r so that we know when the series diverges and when we can stop iterating. In the very first article, we saw that:

[...] for a given polynomial series, there exists a value r such that if for any n, |z_n|>r, then the series diverges.

This only applies to polynomials. The series above is not a polynomial series and its condition for divergence is that the imaginary part of z_n be greater than 50. That doesn't work well with our current code that expects a simple integer. We could replace the integer parameter expected by the mandelbrot and julia functions with a function handle but that would make their usage more complex when using polynomial series. Luckily, we can have the best of both worlds because Octave uses a weakly typed language: the type of a parameter passed to a function is not specified in the function definition, meaning that you can pass anything, you just have to hope that it is what the function expects. While this can be seen as a drawback for long and complex programs where validation can be cumbersome, it is actually a useful feature for a high level language like Octave where functions tend to be short. To make it work, we need to modify the mjcore function so that it can accept a function handle or an integer as its fifth parameter and act accordingly. For this, we will use the Octave isa built-in function to identify whether the parameter is a function handle. Here's how to modify the mjcore function:

mjcore.m

function M=mjcore(z,c,niter,f,r)
  if(isa(r,"function_handle"))
    rf=r;
  else
    rf=@(z) abs(z)<r;
  endif
  M=zeros(length(z(:,1)),length(z(1,:)));
  for s=1:niter
    mask=abs(z)<rrf(z);
    M(mask)=M(mask)+1;
    z(mask)=f(z(mask),c(mask));
  endfor
  M(mask)=0;
endfunction

That's all great and now we should be able to create our first sine Mandelbrot. We need to make sure that we reverse the comparison operator in the last parameter because our code is built so that we provide a convergence condition rather than a divergence one:

octave-3.0.1:1> Msin=mandelbrot(-1.2*pi+0.9*pi*i,1.2*pi-0.9*pi*i,
> 320,64,
> @(z,c) c.*sin(z), @(z) abs(imag(z))<=50);
octave-3.0.1:2> imagesc(Msin)

And we end up with... a dark blue rectangle, not quite what was expected. This is because in this form, the value of z₀ is 0, therefore sin(z₀) is also 0 and so is c sin(z₀): the series is always null and never diverges. To correct this, we would need to have z₀=c, that is initialise the series to the values of c. But if we do that in the code of the mandelbrot function, it will break some of our previous examples. That is unless we can specify a special value of z₀ to the mandelbrot function to force it to initialise to z₀=c. A similar trick to what we did for r above can do that: if the given parameter is the special string "c", initialise to c, otherwise assume the parameter is an integer and initialise as before. Here is the modified mandelbrot function:

mandelbrot.m

function M=mandelbrot(cmin,cmax,hpx,niter,
                      f=@(z,c) z.^2.+c,r=2,
                      z0=0)
  vpx=round(hpx*abs(imag(cmax-cmin)/real(cmax-cmin)));
  z=zeros(vpx,hpx).+z0;
  [cRe,cIm]=meshgrid(linspace(real(cmin),real(cmax),hpx),
                     linspace(imag(cmin),imag(cmax),vpx));
  c=cRe+i*cIm;
  if(strcmp("c",z0))
    z=c;
  else
    z=zeros(vpx,hpx).+z0;
  endif
  M=mjcore(z,c,niter,f,r);
endfunction

After that modification, we can run the mandelbrot function again with our new special parameter:

octave-3.0.1:1> Msin=mandelbrot(-1.2*pi+0.9*pi*i,1.2*pi-0.9*pi*i,
> 320,64,
> @(z,c) c.*sin(z), @(z) abs(imag(z))<=50, "c");
octave-3.0.1:2> imagesc(Msin)

Octave should open the Gnuplot window and display an image similar to this:

Sine Mandelbrot Set

We can also generate a sine Julia set without any change to the julia function:

octave-3.0.1:1> Jsin=julia(-2.4*pi+1.8*pi*i,2.4*pi-1.8*pi*i,
> 320,64,1+0.5i,
> @(z,c) c.*sin(z), @(z) abs(imag(z))<=50);
octave-3.0.1:2> imagesc(Jsin)

Sine Julia Set

Cool stuff! So let's have a look at other non-polynomial functions and their output.

Cosine

This series is defined by:

z_n+1=i c cos(z_n)

And has the same convergence condition than the previous series.

octave-3.0.1:1> Mcos=mandelbrot(-1.2*pi+0.9*pi*i,1.2*pi-0.9*pi*i,
> 320,64,
> @(z,c) i.*c.*cos(z), @(z) abs(imag(z))<=50, "c");
octave-3.0.1:2> imagesc(Mcos)

Cosine Mandelbrot Set

octave-3.0.1:1> Jcos=julia(-2.4*pi+1.8*pi*i,2.4*pi-1.8*pi*i,
> 320,64,0.55+1.195i,
> @(z,c) i.*c.*cos(z), @(z) abs(imag(z))<=50);
octave-3.0.1:2> imagesc(Jcos)

Cosine Julia Set

Tangent

This series is defined by:

z_n+1=c tan(z_n)

And uses the same convergence function as a polynomial with the default value of r.

octave-3.0.1:1> Mtan=mandelbrot(-1.4+1.4i,1.4-1.4i,
> 320,64,
> @(z,c) c.*tan(z), :, "c");
octave-3.0.1:2> imagesc(Mtan)

Tangent Mandelbrot Set

octave-3.0.1:1> Jtan=julia(-1.4+1.4i,1.4-1.4i,
> 320,64,(1+i)*sqrt(2)/2,
> @(z,c) c.*tan(z));
octave-3.0.1:2> imagesc(Jtan)

Tangent Julia Set

Exponential

First, let's multiply the exponential using the following series:

z_n+1=c exp(z_n²)

This series uses the default values for r and z₀.

octave-3.0.1:1> Mexp1=mandelbrot(-1.4+2i,1.4-2i,
> 240,64,
> @(z,c) c.*exp(z.^2));
octave-3.0.1:2> imagesc(Mexp1)

Exponential Mandelbrot Set #1

Then, let's add to the exponential using the following series:

z_n+1=exp(z_n²)+c

This series also uses the default values for r and z₀.

octave-3.0.1:1> Mexp2=mandelbrot(-2.2+2i,1-2i,
> 240,64,
> @(z,c) exp(z.^2).+c);
octave-3.0.1:2> imagesc(Mexp2)

Exponential Mandelbrot Set #2

Cactus

Finally, let's draw a cactus using the following series:

z_n+1=z_n³+(z₀-1)z_n-z₀

This series uses the default values for r and the special value z₀=c.

octave-3.0.1:1> Mcactus=mandelbrot(-0.8+0.6i,0.8-0.6i,
> 320,64,
> @(z,c) z.^3.+(c.-1).*z.-c,:,"c");
octave-3.0.1:2> imagesc(Mcactus)

Cactus Mandelbrot Set

Next

Next in this series, we will look at a fractal called the burning ship.

Lloyds TSB phishing scam

I just received an email which claims to come from Lloyds TSB. As I am not a customer of the bank, it is obviously a scam. I've been looking on the bank's web site in an attempt to report it to them but there doesn't seem to be any email address advertised where you can do so. There are a couple of telephone numbers for their customers but that's it. Surely, it would be useful to them to have copies of those emails, if only to know they exist and to warn their customers. But it seems they are not interested. Maybe it's a good thing I don't bank with them then.

Photos from Peru

I came back from Peru a week ago and spent a good part of last week tweaking and cleaning the best of them. They are now visible on flickr. The tweaking and cleaning I do on photos is very simple and consists in:

• stretching the contrast if required;
• removing fluff and dirt mark left on the lens or the sensor;
• very occasionally, cropping the shot.

So no heavy Photoshop (or in that instance GIMP) processing, just the bare minimum.

It's all relative

I'm in Peru. We just spent a couple of days around the Cañon de Colca, one of the deepest canyons in the world. The place where we stayed last night, Chivay, is at 3700 metres above sea level and on the way there we stopped at a viewpoint that was 4900 metres above sea level, higher than Mont Blanc. The area around Chivay and the canyon is lush and green and if it wasn't for the thin air, you'd easily forget that you are at an altitude that is higher than most of the Alps' ski resorts.

Nobody told the locals either, as explained by one of the girls in the group this morning: she asked the landlady of the hotel where she stayed overnight whether she kept llamas or alpacas, to which the landlady replied "we don't do that here, only in the mountains". That'd be the steep bits above 4000 metres that we see in the distance then.

Bootnote

This was written on the bus between Chivay and Puno yesterday and uploaded today.

ExifTool on Ubuntu

A while ago, I blogged about installing ExifTool on Ubuntu. There's actually a much simpler way to do this. ExifTool is part of the standard packages directly available from the Ubuntu repositories so it can be installed in one line using apt-get, no need for that make malarkey I mentioned last time:

$ sudo apt-get install libimage-exiftool-perl

Et voilà, ExifTool installed and ready to go!

Making Skype work with PulseAudio in Ubuntu Intrepid and Jaunty

Following my sound issues with Skype on Ubuntu, I did a bit more research and eventually found an excellent how-to article on PulseAudio in the Ubuntu forums. Appendix C explains how to get Skype to work properly and indeed it is a doodle and means there is no need to kill PulseAudio when using Skype! I've tested it in Intrepid and Jaunty and it works a treat. To quote the article, here's how to do it:

Open Skype's Options, then go to Sound Devices. You need to set "Sound Out" and "Ringing" to the "pulse" device, and set "Sound In" to the hardware definition of your microphone. For example, my laptop's microphone is defined as "plughw:I82801DBICH4,0".

You will have to experiment with the different options for "Sound In" until you find the correct one: choose an option, click on "Make a test call" until you find the option that works. On my machine, here's what the option window looks like:

Skype Options, Sound Devices Tab

Google Street View

Like a number of people, I was a bit dubious about Google Street View. If you want to see what a place looks like, there are enough photographs on the web for you to find, why would you want to go look at bad photographs taken by a passing car? However, where Street View does come into its own is when coupled with driving directions in Google Maps: for each intersection, you can have a picture of what it looks like with an arrow that tells you where to go. This is a fantastic visual aid when going to a place you are not familiar with.

Driving directions with Street View

Adding a Static Address to the DHCP Server

I've got a Lexmark X9575 all-in-one printer. Up to now, it was connected to the Mac via USB. I wanted to reconfigure it to connect it to the network instead so that I can use it with other computers. But I wanted to do so in such a way that it would get its network configuration via DHCP while keeping a fixed address and be resolvable via DNS.

The first thing I did was to remove the printer setup from the Mac (no, you can't change the configuration, you have to un-install and re-install, thank you Lexmark!). Then I connected the printer to the network switch with a standard RJ45 cable. It requested a DHCP lease form the server immediately and the DNS entries got updated. That was a good start. However, the IP address was part of the dynamic range and therefore unpredictable and the name of the printer was hard-coded to the very non-intuitive ET0020002C5B7A. I wanted something more memorable like lexmark-printer.

To force a fixed IP address is quite simple, as explained in this DHCP mini-howto. I just added the following to the /etc/dhcp3/dhcpd.conf file:

host lexmark-printer {
 hardware ethernet 00:20:00:2c:5b:7a;
 fixed-address 192.168.1.250;
}

Restarted the DHCP server:

$ sudo /etc/init.d/dhcp3-server restart

And it didn't work, for two reasons:

• The DHCP server already had a lease assigned so wasn't going to renew it;
• The printer already has an IP address configured so wasn't going to renew it either.

To solve the first problem, I needed to force the expiry of the lease on the DHCP server. The ISC DHCP server stores leases between restarts in a file called /var/lib/dhcp3/dhcpd.leases. So revoking the lease can be done by updating that file, removing the entry for the given lease and restarting the DHCP server.

The second problem was of the same ilk: go to the printer's TCP/IP configuration, set DHCP to no, save the settings, go back into them, set DHCP back to yes and that forced the printer to request a new lease.

That gave me a static address for the printer while still keeping it managed by DHCP. However, in such a case, the DHCP server doesn't update the DNS database. This has to be done manually, by adding the following in the forward lookup file:

lexmark-printer.home.   IN A    192.168.1.250

And this in the reverse lookup file:

250.1.168.192.in-addr.arpa. IN PTR  lexmark-printer.home.

Then of course, a restart of the DNS server is required:

$ sudo /etc/init.d/bind9 restart

Skype Sound Issues

Skype is not available in the Ubuntu repositories but you can get hold of it (and a few other things) through medibuntu. However, it looks like that version of Skype has a problem with sound on Ubuntu 8.10 Intrepid Ibex, as detailed in this bug report. The workaround suggested in the discussion that consists in killing pulseaudio before launching Skype works for me so that's good news. Even better news, it looks like the problem is fully sorted in Ubuntu 9.04, which is currently in beta. So, until then, I will use the workaround.

One lesson to learn from this though, is that if you ever have a problem with a piece of software on Linux, a good practice is to start that software from the command line: the output in the terminal window is invaluable for developers and maintainers to understand the problem and is a good place to start in order to find a solution to it. Quite often, just putting the error message in Google will return a number of articles and bug reports on that exact problem.