RandToCircle
// convert 16-bit random value in [0,0xffff] into a position in the unit circle
// based on the approximate equal area transform from
// http://marc-b-reynolds.github.io/math/2017/01/08/SquareDisc.html
float2 RandToCircle(uint2 Rand)
{
float2 sf = float2(Rand) * (sqrt(2.) / 0xffff) - sqrt(0.5); // map 0..0xffff to -sqrt(0.5)..sqrt(0.5)
float2 sq = sf*sf;
float root = sqrt(2.*max(sq.x, sq.y) - min(sq.x, sq.y));
if (sq.x > sq.y)
{
sf.x = sf.x > 0 ? root : -root;
}
else
{
sf.y = sf.y > 0 ? root : -root;
}
return sf;
}
參考
Square/Disc mappings
http://marc-b-reynolds.github.io/math/2017/01/08/SquareDisc.html
Square/Disc mappings
January 8th, 2017
This note is for some auxilary information on some maps between the disc and square. A fair amount of the math is convered by Fong1 2 and the survey by Lambers3. Specifically this contains:
- A pair of animated plots between the various maps
- Show the Jacobian matrices and determinates and visualize with interactive plots
- Reduced complexity radial stretch in both directions
- Reduced complexity concentric disc to square
- Show a (perhaps) new radial approximate area preserving map
I will use the convention that (x, y)
is a coordinate on the square and (u, v)
a coordinate on the disc.
Define a signum like function:
sgn(x)={1−1x≥0x<0
Area distortion
This shadertoy is intended to visualize how shapes and angles are transformed under the various maps. Clicking through to the site allows interactivity.
Point set
As an alternate visualization this plot take uniform points on the square and shows how the various square to disc mappings transform the point set.
−1−0.500.51−1−0.500.51
Radial Stretching
Simplest square to disc map is to simply stretch (scale factor of L1
over L2
norm):
(u, v)=max(|x|, |y|)x2+y2−−−−−−√(x, y)
Fong followed by Lambers provide an alternate formulation intended to be computationally friendly. Note that we can reformulate the scale factor s
applied to the input coordinate as the following pseudo-code:
float x2 = x*x;
float y2 = y*y;
float m = x2 >= y2 ? x : y;
float s = abs(m)*inversesqrt(x2+y2+epsilon);
// return s*(x,y)
where epsilon
is some sufficiently small constant4. Formulated in this manner only needs a single select/cmov like operation and no branching to handle all cases including degenerate.
A method to measure area distortion (growth/shrinkage) at infinitesimal around a point is to compute the Jacobian determinate5 of the function.
First we need the Jacobian, for x2>y2
:
(x2+y2)−32[x3+2y2xy3−x2yx3]
and for x2≤y2
(x2+y2)−32[y3−xy2x3y3+2x2y]
The determinates of these are respectively:
x2x2+y2y2x2+y2
Combining these into a single expression gives:
max(x2,y2)x2+y2
The area of the square is 4 and that of the disc is π
, so the value where no area distortion occurs is π4≈0.785398
. Larger/smaller values are the local stretch/compressing respectively. The plot of the area distortion (determinate):
−1−0.500.5−1−0.500.500.20.40.60.81
The disc to square transform is simply inverting the scale factor:
(x, y)=u2+v2−−−−−−√max(|u|, |v|)(u, v)
and like above all cases can be handled with no branches and a single select.
Concentric
Peter Shirley and Kenneth Chiu in 19976 derived an area preserving map between square to disc.
The square to disc map as given in Shirley’s blog post7 (modifications noted by Dave Cline):
(u, v)=⎧⎩⎨⎪⎪(xcos(π4yx),xsin(π4yx))(ycos(π2−π4xy),ysin(π2−π4xy))|x|≥|y||x|<|y|
Phase shift and pulling out signs:
(u, v)=⎧⎩⎨⎪⎪(xcos(π4yx),xsin(π4yx))(ysin(π4xy),ycos(π4xy))|x|≥|y||x|<|y|
The more interesting form will be architecture/code specific. The second form trig inputs are on ±π4
, so the cosine terms are always positive and removable without sign fixup and a narrow enough range for light weight approximation.
The Jacobian for |x|≥|y|
:
[cos(πy4x)+πy4xsin(πy4x)sin(πy4x)−πy4xcos(πy4x)−π4sin(πy4x)π4cos(πy4x)]
The Jacobian for |x|<|y|
:
⎡⎣⎢π4cos(π4y(x+y))−π4sin(π4y(x+y))sin(π4y(x+y))−πx4ycos(π4y(x+y))cos(π4y(x+y))+πx4ysin(π4y(x+y))⎤⎦⎥
The determinant of both are constant (π4)
as expected for an area preserving map.
The disc to square map:
(x, y)=⎧⎩⎨⎪⎪sgn(u)u2+v2−−−−−−√(1,4πatan(vu))sgn(v)u2+v2−−−−−−√(4πatan(uv),1)u2>v2u2≤v2
Where atan
is single parameter. This is degenerate when v is approaching zero. This could be corrected by a select or since atan
is an odd function one possible rewrite choice is:
(x, y)=⎧⎩⎨⎪⎪⎪⎪u2+v2−−−−−−√(sgn(u),4πatan(v|u|))u2+v2−−−−−−√(4πatan(u|v|+ϵ),sgn(v))u2>v2u2≤v2
In either case the range of atan
in on ±1 so can be lightweight (for it) approximated. (ϵ
is some small constant4 as above)
F.G. Squircle
Fong1 derivated mappings based on Fernandez-Guasti’s squircle.
(u, v)=x2+y2−x2y2−−−−−−−−−−−−√x2+y2−−−−−−√(x,y)
The Jacobian:
1(x2+y2)3/2y2−x2(y2−1)−−−−−−−−−−−−−√[y4−(y2−1)x4−2y2(y2−1)x2−xy5−x5yy4+(1−2y2)x4−y2(y2−2)x2]
The Jacobian determinate:
1−2x2y2x2+y2
the plot of the area distortion (determinate):
−1−0.500.5−1−0.500.500.20.40.60.81
The disc to square map:
n=u2+v2s=sgn(uv)2–√n−n(n−4u2v2)−−−−−−−−−−−√−−−−−−−−−−−−−−−−√
(x, y)=s(1u, 1v)
Elliptical
The square to disc mapping was derived by Nowell in a blog post8:
(u, v)=(x1−y22−−−−−−√, y1−x22−−−−−−√)
The Jacobian:
⎛⎝⎜⎜1−y22−−−−−√−xy4−2x2√−xy4−2y2√1−x22−−−−−√⎞⎠⎟⎟
The Jacobian determinate:
2−(x2+y2)2−x2−−−−−√2−y2−−−−−√
the plot of the area distortion (determinate):
−1−0.500.5−1−0.500.500.20.40.60.81
Fong1 provides two derivations of the inverse map, one of which:
(x, y)=12(2+t+22–√u−−−−−−−−−−−√−2+t−22–√u−−−−−−−−−−−√,2−t+22–√v−−−−−−−−−−√−2−t−22–√v−−−−−−−−−−√)
where:
t=u2−v2
Approximate equal area
I haven’t been able to find this method referenced anywhere and a twitter query came up negative as well.
The basic idea here is to form a low complexity middle ground between radial stretch and concentric. See the point set animation above. As such it is only of interest if computationally cheaper than concentric and higher performance is a priority.
The disc to square map:
(x, y)=⎧⎩⎨⎪⎪⎪⎪(sgn(u)u2+v2−−−−−−√, 2–√v)(2–√u, sgn(v)u2+v2−−−−−−√)u2>v2u2≤v2
The square to disc map:
(u, v)=⎧⎩⎨⎪⎪⎪⎪⎪⎪(x1−y22x2−−−−−−√, y2√)(x2√, y1−x22y2−−−−−−√)x2>y2x2≤y2
Which is degenerate when y≤x
and y
approaching zero. Since the division term is positive we can simply add a bias as before. More likely to be interesting is we can multiply through by the denominator and get:
(u, v)=⎧⎩⎨⎪⎪12√(sgn(x)2x2−y2−−−−−−−√, y)12√(x, sgn(y)2y2−x2−−−−−−−√)x2>y2x2≤y2
or leave the half inside:
(u, v)=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪(sgn(x)x2−y22−−−−−−√, y2√)(x2√, sgn(y)y2−x22−−−−−−√)x2>y2x2≤y2
The Jacobian for |x|≥|y|
:
⎛⎝⎜⎜11−y22x2√0−yx4−2y2x2√12√⎞⎠⎟⎟
and the determinate:
12−y2x2−−−−−√
The Jacobian for |x|<|y|
:
⎛⎝⎜⎜12√−x4−2x2y2√y011−x22y2√⎞⎠⎟⎟
and the determinate:
12−x2y2−−−−−√
and the heat map of area distortion. Recall equal area value is π4≈0.785398
−1−0.500.5−1−0.500.500.20.40.60.81
Approximate area preserving (variant 2)
In another post9 we have a volume preserving map between cylinder and sphere which breaks the space into two parts: point inside and outside of the embedded conic. We can take a planar slice through the cylinder’s axis to form a square to disc map. Reducing the 3D equation to the slice gives:
(u, v)=⎧⎩⎨⎪⎪⎪⎪⎪⎪(x1−49y2x2−−−−−−−√, 23y)(x23−19x2y2−−−−−−−√, y−13x2y)x2≥y2x2<y2
We can rewrite the x2≥y2
case to eliminate division (which can be zero) as:
t=23y(sgn(x)x2−t2−−−−−−√, t)
and re-express the x2<y2
case as:
t=x3y(x23−t2−−−−−−√, y−xt)
The Jacobian for for the cases are:
⎛⎝⎜⎜39−4y2x2√0−4y3x9−4y2x2√23⎞⎠⎟⎟
⎛⎝⎜⎜⎜6y2−2x236−y2x2y2√−2x3yx33y36−x2y2√1+x23y2⎞⎠⎟⎟⎟
with the determinates:
29−4y2x2−−−−−−√26−x2y2−−−−−√
The area distortion:
−1−0.500.5−1−0.500.500.20.40.60.81
For the disc to square map, start with a common term:
t=u2+v2−−−−−−√
then the map is:
(x, y)=⎧⎩⎨⎪⎪⎪⎪⎪⎪(u3tt+|v|−−−−√, sgn(v)t)(sgn(u)t, 32v)54v2≥u254v2<u2
References and Footnotes
-
“Analytical Methods for Squaring the Disc”, Chamberlain Fong, 2015. (PDF) ↩ ↩2 ↩3
-
“Analytical Methods for Squaring the Disc (presentation slides)”, Chamberlain Fong, 2014. (link) ↩
-
“Mappings between Sphere, Disc, and Square”, Martin Lambers, 2016. (link) ↩
-
A reasonable(ish) epsilon for the funcs here under generic usage is the smallest normal fp value (single: 2−126