ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/sm_qtree.c
(Generate patch)

Comparing ray/src/hd/sm_qtree.c (file contents):
Revision 3.5 by gwlarson, Mon Sep 14 10:33:47 1998 UTC vs.
Revision 3.6 by gwlarson, Wed Sep 16 18:16:29 1998 UTC

# Line 174 | Line 174 | FVECT t0,t1,t2;
174        tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
175        intersect_vector_plane(pt,n,pd,NULL,i_pt);
176  
177 <      i = max_index(n);
177 >      i = max_index(n,NULL);
178        x = (i+1)%3;
179        y = (i+2)%3;
180        /* Calculate barycentric coordinates of i_pt */
# Line 374 | Line 374 | int n;
374                for(optr = QT_SET_PTR(os),i = QT_SET_CNT(os); i > 0; i--)
375                  {
376                    id = QT_SET_NEXT_ELEM(optr);
377 <                  qtTri_from_id(id,NULL,NULL,NULL,r0,r1,r2,NULL,NULL,NULL);
377 >                  qtTri_from_id(id,r0,r1,r2,NULL,NULL,NULL,NULL,NULL,NULL);
378                    found=qtAdd_tri(qtptr,q0,q1,q2,r0,r1,r2,id,n);
379   #ifdef DEBUG
380                    if(!found)
# Line 527 | Line 527 | double *tptr;
527   }
528  
529   int
530 – qtTrace_edge(qtptr,b,db,orig,dir,max_t,func,arg1,arg2)
531 –   QUADTREE *qtptr;
532 –   double b[3],db[3];
533 –   FVECT orig,dir;
534 –   double max_t;
535 –   int (*func)();
536 –   int *arg1,arg2;
537 – {
538 –
539 –    int i,found;
540 –    QUADTREE *child;
541 –    int nbr,next;
542 –    double t;
543 – #ifdef DEBUG_TEST_DRIVER
544 –    
545 –    FVECT a1,b1,c1;
546 –    int Pick_parent = Pick_cnt-1;
547 –    qtSubdivide_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
548 –                           Pick_v2[Pick_parent],a1,b1,c1);
549 –
550 – #endif
551 –    if(QT_IS_TREE(*qtptr))
552 –    {
553 –        /* Find the appropriate child and reset the coord */
554 –        i = bary_child(b);
555 –
556 –        QT_SET_FLAG(*qtptr);
557 –
558 –        for(;;)
559 –        {
560 –           child = QT_NTH_CHILD_PTR(*qtptr,i);
561 –
562 –           if(i != 3)
563 –             {
564 –
565 –               db[0] *= 2.0;db[1] *= 2.0; db[2] *= 2.0;
566 –              nbr = qtTrace_edge(child,b,db,orig,dir,max_t,func,arg1,arg2);
567 –               db[0] *= 0.5;db[1] *= 0.5; db[2] *= 0.5;
568 –             }
569 –           else
570 –             {
571 –               db[0] *=-2.0;db[1] *= -2.0; db[2] *= -2.0;
572 –              /* If the center cell- must flip direction signs */
573 –              nbr = qtTrace_edge(child,b,db,orig,dir,max_t,func,arg1,arg2);
574 –              db[0] *=-0.5;db[1] *= -0.5; db[2] *= -0.5;
575 –             }
576 –           if(nbr == QT_DONE)
577 –              return(nbr);
578 –
579 –           /* If in same block: traverse */
580 –           if(i==3)
581 –              next = nbr;
582 –             else
583 –                if(nbr == i)
584 –                   next = 3;
585 –             else
586 –               {
587 –                 /* reset the barycentric coordinates in the parents*/
588 –                 bary_parent(b,i);
589 –                 /* Else pop up to parent and traverse from there */
590 –                 return(nbr);
591 –               }
592 –             bary_from_child(b,i,next);
593 –             i = next;
594 –         }
595 –    }
596 –    else
597 –   {
598 – #ifdef DEBUG_TEST_DRIVER
599 –           qtNth_child_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
600 –                           Pick_v2[Pick_parent],a1,b1,c1,i,
601 –                           Pick_v0[Pick_cnt],Pick_v1[Pick_cnt],
602 –                           Pick_v2[Pick_cnt]);
603 –           Pick_cnt++;
604 – #endif
605 –
606 –       if(func(qtptr,arg1,arg2) == QT_DONE)
607 –          return(QT_DONE);
608 –
609 –       /* Advance to next node */
610 –       /* NOTE: Optimize: should only have to check 1/2 */
611 –       nbr = move_to_nbr(b,db[0],db[1],db[2],&t);
612 –
613 –       if(t >= max_t)
614 –          return(QT_DONE);
615 –       if(nbr != -1)
616 –       {
617 –         b[0] += t * db[0];
618 –         b[1] += t * db[1];
619 –         b[2] += t * db[2];
620 –         db[0] *= (1.0 - t);
621 –         db[1] *= (1.0 - t);
622 –         db[2] *= (1.0 - t);
623 –       }
624 –       return(nbr);
625 –   }
626 –    
627 – }
628 –
629 –
630 – int
530   qtTrace_ray(qtptr,b,db0,db1,db2,orig,dir,func,arg1,arg2)
531     QUADTREE *qtptr;
532     double b[3],db0,db1,db2;
# Line 746 | Line 645 | qtRoot_trace_ray(qtptr,q0,q1,q2,orig,dir,func,arg1,arg
645    intersect_vector_plane(c,n,pd,&t,c);
646  
647      /* map to 2d by dropping maximum magnitude component of normal */
648 <  i = max_index(n);
648 >  i = max_index(n,NULL);
649    x = (i+1)%3;
650    y = (i+2)%3;
651    /* Calculate barycentric coordinates of orig */
# Line 772 | Line 671 | qtRoot_trace_ray(qtptr,q0,q1,q2,orig,dir,func,arg1,arg
671      
672   }
673  
674 <
776 < int
777 < qtRoot_trace_edge(qtptr,q0,q1,q2,orig,dir,max_t,func,arg1,arg2)
674 > qtVisit_tri_interior(qtptr,q0,q1,q2,t0,t1,t2,n,func,arg1,arg2,arg3)
675     QUADTREE *qtptr;
676     FVECT q0,q1,q2;
780 –   FVECT orig,dir;
781 –   double max_t;
782 –   int (*func)();
783 –   int *arg1,arg2;
784 – {
785 –  int i,x,y,nbr;
786 –  QUADTREE *child;
787 –  FVECT n,c,i_pt,d;
788 –  double pd,b[3],db[3],t;
789 –    /* Determine if point lies within pyramid (and therefore
790 –       inside a spherical quadtree cell):GT_INTERIOR, on one of the
791 –       pyramid sides (and on cell edge):GT_EDGE(1,2 or 3),
792 –       or on pyramid vertex (and on cell vertex):GT_VERTEX(1,2, or 3).
793 –       For each triangle edge: compare the
794 –       point against the plane formed by the edge and the view center
795 –    */
796 –  i = point_in_stri(q0,q1,q2,orig);  
797 –    
798 –  /* Not in this triangle */
799 –  if(!i)
800 –     return(-1);
801 –  /* Project the origin onto the root node plane */
802 –
803 –  /* Find the intersection point of the origin */
804 –  tri_plane_equation(q0,q1,q2,n,&pd,FALSE);
805 –  intersect_vector_plane(orig,n,pd,NULL,i_pt);
806 –  /* project the dir as well */
807 –  VADD(c,orig,dir);
808 –  intersect_vector_plane(c,n,pd,&t,c);
809 –
810 –    /* map to 2d by dropping maximum magnitude component of normal */
811 –  i = max_index(n);
812 –  x = (i+1)%3;
813 –  y = (i+2)%3;
814 –  /* Calculate barycentric coordinates of orig */
815 –  bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b);
816 –  /* Calculate barycentric coordinates of dir */
817 –  bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],c[x],c[y],db);
818 –  if(t < 0.0)
819 –     VSUB(db,b,db);
820 –  else
821 –     VSUB(db,db,b);
822 –
823 –
824 – #ifdef DEBUG_TEST_DRIVER
825 –    VCOPY(Pick_v0[Pick_cnt],q0);
826 –    VCOPY(Pick_v1[Pick_cnt],q1);
827 –    VCOPY(Pick_v2[Pick_cnt],q2);
828 –    Pick_cnt++;
829 – #endif
830 –      /* trace the ray starting with this node */
831 –    nbr = qtTrace_edge(qtptr,b,db,orig,d,max_t,func,arg1,arg2);
832 –    return(nbr);
833 –    
834 – }
835 –
836 –
837 – qtVisit_tri_interior(qtptr,q0,q1,q2,t0,t1,t2,n,func,arg1,arg2)
838 –   QUADTREE *qtptr;
839 –   FVECT q0,q1,q2;
677     FVECT t0,t1,t2;
678     int n;
679     int (*func)();
680 <   int *arg1,arg2;
680 >   int *arg1,arg2,*arg3;
681   {
682      int i,found,test;
683      QUADTREE *child;
# Line 863 | Line 700 | tree_modified:
700                  child = QT_NTH_CHILD_PTR(*qtptr,i);
701                  qtNth_child_tri(q0,q1,q2,a,b,c,i,c0,c1,c2);
702                  qtVisit_tri_interior(child,c0,c1,c2,t0,t1,t2,n+1,
703 <                                     func,arg1,arg2);
703 >                                     func,arg1,arg2,arg3);
704              }
705          }
706      }
# Line 873 | Line 710 | tree_modified:
710        if(!QT_IS_EMPTY(*qtptr))
711        {
712           if(qtinset(*qtptr,arg2))
713 <           if(func(qtptr,q0,q1,q2,t0,t1,t2,n,arg1,arg2)==QT_MODIFIED)
713 >           if(func(qtptr,q0,q1,q2,t0,t1,t2,n,arg1,arg2,arg3)==QT_MODIFIED)
714               goto tree_modified;
715             else
716               return;
717        }
718        if(point_in_stri(t0,t1,t2,q0) )
719 <          if(func(qtptr,q0,q1,q2,t0,t1,t2,n,arg1,arg2)==QT_MODIFIED)
719 >          if(func(qtptr,q0,q1,q2,t0,t1,t2,n,arg1,arg2,arg3)==QT_MODIFIED)
720              goto tree_modified;
721      }
722   }
# Line 889 | Line 726 | tree_modified:
726  
727  
728  
729 + /* NOTE: SINCE DIR could be unit: then we could use integer math */
730   int
731 < qtVisit_tri_edges(qtptr,b,db,wptr,sfactor,func,arg1,arg2)
731 > qtVisit_tri_edges(qtptr,b,db0,db1,db2,
732 >                   db,wptr,t,sign,sfactor,func,arg1,arg2,arg3)
733     QUADTREE *qtptr;
734 <   double b[3],db[3][3];
734 >   double b[3],db0,db1,db2,db[3][3];
735     int *wptr;
736 +   double t[3];
737 +   int sign;
738     double sfactor;
739     int (*func)();
740 <   int *arg1,arg2;
740 >   int *arg1,arg2,*arg3;
741   {
742      int i,found;
743      QUADTREE *child;
744      int nbr,next,w;
745 <    double t;
745 >    double t_l,t_g;
746   #ifdef DEBUG_TEST_DRIVER
747      FVECT a1,b1,c1;
748      int Pick_parent = Pick_cnt-1;
749      qtSubdivide_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
750                             Pick_v2[Pick_parent],a1,b1,c1);
751   #endif
911 –
752      if(QT_IS_TREE(*qtptr))
753      {
754          /* Find the appropriate child and reset the coord */
# Line 922 | Line 762 | qtVisit_tri_edges(qtptr,b,db,wptr,sfactor,func,arg1,ar
762             child = QT_NTH_CHILD_PTR(*qtptr,i);
763  
764             if(i != 3)
765 <             {
766 <
767 <               db[w][0] *= 2.0;db[w][1] *= 2.0; db[w][2] *= 2.0;
928 <               nbr = qtVisit_tri_edges(child,b,db,wptr,sfactor*2.0,
929 <                                      func,arg1,arg2);
930 <               w = *wptr;
931 <               db[w][0] *= 0.5;db[w][1] *= 0.5; db[w][2] *= 0.5;
932 <             }
765 >               nbr = qtVisit_tri_edges(child,b,db0,db1,db2,
766 >                                        db,wptr,t,sign,
767 >                                        sfactor*2.0,func,arg1,arg2,arg3);
768             else
769 <             {
770 <               db[w][0] *=-2.0;db[w][1] *= -2.0; db[w][2] *= -2.0;
771 <              /* If the center cell- must flip direction signs */
772 <               nbr = qtVisit_tri_edges(child,b,db,wptr,sfactor*(-2.0),
938 <                                 func,arg1,arg2);
939 <               w = *wptr;
940 <               db[w][0] *=-0.5;db[w][1] *= -0.5; db[w][2] *= -0.5;
941 <             }
942 <           if(nbr == QT_DONE)
943 <              return(nbr);
769 >             /* If the center cell- must flip direction signs */
770 >             nbr = qtVisit_tri_edges(child,b,-db0,-db1,-db2,
771 >                                     db,wptr,t,1-sign,
772 >                                     sfactor*2.0,func,arg1,arg2,arg3);
773  
774 +           if(nbr == QT_DONE)
775 +             return(nbr);
776 +           if(*wptr != w)
777 +           {
778 +             w = *wptr;
779 +             db0 = db[w][0];db1 = db[w][1];db2 = db[w][2];
780 +             if(sign)
781 +               {  db0 *= -1.0;db1 *= -1.0; db2 *= -1.0;}
782 +           }
783             /* If in same block: traverse */
784             if(i==3)
785                next = nbr;
# Line 955 | Line 793 | qtVisit_tri_edges(qtptr,b,db,wptr,sfactor,func,arg1,ar
793                   /* Else pop up to parent and traverse from there */
794                   return(nbr);
795                 }
796 <             bary_from_child(b,i,next);
797 <             i = next;
798 <         }
796 >           bary_from_child(b,i,next);
797 >           i = next;
798 >        }
799      }
800      else
801     {
# Line 968 | Line 806 | qtVisit_tri_edges(qtptr,b,db,wptr,sfactor,func,arg1,ar
806             Pick_cnt++;
807   #endif
808  
809 <       if(func(qtptr,arg1,arg2) == QT_DONE)
809 >       if(func(qtptr,arg1,arg2,arg3) == QT_DONE)
810            return(QT_DONE);
811  
812         /* Advance to next node */
975 –       /* NOTE: Optimize: should only have to check 1/2 */
813         w = *wptr;
814         while(1)
815         {
816 <         nbr = move_to_nbr(b,db[w][0],db[w][1],db[w][2],&t);
816 >         nbr = move_to_nbr(b,db0,db1,db2,&t_l);
817  
818 <         if(t >= 1.0)
819 <         {
818 >         t_g = t_l/sfactor;
819 > #ifdef DEBUG
820 >         if(t[w] <= 0.0)
821 >           eputs("qtVisit_tri_edges():negative t\n");
822 > #endif
823 >         if(t_g >= t[w])
824 >         {
825             if(w == 2)
826               return(QT_DONE);
827 <           b[0] += db[w][0];
828 <           b[1] += db[w][1];
829 <           b[2] += db[w][2];
827 >
828 >           b[0] += (t[w])*sfactor*db0;
829 >           b[1] += (t[w])*sfactor*db1;
830 >           b[2] += (t[w])*sfactor*db2;
831             w++;
832 <           db[w][0] *= sfactor;
833 <           db[w][1] *= sfactor;
834 <           db[w][2] *= sfactor;
832 >           db0 = db[w][0];
833 >           db1 = db[w][1];
834 >           db2 = db[w][2];
835 >           if(sign)
836 >          {  db0 *= -1.0;db1 *= -1.0; db2 *= -1.0;}
837           }
838         else
839           if(nbr != INVALID)
840           {
841 <           b[0] += t * db[w][0];
842 <           b[1] += t * db[w][1];
843 <           b[2] += t * db[w][2];
844 <           db[w][0] *= (1.0 - t);
845 <           db[w][1] *= (1.0 - t);
1001 <           db[w][2] *= (1.0 - t);
841 >           b[0] += t_l * db0;
842 >           b[1] += t_l * db1;
843 >           b[2] += t_l * db2;
844 >
845 >           t[w] -= t_g;
846             *wptr = w;
847             return(nbr);
848           }
# Line 1011 | Line 855 | qtVisit_tri_edges(qtptr,b,db,wptr,sfactor,func,arg1,ar
855  
856  
857   int
858 < qtRoot_visit_tri_edges(qtptr,q0,q1,q2,tri,dir,wptr,func,arg1,arg2)
858 > qtRoot_visit_tri_edges(qtptr,q0,q1,q2,tri,i_pt,wptr,func,arg1,arg2,arg3)
859     QUADTREE *qtptr;
860     FVECT q0,q1,q2;
861 <   FVECT tri[3],dir[3];
861 >   FVECT tri[3],i_pt;
862     int *wptr;
863     int (*func)();
864 <   int *arg1,arg2;
864 >   int *arg1,arg2,*arg3;
865   {
866 <  int i,x,y,nbr,w;
866 >  int x,y,z,nbr,w,i,j;
867    QUADTREE *child;
868 <  FVECT n,c,i_pt,d;
869 <  double pd,b[3][3],db[3][3],t;
868 >  FVECT n,c,d,v[3];
869 >  double pd,b[4][3],db[3][3],et[3],t[3],exit_pt;
870      
871    w = *wptr;
872  
# Line 1031 | Line 875 | qtRoot_visit_tri_edges(qtptr,q0,q1,q2,tri,dir,wptr,fun
875    /* Find the intersection point of the origin */
876    tri_plane_equation(q0,q1,q2,n,&pd,FALSE);
877    /* map to 2d by dropping maximum magnitude component of normal */
878 <  i = max_index(n);
879 <  x = (i+1)%3;
880 <  y = (i+2)%3;
878 >  z = max_index(n,NULL);
879 >  x = (z+1)%3;
880 >  y = (z+2)%3;
881    /* Calculate barycentric coordinates for current vertex */
882 <  
1039 <  for(i=0;i < 3; i++)
882 >  if(w != -1)
883    {
884 <    /* If processing 3rd edge-dont need info for t1 */
885 <    if(i==1 && w==2)
886 <      continue;
1044 <    /* project the dir as well */
1045 <    intersect_vector_plane(tri[i],n,pd,NULL,i_pt);
1046 <    bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b[i]);  
1047 <    VADD(c,tri[i],dir[i]);
1048 <    intersect_vector_plane(c,n,pd,&t,c);
1049 <    /* Calculate barycentric coordinates of dir */
1050 <    bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],c[x],c[y],db[i]);
1051 <    if(t < 0.0)
1052 <      VSUB(db[i],b[i],db[i]);
1053 <    else
1054 <      VSUB(db[i],db[i],b[i]);
884 >    bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b[3]);  
885 >    intersect_vector_plane(tri[w],n,pd,&(et[w]),v[w]);
886 >    bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[w][x],v[w][y],b[w]);  
887    }
888 < #ifdef DEBUG_TEST_DRIVER
889 <    VCOPY(Pick_v0[Pick_cnt],q0);
890 <    VCOPY(Pick_v1[Pick_cnt],q1);
891 <    VCOPY(Pick_v2[Pick_cnt],q2);
892 <    Pick_cnt++;
893 < #endif
894 <      /* trace the ray starting with this node */
895 <    nbr = qtVisit_tri_edges(qtptr,b[w],db,wptr,1.0,func,arg1,arg2);
888 >  else
889 >  /* Just starting: b[0] is the origin point: guaranteed to be valid b*/
890 >  {
891 >    w = 0;
892 >    intersect_vector_plane(tri[0],n,pd,&(et[0]),v[0]);
893 >    bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[0][x],v[0][y],b[0]);  
894 >    VCOPY(b[3],b[0]);
895 >  }
896 >
897 >
898 >  j = (w+1)%3;
899 >  intersect_vector_plane(tri[j],n,pd,&(et[j]),v[j]);
900 >  bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[j][x],v[j][y],b[j]);  
901 >  if(et[j] < 0.0)
902 >  {
903 >      VSUB(db[w],b[3],b[j]);
904 >      t[w] = FHUGE;
905 >  }
906 >  else
907 > {
908 >   /* NOTE: for stability: do not increment with ipt- use full dir and
909 >      calculate t: but for wrap around case: could get same problem?
910 >      */
911 >   VSUB(db[w],b[j],b[3]);
912 >   t[w] = 1.0;
913 >   move_to_nbr(b[3],db[w][0],db[w][1],db[w][2],&exit_pt);
914 >   if(exit_pt >= 1.0)
915 >   {
916 >     for(;j < 3;j++)
917 >      {
918 >          i = (j+1)%3;
919 >          if(i!= w)
920 >            {
921 >              intersect_vector_plane(tri[i],n,pd,&(et[i]),v[i]);
922 >              bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[i][x],
923 >                     v[i][y],b[i]);
924 >            }
925 >          if(et[i] < 0.0)
926 >             {
927 >                 VSUB(db[j],b[j],b[i]);
928 >                 t[j] = FHUGE;
929 >                 break;
930 >             }
931 >          else
932 >             {
933 >                 VSUB(db[j],b[i],b[j]);
934 >                 t[j] = 1.0;
935 >             }
936 >          move_to_nbr(b[j],db[j][0],db[j][1],db[j][2],&exit_pt);
937 >          if(exit_pt < 1.0)
938 >             break;
939 >      }
940 >   }
941 > }
942 >  *wptr = w;
943 >  /* trace the ray starting with this node */
944 >    nbr = qtVisit_tri_edges(qtptr,b[3],db[w][0],db[w][1],db[w][2],
945 >                             db,wptr,t,0,1.0,func,arg1,arg2,arg3);
946 >    if(nbr != INVALID && nbr != QT_DONE)
947 >    {
948 >      i_pt[x] = b[3][0]*q0[x] + b[3][1]*q1[x] + b[3][2]*q2[x];
949 >      i_pt[y] = b[3][0]*q0[y] + b[3][1]*q1[y] + b[3][2]*q2[y];
950 >      i_pt[z] = (-n[x]*i_pt[x] - n[y]*i_pt[y] -pd)/n[z];
951 >    }
952      return(nbr);
953      
954   }
955  
956 <
957 <
958 <
956 > int
957 > move_to_nbri(b,db0,db1,db2,tptr)
958 > BCOORD b[3];
959 > BDIR db0,db1,db2;
960 > TINT *tptr;
961 > {
962 >  TINT t,dt;
963 >  BCOORD bc;
964 >  int nbr;
965 >  
966 >  nbr = -1;
967 >  /* Advance to next node */
968 >  if(db0 < 0)
969 >   {
970 >     bc = MAXBCOORD*b[0];
971 >     t = bc/-db0;
972 >     nbr = 0;
973 >   }
974 >  else
975 >    t = HUGET;
976 >  if(db1 < 0)
977 >  {
978 >     bc = MAXBCOORD*b[1];
979 >     dt = bc/-db1;
980 >    if( dt < t)
981 >    {
982 >      t = dt;
983 >      nbr = 1;
984 >    }
985 >  }
986 >  if(db2 < 0)
987 >  {
988 >     bc = MAXBCOORD*b[2];
989 >     dt = bc/-db2;
990 >       if( dt < t)
991 >      {
992 >        t = dt;
993 >        nbr = 2;
994 >      }
995 >    }
996 >  *tptr = t;
997 >  return(nbr);
998 > }
999   /* NOTE: SINCE DIR could be unit: then we could use integer math */
1000   int
1001 < qtVisit_tri_edges2(qtptr,b,db0,db1,db2,
1002 <                   db,wptr,t,sign,sfactor,func,arg1,arg2)
1001 > qtVisit_tri_edgesi(qtptr,b,db0,db1,db2,
1002 >                   db,wptr,t,sign,sfactor,func,arg1,arg2,arg3)
1003     QUADTREE *qtptr;
1004 <   double b[3],db0,db1,db2,db[3][3];
1004 >   BCOORD b[3];
1005 >   BDIR db0,db1,db2,db[3][3];
1006     int *wptr;
1007 <   double t[3];
1007 >   TINT t[3];
1008     int sign;
1009 <   double sfactor;
1009 >   int sfactor;
1010     int (*func)();
1011 <   int *arg1,arg2;
1011 >   int *arg1,arg2,*arg3;
1012   {
1013      int i,found;
1014      QUADTREE *child;
1015      int nbr,next,w;
1016 <    double t_l,t_g;
1016 >    TINT t_g,t_l,t_i;
1017   #ifdef DEBUG_TEST_DRIVER
1018      FVECT a1,b1,c1;
1019      int Pick_parent = Pick_cnt-1;
# Line 1094 | Line 1023 | qtVisit_tri_edges2(qtptr,b,db0,db1,db2,
1023      if(QT_IS_TREE(*qtptr))
1024      {
1025          /* Find the appropriate child and reset the coord */
1026 <        i = bary_child(b);
1026 >        i = baryi_child(b);
1027  
1028          QT_SET_FLAG(*qtptr);
1029  
# Line 1104 | Line 1033 | qtVisit_tri_edges2(qtptr,b,db0,db1,db2,
1033             child = QT_NTH_CHILD_PTR(*qtptr,i);
1034  
1035             if(i != 3)
1036 <               nbr = qtVisit_tri_edges2(child,b,db0,db1,db2,
1036 >               nbr = qtVisit_tri_edgesi(child,b,db0,db1,db2,
1037                                          db,wptr,t,sign,
1038 <                                        sfactor*2.0,func,arg1,arg2);
1038 >                                        sfactor+1,func,arg1,arg2,arg3);
1039             else
1040               /* If the center cell- must flip direction signs */
1041 <             nbr = qtVisit_tri_edges2(child,b,-db0,-db1,-db2,
1041 >             nbr = qtVisit_tri_edgesi(child,b,-db0,-db1,-db2,
1042                                       db,wptr,t,1-sign,
1043 <                                     sfactor*2.0,func,arg1,arg2);
1043 >                                     sfactor+1,func,arg1,arg2,arg3);
1044  
1045             if(nbr == QT_DONE)
1046               return(nbr);
# Line 1120 | Line 1049 | qtVisit_tri_edges2(qtptr,b,db0,db1,db2,
1049               w = *wptr;
1050               db0 = db[w][0];db1 = db[w][1];db2 = db[w][2];
1051               if(sign)
1052 <               {  db0 *= -1.0;db1 *= -1.0; db2 *= -1.0;}
1052 >               {  db0 *= -1;db1 *= -1; db2 *= -1;}
1053             }
1054             /* If in same block: traverse */
1055             if(i==3)
# Line 1131 | Line 1060 | qtVisit_tri_edges2(qtptr,b,db0,db1,db2,
1060               else
1061                 {
1062                   /* reset the barycentric coordinates in the parents*/
1063 <                 bary_parent(b,i);
1063 >                 baryi_parent(b,i);
1064                   /* Else pop up to parent and traverse from there */
1065                   return(nbr);
1066                 }
1067 <           bary_from_child(b,i,next);
1067 >           baryi_from_child(b,i,next);
1068             i = next;
1069          }
1070      }
# Line 1148 | Line 1077 | qtVisit_tri_edges2(qtptr,b,db0,db1,db2,
1077             Pick_cnt++;
1078   #endif
1079  
1080 <       if(func(qtptr,arg1,arg2) == QT_DONE)
1080 >       if(func(qtptr,arg1,arg2,arg3) == QT_DONE)
1081            return(QT_DONE);
1082  
1083         /* Advance to next node */
1084         w = *wptr;
1085         while(1)
1086         {
1087 <         nbr = move_to_nbr(b,db0,db1,db2,&t_l);
1087 >           nbr = move_to_nbri(b,db0,db1,db2,&t_i);
1088  
1089 <         t_g = t_l/sfactor;
1090 < #ifdef DEBUG
1162 <         if(t[w] <= 0.0)
1163 <           eputs("qtVisit_tri_edges2():negative t\n");
1164 < #endif
1089 >         t_g = t_i >> sfactor;
1090 >                
1091           if(t_g >= t[w])
1092           {
1093             if(w == 2)
1094               return(QT_DONE);
1095  
1096 <           b[0] += (t[w])*sfactor*db0;
1097 <           b[1] += (t[w])*sfactor*db1;
1098 <           b[2] += (t[w])*sfactor*db2;
1096 >           /* The edge fits in the cell- therefore the result of shifting db by
1097 >              sfactor is guaranteed to be less than MAXBCOORD
1098 >            */
1099 >           /* Caution: (t[w]*db) must occur before divide by MAXBCOORD
1100 >              since t[w] will always be < MAXBCOORD
1101 >           */
1102 >           b[0] = b[0] + (t[w]*db0*(1 << sfactor))/MAXBCOORD;
1103 >           b[1] = b[1] + (t[w]*db1*(1 << sfactor))/MAXBCOORD;
1104 >           b[2] = b[2] + (t[w]*db2*(1 << sfactor))/MAXBCOORD;
1105             w++;
1106 <           db0 = db[w][0];
1175 <           db1 = db[w][1];
1176 <           db2 = db[w][2];
1106 >           db0 = db[w][0]; db1 = db[w][1]; db2 = db[w][2];
1107             if(sign)
1108 <          {  db0 *= -1.0;db1 *= -1.0; db2 *= -1.0;}
1108 >           {  db0 *= -1;db1 *= -1; db2 *= -1;}
1109           }
1110         else
1111           if(nbr != INVALID)
1112           {
1113 <           b[0] += t_l * db0;
1114 <           b[1] += t_l * db1;
1115 <           b[2] += t_l * db2;
1113 >           /* Caution: (t_i*db) must occur before divide by MAXBCOORD
1114 >              since t_i will always be < MAXBCOORD
1115 >           */
1116 >           b[0] = b[0] + (t_i *db0) / MAXBCOORD;
1117 >           b[1] = b[1] + (t_i *db1) / MAXBCOORD;
1118 >           b[2] = b[2] + (t_i *db2) / MAXBCOORD;
1119  
1120             t[w] -= t_g;
1121             *wptr = w;
# Line 1191 | Line 1124 | qtVisit_tri_edges2(qtptr,b,db0,db1,db2,
1124           else
1125             return(INVALID);
1126         }
1127 <   }
1195 <    
1127 >   }    
1128   }
1129  
1130  
1131   int
1132 < qtRoot_visit_tri_edges2(qtptr,q0,q1,q2,tri,i_pt,wptr,func,arg1,arg2)
1132 > qtRoot_visit_tri_edgesi(qtptr,q0,q1,q2,tri,i_pt,wptr,func,arg1,arg2,arg3)
1133     QUADTREE *qtptr;
1134     FVECT q0,q1,q2;
1135     FVECT tri[3],i_pt;
1136     int *wptr;
1137     int (*func)();
1138 <   int *arg1,arg2;
1138 >   int *arg1,arg2,*arg3;
1139   {
1140    int x,y,z,nbr,w,i,j;
1141    QUADTREE *child;
1142    FVECT n,c,d,v[3];
1143 <  double pd,b[4][3],db[3][3],et[3],t[3],exit_pt;
1144 <    
1143 >  double pd,b[4][3],db[3][3],et[3],exit_pt;
1144 >  BCOORD bi[3];
1145 >  TINT t[3];
1146 >  BDIR dbi[3][3];
1147    w = *wptr;
1148  
1149    /* Project the origin onto the root node plane */
1150  
1151 +  t[0] = t[1] = t[2] = 0;
1152    /* Find the intersection point of the origin */
1153    tri_plane_equation(q0,q1,q2,n,&pd,FALSE);
1154    /* map to 2d by dropping maximum magnitude component of normal */
1155 <  z = max_index(n);
1155 >  z = max_index(n,NULL);
1156    x = (z+1)%3;
1157    y = (z+2)%3;
1158    /* Calculate barycentric coordinates for current vertex */
# Line 1243 | Line 1178 | qtRoot_visit_tri_edges2(qtptr,q0,q1,q2,tri,i_pt,wptr,f
1178    if(et[j] < 0.0)
1179    {
1180        VSUB(db[w],b[3],b[j]);
1181 <      t[w] = FHUGE;
1181 >      t[w] = HUGET;
1182    }
1183    else
1184   {
# Line 1251 | Line 1186 | qtRoot_visit_tri_edges2(qtptr,q0,q1,q2,tri,i_pt,wptr,f
1186        calculate t: but for wrap around case: could get same problem?
1187        */
1188     VSUB(db[w],b[j],b[3]);
1189 <   t[w] = 1.0;
1190 <   move_to_nbr(b[3],db[w][0],db[w][1],db[w][2],&exit_pt);
1191 <   if(exit_pt >= 1.0)
1189 >   /* Check if the point is out side of the triangle: if so t[w] =HUGET */
1190 >   if((fabs(b[j][0])>(FTINY+1.0)) ||(fabs(b[j][1])>(FTINY+1.0)) ||
1191 >      (fabs(b[j][2])>(FTINY+1.0)))
1192 >      t[w] = HUGET;
1193 >   else
1194     {
1195 <     for(;j < 3;j++)
1196 <      {
1197 <          i = (j+1)%3;
1198 <          if(i!= w)
1199 <            {
1200 <              intersect_vector_plane(tri[i],n,pd,&(et[i]),v[i]);
1264 <              bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[i][x],
1265 <                     v[i][y],b[i]);
1266 <            }
1267 <          if(et[i] < 0.0)
1195 >       /* The next point is in the triangle- so db values will be in valid
1196 >          range and t[w]= MAXT
1197 >        */  
1198 >       t[w] = MAXT;
1199 >       if(j != 0)
1200 >          for(;j < 3;j++)
1201               {
1202 <                 VSUB(db[j],b[j],b[i]);
1203 <                 t[j] = FHUGE;
1204 <                 break;
1202 >                 i = (j+1)%3;
1203 >                 if(i!= w)
1204 >                 {
1205 >                     intersect_vector_plane(tri[i],n,pd,&(et[i]),v[i]);
1206 >                     bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[i][x],
1207 >                            v[i][y],b[i]);
1208 >                 }
1209 >                 if(et[i] < 0.0)
1210 >                 {
1211 >                     VSUB(db[j],b[j],b[i]);
1212 >                     t[j] = HUGET;
1213 >                     break;
1214 >                 }
1215 >                 else
1216 >                 {
1217 >                     VSUB(db[j],b[i],b[j]);
1218 >                     if((fabs(b[j][0])>(FTINY+1.0)) ||
1219 >                        (fabs(b[j][1])>(FTINY+1.0)) ||
1220 >                        (fabs(b[j][2])>(FTINY+1.0)))
1221 >                        {
1222 >                            t[j] = HUGET;
1223 >                            break;
1224 >                        }
1225 >                     else
1226 >                        t[j] = MAXT;
1227 >                 }
1228               }
1273 –          else
1274 –             {
1275 –                 VSUB(db[j],b[i],b[j]);
1276 –                 t[j] = 1.0;
1277 –             }
1278 –          move_to_nbr(b[j],db[j][0],db[j][1],db[j][2],&exit_pt);
1279 –          if(exit_pt < 1.0)
1280 –             break;
1281 –      }
1229     }
1230 < }
1230 > }
1231    *wptr = w;
1232 +  bary_dtol(b[3],db,bi,dbi,t);
1233 +
1234    /* trace the ray starting with this node */
1235 <    nbr = qtVisit_tri_edges2(qtptr,b[3],db[w][0],db[w][1],db[w][2],
1236 <                             db,wptr,t,0,1.0,func,arg1,arg2);
1235 >    nbr = qtVisit_tri_edgesi(qtptr,bi,dbi[w][0],dbi[w][1],dbi[w][2],
1236 >                             dbi,wptr,t,0,0,func,arg1,arg2,arg3);
1237      if(nbr != INVALID && nbr != QT_DONE)
1238      {
1239 <      i_pt[x] = b[3][0]*q0[x] + b[3][1]*q1[x] + b[3][2]*q2[x];
1240 <      i_pt[y] = b[3][0]*q0[y] + b[3][1]*q1[y] + b[3][2]*q2[y];
1241 <      i_pt[z] = (-n[x]*i_pt[x] - n[y]*i_pt[y] -pd)/n[z];
1239 >        b[3][0] = (double)bi[0]/(double)MAXBCOORD;
1240 >        b[3][1] = (double)bi[1]/(double)MAXBCOORD;
1241 >        b[3][2] = (double)bi[2]/(double)MAXBCOORD;
1242 >        i_pt[x] = b[3][0]*q0[x] + b[3][1]*q1[x] + b[3][2]*q2[x];
1243 >        i_pt[y] = b[3][0]*q0[y] + b[3][1]*q1[y] + b[3][2]*q2[y];
1244 >        i_pt[z] = (-n[x]*i_pt[x] - n[y]*i_pt[y] -pd)/n[z];
1245      }
1246      return(nbr);
1247      
1248   }
1297 –
1298 –
1249  
1250  
1251  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines