home *** CD-ROM | disk | FTP | other *** search
/ PC Pro 2002 April / pcpro0402.iso / essentials / graphics / Gimp / gimp-src-20001226.exe / src / gimp / plug-ins / MapObject / mapobject_shade.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-06-06  |  34.7 KB  |  1,261 lines
  1. /*****************/
  2. /* Shading stuff */
  3. /*****************/
  4.  
  5. #include <string.h>
  6.  
  7. #include <libgimp/gimp.h>
  8.  
  9. #include <gck/gck.h>
  10.  
  11. #include "mapobject_apply.h"
  12. #include "mapobject_main.h"
  13. #include "mapobject_image.h"
  14. #include "mapobject_shade.h"
  15.  
  16. gdouble            bx1, by1, bx2, by2;
  17. get_ray_color_func get_ray_color;
  18.  
  19. typedef struct
  20. {
  21.   gdouble     u, v;
  22.   gdouble     t;
  23.   GimpVector3 s;
  24.   GimpVector3 n;
  25.   gint        face;
  26. } FaceIntersectInfo;
  27.  
  28. /*****************/
  29. /* Phong shading */
  30. /*****************/
  31.  
  32. static GckRGB
  33. phong_shade (GimpVector3 *pos,
  34.          GimpVector3 *viewpoint,
  35.          GimpVector3 *normal,
  36.          GimpVector3 *light,
  37.          GckRGB      *diff_col,
  38.          GckRGB      *spec_col,
  39.          gint         type)
  40. {
  41.   GckRGB ambientcolor, diffusecolor, specularcolor;
  42.   gdouble NL, RV, dist;
  43.   GimpVector3 L, NN, V, N;
  44.  
  45.   /* Compute ambient intensity */
  46.   /* ========================= */
  47.  
  48.   N = *normal;
  49.   ambientcolor = *diff_col;
  50.   gck_rgb_mul (&ambientcolor, mapvals.material.ambient_int);
  51.  
  52.   /* Compute (N*L) term of Phong's equation */
  53.   /* ====================================== */
  54.  
  55.   if (type == POINT_LIGHT)
  56.     gimp_vector3_sub (&L, light, pos);
  57.   else
  58.     L = *light;
  59.  
  60.   dist = gimp_vector3_length (&L);
  61.  
  62.   if (dist != 0.0)
  63.     gimp_vector3_mul (&L, 1.0 / dist);
  64.  
  65.   NL = 2.0 * gimp_vector3_inner_product (&N, &L);
  66.  
  67.   if (NL >= 0.0)
  68.     {
  69.       /* Compute (R*V)^alpha term of Phong's equation */
  70.       /* ============================================ */
  71.  
  72.       gimp_vector3_sub (&V, viewpoint, pos);
  73.       gimp_vector3_normalize (&V);
  74.  
  75.       gimp_vector3_mul (&N, NL);
  76.       gimp_vector3_sub (&NN, &N, &L);
  77.       RV = gimp_vector3_inner_product (&NN, &V);
  78.       RV = pow (RV, mapvals.material.highlight);
  79.  
  80.       /* Compute diffuse and specular intensity contribution */
  81.       /* =================================================== */
  82.  
  83.       diffusecolor = *diff_col;
  84.       gck_rgb_mul (&diffusecolor, mapvals.material.diffuse_ref);
  85.       gck_rgb_mul (&diffusecolor, NL);
  86.  
  87.       specularcolor = *spec_col;
  88.       gck_rgb_mul (&specularcolor, mapvals.material.specular_ref);
  89.       gck_rgb_mul (&specularcolor, RV);
  90.  
  91.       gck_rgb_add (&diffusecolor, &specularcolor);
  92.       gck_rgb_mul (&diffusecolor, mapvals.material.diffuse_int);
  93.       gck_rgb_clamp (&diffusecolor);
  94.  
  95.       gck_rgb_add (&ambientcolor, &diffusecolor);
  96.     }
  97.  
  98.   return ambientcolor;
  99. }
  100.  
  101. static gint
  102. plane_intersect (GimpVector3 *dir,
  103.          GimpVector3 *viewp,
  104.          GimpVector3 *ipos,
  105.          gdouble     *u,
  106.          gdouble     *v)
  107. {
  108.   static gdouble det, det1, det2, det3, t;
  109.  
  110.   imat[0][0] = dir->x;
  111.   imat[1][0] = dir->y;
  112.   imat[2][0] = dir->z;
  113.  
  114.   /* Compute determinant of the first 3x3 sub matrix (denominator) */
  115.   /* ============================================================= */
  116.  
  117.   det = (imat[0][0] * imat[1][1] * imat[2][2] +
  118.      imat[0][1] * imat[1][2] * imat[2][0] +
  119.      imat[0][2] * imat[1][0] * imat[2][1] -
  120.      imat[0][2] * imat[1][1] * imat[2][0] -
  121.      imat[0][0] * imat[1][2] * imat[2][1] -
  122.      imat[2][2] * imat[0][1] * imat[1][0]);
  123.  
  124.   /* If the determinant is non-zero, a intersection point exists */
  125.   /* =========================================================== */
  126.  
  127.   if (det != 0.0)
  128.     {
  129.       /* Now, lets compute the numerator determinants (wow ;) */
  130.       /* ==================================================== */
  131.  
  132.       det1 = (imat[0][3] * imat[1][1] * imat[2][2] +
  133.           imat[0][1] * imat[1][2] * imat[2][3] +
  134.           imat[0][2] * imat[1][3] * imat[2][1] -
  135.           imat[0][2] * imat[1][1] * imat[2][3] -
  136.           imat[1][2] * imat[2][1] * imat[0][3] -
  137.           imat[2][2] * imat[0][1] * imat[1][3]);
  138.  
  139.       det2 = (imat[0][0] * imat[1][3] * imat[2][2] +
  140.           imat[0][3] * imat[1][2] * imat[2][0] +
  141.           imat[0][2] * imat[1][0] * imat[2][3] -
  142.           imat[0][2] * imat[1][3] * imat[2][0] -
  143.           imat[1][2] * imat[2][3] * imat[0][0] -
  144.           imat[2][2] * imat[0][3] * imat[1][0]);
  145.  
  146.       det3 = (imat[0][0] * imat[1][1] * imat[2][3] +
  147.           imat[0][1] * imat[1][3] * imat[2][0] +
  148.           imat[0][3] * imat[1][0] * imat[2][1] -
  149.           imat[0][3] * imat[1][1] * imat[2][0] -
  150.           imat[1][3] * imat[2][1] * imat[0][0] -
  151.           imat[2][3] * imat[0][1] * imat[1][0]);
  152.  
  153.       /* Now we have the simultanous solutions. Lets compute the unknowns */
  154.       /* (skip u&v if t is <0, this means the intersection is behind us)  */
  155.       /* ================================================================ */
  156.  
  157.       t = det1 / det;
  158.  
  159.       if (t > 0.0)
  160.         {
  161.           *u = 1.0 + ((det2 / det) - 0.5);
  162.           *v = 1.0 + ((det3 / det) - 0.5);
  163.  
  164.             ipos->x = viewp->x + t * dir->x;
  165.           ipos->y = viewp->y + t * dir->y;
  166.           ipos->z = viewp->z + t * dir->z;
  167.  
  168.           return TRUE;
  169.         }
  170.     }
  171.  
  172.   return FALSE;
  173. }
  174.  
  175. /*****************************************************************************
  176.  * These routines computes the color of the surface
  177.  * of the plane at a given point
  178.  *****************************************************************************/
  179.  
  180. GckRGB
  181. get_ray_color_plane (GimpVector3 *pos)
  182. {
  183.   GckRGB color = background;
  184.   static gint inside = FALSE;
  185.   static GimpVector3 ray, spos;
  186.   static gdouble vx, vy;
  187.  
  188.   /* Construct a line from our VP to the point */
  189.   /* ========================================= */
  190.  
  191.   gimp_vector3_sub (&ray, pos, &mapvals.viewpoint);
  192.   gimp_vector3_normalize (&ray);
  193.  
  194.   /* Check for intersection. This is a quasi ray-tracer. */
  195.   /* =================================================== */
  196.  
  197.   if (plane_intersect (&ray, &mapvals.viewpoint, &spos, &vx, &vy) == TRUE)
  198.     {
  199.       color = get_image_color (vx, vy, &inside);
  200.  
  201.       if (color.a!=0.0 && inside == TRUE &&
  202.       mapvals.lightsource.type != NO_LIGHT)
  203.         {
  204.           /* Compute shading at this point */
  205.           /* ============================= */
  206.  
  207.           color = phong_shade (&spos,
  208.                    &mapvals.viewpoint,
  209.                    &mapvals.normal,
  210.                    &mapvals.lightsource.position,
  211.                    &color,
  212.                    &mapvals.lightsource.color,
  213.                    mapvals.lightsource.type);
  214.  
  215.           gck_rgb_clamp (&color);
  216.         }
  217.     }
  218.   
  219.   if (color.a == 0.0)
  220.     color = background;
  221.  
  222.   return color;
  223. }
  224.  
  225. /***********************************************************************/
  226. /* Given the NorthPole, Equator and a third vector (normal) compute    */
  227. /* the conversion from spherical oordinates to image space coordinates */
  228. /***********************************************************************/
  229.  
  230. static void
  231. sphere_to_image (GimpVector3 *normal,
  232.          gdouble     *u,
  233.          gdouble     *v)
  234. {
  235.   static gdouble alpha, fac;
  236.   static GimpVector3 cross_prod;
  237.  
  238.   alpha = acos (-gimp_vector3_inner_product (&mapvals.secondaxis, normal));
  239.  
  240.   *v = alpha / G_PI;
  241.  
  242.   if (*v == 0.0 || *v == 1.0)
  243.     *u = 0.0;
  244.   else
  245.     {
  246.       fac = (gimp_vector3_inner_product (&mapvals.firstaxis, normal) /
  247.          sin (alpha));
  248.  
  249.       /* Make sure that we map to -1.0..1.0 (take care of rounding errors) */
  250.       /* ================================================================= */
  251.  
  252.       fac = CLAMP (fac, -1.0, 1.0);
  253.  
  254.       *u = acos (fac) / (2.0 * G_PI);
  255.  
  256.       cross_prod = gimp_vector3_cross_product (&mapvals.secondaxis,
  257.                            &mapvals.firstaxis);
  258.  
  259.       if (gimp_vector3_inner_product (&cross_prod, normal) < 0.0)
  260.         *u = 1.0 - *u;
  261.     }
  262. }
  263.  
  264. /***************************************************/
  265. /* Compute intersection point with sphere (if any) */
  266. /***************************************************/
  267.  
  268. static gint
  269. sphere_intersect (GimpVector3 *dir,
  270.           GimpVector3 *viewp,
  271.           GimpVector3 *spos1,
  272.           GimpVector3 *spos2)
  273. {
  274.   static gdouble alpha, beta, tau, s1, s2, tmp;
  275.   static GimpVector3 t;
  276.  
  277.   gimp_vector3_sub (&t, &mapvals.position, viewp);
  278.  
  279.   alpha = gimp_vector3_inner_product (dir, &t);
  280.   beta  = gimp_vector3_inner_product (&t, &t);
  281.  
  282.   tau = alpha * alpha - beta + mapvals.radius * mapvals.radius;
  283.  
  284.   if (tau >= 0.0)
  285.     {
  286.       tau = sqrt (tau);
  287.       s1 = alpha + tau;
  288.       s2 = alpha - tau;
  289.  
  290.       if (s2 < s1)
  291.         {
  292.           tmp = s1;
  293.           s1 = s2;
  294.           s2 = tmp;
  295.         }
  296.  
  297.       spos1->x = viewp->x + s1 * dir->x;
  298.       spos1->y = viewp->y + s1 * dir->y;
  299.       spos1->z = viewp->z + s1 * dir->z;
  300.       spos2->x = viewp->x + s2 * dir->x;
  301.       spos2->y = viewp->y + s2 * dir->y;
  302.       spos2->z = viewp->z + s2 * dir->z;
  303.  
  304.       return TRUE ;
  305.     }
  306.  
  307.   return FALSE;
  308. }
  309.  
  310. /*****************************************************************************
  311.  * These routines computes the color of the surface
  312.  * of the sphere at a given point
  313.  *****************************************************************************/
  314.  
  315. GckRGB
  316. get_ray_color_sphere (GimpVector3 *pos)
  317. {
  318.   GckRGB color=background;
  319.   static GckRGB color2;
  320.   static gint inside=FALSE;
  321.   static GimpVector3 normal,ray,spos1,spos2;
  322.   static gdouble vx,vy;
  323.  
  324.   /* Check if ray is within the bounding box */
  325.   /* ======================================= */
  326.   
  327.   if (pos->x<bx1 || pos->x>bx2 || pos->y<by1 || pos->y>by2)
  328.     return(color);
  329.  
  330.   /* Construct a line from our VP to the point */
  331.   /* ========================================= */
  332.  
  333.   gimp_vector3_sub (&ray, pos, &mapvals.viewpoint);
  334.   gimp_vector3_normalize (&ray);
  335.  
  336.   /* Check for intersection. This is a quasi ray-tracer. */
  337.   /* =================================================== */
  338.  
  339.   if (sphere_intersect (&ray, &mapvals.viewpoint, &spos1, &spos2) == TRUE)
  340.     {
  341.       /* Compute spherical to rectangular mapping */
  342.       /* ======================================== */
  343.  
  344.       gimp_vector3_sub (&normal, &spos1, &mapvals.position);
  345.       gimp_vector3_normalize (&normal);
  346.       sphere_to_image (&normal, &vx, &vy);
  347.       color=get_image_color (vx, vy, &inside);
  348.  
  349.       /* Check for total transparency... */
  350.       /* =============================== */
  351.  
  352.       if (color.a < 1.0)
  353.         {
  354.           /* Hey, we can see  through here!      */
  355.           /* Lets see what's on the other side.. */
  356.           /* =================================== */
  357.  
  358.           color = phong_shade (&spos1,
  359.                    &mapvals.viewpoint,
  360.                    &normal,
  361.                    &mapvals.lightsource.position,
  362.                    &color,
  363.                    &mapvals.lightsource.color,
  364.                    mapvals.lightsource.type);
  365.  
  366.           gck_rgba_clamp (&color);
  367.  
  368.           gimp_vector3_sub (&normal, &spos2, &mapvals.position);
  369.           gimp_vector3_normalize (&normal);
  370.           sphere_to_image (&normal, &vx, &vy); 
  371.           color2 = get_image_color (vx, vy, &inside);
  372.  
  373.           /* Make the normal point inwards */
  374.           /* ============================= */
  375.  
  376.           gimp_vector3_mul (&normal, -1.0);
  377.  
  378.           color2=phong_shade (&spos2,
  379.                   &mapvals.viewpoint,
  380.                   &normal,
  381.                   &mapvals.lightsource.position,
  382.                   &color2,
  383.                   &mapvals.lightsource.color,
  384.                   mapvals.lightsource.type);
  385.  
  386.           gck_rgba_clamp (&color2);
  387.  
  388.           if (mapvals.transparent_background == FALSE && color2.a < 1.0)
  389.             {
  390.               color2.r = (color2.r*color2.a)+(background.r*(1.0-color2.a));
  391.               color2.g = (color2.g*color2.a)+(background.g*(1.0-color2.a));
  392.               color2.b = (color2.b*color2.a)+(background.b*(1.0-color2.a));
  393.               color2.a = 1.0;
  394.             }
  395.  
  396.           /* Compute a mix of the first and second colors */
  397.           /* ============================================ */
  398.  
  399.           color.r = color.r*color.a+(1.0-color.a)*color2.r;
  400.           color.g = color.g*color.a+(1.0-color.a)*color2.g;
  401.           color.b = color.b*color.a+(1.0-color.a)*color2.b; 
  402.           color.a = color.a+color2.a;
  403.  
  404.           gck_rgba_clamp (&color);
  405.         }
  406.       else if (color.a!=0.0 &&
  407.            inside==TRUE &&
  408.            mapvals.lightsource.type!=NO_LIGHT)
  409.         {
  410.           /* Compute shading at this point */
  411.           /* ============================= */
  412.  
  413.           color = phong_shade (&spos1,
  414.                    &mapvals.viewpoint,
  415.                    &normal,
  416.                    &mapvals.lightsource.position,
  417.                    &color,
  418.                    &mapvals.lightsource.color,
  419.                    mapvals.lightsource.type);
  420.  
  421.           gck_rgba_clamp (&color);
  422.         }
  423.     }
  424.   
  425.   if (color.a == 0.0)
  426.     color = background;
  427.   
  428.   return color;
  429. }
  430.  
  431. /***************************************************/
  432. /* Transform the corners of the bounding box to 2D */
  433. /***************************************************/
  434.  
  435. void
  436. compute_bounding_box (void)
  437. {
  438.   GimpVector3 p1,p2;
  439.   gdouble t;
  440.   GimpVector3 dir;
  441.  
  442.   p1=mapvals.position;
  443.   p1.x-=(mapvals.radius+0.01);
  444.   p1.y-=(mapvals.radius+0.01);
  445.  
  446.   p2=mapvals.position;
  447.   p2.x+=(mapvals.radius+0.01);
  448.   p2.y+=(mapvals.radius+0.01);
  449.  
  450.   gimp_vector3_sub (&dir, &p1, &mapvals.viewpoint);
  451.   gimp_vector3_normalize (&dir);
  452.  
  453.   if (dir.z!=0.0)
  454.     {
  455.       t=(-1.0*mapvals.viewpoint.z)/dir.z;
  456.       p1.x=(mapvals.viewpoint.x+t*dir.x);
  457.       p1.y=(mapvals.viewpoint.y+t*dir.y);
  458.     }
  459.  
  460.   gimp_vector3_sub (&dir, &p2, &mapvals.viewpoint);
  461.   gimp_vector3_normalize (&dir);
  462.  
  463.   if (dir.z!=0.0)
  464.     {
  465.       t=(-1.0*mapvals.viewpoint.z)/dir.z;
  466.       p2.x=(mapvals.viewpoint.x+t*dir.x);
  467.       p2.y=(mapvals.viewpoint.y+t*dir.y);
  468.     }
  469.  
  470.   bx1=p1.x;
  471.   by1=p1.y;
  472.   bx2=p2.x;
  473.   by2=p2.y;
  474. }
  475.  
  476. /* These two were taken from the Mesa source. Mesa is written   */
  477. /* and is (C) by Brian Paul. vecmulmat() performs a post-mul by */
  478. /* a 4x4 matrix to a 1x4(3) vector. rotmat() creates a matrix   */
  479. /* that by post-mul will rotate a 1x4(3) vector the given angle */
  480. /* about the given axis.                                        */
  481. /* ============================================================ */
  482.  
  483. void
  484. vecmulmat (GimpVector3 *u,
  485.        GimpVector3 *v,
  486.        gfloat       m[16])
  487. {
  488.    gfloat v0=v->x, v1=v->y, v2=v->z;
  489. #define M(row,col)  m[col*4+row]
  490.    u->x = v0 * M(0,0) + v1 * M(1,0) + v2 * M(2,0) + M(3,0);
  491.    u->y = v0 * M(0,1) + v1 * M(1,1) + v2 * M(2,1) + M(3,1);
  492.    u->z = v0 * M(0,2) + v1 * M(1,2) + v2 * M(2,2) + M(3,2);
  493. #undef M
  494. }
  495.  
  496. void
  497. rotatemat (gfloat       angle,
  498.        GimpVector3 *v,
  499.        gfloat       m[16])
  500. {
  501.    /* This function contributed by Erich Boleyn (erich@uruk.org) */
  502.    gfloat mag, s, c;
  503.    gfloat xx, yy, zz, xy, yz, zx, xs, ys, zs, one_c;
  504.    gfloat IdentityMat[16];
  505.    gint cnt;
  506.  
  507.    s = sin( angle * (G_PI / 180.0) );
  508.    c = cos( angle * (G_PI / 180.0) );
  509.  
  510.    mag = sqrt( v->x*v->x + v->y*v->y + v->z*v->z );
  511.  
  512.    if (mag == 0.0) {
  513.       /* generate an identity matrix and return */
  514.  
  515.       for (cnt=0;cnt<16;cnt++)
  516.         IdentityMat[cnt]=0.0;
  517.     
  518.       IdentityMat[0] = 1.0;
  519.       IdentityMat[5] = 1.0;
  520.       IdentityMat[10] = 1.0;
  521.       IdentityMat[15] = 1.0;
  522.  
  523.       memcpy(m, IdentityMat, sizeof(gfloat)*16);
  524.       return;
  525.    }
  526.  
  527.    v->x /= mag;
  528.    v->y /= mag;
  529.    v->z /= mag;
  530.  
  531. #define M(row,col)  m[col*4+row]
  532.  
  533.    xx = v->x * v->x;
  534.    yy = v->y * v->y;
  535.    zz = v->z * v->z;
  536.    xy = v->x * v->y;
  537.    yz = v->y * v->z;
  538.    zx = v->z * v->x;
  539.    xs = v->x * s;
  540.    ys = v->y * s;
  541.    zs = v->z * s;
  542.    one_c = 1.0F - c;
  543.  
  544.    M(0,0) = (one_c * xx) + c;
  545.    M(0,1) = (one_c * xy) - zs;
  546.    M(0,2) = (one_c * zx) + ys;
  547.    M(0,3) = 0.0F;
  548.  
  549.    M(1,0) = (one_c * xy) + zs;
  550.    M(1,1) = (one_c * yy) + c;
  551.    M(1,2) = (one_c * yz) - xs;
  552.    M(1,3) = 0.0F;
  553.  
  554.    M(2,0) = (one_c * zx) - ys;
  555.    M(2,1) = (one_c * yz) + xs;
  556.    M(2,2) = (one_c * zz) + c;
  557.    M(2,3) = 0.0F;
  558.  
  559.    M(3,0) = 0.0F;
  560.    M(3,1) = 0.0F;
  561.    M(3,2) = 0.0F;
  562.    M(3,3) = 1.0F;
  563.  
  564. #undef M
  565. }
  566.  
  567. /* Transpose the matrix m. If m is orthogonal (like a rotation matrix), */
  568. /* this is equal to the inverse of the matrix.                          */
  569. /* ==================================================================== */
  570.  
  571. void
  572. transpose_mat (gfloat m[16])
  573. {
  574.   gint i,j;
  575.   gfloat t;
  576.  
  577.   for (i=0;i<4;i++)
  578.     {
  579.       for (j=0;j<i;j++)
  580.         {
  581.           t = m[j*4+i];
  582.           m[j*4+i] = m[i*4+j];
  583.           m[i*4+j] = t;
  584.         }
  585.     }
  586. }
  587.  
  588. /* Compute the matrix product c=a*b */
  589. /* ================================ */
  590.  
  591. void
  592. matmul (gfloat a[16],
  593.     gfloat b[16],
  594.     gfloat c[16])
  595. {
  596.   gint   i,j,k;
  597.   gfloat value;
  598.   
  599. #define A(row,col)  a[col*4+row]
  600. #define B(row,col)  b[col*4+row]
  601. #define C(row,col)  c[col*4+row]
  602.  
  603.   for (i=0;i<4;i++)
  604.     {
  605.       for (j=0;j<4;j++)
  606.         {
  607.       value = 0.0;
  608.       for (k=0;k<4;k++)
  609.         value += A(i,k)*B(k,j);
  610.           
  611.           C(i,j) = value;
  612.         }
  613.     }
  614.  
  615. #undef A
  616. #undef B
  617. #undef C
  618. }
  619.  
  620. void
  621. ident_mat (gfloat m[16])
  622. {
  623.   gint   i,j;
  624.   
  625. #define M(row,col)  m[col*4+row]
  626.   
  627.   for (i=0;i<4;i++)
  628.     {
  629.       for (j=0;j<4;j++)
  630.         {
  631.           if (i==j)
  632.             M(i,j) = 1.0;
  633.           else
  634.             M(i,j) = 0.0;
  635.         }
  636.     }  
  637.  
  638. #undef M
  639. }
  640.  
  641. static gboolean
  642. intersect_rect (gdouble            u,
  643.         gdouble            v,
  644.         gdouble            w,
  645.         GimpVector3        viewp,
  646.         GimpVector3        dir,
  647.         FaceIntersectInfo *face_info)
  648. {
  649.   gboolean result = FALSE;
  650.   gdouble u2,v2;
  651.   
  652.   if (dir.z!=0.0)
  653.     {
  654.       u2 = u / 2.0;
  655.       v2 = v / 2.0;
  656.  
  657.       face_info->t = (w-viewp.z)/dir.z;
  658.       face_info->s.x = viewp.x + face_info->t*dir.x;
  659.       face_info->s.y = viewp.y + face_info->t*dir.y;
  660.       face_info->s.z = w;
  661.  
  662.       if (face_info->s.x>=-u2 && face_info->s.x<=u2 && 
  663.           face_info->s.y>=-v2 && face_info->s.y<=v2)
  664.         {
  665.           face_info->u = (face_info->s.x + u2)/u;
  666.           face_info->v = (face_info->s.y + v2)/v;
  667.           result = TRUE;
  668.         }
  669.     }
  670.   
  671.   return(result);
  672. }
  673.  
  674. static gboolean
  675. intersect_box (GimpVector3        scale,
  676.            GimpVector3        viewp,
  677.            GimpVector3        dir,
  678.            FaceIntersectInfo *face_intersect)
  679. {
  680.   GimpVector3 v,d,tmp,axis[3];
  681.   FaceIntersectInfo face_tmp;
  682.   gboolean result = FALSE;
  683.   gfloat m[16];
  684.   gint i = 0;
  685.  
  686.   gimp_vector3_set(&axis[0], 1.0,0.0,0.0);
  687.   gimp_vector3_set(&axis[1], 0.0,1.0,0.0);
  688.   gimp_vector3_set(&axis[2], 0.0,0.0,1.0);
  689.  
  690.   /* Front side */
  691.   /* ========== */
  692.  
  693.   if (intersect_rect(scale.x,scale.y,scale.z/2.0,viewp,dir,&face_intersect[i])==TRUE)
  694.     {
  695.       face_intersect[i].face = 0;
  696.       gimp_vector3_set(&face_intersect[i++].n, 0.0,0.0,1.0);
  697.       result = TRUE;
  698.     }
  699.  
  700.   /* Back side */
  701.   /* ========= */
  702.  
  703.   if (intersect_rect(scale.x,scale.y,-scale.z/2.0,viewp,dir,&face_intersect[i])==TRUE)
  704.     {
  705.       face_intersect[i].face = 1;
  706.       face_intersect[i].u = 1.0 - face_intersect[i].u;
  707.       gimp_vector3_set(&face_intersect[i++].n, 0.0,0.0,-1.0);
  708.       result = TRUE;
  709.     }
  710.  
  711.   /* Check if we've found the two possible intersection points */
  712.   /* ========================================================= */
  713.   
  714.   if (i<2)
  715.     {
  716.       /* Top: Rotate viewpoint and direction into rectangle's local coordinate system */
  717.       /* ============================================================================ */
  718.   
  719.       rotatemat(90, &axis[0], m);
  720.       vecmulmat(&v,&viewp,m);
  721.       vecmulmat(&d,&dir,m);
  722.     
  723.       if (intersect_rect(scale.x,scale.z,scale.y/2.0,v,d,&face_intersect[i])==TRUE)
  724.         {
  725.           face_intersect[i].face = 2;
  726.  
  727.           transpose_mat(m);
  728.           vecmulmat(&tmp, &face_intersect[i].s, m);
  729.           face_intersect[i].s = tmp;
  730.           
  731.           gimp_vector3_set(&face_intersect[i++].n, 0.0,-1.0,0.0);
  732.           result = TRUE;
  733.         } 
  734.     }
  735.  
  736.   /* Check if we've found the two possible intersection points */
  737.   /* ========================================================= */
  738.   
  739.   if (i<2)
  740.     {
  741.       /* Bottom: Rotate viewpoint and direction into rectangle's local coordinate system */
  742.       /* =============================================================================== */
  743.     
  744.       rotatemat(90, &axis[0], m);  
  745.       vecmulmat(&v,&viewp,m);
  746.       vecmulmat(&d,&dir,m);
  747.     
  748.       if (intersect_rect(scale.x,scale.z,-scale.y/2.0,v,d,&face_intersect[i])==TRUE)
  749.         {
  750.           face_intersect[i].face = 3;
  751.  
  752.           transpose_mat(m);
  753.  
  754.           vecmulmat(&tmp, &face_intersect[i].s, m);
  755.           face_intersect[i].s = tmp;
  756.  
  757.           face_intersect[i].v = 1.0 - face_intersect[i].v;
  758.  
  759.           gimp_vector3_set(&face_intersect[i++].n, 0.0,1.0,0.0);
  760.  
  761.           result = TRUE;
  762.         }
  763.     }
  764.  
  765.   /* Check if we've found the two possible intersection points */
  766.   /* ========================================================= */
  767.   
  768.   if (i<2)
  769.     {
  770.       /* Left side: Rotate viewpoint and direction into rectangle's local coordinate system */
  771.       /* ================================================================================== */
  772.     
  773.       rotatemat(90, &axis[1], m);  
  774.       vecmulmat(&v,&viewp,m);
  775.       vecmulmat(&d,&dir,m);
  776.  
  777.       if (intersect_rect(scale.z,scale.y,scale.x/2.0,v,d,&face_intersect[i])==TRUE)
  778.         {
  779.           face_intersect[i].face = 4;
  780.     
  781.           transpose_mat(m);
  782.           vecmulmat(&tmp, &face_intersect[i].s, m);
  783.           face_intersect[i].s = tmp;
  784.  
  785.           gimp_vector3_set(&face_intersect[i++].n, 1.0,0.0,0.0);
  786.           result = TRUE;
  787.         }
  788.     }
  789.  
  790.   /* Check if we've found the two possible intersection points */
  791.   /* ========================================================= */
  792.   
  793.   if (i<2)
  794.     {
  795.       /* Right side: Rotate viewpoint and direction into rectangle's local coordinate system */
  796.       /* =================================================================================== */
  797.   
  798.       rotatemat(90, &axis[1], m);  
  799.       vecmulmat(&v,&viewp,m);
  800.       vecmulmat(&d,&dir,m);
  801.     
  802.       if (intersect_rect(scale.z,scale.y,-scale.x/2.0,v,d,&face_intersect[i])==TRUE)
  803.         {
  804.           face_intersect[i].face = 5;
  805.  
  806.           transpose_mat(m);
  807.           vecmulmat(&tmp, &face_intersect[i].s, m);
  808.  
  809.           face_intersect[i].u = 1.0 - face_intersect[i].u;
  810.  
  811.           gimp_vector3_set(&face_intersect[i++].n, -1.0,0.0,0.0);
  812.           result = TRUE;
  813.         }
  814.     }
  815.  
  816.   /* Sort intersection points */
  817.   /* ======================== */
  818.   
  819.   if (face_intersect[0].t>face_intersect[1].t)
  820.     {
  821.       face_tmp = face_intersect[0];
  822.       face_intersect[0] = face_intersect[1];
  823.       face_intersect[1] = face_tmp;
  824.     }
  825.  
  826.   return(result);  
  827. }
  828.  
  829. GckRGB
  830. get_ray_color_box (GimpVector3 *pos)
  831. {
  832.   GimpVector3 lvp,ldir,vp,p,dir,ns,nn;
  833.   GckRGB color, color2;
  834.   gfloat m[16];
  835.   gint i;
  836.   FaceIntersectInfo face_intersect[2];
  837.   
  838.   color=background;
  839.   vp = mapvals.viewpoint;
  840.   p = *pos;
  841.  
  842.   /* Translate viewpoint so that the box has its origin */
  843.   /* at its lower left corner.                          */
  844.   /* ================================================== */
  845.  
  846.   vp.x = vp.x - mapvals.position.x;
  847.   vp.y = vp.y - mapvals.position.y;
  848.   vp.z = vp.z - mapvals.position.z;
  849.  
  850.   p.x = p.x - mapvals.position.x;
  851.   p.y = p.y - mapvals.position.y;
  852.   p.z = p.z - mapvals.position.z;
  853.  
  854.   /* Compute direction */
  855.   /* ================= */
  856.   
  857.   gimp_vector3_sub(&dir,&p,&vp);
  858.   gimp_vector3_normalize(&dir);
  859.  
  860.   /* Compute inverse of rotation matrix and apply it to   */
  861.   /* the viewpoint and direction. This transforms the     */
  862.   /* observer into the local coordinate system of the box */
  863.   /* ==================================================== */
  864.  
  865.   memcpy(m,rotmat,sizeof(gfloat)*16);
  866.   
  867.   transpose_mat(m);
  868.   
  869.   vecmulmat(&lvp,&vp,m);
  870.   vecmulmat(&ldir,&dir,m);
  871.  
  872.   /* Ok. Now the observer is in the space where the box is located */
  873.   /* with its lower left corner at the origin and its axis aligned */
  874.   /* to the cartesian basis. Check if the transformed ray hits it. */
  875.   /* ============================================================= */
  876.  
  877.   face_intersect[0].t = 1000000.0;
  878.   face_intersect[1].t = 1000000.0;
  879.   
  880.   if (intersect_box(mapvals.scale,lvp,ldir,face_intersect)==TRUE)
  881.     {
  882.       /* We've hit the box. Transform the hit points and */
  883.       /* normals back into the world coordinate system   */
  884.       /* =============================================== */
  885.       
  886.       for (i=0;i<2;i++)
  887.         {
  888.           vecmulmat(&ns,&face_intersect[i].s,rotmat);
  889.           vecmulmat(&nn,&face_intersect[i].n,rotmat);
  890.  
  891.           ns.x = ns.x + mapvals.position.x;
  892.           ns.y = ns.y + mapvals.position.y;
  893.           ns.z = ns.z + mapvals.position.z;
  894.           
  895.           face_intersect[i].s = ns;
  896.           face_intersect[i].n = nn;
  897.         }
  898.  
  899.       color = get_box_image_color(face_intersect[0].face,
  900.                                   face_intersect[0].u,face_intersect[0].v);
  901.  
  902.       /* Check for total transparency... */
  903.       /* =============================== */
  904.  
  905.       if (color.a<1.0)
  906.         {
  907.           /* Hey, we can see  through here!      */
  908.           /* Lets see what's on the other side.. */
  909.           /* =================================== */
  910.           
  911.           color = phong_shade (&face_intersect[0].s,
  912.                    &mapvals.viewpoint,
  913.                    &face_intersect[0].n,
  914.                    &mapvals.lightsource.position,
  915.                    &color,
  916.                    &mapvals.lightsource.color,
  917.                    mapvals.lightsource.type);      
  918.     
  919.           gck_rgba_clamp(&color);
  920.  
  921.           color2 = get_box_image_color(face_intersect[1].face,
  922.                                        face_intersect[1].u,face_intersect[1].v);
  923.  
  924.           /* Make the normal point inwards */
  925.           /* ============================= */
  926.  
  927.           gimp_vector3_mul(&face_intersect[1].n,-1.0);
  928.  
  929.           color2 = phong_shade (&face_intersect[1].s,
  930.                 &mapvals.viewpoint,
  931.                 &face_intersect[1].n,
  932.                 &mapvals.lightsource.position,
  933.                 &color2,
  934.                 &mapvals.lightsource.color,
  935.                 mapvals.lightsource.type);      
  936.  
  937.           gck_rgba_clamp(&color2);
  938.           
  939.           if (mapvals.transparent_background==FALSE && color2.a<1.0)
  940.             {
  941.               color2.r = (color2.r*color2.a)+(background.r*(1.0-color2.a));
  942.               color2.g = (color2.g*color2.a)+(background.g*(1.0-color2.a));
  943.               color2.b = (color2.b*color2.a)+(background.b*(1.0-color2.a));
  944.               color2.a = 1.0;
  945.             }
  946.  
  947.           /* Compute a mix of the first and second colors */
  948.           /* ============================================ */
  949.           
  950.           color.r = color.r*color.a+(1.0-color.a)*color2.r;
  951.           color.g = color.g*color.a+(1.0-color.a)*color2.g;
  952.           color.b = color.b*color.a+(1.0-color.a)*color2.b; 
  953.           color.a = color.a+color2.a;
  954.  
  955.           gck_rgba_clamp(&color);
  956.         }
  957.       else if (color.a!=0.0 && mapvals.lightsource.type!=NO_LIGHT)
  958.         {
  959.            color = phong_shade (&face_intersect[0].s,
  960.                 &mapvals.viewpoint,
  961.                 &face_intersect[0].n,
  962.                 &mapvals.lightsource.position,
  963.                 &color,
  964.                 &mapvals.lightsource.color,
  965.                 mapvals.lightsource.type);      
  966.     
  967.            gck_rgba_clamp (&color);
  968.         }
  969.     }
  970.   else
  971.     {
  972.       if (mapvals.transparent_background == TRUE)
  973.         color.a = 0.0;
  974.     }
  975.  
  976.   return color;
  977. }
  978.  
  979. static gboolean
  980. intersect_circle (GimpVector3        vp,
  981.           GimpVector3        dir,
  982.           gdouble            w,
  983.           FaceIntersectInfo *face_info)
  984. {
  985.   gboolean result = FALSE;
  986.   gdouble r,d;
  987.   
  988. #define sqr(a) a*a
  989.  
  990.   if (dir.y!=0.0)
  991.     {
  992.       face_info->t = (w-vp.y)/dir.y;
  993.       face_info->s.x = vp.x + face_info->t*dir.x;
  994.       face_info->s.y = w;
  995.       face_info->s.z = vp.z + face_info->t*dir.z;
  996.  
  997.       r = sqrt(sqr(face_info->s.x) + sqr(face_info->s.z));
  998.       
  999.       if (r<=mapvals.cylinder_radius)
  1000.         {
  1001.           d = 2.0*mapvals.cylinder_radius;
  1002.           face_info->u = (face_info->s.x+mapvals.cylinder_radius)/d;
  1003.           face_info->v = (face_info->s.z+mapvals.cylinder_radius)/d;
  1004.           result = TRUE;
  1005.         }
  1006.     }
  1007.  
  1008. #undef sqr
  1009.   
  1010.   return(result);  
  1011. }
  1012.  
  1013. static gboolean
  1014. intersect_cylinder (GimpVector3        vp,
  1015.             GimpVector3        dir,
  1016.             FaceIntersectInfo *face_intersect)
  1017. {
  1018.   gdouble a,b,c,d,e,f,tmp,l;
  1019.   gboolean result = FALSE;
  1020.   gint i;
  1021.  
  1022. #define sqr(a) a*a
  1023.  
  1024.   a = sqr(dir.x) + sqr(dir.z);
  1025.   b = 2.0*(vp.x*dir.x+vp.z*dir.z);
  1026.   c = sqr(vp.x)+sqr(vp.z)-sqr(mapvals.cylinder_radius);
  1027.  
  1028.   d = sqr(b)-4.0*a*c;
  1029.   if (d>=0.0)
  1030.     {
  1031.       e = sqrt(d);
  1032.       f = 2.0*a;
  1033.     
  1034.       if (f!=0.0)
  1035.         {
  1036.           result = TRUE;      
  1037.     
  1038.           face_intersect[0].t = (-b+e)/f;
  1039.           face_intersect[1].t = (-b-e)/f;
  1040.     
  1041.           if (face_intersect[0].t>face_intersect[1].t)
  1042.             {
  1043.               tmp = face_intersect[0].t;
  1044.               face_intersect[0].t = face_intersect[1].t;
  1045.               face_intersect[1].t = tmp;
  1046.             }
  1047.           
  1048.           for (i=0;i<2;i++)
  1049.             {
  1050.               face_intersect[i].s.x = vp.x + face_intersect[i].t * dir.x;
  1051.               face_intersect[i].s.y = vp.y + face_intersect[i].t * dir.y;
  1052.               face_intersect[i].s.z = vp.z + face_intersect[i].t * dir.z;
  1053.  
  1054.               face_intersect[i].n = face_intersect[i].s;
  1055.               face_intersect[i].n.y = 0.0;
  1056.               gimp_vector3_normalize(&face_intersect[i].n);
  1057.  
  1058.               l = mapvals.cylinder_length/2.0;
  1059.  
  1060.               face_intersect[i].u = (atan2(face_intersect[i].s.x,face_intersect[i].s.z)+G_PI)/(2.0*G_PI);
  1061.               face_intersect[i].v = (face_intersect[i].s.y+l)/mapvals.cylinder_length;
  1062.  
  1063.               /* Mark hitpoint as on the cylinder hull */
  1064.               /* ===================================== */
  1065.               
  1066.               face_intersect[i].face = 0;
  1067.  
  1068.               /* Check if we're completely off the cylinder axis */
  1069.               /* =============================================== */
  1070.     
  1071.               if (face_intersect[i].s.y>l || face_intersect[i].s.y<-l)
  1072.                 {
  1073.                   /* Check if we've hit a cap */
  1074.                   /* ======================== */
  1075.     
  1076.                   if (face_intersect[i].s.y>l)
  1077.                     {
  1078.                       if (intersect_circle(vp,dir,l,&face_intersect[i])==FALSE)
  1079.                         result = FALSE;
  1080.                       else
  1081.                         {
  1082.                           face_intersect[i].face = 2;
  1083.                           face_intersect[i].v = 1 - face_intersect[i].v;
  1084.                           gimp_vector3_set(&face_intersect[i].n, 0.0, 1.0, 0.0);
  1085.                         }
  1086.                     }
  1087.                   else
  1088.                     {
  1089.                       if (intersect_circle(vp,dir,-l,&face_intersect[i])==FALSE)
  1090.                         result = FALSE;
  1091.                       else
  1092.                         {
  1093.                           face_intersect[i].face = 1;
  1094.                           gimp_vector3_set(&face_intersect[i].n, 0.0, -1.0, 0.0);
  1095.                         }
  1096.                     }
  1097.                 }
  1098.             }
  1099.         }
  1100.     }
  1101.  
  1102. #undef sqr
  1103.   
  1104.   return(result);
  1105. }
  1106.  
  1107. static GckRGB
  1108. get_cylinder_color (gint    face,
  1109.             gdouble u,
  1110.             gdouble v)
  1111. {
  1112.   GckRGB color;
  1113.   gint inside;
  1114.   
  1115.   if (face==0)
  1116.     color = get_image_color(u,v,&inside);
  1117.   else
  1118.     color = get_cylinder_image_color (face-1,u,v);
  1119.   
  1120.   return(color);
  1121. }
  1122.  
  1123. GckRGB
  1124. get_ray_color_cylinder (GimpVector3 *pos)
  1125. {
  1126.   GimpVector3 lvp,ldir,vp,p,dir,ns,nn;
  1127.   GckRGB color, color2;
  1128.   gfloat m[16];
  1129.   gint i;
  1130.   FaceIntersectInfo face_intersect[2];
  1131.   
  1132.   color=background;
  1133.   vp = mapvals.viewpoint;
  1134.   p = *pos;
  1135.  
  1136.   vp.x = vp.x - mapvals.position.x;
  1137.   vp.y = vp.y - mapvals.position.y;
  1138.   vp.z = vp.z - mapvals.position.z;
  1139.  
  1140.   p.x = p.x - mapvals.position.x;
  1141.   p.y = p.y - mapvals.position.y;
  1142.   p.z = p.z - mapvals.position.z;
  1143.  
  1144.   /* Compute direction */
  1145.   /* ================= */
  1146.   
  1147.   gimp_vector3_sub(&dir,&p,&vp);
  1148.   gimp_vector3_normalize(&dir);
  1149.  
  1150.   /* Compute inverse of rotation matrix and apply it to   */
  1151.   /* the viewpoint and direction. This transforms the     */
  1152.   /* observer into the local coordinate system of the box */
  1153.   /* ==================================================== */
  1154.  
  1155.   memcpy(m,rotmat,sizeof(gfloat)*16);
  1156.   
  1157.   transpose_mat(m);
  1158.   
  1159.   vecmulmat(&lvp,&vp,m);
  1160.   vecmulmat(&ldir,&dir,m);
  1161.  
  1162.   if (intersect_cylinder(lvp,ldir,face_intersect)==TRUE)
  1163.     {
  1164.       /* We've hit the cylinder. Transform the hit points and */
  1165.       /* normals back into the world coordinate system        */
  1166.       /* ==================================================== */
  1167.       
  1168.       for (i=0;i<2;i++)
  1169.         {
  1170.           vecmulmat(&ns,&face_intersect[i].s,rotmat);
  1171.           vecmulmat(&nn,&face_intersect[i].n,rotmat);
  1172.  
  1173.           ns.x = ns.x + mapvals.position.x;
  1174.           ns.y = ns.y + mapvals.position.y;
  1175.           ns.z = ns.z + mapvals.position.z;
  1176.           
  1177.           face_intersect[i].s = ns;
  1178.           face_intersect[i].n = nn;
  1179.         }
  1180.  
  1181.       color = get_cylinder_color(face_intersect[0].face,
  1182.                                  face_intersect[0].u,face_intersect[0].v);
  1183.  
  1184.       /* Check for total transparency... */
  1185.       /* =============================== */
  1186.  
  1187.       if (color.a<1.0)
  1188.         {
  1189.           /* Hey, we can see  through here!      */
  1190.           /* Lets see what's on the other side.. */
  1191.           /* =================================== */
  1192.           
  1193.           color = phong_shade (&face_intersect[0].s,
  1194.                    &mapvals.viewpoint,
  1195.                    &face_intersect[0].n,
  1196.                    &mapvals.lightsource.position,
  1197.                    &color,
  1198.                    &mapvals.lightsource.color,
  1199.                    mapvals.lightsource.type);      
  1200.     
  1201.           gck_rgba_clamp (&color);
  1202.  
  1203.           color2 = get_cylinder_color (face_intersect[1].face,
  1204.                        face_intersect[1].u,
  1205.                        face_intersect[1].v);
  1206.  
  1207.           /* Make the normal point inwards */
  1208.           /* ============================= */
  1209.  
  1210.           gimp_vector3_mul (&face_intersect[1].n, -1.0);
  1211.  
  1212.           color2 = phong_shade (&face_intersect[1].s,
  1213.                 &mapvals.viewpoint,
  1214.                 &face_intersect[1].n,
  1215.                 &mapvals.lightsource.position,
  1216.                 &color2,
  1217.                 &mapvals.lightsource.color,
  1218.                 mapvals.lightsource.type);
  1219.  
  1220.           gck_rgba_clamp (&color2);
  1221.           
  1222.           if (mapvals.transparent_background == FALSE && color2.a < 1.0)
  1223.             {
  1224.               color2.r = (color2.r*color2.a) + (background.r*(1.0-color2.a));
  1225.               color2.g = (color2.g*color2.a) + (background.g*(1.0-color2.a));
  1226.               color2.b = (color2.b*color2.a) + (background.b*(1.0-color2.a));
  1227.               color2.a = 1.0;
  1228.             }
  1229.  
  1230.           /* Compute a mix of the first and second colors */
  1231.           /* ============================================ */
  1232.           
  1233.           color.r = color.r*color.a + (1.0-color.a)*color2.r;
  1234.           color.g = color.g*color.a + (1.0-color.a)*color2.g;
  1235.           color.b = color.b*color.a + (1.0-color.a)*color2.b; 
  1236.           color.a = color.a+color2.a;
  1237.  
  1238.           gck_rgba_clamp (&color);
  1239.         }
  1240.       else if (color.a != 0.0 && mapvals.lightsource.type != NO_LIGHT)
  1241.         {
  1242.            color = phong_shade (&face_intersect[0].s,
  1243.                 &mapvals.viewpoint,
  1244.                 &face_intersect[0].n,
  1245.                 &mapvals.lightsource.position,
  1246.                 &color,
  1247.                 &mapvals.lightsource.color,
  1248.                 mapvals.lightsource.type);
  1249.  
  1250.            gck_rgba_clamp (&color);
  1251.         }
  1252.     }
  1253.   else
  1254.     {
  1255.       if (mapvals.transparent_background == TRUE)
  1256.         color.a = 0.0;
  1257.     }
  1258.   
  1259.   return color;
  1260. }
  1261.