1 /* Copyright (C) 1992-1998 The Geometry Center
2 * Copyright (C) 1998-2000 Stuart Levy, Tamara Munzner, Mark Phillips
3 *
4 * This file is part of Geomview.
5 *
6 * Geomview is free software; you can redistribute it and/or modify it
7 * under the terms of the GNU Lesser General Public License as published
8 * by the Free Software Foundation; either version 2, or (at your option)
9 * any later version.
10 *
11 * Geomview is distributed in the hope that it will be useful, but
12 * WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with Geomview; see the file COPYING. If not, write
18 * to the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139,
19 * USA, or visit http://www.gnu.org.
20 */
21
22 #if 0
23 static char copyright[] = "Copyright (C) 1992-1998 The Geometry Center\n\
24 Copyright (C) 1998-2000 Stuart Levy, Tamara Munzner, Mark Phillips";
25 #endif
26
27 #include <math.h>
28 #include <stdio.h>
29 #include "transform.h"
30 #include "hpoint3.h"
31 #include "vec4.h"
32
33 int make_rotation(HPoint3 *v1, HPoint3 *v2, Transform m);
34
make_tform(HPoint3 * p1,HPoint3 * p2,HPoint3 * p3,Transform m)35 int make_tform(HPoint3 *p1, HPoint3 *p2, HPoint3 *p3, Transform m)
36 /* Generate a Euclidean isometry which moves (0,0,0) -> p1.
37 Furthermore,
38 the vector (1,0,0) will go a vector tx based at p1 and ending at p2.
39 the vector (0,1,0) will go a vector ty based at p1, perpendicular to tx,
40 and lying in the plane of p1-p2-p3, in the "direction" of p3.
41 the vector (0,0,1) will go a vector tz based at p1, perpendicular to tx and ty.
42 The rotational part is computed by make_rotation. See comments there */
43 {
44 HPoint3 nv1, nv2;
45 Transform trans;
46
47 VSUB3(p2, p1, &nv1);
48 VSUB3(p3, p2, &nv2);
49 nv1.w = 1.0; nv2.w = 1.0;
50 NORMALIZE3(&nv1);
51 NORMALIZE3(&nv2);
52 if ( !make_rotation(&nv1, &nv2, m)) /* if collinear, use most recent rotation */
53 {
54 #ifdef DEBUG
55 fprintf(stderr,"Non-independent vectors: make_rotation\n");
56 #endif
57 return(-1);
58 }
59 /* compute translation to take p1->(0,0,0) */
60 TmTranslate(trans, p1->x, p1->y, p1->z);
61 /* first do the rotation, then the translation */
62 TmConcat(m, trans, m);
63 return(0);
64 }
65
66 /* make rotation that takes (1,0,0) to v1; and (0,1,0) to the orthog. proj of v2 onto
67 the perpendicular subspace of v1. Return (0) if v1 and v2 are collinear.
68 */
make_rotation(v1,v2,m)69 int make_rotation(v1, v2, m)
70 HPoint3 *v1, *v2;
71 Transform m;
72 {
73 double a;
74 HPoint3 t1, t2;
75
76 a = VDOT3(v1, v2);
77 if (a > .9999 || a < -.9999) return(0);
78 t1.x = v1->x * a;
79 t1.y = v1->y * a;
80 t1.z = v1->z * a;
81 VSUB3(v2, &t1, &t2) /* now t2 is orthogonal to v1 */
82 NORMALIZE3(&t2);
83 TmIdentity(m);
84 m[0][0] = v1->x; m[0][1] = v1->y; m[0][2] = v1->z;
85 m[1][0] = t2.x; m[1][1] = t2.y; m[1][2] = t2.z;
86 /* last row is cross product of the first two rows */
87 m[2][0] = v1->y*t2.z - v1->z*t2.y;
88 m[2][1] = v1->z*t2.x - v1->x*t2.z;
89 m[2][2] = v1->x*t2.y - v1->y*t2.x;
90
91 return(1);
92
93 }
94