home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Graphics / Graphics.zip / povsrc31.zip / bezier2.c < prev    next >
C/C++ Source or Header  |  2000-11-26  |  22KB  |  912 lines
  1. #include "frame.h"
  2. #include "vector.h"
  3. #include "povproto.h"
  4. #include "bezier.h"
  5. #include "matrices.h"
  6. #include "objects.h"
  7. #include "povray.h"
  8.  
  9. #ifdef RBezierPatch
  10.  
  11.    /* following magical constant serves to detection of multiple intesections.
  12.       according do NSK (Nishita, Sederberg,Kakimoto) set to 20% */
  13.  
  14. #define MAGICAL_CONSTANT    (1.0 -  0.2)
  15.  
  16. /*YS*/
  17. #ifndef INFINITY
  18. #define INFINITY        1.0e100
  19. #endif
  20. /*YS*/
  21. #define VERY_SMALL_EPSILON    1.0e-20
  22.  
  23. /*YS*/
  24. #ifndef MACOS
  25. /* just for debugging */
  26. #define FIXME            { printf("FIXME: file %s: line %d", \
  27.                                         __FILE__,                   \
  28.                                         __LINE__); exit(0); }   
  29.  
  30. #else
  31.     #define FIXME
  32. #endif
  33. /*YS*/
  34. static RAY        *Current_Ray;
  35. static BICUBIC_PATCH     *Current_Shape;
  36. static ISTACK         *Current_Depth_Stack;
  37.  
  38. static int bezier_compute_plane_points_distances(DBL     p[4][4],VECTOR n,DBL d, VECTOR cp[4][4]);
  39. static int rbezier_compute_plane_points_distances(DBL    p[4][4], VECTOR    n, DBL d, VECTOR cp[4][4], WEIGHTS w);
  40. static void compute_line_points_distance( DISTANCES d, DBL vec[2], DISTANCES points[2]);
  41. static int new_line0(DISTANCES p[2], DBL Line[2]);
  42. static int new_line1(DISTANCES p[2], DBL Line[2]);
  43. static int bezier_bounds(DISTANCES p[2], DBL interval[2][2]);
  44. static void split_patch_0(DISTANCES p[2], DBL where, DISTANCES first[2],DISTANCES second[2] );
  45. static void split_patch_1(DISTANCES p[2], DBL where, DISTANCES first[2],DISTANCES second[2] ) ;
  46. static int bezier_get_minmax_0(DISTANCES d,DBL *min,DBL *max);
  47. static int bezier_get_minmax_1(DISTANCES d,DBL *min,DBL *max);
  48. static void decasteljau (VECTOR  cpoints[4][4], DBL u, DBL v, VECTOR  point, VECTOR normal);
  49. static void rb_decasteljau (VECTOR  cpoints[4][4], WEIGHTS w, DBL u, DBL v, VECTOR  point, VECTOR normal);
  50. static int add_solution(DBL u,DBL v);
  51. static int find_solution( DISTANCES    p[2], DBL interval[2][2], DBL bounds[2][2]);
  52.  
  53. void find_planes (VECTOR point,VECTOR dir,VECTOR n1,DBL * d1,VECTOR n2,DBL * d2)        
  54. {
  55.   DBL tmp1, tmp2;
  56.  
  57.   if ( fabs(dir[X] < 0.5 ))
  58.     {
  59.             n1[X]=     0.0;
  60.       tmp1= n1[Y]=  dir[Z];
  61.       tmp2= n1[Z]= -dir[Y];
  62.       tmp1= sqrt(tmp1*tmp1+tmp2*tmp2);
  63.       n1[Y]/= tmp1;
  64.       n1[Z]/= tmp1;
  65.     }
  66.  else
  67.    {    
  68.       tmp1= n1[X]= -dir[Z];
  69.             n1[Y]=     0.0;
  70.       tmp2= n1[Z]=  dir[X];
  71.       tmp1= sqrt(tmp1*tmp1+tmp2*tmp2);
  72.       n1[X]/= tmp1;
  73.       n1[Z]/= tmp1;
  74.    }
  75.  
  76.   VCross(n2, dir, n1);
  77.  
  78.   VDot( *d1, point, n1);
  79.   *d1= - *d1;
  80.  
  81.   VDot( *d2, point, n2);
  82.   *d2= - *d2;
  83. }
  84.  
  85. /*
  86. this funciton computes distace between control points (cp) and plane (normal n, distance d)
  87. if there can not be intersection between simple convex hul and ray (ie all distances are greater than 0 or less than 0)
  88. returns 0
  89.  
  90. */ 
  91.  
  92. static int bezier_compute_plane_points_distances(DBL     p[4][4],VECTOR n,DBL d, VECTOR cp[4][4])
  93. {
  94.  int         i,j;
  95.  VECTOR      *cptmp;
  96.  DBL        *ptmp;
  97.  DBL        a,b,c;
  98.  DBL        min,max;
  99.  
  100.  ptmp=  *p;
  101.  cptmp=  *cp;
  102.  a= n[X]; b= n[Y]; c= n[Z];
  103.  min=   INFINITY;
  104.  max= - INFINITY;
  105.  
  106.  for (i=0; i<4; i++)
  107.    for (j=0; j<4; j++, ptmp++, cptmp++)
  108.      {
  109.        *ptmp= (*cptmp)[X]*a + (*cptmp)[Y]*b + (*cptmp)[Z]*c +d;
  110.        if (*ptmp < min) min=*ptmp;
  111.        if (*ptmp > max) max=*ptmp;
  112.      }
  113.  
  114.  if ((min>0.0)||(max<0.0))
  115.    return 0;
  116.  else
  117.    return 1;
  118. }
  119.  
  120. static int rbezier_compute_plane_points_distances(DBL    p[4][4], VECTOR    n, DBL d, VECTOR cp[4][4], WEIGHTS w)
  121. {
  122.  int         i,j;
  123.  VECTOR      *cptmp;
  124.  DBL        *ptmp;
  125.  DBL        *wptr;
  126.  
  127.  DBL        a,b,c;
  128.  DBL        min,max;
  129.  
  130.  ptmp=   *p;
  131.  cptmp=  *cp;
  132.  wptr=      *w;
  133.  
  134.  a= n[X]; b= n[Y]; c= n[Z];
  135.  min=   INFINITY;
  136.  max= - INFINITY;
  137.  
  138.  for (i=0; i<4; i++)
  139.    for (j=0; j<4; j++, ptmp++, cptmp++, wptr++)
  140.      {
  141.        *ptmp= (*wptr) * ((*cptmp)[X]*a + (*cptmp)[Y]*b + (*cptmp)[Z]*c +d );
  142.        if (*ptmp < min) min=*ptmp;
  143.        if (*ptmp > max) max=*ptmp;
  144.      }
  145.  
  146.  if ((min>0.0)||(max<0.0))
  147.    return 0;
  148.  else
  149.    return 1;
  150. }
  151.  
  152.  
  153. static void compute_line_points_distance( DISTANCES d, DBL vec[2], DISTANCES points[2])
  154. {
  155.   int i,j;
  156.   DBL a,b;
  157.   DBL     *ptmp, *xtmp, *ytmp;
  158.  
  159.   a=vec[0];
  160.   b=vec[1];  
  161.  
  162.   ptmp= *d;
  163.   xtmp= **points;
  164.   ytmp= *(points[1]);
  165.  
  166.   for (i=0; i<4;i++)
  167.     for (j=0; j<4; j++, ptmp++, xtmp++, ytmp++)
  168.       *ptmp= (*xtmp)*a + (*ytmp)*b;
  169. }
  170.  
  171. static int new_line0(DISTANCES p[2], DBL Line[2])
  172. {
  173.   DBL tmp;
  174.  
  175.   Line[0]=    p[1][3][0] - p[1][0][0] + p[1][3][3] - p[1][0][3];
  176.   Line[1]= - (p[0][3][0] - p[0][0][0] + p[0][3][3] - p[0][0][3]);
  177.  
  178.   tmp= Line[0]*Line[0]+Line[1]*Line[1];
  179.  
  180.   if (tmp < VERY_SMALL_EPSILON)
  181.     return 0;
  182.  
  183.   tmp= sqrt(tmp);
  184.  
  185.   Line[0]/= tmp;
  186.   Line[1]/= tmp;    
  187.   return 1;
  188. }
  189.  
  190. static int new_line1(DISTANCES p[2], DBL Line[2])
  191. {
  192.   DBL tmp;
  193.  
  194.   Line[0]=    p[1][0][3] -  p[1][0][0]  + p[1][3][3] - p[1][3][0];
  195.   Line[1]= -( p[0][0][3] -  p[0][0][0]  + p[0][3][3] - p[0][3][0] );
  196.  
  197.   tmp= Line[0]*Line[0]+Line[1]*Line[1];
  198.  
  199.   if (tmp < VERY_SMALL_EPSILON)
  200.     return 0;
  201.  
  202.   tmp= sqrt (tmp);
  203.   Line[0]/= tmp;
  204.   Line[1]/= tmp;
  205.   return 1;
  206. }
  207.  
  208.  
  209. /* conv. hull does not contain [0,0] */
  210. static int bezier_bounds(DISTANCES p[2], DBL interval[2][2])
  211. {
  212.   int i,j;
  213.   
  214.   DBL *itmp;
  215.   DBL *min, *max;
  216.  
  217.   int ret=0;
  218.   
  219.   for (i=0; i<2; i++)
  220.     {
  221.       itmp = &(p[i][0][0]);
  222.       min= & (interval[i][0]);
  223.       max= & (interval[i][1]);
  224.       *min=*max= *(itmp++);
  225.       for(j=15; j  ; j--, itmp++ )
  226.     {
  227.       if ( *itmp < *min )
  228.         *min= *itmp;
  229.       else
  230.         if ( *itmp > *max )
  231.           *max= *itmp;
  232.     }
  233.       ret|= (*min > EPSILON) || (*max < -EPSILON);
  234.     }
  235.   return !ret;
  236. }
  237.  
  238. static void split_patch_0(DISTANCES p[2], DBL where, DISTANCES first[2],DISTANCES second[2] )  
  239. {
  240.   int i,l;
  241.   DISTANCES tmp;
  242.  
  243.   DBL     *p_i;
  244.   DBL      *fi_i;    
  245.   DBL      *se_i;
  246.   DBL   where2;
  247.  
  248.   where2= 1.0 - where;
  249.  
  250.   for (l=0; l<2; l++)
  251.     {
  252.       p_i=   & (p[l][0][0]);
  253.       fi_i=  & (first[l][0][0]);
  254.       se_i=  & (second[l][0][0]);
  255.       se_i+= 3;
  256.  
  257.       for (i=0; 
  258.        i<4; 
  259.        i++
  260.        )
  261.     {
  262.       tmp[1][0]= where2 * (* fi_i= * p_i); p_i++; fi_i++; 
  263.                                                *fi_i= tmp[1][0]+= where * (* p_i);
  264.                            fi_i++;
  265.       tmp[1][1]= where2 * (* p_i);            p_i+= 1; tmp[1][1]+= where * (* p_i);
  266.       tmp[1][2]= where2 * (* p_i);            p_i+= 1; tmp[1][2]+= where * (*se_i = * p_i);
  267.                                           se_i--;
  268.                           p_i++;
  269.  
  270.       *fi_i= tmp[2][0]= where2 * tmp[1][0] + where * tmp[1][1];
  271.       fi_i++;
  272.                  tmp[2][1]= where2 * tmp[1][1] + where * (*se_i= tmp[1][2]);
  273.       se_i--;
  274.  
  275.       *fi_i= where2 * tmp[2][0] + where * (*se_i= tmp[2][1]);
  276.       se_i--;
  277.       *se_i= *fi_i;
  278.       fi_i++;
  279.       se_i+= 7;
  280.     }
  281.     }
  282. }
  283.  
  284. static void split_patch_1(DISTANCES p[2], DBL where, DISTANCES first[2],DISTANCES second[2] ) 
  285. {
  286.   int i,l;
  287.   DISTANCES tmp;
  288.  
  289.   DBL     *p_i;
  290.   DBL      *fi_i;    
  291.   DBL      *se_i;
  292.   DBL   where2;
  293.  
  294.   where2= 1.0 - where;
  295.  
  296.   for (l=0; l<2; l++)
  297.     {
  298.       p_i=   & (p[l][0][0]);
  299.       fi_i=  & (first[l][0][0]);
  300.       se_i=  & (second[l][0][0]);
  301.       se_i+= 3* 4;
  302.  
  303.       for (i=0; 
  304.        i<4; 
  305.        i++
  306.        )
  307.     {
  308.       tmp[1][0]= where2 * (* fi_i= * p_i); p_i+= 4; fi_i+= 4; 
  309.                                                *fi_i= tmp[1][0]+= where * (* p_i);
  310.                            fi_i+= 4;
  311.       tmp[1][1]= where2 * (* p_i);      p_i+= 4; tmp[1][1]+= where * (* p_i);
  312.       tmp[1][2]= where2 * (* p_i);      p_i+= 4; tmp[1][2]+= where * (*se_i = * p_i);
  313.                                           se_i-= 4;
  314.                           
  315.       *fi_i= tmp[2][0]= where2 * tmp[1][0] + where * tmp[1][1];
  316.       fi_i+= 4;
  317.                  tmp[2][1]= where2 * tmp[1][1] + where * (*se_i= tmp[1][2]);
  318.       se_i-= 4;
  319.  
  320.       *fi_i= where2 * tmp[2][0] + where * (*se_i= tmp[2][1]);
  321.       se_i-= 4;
  322.       *se_i= *fi_i;
  323.       fi_i-= 3*4 - 1;
  324.       p_i -= 3*4 - 1;
  325.       se_i+= 1 + 3* 4;
  326.     }
  327.     }
  328. }
  329.  
  330. static void (*split_patch[2])(DISTANCES *, DBL, DISTANCES *, DISTANCES *) =
  331.   { split_patch_0, split_patch_1 };
  332.  
  333.  
  334. static int bezier_get_minmax_0(DISTANCES d,DBL *min,DBL *max)
  335. {
  336.  DBL     *p, *n, *pp, tmp;
  337.  int     f,l;
  338.  
  339.  int     i;
  340.  
  341.  
  342.  *min =  2.0;
  343.  *max = -2.0;
  344.  
  345.  f=l= 0;
  346.  
  347.  for (i=0; i<4; i++)
  348.    {
  349.      n= 1 + (p= pp= &(d[i][0]) );
  350.      if ( *p < 0.0 ) f++;
  351.     
  352.      if ( (*p <= 0.0) ^ (*n  <= 0.0) )  /* 0,1 */
  353.        {
  354.      tmp= /* 0.0 / 3.0  + */  (1.0 / 3.0) * (*p) / ( (*p) - (*n) );
  355.  
  356.      if (tmp < *min) *min= tmp;  
  357.      if (tmp > *max) *max= tmp;
  358.        }
  359.      p++; n++;
  360.  
  361.      if ( (*p <= 0.0) ^ (*n  <= 0.0) )  /* 1,2 */
  362.        {
  363.      tmp=  1.0 / 3.0 +  (1.0 / 3.0) * (*p) / ( (*p) - (*n) );
  364.  
  365.      if (tmp < *min) *min= tmp;  
  366.      if (tmp > *max) *max= tmp;
  367.        }
  368.  
  369.      if ( (*pp <= 0.0) ^ (*n  <= 0.0) )  /* 0,2 */
  370.        {
  371.      tmp=  /* 0.0 / 3.0 + */  (2.0 / 3.0) * (*pp) / ( (*pp) - (*n) );
  372.  
  373.      if (tmp < *min) *min= tmp;  
  374.      if (tmp > *max) *max= tmp;
  375.        }
  376.  
  377.      p++; n++;
  378.  
  379.      if ( (*p <= 0.0) ^ (*n  <= 0.0) )  /* 2,3 */
  380.        {
  381.      tmp=  2.0 / 3.0 +  (1.0 / 3.0) * (*p) / ( (*p) - (*n) );
  382.  
  383.      if (tmp < *min) *min= tmp;  
  384.      if (tmp > *max) *max= tmp;
  385.        }
  386.  
  387.      if ( (*pp <= 0.0) ^ (*n  <= 0.0) )  /* 0,3 */
  388.        {
  389.      tmp= /* 0.0 / 3.0 +  (3.0 / 3.0) */  (*pp) / ( (*pp) - (*n) );
  390.  
  391.      if (tmp < *min) *min= tmp;  
  392.      if (tmp > *max) *max= tmp;
  393.        }
  394.  
  395.      pp++;
  396.  
  397.      if ( (*pp <= 0.0) ^ (*n  <= 0.0) )  /* 1,3 */
  398.        {
  399.      tmp=  1.0 / 3.0 +  (2.0 / 3.0) * (*pp) / ( (*pp) - (*n) );
  400.  
  401.      if (tmp < *min) *min= tmp;  
  402.      if (tmp > *max) *max= tmp;
  403.        }
  404.  
  405.      if ( *n  < 0.0 ) l++;
  406.    }
  407.  
  408.  if ( (f & 3) != 0) *min = 0.0;
  409.  if ( (l & 3) != 0) *max = 1.0;
  410.  
  411.  if (*min > *max) return 0;  /* no solution found */
  412.  else return 1;
  413. }
  414.  
  415. static int bezier_get_minmax_1(DISTANCES d,DBL *min,DBL *max)
  416. {
  417.  DBL     *p, *n, *pp, tmp;
  418.  int    f,l;
  419.  
  420.  int     i;
  421.  
  422.  f=l= 0;
  423.  *min =  2.0;
  424.  *max = -2.0;
  425.  
  426.  
  427.  for (i=0; i<4; i++)
  428.    {
  429.      n= 4 + (p= pp= &(d[0][i]) );
  430.      if ( *p < 0.0) f++;
  431.  
  432.      if ( (*p <= 0.0) ^ (*n  <= 0.0) )  /* 0,1 */
  433.        {
  434.      tmp= /* 0.0 / 3.0  + */  (1.0 / 3.0) * (*p) / ( (*p) - (*n) );
  435.  
  436.      if (tmp < *min) *min= tmp;  
  437.      if (tmp > *max) *max= tmp;
  438.        }
  439.      p+= 4; n+= 4;
  440.  
  441.      if ( (*p <= 0.0) ^ (*n  <= 0.0) )  /* 1,2 */
  442.        {
  443.      tmp=  1.0 / 3.0 +  (1.0 / 3.0) * (*p) / ( (*p) - (*n) );
  444.  
  445.      if (tmp < *min) *min= tmp;  
  446.      if (tmp > *max) *max= tmp;
  447.        }
  448.  
  449.      if ( (*pp <= 0.0) ^ (*n  <= 0.0) )  /* 0,2 */
  450.        {
  451.      tmp=  /* 0.0 / 3.0 */ +  (2.0 / 3.0) * (*pp) / ( (*pp) - (*n) );
  452.  
  453.      if (tmp < *min) *min= tmp;  
  454.      if (tmp > *max) *max= tmp;
  455.        }
  456.  
  457.      p+= 4; n+= 4;
  458.  
  459.      if ( (*p <= 0.0) ^ (*n  <= 0.0) )  /* 2,3 */
  460.        {
  461.      tmp=  2.0 / 3.0 +  (1.0 / 3.0) * (*p) / ( (*p) - (*n) );
  462.  
  463.      if (tmp < *min) *min= tmp;  
  464.      if (tmp > *max) *max= tmp;
  465.        }
  466.  
  467.      if ( (*pp <= 0.0) ^ (*n  <= 0.0) )  /* 0,3 */
  468.        {
  469.      tmp= /* 0.0 / 3.0 +  (3.0 / 3.0) */  (*pp) / ( (*pp) - (*n) );
  470.  
  471.      if (tmp < *min) *min= tmp;  
  472.      if (tmp > *max) *max= tmp;
  473.        }
  474.      pp+=4;
  475.  
  476.      if ( (*pp <= 0.0) ^ (*n  <= 0.0) )  /* 1,3 */
  477.        {
  478.      tmp=    1.0 / 3.0 +  (2.0 / 3.0) * (*pp) / ( (*pp) - (*n) );
  479.  
  480.      if (tmp < *min) *min= tmp;  
  481.      if (tmp > *max) *max= tmp;
  482.        }
  483.  
  484.      if ( *n < 0.0 ) l++;
  485.    }
  486.  
  487.  if ( (f & 3) != 0) *min = 0.0;
  488.  if ( (l & 3) != 0) *max = 1.0;
  489.  
  490.  if (*min > *max) return 0;  /* no solution found */
  491.  else return 1;
  492. }
  493.  
  494. static int (*bezier_get_minmax[2])(DISTANCES ,DBL *, DBL *) =
  495.   { bezier_get_minmax_0, bezier_get_minmax_1 };
  496.  
  497. /*
  498. static void
  499. solve_as_curve(p, iii, where)
  500. DISTANCES p[2];
  501. int iii;
  502. DBL where;
  503. {
  504.  FIXME;
  505. }
  506. */
  507.  
  508. static void decasteljau (VECTOR  cpoints[4][4], DBL u, DBL v, VECTOR  point, VECTOR normal)
  509. {
  510.  int      i,j;
  511.  DBL     m, n;
  512.  VECTOR  v0, v1;
  513.  DBL      *pp, *nn;
  514.  DBL      q[4];
  515.  DBL      ftmp;
  516.  
  517.  DISTANCES tmp;
  518.  
  519.  m= 1.0 - u; 
  520.  n= 1.0 - v;
  521.  
  522.  for (i=0; i<3; i++)
  523.    {
  524.      nn= 3 + (pp= &(cpoints [0][0][i]) );
  525.      for (j=0; j < 4; j++)
  526.        {
  527.      tmp[0][0]= m * (*pp) + u * (*nn);        pp+=3; nn+=3;
  528.      tmp[0][1]= m * (*pp) + u * (*nn);        pp+=3; nn+=3;
  529.      tmp[0][2]= m * (*pp) + u * (*nn);        pp+=6; nn+=6;
  530.  
  531.      tmp[1][0]= m * tmp[0][0] + u * tmp[0][1];
  532.      tmp[1][1]= m * tmp[0][1] + u * tmp[0][2];
  533.  
  534.      q[j]= m * tmp[1][0] + u * tmp[1][1];
  535.        }
  536.      
  537.      tmp[0][0]= n * q[0] + v * q[1];
  538.      tmp[0][1]= n * q[1] + v * q[2];
  539.      tmp[0][2]= n * q[2] + v * q[3];
  540.  
  541.      tmp[1][0] = n * tmp[0][0] + v * tmp[0][1];
  542.      tmp[1][1] = n * tmp[0][1] + v * tmp[0][2];
  543.  
  544.      v0[i]= tmp[1][1] - tmp[1][0];
  545.      point[i]= n * tmp[1][0] + v * tmp[1][1];
  546.    }
  547.  
  548.  for (i=0; i<3; i++)
  549.    {
  550.      nn= 3 * 4 + ( pp= &( cpoints [0][0][i] ));
  551.      for (j=0; j < 4; j++)
  552.        {
  553.      tmp[0][0]= n * (*pp) + v * (*nn);        pp+=3*4;   nn+=3*4;
  554.      tmp[0][1]= n * (*pp) + v * (*nn);        pp+=3*4;   nn+=3*4;
  555.      tmp[0][2]= n * (*pp) + v * (*nn);        pp-=2*3*4-3; nn-=2*3*4-3;
  556.  
  557.      tmp[1][0]= n * tmp[0][0] + v * tmp[0][1];
  558.      tmp[1][1]= n * tmp[0][1] + v * tmp[0][2];
  559.  
  560.      q[j]= n * tmp[1][0] + v * tmp[1][1];
  561.        }
  562.      
  563.      tmp[0][0]= m * q[0] + u * q[1];
  564.      tmp[0][1]= m * q[1] + u * q[2];
  565.      tmp[0][2]= m * q[2] + u * q[3];
  566.  
  567.      tmp[1][0] = m * tmp[0][0] + u * tmp[0][1];
  568.      tmp[1][1] = m * tmp[0][1] + u * tmp[0][2];
  569.  
  570.      v1[i]= tmp[1][1] - tmp[1][0];
  571.    }
  572.  
  573.  VCross(normal,v0,v1);
  574.  VLength(ftmp, normal);
  575.  
  576.  if (ftmp < VERY_SMALL_EPSILON)
  577.    { normal[X]=normal[Y]=0.0; normal[Z]=1.0; }
  578.  else
  579.    { normal[X]/= ftmp; normal[Y]/= ftmp; normal[Z]/= ftmp; }
  580. }
  581.  
  582. static void rb_decasteljau (VECTOR  cpoints[4][4], WEIGHTS w, DBL u, DBL v, VECTOR  point, VECTOR normal)
  583. {
  584.  int      i,j;
  585.  DBL     m, n;
  586.  
  587.  VECTOR  v0, v1;
  588.  DBL      *pp, *nn, *wptr;
  589.  
  590.  DBL      q[4];
  591.  DBL      ftmp;
  592.  DBL       wpoints[4][4][4];
  593.  DBL       w0,w1,w00;
  594.  
  595.  DBL    tmp[2][4];
  596.  
  597.  
  598.  wptr= & ( w[0][0] );
  599.  pp= & ( cpoints[0][0][0] );
  600.  nn= & ( wpoints[0][0][0] ) ;
  601.    
  602.  for (i=0; i<16; i++)
  603.    {
  604.      *nn= (*pp) * (*wptr); nn++; pp++;
  605.      *nn= (*pp) * (*wptr); nn++; pp++;
  606.      *nn= (*pp) * (*wptr); nn++; pp++;
  607.      *nn=       (*wptr); nn++; wptr++;
  608.    }
  609.  
  610.  m= 1.0 - u; 
  611.  n= 1.0 - v;
  612.  
  613.  for (i=3; i>=0; i--)
  614.    {
  615.      nn= 4 + (pp= &(wpoints [0][0][i]) );
  616.      for (j=0; j < 4; j++)
  617.        {
  618.      tmp[0][0]= m * (*pp) + u * (*nn);        pp+=4; nn+=4;
  619.      tmp[0][1]= m * (*pp) + u * (*nn);        pp+=4; nn+=4;
  620.      tmp[0][2]= m * (*pp) + u * (*nn);        pp+=8; nn+=8;
  621.  
  622.      tmp[1][0]= m * tmp[0][0] + u * tmp[0][1];
  623.      tmp[1][1]= m * tmp[0][1] + u * tmp[0][2];
  624.  
  625.      q[j]= m * tmp[1][0] + u * tmp[1][1];
  626.        }
  627.  
  628.      tmp[0][0]= n * q[0] + v * q[1];
  629.      tmp[0][1]= n * q[1] + v * q[2];
  630.      tmp[0][2]= n * q[2] + v * q[3];
  631.  
  632.      tmp[1][0] = n * tmp[0][0] + v * tmp[0][1];
  633.      tmp[1][1] = n * tmp[0][1] + v * tmp[0][2];
  634.  
  635.      if (i!=3)
  636.        {
  637.      v0[i]= tmp[1][1] / w1 - tmp[1][0] / w0; 
  638.      point[i]= (n * tmp[1][0] + v * tmp[1][1]) / w00; 
  639.        }
  640.      else
  641.        {
  642.      w00= (n *  (w0= tmp[1][0]) + v * (w1= tmp[1][1]) );
  643.        }
  644.    }
  645.  
  646.  for (i=3; i>=0; i--)
  647.    {
  648.      nn= 4 * 4 + ( pp= &( wpoints [0][0][i] ));
  649.      for (j=0; j < 4; j++)
  650.        {
  651.      tmp[0][0]= n * (*pp) + v * (*nn);        pp+=4*4;   nn+=4*4;
  652.      tmp[0][1]= n * (*pp) + v * (*nn);        pp+=4*4;   nn+=4*4;
  653.      tmp[0][2]= n * (*pp) + v * (*nn);        pp-=2*4*4-4; nn-=2*4*4-4;
  654.  
  655.      tmp[1][0]= n * tmp[0][0] + v * tmp[0][1];
  656.      tmp[1][1]= n * tmp[0][1] + v * tmp[0][2];
  657.  
  658.      q[j]= n * tmp[1][0] + v * tmp[1][1];
  659.        }
  660.      
  661.      tmp[0][0]= m * q[0] + u * q[1];
  662.      tmp[0][1]= m * q[1] + u * q[2];
  663.      tmp[0][2]= m * q[2] + u * q[3];
  664.  
  665.      tmp[1][0] = m * tmp[0][0] + u * tmp[0][1];
  666.      tmp[1][1] = m * tmp[0][1] + u * tmp[0][2];
  667.  
  668.      if (i!=3)
  669.        {
  670.      v1[i]= tmp[1][1] / w1 - tmp[1][0] / w0;
  671.        }
  672.      else
  673.        {
  674.      w0= tmp[1][0];     w1= tmp[1][1];
  675.        }
  676.    }
  677.  
  678.  VCross(normal,v0,v1);
  679.  VLength(ftmp, normal);
  680.  
  681.  if (ftmp < VERY_SMALL_EPSILON)
  682.    { normal[X]=normal[Y]=0.0; normal[Z]=1.0; }
  683.  else
  684.    { normal[X]/= ftmp; normal[Y]/= ftmp; normal[Z]/= ftmp; }
  685.  
  686. }
  687.  
  688.  
  689. static int add_solution(DBL u,DBL v)
  690. {
  691.   VECTOR point, normal;
  692.   DBL  depth;
  693. #ifdef UpdatedUVMappingPatch
  694.   UV_VECT UV;  /* MH */
  695.   double uv_point[2], tpoint[3]; /*MH */
  696. #endif
  697.   if ((u<0.0) || (u> 1.0) || (v < 0.0) || ( v > 1.0)) return 0;
  698.  
  699.   if ( Current_Shape ->  Patch_Type ==  RATIONAL_BEZIER_NSK)
  700.     rb_decasteljau( Current_Shape -> Control_Points,*Current_Shape -> Weights,
  701.             u, v, point,  normal );
  702.   else
  703.     decasteljau(  Current_Shape -> Control_Points, u, v, point,  normal );
  704.  
  705.   if (fabs(Current_Ray->Direction[X]) > 0.5 ) 
  706.     depth= (point[X]- Current_Ray->Initial[X]) / Current_Ray->Direction[X];
  707.   else
  708.     if (fabs(Current_Ray->Direction[Y]) > 0.5 ) 
  709.       depth= (point[Y]- Current_Ray->Initial[Y]) / Current_Ray->Direction[Y];
  710.     else
  711.       depth= (point[Z]- Current_Ray->Initial[Z]) / Current_Ray->Direction[Z];
  712.  
  713.   if  ( depth < 4.0 * Current_Shape-> accuracy) 
  714.     return 0; 
  715.  
  716. #ifndef UpdatedUVMappingPatch
  717.   push_normal_entry(depth, point, normal, (OBJECT *)Current_Shape, Current_Depth_Stack);
  718. #else
  719. /* MH */
  720.   /* transform current point from uv space to texture space */
  721.    uv_point[0] = u; 
  722.    uv_point[1] = v; 
  723.    MTransUVPoint(uv_point, Current_Shape->Mapping, tpoint); 
  724.    
  725.    UV[U] = tpoint[0];
  726.    UV[V] = tpoint[1];
  727.  
  728.    push_normal_uv_entry(depth, point, normal, UV, (OBJECT *)Current_Shape, Current_Depth_Stack);
  729. /* MH */
  730. #endif
  731.   return 1;
  732. }
  733.  
  734.  
  735. static int find_solution( DISTANCES    p[2], 
  736.     DBL        interval[2][2],/* attention, this is strange  interval,  first is begin, second is length */
  737.      DBL        bounds[2][2] /* this is bounding box in trans.  coordinates */
  738.      )
  739. {
  740.   int         iii;
  741.   DBL     Line[2];       
  742.   DBL        tmp;
  743.   
  744.   int        in_half=0;
  745.  
  746.   DISTANCES    d;
  747.   DBL        min,max;
  748.  
  749.  
  750.  
  751.   if ( (bounds[1][1]-bounds[1][0] < Current_Shape->accuracy)
  752.        && (bounds[0][1]-bounds[0][0] < Current_Shape->accuracy))
  753.     {   /* end of iteration ...*/
  754.       return add_solution( interval[0][0]+ interval[0][1]/2.0, interval[1][0] + interval[1][1]/2.0 );
  755.       FIXME  /* soupni sem presnejsi resitko ...!!! */
  756.     }
  757.  
  758.   if ( interval[0][1] > interval[1][1] )
  759.     {
  760.       iii=0;
  761.       if ( new_line0(p, Line) == 0)
  762.     if ( new_line1(p, Line) == 0)
  763.       in_half=1;   /* degenerated patch -=> split in half */
  764.     else
  765.       {
  766.         tmp= -Line[0];
  767.         Line[0]= Line[1];
  768.         Line[1]= tmp;
  769.       }
  770.     }
  771.   else
  772.     {
  773.       iii=1;
  774.       if ( new_line1(p, Line) == 0)
  775.         if ( new_line0(p, Line) == 0)
  776.           in_half=1;   /* degenerated patch -=> split in half */
  777.         else
  778.           {
  779.             tmp= -Line[0];
  780.             Line[0]= Line[1];
  781.             Line[1]= tmp;
  782.           }
  783.     }
  784.  
  785.   /* FIXME */
  786.   /* in_half=1;  */
  787.  
  788.   if (in_half == 0)
  789.     {
  790.       compute_line_points_distance( d, Line, p );
  791.       if (bezier_get_minmax[iii](d, &min, &max)==0) return 0;
  792.       if ( (tmp= max-min) > MAGICAL_CONSTANT )
  793.     in_half= 1;
  794.       else
  795.     {
  796.  
  797.       /*
  798.       if (tmp < EPSILON)  
  799.         solve_as_surve(p,iii,min);
  800.         FIXME
  801.       */
  802.  
  803.       {
  804.         DISTANCES ptmp1[2],ptmp2[2];
  805.  
  806.         /* is seems that this "speed-up" is "slow-down" :-)
  807.         if (min == 0.0)
  808.           {
  809.         split_patch[iii]( p, max, ptmp1, ptmp2);
  810.         if (bezier_bounds(ptmp1, bounds ) == 0) return 0;  
  811.         interval[iii][1] = interval[iii][1]*max;
  812.         return find_solution(ptmp1, interval, bounds);        
  813.           }
  814.  
  815.         if (max == 1.0)
  816.           {
  817.         split_patch[iii]( p, min, ptmp1, ptmp2);
  818.         if (bezier_bounds(ptmp2, bounds ) == 0) return 0;  
  819.         interval[iii][0]+= interval[iii][1]*min;
  820.         interval[iii][1] = interval[iii][1]*tmp;
  821.         return find_solution(ptmp2, interval, bounds);        
  822.           }
  823.         */
  824.         
  825.         if (min < 0.5) 
  826.           {
  827.         split_patch[iii]( p,min, ptmp1, ptmp2);
  828.         split_patch[iii]( ptmp2, (max-min) / (1.0 - min), p, ptmp1);
  829.           }
  830.         else 
  831.           {
  832.         split_patch[iii]( p,max,ptmp1, ptmp2);
  833.         split_patch[iii]( ptmp1,min/max, ptmp2, p);
  834.           }
  835.       }
  836.  
  837.       if (bezier_bounds(p, bounds ) == 0) return 0;  
  838.      
  839.       interval[iii][0]+= interval[iii][1]*min;
  840.       interval[iii][1] = interval[iii][1]*tmp;
  841.       
  842.       return find_solution(p, interval, bounds);
  843.     }
  844.     }
  845.  
  846.   if (in_half == 1)
  847.     {
  848.       DISTANCES ld[2], rd[2];
  849.       DBL    ni[2][2];
  850.       int    ret=0;
  851.  
  852.       split_patch[iii](p, 0.5, ld, rd);
  853.  
  854.       tmp= interval[iii][1]/= 2.0;
  855.       memcpy( &(ni[0][0]), &(interval[0][0]), 4 * sizeof(DBL));
  856.       ni[iii][0]+= tmp;
  857.  
  858.       if (bezier_bounds(ld,bounds) != 0) 
  859.     ret|= find_solution(ld, interval, bounds);
  860.  
  861.       if (bezier_bounds(rd, bounds) != 0)
  862.     ret|= find_solution(rd, ni, bounds);
  863.       
  864.       return ret;
  865.     }
  866.  
  867.   printf("Bezier 2: This should never ever happend :-)\n");
  868.   return 0;
  869. }
  870.  
  871. int intersect_bicubic_patch_nsk (RAY * Ray,BICUBIC_PATCH *Shape,ISTACK * Depth_Stack)
  872. {
  873.  VECTOR n0,n1;           
  874.  DBL    d0,d1;  /* two perpendicular planes. intersection( pl1, pl2 ) == Ray */
  875.  DISTANCES      p[2];
  876.  DBL        interval[2][2];
  877.  DBL        sti[2][2];
  878.  
  879.  find_planes( Ray->Initial, Ray->Direction, n0, &d0, n1, &d1 );
  880.  
  881.  if ( Shape -> Patch_Type == BEZIER_NSK )
  882.    {
  883.      if (bezier_compute_plane_points_distances(p[0], n0, d0, Shape->Control_Points) == 0)
  884.        return 0;
  885.      if (bezier_compute_plane_points_distances(p[1], n1, d1, Shape->Control_Points) == 0)
  886.        return 0;
  887.    }
  888.  else
  889.    {
  890.      if (rbezier_compute_plane_points_distances(p[0], n0, d0, Shape->Control_Points, *Shape-> Weights) == 0)
  891.        return 0;
  892.      if (rbezier_compute_plane_points_distances(p[1], n1, d1, Shape->Control_Points, *Shape-> Weights) == 0)
  893.        return 0;
  894.    }
  895.  
  896.  Current_Ray=Ray;
  897.  Current_Shape=Shape;
  898.  Current_Depth_Stack= Depth_Stack;
  899.  
  900.  interval[0][0]=interval[1][0]= 0.0;
  901.  interval[0][1]=interval[1][1]= 1.0;
  902.  
  903.  sti[1][1]= sti[0][1]= INFINITY; 
  904.  sti[1][0]= sti[0][0]=-INFINITY;
  905.  
  906.  return find_solution( p, interval, sti );
  907. }
  908. #else
  909.     static char dummy[2];
  910. #endif
  911.  
  912.