.ad
.po -0.4i
.ll +1.1i
.he ''\s16Linear Pixel Shuffling Made Easy\s0''
.2c
.EQ
delim $$
gsize 12
.EN
.SZ12
.sz12
.sh1 Abstract
.pp
We describe how to generate the coordinates of
pixels in an $N times N$ image so that the first
$M$ pixels are evenly spread across the whole 
image, as $M$ grows from 1 to $N sup 2$.
.sh1 "A Useful Sequence"
.pp
We need to make heavy use of the following Fibonacci-like
sequence, $G sub n$, defined as follows:
.EQ
G sub n ^^=^^
left { ^^
lpile { 0 above 1 above 1 above G sub n-1 ^^+^^ G sub n-3 } ^^^^^^
lpile { if^^n^=^0 above if^^n^=^1 above if^^n^=^2 above if^^n^>=^3 }
.EN
The first few terms of this sequence are:
.EQ
0^^1^^1^^1^^2^^3^^4^^6^^9^^13^^19^^28^^41^^60^^88
.EN
.EQ
129^^189^^277^^406^^595^^872^^1278^^1873
.EN
Notice that the growth rate for every other term
is roughly a factor of two, so there are twice as
many $G sub n$'s as there are powers of two.
.sh1 "Enumerating Pixels"
.pp
The usual order of enumerating the pixel
coordinates in an $N times N$ square is to
``visit'' the pixels in the top line, left to
right, then the pixels in the next line, left to
right, and so on.  This is shown in the following
program fragment:
.TS
center,box;
l.
.nf
    I = (1,0);
    J = (0,1);
    for( x = 0; x < N; x++ )
        for( y = 0; y < N; y++ )
            process pixel x*I+y*J;    
.TE
We have expressed the ``process'' line of this code
using a linear combination of two vectors; but this means
the same thing as
.TS
center,box;
l.
.nf
    process pixel (x,y)
.TE
But the vector notation expresses the part we
will change.
.sh1 "Linear Pixel Shuffling"
.pp
We modify the basis vectors, $I$ and $J$, so that
the enumeration of the combinations,
$x^I^^+^^y^J$ still eventually visits all $N sup
2$ pixels, but does so in a manner that evenly
visits all over the image.  The following four
rules show how to change the program:
.np
Choose $N^^=^^G sub n$ as appropriate.
.np
Replace $I$ with
$( G sub n-1 sup 2 ,$ $ G sub n-2 sup 2 ^+^ G sub n-1 sup 2 )$.
.np
Replace $J$ with $(G sub n-2 ,$ $ G sub n-3 )$.
.np
Assure that any generated subscripts are legal
by remaindering them with $N$.
.sh1 "An Example"
.pp
Suppose you have an image that fits within a $129 times 129$ square.
.np
Let $N^^=^^129$.
.np
Let $I^^=^^( 88 sup 2 ,^^ 60 sup 2 ^^+^^ 88 sup 2 )$ 
= $( 7744,$ $ 11344)$
= $( 4,$ $ 121 )$.
(The last expression shows the values remaindered by 129.)
.np
Let $J^^=^^( 60,$ $ 41 )$.
.lp
The program may be transformed to the following:
.TS
center,box;
l.
.nf
   N = 129
   I1 = 4; I2 = 121;
   J1 = 60; J2 = 41;
   for( x = 0; x < N; x++ )
      for( y = 0; y < N; y++ )    
         process pixel
            ( (x*I1+y*J1)%N, (x*I2+y*J2)%N )  
.TE
The expression in the final line of this program can be maintained by
simple additions and conditional subtractions rather than multiplications
and remainderings (i.e., division), because the factors $x$ and $y$
are simple for loop induction variables.
....0 1 1 1 2 3 4 6 9 13 19 28 41 60 88 129 189 277 406 595 872 1278 1873
.lp
Here is the full program.
Two items need to be filled in by the user:
(1) $n$, such that $G[n]^^=^^N$ is an appropriate size
for the image, and
(2) application-specific
code ``process pixel $(i,j)$''.
.TS
center,box;
l.
.nf
   int N, I1, I2, J1, J2;
   int x, y, i, j;
   int G[100];
   
   G[0] = 0; G[1] = 1; G[2] = 1;
   for( i = 3; i <= n; i++ )
      G[i] = G[i-1] + G[i-3];
    
   N = G[n];
   I1 = (G[n - 1] * G[n - 1]) % N;
   I2 = (G[n - 2] * G[n - 2] + I1) % N;   
   J1 = G[n - 2];
   J2 = G[n - 3];
   
   for (x = 0; x < N; x++)
      for (y = 0; y < N; y++) {
         i = (x * I1 + y * J1) % N;
         j = (x * I2 + y * J2) % N;

         ...process pixel (i,j)...

        }
.TE
.sh1 "The Infinite Case."
.pp
Next, we address the problem of generating
an unlimited number of points in the unit cube in
$N$-space; i.e., a set of $N$-tuples
.EQ
(x sub 1 ,
x sub 2 ,
x sub 3 , ... ,
x sub N )
.EN
where $0^<=x sub k ^<1$, for each $k$, and the points
are as smoothly distributed as possible.
.pp
Determine the unique positive root of the equation
.EQ
x sup n+1 ~=~ x sup n ~+~ 1
.EN
using standard numerical techniques.
Let $a sub 1 ~=~ x^-^1$, and
$a sub k ~=~ a sub 0 ( a sub 0 ^+^ 1 ) sup k-1 $
for $k~=~2,^ ... ,^ N$.
.pp
Then, the $m$-th point in the $N$-cube is made up of the
fractional parts of the multiples of the $a sub k$;
that is,
.EQ
p sub k ~=~ (
"{"k a sub 1 "}",~
"{"k a sub 2 "}",~
"{"k a sub 3 "}",~ ... ,^
"{"k a sub N "}")
.EN
The fractional part notation is: $"{"y"}"$ = $y^-^[y]$.
parts of the