geom.c 8.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266
  1. /*
  2. * SGI FREE SOFTWARE LICENSE B (Version 2.0, Sept. 18, 2008)
  3. * Copyright (C) 1991-2000 Silicon Graphics, Inc. All Rights Reserved.
  4. *
  5. * Permission is hereby granted, free of charge, to any person obtaining a
  6. * copy of this software and associated documentation files (the "Software"),
  7. * to deal in the Software without restriction, including without limitation
  8. * the rights to use, copy, modify, merge, publish, distribute, sublicense,
  9. * and/or sell copies of the Software, and to permit persons to whom the
  10. * Software is furnished to do so, subject to the following conditions:
  11. *
  12. * The above copyright notice including the dates of first publication and
  13. * either this permission notice or a reference to
  14. * http://oss.sgi.com/projects/FreeB/
  15. * shall be included in all copies or substantial portions of the Software.
  16. *
  17. * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  18. * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  19. * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  20. * SILICON GRAPHICS, INC. BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
  21. * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF
  22. * OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  23. * SOFTWARE.
  24. *
  25. * Except as contained in this notice, the name of Silicon Graphics, Inc.
  26. * shall not be used in advertising or otherwise to promote the sale, use or
  27. * other dealings in this Software without prior written authorization from
  28. * Silicon Graphics, Inc.
  29. */
  30. /*
  31. ** Author: Eric Veach, July 1994.
  32. **
  33. */
  34. #include "../prboom/SDL_opengl.h" // JDC
  35. //#include "gluos.h"
  36. #include <assert.h>
  37. #include "mesh.h"
  38. #include "geom.h"
  39. int __gl_vertLeq( GLUvertex *u, GLUvertex *v )
  40. {
  41. /* Returns TRUE if u is lexicographically <= v. */
  42. return VertLeq( u, v );
  43. }
  44. GLdouble __gl_edgeEval( GLUvertex *u, GLUvertex *v, GLUvertex *w )
  45. {
  46. /* Given three vertices u,v,w such that VertLeq(u,v) && VertLeq(v,w),
  47. * evaluates the t-coord of the edge uw at the s-coord of the vertex v.
  48. * Returns v->t - (uw)(v->s), ie. the signed distance from uw to v.
  49. * If uw is vertical (and thus passes thru v), the result is zero.
  50. *
  51. * The calculation is extremely accurate and stable, even when v
  52. * is very close to u or w. In particular if we set v->t = 0 and
  53. * let r be the negated result (this evaluates (uw)(v->s)), then
  54. * r is guaranteed to satisfy MIN(u->t,w->t) <= r <= MAX(u->t,w->t).
  55. */
  56. GLdouble gapL, gapR;
  57. assert( VertLeq( u, v ) && VertLeq( v, w ));
  58. gapL = v->s - u->s;
  59. gapR = w->s - v->s;
  60. if( gapL + gapR > 0 ) {
  61. if( gapL < gapR ) {
  62. return (v->t - u->t) + (u->t - w->t) * (gapL / (gapL + gapR));
  63. } else {
  64. return (v->t - w->t) + (w->t - u->t) * (gapR / (gapL + gapR));
  65. }
  66. }
  67. /* vertical line */
  68. return 0;
  69. }
  70. GLdouble __gl_edgeSign( GLUvertex *u, GLUvertex *v, GLUvertex *w )
  71. {
  72. /* Returns a number whose sign matches EdgeEval(u,v,w) but which
  73. * is cheaper to evaluate. Returns > 0, == 0 , or < 0
  74. * as v is above, on, or below the edge uw.
  75. */
  76. GLdouble gapL, gapR;
  77. assert( VertLeq( u, v ) && VertLeq( v, w ));
  78. gapL = v->s - u->s;
  79. gapR = w->s - v->s;
  80. if( gapL + gapR > 0 ) {
  81. return (v->t - w->t) * gapL + (v->t - u->t) * gapR;
  82. }
  83. /* vertical line */
  84. return 0;
  85. }
  86. /***********************************************************************
  87. * Define versions of EdgeSign, EdgeEval with s and t transposed.
  88. */
  89. GLdouble __gl_transEval( GLUvertex *u, GLUvertex *v, GLUvertex *w )
  90. {
  91. /* Given three vertices u,v,w such that TransLeq(u,v) && TransLeq(v,w),
  92. * evaluates the t-coord of the edge uw at the s-coord of the vertex v.
  93. * Returns v->s - (uw)(v->t), ie. the signed distance from uw to v.
  94. * If uw is vertical (and thus passes thru v), the result is zero.
  95. *
  96. * The calculation is extremely accurate and stable, even when v
  97. * is very close to u or w. In particular if we set v->s = 0 and
  98. * let r be the negated result (this evaluates (uw)(v->t)), then
  99. * r is guaranteed to satisfy MIN(u->s,w->s) <= r <= MAX(u->s,w->s).
  100. */
  101. GLdouble gapL, gapR;
  102. assert( TransLeq( u, v ) && TransLeq( v, w ));
  103. gapL = v->t - u->t;
  104. gapR = w->t - v->t;
  105. if( gapL + gapR > 0 ) {
  106. if( gapL < gapR ) {
  107. return (v->s - u->s) + (u->s - w->s) * (gapL / (gapL + gapR));
  108. } else {
  109. return (v->s - w->s) + (w->s - u->s) * (gapR / (gapL + gapR));
  110. }
  111. }
  112. /* vertical line */
  113. return 0;
  114. }
  115. GLdouble __gl_transSign( GLUvertex *u, GLUvertex *v, GLUvertex *w )
  116. {
  117. /* Returns a number whose sign matches TransEval(u,v,w) but which
  118. * is cheaper to evaluate. Returns > 0, == 0 , or < 0
  119. * as v is above, on, or below the edge uw.
  120. */
  121. GLdouble gapL, gapR;
  122. assert( TransLeq( u, v ) && TransLeq( v, w ));
  123. gapL = v->t - u->t;
  124. gapR = w->t - v->t;
  125. if( gapL + gapR > 0 ) {
  126. return (v->s - w->s) * gapL + (v->s - u->s) * gapR;
  127. }
  128. /* vertical line */
  129. return 0;
  130. }
  131. int __gl_vertCCW( GLUvertex *u, GLUvertex *v, GLUvertex *w )
  132. {
  133. /* For almost-degenerate situations, the results are not reliable.
  134. * Unless the floating-point arithmetic can be performed without
  135. * rounding errors, *any* implementation will give incorrect results
  136. * on some degenerate inputs, so the client must have some way to
  137. * handle this situation.
  138. */
  139. return (u->s*(v->t - w->t) + v->s*(w->t - u->t) + w->s*(u->t - v->t)) >= 0;
  140. }
  141. /* Given parameters a,x,b,y returns the value (b*x+a*y)/(a+b),
  142. * or (x+y)/2 if a==b==0. It requires that a,b >= 0, and enforces
  143. * this in the rare case that one argument is slightly negative.
  144. * The implementation is extremely stable numerically.
  145. * In particular it guarantees that the result r satisfies
  146. * MIN(x,y) <= r <= MAX(x,y), and the results are very accurate
  147. * even when a and b differ greatly in magnitude.
  148. */
  149. #define RealInterpolate(a,x,b,y) \
  150. (a = (a < 0) ? 0 : a, b = (b < 0) ? 0 : b, \
  151. ((a <= b) ? ((b == 0) ? ((x+y) / 2) \
  152. : (x + (y-x) * (a/(a+b)))) \
  153. : (y + (x-y) * (b/(a+b)))))
  154. #ifndef FOR_TRITE_TEST_PROGRAM
  155. #define Interpolate(a,x,b,y) RealInterpolate(a,x,b,y)
  156. #else
  157. /* Claim: the ONLY property the sweep algorithm relies on is that
  158. * MIN(x,y) <= r <= MAX(x,y). This is a nasty way to test that.
  159. */
  160. #include <stdlib.h>
  161. extern int RandomInterpolate;
  162. GLdouble Interpolate( GLdouble a, GLdouble x, GLdouble b, GLdouble y)
  163. {
  164. printf("*********************%d\n",RandomInterpolate);
  165. if( RandomInterpolate ) {
  166. a = 1.2 * drand48() - 0.1;
  167. a = (a < 0) ? 0 : ((a > 1) ? 1 : a);
  168. b = 1.0 - a;
  169. }
  170. return RealInterpolate(a,x,b,y);
  171. }
  172. #endif
  173. #define Swap(a,b) if (1) { GLUvertex *t = a; a = b; b = t; } else
  174. void __gl_edgeIntersect( GLUvertex *o1, GLUvertex *d1,
  175. GLUvertex *o2, GLUvertex *d2,
  176. GLUvertex *v )
  177. /* Given edges (o1,d1) and (o2,d2), compute their point of intersection.
  178. * The computed point is guaranteed to lie in the intersection of the
  179. * bounding rectangles defined by each edge.
  180. */
  181. {
  182. GLdouble z1, z2;
  183. /* This is certainly not the most efficient way to find the intersection
  184. * of two line segments, but it is very numerically stable.
  185. *
  186. * Strategy: find the two middle vertices in the VertLeq ordering,
  187. * and interpolate the intersection s-value from these. Then repeat
  188. * using the TransLeq ordering to find the intersection t-value.
  189. */
  190. if( ! VertLeq( o1, d1 )) { Swap( o1, d1 ); }
  191. if( ! VertLeq( o2, d2 )) { Swap( o2, d2 ); }
  192. if( ! VertLeq( o1, o2 )) { Swap( o1, o2 ); Swap( d1, d2 ); }
  193. if( ! VertLeq( o2, d1 )) {
  194. /* Technically, no intersection -- do our best */
  195. v->s = (o2->s + d1->s) / 2;
  196. } else if( VertLeq( d1, d2 )) {
  197. /* Interpolate between o2 and d1 */
  198. z1 = EdgeEval( o1, o2, d1 );
  199. z2 = EdgeEval( o2, d1, d2 );
  200. if( z1+z2 < 0 ) { z1 = -z1; z2 = -z2; }
  201. v->s = Interpolate( z1, o2->s, z2, d1->s );
  202. } else {
  203. /* Interpolate between o2 and d2 */
  204. z1 = EdgeSign( o1, o2, d1 );
  205. z2 = -EdgeSign( o1, d2, d1 );
  206. if( z1+z2 < 0 ) { z1 = -z1; z2 = -z2; }
  207. v->s = Interpolate( z1, o2->s, z2, d2->s );
  208. }
  209. /* Now repeat the process for t */
  210. if( ! TransLeq( o1, d1 )) { Swap( o1, d1 ); }
  211. if( ! TransLeq( o2, d2 )) { Swap( o2, d2 ); }
  212. if( ! TransLeq( o1, o2 )) { Swap( o1, o2 ); Swap( d1, d2 ); }
  213. if( ! TransLeq( o2, d1 )) {
  214. /* Technically, no intersection -- do our best */
  215. v->t = (o2->t + d1->t) / 2;
  216. } else if( TransLeq( d1, d2 )) {
  217. /* Interpolate between o2 and d1 */
  218. z1 = TransEval( o1, o2, d1 );
  219. z2 = TransEval( o2, d1, d2 );
  220. if( z1+z2 < 0 ) { z1 = -z1; z2 = -z2; }
  221. v->t = Interpolate( z1, o2->t, z2, d1->t );
  222. } else {
  223. /* Interpolate between o2 and d2 */
  224. z1 = TransSign( o1, o2, d1 );
  225. z2 = -TransSign( o1, d2, d1 );
  226. if( z1+z2 < 0 ) { z1 = -z1; z2 = -z2; }
  227. v->t = Interpolate( z1, o2->t, z2, d2->t );
  228. }
  229. }