How can I improve perfomace of Hilbert scan of image?

This method of image scanning based on Hilbert Curve. Curve looks like (from 1 to 6 order):

It can be used for image scan. So, for example, my code for 3-order curve is:

Hilbert=[C(1,1) C(1,2) C(2,2) C(2,1) C(3,1) C(4,1) C(4,2) C(3,2) C(3,3) C(4,3) C(4,4) C(3,4)...
      C(2,4) C(2,3) C(1,3) C(1,4) C(1,5) C(2,5) C(2,6) C(1,6) C(1,7) C(1,8) C(2,8) C(2,7)...
      C(3,7) C(3,8) C(4,8) C(4,7) C(4,6) C(3,6) C(3,5) C(4,5) C(5,5) C(6,5) C(6,6) C(5,6)...
      C(5,7) C(5,8) C(6,8) C(6,7) C(7,7) C(7,8) C(8,8) C(8,7) C(8,6) C(7,6) C(7,5) C(8,5)...
      C(8,4) C(8,3) C(7,3) C(7,4) C(6,4) C(5,4) C(5,3) C(6,3) C(6,2) C(5,2) C(5,1) C(6,1)...
      C(7,1) C(7,2) C(8,2) C(8,1)];

And it works and works fast. I made the same functions for 8- and 9-order curves, but it works very-very slow. 9-order, perhaps, will never end. At least, I have not had the patience to wait for the end - after 2 hours I just turned off the program. But the 7-order curve runs for 15 seconds. What's the matter? Can I do the same, but faster? Yes, the program needs to read 512 * 512 array elements, but it can not be impossible to make it faster.

So, what exactly I need - I have the coordinates of the array elements, and they are arranged in the order in which should be read. I need for acceptable time to read them and write in the new array. How to do it?

p.s. English is still hard for me, if something is unclear - ask me, please.

Doing a quick search online, you can find a post about Hilbert curves on Steve Eddins blog. Here is his implementation to generate the curve:

function [x,y] = hilbert_curve(order)
    A = zeros(0,2);
    B = zeros(0,2);
    C = zeros(0,2);
    D = zeros(0,2);

    north = [ 0  1];
    east  = [ 1  0];
    south = [ 0 -1];
    west  = [-1  0];

    for i=1:order
        AA = [B ; north ; A ; east  ; A ; south ; C];
        BB = [A ; east  ; B ; north ; B ; west  ; D];
        CC = [D ; west  ; C ; south ; C ; east  ; A];
        DD = [C ; south ; D ; west  ; D ; north ; B];

        A = AA;
        B = BB;
        C = CC;
        D = DD;
    end

    subs = [0 0; cumsum(A)] + 1;
    x = subs(:,1);
    y = subs(:,2);
end

The returned xy-coordinates are integers in the range [1,2^order]. As you can see below, the function is sufficiently fast:

>> for order=1:10, tic, [x,y] = hilbert_curve(order); toc; end
Elapsed time is 0.001478 seconds.
Elapsed time is 0.000603 seconds.
Elapsed time is 0.000228 seconds.
Elapsed time is 0.000251 seconds.
Elapsed time is 0.000361 seconds.
Elapsed time is 0.000623 seconds.
Elapsed time is 0.001288 seconds.
Elapsed time is 0.007269 seconds.
Elapsed time is 0.029440 seconds.
Elapsed time is 0.117728 seconds.

Now let's test it with an image with the curve overlayed. We resize the image down to 128x128 so that we can see the pattern without getting overcrowded, but you can definitely do 512x512 for your case:

%// some grayscale square image
img = imread('cameraman.tif');

%// scale it down for better visualization
N = 128
assert(N > 0 && mod(N,2)==0);
img = imresize(img, [N N]);

%// space-filling Hilbert curve
order = log2(N)
[x,y] = hilbert_curve(order);

%// show image with curve overlayed
imshow(img, 'InitialMagnification',400)
h = line(x, y);

Let's zoom in a bit to better see how the curve covers all the pixels:

>> zoom(10)
>> set(h, 'Marker','.')

Finally you can use the subscripts to index into the image matrix:

>> ind = sub2ind([N N], x, y);
>> pix = img(ind);    %// linear indexing

where:

>> whos ind
  Name          Size             Bytes  Class     Attributes

  ind       16384x1             131072  double