home *** CD-ROM | disk | FTP | other *** search
/ Photo CD Demo 1 / Demo.bin / gems / gemsiii / hemis.c < prev    next >
C/C++ Source or Header  |  1992-04-07  |  7KB  |  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, j;
  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 */double projectarea(vect)
  74. Vector3 *vect[PONPATCH];
  75. {
  76.     int i;
  77.     double out= 0, twopntarea();
  78.     for (i = 0; i < PONPATCH; i++)
  79.     out += twopntarea(vect[i], vect[(i+1)%PONPATCH]);
  80.     return(out);
  81. }
  82. /* Calculate the area of the ellipsoid triangle */double twopntarea(vect1, vect2)
  83. Vector3 *vect1, *vect2;
  84. {
  85.     double a, b;
  86.     double sum, cos_sita, halfsita, halfa, d1, d2;
  87.     double arccos(), arctan();
  88.     Vector3 *major, *norm, *v;
  89. /* Find out the normal vector of the great circle */
  90.     major = V3New(0, 0, 1.0);
  91.     norm = V3Duplicate(vect1);
  92.     v = NEWTYPE(Vector3);
  93.     norm = V3Normalize(V3Cross(norm, vect2, v));
  94.     cos_sita = ABS(V3Dot(major, norm)); 
  95.     if (cos_sita < ERR){
  96. /* if the normal vector is perpendicular  to z axis, form factor is zero */
  97.         sum = 0.0;
  98.         free(norm);
  99.         free(major);
  100.         return(sum);
  101.     } 
  102.     halfsita = cos_sita / 2.0;
  103.     halfa = PI * halfsita;
  104.     if (cos_sita < ONE){
  105. /* project the great circle on to the equatorial plane and calculate the form factor */
  106.         major = V3Normalize(V3Cross(major, norm, v));
  107.         norm->x *= -1; norm->y *= -1; norm->z = 0;
  108.         d1 = sqrt(V3Dot(vect1, vect1));
  109.         d2 = sqrt(V3Dot(vect2, vect2));
  110.         a  = V3Dot(major, vect1) / d1;
  111.         if (ABS(a) > 1.0) a /= ABS(a);
  112.         b = V3Dot(norm, vect1);
  113.         sum  = arccos(a, b);
  114.         a = V3Dot(major, vect2) / d2;
  115.         if (ABS(a) > 1.0) a /= ABS(a);
  116.         b  = V3Dot(norm, vect2);
  117.         sum  -= arccos(a, b);
  118.         }
  119.     else{
  120. /* if the normal vector is parallel to z axis, form factor is maximum */
  121.         sum  = arctan(vect1->x, vect1->y);
  122.         sum -= arctan(vect2->x, vect2->y); 
  123.         }
  124.     sum *= halfsita;
  125.     free(norm); free(major); free(v);
  126.     if (sum > halfa) sum -= 2.0 * halfa;
  127.     else if (sum < -halfa) sum += 2.0 * halfa;
  128.     return(sum);
  129. }
  130.  
  131. double arccos(x, y)
  132. double x, y;
  133. {
  134.     double angle;
  135.     angle = acos(x);
  136.     if (y < 0) angle = 2 * PI - angle;
  137.     return(angle);
  138. }
  139.  
  140. double arctan(x, y)
  141. double x, y;
  142. {
  143.     double angle;
  144.     angle = atan2(y, x);
  145.     if (y < 0) angle = 2 * PI + angle;
  146.     return(angle);
  147. }
  148.  
  149. /* Derive the homogenous transformation matrix for the given patch */
  150. void ghmgn(suf, hm) 
  151. Patch *suf;
  152. double hm[HOMON][HOMON];
  153. {
  154.     int i, j;
  155.     Vector3 *v1;
  156.     Point3 *pnt = NEWTYPE(Vector3);
  157.     double a, b, c, d, x1, y1, z1;
  158.     double tt[HOMON][HOMON], tr1[HOMON][HOMON];
  159.     double tr2[HOMON][HOMON], rr[HOMON][HOMON];
  160.     void mtxprd();
  161.  
  162.     pnt->x = (suf->pt[0].x+suf->pt[1].x+suf->pt[2].x)/PONPATCH;
  163.     pnt->y = (suf->pt[0].y+suf->pt[1].y+suf->pt[2].y)/PONPATCH;
  164.     pnt->z = (suf->pt[0].z+suf->pt[1].z+suf->pt[2].z)/PONPATCH;
  165.     for (i = 0; i < HOMON; i++) 
  166.         for (j = 0; j < HOMON; j++)
  167.             if (i == j){
  168.                 tt[i][j] = 1.0; tr1[i][j] = 1.0; tr2[i][j] = 1.0; }
  169.             else{
  170.                 tt[i][j] = 0; tr1[i][j] = 0; tr2[i][j] = 0; }
  171.     v1 = patchnorm(suf);
  172.     tt[3][0] = -pnt->x; tt[3][1] = -pnt->y; tt[3][2] = -pnt->z;
  173.     free(pnt);
  174.     d = sqrt(SQR(v1->y)+SQR(v1->z));
  175.     if (d > ERR){
  176.         tr1[0][0] = 1;
  177.         tr1[3][3] = 1;
  178.         tr1[1][1] = v1->z/d;
  179.         tr1[2][2] = v1->z/d;
  180.         tr1[2][1] = -v1->y/d;
  181.         tr1[1][2] = v1->y/d;
  182.     }
  183.     mtxprd(tt, tr1, rr);
  184.     tr2[0][0] = d; tr2[2][2] = d; tr2[2][0] = -v1->x; tr2[0][2] = v1->x;
  185.     free(v1);
  186.     mtxprd(rr, tr2, hm);
  187. }
  188.  
  189.  
  190. /* Transform a given patch */
  191. void patchtransf(hm, suf)
  192. double hm[HOMON][HOMON];
  193. Patch *suf;
  194. {
  195.     int i;
  196.     void pnttransf();
  197.     for (i = 0; i < PONPATCH; i++)
  198.         pnttransf(hm, &suf->pt[i]);
  199. }
  200.  
  201. /* Transform a given point */
  202. void pnttransf(hm, pnt)
  203. double hm[HOMON][HOMON];
  204. Point3 *pnt;
  205. {
  206.     int i;
  207.     double hi[HOMON], ho[HOMON];
  208.     void transform();
  209.     hi[0] = pnt->x; hi[1] = pnt->y; hi[2] = pnt->z; hi[3] = 1.;
  210.     transform(hm, hi, ho);
  211.     pnt->x = ho[0] / ho[3]; pnt->y = ho[1] / ho[3]; pnt->z = ho[2] / ho[3];
  212. }
  213.  
  214. /* Transform a given vector */
  215. void vecttransf(hm, vect)
  216. double hm[HOMON][HOMON];
  217. Vector3 *vect;
  218. {
  219.     int i;
  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.