1 /************************************/
2 /* ArcBall.c (c) Ken Shoemake, 1993 */
3 /* Modified by Tom Bech, 1996 */
4 /************************************/
5
6 #include "config.h"
7
8 #include <libgimp/gimp.h>
9
10 #include "arcball.h"
11
12 /* Global variables */
13 /* ================ */
14
15 Quat qOne = { 0, 0, 0, 1 };
16
17 static HVect center;
18 static double radius;
19 static Quat qNow, qDown, qDrag;
20 static HVect vNow, vDown, vFrom, vTo, vrFrom, vrTo;
21 static HMatrix mNow, mDown;
22 static unsigned int showResult, dragging;
23 static ConstraintSet sets[NSets];
24 static int setSizes[NSets];
25 static AxisSet axisSet;
26 static int axisIndex;
27
28 static HMatrix mId =
29 {
30 { 1, 0, 0, 0 },
31 { 0, 1, 0, 0 },
32 { 0, 0, 1, 0 },
33 { 0, 0, 0, 1 }
34 };
35
36 static double otherAxis[][4] =
37 {
38 {-0.48, 0.80, 0.36, 1}
39 };
40
41 /* Internal methods */
42 /* ================ */
43
44 static void Qt_ToMatrix(Quat q,HMatrix out);
45 static Quat Qt_Conj(Quat q);
46 static Quat Qt_Mul(Quat qL, Quat qR);
47 static Quat Qt_FromBallPoints(HVect from, HVect to);
48 static void Qt_ToBallPoints(Quat q, HVect *arcFrom, HVect *arcTo);
49
50 static HVect V3_(double x, double y, double z);
51 static double V3_Norm(HVect v);
52 static HVect V3_Unit(HVect v);
53 static HVect V3_Scale(HVect v, double s);
54 static HVect V3_Negate(HVect v);
55 /*
56 static HVect V3_Add(HVect v1, HVect v2);
57 */
58 static HVect V3_Sub(HVect v1, HVect v2);
59 static double V3_Dot(HVect v1, HVect v2);
60 /*
61 static HVect V3_Cross(HVect v1, HVect v2);
62 static HVect V3_Bisect(HVect v0, HVect v1);
63 */
64
65 static HVect MouseOnSphere(HVect mouse, HVect ballCenter, double ballRadius);
66 static HVect ConstrainToAxis(HVect loose, HVect axis);
67 static int NearestConstraintAxis(HVect loose, HVect *axes, int nAxes);
68
69 /* Establish reasonable initial values for controller. */
70 /* =================================================== */
71
72 void
ArcBall_Init(void)73 ArcBall_Init (void)
74 {
75 int i;
76
77 center = qOne;
78 radius = 1.0;
79 vDown = vNow = qOne;
80 qDown = qNow = qOne;
81 for (i=15; i>=0; i--)
82 ((double *)mNow)[i] = ((double *)mDown)[i] = ((double *)mId)[i];
83
84 showResult = dragging = FALSE;
85 axisSet = NoAxes;
86 sets[CameraAxes] = mId[X];
87 setSizes[CameraAxes] = 3;
88 sets[BodyAxes] = mDown[X];
89 setSizes[BodyAxes] = 3;
90 sets[OtherAxes] = otherAxis[X];
91 setSizes[OtherAxes] = 1;
92 }
93
94 /* Set the center and size of the controller. */
95 /* ========================================== */
96
97 void
ArcBall_Place(HVect Center,double Radius)98 ArcBall_Place (HVect Center,
99 double Radius)
100 {
101 center = Center;
102 radius = Radius;
103 }
104
105 /* Incorporate new mouse position. */
106 /* =============================== */
107
108 void
ArcBall_Mouse(HVect v_Now)109 ArcBall_Mouse (HVect v_Now)
110 {
111 vNow = v_Now;
112 }
113
114 /* Choose a constraint set, or none. */
115 /* ================================= */
116
117 void
ArcBall_UseSet(AxisSet axis_Set)118 ArcBall_UseSet (AxisSet axis_Set)
119 {
120 if (!dragging) axisSet = axis_Set;
121 }
122
123 /* Using vDown, vNow, dragging, and axisSet, compute rotation etc. */
124 /* =============================================================== */
125
126 void
ArcBall_Update(void)127 ArcBall_Update (void)
128 {
129 int setSize = setSizes[axisSet];
130 HVect *set = (HVect *)(sets[axisSet]);
131
132 vFrom = MouseOnSphere(vDown, center, radius);
133 vTo = MouseOnSphere(vNow, center, radius);
134 if (dragging)
135 {
136 if (axisSet!=NoAxes)
137 {
138 vFrom = ConstrainToAxis(vFrom, set[axisIndex]);
139 vTo = ConstrainToAxis(vTo, set[axisIndex]);
140 }
141 qDrag = Qt_FromBallPoints(vFrom, vTo);
142 qNow = Qt_Mul(qDrag, qDown);
143 }
144 else
145 {
146 if (axisSet!=NoAxes) axisIndex = NearestConstraintAxis(vTo, set, setSize);
147 }
148 Qt_ToBallPoints(qDown, &vrFrom, &vrTo);
149 Qt_ToMatrix(Qt_Conj(qNow), mNow); /* Gives transpose for GL. */
150 }
151
152 /* Return rotation matrix defined by controller use. */
153 /* ================================================= */
154
155 void
ArcBall_Value(HMatrix m_Now)156 ArcBall_Value (HMatrix m_Now)
157 {
158 ArcBall_CopyMat (mNow, m_Now);
159 }
160
161 /* Extract rotation angles from matrix */
162 /* =================================== */
163
164 void
ArcBall_Values(double * alpha,double * beta,double * gamma)165 ArcBall_Values (double *alpha,
166 double *beta,
167 double *gamma)
168 {
169 if ((*beta=asin(-mNow[0][2]))!=0.0)
170 {
171 *gamma=atan2(mNow[1][2],mNow[2][2]);
172 *alpha=atan2(mNow[0][1],mNow[0][0]);
173 }
174 else
175 {
176 *gamma=atan2(mNow[1][0],mNow[1][1]);
177 *alpha=0.0;
178 }
179 }
180
181 /* Begin drag sequence. */
182 /* ==================== */
183
184 void
ArcBall_BeginDrag(void)185 ArcBall_BeginDrag (void)
186 {
187 dragging = TRUE;
188 vDown = vNow;
189 }
190
191 /* Stop drag sequence. */
192 /* =================== */
193
194 void
ArcBall_EndDrag(void)195 ArcBall_EndDrag (void)
196 {
197 dragging = FALSE;
198 qDown = qNow;
199
200 ArcBall_CopyMat (mNow, mDown);
201 }
202
203 /*===================*/
204 /***** BallAux.c *****/
205 /*===================*/
206
207 /* Return quaternion product qL * qR. Note: order is important! */
208 /* To combine rotations, use the product Mul(qSecond, qFirst), */
209 /* which gives the effect of rotating by qFirst then qSecond. */
210 /* ============================================================= */
211
212 static Quat
Qt_Mul(Quat qL,Quat qR)213 Qt_Mul (Quat qL,
214 Quat qR)
215 {
216 Quat qq;
217 qq.w = qL.w*qR.w - qL.x*qR.x - qL.y*qR.y - qL.z*qR.z;
218 qq.x = qL.w*qR.x + qL.x*qR.w + qL.y*qR.z - qL.z*qR.y;
219 qq.y = qL.w*qR.y + qL.y*qR.w + qL.z*qR.x - qL.x*qR.z;
220 qq.z = qL.w*qR.z + qL.z*qR.w + qL.x*qR.y - qL.y*qR.x;
221 return (qq);
222 }
223
224 /* Construct rotation matrix from (possibly non-unit) quaternion. */
225 /* Assumes matrix is used to multiply column vector on the left: */
226 /* vnew = mat vold. Works correctly for right-handed coordinate */
227 /* system and right-handed rotations. */
228 /* ============================================================== */
229
230 static void
Qt_ToMatrix(Quat q,HMatrix out)231 Qt_ToMatrix (Quat q,
232 HMatrix out)
233 {
234 double Nq = q.x*q.x + q.y*q.y + q.z*q.z + q.w*q.w;
235 double s = (Nq > 0.0) ? (2.0 / Nq) : 0.0;
236 double xs = q.x*s, ys = q.y*s, zs = q.z*s;
237 double wx = q.w*xs, wy = q.w*ys, wz = q.w*zs;
238 double xx = q.x*xs, xy = q.x*ys, xz = q.x*zs;
239 double yy = q.y*ys, yz = q.y*zs, zz = q.z*zs;
240 out[X][X] = 1.0 - (yy + zz); out[Y][X] = xy + wz; out[Z][X] = xz - wy;
241 out[X][Y] = xy - wz; out[Y][Y] = 1.0 - (xx + zz); out[Z][Y] = yz + wx;
242 out[X][Z] = xz + wy; out[Y][Z] = yz - wx; out[Z][Z] = 1.0 - (xx + yy);
243 out[X][W] = out[Y][W] = out[Z][W] = out[W][X] = out[W][Y] = out[W][Z] = 0.0;
244 out[W][W] = 1.0;
245 }
246
247 /* Return conjugate of quaternion. */
248 /* =============================== */
249
250 static Quat
Qt_Conj(Quat q)251 Qt_Conj (Quat q)
252 {
253 Quat qq;
254 qq.x = -q.x; qq.y = -q.y; qq.z = -q.z; qq.w = q.w;
255 return (qq);
256 }
257
258 /* Return vector formed from components */
259 /* ==================================== */
260
261 static HVect
V3_(double x,double y,double z)262 V3_ (double x,
263 double y,
264 double z)
265 {
266 HVect v;
267 v.x = x; v.y = y; v.z = z; v.w = 0;
268 return (v);
269 }
270
271 /* Return norm of v, defined as sum of squares of components */
272 /* ========================================================= */
273
274 static double
V3_Norm(HVect v)275 V3_Norm (HVect v)
276 {
277 return ( v.x*v.x + v.y*v.y + v.z*v.z );
278 }
279
280 /* Return unit magnitude vector in direction of v */
281 /* ============================================== */
282
283 static HVect
V3_Unit(HVect v)284 V3_Unit (HVect v)
285 {
286 static HVect u = {0, 0, 0, 0};
287 double vlen = sqrt(V3_Norm(v));
288
289 if (vlen != 0.0)
290 {
291 u.x = v.x/vlen;
292 u.y = v.y/vlen;
293 u.z = v.z/vlen;
294 }
295 return (u);
296 }
297
298 /* Return version of v scaled by s */
299 /* =============================== */
300
301 static HVect
V3_Scale(HVect v,double s)302 V3_Scale (HVect v,
303 double s)
304 {
305 HVect u;
306 u.x = s*v.x; u.y = s*v.y; u.z = s*v.z; u.w = v.w;
307 return (u);
308 }
309
310 /* Return negative of v */
311 /* ==================== */
312
313 static HVect
V3_Negate(HVect v)314 V3_Negate (HVect v)
315 {
316 static HVect u = {0, 0, 0, 0};
317 u.x = -v.x; u.y = -v.y; u.z = -v.z;
318 return (u);
319 }
320
321 /* Return sum of v1 and v2 */
322 /* ======================= */
323 /*
324 static HVect
325 V3_Add (HVect v1,
326 HVect v2)
327 {
328 static HVect v = {0, 0, 0, 0};
329 v.x = v1.x+v2.x; v.y = v1.y+v2.y; v.z = v1.z+v2.z;
330 return (v);
331 }
332 */
333 /* Return difference of v1 minus v2 */
334 /* ================================ */
335
336 static HVect
V3_Sub(HVect v1,HVect v2)337 V3_Sub (HVect v1,
338 HVect v2)
339 {
340 static HVect v = {0, 0, 0, 0};
341 v.x = v1.x-v2.x; v.y = v1.y-v2.y; v.z = v1.z-v2.z;
342 return (v);
343 }
344
345 /* Halve arc between unit vectors v0 and v1. */
346 /* ========================================= */
347 /*
348 static HVect
349 V3_Bisect (HVect v0,
350 HVect v1)
351 {
352 HVect v = {0, 0, 0, 0};
353 double Nv;
354
355 v = V3_Add(v0, v1);
356 Nv = V3_Norm(v);
357 if (Nv < 1.0e-5) v = V3_(0, 0, 1);
358 else v = V3_Scale(v, 1/sqrt(Nv));
359 return (v);
360 }
361 */
362
363 /* Return dot product of v1 and v2 */
364 /* =============================== */
365
366 static double
V3_Dot(HVect v1,HVect v2)367 V3_Dot (HVect v1,
368 HVect v2)
369 {
370 return (v1.x*v2.x + v1.y*v2.y + v1.z*v2.z);
371 }
372
373
374 /* Return cross product, v1 x v2 */
375 /* ============================= */
376 /*
377 static HVect
378 V3_Cross (HVect v1,
379 HVect v2)
380 {
381 static HVect v = {0, 0, 0, 0};
382 v.x = v1.y*v2.z-v1.z*v2.y;
383 v.y = v1.z*v2.x-v1.x*v2.z;
384 v.z = v1.x*v2.y-v1.y*v2.x;
385 return (v);
386 }
387 */
388
389 void
ArcBall_CopyMat(HMatrix inm,HMatrix outm)390 ArcBall_CopyMat (HMatrix inm,
391 HMatrix outm)
392 {
393 int x=0,y=0;
394
395 for (x=0;x<4;x++)
396 {
397 for (y=0;y<4;y++)
398 {
399 outm[y][x]=inm[y][x];
400 }
401 }
402 }
403
404 /*=====================================================*/
405 /**** BallMath.c - Essential routines for ArcBall. ****/
406 /*=====================================================*/
407
408 /* Convert window coordinates to sphere coordinates. */
409 /* ================================================= */
410
411 static HVect
MouseOnSphere(HVect mouse,HVect ballCenter,double ballRadius)412 MouseOnSphere (HVect mouse,
413 HVect ballCenter,
414 double ballRadius)
415 {
416 HVect ballMouse;
417 register double mag;
418
419 ballMouse.x = (mouse.x - ballCenter.x) / ballRadius;
420 ballMouse.y = (mouse.y - ballCenter.y) / ballRadius;
421 mag = ballMouse.x*ballMouse.x + ballMouse.y*ballMouse.y;
422 if (mag > 1.0)
423 {
424 register double scale = 1.0/sqrt(mag);
425 ballMouse.x *= scale; ballMouse.y *= scale;
426 ballMouse.z = 0.0;
427 }
428 else ballMouse.z = sqrt(1 - mag);
429 ballMouse.w = 0.0;
430 return (ballMouse);
431 }
432
433 /* Construct a unit quaternion from two points on unit sphere */
434 /* ========================================================== */
435
436 static Quat
Qt_FromBallPoints(HVect from,HVect to)437 Qt_FromBallPoints (HVect from,
438 HVect to)
439 {
440 Quat qu;
441 qu.x = from.y*to.z - from.z*to.y;
442 qu.y = from.z*to.x - from.x*to.z;
443 qu.z = from.x*to.y - from.y*to.x;
444 qu.w = from.x*to.x + from.y*to.y + from.z*to.z;
445 return (qu);
446 }
447
448 /* Convert a unit quaternion to two points on unit sphere */
449 /* ====================================================== */
450
451 static void
Qt_ToBallPoints(Quat q,HVect * arcFrom,HVect * arcTo)452 Qt_ToBallPoints (Quat q,
453 HVect *arcFrom,
454 HVect *arcTo)
455 {
456 double s;
457
458 s = sqrt(q.x*q.x + q.y*q.y);
459 if (s == 0.0) *arcFrom = V3_(0.0, 1.0, 0.0);
460 else *arcFrom = V3_(-q.y/s, q.x/s, 0.0);
461 arcTo->x = q.w*arcFrom->x - q.z*arcFrom->y;
462 arcTo->y = q.w*arcFrom->y + q.z*arcFrom->x;
463 arcTo->z = q.x*arcFrom->y - q.y*arcFrom->x;
464 if (q.w < 0.0) *arcFrom = V3_(-arcFrom->x, -arcFrom->y, 0.0);
465 }
466
467 /* Force sphere point to be perpendicular to axis. */
468 /* =============================================== */
469
470 static HVect
ConstrainToAxis(HVect loose,HVect axis)471 ConstrainToAxis (HVect loose,
472 HVect axis)
473 {
474 HVect onPlane;
475 register double norm;
476
477 onPlane = V3_Sub(loose, V3_Scale(axis, V3_Dot(axis, loose)));
478 norm = V3_Norm(onPlane);
479 if (norm > 0.0)
480 {
481 if (onPlane.z < 0.0) onPlane = V3_Negate(onPlane);
482 return ( V3_Scale(onPlane, 1/sqrt(norm)) );
483 }
484 /* else drop through */
485 /* ================= */
486
487 if (axis.z == 1) onPlane = V3_(1.0, 0.0, 0.0);
488 else onPlane = V3_Unit(V3_(-axis.y, axis.x, 0.0));
489 return (onPlane);
490 }
491
492 /* Find the index of nearest arc of axis set. */
493 /* ========================================== */
494
495 static int
NearestConstraintAxis(HVect loose,HVect * axes,int nAxes)496 NearestConstraintAxis (HVect loose,
497 HVect *axes,
498 int nAxes)
499 {
500 HVect onPlane;
501 register double max, dot;
502 register int i, nearest;
503 max = -1; nearest = 0;
504
505 for (i=0; i<nAxes; i++)
506 {
507 onPlane = ConstrainToAxis(loose, axes[i]);
508 dot = V3_Dot(onPlane, loose);
509 if (dot>max)
510 {
511 max = dot; nearest = i;
512 }
513 }
514 return (nearest);
515 }
516