Simple image processing very slow

I have compared the speed of some simple image processing routines between Mathematica an IDL.

1. Reading of a grey scale png image of (720,577) pixels:

The image is: Mathematica:

t1 = AbsoluteTime[];
img = Import[filename];
Print[AbsoluteTime[] - t1]
0.029513

IDL:

t1=systime(1)
print,systime(1)-t1
0.013546944

IDL is giving the byte array of pixel values of the image.

To do that in addition with Mathematica:

t1 = AbsoluteTime[];
imgData = ImageData[img, "Byte"];
Print[AbsoluteTime[] - t1]
0.002133

costs another few msec.

---> IDL is about 2 times faster

2. Finding positions of pixels at which the image pixels are exceeding a threshold of 50:

Mathematica:

t1 = AbsoluteTime[];
is = Position[imgData, n_ /; n >= 50];
Print[AbsoluteTime[] - t1]
0.222408

IDL:

t1=systime(1)
is=where(img GE 50)
print,systime(1)-t1
0.00085401535

---> IDL is more than 250 times faster

3. Replacing the pixel values that exceed the threshold of 50 with 255:

Mathematica:

t1 = AbsoluteTime[];
img = ReplacePart[imgData, is -> 255];
Print[AbsoluteTime[] - t1]
0.080377

IDL:

t1=systime(1)
img(is)=255
print,systime(1)-t1
3.2901764e-05

---> IDL is nearly 2500 times faster

These differences are really disappointing.

I would be happy if somebody has a solution how to do the things faster.

• If you are working with grayscale, then you process each pixel three times for RGB - same with replace. – Yves Klett Dec 15 '14 at 11:59
• For finding the positions this should be quite fast: SparseArray[UnitStep[imgData - 50]]["NonzeroPositions"] – Simon Woods Dec 15 '14 at 12:47
• ImageApply[Min[{#, 255}] &, img] // AbsoluteTiming can do the last two quickly. – C. E. Dec 15 '14 at 13:39
• I think you should change the title/description, to be honest. Neither Import, Position or ReplacePart are image processing routines. You should rather be asking, "how can I do this with image processing routines?" or "how can I make this code faster?" – C. E. Dec 15 '14 at 16:43
• Thank you for this hint and for your help ... – mrz Dec 15 '14 at 21:17

All three parts of this operation can be done in 0.026 seconds:

In:= AbsoluteTiming[ImageAdd[img, Binarize[img = Import[filename, "PNG"], {50/255, 1}]];]
Out= {0.026024, Null}

PixelValuePositions can be used for extracting pixel positions:

In:= img = Import["http://i.stack.imgur.com/uye1v.png"];
AbsoluteTiming[bimg = Binarize[img, {50/255, 1}];]
AbsoluteTiming[pos = PixelValuePositions[bimg, 1];]
Out= {0.001001, Null}
Out= {0.004005, Null}
• I think this only addresses two out of the three elements: reading the image and scaling the pixel values. It does not address the retrieval of the pixel positions. – Mac Dec 16 '14 at 10:38
• Pixel positions are not necessary for the problem of replacing pixel values that exceed the threshold of 50 with 255. Pixel positions can be retrieved with PixelValuePositions. – Piotr Wendykier Dec 16 '14 at 14:12
• You could amend your answer to address and benchmark the position part as well. – Yves Klett Dec 16 '14 at 14:20

As for part three of your operation, this seems to be rather zippy (about 800x faster than the ReplacePart line alone, with the additional benefit of shedding the Position part):

img = Clip[imgData, {0, 49}, {0, 255}];
• Dear Yves, I have a MacBook Pro (Retina, Mid 2012), 2,7 GHz Intel Core i7, 16 GB 1600 MHz DDR3, NVIDIA GeForce GT 650M 1024 MB. This is great with "Clip", but I also would need for the further analysis the positions of the pixels above the threshold ... is there also a trick ... – mrz Dec 15 '14 at 9:59
• You might want to look at PixelValuePositions and ReplacePixelValue and related for a really performant solution. – Yves Klett Dec 15 '14 at 11:51

We focus on retrieving pixel positions meeting a criteria using Yves Klett and Simon Wood's answer. Let's run a test case using pixel values of 1 (this can easily be generalised to conditions such as pixel value >= 50). First generate some random coordinates in a 900 x 900 image.

xc = RandomSample[Range, 500];
yc = RandomSample[Range, 500];

Then generate a black image and colour the random points in white (i.e. set to 1).

img2 = ReplacePart[ConstantArray[0, {900, 900}], Transpose[{xc, yc}] -> 1];

Now time the various options.

t1 = Timing[ind1 = Position[img2, pix_ /; pix == 1];] // First
t2 = Timing[ind2 = PixelValuePositions[Image[img2], 1];] // First
t3 = Timing[ind3 = SparseArray[img2]["NonzeroPositions"];] // First
ind2 = Map[{901 - #[], #[]} &, Reverse[ind2, {2}]];

Notice that the image processing routine PixelValuePositions[] used image based coordinates rather than array based as for the other routines, hence the extra step to compute ind2.

First we check that all answers are now the same (the coordinates of the white pixels in array coordinates).

ind1 == ind2 == ind3

which gives

True

and then the timing improvement on my MacBookPro.

t1/{t2, t3}

giving

{28.96, 171.9}

So if the main information of interest is pixel positions you can get back about a factor of 30 in speed using the dedicated image processing routines and a factor of 170 using the SparseArray technique (which I believe is undocumented). Considering that PixelValuePositions[] is a dedicated function for this purpose there is definitely room for improvement here.

• The only thing to add is the range for values >=50. Could probably be done with the third argument for PixelValuePositions. – Yves Klett Dec 15 '14 at 14:10