|
16 | 16 | import xarray as xr
|
17 | 17 | import matplotlib
|
18 | 18 | import quaternion as quat
|
| 19 | +import math |
19 | 20 |
|
20 | 21 | from utils import is_notebook, download, pushd
|
21 | 22 |
|
@@ -1004,6 +1005,100 @@ def vec_rotate(v, theta, axis):
|
1004 | 1005 | return v_prime.imag
|
1005 | 1006 |
|
1006 | 1007 |
|
| 1008 | +def magnitude(vec): |
| 1009 | + return np.linalg.norm(vec) |
| 1010 | + |
| 1011 | +def normalise(vec): |
| 1012 | + norm = np.linalg.norm(vec) |
| 1013 | + if norm == 0: |
| 1014 | + vn = vec |
| 1015 | + else: |
| 1016 | + vn = vec / norm |
| 1017 | + return vn |
| 1018 | + |
| 1019 | +def vector_align(v1, v2, up=[0,1,0], lvformat=True): |
| 1020 | + """ |
| 1021 | + Get a rotation quaterion to align vectors v1 with v2 |
| 1022 | +
|
| 1023 | + Parameters |
| 1024 | + ---------- |
| 1025 | + v1 : list/numpy.ndarray |
| 1026 | + First 3 component vector |
| 1027 | + v2 : list/numpy.ndarray |
| 1028 | + Second 3 component vector to align the first to |
| 1029 | +
|
| 1030 | + Returns |
| 1031 | + ------- |
| 1032 | + list: quaternion to rotate v1 to v2 (in lavavu format) |
| 1033 | + """ |
| 1034 | + |
| 1035 | + #Check for parallel or opposite |
| 1036 | + v1 = normalise(np.array(v1)) |
| 1037 | + v2 = normalise(np.array(v2)) |
| 1038 | + epsilon = np.finfo(np.float32).eps |
| 1039 | + one_minus_eps = 1.0 - epsilon |
| 1040 | + if np.dot(v1, v2) > one_minus_eps: # 1.0 |
| 1041 | + #No rotation |
| 1042 | + return [0,0,0,1] |
| 1043 | + elif np.dot(v1, v2) < -one_minus_eps: # -1.0 |
| 1044 | + #180 rotation about Y |
| 1045 | + return [0,1,0,1] |
| 1046 | + xyz = np.cross(v1, v2) |
| 1047 | + l1 = np.linalg.norm(v1) |
| 1048 | + l2 = np.linalg.norm(v2) |
| 1049 | + w = math.sqrt((l1*l1) * (l2*l2)) + np.dot(v1, v2) |
| 1050 | + qr = quat.quaternion(w, xyz[0], xyz[1], xyz[2]) |
| 1051 | + qr = qr.normalized() |
| 1052 | + #Return in LavaVu quaternion format |
| 1053 | + if lvformat: |
| 1054 | + return [qr.x, qr.y, qr.z, qr.w] |
| 1055 | + else: |
| 1056 | + return qr |
| 1057 | + |
| 1058 | +def lookat(lv, pos, lookat=None, up=None): |
| 1059 | + """ |
| 1060 | + Set the camera with a position coord and lookat coord |
| 1061 | +
|
| 1062 | + Parameters |
| 1063 | + ---------- |
| 1064 | + lv : lavavu.Viewer |
| 1065 | + The viewer object |
| 1066 | + pos : list/numpy.ndarray |
| 1067 | + Camera position in world coords |
| 1068 | + lookat : list/numpy.ndarray |
| 1069 | + Look at position in world coords, defaults to model origin |
| 1070 | + up : list/numpy.ndarray |
| 1071 | + Up vector, defaults to Y axis [0,1,0] |
| 1072 | + """ |
| 1073 | + |
| 1074 | + #Use the origin from viewer if no target provided |
| 1075 | + if lookat is None: |
| 1076 | + lookat = lv['focus'] |
| 1077 | + else: |
| 1078 | + lv['focus'] = lookat |
| 1079 | + |
| 1080 | + #Default to Y-axis up vector |
| 1081 | + if up is None: |
| 1082 | + up = np.array([0,1,0]) |
| 1083 | + |
| 1084 | + #Calculate the lookat rotation matrix |
| 1085 | + heading = np.array(pos) - np.array(lookat) # inverse line of sight |
| 1086 | + zd = normalise((heading)) |
| 1087 | + xd = normalise(np.cross(up, zd)) |
| 1088 | + yd = normalise(np.cross(zd, xd)) |
| 1089 | + q = quat.from_rotation_matrix(np.array([xd, yd, zd])) |
| 1090 | + q = q.normalized() |
| 1091 | + qr = [q.x, q.y, q.z, q.w] |
| 1092 | + |
| 1093 | + #Apply the rotation |
| 1094 | + lv.rotation(*qr) |
| 1095 | + |
| 1096 | + #Calc translation, just zoom back by heading vector length in Z |
| 1097 | + tr = [0, 0, -magnitude(np.array(pos) - np.array(lookat))] |
| 1098 | + |
| 1099 | + #Apply translation |
| 1100 | + lv.translation(tr) |
| 1101 | + |
1007 | 1102 | def sun_light(time=None, now=False, local=True, tz=None, hour=None, minute=None, xyz=[0.4, 0.4, 1.0]):
|
1008 | 1103 | """
|
1009 | 1104 | Setup a sun light for earth model illumination
|
|
0 commit comments