1 /* $Id: fixseg.c,v 1.3 2004-12-19 15:21:11 btb Exp $ */
3 THE COMPUTER CODE CONTAINED HEREIN IS THE SOLE PROPERTY OF PARALLAX
4 SOFTWARE CORPORATION ("PARALLAX"). PARALLAX, IN DISTRIBUTING THE CODE TO
5 END-USERS, AND SUBJECT TO ALL OF THE TERMS AND CONDITIONS HEREIN, GRANTS A
6 ROYALTY-FREE, PERPETUAL LICENSE TO SUCH END-USERS FOR USE BY SUCH END-USERS
7 IN USING, DISPLAYING, AND CREATING DERIVATIVE WORKS THEREOF, SO LONG AS
8 SUCH USE, DISPLAY OR CREATION IS FOR NON-COMMERCIAL, ROYALTY OR REVENUE
9 FREE PURPOSES. IN NO EVENT SHALL THE END-USER USE THE COMPUTER CODE
10 CONTAINED HEREIN FOR REVENUE-BEARING PURPOSES. THE END-USER UNDERSTANDS
11 AND AGREES TO THE TERMS HEREIN AND ACCEPTS THE SAME BY USE OF THIS FILE.
12 COPYRIGHT 1993-1998 PARALLAX SOFTWARE CORPORATION. ALL RIGHTS RESERVED.
17 * Functions to make faces planar and probably other things.
22 static char rcsid[] = "$Id: fixseg.c,v 1.3 2004-12-19 15:21:11 btb Exp $";
40 //#include "segment2.h"
45 #define SWAP(a,b) {temp = (a); (a) = (b); (b) = temp;}
47 // -----------------------------------------------------------------------------------------------------------------
48 // Gauss-Jordan elimination solution of a system of linear equations.
49 // a[1..n][1..n] is the input matrix. b[1..n][1..m] is input containing the m right-hand side vectors.
50 // On output, a is replaced by its matrix inverse and b is replaced by the corresponding set of solution vectors.
51 void gaussj(fix **a, int n, fix **b, int m)
53 int indxc[4], indxr[4], ipiv[4];
54 int i, icol=0, irow=0, j, k, l, ll;
55 fix big, dum, pivinv, temp;
58 mprintf((0,"Error -- array too large in gaussj.\n"));
65 for (i=1; i<=n; i++) {
69 for (k=1; k<=n; k++) {
71 if (abs(a[j][k]) >= big) {
76 } else if (ipiv[k] > 1) {
77 mprintf((0,"Error: Singular matrix-1\n"));
84 // We now have the pivot element, so we interchange rows, if needed, to put the pivot
85 // element on the diagonal. The columns are not physically interchanged, only relabeled:
86 // indxc[i], the column of the ith pivot element, is the ith column that is reduced, while
87 // indxr[i] is the row in which that pivot element was originally located. If indxr[i] !=
88 // indxc[i] there is an implied column interchange. With this form of bookkeeping, the
89 // solution b's will end up in the correct order, and the inverse matrix will be scrambled
94 SWAP(a[irow][l], a[icol][l]);
96 SWAP(b[irow][l], b[icol][l]);
101 if (a[icol][icol] == 0) {
102 mprintf((0,"Error: Singular matrix-2\n"));
105 pivinv = fixdiv(F1_0, a[icol][icol]);
106 a[icol][icol] = F1_0;
109 a[icol][l] = fixmul(a[icol][l], pivinv);
111 b[icol][l] = fixmul(b[icol][l], pivinv);
113 for (ll=1; ll<=n; ll++)
118 a[ll][l] -= a[icol][l]*dum;
120 b[ll][l] -= b[icol][l]*dum;
124 // This is the end of the main loop over columns of the reduction. It only remains to unscramble
125 // the solution in view of the column interchanges. We do this by interchanging pairs of
126 // columns in the reverse order that the permutation was built up.
127 for (l=n; l>=1; l--) {
128 if (indxr[l] != indxc[l])
130 SWAP(a[k][indxr[l]], a[k][indxc[l]]);
136 // -----------------------------------------------------------------------------------------------------------------
137 // Return true if side is planar, else return false.
138 int side_is_planar_p(segment *sp, int side)
141 vms_vector *v0,*v1,*v2,*v3;
144 vp = Side_to_verts[side];
145 v0 = &Vertices[sp->verts[vp[0]]];
146 v1 = &Vertices[sp->verts[vp[1]]];
147 v2 = &Vertices[sp->verts[vp[2]]];
148 v3 = &Vertices[sp->verts[vp[3]]];
150 vm_vec_normalize(vm_vec_normal(&va,v0,v1,v2));
151 vm_vec_normalize(vm_vec_normal(&vb,v0,v2,v3));
153 // If the two vectors are very close to being the same, then generate one quad, else generate two triangles.
154 return (vm_vec_dist(&va,&vb) < F1_0/1000);
157 // -------------------------------------------------------------------------------------------------
158 // Return coordinates of a vertex which is vertex v moved so that all sides of which it is a part become planar.
159 void compute_planar_vert(segment *sp, int side, int v, vms_vector *vp)
161 if ((sp) && (side > -3))
165 // -------------------------------------------------------------------------------------------------
166 // Making Cursegp:Curside planar.
167 // If already planar, return.
168 // for each vertex v on side, not part of another segment
169 // choose the vertex v which can be moved to make all sides of which it is a part planar, minimizing distance moved
170 // if there is no vertex v on side, not part of another segment, give up in disgust
172 // 0 curside made planar (or already was)
173 // 1 did not make curside planar
174 int make_curside_planar(void)
178 vms_vector planar_verts[4]; // store coordinates of up to 4 vertices which will make Curside planar, corresponding to each of 4 vertices on side
179 int present_verts[4]; // set to 1 if vertex is present
181 if (side_is_planar_p(Cursegp, Curside))
184 // Look at all vertices in side to find a free one.
186 present_verts[v] = 0;
188 vp = Side_to_verts[Curside];
190 for (v=0; v<4; v++) {
191 int v1 = vp[v]; // absolute vertex id
192 if (is_free_vertex(Cursegp->verts[v1])) {
193 compute_planar_vert(Cursegp, Curside, Cursegp->verts[v1], &planar_verts[v]);
194 present_verts[v] = 1;
198 // Now, for each v for which present_verts[v] == 1, there is a vector (point) in planar_verts[v].
199 // See which one is closest to the plane defined by the other three points.
200 // Nah...just use the first one we find.
202 if (present_verts[v]) {
203 med_set_vertex(vp[v],&planar_verts[v]);
204 validate_segment(Cursegp);
205 // -- should propagate tmaps to segments or something here...
209 // We tried, but we failed, to make Curside planer.