#Octave 234 Bytes

I=imread(input(''));w=rows(I);X=[0,3;1,2];for k=2:log2(w);n=numel(X);X=[X',rot90(X',2)+3*n;X+n,X+2*n];end;for k = 1:3;I(2:2:end,:,k)=fliplr(I(2:2:end,:,k));end[~,S]=sort(X(:));I(S+(0:w^2:2*w^2))=permute(I,[2 1 3]);imwrite(I,input(''))

File names of input and output images should be provided form standard input. size of code without input/output is 194 bytes.
Explanation:

The base pattern of indices is:

X =
 0 3
 1 2

In each iteration 4 copies from result of the previous iteration made and some transformation applied to each copy then all blocks concatenated to form the current result.

X =[0,3;1,2];
for k = 2:log2(s)
 n=numel(X);
 X = [X',rot90(X',2)+3*n;X+n,X+2*n];
end

so we have:

block(1,1): X' 
block(1,2): rot90(X',2)+3*n 
block(2,1): X+n
block(2,2): X+2*n
0 1 | 14 15
3 2 | 13 12
--------|--------
4 7 | 8 11
5 6 | 9 10

Hilbert indices sorted and indexes of sorted elements returned:

[~,S]=sort(X(:));

Unraveling applied flipping all even rows:

for k = 1:3
 I(2:2:end,:,k) = fliplr(I(2:2:end,:,k));
end

Reraveling applied :
-S repeted for each channel
-permutation applied since in Octave data arranged column wise

I(S+(0:w^2:2*w^2))=permute(I,[2 1 3]);

Example images:

enter image description here

enter image description here

enter image description here

enter image description here

Octave 234 Bytes

I=imread(input(''));w=rows(I);X=[0,3;1,2];for k=2:log2(w);n=numel(X);X=[X',rot90(X',2)+3*n;X+n,X+2*n];end;for k = 1:3;I(2:2:end,:,k)=fliplr(I(2:2:end,:,k));end[~,S]=sort(X(:));I(S+(0:w^2:2*w^2))=permute(I,[2 1 3]);imwrite(I,input(''))

File names of input and output images should be provided form standard input. size of code without input/output is 194 bytes.
Explanation:

The base pattern of indices is:

X =
 0 3
 1 2

In each iteration 4 copies from result of the previous iteration made and some transformation applied to each copy then all blocks concatenated to form the current result.

X =[0,3;1,2];
for k = 2:log2(s)
 n=numel(X);
 X = [X',rot90(X',2)+3*n;X+n,X+2*n];
end

so we have:

block(1,1): X' 
block(1,2): rot90(X',2)+3*n 
block(2,1): X+n
block(2,2): X+2*n
0 1 | 14 15
3 2 | 13 12
--------|--------
4 7 | 8 11
5 6 | 9 10

Hilbert indices sorted and indexes of sorted elements returned:

[~,S]=sort(X(:));

Unraveling applied flipping all even rows:

for k = 1:3
 I(2:2:end,:,k) = fliplr(I(2:2:end,:,k));
end

Reraveling applied :
-S repeted for each channel
-permutation applied since in Octave data arranged column wise

I(S+(0:w^2:2*w^2))=permute(I,[2 1 3]);

Example images:

enter image description here

enter image description here

enter image description here

enter image description here

#Octave 234 Bytes

I=imread(input(''));w=rows(I);X=[0,3;1,2];for k=2:log2(w);n=numel(X);X=[X',rot90(X',2)+3*n;X+n,X+2*n];end;for k = 1:3;I(2:2:end,:,k)=fliplr(I(2:2:end,:,k));end[~,S]=sort(X(:));I(S+(0:w^2:2*w^2))=permute(I,[2 1 3]);imwrite(I,input(''))

File names of input and output images should be provided form standard input. size of code without input/output is 194 bytes.
Explanation:

The base pattern of indices is:

X =
 0 3
 1 2

In each iteration 4 copies from result of the previous iteration made and some transformation applied to each copy then all blocks concatenated to form the current result.

X =[0,3;1,2];
for k = 2:log2(s)
 n=numel(X);
 X = [X',rot90(X',2)+3*n;X+n,X+2*n];
end

so we have:

block(1,1): X' 
block(1,2): rot90(X',2)+3*n 
block(2,1): X+n
block(2,2): X+2*n
0 1 | 14 15
3 2 | 13 12
--------|--------
4 7 | 8 11
5 6 | 9 10

Hilbert indices sorted ans indexes of sorted elements returned:

[~,S]=sort(X(:));

Unraveling applied flipping all even rows:

for k = 1:3
 I(2:2:end,:,k) = fliplr(I(2:2:end,:,k));
end

Reraveling applied :
-S repeted for each channel
-permutation applied since in Octave data arranged column wise

I(S+(0:w^2:2*w^2))=permute(I,[2 1 3]);

Example images:

enter image description here

enter image description here

enter image description here

enter image description here

#Octave 234 Bytes

I=imread(input(''));w=rows(I);X=[0,3;1,2];for k=2:log2(w);n=numel(X);X=[X',rot90(X',2)+3*n;X+n,X+2*n];end;for k = 1:3;I(2:2:end,:,k)=fliplr(I(2:2:end,:,k));end[~,S]=sort(X(:));I(S+(0:w^2:2*w^2))=permute(I,[2 1 3]);imwrite(I,input(''))

File names of input and output images should be provided form standard input. size of code without input/output is 194 bytes.
Explanation:

The base pattern of indices is:

X =
 0 3
 1 2

In each iteration 4 copies from result of the previous iteration made and some transformation applied to each copy then all blocks concatenated to form the current result.

X =[0,3;1,2];
for k = 2:log2(s)
 n=numel(X);
 X = [X',rot90(X',2)+3*n;X+n,X+2*n];
end

so we have:

block(1,1): X' 
block(1,2): rot90(X',2)+3*n 
block(2,1): X+n
block(2,2): X+2*n
0 1 | 14 15
3 2 | 13 12
--------|--------
4 7 | 8 11
5 6 | 9 10

Hilbert indices sorted and indexes of sorted elements returned:

[~,S]=sort(X(:));

Unraveling applied flipping all even rows:

for k = 1:3
 I(2:2:end,:,k) = fliplr(I(2:2:end,:,k));
end

Reraveling applied :
-S repeted for each channel
-permutation applied since in Octave data arranged column wise

I(S+(0:w^2:2*w^2))=permute(I,[2 1 3]);

Example images:

enter image description here

enter image description here

enter image description here

enter image description here