Friday, December 2, 2011

The sieve of Eratosthenes



Finding prime numbers is a good exercise to improve your skills in the programming scripting language of ImageJ. Moreover, that's a good opportunity to use a 2D image as a temporary array.


 One of the oldest algorithm is the sieve of Erathosthenes [Wikipedia] allowing to find prime numbers up to an integer number num. The principle consists of filling an array from 2 to num and for each number, to search and to remove all the multiples. Thus, the 'survivors' are the prime numbers.
Here, rather than using the built-in function newArray(), we'll try to use a 32-bit image.

1- Data initialization

 The data initialization is done in two steps:
  1. Creation of a 32-bit image using newImage(...). For sake of convenience, the width is set to 10 and the height is equal to num/10 where num is the max number.
  2. Data are filled with setPixel(...)

newImage("prime numbers", "32-bit Black", 10, height, 1);
// Fill data from 0 to num

for (i=0;i<num;i++)
{
   x = i%10;  y = floor(i/10); // Extract X and Y coordinates from the index
   setPixel(x,y,i);
}

Better than using this loop, the same job can be done in one single line of code with Process > Math >Macro... as shown in the script (at the end of this post).

2- Loop

The core of this script comprises a nested loop. For each ith non-zero pixel value, scan the array from i+1 to num, and check if the remainder of the division by the ith pixel value is equal to 0 (zero). If yes, the pixel isn't a prime number and it sets to 0 (ie removed).

3- The script

After executing the script, the image "prime numbers" is created and it only remains the prime numbers (the non-prime numbers were set to 0 (black)) .
Fig. 1: The prime numbers correspond to the non-zero pixels in the image.

In conclusion, keep in mind that using images as data storage (arrays) are really powerful and very convenient specially for numbers.

Here is the whole macro/script. Just modify the variable 'num' corresponding to the maximum number (by default, equal to 217).

+++ IJ snippet +++ +++End of IJ snippet +++

A small optimization ... you can replace the inner loop by the following line:

run("Macro...","code=[if (x+w*y >"+primeNumber+ "&& v%"+primeNumber+"==0) v=0;]");

Another optimization, if you want to see the evolution of your image like in Wikipedia, add the two following lines just before the if statement if (getPixel(x,y)!=0)...

run("Select All");run("Copy");
run("Add Slice");run("Paste");

... and the result is ...
Fig. 2: Evolution of the array contents according to the 14 first numbers

No comments:

Post a Comment