home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Graphics / Graphics.zip / povsrc31.zip / rbezier.c < prev    next >
C/C++ Source or Header  |  2000-11-26  |  40KB  |  1,671 lines
  1.  
  2. /* 
  3.   overit rychlost:
  4.    zda add_solution pocita spravne nejdrive ty s mensim orderem
  5. */
  6.  
  7. #include "frame.h"
  8. #include "vector.h"
  9. #include "povproto.h"
  10. #include "bezier.h"
  11. #include "matrices.h"
  12. #include "objects.h"
  13. #include "povray.h"
  14.  
  15. #include "bezier.h"
  16. #include "rbezier.h"
  17.  
  18. #ifdef RBezierPatch
  19.  
  20.  
  21.            /* Yet another epsilon */    
  22. #define   YAEP            (1.0e-8)
  23. #define   MAGICAL_CONSTANT  (1.0 - 0.2)
  24. /*YS*/
  25. #ifndef INFINITY
  26. #define INFINITY        (1.0e100)
  27. #endif
  28. /*YS*/
  29. #define   VERY_SMALL_EPSILON (1e-20)
  30.  
  31. /*YS*/
  32. #ifndef MACOS
  33. /* just for debugging */
  34. #define FIXME            { printf("FIXME: file %s: line %d", \
  35.                                         __FILE__,                   \
  36.                                         __LINE__); exit(0); }   
  37.  
  38. #else
  39.     #define FIXME
  40. #endif
  41. /*YS*/
  42.  
  43. static void RBezier_Patch_Normal (VECTOR, OBJECT *, INTERSECTION *);
  44. static int Inside_RBezier_Patch (VECTOR, OBJECT *);
  45. static void Translate_RBezier_Patch (OBJECT *, VECTOR, TRANSFORM *);
  46. static void Scale_RBezier_Patch (OBJECT *, VECTOR, TRANSFORM *);
  47. static void Rotate_RBezier_Patch  (OBJECT *, VECTOR, TRANSFORM *);
  48. static void Transform_RBezier_Patch (OBJECT *, TRANSFORM *);
  49. static void Invert_RBezier_Patch (OBJECT *);
  50. static void *Copy_RBezier_Patch (OBJECT *);
  51. static void Destroy_RBezier_Patch (OBJECT *);
  52. void Compute_RBezier_Patch_BBox (RBEZIER_PATCH *);
  53. static int All_RBezier_Patch_Intersections (OBJECT *, RAY *, ISTACK *);
  54. static int rbezier_check_trims (DBL,DBL);
  55.  
  56. static RAY        *Current_Ray;
  57. static RBEZIER_PATCH     *Current_Shape;
  58. static ISTACK         *Current_Depth_Stack;
  59. static TRIM_SHAPE    *Current_Trim;
  60.  
  61. static int        cu_order, cv_order;
  62. static int        cu_order_1, cv_order_1;
  63. static int        c4_u_order;
  64.  
  65. static METHODS RBezier_Patch_Methods =
  66. {
  67.   All_RBezier_Patch_Intersections,
  68.   Inside_RBezier_Patch, RBezier_Patch_Normal, Default_UVCoord,
  69.   Copy_RBezier_Patch,
  70.   Translate_RBezier_Patch, Rotate_RBezier_Patch,
  71.   Scale_RBezier_Patch, Transform_RBezier_Patch, Invert_RBezier_Patch,
  72.   Destroy_RBezier_Patch
  73. };
  74.  
  75. /*******************************************************
  76.  
  77.    (r|)bezier_compute_plane_points_distances functions
  78.  
  79. ********************************************************/
  80.  
  81. static int bezier_plane_points_distances(DBL         p[4][4],VECTOR n,DBL d,RBEZIER_PATCH *Patch)
  82. {
  83.   int        i, j;
  84.   DBL        *x,*y,*z;
  85.   DBL        *ptmp;
  86.   DBL        a,b,c;
  87.   DBL        min, max;
  88.  
  89.   ptmp=  *p;
  90.   x= Patch->pts[X];
  91.   y= Patch->pts[Y];
  92.   z= Patch->pts[Z];
  93.  
  94.   a= n[X]; b= n[Y]; c=n[Z];
  95.   min=     INFINITY;
  96.   max= -INFINITY;
  97.  
  98.   for ( i= cv_order; i; i--, ptmp+= c4_u_order, 
  99.                          x+=c4_u_order, y+=c4_u_order, z+=c4_u_order)
  100.     for ( j= cu_order; j; j--, ptmp++,  x++, y++, z++)
  101.       {
  102.     *ptmp = (*x)*a + (*y)*b + (*z)*c + d;
  103.     if (*ptmp < min) min= *ptmp;
  104.     if (*ptmp > max) max= *ptmp;
  105.       }    
  106.  
  107.   if ((min>0.0)||(max<0.0))
  108.     return 0;
  109.   else 
  110.     return 1;
  111. }
  112.  
  113. static int rbezier_plane_points_distances(DBL    p[4][4],VECTOR n,DBL d,RBEZIER_PATCH *Patch)
  114. {
  115.   int        i, j;
  116.   DBL        *x,*y,*z;
  117.   DBL        *ptmp, *wptr;
  118.   DBL        a,b,c;
  119.   DBL        min, max;
  120.  
  121.   ptmp=  *p;
  122.   x= Patch->pts[X];
  123.   y= Patch->pts[Y];
  124.   z= Patch->pts[Z];
  125.   wptr= Patch->Weights;
  126.   a= n[X]; b= n[Y]; c=n[Z];
  127.   min=     INFINITY;
  128.   max= -INFINITY;
  129.  
  130.   for ( i=  cv_order; i; i--, ptmp+= c4_u_order, wptr+=c4_u_order,
  131.                           x+=c4_u_order, y+=c4_u_order, z+=c4_u_order)
  132.     for ( j=cu_order; j; j--, wptr++, ptmp++, x++, y++, z++)
  133.       {
  134.     *ptmp = (*x)*a + (*y)*b + (*z)*c + (*wptr)*d;
  135.     if (*ptmp < min) min= *ptmp;
  136.     if (*ptmp > max) max= *ptmp;
  137.       }    
  138.  
  139.   if ((min>0.0)||(max<0.0))
  140.     return 0;
  141.   else 
  142.     return 1;
  143. }
  144.  
  145.  
  146. static void
  147. rbezier_precompute_weights(RBEZIER_PATCH *Patch)
  148. {
  149.   int     i,j;
  150.   DBL    *x,*y,*z, *wptr;
  151.  
  152.   Patch-> w_precomp=1; 
  153.   wptr= Patch->Weights;
  154.   x= Patch->pts[X];
  155.   y= Patch->pts[Y];
  156.   z= Patch->pts[Z];
  157.   for (i=cv_order; i; i--, x+=c4_u_order, y+=c4_u_order, z+=c4_u_order, wptr+=c4_u_order)
  158.     for( j= cu_order; j; j--, wptr++, x++, y++, z++)
  159.       {
  160.     *x *= *wptr;
  161.     *y *= *wptr;
  162.     *z *= *wptr;
  163.       }
  164. }
  165.  
  166. /*********************************************************
  167.  
  168.    rbezier_new_line[0|1] functions          
  169.  
  170. **********************************************************/
  171.  
  172. static int rbezier_new_line0(DISTANCES p[2],DBL      Line[2])
  173. {
  174.   DBL tmp;
  175.  
  176.   Line[0]=    p[1][cv_order_1][0] - p[1][0][0] + p[1][cv_order_1][cu_order_1] - p[1][0][cu_order_1];
  177.   Line[1]= - (p[0][cv_order_1][0] - p[0][0][0] + p[0][cv_order_1][cu_order_1] - p[0][0][cu_order_1]);
  178.  
  179.   tmp= Line[0]*Line[0]+Line[1]*Line[1];
  180.  
  181.   if (tmp < VERY_SMALL_EPSILON)
  182.     return 0;
  183.  
  184.   tmp= sqrt(tmp);
  185.  
  186.   Line[0]/= tmp;
  187.   Line[1]/= tmp;    
  188.   return 1;
  189. }
  190.  
  191. static int rbezier_new_line1(DISTANCES p[2],DBL  Line[2])
  192. {
  193.   DBL tmp;
  194.  
  195.   Line[0]=    p[1][0][cu_order_1] -  p[1][0][0]  + p[1][cv_order_1][cu_order_1] - p[1][cv_order_1][0];
  196.   Line[1]= -( p[0][0][cu_order_1] -  p[0][0][0]  + p[0][cv_order_1][cu_order_1] - p[0][cv_order_1][0] );
  197.  
  198.   tmp= Line[0]*Line[0]+Line[1]*Line[1];
  199.  
  200.   if (tmp < VERY_SMALL_EPSILON)
  201.     return 0;
  202.  
  203.   tmp= sqrt (tmp);
  204.   Line[0]/= tmp;
  205.   Line[1]/= tmp;
  206.   return 1;
  207. }
  208.  
  209.  
  210. /*********************************************************
  211.  
  212.    rbezier_line_points_distances
  213.  
  214. **********************************************************/
  215.  
  216.  
  217. static void rbezier_line_points_distances(DISTANCES d,DBL vec[2],DISTANCES points[2])
  218. {
  219.   int i,j;
  220.   DBL a,b;
  221.   DBL     *ptmp, *xtmp, *ytmp;
  222.  
  223.   a=vec[0];
  224.   b=vec[1];  
  225.  
  226.   ptmp= *d;
  227.   xtmp= **points;
  228.   ytmp= *(points[1]);
  229.  
  230.   for (i=cv_order; i;i--,ptmp+=c4_u_order, xtmp+= c4_u_order, ytmp+= c4_u_order)
  231.     for (j=cu_order; j; j--, ptmp++, xtmp++, ytmp++)
  232.       *ptmp= (*xtmp)*a + (*ytmp)*b;
  233. }
  234.  
  235.  
  236. /*********************************************************
  237.  
  238.    rbezier_bounds          
  239.  
  240. **********************************************************/
  241.  
  242. static int rbezier_bounds( DISTANCES p[2],DBL      interval[2][2])
  243. {
  244.   int i,j,k;
  245.   
  246.   DBL *itmp;
  247.   DBL *min, *max;
  248.  
  249.   int ret=0;
  250.   
  251.   for (i=0; i<2; i++)
  252.     {
  253.       itmp = &(p[i][0][0]);
  254.       min= & (interval[i][0]);
  255.       max= & (interval[i][1]);
  256.       *min=*max= *(itmp);
  257.       for (j= cv_order; j; j--, itmp+= c4_u_order)
  258.     for(k= cu_order; k  ; k--, itmp++ )
  259.       {    
  260.         if ( *itmp < *min )
  261.           *min= *itmp;
  262.         else
  263.           if ( *itmp > *max )
  264.         *max= *itmp;
  265.       }
  266.       ret|= (*min > EPSILON) || (*max < -EPSILON);
  267.     }
  268.   return !ret;
  269. }
  270.  
  271. /*********************************************************
  272.  
  273.   rb_decasteljau_D1_D2 functions
  274.  
  275.     D1 - direction
  276.     D2 - order + 1 (ie number of control points
  277.   
  278.          p  -  points
  279.      f  -  first
  280.      s  -  second
  281.      w  -  where to split
  282.     dw  -  1.0 - dw     
  283. **********************************************************/
  284.  
  285. #define NEXT_0(T)    ((T)++)
  286. #define NEXT_1(T)    ((T)+=4)
  287.  
  288. #define PREV_0(T)    ((T)--)
  289. #define PREV_1(T)    ((T)-=4)
  290.  
  291. #define ADD_0(T,X)    ((T)+=(X))
  292. #define ADD_1(T,X)    ((T)+=(X)*4)
  293.  
  294. #define AV(F,S)        ( dw * (F) + w * (S)) 
  295.  
  296. /*
  297.     2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
  298.  
  299.                   (swans)
  300.               
  301. */
  302.  
  303. static void rb_decasteljau_0_2( DBL *p,DBL *f,DBL *s,DBL w,DBL dw,int i)
  304. {
  305.   NEXT_0(s);
  306.   for(; i; i--)
  307.     {
  308.       *f= *p;
  309.       NEXT_0(p);
  310.       *s= *p;
  311.       PREV_0(s);
  312.       *s= AV( *f,*p );
  313.       NEXT_0(f);
  314.       *f= *s;
  315.       s+= 5;
  316.       f+= 3;
  317.       p+= 3;
  318.     }
  319. }
  320.  
  321. static void rb_decasteljau_1_2( DBL *p,DBL *f,DBL *s,DBL w,DBL dw,int i)
  322. {
  323.   NEXT_1(s);
  324.   for(; i; i--)
  325.     {
  326.       *f= *p;
  327.       NEXT_1(p);
  328.       *s= *p;
  329.       PREV_1(s);
  330.       *s= AV( *f,*p );
  331.       NEXT_1(f);
  332.       *f= *s;
  333.  
  334.       s+= 5;
  335.       f-= 3;
  336.       p-= 3;
  337.     }
  338. }
  339.  
  340.  
  341. /*
  342.   3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
  343.   
  344.                      .....
  345.  */
  346.  
  347. static void rb_decasteljau_0_3(  DBL *p,DBL *f,DBL *s,DBL w,DBL dw,int i)
  348. {
  349.   ADD_0( s , 2);
  350.  
  351.   for (; i ; i--)
  352.     {
  353.       *f= *p;
  354.  
  355.       NEXT_0(f);
  356.       *f= AV( *p, *(p+1));
  357.       ADD_0( p, 2);
  358.  
  359.       *s= *p;
  360.       PREV_0(s);
  361.       *s= AV( *(p-1), *p);
  362.  
  363.       PREV_0(s);
  364.       *s= AV( *f, *(s+1) );
  365.       NEXT_0(f);
  366.       *f= *s;
  367.  
  368.       f+= 2;
  369.       p+= 2;
  370.       s+= 6;
  371.     }
  372. }
  373.  
  374. static void rb_decasteljau_1_3(  DBL *p,DBL *f,DBL *s,DBL w,DBL dw,int i)
  375. {
  376.   ADD_1( s , 2 );
  377.  
  378.   for (; i ; i--)
  379.     {
  380.       *f= *p;
  381.  
  382.       NEXT_1(f);
  383.       *f= AV( *p, *(p+4));
  384.       ADD_1( p, 2);
  385.  
  386.       *s= *p;
  387.       PREV_1(s);
  388.       *s= AV( *(p-4), *p);
  389.  
  390.       PREV_1(s);
  391.       *s= AV( *f, *(s+4) );
  392.       NEXT_1(f);
  393.       *f= *s;
  394.  
  395.       f+= 1 - 2 * 4;
  396.       p+= 1 - 2 * 4;
  397.       s+= 1 + 2 * 4;
  398.     }
  399. }
  400.  
  401. /*
  402.    4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
  403.  
  404.                    (chairs)
  405. */
  406.  
  407. static void rb_decasteljau_0_4( DBL *p,DBL *f,DBL *s,DBL w,DBL dw,int i)
  408. {
  409.   DBL     tmp;
  410.   ADD_0( s, 3 );
  411.  
  412.   for (; i; i--)
  413.     {
  414.       *f= *p;
  415.       NEXT_0(f);
  416.       *f= AV(*p, *(p+1));
  417.  
  418.       NEXT_0(f);
  419.       NEXT_0(p);
  420.       tmp= AV(*p, *(p+1));
  421.       *f= AV(*(f-1), tmp );
  422.  
  423.       NEXT_0(p);
  424.       *s= *(p+1);
  425.  
  426.       PREV_0( s );
  427.       *s= AV( *p, *(s+1) );
  428.       PREV_0( s );
  429.  
  430.       *s= AV( tmp, *(s+1));
  431.  
  432.       PREV_0(s);
  433.       *s= AV( *f, *(s+1) );
  434.       NEXT_0(f);
  435.       *f=*s;
  436.  
  437.       f++;
  438.       p+= 2;   /* p was on [.][2] */
  439.       s+= 7;
  440.     }
  441. }
  442.  
  443. static void rb_decasteljau_1_4(  DBL *p,DBL *f,DBL *s,DBL w,DBL dw,int i)
  444. {
  445.   DBL     tmp;
  446.   ADD_1( s, 3);
  447.  
  448.   for (; i; i--)
  449.     {
  450.       *f= *p;
  451.       NEXT_1(f);
  452.       *f= AV(*p, *(p+4));
  453.  
  454.       NEXT_1(f);
  455.       NEXT_1(p);
  456.       tmp= AV(*p, *(p+4));
  457.       *f= AV(*(f-4), tmp);
  458.  
  459.       NEXT_1(p);
  460.       *s= *(p+4);
  461.  
  462.       PREV_1( s ); 
  463.       *s= AV( *p, *(s+4) );
  464.       PREV_1( s );
  465.  
  466.       *s= AV( tmp, *(s+4));
  467.  
  468.       PREV_1(s);
  469.       *s= AV( *f, *(s+4) );
  470.       NEXT_1(f);
  471.       *f=*s;
  472.  
  473.       f+= 1 - 3 * 4;
  474.       p+= 1 - 2 * 4;   /* p was on [2][.] */
  475.       s+= 1 + 3 * 4;
  476.     }
  477. }
  478.  
  479. /*********************************************************
  480.  
  481.   rbezier_get_minmax_D1_D2 functions    
  482.   
  483.  **********************************************************/
  484.  
  485. #define FIND_INTER( P, N, P0, N0 )             \
  486.   if ( ( (P) <= 0.0) ^ ( (N)  <= 0.0) )            \
  487.    {                            \
  488.      tmp=  (P0) + ((N0)-(P0)) * (P) / ( (P) - (N) );    \
  489.                             \
  490.      if (tmp < *min) *min= tmp;              \
  491.      if (tmp > *max) *max= tmp;                \
  492.    }                            
  493.  
  494. static int rb_get_minmax_0_2(DISTANCES d,DBL *min, DBL *max)
  495. {
  496.   DBL *n,*p,tmp;
  497.  
  498.   int  f,l;    
  499.   int  i;
  500.  
  501.   *min = 2.0;
  502.   *max =-2.0;
  503.  
  504.   n= 1 + ( p= *d );
  505.   f=l=0;
  506.  
  507.   for (i= cv_order; i; i--, NEXT_1(p), NEXT_1(n))
  508.     {
  509.       if ( *p < 0.0) f++;
  510.       if ( *n < 0.0) l++;
  511.       FIND_INTER( *p, *n, 0.0, 1.0 );
  512.     }
  513.  
  514.   if (( f != 0 ) && ( f != cv_order )) *min = 0.0;
  515.   if (( l != 0 ) && ( l != cv_order )) *max = 1.0;
  516.  
  517.   return ( *min <=  * max );
  518. }
  519.  
  520. static int rb_get_minmax_1_2(DISTANCES d,DBL *min, DBL *max)
  521. {
  522.   DBL *n,*p,tmp;
  523.  
  524.   int  f,l;    
  525.   int  i;
  526.  
  527.   *min = 2.0;
  528.   *max =-2.0;
  529.  
  530.   n= 4 + ( p= *d );
  531.   f=l=0;
  532.  
  533.   for (i= cu_order; i; i--, NEXT_0(p), NEXT_0(n))
  534.     {
  535.       if ( *p < 0.0) f++;
  536.       if ( *n < 0.0) l++;
  537.       FIND_INTER( *p, *n, 0.0, 1.0 );
  538.     }
  539.  
  540.   if (( f != 0 ) && ( f != cu_order )) *min = 0.0;
  541.   if (( l != 0 ) && ( l != cu_order )) *max = 1.0;
  542.  
  543.   return ( *min <=  * max );
  544. }
  545.  
  546. static int rb_get_minmax_0_3(DISTANCES d,DBL *min, DBL *max)
  547. {
  548.   DBL *n,*p, *m,tmp;
  549.  
  550.   int  f,l;    
  551.   int  i;
  552.  
  553.   *min = 2.0;
  554.   *max =-2.0;
  555.  
  556.   n= ( 1+ (m= 1 + ( p= *d )));
  557.   f=l=0;
  558.  
  559.   for (i= cv_order; i; i--, NEXT_1(p), NEXT_1(m), NEXT_1(n))
  560.     {
  561.       if ( *p < 0.0) f++;
  562.       if ( *n < 0.0) l++;
  563.  
  564.       FIND_INTER( *p, *m, 0.0, 0.5 );
  565.       FIND_INTER( *m, *n, 0.5, 1.0 );
  566.       FIND_INTER( *p, *n, 0.0, 1.0 );
  567.     }
  568.  
  569.   if (( f != 0 ) && ( f != cv_order )) *min = 0.0;
  570.   if (( l != 0 ) && ( l != cv_order )) *max = 1.0;
  571.  
  572.   return ( *min <=  * max );
  573. }
  574.  
  575. static int rb_get_minmax_1_3(DISTANCES d,DBL *min, DBL *max)
  576. {
  577.   DBL *n,*p, *m,tmp;
  578.  
  579.   int  f,l;    
  580.   int  i;
  581.  
  582.   *min = 2.0;
  583.   *max =-2.0;
  584.  
  585.   n= ( 4 + (m= 4 + ( p= *d )));
  586.   f=l=0;
  587.  
  588.   for (i= cu_order; i; i--, NEXT_0(p), NEXT_0(m), NEXT_0(n))
  589.     {
  590.       if ( *p < 0.0) f++;
  591.       if ( *n < 0.0) l++;
  592.  
  593.       FIND_INTER( *p, *m, 0.0, 0.5 );
  594.       FIND_INTER( *m, *n, 0.5, 1.0 );
  595.       FIND_INTER( *p, *n, 0.0, 1.0 );
  596.     }
  597.  
  598.   if (( f != 0 ) && ( f != cu_order )) *min = 0.0;
  599.   if (( l != 0 ) && ( l != cu_order )) *max = 1.0;
  600.  
  601.   return ( *min <=  * max );
  602. }
  603.  
  604. static int rb_get_minmax_0_4(DISTANCES d,DBL *min, DBL *max)
  605. {
  606.   DBL *n,*p, *mn, *mp,tmp;
  607.  
  608.   int  f,l;    
  609.   int  i;
  610.  
  611.   *min = 2.0;
  612.   *max =-2.0;
  613.  
  614.   n= 1 + (mn=  1+ (mp= 1 + ( p= *d )));
  615.   f=l=0;
  616.  
  617.   for (i= cv_order; i; i--, NEXT_1(p), NEXT_1(mp), NEXT_1(mn), NEXT_1(n))
  618.     {
  619.       if ( *p < 0.0) f++;
  620.       if ( *n < 0.0) l++;
  621.  
  622.       FIND_INTER( * p, *mp,       0.0, (1.0/3.0) );
  623.       FIND_INTER( *mp, *mn, (1.0/3.0), (2.0/3.0) );
  624.       FIND_INTER( * p, *mn,       0.0, (2.0/3.0) );
  625.  
  626.       FIND_INTER( *mn, * n, (2.0/3.0),      1.0  );
  627.       FIND_INTER( *mp, * n, (1.0/3.0),      1.0  );
  628.  
  629.       FIND_INTER( * p, * n,      0.0 ,      1.0  );
  630.     }
  631.  
  632.   if (( f != 0 ) && ( f != cv_order )) *min = 0.0;
  633.   if (( l != 0 ) && ( l != cv_order )) *max = 1.0;
  634.  
  635.   return ( *min <=  * max );
  636. }
  637.  
  638. static int rb_get_minmax_1_4(DISTANCES d,DBL *min, DBL *max)
  639. {
  640.   DBL *n,*p, *mn, *mp,tmp;
  641.  
  642.   int  f,l;    
  643.   int  i;
  644.  
  645.   *min = 2.0;
  646.   *max =-2.0;
  647.  
  648.   n= 4 + (mn=  4+ (mp= 4 + ( p= *d )));
  649.   f=l=0;
  650.  
  651.   for (i= cu_order; i; i--, NEXT_0(p), NEXT_0(mp), NEXT_0(mn), NEXT_0(n))
  652.     {
  653.       if ( *p < 0.0) f++;
  654.       if ( *n < 0.0) l++;
  655.  
  656.       FIND_INTER( * p, *mp,       0.0, (1.0/3.0) );
  657.       FIND_INTER( *mp, *mn, (1.0/3.0), (2.0/3.0) );
  658.       FIND_INTER( * p, *mn,       0.0, (2.0/3.0) );
  659.  
  660.       FIND_INTER( *mn, * n, (2.0/3.0),      1.0  );
  661.       FIND_INTER( *mp, * n, (1.0/3.0),      1.0  );
  662.  
  663.       FIND_INTER( * p, * n,      0.0 ,      1.0  );
  664.     }
  665.  
  666.   if (( f != 0 ) && ( f != cu_order )) *min = 0.0;
  667.   if (( l != 0 ) && ( l != cu_order )) *max = 1.0;
  668.  
  669.   return ( *min <=  * max );
  670. }
  671.  
  672. static
  673.   int (*rbezier_get_minmax[2][4])(DISTANCES, DBL*, DBL*)= {
  674.     { NULL, rb_get_minmax_0_2, rb_get_minmax_0_3, rb_get_minmax_0_4 },
  675.     { NULL, rb_get_minmax_1_2, rb_get_minmax_1_3, rb_get_minmax_1_4 }
  676.   };
  677.  
  678. static
  679.   void (*rbezier_decasteljau[2][4])(DBL*, DBL*, DBL *, DBL, DBL, int ) = {
  680.     { NULL, rb_decasteljau_0_2, rb_decasteljau_0_3, rb_decasteljau_0_4 },
  681.     { NULL, rb_decasteljau_1_2, rb_decasteljau_1_3, rb_decasteljau_1_4 }
  682.   };
  683.  
  684.  
  685. /**************************************************
  686.  
  687.   rbezier_point 
  688.  
  689. **************************************************/
  690.  
  691. static void rbezier_point(DBL u,DBL v,int iii,DBL *cp,DBL *p0,DBL *p10,DBL *p11 )
  692. {
  693.   DISTANCES  d1, d2;
  694.  
  695.   if (iii == 0)
  696.     {
  697.       rbezier_decasteljau[0][cu_order_1]( cp, &(d1[0][0]), &(d2[0][0]),u,1.0-u, cv_order);
  698.       rbezier_decasteljau[1][cv_order_1]( &(d2[0][0]), &(d2[0][1]), &(d2[0][2]), v, 1.0-v, 1);
  699.  
  700.       *p0=  d2[0][2];
  701.       *p10= d2[1][2];
  702.       *p11= d2[cv_order - 2][1];
  703.     }
  704.   else
  705.     {
  706.       rbezier_decasteljau[1][cv_order_1]( cp, &(d1[0][0]), &(d2[0][0]),v,1.0-v, cu_order);
  707.       rbezier_decasteljau[0][cu_order_1]( &(d2[0][0]), &(d2[1][0]), &(d2[2][0]), u, 1.0-u, 1);
  708.  
  709.       *p0=  d2[2][0];
  710.       *p10= d2[2][1];
  711.       *p11= d2[1][cu_order - 2];
  712.     }
  713. }
  714.  
  715. /**************************************************
  716.  
  717.   rbezier_add_solution
  718.  
  719. **************************************************/
  720.  
  721. static int rbezier_add_solution(DBL u,DBL v)
  722. {
  723.   int     coord;
  724.   DBL     *cp, d;
  725.  
  726.   int        iii;
  727.   DBL       c1,c2,c3;
  728.   DBL        w1,w2,w3;
  729.   VECTOR     IPoint;
  730.   ISTACK    *i;
  731.  
  732.   if ((u<0.0) || (u> 1.0) || (v < 0.0) || ( v > 1.0)) return 0;
  733.  
  734.   if (fabs(Current_Ray->Direction[X]) > 0.5 ) 
  735.      { coord= X; cp = Current_Shape -> pts[X]; }
  736.    else
  737.      if (fabs(Current_Ray->Direction[Y]) > 0.5 ) 
  738.        { coord= Y; cp = Current_Shape -> pts[Y]; }
  739.      else
  740.        { coord= Z; cp = Current_Shape -> pts[Z]; }
  741.    
  742.    iii= cu_order > cv_order;
  743.  
  744.    rbezier_point( u, v, iii,  cp, &c1, &c2, &c3);
  745.  
  746.    if ( Current_Shape -> Weights )
  747.      {
  748.        rbezier_point( u, v, iii, Current_Shape->Weights, &w1, &w2, &w3 );
  749.        c1/= w1;
  750.        c2= (c2 / w2 - c3 / w3);
  751.      }
  752.    else
  753.      c2= c2 -c3;
  754.  
  755.    d= (c1 - Current_Ray->Initial[coord]) / (Current_Ray->Direction[coord]);
  756.  
  757.    if ( d  < 4.0 * Current_Shape -> accuracy) return 0;
  758.  
  759.    if (Current_Shape->num_of_trims != 0)
  760.      if (rbezier_check_trims(u,v) == 0) return 0;
  761.  
  762.    VScale (IPoint, Current_Ray->Direction, d);
  763.    VAddEq (IPoint, Current_Ray->Initial);
  764.  
  765.  
  766.    i= Current_Depth_Stack;
  767.    /* 
  768.       I was too lazy to create own macro in frame.h :-) 
  769.       so you can imagine this is macro push_entry_blah_blah 
  770.     */
  771.  
  772.   itop(i).Depth  = d;                     
  773.   itop(i).Object = (OBJECT *) Current_Shape;                      
  774.   itop(i).i1 = coord;     
  775.   itop(i).i2 = iii;
  776.  
  777.   itop(i).u = u;                          
  778.   itop(i).v = v;
  779.                        
  780.   itop(i).d1 = c1;                          
  781.   itop(i).d2 = c2;
  782.  
  783.   if (Current_Shape -> Weights)
  784.     {
  785.       itop(i).w1 = w1;                          
  786.       itop(i).w2 = w2;                          
  787.       itop(i).w3 = w3;                          
  788.     }
  789.  
  790.   Assign_Vector(itop(i).IPoint,IPoint);         
  791.   incstack(i);                       
  792.  
  793.   /* end of `virtual' macro :-) */
  794.  
  795.   return 1; 
  796. }
  797.  
  798.  
  799. /*********************************************************
  800.  
  801.    rbezier_find_solution          
  802.  
  803. **********************************************************/
  804.  
  805. static int rbezier_find_solution( DISTANCES    p[2],DBL interval[2][2], DBL bounds[2][2])
  806. {
  807.   int        iii, yyy,ppp;
  808.   DBL        Line[2];
  809.   DBL        tmp;
  810.  
  811.   char        in_half=0;
  812.  
  813.   DISTANCES     d;
  814.   DBL        min,max;
  815.  
  816.   
  817.   if ( (bounds[1][1]-bounds[1][0] < Current_Shape->accuracy)
  818.        && (bounds[0][1]-bounds[0][0] < Current_Shape->accuracy))
  819.     {   /* end of iteration ...*/
  820.       return rbezier_add_solution( interval[0][0]+ interval[0][1]/2.0, interval[1][0] + interval[1][1]/2.0 );
  821.       FIXME  /* soupni sem presnejsi resitko ...!!! */
  822.     }
  823.  
  824.   if ( interval[0][1] > interval[1][1] )
  825.     {
  826.       iii=0;
  827.       if ( rbezier_new_line0(p, Line) == 0)
  828.     if ( rbezier_new_line1(p, Line) == 0)
  829.       in_half=1;   /* degenerated patch -=> split in half */
  830.     else
  831.       {
  832.         tmp= -Line[0];
  833.         Line[0]= Line[1];
  834.         Line[1]= tmp;
  835.       }
  836.     }
  837.   else
  838.     {
  839.       iii=1;
  840.       if ( rbezier_new_line1(p, Line) == 0)
  841.         if ( rbezier_new_line0(p, Line) == 0)
  842.           in_half=1;   /* degenerated patch -=> split in half */
  843.         else
  844.           {
  845.             tmp= -Line[0];
  846.             Line[0]= Line[1];
  847.             Line[1]= tmp;
  848.           }
  849.     }
  850.  
  851.   yyy= (iii == 0) ? cu_order_1 : cv_order_1;
  852.   ppp= (iii != 0) ? cu_order   : cv_order;
  853.  
  854.   if (in_half == 0)
  855.     {
  856.       rbezier_line_points_distances( d, Line, p);
  857.       if ( rbezier_get_minmax[iii][yyy](d,&min, &max)==0) return 0;
  858.       if ( (tmp= max-min) > MAGICAL_CONSTANT )
  859.     in_half= 1;
  860.       else
  861.     {
  862.       {
  863.         DISTANCES ptmp1[2],ptmp2[2];
  864.  
  865.  
  866.         if (min < 0.5) 
  867.           {
  868.         rbezier_decasteljau[iii][yyy]( &(p[0][0][0]), &(ptmp1[0][0][0]), &(ptmp2[0][0][0]), min, 1.0 - min, ppp);
  869.         rbezier_decasteljau[iii][yyy]( &(p[1][0][0]), &(ptmp1[1][0][0]), &(ptmp2[1][0][0]), min, 1.0 - min, ppp);
  870.  
  871.         rbezier_decasteljau[iii][yyy]( &(ptmp2[0][0][0]), &(p[0][0][0]), &(ptmp1[0][0][0]), 
  872.                            (max-min) / (1.0 - min), (1.0 - max) / (1.0 - min),  ppp);
  873.         rbezier_decasteljau[iii][yyy]( &(ptmp2[1][0][0]), &(p[1][0][0]), &(ptmp1[1][0][0]), 
  874.                            (max-min) / (1.0 - min), (1.0 - max) / (1.0 - min),  ppp);
  875.           }
  876.         else 
  877.           {
  878.         rbezier_decasteljau[iii][yyy]( &(p[0][0][0]), &(ptmp1[0][0][0]), &(ptmp2[0][0][0]), max, 1.0 - max, ppp);
  879.         rbezier_decasteljau[iii][yyy]( &(p[1][0][0]), &(ptmp1[1][0][0]), &(ptmp2[1][0][0]), max, 1.0 - max, ppp);
  880.  
  881.         rbezier_decasteljau[iii][yyy]( &(ptmp1[0][0][0]), &(ptmp2[0][0][0]), &(p[0][0][0]), 
  882.                            min / max, 1.0 - min / max, ppp);
  883.         rbezier_decasteljau[iii][yyy]( &(ptmp1[1][0][0]), &(ptmp2[1][0][0]), &(p[1][0][0]), 
  884.                            min / max, 1.0 - min/max, ppp);
  885.           }
  886.       }
  887.  
  888.       if (rbezier_bounds(p, bounds ) == 0) return 0;  
  889.      
  890.       interval[iii][0]+= interval[iii][1]*min;
  891.       interval[iii][1] = interval[iii][1]*tmp;
  892.       
  893.       return rbezier_find_solution(p, interval, bounds);
  894.     }
  895.     }
  896.  
  897.       /*  (in_half == 1) */
  898.   {
  899.     DISTANCES ld[2], rd[2];
  900.     DBL    ni[2][2];
  901.     int    ret=0;
  902.     
  903.     rbezier_decasteljau[iii][yyy]( &(p[0][0][0]), &(ld[0][0][0]), &(rd[0][0][0]), 0.5, 0.5 ,ppp);
  904.     rbezier_decasteljau[iii][yyy]( &(p[1][0][0]), &(ld[1][0][0]), &(rd[1][0][0]), 0.5, 0.5 ,ppp);
  905.     
  906.     tmp= interval[iii][1]/= 2.0;
  907.     memcpy( &(ni[0][0]), &(interval[0][0]), 4 * sizeof(DBL));
  908.     ni[iii][0]+= tmp;
  909.     
  910.     if (rbezier_bounds(ld,bounds) != 0) 
  911.       ret|= rbezier_find_solution(ld, interval, bounds);
  912.  
  913.     if (rbezier_bounds(rd, bounds) != 0)
  914.       ret|= rbezier_find_solution(rd, ni, bounds);
  915.       
  916.     return ret;
  917.   }
  918. }
  919.  
  920.  
  921. /*********************************************************
  922.  
  923.    All_RBezier_Patch_Intersections
  924.  
  925. **********************************************************/
  926.  
  927. static int All_RBezier_Patch_Intersections(OBJECT *Object,RAY *Ray,ISTACK *Depth_Stack)
  928. {
  929.   VECTOR     n0,n1;
  930.   DBL         d0, d1;
  931.  
  932.   DISTANCES    p[2];
  933.   DBL        interval[2][2];
  934.   DBL        sti[2][2];
  935.   RBEZIER_PATCH *Patch;
  936.   int        ret;
  937.  
  938. #ifdef DEBUG_RBEZIER
  939.   if (((RBEZIER_PATCH *)Object)->debug)
  940.     if (Ray->Direction[Z] != 1.0) return 0; 
  941. #endif
  942.   
  943.   Increase_Counter(stats[Ray_RBezier_Tests]);    
  944.  
  945.   Patch= (RBEZIER_PATCH *) Object;
  946.   
  947.   find_planes( Ray->Initial, Ray->Direction, n0,&d0,n1,&d1);
  948.  
  949.   cu_order= Patch->u_order;
  950.   cv_order= Patch->v_order;
  951.   
  952.   cu_order_1= cu_order-1;
  953.   cv_order_1= cv_order-1;
  954.  
  955.   c4_u_order= 4 - cu_order;
  956.  
  957.   if ( Patch -> Weights == NULL)
  958.     {
  959.      if (bezier_plane_points_distances(p[0], n0, d0, Patch) == 0)
  960.        return 0;
  961.      if (bezier_plane_points_distances(p[1], n1, d1, Patch) == 0)
  962.        return 0;
  963.    }
  964.  else
  965.    {
  966.      if (Patch-> w_precomp == 0)
  967.        {
  968.      rbezier_precompute_weights(Patch);
  969.        }
  970.      if (rbezier_plane_points_distances(p[0], n0, d0, Patch) == 0)
  971.        return 0;
  972.      if (rbezier_plane_points_distances(p[1], n1, d1, Patch) == 0)
  973.        return 0;
  974.    }
  975.  
  976.   interval[0][0]=interval[1][0]= 0.0;
  977.   if ( cu_order > cv_order)
  978.     { 
  979.       /* small trick to ensure that we first split patch in
  980.      smaller order */
  981.       interval[0][1]= 1.0 + YAEP; 
  982.       interval[1][1]= 1.0;
  983.     }
  984.   else
  985.     {
  986.       interval[0][1]= 1.0;
  987.       interval[1][1]= 1.0 + YAEP; 
  988.     }
  989.  
  990.   sti[1][1]= sti[0][1]=  INFINITY;
  991.   sti[1][0]= sti[0][0]= -INFINITY;
  992.  
  993.   Current_Ray= Ray;
  994.   Current_Shape= Patch;
  995.   Current_Depth_Stack= Depth_Stack; 
  996.  
  997.   ret= rbezier_find_solution( p, interval, sti );
  998.   if (ret) 
  999.     Increase_Counter(stats[Ray_RBezier_Tests_Succeeded]);
  1000.   return ret;
  1001. }
  1002.  
  1003. static void RBezier_Patch_Normal(VECTOR Result, OBJECT *Object, INTERSECTION *Inter)
  1004. {
  1005.   VECTOR v1,v2;
  1006.   int      i, rat;
  1007.   DBL     u,v;
  1008.   DBL     c1,c2,c3;
  1009.   DBL     w1,w2,w3;
  1010.   RBEZIER_PATCH *Patch;
  1011.   int     iii;
  1012.  
  1013.   Patch= (RBEZIER_PATCH *) Object;
  1014.  
  1015.   cu_order= Patch-> u_order;
  1016.   cv_order= Patch-> v_order;
  1017.  
  1018.   cu_order_1= cu_order-1;
  1019.   cv_order_1= cv_order-1;
  1020.  
  1021.  
  1022.   iii= Inter-> i2;
  1023.   u= Inter->u;
  1024.   v= Inter->v;
  1025.  
  1026.   if ( (rat= (Patch ->Weights != NULL))) 
  1027.      {  w2= Inter->w2;  w3= Inter->w3; } 
  1028.     
  1029.   for ( i = 0; i <=Z; i++)
  1030.     if ( i == Inter -> i1 ) /* computed in rbezier_add_solution */
  1031.       v1[i]= Inter->d2;
  1032.     else
  1033.       {
  1034.     rbezier_point( u, v, iii, Patch->pts[i], &c1, &c2, &c3);
  1035.     if (rat)
  1036.       v1[i]= (c2 / w2 - c3 / w3 );
  1037.     else
  1038.       v1[i]= c2 - c3;
  1039.       }
  1040.  
  1041.   iii^= 1;
  1042.  
  1043.   if (rat)
  1044.     rbezier_point( u, v, iii,  Patch->Weights, &w1, &w2, &w3);
  1045.  
  1046.   for ( i = 0; i <=Z; i++)
  1047.       {
  1048.     rbezier_point( u, v, iii, Patch->pts[i], &c1, &c2, &c3);
  1049.     if (rat)
  1050.       v2[i]= (c2 / w2 - c3 / w3 );
  1051.     else
  1052.       v2[i]= c2 - c3;
  1053.       }
  1054.     
  1055.   VCross(Result, v1,v2);
  1056.   VLength(w2, Result);
  1057.  
  1058.   if (w2 < VERY_SMALL_EPSILON)
  1059.     { Result[X]= Result[Y]= 0.0; Result[Z]= 1.0; }
  1060.   else
  1061.     { Result[X]/= w2; Result[Y]/= w2; Result[Z]/= w2; }
  1062. }
  1063.  
  1064. /***********************************************************
  1065.  
  1066.    other "not so funny functions" (like transform, create etc)
  1067.  
  1068. ************************************************************/
  1069.  
  1070. static int Inside_RBezier_Patch(VECTOR IPoint,OBJECT *Object)
  1071. {
  1072.   return (0);
  1073. }
  1074.  
  1075. static void Translate_RBezier_Patch(OBJECT *Object,VECTOR Vector,TRANSFORM *Trans)
  1076. {
  1077.   int i,n;
  1078.  
  1079.   RBEZIER_PATCH *Patch = (RBEZIER_PATCH *) Object;
  1080.   n= 4 * Patch -> v_order;
  1081.  
  1082.   for (i=0; i < n; i++)
  1083.     {
  1084.       Patch->Points[X*n+i]+= Vector[X];
  1085.       Patch->Points[Y*n+i]+= Vector[Y];
  1086.       Patch->Points[Z*n+i]+= Vector[Z];
  1087.     }          
  1088.  
  1089.   Compute_RBezier_Patch_BBox(Patch);
  1090. }
  1091.  
  1092. static void Scale_RBezier_Patch(OBJECT *Object,VECTOR Vector,TRANSFORM *Trans)
  1093. {
  1094.   int i,n;
  1095.  
  1096.   RBEZIER_PATCH *Patch = (RBEZIER_PATCH *) Object;
  1097.   n= 4 * Patch-> v_order;
  1098.  
  1099.   for (i=0; i < n; i++)
  1100.     {      
  1101.       Patch->Points[X*n+i]*= Vector[X];
  1102.       Patch->Points[Y*n+i]*= Vector[Y];
  1103.       Patch->Points[Z*n+i]*= Vector[Z];
  1104.     }
  1105.  
  1106.   Compute_RBezier_Patch_BBox(Patch);
  1107. }
  1108.  
  1109. static void Rotate_RBezier_Patch (OBJECT *Object,VECTOR Vector,TRANSFORM *Trans)
  1110. {
  1111.   Transform_RBezier_Patch(Object, Trans);
  1112. }
  1113.  
  1114. static void Transform_RBezier_Patch(OBJECT *Object,TRANSFORM *Trans)
  1115. {
  1116.   int i,n;
  1117.   VECTOR v;
  1118.   RBEZIER_PATCH *Patch = (RBEZIER_PATCH *) Object;
  1119.  
  1120.   n= 4 * Patch->v_order;
  1121.  
  1122.   for (i = 0; i < n; i++)
  1123.     {
  1124.       v[X]= Patch->Points[X*n+i];
  1125.       v[Y]= Patch->Points[Y*n+i];
  1126.       v[Z]= Patch->Points[Z*n+i];
  1127.  
  1128.       MTransPoint(v, v, Trans);
  1129.       
  1130.       Patch->Points[X*n+i]= v[X];      
  1131.       Patch->Points[Y*n+i]= v[Y];      
  1132.       Patch->Points[Z*n+i]= v[Z];      
  1133.     }
  1134.  
  1135.   Compute_RBezier_Patch_BBox (Patch);
  1136. }
  1137.  
  1138.  
  1139. static void Invert_RBezier_Patch(OBJECT *Object)
  1140. {
  1141. }
  1142.  
  1143. RBEZIER_PATCH *Create_RBezier_Patch(void)
  1144. {
  1145.   RBEZIER_PATCH *New;
  1146.   
  1147.   New= (RBEZIER_PATCH *) POV_MALLOC ( sizeof (RBEZIER_PATCH), "bezier patch" );
  1148.   
  1149.  INIT_OBJECT_FIELDS(New, RBEZIER_PATCH_OBJECT, &RBezier_Patch_Methods)
  1150.  
  1151.   New-> Points= NULL;
  1152.   New-> Weights= NULL;
  1153.   New-> trims= NULL;
  1154.   New-> num_of_trims= 0;
  1155.   New-> w_precomp= 0;
  1156.   New-> u_order= 4;  
  1157.   New-> v_order= 4;
  1158.  
  1159.   New-> accuracy= 0.01;
  1160.  
  1161. #ifdef DEBUG_RBEZIER
  1162.   New-> debug= 0;
  1163. #endif 
  1164.  
  1165.   return (New);
  1166. }
  1167.  
  1168. static void *Copy_RBezier_Patch(OBJECT *Object)
  1169. {
  1170.   RBEZIER_PATCH     *New, *Old;
  1171.   int             i;
  1172.  
  1173.   New= Create_RBezier_Patch();
  1174.   Old= (RBEZIER_PATCH *) Object;
  1175.  
  1176.   New->UV_Trans = Copy_Transform( Old->UV_Trans );
  1177.  
  1178.   New->accuracy= Old->accuracy;
  1179.   New->u_order= Old-> u_order;
  1180.   New->v_order= Old-> v_order;
  1181.  
  1182.   i= 4 * Old-> v_order;
  1183.   
  1184.   if (Old -> Points)
  1185.     {
  1186.       New->Points= (DBL *) POV_MALLOC( sizeof(DBL) * (i * 3), "bezier patch");
  1187.       memcpy ( New->Points, Old->Points, sizeof(DBL) * (i * 3) );
  1188.     }
  1189.  
  1190.   if (Old -> Weights)
  1191.     {
  1192.       New->Weights= (DBL *) POV_MALLOC ( sizeof(DBL) * i, "bezier patch");
  1193.       memcpy( New->Points, Old->Points, sizeof(DBL) * i);
  1194.     }
  1195.  
  1196.   New->pts[X]= New->Points;
  1197.   New->pts[Y]= &(New->Points[4*New->v_order]);
  1198.   New->pts[Z]= &(New->Points[2*4*New->v_order]);
  1199.  
  1200.   if (New-> num_of_trims != 0)
  1201.     {
  1202.       for (i=New->num_of_trims; i; i--)
  1203.     {
  1204.       Copy_Trim( New->trims[i] );
  1205.       /* just to incremetn counter in trim */
  1206.     }
  1207.     }
  1208.  
  1209.   return (New);
  1210. }
  1211.  
  1212. static void Destroy_RBezier_Patch(OBJECT *Object)
  1213. {
  1214.   RBEZIER_PATCH *Patch;
  1215.   int        i;
  1216.  
  1217.   Patch = (RBEZIER_PATCH *)Object;
  1218.  
  1219.   if ( Patch -> Points ) POV_FREE( Patch-> Points );
  1220.   if ( Patch -> Weights ) POV_FREE( Patch -> Weights );
  1221.  
  1222.   if (Patch -> trims)
  1223.     {
  1224.       for (i= 0; i < Patch->num_of_trims; i++)
  1225.     Destroy_Trim( Patch -> trims[i] );
  1226.       POV_FREE( Patch->trims );
  1227.     }
  1228.   POV_FREE( Patch );
  1229. }
  1230.  
  1231. void Compute_RBezier_Patch_BBox(RBEZIER_PATCH *Patch)
  1232. {
  1233.   int i,n;
  1234.   VECTOR Min, Max;
  1235.  
  1236.   Make_Vector(Min, BOUND_HUGE, BOUND_HUGE, BOUND_HUGE);
  1237.   Make_Vector(Max, -BOUND_HUGE, -BOUND_HUGE, -BOUND_HUGE);
  1238.  
  1239.   n= 4 * (Patch-> v_order);
  1240.  
  1241.   for (i=0; i < n; i++)
  1242.     {
  1243.       Min[X] = min(Min[X], Patch->Points[X*n+i]);    
  1244.       Min[Y] = min(Min[Y], Patch->Points[Y*n+i]);    
  1245.       Min[Z] = min(Min[Z], Patch->Points[Z*n+i]);    
  1246.       Max[X] = max(Max[X], Patch->Points[X*n+i]);    
  1247.       Max[Y] = max(Max[Y], Patch->Points[Y*n+i]);    
  1248.       Max[Z] = max(Max[Z], Patch->Points[Z*n+i]);    
  1249.     }
  1250.   Make_BBox_from_min_max( Patch-> BBox, Min, Max);
  1251. }
  1252.  
  1253. /*******************************************************
  1254.  
  1255.     functon for handling trimmed surface
  1256.  
  1257.  *******************************************************/
  1258.  
  1259. /**************************************************
  1260.  
  1261.   rbezier_trim
  1262.  
  1263. **************************************************/
  1264.  
  1265. static int trim_choose_direction( DBL u,DBL v )
  1266. {
  1267.   char    dir;
  1268.   DBL    min, tmp;
  1269.   
  1270.   dir=0;
  1271.   min= Current_Trim-> bounds[0][1] - u ;
  1272.   
  1273.   tmp= u - Current_Trim-> bounds[0][0];
  1274.   if ( min > tmp ) { min= tmp; dir = 2; }
  1275.  
  1276.   tmp= Current_Trim-> bounds[1][1] - v;
  1277.   if ( min > tmp ) { min= tmp; dir = 1; }
  1278.  
  1279.   tmp= v - Current_Trim-> bounds[1][0];
  1280.   if ( min > tmp ) { min= tmp; dir = 3; }
  1281.   
  1282.   if (min < 0.0) return -1;
  1283.   else return dir;
  1284. }
  1285.  
  1286. static void trim_translate_points(char     dir,DBL    u,DBL v,DBL *newu,DBL *newv)
  1287. {
  1288.   int     i;
  1289.   DBL    *uu, *vv;
  1290.  
  1291.   uu= Current_Trim-> u;
  1292.   vv= Current_Trim-> v;
  1293.  
  1294.   if ((dir & 1) == 0)
  1295.     if (dir==0)
  1296.       for (i=Current_Trim->num_of_cp; i; *newu=(*uu)-u, (*newv)=(*vv)-v, uu++, vv++, newu++, newv++, i--);
  1297.     else                                 
  1298.       for (i=Current_Trim->num_of_cp; i; *newu=u-(*uu), (*newv)=v-(*vv), uu++, vv++, newu++, newv++, i--);
  1299.   else                                     
  1300.     if (dir==1)                                 
  1301.       for (i=Current_Trim->num_of_cp; i; *newu=(*vv)-v, (*newv)=(*uu)-u, uu++, vv++, newu++, newv++, i--);
  1302.     else                                 
  1303.       for (i=Current_Trim->num_of_cp; i; *newu=v-(*vv), (*newv)=u-(*uu), uu++, vv++, newu++, newv++, i--);
  1304. }
  1305.  
  1306.  
  1307. #define TRIM_FIND_INTER( P, N, P0, N0 )            \
  1308.    {                            \
  1309.      tmp=  (P0) + ((N0)-(P0)) * (P) / ( (P) - (N) );    \
  1310.                             \
  1311.      if (tmp < min) min= tmp;                  \
  1312.      if (tmp > max) max= tmp;                \
  1313.    }                            
  1314.  
  1315. static int trim_check_3 (DBL *uu, DBL *vv,DBL len)
  1316. {
  1317.   int     vc,vcc, uc;
  1318.   DBL    min, max, tmp;
  1319.  
  1320.   vc= (vv[0]<0.0) + (vv[2]<0.0);
  1321.   vcc= vc + (vv[1]<0.0);
  1322.   if ((vcc==0)||(vcc==3)) return 0;
  1323.  
  1324.   uc= (uu[0]<0.0) + (uu[1]<0.0) + (uu[2]<0.0);
  1325.   if (uc==3) return 0;
  1326.   if (uc ==0) 
  1327.     if (vc==1)
  1328.       return 1;
  1329.     else  /* vc from {0,2} */
  1330.       return 0;
  1331.  
  1332.   /*
  1333.     if (len < Current_Shape->accuracy) return 0;
  1334.   */
  1335.  
  1336.   min= 2; max= -2;
  1337.  
  1338.   if ((vv[0]<0.0)^(vv[1]<0.0)) 
  1339.     TRIM_FIND_INTER( vv[0], vv[1],     0.0,  0.5);
  1340.  
  1341.   if ((vv[1]<0.0)^(vv[2]<0.0)) 
  1342.     TRIM_FIND_INTER( vv[1], vv[2],     0.5,  1.0);
  1343.  
  1344.   if (vc == 1)
  1345.     TRIM_FIND_INTER( vv[0], vv[2],     0.0,  1.0);
  1346.  
  1347.   if (min > max) return 0;
  1348.  
  1349.   if ((max-min) < Current_Shape->accuracy)
  1350.     { /* size of new part would be too small -=> stop */
  1351.       DBL uu1[3], uu2[3], vv1[3], vv2[3];
  1352.       rb_decasteljau_0_3( uu, uu1, uu2, min, 1.0-min, 1);
  1353.       rb_decasteljau_0_3( vv, vv1, vv2, min, 1.0-min, 1);
  1354.  
  1355.       if ((uu2[0] >= 0.0) && (fabs(vv2[0]) < Current_Shape-> accuracy)) return 1;
  1356.       else return 0;
  1357.     }
  1358.   if ((max-min) > MAGICAL_CONSTANT)
  1359.     {  /* split in half */
  1360.       DBL    uu1[3], uu2[3], vv1[3],vv2[3];
  1361.  
  1362.       rb_decasteljau_0_3( uu, uu1, uu2, 0.5, 0.5, 1);
  1363.       rb_decasteljau_0_3( vv, vv1, vv2, 0.5, 0.5, 1);
  1364.  
  1365.       return     trim_check_3(uu1, vv1, 0.5 * len) 
  1366.           + trim_check_3(uu2, vv2, 0.5 * len);
  1367.     }
  1368.   else
  1369.     {
  1370.       DBL uu1[3], uu2[3], uu3[3], vv1[3], vv2[3], vv3[3];
  1371.       if (min < 0.5)
  1372.     {
  1373.       rb_decasteljau_0_3( uu, uu1, uu2, min, 1.0-min, 1);
  1374.       rb_decasteljau_0_3( vv, vv1, vv2, min, 1.0-min, 1);
  1375.         
  1376.       rb_decasteljau_0_3( uu2, uu1, uu3, (max-min)/(1.0-min), (1.0-max)/(1.0-min), 1);
  1377.       rb_decasteljau_0_3( vv2, vv1, vv3, (max-min)/(1.0-min), (1.0-max)/(1.0-min), 1);
  1378.  
  1379.       return trim_check_3( uu1, vv1, len * (max-min) );
  1380.     }
  1381.       else
  1382.     {
  1383.       rb_decasteljau_0_3( uu, uu1, uu2, max, 1.0-max, 1);
  1384.       rb_decasteljau_0_3( vv, vv1, vv2, max, 1.0-max, 1);
  1385.         
  1386.       rb_decasteljau_0_3( uu1, uu2, uu3, min / max, 1.0 - min / max, 1);
  1387.       rb_decasteljau_0_3( vv1, vv2, vv3, min / max, 1.0 - min / max, 1);
  1388.  
  1389.       return trim_check_3( uu3, vv3, len * (max-min) );
  1390.     }
  1391.     }
  1392. }
  1393.  
  1394. static int  trim_check_4 (DBL *uu, DBL *vv,DBL len)
  1395. {
  1396.   int     vc,vcc, uc;
  1397.   DBL    min, max, tmp;
  1398.  
  1399.   vc= (vv[0]<0.0) + (vv[3]<0.0);
  1400.   vcc= vc + (vv[1]<0.0) + (vv[2]<0.0);
  1401.   if ((vcc==0)||(vcc==4)) return 0;
  1402.  
  1403.   uc= (uu[0]<0.0) + (uu[1]<0.0) + (uu[2]<0.0) + (uu[3]<0.0);
  1404.   if (uc==4) return 0;
  1405.   if (uc ==0) 
  1406.     if (vc==1)
  1407.       return 1;
  1408.     else  /* vc from {0,2} */
  1409.       return 0;
  1410.  
  1411.   /* seems that following check is useless. Moreover it seems that it's originator of some bad results
  1412.     if (len < Current_Shape->accuracy) return 0;
  1413.   */
  1414.  
  1415.   min= 2; max= -2;
  1416.  
  1417.   if ((vv[0]<0.0)^(vv[1]<0.0)) 
  1418.     TRIM_FIND_INTER( vv[0], vv[1],     0.0, 1.0/3.0);
  1419.   if ((vv[1]<0.0)^(vv[2]<0.0)) 
  1420.     TRIM_FIND_INTER( vv[1], vv[2], 1.0/3.0, 2.0/3.0);
  1421.   if ((vv[2]<0.0)^(vv[3]<0.0)) 
  1422.     TRIM_FIND_INTER( vv[2], vv[3], 2.0/3.0,     1.0);
  1423.  
  1424.   if ((vv[0]<0.0)^(vv[2]<0.0)) 
  1425.     TRIM_FIND_INTER( vv[0], vv[2],     0.0, 2.0/3.0);
  1426.   if ((vv[1]<0.0)^(vv[3]<0.0)) 
  1427.     TRIM_FIND_INTER( vv[1], vv[3], 1.0/3.0,     1.0);
  1428.  
  1429.   if (vc == 1)
  1430.     TRIM_FIND_INTER( vv[0], vv[3],     0.0,     1.0);
  1431.  
  1432.   if (min > max) return 0;
  1433.  
  1434.   if ((max-min) < Current_Shape->accuracy)
  1435.     { /* size of new part would be too small -=> stop */
  1436.       DBL uu1[4], uu2[4],vv1[4], vv2[4];
  1437.  
  1438.       rb_decasteljau_0_4( uu, uu1, uu2, min, 1.0-min, 1);
  1439.       rb_decasteljau_0_4( vv, vv1, vv2, min, 1.0-min, 1);
  1440.  
  1441.       if ((uu2[0] >= 0.0) && (fabs(vv2[0]) < Current_Shape-> accuracy)) return 1;
  1442.       else return 0;
  1443.     }
  1444.   if ((max-min) > MAGICAL_CONSTANT)
  1445.     {  /* split in half */
  1446.       DBL    uu1[4], uu2[4], vv1[4],vv2[4];
  1447.  
  1448.       rb_decasteljau_0_4( uu, uu1, uu2, 0.5, 0.5, 1);
  1449.       rb_decasteljau_0_4( vv, vv1, vv2, 0.5, 0.5, 1);
  1450.  
  1451.       return     trim_check_4(uu1, vv1, 0.5 * len) 
  1452.           + trim_check_4(uu2, vv2, 0.5 * len);
  1453.     }
  1454.   else
  1455.     {
  1456.       DBL uu1[4], uu2[4], uu3[4], vv1[4], vv2[4], vv3[4];
  1457.       if (min < 0.5)
  1458.     {
  1459.       rb_decasteljau_0_4( uu, uu1, uu2, min, 1.0-min, 1);
  1460.       rb_decasteljau_0_4( vv, vv1, vv2, min, 1.0-min, 1);
  1461.         
  1462.       rb_decasteljau_0_4( uu2, uu1, uu3, (max-min)/(1.0-min), (1.0-max)/(1.0-min), 1);
  1463.       rb_decasteljau_0_4( vv2, vv1, vv3, (max-min)/(1.0-min), (1.0-max)/(1.0-min), 1);
  1464.  
  1465.       return trim_check_4( uu1, vv1, len * (max-min) );
  1466.     }
  1467.       else
  1468.     {
  1469.       rb_decasteljau_0_4( uu, uu1, uu2, max, 1.0-max, 1);
  1470.       rb_decasteljau_0_4( vv, vv1, vv2, max, 1.0-max, 1);
  1471.         
  1472.       rb_decasteljau_0_4( uu1, uu2, uu3, min / max, 1.0 - min / max, 1);
  1473.       rb_decasteljau_0_4( vv1, vv2, vv3, min / max, 1.0 - min / max, 1);
  1474.  
  1475.       return trim_check_4( uu3, vv3, len * (max-min) );
  1476.     }
  1477.     }
  1478. }
  1479.  
  1480. #define CONTINUE_2    uu++; vv++; continue; 
  1481. #define CONTINUE_3    uu+=2; vv+=2;  continue; 
  1482. #define CONTINUE_4    uu+=3; vv+=3;  continue; 
  1483.  
  1484. static int trim_check (TRIM_SHAPE *trim, DBL u,DBL v)
  1485. {
  1486.   char         dir, *cc;
  1487.   DBL        up[TRIM_MAX_CONST], vp[TRIM_MAX_CONST];
  1488.   DBL        *uu,*vv, *ww;
  1489.  
  1490.   int        i, inter, tmpi;
  1491.   
  1492.   Current_Trim=    trim;
  1493.   inter= trim->ttype;
  1494.   if ((dir= trim_choose_direction(u,v)) < 0) return inter;
  1495.  
  1496.   trim_translate_points(dir, u,v, up, vp);
  1497.  
  1498.   uu= up;   vv= vp;   ww= trim->w;   cc= trim->ctypes;
  1499.  
  1500.   for (i=trim->num_of_parts; i; i--, cc++)
  1501.     switch ( (*cc) & 0x0f )
  1502.       {
  1503.         case 2:  /* 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1  */
  1504.        tmpi= (vv[0] < 0.0) + (vv[1] < 0.0);
  1505.        if ((tmpi & 1)==0) { CONTINUE_2 }
  1506.  
  1507.        tmpi= (uu[0] < 0.0) + (uu[1] < 0.0);
  1508.        if (tmpi == 2) { CONTINUE_2 }
  1509.        if (tmpi == 0) { inter++; CONTINUE_2 }
  1510.  
  1511.        if  ( 0.0 >= (
  1512.               uu[0] + ((uu[1]-uu[0]) * (vv[0])) / (vv[0] - vv[1])
  1513.                 )) { CONTINUE_2 }
  1514.        inter++;
  1515.        CONTINUE_2
  1516.         case 3: 
  1517.       if (((*cc) & 0xf0) == 0)
  1518.         {
  1519.           inter+= trim_check_3(uu, vv,1.0);
  1520.           CONTINUE_3;
  1521.         }
  1522.        else
  1523.          {
  1524.            DBL u1[3], v1[3];
  1525.            u1[0]= (*ww) * (*uu);    v1[0]= (*ww) * (*vv);    ww++; uu++; vv++;
  1526.            u1[1]= (*ww) * (*uu);    v1[1]= (*ww) * (*vv);    ww++; uu++; vv++;
  1527.            u1[2]= (*ww) * (*uu);    v1[2]= (*ww) * (*vv);    ww++;
  1528.            inter += trim_check_3(u1, v1,1.0);
  1529.            continue;
  1530.          }
  1531.         case 4: 
  1532.       if (((*cc) & 0xf0) == 0)
  1533.         {
  1534.           inter+= trim_check_4(uu, vv,1.0);
  1535.           CONTINUE_4;
  1536.         }
  1537.        else
  1538.          {
  1539.            DBL u1[4], v1[4];
  1540.            u1[0]= (*ww) * (*uu);    v1[0]= (*ww) * (*vv);    ww++; uu++; vv++;
  1541.            u1[1]= (*ww) * (*uu);    v1[1]= (*ww) * (*vv);    ww++; uu++; vv++;
  1542.            u1[2]= (*ww) * (*uu);    v1[2]= (*ww) * (*vv);    ww++; uu++; vv++;
  1543.            u1[3]= (*ww) * (*uu);    v1[3]= (*ww) * (*vv);    ww++;
  1544.            inter += trim_check_4(u1, v1,1.0);
  1545.            continue;
  1546.          }
  1547.         default:
  1548.            Error("This trim is not implemented yet");
  1549.       }
  1550.   return inter;
  1551. }
  1552.  
  1553. static int rbezier_check_trims(DBL u,DBL v)
  1554. {
  1555.   int     i;
  1556.  
  1557.  
  1558.   for (i=0; i < Current_Shape->num_of_trims; i++)
  1559.       if ( (trim_check((Current_Shape->trims)[i], u, v) & 1) == 0)
  1560.     return 0;
  1561.  
  1562.   return  1;
  1563. }
  1564.  
  1565. void Destroy_Trim( TRIM_SHAPE *ptr)
  1566. {
  1567.   (ptr->refs)--;
  1568.   if (ptr->refs == 0)
  1569.     {
  1570.       POV_FREE( ptr-> u );
  1571.       POV_FREE( ptr );
  1572.     }  
  1573. }    
  1574.  
  1575. TRIM_SHAPE *Copy_Trim (TRIM_SHAPE *trim)
  1576. {
  1577.   if (trim==NULL) return NULL;
  1578.   else
  1579.     {
  1580.       (trim-> refs)++;
  1581.       return trim;
  1582.     }
  1583. }
  1584.  
  1585. void trim_rotate_cp (DBL angle,DBL *up,DBL *vp,int num)
  1586. {
  1587.   DBL s,c, x,y;
  1588.  
  1589.   angle*= M_PI_180;
  1590.   s= sin(angle);
  1591.   c= cos(angle);
  1592.  
  1593.   for (; num; num--, up++, vp++)
  1594.     {
  1595.       x= *up;
  1596.       y= *vp;
  1597.       *up= x * c - y * s;
  1598.       *vp= x * s + y * c;
  1599.     }
  1600. }
  1601.  
  1602. void trim_scale_cp ( DBL *vec,DBL *up,DBL *vp,int num)
  1603. {
  1604.   for (; num; num--, up++, vp++)
  1605.     {    
  1606.       *up *= vec[0];
  1607.       *vp *= vec[1];
  1608.     }
  1609. }
  1610.  
  1611. void trim_translate_cp  ( DBL *vec,DBL *up,DBL *vp,int num)
  1612. {
  1613.   for (; num; num--, up++, vp++)
  1614.     {    
  1615.       *up += vec[0];
  1616.       *vp += vec[1];
  1617.     }
  1618. }
  1619.  
  1620. void trim_compute_bounds (TRIM_SHAPE *trim)
  1621. {
  1622.   int tmp;
  1623.   DBL *u,*v;
  1624.  
  1625.   trim->bounds[0][0]=trim->bounds[0][1]= trim-> u[0];
  1626.   trim->bounds[1][0]=trim->bounds[1][1]= trim->v[0];
  1627.   
  1628.   u= trim->u;
  1629.   v= trim->v;
  1630.  
  1631.   for (tmp=1; tmp < trim->num_of_cp; tmp++)
  1632.     {
  1633.       if (trim->bounds[0][0] > u[tmp]) trim->bounds[0][0]= u[tmp];
  1634.       if (trim->bounds[0][1] < u[tmp]) trim->bounds[0][1]= u[tmp];
  1635.       if (trim->bounds[1][0] > v[tmp]) trim->bounds[1][0]= v[tmp];
  1636.       if (trim->bounds[1][1] < v[tmp]) trim->bounds[1][1]= v[tmp];
  1637.     }
  1638. }
  1639.  
  1640.  
  1641. TRIM_SHAPE *trim_deep_copy( TRIM_SHAPE *trim)
  1642. {
  1643.   char        *b;
  1644.   int         len;
  1645.   TRIM_SHAPE     *New;
  1646.  
  1647.   New= (TRIM_SHAPE *) POV_MALLOC (sizeof(TRIM_SHAPE), "trimming shape");
  1648.  
  1649.   memcpy( New, trim, sizeof(TRIM_SHAPE));
  1650.  
  1651.   b= (char *) trim -> u;
  1652.   len= (trim->ctypes - b);
  1653.   len+= trim->num_of_parts;
  1654.  
  1655.   New-> u= (DBL *) POV_MALLOC ( len, "trimming shape");
  1656.   memcpy (New->u, trim->u, len );
  1657.  
  1658.   /* what a horrible adress arithmetic */
  1659.   New-> v= New->u + ((trim->v) - (trim->u));
  1660.   New-> w= New->u + ((trim->w) - (trim->u));
  1661.  
  1662.   New-> ctypes= (char *) New->u;
  1663.   New-> ctypes+= (trim->ctypes - b);
  1664.   New-> refs= 1;
  1665.  
  1666.   return New;
  1667. } /* End of really ugly function */
  1668. #else
  1669. static char dummy[2];
  1670. #endif /*#ifdef RBezierPatch */
  1671.