Skip to content

Instantly share code, notes, and snippets.

@k3a
Last active May 30, 2022 13:26
Show Gist options
  • Star 2 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save k3a/2903719bb42b48c9198d20c2d6f73ac1 to your computer and use it in GitHub Desktop.
Save k3a/2903719bb42b48c9198d20c2d6f73ac1 to your computer and use it in GitHub Desktop.
#!/usr/bin/env python3
#
# Test script to convert between xtilt/ytilt (Qt, Windows) and azimuth/altitude (iOS)
# Free to use for any purpose
import math
import secrets
# [azimuthRad, altitudeRad] => [tiltXrad, tiltYrad]
# azimuth in [0, 2*pi] range, altitude in [0, pi/2]
# X = sin(azimuth) * cos(altitude)
# Y = cos(azimuth) * cos(altitude)
# Z = sin(altitude)
# X Tilt = arctan(X / Z)
# Y Tilt = arctan(Y / Z)
# Source: Qt https://code.woboq.org/qt5/qtbase/src/plugins/platforms/windows/qwindowstabletsupport.cpp.html#527
def spherical2tilt(inVec):
azimuthRad, altitudeRad = inVec[0], inVec[1]
if altitudeRad == math.pi/2.0:
# perpendicular to the pad
return [0, 0]
elif altitudeRad == 0:
# when pen is laying on the pad it is impossible to precisely encode but at least approximate for 4 cases
if azimuthRad > 7 * math.pi/4 or azimuthRad <= math.pi/4:
# for azimuthRad == 0, the pen is on the positive Y axis
return [0, math.pi/2]
elif azimuthRad > math.pi/4 and azimuthRad <= 3 * math.pi/4:
# for azimuthRad == math.pi/2 the pen is on the positive X axis
return [math.pi/2, 0]
elif azimuthRad > 3 * math.pi/4 and azimuthRad <= 5 * math.pi/4:
# for azimuthRad == math.pi, the pen is on the negative Y axis
return [0, -math.pi/2]
elif azimuthRad > 5 * math.pi/4 and azimuthRad <= 7 * math.pi/4:
# for azimuthRad == math.pi + math.pi/2 pen on negative X axis
return [-math.pi/2, 0]
tanAlt = math.tan(altitudeRad) # tan(x) = sin(x)/cos(x)
tiltXrad = math.atan(math.sin(azimuthRad) / tanAlt)
tiltYrad = math.atan(math.cos(azimuthRad) / tanAlt)
return [tiltXrad, tiltYrad]
# [tiltXrad, tiltYrad] => [azimuthRad, altitudeRad]
# tiltX and tiltY should be in the range [0, pi/2]
def tilt2spherical(inVec):
tiltXrad, tiltYrad = inVec[0], inVec[1]
if tiltXrad == 0 and tiltYrad == 0:
# pen perpendicular to the pad
return [0, math.pi/2]
# X and Y of a vector created by the intersection of tilt planes
# first normal vectors of tiltX and tiltY planes are defined
# from that cross product is done to find this vector perpendicular to both plane's normal vector
# in this unit vector Z is ignored to get x y coords projected on the pad
y = math.cos(tiltXrad) * math.sin(tiltYrad)
x = -math.sin(tiltXrad) * -math.cos(tiltYrad)
z = -math.cos(tiltXrad)*-math.cos(tiltYrad)
# compute angle of the projected 2D vector to get azimuth in the proper direction
azimuthRad = -math.atan2(y, x) + math.pi/2
if azimuthRad < 0:
# make always positive in range from 0 to 2*pi
azimuthRad += 2*math.pi
vecLenOn2DPad = math.sqrt(x*x+y*y)
altitudeRad = math.atan(z / vecLenOn2DPad)
# other possible, simpler way to get altitudeRad which is not 100% correct:
# deviation: max(7.96°) / avg(2.00°) / median(0.91°)
# not derived from anything, just two 2D situations combined by a multiplication
# altitudeRad = math.pi/2-math.acos(math.cos(tiltXrad) * math.cos(tiltYrad))
return [azimuthRad, altitudeRad]
def randfloat():
return secrets.randbelow(2**32)/2**32
azMax = 10
altMax = 10
n = 0
totalBugs = 0
for az in range(azMax+1):
for alt in range(azMax+1):
azimuthRad = (2.0*math.pi) / azMax * az
altitudeRad = (math.pi/2.0) / altMax * alt
azimAltVec = [azimuthRad, altitudeRad]
tiltVec = spherical2tilt(azimAltVec)
azimAltVec2 = tilt2spherical(tiltVec)
# correction check
res = ""
tolerance = math.pi/4 if altitudeRad == 0 else 0.1
# normalize azimuth before making difference comparison
# by snapping azimuth to the 0 azimuth for azimuths between [2*pi-tolerance, 2*pi]
azim = azimAltVec[0]
azim2 = azimAltVec2[0]
if azim >= 2*math.pi - tolerance:
azim = 0
tolerance = 0.000001
if azim2 >= 2*math.pi - tolerance:
azim2 = 0
tolerance = 0.000001
diffAz = azim2-azim
diffAlt = azimAltVec2[1]-azimAltVec[1]
if altitudeRad < math.pi/2:
if abs(diffAz) > tolerance:
res = f"!! azimuth over tolerance {diffAz} at {n} !!"
totalBugs += 1
if abs(diffAlt) > tolerance:
res = f"!! altitude over tolerance {diffAlt} at {n} !!"
totalBugs += 1
print("az,alt: %6.2f, %6.2f => tx,tY: %6.2f, %6.2f => az,alt: %6.2f, %6.2f (diff %6.2f, %6.2f rad) %s"
% (azimAltVec[0], azimAltVec[1],
tiltVec[0], tiltVec[1],
azimAltVec2[0], azimAltVec2[1],
diffAz, diffAlt,
res
)
)
n += 1
if totalBugs > 0:
print(f"TOTAL {totalBugs} BUGS")

spherical2tilt

  • it is copy & paste from Qt code mentioned there
  • assumes X Y Z position of the pen end ("pen rubber") like this: X = sin(azimuth) * cos(altitude) Y = cos(azimuth) * cos(altitude) Z = sin(altitude)
  • from the X, Y, Z formulas there it is clear that: X is 1 when azimuth is pi/2 rad (90 deg) and altitude is 0 Y is 1 when azimuth is 0 and altitude is 0 Z is 1 when altitude is pi/2 (90 deg) regardless of azimuth (pen being perpendicular to the plane) thus when altitude is 0 (pen laying on the pad) and azimuth slowly rising from 0 into pi/2 you can see pen end leaving positive Y axis and going towards positive X axis still laying (altitude is 0) and azimuth going from pi/2 into pi, you can see pen end leaving positive X axis, going towards negative Y axis still laying with azimuth going from pi into 3pi/2 you can see pen end leaving negative Y axis, going towards negative X still laying when azimuth going from 3pi/2 into 2*pi (=0) you can see pen end levaing legative X, going towards to the positive Y axis again
  • position Z of the pen end is easy sin(altitude) because it is 1 when altitude = pi/2 (90 deg) and pen is perpendicular to the pad and Z=0 when pen parallel with the pad
  • then Qt gues say this:
    • XTilt = arctan(X / Z) arctan is reverse of tan (you ask for which angle the tan function would return the number you pass into arctan)
    • XTilt is actually an angle between two planes where intersection of those planes is equal to Y axis (see https://www.w3.org/Submission/pointer-events/tiltX_600px.png)
    • tan(x) in 2D for a triangle is defined as = opposite / adjacent
    • When you add an imaginary line on top corners on the planes in https://www.w3.org/Submission/pointer-events/tiltX_600px.png, you see such triangle formed and for that triangle tan(tiltX) = "pen end position in the X axis" (opposite) / "pen end position in Z axis" (adjacent) . Where "pen end in X axis" is actually the same length as the length of the added imaginary line.
    • So to get the angle XTilt when we know tan, we do atan, thus XTilt = arctan(X / Z)
    • It works like that but Qt guys further optimized the code this way:
    • Xtilt = sin(radAzim) / tanAlt why is it equal to the previous formula for Xtilt? see:
    • fill in X and Z formulas into the XTilt = arctan(X / Z) = arctan( (sin(azimuth) * cos(altitude)) / sin(altitude) )
    • another definition of tan is: tan(x) = sin(x)/cos(x) so compute tanAlt = sin(altitude) / cos(altitude)
    • as you can see tanAlt = sin(altitude) / cos(altitude) which is opposite of what we have in arctan argument of XTilt formula: (sin(azimuth) * cos(altitude)) / sin(altitude
    • thus we divide XTilt = sin(radAzim) / tanAlt
    • Ytilt works the same way, just for the other axis

tilt2spherical

  • this is my code (I never needed this conversion, I made it just because I was curious whether apple azimuth/atitude encoding has the same or better encoding precision as the older tiltX/tiltY method) and maybe to learn something new. Later I discovered some possible inaccuracies but I never had time to devote so much time again to fix it up. I don't have time now either but I decided it must be done finally.
  • Result: there were wrong signs in some cases and it missed couple of special cases. Read on to understand it fully:

azimuthRad:

  • I compute azimuthRad using atan2
  • now atan2(y, x) is a function which returns the angle between X axis and a vector with coordinates X,Y with the angle rising counter-clockwise see https://upload.wikimedia.org/wikipedia/commons/thumb/7/72/Sinus_und_Kosinus_am_Einheitskreis_1.svg/375px-Sinus_und_Kosinus_am_Einheitskreis_1.svg.png works like this:
    • for I. quadrant (X>0, Y>0) it returns angles from 0 to pi/2
    • for II. quadrant (X<0, Y>0) it returns angles from pi/2 to pi
    • for III. quadrant (X<0, Y<0) it returns angles from -pi to -pi/2
    • for IV. quandrant (X>0, Y<0) it returns angles -pi/2 to 0 The vector here is pen "rubber" position on 2D pad (ground) plane. But first I need to get 3D position of the pen "rubber" by doing intersection of tilt planes:

tiltX plane normal vector when tiltX=0: X=-1 Y=0 Z=0 with tiltX rising positively: X must go towards 0 so X=-cos(tiltX), Z must rise positively so Z=sin(tiltX) now we know normal vector of tiltX plane is: X=-cos(tiltX) Y=0 Z=sin(tiltX)

tiltY plane normal vector when tiltY=0: X=0 Y=-1 Z=0 with tiltY rising positively: Y must go towards 0 so Y=-cos(tiltY) X=0 Z must rise so Z=sin(tiltY) not we know normal vector of tiltY plane: X=0 Y=-cos(tiltY) Z=sin(tiltY)

Now we just need to compute intersection of two planes to find the intersecting line vector. This line is perpendicular to normal vectors of those planes so we need to do cross product: vec = tiltYnormal (=a) x tiltYnormal (=b) Cross product is defined as: vec.x = a.yb.z - a.zb.y vec.y = a.zb.x - a.xb.z vec.z = a.xb.y - a.yb.x thus: vec.x = 0sin(tiltY) - sin(tiltX)-cos(tiltY) vec.y = sin(tiltX)0 - -cos(tiltX)sin(tiltY) vec.z = -cos(tiltX)-cos(tiltY) - 00 remove zeroes: vec.x = -sin(tiltX)*-cos(tiltY) vec.y = cos(tiltX)sin(tiltY) vec.z = -cos(tiltX)-cos(tiltY)

Now ve heve vector representing direction from the 0,0,0 origin to the pen end (rubber). We need to project it to the ground to get 2D representation so we ignore Z and use X and Y only for atan2(y,x). Result of atan2 needs to be negated, moved by pi/2 and normalized into 0 - 2pi range (this works for all 4 quadrants). Finally that gives us "azimuth" angle.

The last thing is to compute altitude angle. Vector vec which we got as a result of cross product may not always be of unit length (if it was the case we could optimize). Altitude angle is angle between the 3D position of the pen end [vec.x, vec.y, vec.z] and its projection to XY 2d plane [X, Y, 0] which can be computed from atan: altitudeAngle = atan( vec.z / (length-of-2d-projection-of-vec-on-the-pad) )

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment