1 /*
2  * Adapted to D by Antonio Monteiro
3  */
4 /*
5  * (c) Copyright 1993, 1994, Silicon Graphics, Inc.
6  * ALL RIGHTS RESERVED
7  * Permission to use, copy, modify, and distribute this software for
8  * any purpose and without fee is hereby granted, provided that the above
9  * copyright notice appear in all copies and that both the copyright notice
10  * and this permission notice appear in supporting documentation, and that
11  * the name of Silicon Graphics, Inc. not be used in advertising
12  * or publicity pertaining to distribution of the software without specific,
13  * written prior permission.
14  *
15  * THE MATERIAL EMBODIED ON THIS SOFTWARE IS PROVIDED TO YOU "AS-IS"
16  * AND WITHOUT WARRANTY OF ANY KIND, EXPRESS, IMPLIED OR OTHERWISE,
17  * INCLUDING WITHOUT LIMITATION, ANY WARRANTY OF MERCHANTABILITY OR
18  * FITNESS FOR A PARTICULAR PURPOSE.  IN NO EVENT SHALL SILICON
19  * GRAPHICS, INC.  BE LIABLE TO YOU OR ANYONE ELSE FOR ANY DIRECT,
20  * SPECIAL, INCIDENTAL, INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY
21  * KIND, OR ANY DAMAGES WHATSOEVER, INCLUDING WITHOUT LIMITATION,
22  * LOSS OF PROFIT, LOSS OF USE, SAVINGS OR REVENUE, OR THE CLAIMS OF
23  * THIRD PARTIES, WHETHER OR NOT SILICON GRAPHICS, INC.  HAS BEEN
24  * ADVISED OF THE POSSIBILITY OF SUCH LOSS, HOWEVER CAUSED AND ON
25  * ANY THEORY OF LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE
26  * POSSESSION, USE OR PERFORMANCE OF THIS SOFTWARE.
27  *
28  * US Government Users Restricted Rights
29  * Use, duplication, or disclosure by the Government is subject to
30  * restrictions set forth in FAR 52.227.19(c)(2) or subparagraph
31  * (c)(1)(ii) of the Rights in Technical Data and Computer Software
32  * clause at DFARS 252.227-7013 and/or in similar or successor
33  * clauses in the FAR or the DOD or NASA FAR Supplement.
34  * Unpublished-- rights reserved under the copyright laws of the
35  * United States.  Contractor/manufacturer is Silicon Graphics,
36  * Inc., 2011 N.  Shoreline Blvd., Mountain View, CA 94039-7311.
37  *
38  * OpenGL(TM) is a trademark of Silicon Graphics, Inc.
39  */
40 /*
41  * trackball.h
42  * A virtual trackball implementation
43  * Written by Gavin Bell for Silicon Graphics, November 1988.
44  */
45 
46 /*
47  * This size should really be based on the distance from the center of
48  * rotation to the point on the object underneath the mouse.  That
49  * point would then track the mouse as closely as possible.  This is a
50  * simple example, though, so that is left as an Exercise for the
51  * Programmer.
52  */
53 
54 module gl.TrackBall;
55 
56 //private import tango.stdc.math;//std.math;
57 version (Tango)
58 {
59     private import tango.math.Math;
60 }
61 else
62 {
63     private import std.math;
64 }
65 
66 public:
67 const float TRACKBALLSIZE = 0.8;
68 const int RENORMCOUNT = 97;
69  
70 /**
71  * Pass the x and y coordinates of the last and current positions of
72  * the mouse, scaled so they are from (-1.0 ... 1.0).
73  *
74  * The resulting rotation is returned as a quaternion rotation in the
75  * first paramater.
76  */
77 void
78 trackball(float q[4], float p1x, float p1y, float p2x, float p2y)
79 {
80 /*
81  * Ok, simulate a track-ball.  Project the points onto the virtual
82  * trackball, then figure out the axis of rotation, which is the cross
83  * product of P1 P2 and O P1 (O is the center of the ball, 0,0,0)
84  * Note:  This is a deformed trackball-- is a trackball in the center,
85  * but is deformed into a hyperbolic sheet of rotation away from the
86  * center.  This particular function was chosen after trying out
87  * several variations.
88  *
89  * It is assumed that the arguments to this routine are in the range
90  * (-1.0 ... 1.0)
91  */
92     float[3] a; /* Axis of rotation */
93     float phi;  /* how much to rotate about axis */
94     float[3] p1;
95 	float[3] p2;
96 	float[3] d;
97     float t;
98 
99     if (p1x == p2x && p1y == p2y) {
100         /* Zero rotation */
101         vzero(q.ptr);
102         q[3] = 1.0;
103         return;
104     }
105 
106     /*
107      * First, figure out z-coordinates for projection of P1 and P2 to
108      * deformed sphere
109      */
110     vset(p1,p1x,p1y,tb_project_to_sphere(TRACKBALLSIZE,p1x,p1y));
111     vset(p2,p2x,p2y,tb_project_to_sphere(TRACKBALLSIZE,p2x,p2y));
112 
113     /*
114      *  Now, we want the cross product of P1 and P2
115      */
116     vcross(p2.ptr,p1.ptr,a.ptr);
117 
118     /*
119      *  Figure out how much to rotate around that axis.
120      */
121     vsub(p1,p2,d);
122     t = vlength(d) / (2.0*TRACKBALLSIZE);
123 
124     /*
125      * Avoid problems with out-of-control values...
126      */
127     if (t > 1.0) t = 1.0;
128     if (t < -1.0) t = -1.0;
129     phi = 2.0 * asin(t);//std.math.asin(t);
130 
131     axis_to_quat(a,phi,q);
132 }
133 
134 
135 /*
136  * Given two quaternions, add them together to get a third quaternion.
137  * Adding quaternions to get a compound rotation is analagous to adding
138  * translations to get a compound translation.  When incrementally
139  * adding rotations, the first argument here should be the new
140  * rotation, the second and third the total rotation (which will be
141  * over-written with the resulting new total rotation).
142  */
143 void
144 add_quats(float[4] q1, float[4] q2, float[4] dest)
145 {
146 /*
147  * Given two rotations, e1 and e2, expressed as quaternion rotations,
148  * figure out the equivalent single rotation and stuff it into dest.
149  *
150  * This routine also normalizes the result every RENORMCOUNT times it is
151  * called, to keep error from creeping in.
152  *
153  * NOTE: This routine is written so that q1 or q2 may be the same
154  * as dest (or each other).
155  */
156     static int count=0;
157     float[4] t1;
158     float[4] t2;
159     float[4] t3;
160     float tf[4];
161 
162     vcopy(q1.ptr,t1.ptr);
163     vscale(t1.ptr,q2[3]);
164 
165     vcopy(q2.ptr,t2.ptr);
166     vscale(t2.ptr,q1[3]);
167 
168     vcross(q2.ptr,q1.ptr,t3.ptr);
169     vadd(t1.ptr,t2.ptr,tf.ptr);
170     vadd(t3.ptr,tf.ptr,tf.ptr);
171     tf[3] = q1[3] * q2[3] - vdot(q1.ptr,q2.ptr);
172 
173     dest[0] = tf[0];
174     dest[1] = tf[1];
175     dest[2] = tf[2];
176     dest[3] = tf[3];
177 
178     if (++count > RENORMCOUNT) {
179         count = 0;
180         normalize_quat(dest);
181     }
182 }
183 
184 
185 /*
186  * A useful function, builds a rotation matrix in Matrix based on
187  * given quaternion.
188  */
189 void
190 build_rotmatrix(float[4][4] m, float[4] q)
191 {
192     m[0][0] = 1.0 - 2.0 * (q[1] * q[1] + q[2] * q[2]);
193     m[0][1] = 2.0 * (q[0] * q[1] - q[2] * q[3]);
194     m[0][2] = 2.0 * (q[2] * q[0] + q[1] * q[3]);
195     m[0][3] = 0.0;
196 
197     m[1][0] = 2.0 * (q[0] * q[1] + q[2] * q[3]);
198     m[1][1]= 1.0 - 2.0 * (q[2] * q[2] + q[0] * q[0]);
199     m[1][2] = 2.0 * (q[1] * q[2] - q[0] * q[3]);
200     m[1][3] = 0.0;
201 
202     m[2][0] = 2.0 * (q[2] * q[0] - q[1] * q[3]);
203     m[2][1] = 2.0 * (q[1] * q[2] + q[0] * q[3]);
204     m[2][2] = 1.0 - 2.0 * (q[1] * q[1] + q[0] * q[0]);
205     m[2][3] = 0.0;
206 
207     m[3][0] = 0.0;
208     m[3][1] = 0.0;
209     m[3][2] = 0.0;
210     m[3][3] = 1.0;
211 }
212 
213 
214 /*
215  * This function computes a quaternion based on an axis (defined by
216  * the given vector) and an angle about which to rotate.  The angle is
217  * expressed in radians.  The result is put into the third argument.
218  */
219 void
220 axis_to_quat(float[3] a, float phi, float[4] q)
221 {
222     vnormal(a);
223     vcopy(a.ptr,q.ptr);
224     vscale(q.ptr,sin(phi/2.0));
225     q[3] = cos(phi/2.0);
226 }
227 
228 /*
229  * Quaternions always obey:  a^2 + b^2 + c^2 + d^2 = 1.0
230  * If they don't add up to 1.0, dividing by their magnitued will
231  * renormalize them.
232  *
233  * Note: See the following for more information on quaternions:
234  *
235  * - Shoemake, K., Animating rotation with quaternion curves, Computer
236  *   Graphics 19, No 3 (Proc. SIGGRAPH'85), 245-254, 1985.
237  * - Pletinckx, D., Quaternion calculus as a basic tool in computer
238  *   graphics, The Visual Computer 5, 2-13, 1989.
239  */
240 static void
241 normalize_quat(float[4] q)
242 {
243     int i;
244     float mag;
245 
246     mag = (q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + q[3]*q[3]);
247     for (i = 0; i < 4; i++) q[i] /= mag;
248 }
249 
250 /*
251  * Project an x,y pair onto a sphere of radius r OR a hyperbolic sheet
252  * if we are away from the center of the sphere.
253  */
254 static float
255 tb_project_to_sphere(float r, float x, float y)
256 {
257     float d, t, z;
258 
259     d = sqrt(x*x + y*y);
260     if (d < r * 0.70710678118654752440) {    /* Inside sphere */
261         z = sqrt(r*r - d*d);
262     } else {           /* On hyperbola */
263         t = r / 1.41421356237309504880;
264         z = t*t / d;
265     }
266     return z;
267 }
268 
269 void
270 vzero(float* v)
271 {
272     v[0] = 0.0;
273     v[1] = 0.0;
274     v[2] = 0.0;
275 }
276 
277 void
278 vset(float[3] v, float x, float y, float z)
279 {
280     v[0] = x;
281     v[1] = y;
282     v[2] = z;
283 }
284 
285 void
286 vsub(float[3] src1, float[3] src2, float[3] dst)
287 {
288     dst[0] = src1[0] - src2[0];
289     dst[1] = src1[1] - src2[1];
290     dst[2] = src1[2] - src2[2];
291 }
292 
293 void
294 vcopy(float* v1, float* v2)
295 {
296     int i;
297     for (i = 0 ; i < 3 ; i++)
298         v2[i] = v1[i];
299 }
300 
301 void
302 vcross(float* v1, float* v2, float* cross)
303 {
304     float[3] temp;
305 
306     temp[0] = (v1[1] * v2[2]) - (v1[2] * v2[1]);
307     temp[1] = (v1[2] * v2[0]) - (v1[0] * v2[2]);
308     temp[2] = (v1[0] * v2[1]) - (v1[1] * v2[0]);
309     vcopy(temp.ptr, cross);
310 }
311 
312 float
313 vlength(float[3] v)
314 {
315     return sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
316 }
317 
318 void
319 vscale(float* v, float div)
320 {
321     v[0] *= div;
322     v[1] *= div;
323     v[2] *= div;
324 }
325 
326 void
327 vnormal(float[3] v)
328 {
329     vscale(v.ptr,1.0/vlength(v));
330 }
331 
332 float
333 vdot(float* v1, float* v2)
334 {
335     return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
336 }
337 
338 void
339 vadd(float* src1, float* src2, float* dst)
340 {
341     dst[0] = src1[0] + src2[0];
342     dst[1] = src1[1] + src2[1];
343     dst[2] = src1[2] + src2[2];
344 }
345