Auto-contrast: an optimization example where Matlab beats Matlab

For this post, we are glad to have Thierry as a guest blogger. He talks about an optimisation example he recently encountered that is a good indirect example of ‘inlining‘.

As my first post here, I will first describe a standard “auto-contrast” method with some optimizing hints. I will then try a more powerful one, and hit a situation where the buit-in Matlab function (quantile) can be beaten by a simple Matlab script.

Suppose you want to “auto-contrast” an image.
A simple way to do this is to scale the value of your matrix from 0 to 255 (or [0..2^16-1] for a 16 bit image). You can make it with this code :

Of course, this code makes more sense if you start from a matrix of double and want to store it in an image file (TIFF, png…). I am sure many of you know that the “imagesc” Matlab function is doing it for you, in the background.

With such a formula, you are sensitive to “hot pixels” or dust in optic which are out of scale compared to the rest of your data. To further improve the result, you therefore probably want to “saturate” both the highest (hot) pixels values as well as the lowest. This is what the “clims” parameter of imagesc is for. In ImageJ, there is a toot that “enhances contrast with XX% saturation“.

To perform the same calculation in Matlab, we just have to modify the min and max value in the previous formula by the XX% lowest and highest values. We might get these value with this algorithm :

sortedValues=sort(myImage(:));
N=length(sortedValues);
satMin=sortedValues(floor(N*saturation));
satMax=sortedValues(floor(N*(1-saturation)));
[and now, you use the previous autocontrast/scaling formula with this "satMin" and satMax values]

saturation is here a number between 0 and 1. I often choose a saturation of 0.001, which means that 0.1% of all my pixels will have the 0 value and the same amount will have 255.
Please note again the “(:)” to linearize myImage, allowing to sort all the value at the same time.

The above code happened to be what I produced as my first (after some debugging) guess. Using the profiler, I got a noticeable long calculation time on the “sort(” line. So I searched for a dedicated Matlab function (which in principle should have been the first thing to do). And guess what?

First thing to be proud of: I reproduce a Matlab function with no bug!

This is the now the equivalent code using that quantile function.

quantile(myImage(:),[saturation, 1-saturation])

So now, I wanted to compare the calculation time of the two variants. I used tic…toc functions, and found that my code was actually twice faster than the “quantile” built-in function.

I also found that the quantile function has a strange scaling with various input parameters, in term of computational speed, as shown below.

tic;
for i=1:220,quantile(I(:),.9);end;
toc

elapsed time : 11sec

tic;
for i=1:220,quantile(I(:),[.1 .9]);end;
toc

elapsed time : 19 sec

My small algorithm with “sort” does the job in 6 seconds, no matter how many values I want to extract (sorting the image is the bottleneck; extracting the “n-th” value is instantaneous).

To conclude, I used the profiler to find the line which takes the longest to compute. Usually, one can improve the speed by a factor (between 2 and 10) by removing bottlenecks. It was not possible here as even the builtin Matlab function was slower. The lesson to this is the following : if the slowest thing in your code is a basic Matlab function like matrix multiplication, convolution or (in this example) sort, there is usually nothing to do.

I’m sure you’ve already encountered such a situation, feel free to write it down here as a comment. As you can see, by writing simpler function, you can even beat Matlab provided function with a Matlab script of your own.