home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 3 / RISC_DISC_3.iso / resources / etexts / gems / gemsiii / hemis.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-12-14  |  6.9 KB  |  252 lines
  1. /*
  2. Hemispherical Projection of a Triangle
  3.  
  4. Buming Bian
  5.  
  6. UT System---Center for High Performance Computing 
  7.  
  8. Austin, Texas
  9.  
  10. */
  11. #include <math.h>
  12. #include <stdio.h>
  13.  
  14. #include "GraphicsGems.h"
  15. #include "GraphicsGems.c"
  16.  
  17.  
  18. #define HOMON 4
  19. #define PONPATCH 3
  20. #define ERR  1.0e-10
  21. #define ONE  1 - ERR
  22. /* types of structures */
  23. typedef struct { 
  24.         Point3 pt[PONPATCH];   /* Vertices coordinates */
  25.         double ref;            /* Reflectivity of the patch */
  26.         double trs;            /* transmittance of the patch */
  27.    } Patch;
  28.  
  29. /* return a patch normal. */
  30. Vector3 *patchnorm(suf)
  31. Patch *suf;
  32. {
  33.     Vector3 *vt1 = NEWTYPE(Vector3); 
  34.     Vector3 *vt2 = NEWTYPE(Vector3);
  35.     Vector3 *norm = NEWTYPE(Vector3);
  36.     vt1 = V3Sub(&suf->pt[1], &suf->pt[0], vt1);
  37.     vt2 = V3Sub(&suf->pt[2], &suf->pt[0], vt2);
  38.     norm = V3Normalize(V3Cross(vt1, vt2, norm));
  39.     free(vt1); free(vt2);
  40.     return(norm);
  41. }
  42. /******************************************************************************
  43. *  Compute the form factor. Homogeneous transformation is applied to the two  *
  44. *  input Patches such that the center of patch suf1 located at origin with    *
  45. *  its normal coinciding with z axis. A hemisphere is constructed above the xy* 
  46. *  plane, patch suf2 is projected on the hemisphere to form a spherical       * 
  47. *  triangle. This spherical triangle is again projected on the xy plane and   * 
  48. *  the projected area is the form factor of patch suf2 to patch suf1.         * 
  49. ******************************************************************************/
  50. double hemisphere(suf1, suf2)
  51. Patch *suf1, *suf2;
  52. {
  53.     int i;
  54.     double out, hm[HOMON][HOMON], projectarea();
  55.     Vector3 *vect[PONPATCH];
  56.     void patchtransf(), ghmgn();
  57.  
  58. /* use patch suf1 to form the homogeneous transformation matrix */
  59.  
  60.     ghmgn(suf1, hm);
  61. /* apply the transformation to patch suf2 */
  62.     patchtransf(hm, suf2);
  63. /* represent the vertices with vectors */
  64.     for (i = 0; i < PONPATCH; i++)
  65.         vect[i] = V3New(suf2->pt[i].x, suf2->pt[i].y, suf2->pt[i].z);
  66. /* Normalize the form factor so the all sub-total of the form factor */
  67.     out = projectarea(vect);
  68.     for (i = 0; i < PONPATCH; i++)
  69.         free(vect[i]);
  70.     return(ABS(out/PI));
  71. }
  72. /* Calculate the form factor by adding up all the ellipse sectors formed by two 
  73.  * neighbor points */
  74. double projectarea(vect)
  75. Vector3 *vect[PONPATCH];
  76. {
  77.     int i;
  78.     double out= 0, twopntarea();
  79.     for (i = 0; i < PONPATCH; i++)
  80.     out += twopntarea(vect[i], vect[(i+1)%PONPATCH]);
  81.     return(out);
  82. }
  83. /* Calculate the area of the ellipsoid triangle */
  84. double twopntarea(vect1, vect2)
  85. Vector3 *vect1, *vect2;
  86. {
  87.     double a, b;
  88.     double sum, cos_sita, halfsita, halfa, d1, d2;
  89.     double arccos(), arctan();
  90.     Vector3 *major, *norm, *v;
  91. /* Find out the normal vector of the great circle */
  92.     major = V3New(0.0, 0.0, 1.0);
  93.     norm = V3Duplicate(vect1);
  94.     v = NEWTYPE(Vector3);
  95.     norm = V3Normalize(V3Cross(norm, vect2, v));
  96.     cos_sita = ABS(V3Dot(major, norm)); 
  97.     if (cos_sita < ERR){
  98. /* if the normal vector is perpendicular  to z axis, form factor is zero */
  99.         sum = 0.0;
  100.         free(norm);
  101.         free(major);
  102.         return(sum);
  103.     } 
  104.     halfsita = cos_sita / 2.0;
  105.     halfa = PI * halfsita;
  106.     if (cos_sita < ONE){
  107. /* project the great circle on to the equatorial plane and calculate the form factor */
  108.         major = V3Normalize(V3Cross(major, norm, v));
  109.         norm->x *= -1; norm->y *= -1; norm->z = 0;
  110.         d1 = sqrt(V3Dot(vect1, vect1));
  111.         d2 = sqrt(V3Dot(vect2, vect2));
  112.         a  = V3Dot(major, vect1) / d1;
  113.         if (ABS(a) > 1.0) a /= ABS(a);
  114.         b = V3Dot(norm, vect1);
  115.         sum  = arccos(a, b);
  116.         a = V3Dot(major, vect2) / d2;
  117.         if (ABS(a) > 1.0) a /= ABS(a);
  118.         b  = V3Dot(norm, vect2);
  119.         sum  -= arccos(a, b);
  120.         }
  121.     else{
  122. /* if the normal vector is parallel to z axis, form factor is maximum */
  123.         sum  = arctan(vect1->x, vect1->y);
  124.         sum -= arctan(vect2->x, vect2->y); 
  125.         }
  126.     sum *= halfsita;
  127.     free(norm); free(major); free(v);
  128.     if (sum > halfa) sum -= 2.0 * halfa;
  129.     else if (sum < -halfa) sum += 2.0 * halfa;
  130.     return(sum);
  131. }
  132.  
  133. double arccos(x, y)
  134. double x, y;
  135. {
  136.     double angle;
  137.     angle = acos(x);
  138.     if (y < 0) angle = 2 * PI - angle;
  139.     return(angle);
  140. }
  141.  
  142. double arctan(x, y)
  143. double x, y;
  144. {
  145.     double angle;
  146.     angle = atan2(y, x);
  147.     if (y < 0) angle = 2 * PI + angle;
  148.     return(angle);
  149. }
  150.  
  151. /* Derive the homogeneous transformation matrix for the given patch */
  152. void ghmgn(suf, hm) 
  153. Patch *suf;
  154. double hm[HOMON][HOMON];
  155. {
  156.     int i, j;
  157.     Vector3 *v1;
  158.     Point3 *pnt = NEWTYPE(Vector3);
  159.     double d;
  160.     double tt[HOMON][HOMON], tr1[HOMON][HOMON];
  161.     double tr2[HOMON][HOMON], rr[HOMON][HOMON];
  162.     void mtxprd();
  163.  
  164.     pnt->x = (suf->pt[0].x+suf->pt[1].x+suf->pt[2].x)/PONPATCH;
  165.     pnt->y = (suf->pt[0].y+suf->pt[1].y+suf->pt[2].y)/PONPATCH;
  166.     pnt->z = (suf->pt[0].z+suf->pt[1].z+suf->pt[2].z)/PONPATCH;
  167.     for (i = 0; i < HOMON; i++) 
  168.         for (j = 0; j < HOMON; j++)
  169.             if (i == j){
  170.                 tt[i][j] = 1.0; tr1[i][j] = 1.0; tr2[i][j] = 1.0; }
  171.             else{
  172.                 tt[i][j] = 0; tr1[i][j] = 0; tr2[i][j] = 0; }
  173.     v1 = patchnorm(suf);
  174.     tt[3][0] = -pnt->x; tt[3][1] = -pnt->y; tt[3][2] = -pnt->z;
  175.     free(pnt);
  176.     d = sqrt(SQR(v1->y)+SQR(v1->z));
  177.     if (d > ERR){
  178.         tr1[0][0] = 1;
  179.         tr1[3][3] = 1;
  180.         tr1[1][1] = v1->z/d;
  181.         tr1[2][2] = v1->z/d;
  182.         tr1[2][1] = -v1->y/d;
  183.         tr1[1][2] = v1->y/d;
  184.     }
  185.     mtxprd(tt, tr1, rr);
  186.     tr2[0][0] = d; tr2[2][2] = d; tr2[2][0] = -v1->x; tr2[0][2] = v1->x;
  187.     free(v1);
  188.     mtxprd(rr, tr2, hm);
  189. }
  190.  
  191.  
  192. /* Transform a given patch */
  193. void patchtransf(hm, suf)
  194. double hm[HOMON][HOMON];
  195. Patch *suf;
  196. {
  197.     int i;
  198.     void pnttransf();
  199.     for (i = 0; i < PONPATCH; i++)
  200.         pnttransf(hm, &suf->pt[i]);
  201. }
  202.  
  203. /* Transform a given point */
  204. void pnttransf(hm, pnt)
  205. double hm[HOMON][HOMON];
  206. Point3 *pnt;
  207. {
  208.     double hi[HOMON], ho[HOMON];
  209.     void transform();
  210.     hi[0] = pnt->x; hi[1] = pnt->y; hi[2] = pnt->z; hi[3] = 1.;
  211.     transform(hm, hi, ho);
  212.     pnt->x = ho[0] / ho[3]; pnt->y = ho[1] / ho[3]; pnt->z = ho[2] / ho[3];
  213. }
  214.  
  215. /* Transform a given vector */
  216. void vecttransf(hm, vect)
  217. double hm[HOMON][HOMON];
  218. Vector3 *vect;
  219. {
  220.     double hi[HOMON], ho[HOMON];
  221.     void transform();
  222.     hi[0] = vect->x; hi[1] = vect->y; hi[2] = vect->z; hi[3] = 0.;
  223.     transform(hm, hi, ho);
  224.     vect->x = ho[0]; vect->y = ho[1]; vect->z = ho[2];
  225. }
  226.  
  227. /* Transform a given vector */
  228. void transform(hm, hi, ho)
  229. double hm[HOMON][HOMON];
  230. double hi[HOMON], ho[HOMON];
  231. {
  232.     int j, k;
  233.     for (j = 0; j < HOMON; j++){
  234.         ho[j] = 0;
  235.         for (k = 0; k < HOMON; k++)
  236.             ho[j] += hi[k] * hm[k][j];
  237.     }  
  238. }
  239.  
  240. /* Multiple two matrices */
  241. void mtxprd(h1, h2, h3)
  242. double h1[HOMON][HOMON], h2[HOMON][HOMON],h3[HOMON][HOMON];
  243. {
  244.     int i, j, k;
  245.     for (i = 0; i < HOMON; i++)
  246.         for (j = 0; j < HOMON; j++){
  247.             h3[i][j] = 0;
  248.             for (k = 0; k < HOMON; k++)
  249.                 h3[i][j] += h1[i][k] * h2[k][j];
  250.         }
  251. }
  252.