3D Hilbert Curve
restart;
with(plots):
Arithmetisation of a 3D Hilbert Curve
Definition of the operators H0 to H7 (compare Fig. 8.6, right plot):
H[0] := (x,y,z) -> ( x/2, z/2, y/2);
H[1] := (x,y,z) -> ( z/2, y/2+1/2, x/2);
H[2] := (x,y,z) -> ( z/2+1/2, y/2+1/2, x/2);
H[3] := (x,y,z) -> (-y/2+1, -x/2+1/2, z/2);
H[4] := (x,y,z) -> (-y/2+1, -x/2+1/2, z/2+1/2);
H[5] := (x,y,z) -> (-z/2+1, y/2+1/2,-x/2+1);
H[6] := (x,y,z) -> (-z/2+1/2, y/2+1/2,-x/2+1);
H[7] := (x,y,z) -> ( x/2, -z/2+1/2,-y/2+1);
Zio2JUkieEc2IkkieUdGJUkiekdGJUYlNiRJKW9wZXJhdG9yR0YlSSZhcnJvd0dGJUYlNiUsJDkkIyIiIiIiIywkOSZGLiwkOSVGLkYlRiVGJQ==
Zio2JUkieEc2IkkieUdGJUkiekdGJUYlNiRJKW9wZXJhdG9yR0YlSSZhcnJvd0dGJUYlNiUsJDkmIyIiIiIiIywmOSVGLkYuRi8sJDkkRi5GJUYlRiU=
Zio2JUkieEc2IkkieUdGJUkiekdGJUYlNiRJKW9wZXJhdG9yR0YlSSZhcnJvd0dGJUYlNiUsJjkmIyIiIiIiI0YuRi8sJjklRi5GLkYvLCQ5JEYuRiVGJUYl
Zio2JUkieEc2IkkieUdGJUkiekdGJUYlNiRJKW9wZXJhdG9yR0YlSSZhcnJvd0dGJUYlNiUsJjklIyEiIiIiIyIiIkYxLCY5JEYuI0YxRjBGMSwkOSZGNEYlRiVGJQ==
Zio2JUkieEc2IkkieUdGJUkiekdGJUYlNiRJKW9wZXJhdG9yR0YlSSZhcnJvd0dGJUYlNiUsJjklIyEiIiIiIyIiIkYxLCY5JEYuI0YxRjBGMSwmOSZGNEY0RjFGJUYlRiU=
Zio2JUkieEc2IkkieUdGJUkiekdGJUYlNiRJKW9wZXJhdG9yR0YlSSZhcnJvd0dGJUYlNiUsJjkmIyEiIiIiIyIiIkYxLCY5JSNGMUYwRjRGMSwmOSRGLkYxRjFGJUYlRiU=
Zio2JUkieEc2IkkieUdGJUkiekdGJUYlNiRJKW9wZXJhdG9yR0YlSSZhcnJvd0dGJUYlNiUsJjkmIyEiIiIiIyMiIiJGMEYyLCY5JUYxRjFGMiwmOSRGLkYyRjJGJUYlRiU=
Zio2JUkieEc2IkkieUdGJUkiekdGJUYlNiRJKW9wZXJhdG9yR0YlSSZhcnJvd0dGJUYlNiUsJDkkIyIiIiIiIywmOSYjISIiRjBGLkYvLCY5JUYzRi9GL0YlRiVGJQ==
Standard algorithm to compute the resepective mapping:
hilbert3D := proc(t::numeric, depth::integer)
# t: paramter of the Hilbert mapping
# depth: number of considered octal digits of t
global H;
local q,r;
if depth = 0 then return (0,0,0) end if;
q := floor(8*t); # extract the first octal digit
r := 8*t-q; # and compute remainder
return H[q]( hilbert3D(r,depth-1) );
end proc:
hilbert3D(1/2,5);
NiUiIiIjRiMiIiNGJA==
hilbert3D(1/8,9);
NiUiIiEjIiIiIiIjRiM=
hilbert3D(1/3,5);
NiUjIiNEIiNLIyIjSkYlIyIiJCIiKQ==
We will need a parameter that is approximately mapped to the centre (1/2, 1/2, 1/2) of the unit cube.
For the given curve, such a parameter is not easy to find:
hilbert3D( 10193/3^9,10);evalf(%);
NiUjIiRKIiIkYyMjIiQ0JiIlQzUjIiNMIiNr
NiUkIisrdj08XiEjNSQiK0RKcXFcRiUkIisrK0RjXkYl
Plotting the approximating polygons
corners := pointplot3d([[0,0,0],[0,0,1],[0,1,0],[1,0,0]]):
iter := 1:
curve :=
spacecurve( [ seq( [hilbert3D(i/8^iter,8)], i=0..8^iter-1), [0,0,1] ],
axes=BOXED, scaling=CONSTRAINED, thickness=3):
display3d(corners,curve);
LSUnUExPVDNERzYsLSUnUE9JTlRTRzYmNyUkIiIhISIiJCIiISEiIiQiIiEhIiI3JSQiIiEhIiIkIiIhISIiJCIjNSEiIjclJCIiISEiIiQiIzUhIiIkIiIhISIiNyUkIiM1ISIiJCIiISEiIiQiIiEhIiItJSdDVVJWRVNHNiQ3KzclJCIiISEiIiQiIiEhIiIkIiIhISIiNyUkIiIhISIiJCIiJiEiIiQiIiEhIiI3JSQiIiYhIiIkIiImISIiJCIiISEiIjclJCIjNSEiIiQiIiYhIiIkIiIhISIiNyUkIiM1ISIiJCIiJiEiIiQiIiYhIiI3JSQiIzUhIiIkIiImISIiJCIjNSEiIjclJCIiJiEiIiQiIiYhIiIkIiM1ISIiNyUkIiIhISIiJCIiJiEiIiQiIzUhIiI3JSQiIiEhIiIkIiIhISIiJCIjNSEiIi0lKlRISUNLTkVTU0c2IyIiJC0mJSZfQVhJU0c2IyIiIjYmLSUmQ09MT1JHNiYlJFJHQkckIiIhISIiJCIiISEiIiQiIiEhIiItJSpMSU5FU1RZTEVHNiMiIiEtJSpUSElDS05FU1NHNiMiIiEtJS1UUkFOU1BBUkVOQ1lHNiMkIiIhISIiLSYlJl9BWElTRzYjIiIjNiYtJSZDT0xPUkc2JiUkUkdCRyQiIiEhIiIkIiIhISIiJCIiISEiIi0lKkxJTkVTVFlMRUc2IyIiIS0lKlRISUNLTkVTU0c2IyIiIS0lLVRSQU5TUEFSRU5DWUc2IyQiIiEhIiItJiUmX0FYSVNHNiMiIiQ2Ji0lJkNPTE9SRzYmJSRSR0JHJCIiISEiIiQiIiEhIiIkIiIhISIiLSUqTElORVNUWUxFRzYjIiIhLSUqVEhJQ0tORVNTRzYjIiIhLSUtVFJBTlNQQVJFTkNZRzYjJCIiISEiIi0lK0xJR0hUTU9ERUxHNiNRKExJR0hUXzM2Ii0lK1BST0pFQ1RJT05HNiwkIikoKioqKipcISIpJCEpeDFycSEiKSQiKSgqKioqKlwhIikkIikoKioqKipcISIpJCIpeDFycSEiKSQiKSgqKioqKlwhIikkISl4MXJxISIpJCIiISEiIiQiKXgxcnEhIikkIiM1ISIiLSUqQVhFU1NUWUxFRzYjJSRCT1hHLSUoU0NBTElOR0c2IyUsQ09OU1RSQUlORURHLSUlUk9PVEc2Jy0lKUJPVU5EU19YRzYjJCIkKyIhIiItJSlCT1VORFNfWUc2IyQiJCsiISIiLSUtQk9VTkRTX1dJRFRIRzYjJCIlK1IhIiItJS5CT1VORFNfSEVJR0hURzYjJCIlK0ghIiItJSlDSElMRFJFTkc2Ig==
Plotting the iterations
Note that, as h(1/2) is not the centre of the cube, connecting the images of the interval centres will look strange:
iter := 1;
curve :=
spacecurve( [ seq( [hilbert3D((i+1/2)/8^iter,8)], i=0..8^iter-1) ],
scaling=CONSTRAINED, axes=BOXED, thickness=4):
display3d(corners,curve);
IiIi
LSUnUExPVDNERzYsLSUnUE9JTlRTRzYmNyUkIiIhISIiJCIiISEiIiQiIiEhIiI3JSQiIiEhIiIkIiIhISIiJCIjNSEiIjclJCIiISEiIiQiIzUhIiIkIiIhISIiNyUkIiM1ISIiJCIiISEiIiQiIiEhIiItJSdDVVJWRVNHNiQ3KjclJCIiJiEiIiQiI0QhIiMkIiNEISIjNyUkIiNEISIjJCIjdiEiIyQiIiYhIiI3JSQiI3YhIiMkIiN2ISIjJCIiJiEiIjclJCIjdiEiIyQiIiEhIiIkIiNEISIjNyUkIiN2ISIjJCIiISEiIiQiI3YhIiM3JSQiI3YhIiMkIiN2ISIjJCIiJiEiIjclJCIjRCEiIyQiI3YhIiMkIiImISIiNyUkIiImISIiJCIjRCEiIyQiI3YhIiMtJSpUSElDS05FU1NHNiMiIiUtJiUmX0FYSVNHNiMiIiI2Ji0lJkNPTE9SRzYmJSRSR0JHJCIiISEiIiQiIiEhIiIkIiIhISIiLSUqTElORVNUWUxFRzYjIiIhLSUqVEhJQ0tORVNTRzYjIiIhLSUtVFJBTlNQQVJFTkNZRzYjJCIiISEiIi0mJSZfQVhJU0c2IyIiIzYmLSUmQ09MT1JHNiYlJFJHQkckIiIhISIiJCIiISEiIiQiIiEhIiItJSpMSU5FU1RZTEVHNiMiIiEtJSpUSElDS05FU1NHNiMiIiEtJS1UUkFOU1BBUkVOQ1lHNiMkIiIhISIiLSYlJl9BWElTRzYjIiIkNiYtJSZDT0xPUkc2JiUkUkdCRyQiIiEhIiIkIiIhISIiJCIiISEiIi0lKkxJTkVTVFlMRUc2IyIiIS0lKlRISUNLTkVTU0c2IyIiIS0lLVRSQU5TUEFSRU5DWUc2IyQiIiEhIiItJStMSUdIVE1PREVMRzYjUShMSUdIVF8zNiItJStQUk9KRUNUSU9ORzYsJCIpKCoqKioqXCEiKSQhKXgxcnEhIikkIikoKioqKipcISIpJCIpKCoqKioqXCEiKSQiKXgxcnEhIikkIikoKioqKipcISIpJCEpeDFycSEiKSQiIiEhIiIkIil4MXJxISIpJCIjNSEiIi0lKkFYRVNTVFlMRUc2IyUkQk9YRy0lKFNDQUxJTkdHNiMlLENPTlNUUkFJTkVERy0lJVJPT1RHNictJSlCT1VORFNfWEc2IyQiJCsiISIiLSUpQk9VTkRTX1lHNiMkIiQrIiEiIi0lLUJPVU5EU19XSURUSEc2IyQiJStSISIiLSUuQk9VTkRTX0hFSUdIVEc2IyQiJStSISIiLSUpQ0hJTERSRU5HNiI=
iter := 2;curve :=
spacecurve( [ seq( [hilbert3D((i+10193/3^9)/8^iter,8)], i=0..8^iter-1) ],
scaling=CONSTRAINED, axes=BOXED, thickness=4):
display3d(corners,curve);
IiIj
LSUnUExPVDNERzYsLSUnUE9JTlRTRzYmNyUkIiIhISIiJCIiISEiIiQiIiEhIiI3JSQiIiEhIiIkIiIhISIiJCIjNSEiIjclJCIiISEiIiQiIzUhIiIkIiIhISIiNyUkIiM1ISIiJCIiISEiIiQiIiEhIiItJSdDVVJWRVNHNiQ3XG83JSQiKUQxKkciISIpJCIkRCIhIiQkIiREIiEiJDclJCIkRCIhIiQkIilEMSpHIiEiKSQiJHYkISIkNyUkIiR2JCEiJCQiKUQxKkciISIpJCIkdiQhIiQ3JSQiJHYkISIkJCIkRCIhIiQkIil2JDRAIiEiKTclJCIkdiQhIiQkIiR2JCEiJCQiKXYkNEAiISIpNyUkIiR2JCEiJCQiKXYkNHIkISIpJCIkdiQhIiQ3JSQiJEQiISIkJCIpdiQ0ciQhIikkIiR2JCEiJDclJCIpRDEqRyIhIikkIiR2JCEiJCQiJEQiISIkNyUkIiREIiEiJCQiJEQnISIkJCIpRDEqRyIhIik3JSQiKUQxKkciISIpJCIkdikhIiQkIiREIiEiJDclJCIpRDEqRyIhIikkIiR2KSEiJCQiJHYkISIkNyUkIiREIiEiJCQiKXYkNEAnISIpJCIkdiQhIiQ3JSQiJHYkISIkJCIpdiQ0QCchIikkIiR2JCEiJDclJCIpdiQ0ciQhIikkIiR2KSEiJCQiJHYkISIkNyUkIil2JDRyJCEiKSQiJHYpISIkJCIkRCIhIiQ3JSQiJHYkISIkJCIkRCchIiQkIilEMSpHIiEiKTclJCIkRCchIiQkIiREJyEiJCQiKUQxKkciISIpNyUkIilEMSpHJyEiKSQiJHYpISIkJCIkRCIhIiQ3JSQiKUQxKkcnISIpJCIkdikhIiQkIiR2JCEiJDclJCIkRCchIiQkIil2JDRAJyEiKSQiJHYkISIkNyUkIiR2KSEiJCQiKXYkNEAnISIpJCIkdiQhIiQ3JSQiKXYkNHIpISIpJCIkdikhIiQkIiR2JCEiJDclJCIpdiQ0cikhIikkIiR2KSEiJCQiJEQiISIkNyUkIiR2KSEiJCQiJEQnISIkJCIpRDEqRyIhIik3JSQiJHYpISIkJCIpdiQ0ciQhIikkIiREIiEiJDclJCIkRCchIiQkIiR2JCEiJCQiKUQxKkciISIpNyUkIiREJyEiJCQiJEQiISIkJCIpRDEqRyIhIik3JSQiKUQxKnkpISIpJCIkRCIhIiQkIiREIiEiJDclJCIpRDEqeSkhIikkIiREIiEiJCQiJHYkISIkNyUkIiREJyEiJCQiJEQiISIkJCIpdiQ0ciQhIik3JSQiJEQnISIkJCIkdiQhIiQkIil2JDRyJCEiKTclJCIkdikhIiQkIil2JDRyJCEiKSQiJHYkISIkNyUkIiR2KSEiJCQiKXYkNHIkISIpJCIkRCchIiQ3JSQiJEQnISIkJCIkdiQhIiQkIilEMSpHJyEiKTclJCIkRCchIiQkIiREIiEiJCQiKUQxKkcnISIpNyUkIilEMSp5KSEiKSQiJEQiISIkJCIkRCchIiQ3JSQiKUQxKnkpISIpJCIkRCIhIiQkIiR2KSEiJDclJCIkRCchIiQkIiREIiEiJCQiKXYkNHIpISIpNyUkIiREJyEiJCQiJHYkISIkJCIpdiQ0cikhIik3JSQiJHYpISIkJCIpdiQ0ciQhIikkIiR2KSEiJDclJCIkdikhIiQkIiREJyEiJCQiKXYkNHIpISIpNyUkIil2JDRyKSEiKSQiJHYpISIkJCIkdikhIiQ3JSQiKXYkNHIpISIpJCIkdikhIiQkIiREJyEiJDclJCIkdikhIiQkIil2JDRAJyEiKSQiJEQnISIkNyUkIiREJyEiJCQiKXYkNEAnISIpJCIkRCchIiQ3JSQiKUQxKkcnISIpJCIkdikhIiQkIiREJyEiJDclJCIpRDEqRychIikkIiR2KSEiJCQiJHYpISIkNyUkIiREJyEiJCQiJEQnISIkJCIpdiQ0cikhIik3JSQiJHYkISIkJCIkRCchIiQkIil2JDRyKSEiKTclJCIpdiQ0ciQhIikkIiR2KSEiJCQiJHYpISIkNyUkIil2JDRyJCEiKSQiJHYpISIkJCIkRCchIiQ3JSQiJHYkISIkJCIpdiQ0QCchIikkIiREJyEiJDclJCIkRCIhIiQkIil2JDRAJyEiKSQiJEQnISIkNyUkIilEMSpHIiEiKSQiJHYpISIkJCIkRCchIiQ3JSQiKUQxKkciISIpJCIkdikhIiQkIiR2KSEiJDclJCIkRCIhIiQkIiREJyEiJCQiKXYkNHIpISIpNyUkIilEMSpHIiEiKSQiJHYkISIkJCIkdikhIiQ3JSQiJEQiISIkJCIpdiQ0ciQhIikkIiREJyEiJDclJCIkdiQhIiQkIil2JDRyJCEiKSQiJEQnISIkNyUkIiR2JCEiJCQiJHYkISIkJCIpRDEqeSkhIik3JSQiJHYkISIkJCIkRCIhIiQkIilEMSp5KSEiKTclJCIkdiQhIiQkIilEMSpHIiEiKSQiJEQnISIkNyUkIiREIiEiJCQiKUQxKkciISIpJCIkRCchIiQ3JSQiKUQxKkciISIpJCIkRCIhIiQkIiR2KSEiJC0lKlRISUNLTkVTU0c2IyIiJS0mJSZfQVhJU0c2IyIiIjYmLSUmQ09MT1JHNiYlJFJHQkckIiIhISIiJCIiISEiIiQiIiEhIiItJSpMSU5FU1RZTEVHNiMiIiEtJSpUSElDS05FU1NHNiMiIiEtJS1UUkFOU1BBUkVOQ1lHNiMkIiIhISIiLSYlJl9BWElTRzYjIiIjNiYtJSZDT0xPUkc2JiUkUkdCRyQiIiEhIiIkIiIhISIiJCIiISEiIi0lKkxJTkVTVFlMRUc2IyIiIS0lKlRISUNLTkVTU0c2IyIiIS0lLVRSQU5TUEFSRU5DWUc2IyQiIiEhIiItJiUmX0FYSVNHNiMiIiQ2Ji0lJkNPTE9SRzYmJSRSR0JHJCIiISEiIiQiIiEhIiIkIiIhISIiLSUqTElORVNUWUxFRzYjIiIhLSUqVEhJQ0tORVNTRzYjIiIhLSUtVFJBTlNQQVJFTkNZRzYjJCIiISEiIi0lK0xJR0hUTU9ERUxHNiNRKExJR0hUXzM2Ii0lK1BST0pFQ1RJT05HNiwkIigtOmoiISIoJCIoUm8ieSEiKCQiKSd5Jj5nISIpJCEoejBvIyEiKCQiKWQrQmkhIikkISgpXGF0ISIoJCEpYiFcXCohIikkISh2cDglISIpJCIpTWk1SiEiKSQiIzUhIiItJSpBWEVTU1RZTEVHNiMlJEJPWEctJShTQ0FMSU5HRzYjJSxDT05TVFJBSU5FREctJSVST09URzYnLSUpQk9VTkRTX1hHNiMkIiQrIiEiIi0lKUJPVU5EU19ZRzYjJCIkKyIhIiItJS1CT1VORFNfV0lEVEhHNiMkIiUrUiEiIi0lLkJPVU5EU19IRUlHSFRHNiMkIiUrUiEiIi0lKUNISUxEUkVORzYi
Here's a small test function to find a parameter that is mapped to the cube's centre:
test := proc()
local i,pt;
for i from 1 to 3^9-1 do
pt := hilbert3D(i/3^9,16):
if abs(pt[1]-1/2)<1/50 and abs(pt[2]-1/2)<1/50 and abs(pt[3]-1/2)<1/50
then print(i,pt,evalf(pt));
fi;
od;
end proc:
test();
NikiJSFcKiMiJSQ0IyIlJzQlIyIlWiIpIiYlUTsjIiQkKioiJVs/JCIrIkdqKTReISM1JCIrIT1NRChcRi8kIis3R2pbW0Yv
NikiJiQ+NSMiJSQ0IyIlJzQlIyIlWiIpIiYlUTsjIiZmUCQiJk9iJyQiKyJHaik0XiEjNSQiKyE9TUQoXEYvJCIrK1lAXl5GLw==