home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / os2 / rsg_ray / quad.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-02-09  |  11.6 KB  |  450 lines
  1. /*
  2.  
  3. QUAD.C  Arbitrary quadratic logic
  4.  
  5.                2   2   2
  6. Defined solid for :- ax +by +cz +dxy+eyz+fzx+gx+hy+iz+j<=0
  7.  
  8. Ellisoids, spheres, infinitely long cylinders, double cones and planes are
  9. all special case of this equation.
  10.  
  11. */
  12.  
  13. /*...sincludes:0:*/
  14. #include <stdio.h>
  15. #include <stdlib.h>
  16. #include <stddef.h>
  17. #include <malloc.h>
  18. #include <memory.h>
  19. #include <math.h>
  20. #include "standard.h"
  21. #include "rt.h"
  22. #include "vector.h"
  23. #include "sil.h"
  24. #define    _QUAD_
  25. #include "quad.h"
  26.  
  27. /*...vrt\46\h:0:*/
  28. /*...vvector\46\h:0:*/
  29. /*...vsil\46\h:0:*/
  30. /*...vquad\46\h:0:*/
  31. /*...e*/
  32.  
  33. /*...screate_quad        \45\ create any general quadratic:0:*/
  34. QUAD    *create_quad(
  35.     double a, double b, double c,
  36.     double d, double e, double f,
  37.     double g, double h, double i,
  38.     double j
  39.     )
  40.     {
  41.     QUAD    *quad;
  42.  
  43.     if ( (quad = malloc(sizeof(QUAD))) == NULL )
  44.         return ( NULL );
  45.  
  46.     quad -> a = a;
  47.     quad -> b = b;
  48.     quad -> c = c;
  49.     quad -> d = d;
  50.     quad -> e = e;
  51.     quad -> f = f;
  52.     quad -> g = g;
  53.     quad -> h = h;
  54.     quad -> i = i;
  55.     quad -> j = j;
  56.  
  57.     return ( quad );
  58.     }
  59. /*...e*/
  60. /*...screate_ellipsoid   \45\ create an ellipsoid of a given size at the origin:0:*/
  61. QUAD    *create_ellipsoid(double rx, double ry, double rz)
  62.     {
  63.     return ( create_quad(1.0 / (rx * rx), 1.0 / (ry * ry), 1.0 / (rz * rz),
  64.                  0.0, 0.0, 0.0,
  65.                  0.0, 0.0, 0.0,
  66.                  -1.0) );
  67.     }
  68. /*...e*/
  69. /*...screate_x_ell_cyl   \45\ create elliptical cylinder along x axis of given radii:0:*/
  70. QUAD    *create_x_ell_cyl(double ry, double rz)
  71.     {
  72.     return ( create_quad(0.0, 1.0 / (ry * ry), 1.0 / (rz * rz),
  73.                  0.0, 0.0, 0.0,
  74.                  0.0, 0.0, 0.0,
  75.                  -1.0) );
  76.     }
  77. /*...e*/
  78. /*...screate_y_ell_cyl   \45\ create elliptical cylinder along y axis of given radii:0:*/
  79. QUAD    *create_y_ell_cyl(double rx, double rz)
  80.     {
  81.     return ( create_quad(1.0 / (rx * rx), 0.0, 1.0 / (rz * rz),
  82.                  0.0, 0.0, 0.0,
  83.                  0.0, 0.0, 0.0,
  84.                  -1.0) );
  85.     }
  86. /*...e*/
  87. /*...screate_z_ell_cyl   \45\ create elliptical cylinder along z axis of given radii:0:*/
  88. QUAD    *create_z_ell_cyl(double rx, double ry)
  89.     {
  90.     return ( create_quad(1.0 / (rx * rx), 1.0 / (ry * ry), 0.0,
  91.                  0.0, 0.0, 0.0,
  92.                  0.0, 0.0, 0.0,
  93.                  -1.0) );
  94.     }
  95. /*...e*/
  96. /*...screate_x_cyl       \45\ create cylinder along x axis of given radii:0:*/
  97. QUAD    *create_x_cyl(double r)
  98.     {
  99.     return ( create_x_ell_cyl(r, r) );
  100.     }
  101. /*...e*/
  102. /*...screate_y_cyl       \45\ create cylinder along y axis of given radii:0:*/
  103. QUAD    *create_y_cyl(double r)
  104.     {
  105.     return ( create_y_ell_cyl(r, r) );
  106.     }
  107. /*...e*/
  108. /*...screate_z_cyl       \45\ create cylinder along z axis of given radii:0:*/
  109. QUAD    *create_z_cyl(double r)
  110.     {
  111.     return ( create_z_ell_cyl(r, r) );
  112.     }
  113. /*...e*/
  114. /*...screate_x_ell_cone  \45\ create elliptical cone along x axis of given radii:0:*/
  115. QUAD    *create_x_ell_cone(double ky, double kz)
  116.     {
  117.     return ( create_quad(-1.0, ky * ky, kz * kz,
  118.                  0.0, 0.0, 0.0,
  119.                  0.0, 0.0, 0.0,
  120.                  0.0) );
  121.     }
  122. /*...e*/
  123. /*...screate_y_ell_cone  \45\ create elliptical cone along y axis of given radii:0:*/
  124. QUAD    *create_y_ell_cone(double kx, double kz)
  125.     {
  126.     return ( create_quad(kx * kx, -1.0, kz * kz,
  127.                  0.0, 0.0, 0.0,
  128.                  0.0, 0.0, 0.0,
  129.                  0.0) );
  130.     }
  131. /*...e*/
  132. /*...screate_z_ell_cone  \45\ create elliptical cone along z axis of given radii:0:*/
  133. QUAD    *create_z_ell_cone(double kx, double ky)
  134.     {
  135.     return ( create_quad(kx * kx, ky * ky, -1.0,
  136.                  0.0, 0.0, 0.0,
  137.                  0.0, 0.0, 0.0,
  138.                  0.0) );
  139.     }
  140. /*...e*/
  141. /*...screate_x_cone      \45\ create cone along x axis of given radii:0:*/
  142. QUAD    *create_x_cone(double k)
  143.     {
  144.     return ( create_x_ell_cone(k, k) );
  145.     }
  146. /*...e*/
  147. /*...screate_y_cone      \45\ create cone along y axis of given radii:0:*/
  148. QUAD    *create_y_cone(double k)
  149.     {
  150.     return ( create_y_ell_cone(k, k) );
  151.     }
  152. /*...e*/
  153. /*...screate_z_cone      \45\ create cone along z axis of given radii:0:*/
  154. QUAD    *create_z_cone(double k)
  155.     {
  156.     return ( create_z_ell_cone(k, k) );
  157.     }
  158. /*...e*/
  159. /*...scopy_quad          \45\ make a copy of a quad:0:*/
  160. QUAD    *copy_quad(QUAD *quad)
  161.     {
  162.     QUAD    *copy;
  163.  
  164.     if ( (copy = malloc(sizeof(QUAD))) == NULL )
  165.         return ( NULL );
  166.  
  167.     memcpy(copy, quad, sizeof(QUAD));
  168.     return ( copy );
  169.     }
  170. /*...e*/
  171. /*...sdestroy_quad       \45\ destroy a quad:0:*/
  172. void    destroy_quad(QUAD *quad)
  173.     {
  174.     free(quad);
  175.     }
  176. /*...e*/
  177.  
  178. /*...strans_x_quad       \45\ translate by amount in x direction:0:*/
  179. void    trans_x_quad(QUAD *quad, double t)
  180.     {
  181.     quad -> j += (quad -> a * t * t - quad -> g * t);
  182.     quad -> i -= quad -> f * t;
  183.     quad -> h -= quad -> d * t;
  184.     quad -> g -= 2.0 * quad -> a * t;
  185.     }
  186. /*...e*/
  187. /*...strans_y_quad       \45\ translate by amount in y direction:0:*/
  188. void    trans_y_quad(QUAD *quad, double t)
  189.     {
  190.     quad -> j += (quad -> b * t * t - quad -> h * t);
  191.     quad -> i -= quad -> e * t;
  192.     quad -> h -= 2.0 * quad -> b * t;
  193.     quad -> g -= quad -> d * t;
  194.     }
  195. /*...e*/
  196. /*...strans_z_quad       \45\ translate by amount in z direction:0:*/
  197. void    trans_z_quad(QUAD *quad, double t)
  198.     {
  199.     quad -> j += (quad -> c * t * t - quad -> i * t);
  200.     quad -> i -= 2.0 * quad -> c * t;
  201.     quad -> h -= quad -> e * t;
  202.     quad -> g -= quad -> f * t;
  203.     }
  204. /*...e*/
  205. /*...sscale_x_quad       \45\ scale by factor in x direction:0:*/
  206. void    scale_x_quad(QUAD *quad, double factor)
  207.     {
  208.     quad -> a /= (factor * factor);
  209.     quad -> d /= factor;
  210.     quad -> f /= factor;
  211.     quad -> g /= factor;
  212.     }
  213. /*...e*/
  214. /*...sscale_y_quad       \45\ scale by factor in y direction:0:*/
  215. void    scale_y_quad(QUAD *quad, double factor)
  216.     {
  217.     quad -> b /= (factor * factor);
  218.     quad -> d /= factor;
  219.     quad -> e /= factor;
  220.     quad -> h /= factor;
  221.     }
  222. /*...e*/
  223. /*...sscale_z_quad       \45\ scale by factor in z direction:0:*/
  224. void    scale_z_quad(QUAD *quad, double factor)
  225.     {
  226.     quad -> c /= (factor * factor);
  227.     quad -> e /= factor;
  228.     quad -> f /= factor;
  229.     quad -> i /= factor;
  230.     }
  231. /*...e*/
  232. /*...srot_x_quad         \45\ rotate about x axis by given angle:0:*/
  233. void    rot_x_quad(QUAD *quad, double angle)
  234.     {
  235.     double    b    = quad -> b;
  236.     double    c    = quad -> c;
  237.     double    d    = quad -> d;
  238.     double    e    = quad -> e;
  239.     double    f    = quad -> f;
  240.     double    h    = quad -> h;
  241.     double    i    = quad -> i;
  242.     double    ca   = cos(angle);
  243.     double    sa   = sin(angle);
  244.     double    caca = ca * ca;
  245.     double    sasa = sa * sa;
  246.     double    saca = sa * ca;
  247.  
  248.     quad -> b = b * caca + c * sasa - e * saca;
  249.     quad -> c = b * sasa + c * caca + e * saca;
  250.     quad -> d = d * ca - f * sa;
  251.     quad -> e = 2.0 * (b - c) * saca + e * (caca - sasa);
  252.     quad -> f = d * sa + f * ca;
  253.     quad -> h = h * ca - i * sa;
  254.     quad -> i = h * sa + i * ca;
  255.     }
  256. /*...e*/
  257. /*...srot_y_quad         \45\ rotate about y axis by given angle:0:*/
  258. void    rot_y_quad(QUAD *quad, double angle)
  259.     {
  260.     double    a    = quad -> a;
  261.     double    c    = quad -> c;
  262.     double    d    = quad -> d;
  263.     double    e    = quad -> e;
  264.     double    f    = quad -> f;
  265.     double    g    = quad -> g;
  266.     double    i    = quad -> i;
  267.     double    ca   = cos(angle);
  268.     double    sa   = sin(angle);
  269.     double    caca = ca * ca;
  270.     double    sasa = sa * sa;
  271.     double    saca = sa * ca;
  272.  
  273.     quad -> a = a * caca + c * sasa + f * saca;
  274.     quad -> c = a * sasa + c * caca - f * saca;
  275.     quad -> d = e * sa + d * ca;
  276.     quad -> e = e * ca - d * sa;
  277.     quad -> f = 2.0 * (c - a) * saca + f * (caca - sasa);
  278.     quad -> g = i * sa + g * ca;
  279.     quad -> i = i * ca - g * sa;
  280.     }
  281. /*...e*/
  282. /*...srot_z_quad         \45\ rotate about z axis by given angle:0:*/
  283. void    rot_z_quad(QUAD *quad, double angle)
  284.     {
  285.     double    a    = quad -> a;
  286.     double    b    = quad -> b;
  287.     double    d    = quad -> d;
  288.     double    e    = quad -> e;
  289.     double    f    = quad -> f;
  290.     double    g    = quad -> g;
  291.     double    h    = quad -> h;
  292.     double    ca   = cos(angle);
  293.     double    sa   = sin(angle);
  294.     double    caca = ca * ca;
  295.     double    sasa = sa * sa;
  296.     double    saca = sa * ca;
  297.  
  298.     quad -> a = a * caca + b * sasa - d * saca;
  299.     quad -> b = a * sasa + b * caca + d * saca;
  300.     quad -> d = 2.0 * (a - b) * saca + d * (caca - sasa);
  301.     quad -> e = f * sa + e * ca;
  302.     quad -> f = f * ca - e * sa;
  303.     quad -> g = g * ca - h * sa;
  304.     quad -> h = g * sa + h * ca;
  305.     }
  306. /*...e*/
  307.  
  308. /*...sisects_reqd_quad   \45\ max number of isects we will generate:0:*/
  309. int    isects_reqd_quad(QUAD *quad)
  310.     {
  311.     quad=quad; /* Suppress 'unref arg' compiler warning */
  312.  
  313.     return ( 4 );
  314.     }
  315. /*...e*/
  316. /*...sintersect_quad     \45\ intersect quadratic:0:*/
  317. /*
  318.  
  319. Those who don't like maths should page down a few times very rapidly!
  320.  
  321. Any point along the line we are interested in is of the form p + tq.
  322.                                  ~    ~
  323.       2   2   2
  324. Given:    ax +by +cz +dxy+eyz+fzx+gx+hy+iz+j=0
  325.  
  326. Subs:    x = x +tx        y = y +ty        z = z +tz
  327.          p   q             p   q             p   q
  328.  
  329.        2   2   2                    2
  330. Gives:    (ax +by +cz +dx y +ey z +fz x )t  +
  331.        q   q   q   q q   q q   q q
  332.  
  333.     (2[ax x +by y +cz z ]+d[x y +x y ]+e[y z +y z ]+f[z x +z x ]+gx +hy +iz )t +
  334.          p q   p q   p q     p q  q p     p q  q p     p q  q p    q   q   q
  335.  
  336.        2   2   2
  337.     (ax +by +cz +dx y +ey z +fz x +gx +hy +iz +j) = 0
  338.        p   p   p   p p   p p   p p   p   p   p
  339.  
  340.  
  341.       2   2   2                   
  342. Opt:    ax +by +cz +dx y +ey z +fz x            costs 12x's and 5+'s
  343.       q   q   q   q q   q q   q q
  344.  
  345.     (ax +dy )x +(by +ez )y +(cz +fx )z        costs 9x's and 5+'s
  346.        q   q  q    q   q  q    q   q  q
  347.  
  348.  
  349.       2   2   2
  350. Opt:    ax +by +cz +dx y +ey z +fz x +gx +hy +iz +j    costs 15x's and 9+'s
  351.       p   p   p   p p   p p   p p   p   p   p
  352.  
  353.     (ax +dy +g)x +(by +ez +h)y +(cz +fx +i)z +j    costs 9x's and 9+'s
  354.        p   p    p    p   p    p    p   p    p
  355. */
  356.  
  357. /*...svalue_of_quad:0:*/
  358. /*
  359. Use this to determine value of quadratic at given point in space.
  360.  
  361.           2   2   2
  362. Opt:    ax +by +cz +dxy+eyz+fzx+gx+hy+iz+j        costs 15x's and 9+'s
  363.  
  364.     (ax+dy+g)x+(by+ez+h)y+(cz+fx+i)z+j        costs 9x's and 9+'s
  365.  
  366. */
  367.  
  368. static double value_of_quad(void *shapeinfo, VECTOR v)
  369.     {
  370.     QUAD    *quad = (QUAD *) shapeinfo;
  371.     double    x = v.x, y = v.y, z = v.z;
  372.  
  373.     return ( (quad -> a * x + quad -> d * y + quad -> g) * x +
  374.          (quad -> b * y + quad -> e * z + quad -> h) * y +
  375.          (quad -> c * z + quad -> f * x + quad -> i) * z +
  376.           quad -> j );
  377.     }
  378. /*...e*/
  379.  
  380. void    intersect_quad(QUAD *quad, VECTOR p, VECTOR q, SIL *sil)
  381.     {
  382.     double    a  = quad -> a;
  383.     double    b  = quad -> b;
  384.     double    c  = quad -> c;
  385.     double    d  = quad -> d;
  386.     double    e  = quad -> e;
  387.     double    f  = quad -> f;
  388.     double    g  = quad -> g;
  389.     double    h  = quad -> h;
  390.     double    i  = quad -> i;
  391.     double    j  = quad -> j;
  392.     double    qa = (a * q.x + d * q.y) * q.x +
  393.              (b * q.y + e * q.z) * q.y +
  394.              (c * q.z + f * q.x) * q.z;
  395.     double    xx = (a * p.x * q.x + b * p.y * q.y + c * p.z * q.z);
  396.     double    qb = xx + xx +
  397.              d * (p.x * q.y + q.x * p.y) +
  398.              e * (p.y * q.z + q.y * p.z) +
  399.              f * (p.z * q.x + q.z * p.x) +
  400.              g * q.x + h * q.y + i * q.z;
  401.     double    qc = (a * p.x + d * p.y + g) * p.x +
  402.              (b * p.y + e * p.z + h) * p.y +
  403.              (c * p.z + f * p.x + i) * p.z + j;
  404.  
  405.     intersect_quadratic_t(qa, qb, qc, p, q,
  406.                   (void *) quad, value_of_quad, sil);
  407.     }
  408. /*...e*/
  409. /*...snormal_to_quad     \45\ find normal of surface at a given point:0:*/
  410. /*
  411.  
  412. Use partial derivatives to find normal at a given point.
  413.  
  414. Given:      2   2   2
  415.     ax +by +cz +dxy+eyz+fzx+gx+hy+iz+j=0
  416.  
  417. d/dx:    2ax+dy+fz+g
  418.  
  419. d/dy:    2by+dx+ez+h
  420.  
  421. d/dz:    2cz+ey+fx+i
  422.  
  423. */
  424.  
  425. VECTOR    normal_to_quad(QUAD *quad, VECTOR p)
  426.     {
  427.     double    a = quad -> a;
  428.     double    b = quad -> b;
  429.     double    c = quad -> c;
  430.     double    d = quad -> d;
  431.     double    e = quad -> e;
  432.     double    f = quad -> f;
  433.     double    g = quad -> g;
  434.     double    h = quad -> h;
  435.     double    i = quad -> i;
  436.     double    x = p.x;
  437.     double    y = p.y;
  438.     double    z = p.z;
  439.     double    ax = a * x;
  440.     double    by = b * y;
  441.     double    cz = c * z;
  442.     VECTOR    normal;
  443.  
  444.     normal.x = ax + ax + d * y + f * z + g;
  445.     normal.y = by + by + d * x + e * z + h;
  446.     normal.z = cz + cz + e * y + f * x + i;
  447.     return ( normal );
  448.     }
  449. /*...e*/
  450.