---

gunu_intersect.cpp

Go to the documentation of this file.

/* LIBGUL - Geometry Utility Library
 * Copyright (C) 1998-1999 Norbert Irmer
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */
 
#include "stdafx.h"

// #if !defined(_MSC_VER) || defined(NEED_NEW_H)
//  #include <new.h>
// #endif

#include <iostream>

#include "gul_types.h"
#include "gul_float.h"
#include "gul_vector.h" 
#include "guge_normalize.h"
#include "guge_linear.h"
#include "gunu_basics.h"
#include "gunu_divide.h"
#include "gunu_intersect.h"
#include "guar_exact.h"
#include "gugr_basics.h"
#include "guar_intersect.h"
#include "guma_random.h"
#include "guma_rkdtree.h"
#include "gunu_project.h"
#include "gul_io.h"

namespace gunu {

using guge::IsLinear;
using guge::CalcBoundingBoxVerts;
using guar::rational;
using gul::rtr;
using gul::rel_equal;
using gugr::coord2int;
using gul::point2;
using gul::point;
using gul::rel_equal;
using gul::triangle;
using gul::triangle2;
using guar::IntersectTrianglesUV;
using gul::Max;
using gul::Min;
using gul::test;
using gul::List;
using gul::ListNode;
using gugr::cnv2coord_rnd;
using gul::kdpoint;
using guma::kdnnrec;
using guma::kdrec;
using guma::rkdtree;
using guma::random_des;
using gul::line;
using std::cout;

/*--------------------------------------------------------------------
  Check if a NURBS surface can be approximated with two triangles through
  its four corners. If not the surface is subdivided into four parts
----------------------------------------------------------------------*/

template< class T, class HP >
void LinearizeOrDivide( TessInfo<T,HP> *A, const T tol, 
                        bool need_bbox )
{
  int nu,nv,i,pu,pv,j;
  Ptr< T > U,V;
  Ptr< Ptr< HP > > Pw;
  point<T> P00, P01, P10, P11, SP0, SP1, SP2, SP3, CP, minP, maxP;
  bool subdiv;
  surface< T, HP > S[4], surf;
  Ptr< HP > buf;
  gul::Ptr< int > test;
 
  nu = A->org->nu;
  pu = A->org->pu;
  nv = A->org->nv;
  pv = A->org->pv;
  U = A->org->U;
  V = A->org->V;
  Pw = A->org->Pw;

  P00 = A->org->P00;
  P01 = A->org->P01;
  P10 = A->org->P10;
  P11 = A->org->P11;
  
  buf.reserve_place( reserve_stack(hpoint<T>,nv+1), nv+1 );
  
                          /*** Pruefen, ob Randkurven linear approximierbar: ***/
  subdiv = false;
                                                   /* --- obere Randkurve: --- */
  if( !guge::IsLinear< T,HP,point<T> >( nu, Pw[0], tol ) ) 
  {
    subdiv = true; 
    CurvePoint< T, HP, point<T> >( (T)0.5, nu, pu, U, Pw[0], &SP0 );
  }
  else
  {
    SP0 = (T)0.5 * (P00 + P01);
  }
                                                  /* --- untere Randkurve: --- */
  if( !guge::IsLinear<T,HP,point<T> >( nu, Pw[nv], tol ) ) 
  {
    subdiv = true;
    CurvePoint< T, HP, point<T> >( (T)0.5, nu, pu, U, Pw[nv], &SP1 );  
  }
  else
  {
    SP1 = (T)0.5 * (P10 + P11);
  }
                                                   /* --- linke Randkurve: --- */
  for( i = 0; i <= nv; i++ )
    buf[i] = Pw[i][0];
  if( !guge::IsLinear<T,HP,point<T> >( nv, buf, tol ) ) 
  {
    subdiv = true;
    CurvePoint< T, HP, point<T> >( (T)0.5, nv, pv, V, buf, &SP2 );
  }
  else
  {
    SP2 = (T)0.5 * (P00 + P10);
  }
                                                  /* --- rechte Randkurve: --- */
  for( i = 0; i <= nv; i++ )
    buf[i] = Pw[i][nu];
  if( !guge::IsLinear<T,HP,point<T> >( nv, buf, tol ) ) 
  {
    subdiv = true;
    CurvePoint< T, HP, point<T> >( (T)0.5, nv, pv, V, buf, &SP3 );
  }
  else
  {
    SP3 = (T)0.5 * (P01 + P11);
  }  

  #if 0
  /* not necessary if the tesselation tolerance is fine enough and the
     surface doesn't intersect itself, i.e. in normal cases 
  */
  
  // check if iso curves are linear enough 
  if( !subdiv )
  { 
    for( i = 1; i < nv; i++ )
    {
      if( !guge::IsLinear< T,HP,point<T> >( nu, Pw[i], tol, t ) ) 
      {
        subdiv = true;
        break;
      }
    }
  }
  if( !subdiv )
  { 
    for( j = 1; j < nu; j++ )
    {
      for( i = 0; i <= nv; i++ )
        buf[i] = Pw[i][j];
      if( !guge::IsLinear<T,HP,point<T> >( nv, buf, tol, t ) ) 
      {
        subdiv = true;
        break;
      }
    }
  }
  #endif
                          /* --- Pruefen ob Flaeche linear approxmierbar --- */
  if( !subdiv )
  {
    if( !guge::IsPlanar<T,HP>( nu, nv, Pw, tol ) )
      subdiv = true;
  }
                /*** Flaeche nicht linear approximierbar, also unterteilen: ***/
  if( subdiv )
  {
    surf.pu = pu;
    surf.pv = pv;
    surf.net.nu = nu;
    surf.net.nv = nv;  
    surf.net.Pw = Pw;
    surf.knu.m = nu+pu+1;
    surf.knu.U = U;
    surf.knv.m = nv+pv+1;
    surf.knv.U = V;
   
    SplitSurface< T, HP >( &surf, 0.5, 0.5, &S[0], &S[1], &S[2], &S[3] );

    A->divided = 1;
    delete A->org; /* Speicher fuer ungeteilte Flaeche freigeben */
    A->org = 0;  
    for( i = 0; i < 4; i++ ) A->sub[i] = new TessInfo<T,HP>();

    for( i = 0; i < 4; i++ )
    {
      if( need_bbox )
      {
        guge::CalcBoundingBoxVerts< HP, point<T> >( 
                                   S[i].net.nu+1, S[i].net.Pw[0], minP, maxP );
        for( j = 1; j < S[i].net.nv+1; j++ )
          guge::UpdateBoundingBoxVerts< HP, point<T> >( 
                                   S[i].net.nu+1, S[i].net.Pw[j], minP, maxP );
        A->sub[i]->x1 = minP.x; A->sub[i]->x2 = maxP.x; 
        A->sub[i]->y1 = minP.y; A->sub[i]->y2 = maxP.y;
        A->sub[i]->z1 = minP.z; A->sub[i]->z2 = maxP.z;
      }
      A->sub[i]->org->nu = S[i].net.nu;
      A->sub[i]->org->pu = S[i].pu;
      A->sub[i]->org->nv = S[i].net.nv;
      A->sub[i]->org->pv = S[i].pv;
      A->sub[i]->org->U = S[i].knu.U;
      A->sub[i]->org->V = S[i].knv.U;
      A->sub[i]->org->Pw = S[i].net.Pw;
    }
    CP = euclid( A->sub[0]->org->Pw[ A->sub[0]->org->nv ][ A->sub[0]->org->nu ] );

    A->sub[0]->org->P00 = P00;
    A->sub[0]->org->P01 = SP0;
    A->sub[0]->org->P10 = SP2;
    A->sub[0]->org->P11 = CP;
    A->sub[0]->u1 = A->u1;
    A->sub[0]->u2 = (A->u1 + A->u2) / (T)2.0;
    A->sub[0]->v1 = A->v1;
    A->sub[0]->v2 = (A->v1 + A->v2) / (T)2.0;    
    
    A->sub[1]->org->P00 = SP0;
    A->sub[1]->org->P01 = P01;
    A->sub[1]->org->P10 = CP;
    A->sub[1]->org->P11 = SP3;
    A->sub[1]->u1 = (A->u1 + A->u2) / (T)2.0;;
    A->sub[1]->u2 = A->u2;
    A->sub[1]->v1 = A->v1;
    A->sub[1]->v2 = (A->v1 + A->v2) / (T)2.0;    

    A->sub[2]->org->P00 = SP2;
    A->sub[2]->org->P01 = CP;
    A->sub[2]->org->P10 = P10;
    A->sub[2]->org->P11 = SP1;
    A->sub[2]->u1 = A->u1;
    A->sub[2]->u2 = (A->u1 + A->u2) / (T)2.0;
    A->sub[2]->v1 = (A->v1 + A->v2) / (T)2.0;
    A->sub[2]->v2 = A->v2;    

    A->sub[3]->org->P00 = CP;
    A->sub[3]->org->P01 = SP3;
    A->sub[3]->org->P10 = SP1;
    A->sub[3]->org->P11 = P11;
    A->sub[3]->u1 = (A->u1 + A->u2) / (T)2.0;;
    A->sub[3]->u2 = A->u2;
    A->sub[3]->v1 = (A->v1 + A->v2) / (T)2.0;
    A->sub[3]->v2 = A->v2;    
  }
  else
  {                    /****** Flaeche durch 2 Dreiecke approximierbar *****/
    A->linearized = 1;
    A->divided = 0;
  }
  return;
}

/*--------------------------------------------------------------------
  Check if a NURBS surface can be approximated with two triangles through
  its four corners. If not the surface is subdivided into four parts
----------------------------------------------------------------------*/
template
void LinearizeOrDivide( TessInfo<float,hpoint<float> > *A, const float tol, 
                        bool need_bbox );
template
void LinearizeOrDivide( TessInfo<float,point<float> > *A, const float tol, 
                        bool need_bbox );

template 
void LinearizeOrDivide( TessInfo<double,hpoint<double> > *A, const double tol, 
                        bool need_bbox );
template
void LinearizeOrDivide( TessInfo<double,point<double> > *A, const double tol, 
                        bool need_bbox );

/*-------------------------------------------------------------------
  Approximate NURBS-Surface with triangles
--------------------------------------------------------------------*/
template< class T, class HP >
void DoLinearizeSurface( 
                 int current_iter, int max_iter,
                 TessInfo<T,HP> *A,
                 const T tol,
                 void    outfunc( TessInfo<T,HP> *, void * ),                                   
                 void   *outfunc_data )
{
  int nA,i;
  TessInfo<T,HP> **pA;  

  /*
  cout << "iter: " << current_iter << "\n";
      
  cout << gul::dump_surf<T,HP>( 
              A->org->nu, A->org->pu, A->org->U,
              A->org->nv, A->org->pv, A->org->V,
              A->org->Pw ) << "\n";  
  */

/* --- A vierteln falls A nicht durch Dreiecke approximierbar: */    

  if( !A->linearized && !A->divided )
  {      
    if( current_iter < max_iter )
      gunu::LinearizeOrDivide( A, tol, false );
    else
      A->linearized = 1;
  }

  if( !A->linearized )
  {
    nA = 4;
    pA = A->sub;

    for( i = 0; i < nA; i++ )          /* Rekursion */
    {         
      DoLinearizeSurface( current_iter+1, max_iter, pA[i], tol,
                          outfunc, outfunc_data );
      delete pA[i];
      pA[i] = 0;      
    }
  }
  else       /* ------- Dreiecke ausgeben ------------------ */  
  {
    outfunc( A, outfunc_data ); 
  }
}
/*-------------------------------------------------------------------
  Approximate NURBS-Surface with triangles
--------------------------------------------------------------------*/
template void DoLinearizeSurface( 
                 int current_iter, int max_iter,
                 TessInfo<float,hpoint<float> > *A,
                 const float tol,
                 void    outfunc( TessInfo<float,hpoint<float> > *, void * ),                                   
                 void   *outfunc_data );
template void DoLinearizeSurface( 
                 int current_iter, int max_iter,
                 TessInfo<double,hpoint<double> > *A,
                 const double tol,
                 void    outfunc( TessInfo<double,hpoint<double> > *, void * ),                                   
                 void   *outfunc_data );

template void DoLinearizeSurface( 
                 int current_iter, int max_iter,
                 TessInfo<float,point<float> > *A,
                 const float tol,
                 void    outfunc( TessInfo<float,point<float> > *, void * ),                                   
                 void   *outfunc_data );
template void DoLinearizeSurface( 
                 int current_iter, int max_iter,
                 TessInfo<double,point<double> > *A,
                 const double tol,
                 void    outfunc( TessInfo<double,point<double> > *, void * ),                                   
                 void   *outfunc_data );

/*--------------------------------------------------------------------
  stores an intersection line segment
--------------------------------------------------------------------*/
template< class T >
inline void StoreIntersectionSegmentUV(
   point2<rational>                           *Suv,
   int                                        *Sflag,
   List<ListNode<IntersectionLineInfo<T> > >& I )
{
  ListNode<IntersectionLineInfo<T> > *in;
  T u,v;

  in = new ListNode<IntersectionLineInfo<T> >;
  I.Append(in);
  cnv2coord_rnd(Suv[0].x,&u);
  in->el.Suv[0].x = u - (T)1;
  cnv2coord_rnd(Suv[0].y,&v);
  in->el.Suv[0].y = v - (T)1;
  cnv2coord_rnd(Suv[1].x,&u);
  in->el.Suv[1].x = u - (T)1;
  cnv2coord_rnd(Suv[1].y,&v);
  in->el.Suv[1].y = v - (T)1;
  in->el.Sflag[0] = Sflag[0];
  in->el.Sflag[1] = Sflag[1];
}

// the 3D representation of the intersection segments is not really needed, 
// only for debugging
template< class T >
inline void StoreIntersectionSegment(
   point<rational>                            *S,
   IntersectInfo<T>                           *II )
{
  T scale,x,y,z;
  ListNode<line<T> > *in;

  in = new ListNode<line<T> >;
  II->IS.Append(in);

  scale = (T)1/II->scalei;

  cnv2coord_rnd(S[0].x,&x);
  in->el.P1.x = x * scale + II->minx;
  cnv2coord_rnd(S[0].y,&y);
  in->el.P1.y = y * scale + II->miny;
  cnv2coord_rnd(S[0].z,&z);
  in->el.P1.z = z * scale + II->minz;

  cnv2coord_rnd(S[1].x,&x);
  in->el.P2.x = x * scale + II->minx;
  cnv2coord_rnd(S[1].y,&y);
  in->el.P2.y = y * scale + II->miny;
  cnv2coord_rnd(S[1].z,&z);
  in->el.P2.z = z * scale + II->minz;
}

/*---------------------------------------------------------------------
  intersect the linear approximation of two nurbs surfaces
----------------------------------------------------------------------*/
template< class T, class HP >
void DoIntersectSurfaces( 
                 TessInfo<T,HP> *A, TessInfo<T,HP> *B,
                 T tol,
                 IntersectInfo<T> *II )
{
  int nA, nB, i, j;
  TessInfo<T,HP> **pA, **pB;  
  // IntersectionLineInfo<T> linfo;

  // do the bounding boxes overlap ?  
  if( (A->x1 > B->x2) || (A->x2 < B->x1) ||
      (A->y1 > B->y2) || (A->y2 < B->y1) ||
      (A->z1 > B->z2) || (A->z2 < B->z1) )
  {
    return;
  }   

  // split A into 4 parts, if A cannot be approximated by 2 triangles    
  if( !A->linearized && !A->divided )
  {      
    LinearizeOrDivide<T,HP>( A, tol, true );
  }

  // split B into 4 parts, if B cannot be approximated by 2 triangles    
  if( !B->linearized && !B->divided )
  {      
    LinearizeOrDivide<T,HP>( B, tol, true );
  }

  // recursion, if not both surfaces can be approximated by triangles
  if( !A->linearized || !B->linearized )
  {
    if( !A->linearized )
    {
      nA = 4;
      pA = A->sub;
    }
    else
    {
      nA = 1;
      pA = &A;        
    }

    if( !B->linearized )
    {
      nB = 4;
      pB = B->sub;
    }
    else
    {
      nB = 1;
      pB = &B;        
    }

    for( i = 0; i < nA; i++ )          // recursion
    {
      for( j = 0; j < nB; j++ )
      {
        DoIntersectSurfaces( pA[i], pB[j], tol, II );
      }      
    }
  }
  else  // intersect triangles pairwise
  {
    point<T>         PA[2][2],PB[2][2];
    point<rational>  PRA[2][2],PRB[2][2];
    point2<rational> PRA_UV[2][2],PRB_UV[2][2];
    T                x,y,z;
    unsigned long    *cbuf = 
      (unsigned long *)alloca(gul::rtr<T>::mantissa_length());
    triangle<rational> TA[2],TB[2];
    triangle2<rational> TAuv[2],TBuv[2];
    point<rational> S[2];
    point2<rational> S1uv[2], S2uv[2];
    int S1flag[2], S2flag[2];
    bool result;

    // convert A to rational representation
    PA[0][0] = A->org->P00;
    PA[0][1] = A->org->P01;
    PA[1][0] = A->org->P10;
    PA[1][1] = A->org->P11;

    for( i = 0; i < 2; i++ )
    {
      for( j = 0; j < 2; j++ )
      {
        x = (PA[i][j].x - II->minx)*II->scalei;
        y = (PA[i][j].y - II->miny)*II->scalei;
        z = (PA[i][j].z - II->minz)*II->scalei;

        PRA[i][j].x = rational(coord2int(x,cbuf),cbuf);
        PRA[i][j].y = rational(coord2int(y,cbuf),cbuf);
        PRA[i][j].z = rational(coord2int(z,cbuf),cbuf);
      }
    }
    // domain coordinates
    PRA_UV[0][0].x = rational(coord2int(A->u1+(T)1,cbuf),cbuf);
    PRA_UV[0][0].y = rational(coord2int(A->v1+(T)1,cbuf),cbuf);

    PRA_UV[0][1].x = rational(coord2int(A->u2+(T)1,cbuf),cbuf);
    PRA_UV[0][1].y = PRA_UV[0][0].y ;

    PRA_UV[1][0].x = PRA_UV[0][0].x;
    PRA_UV[1][0].y = rational(coord2int(A->v2+(T)1,cbuf),cbuf);

    PRA_UV[1][1].x = PRA_UV[0][1].x;
    PRA_UV[1][1].y = PRA_UV[1][0].y;

    // convert B to rational representation
    PB[0][0] = B->org->P00;
    PB[0][1] = B->org->P01;
    PB[1][0] = B->org->P10;
    PB[1][1] = B->org->P11;

    for( i = 0; i < 2; i++ )
    {
      for( j = 0; j < 2; j++ )
      {
        x = (PB[i][j].x - II->minx)*II->scalei;
        y = (PB[i][j].y - II->miny)*II->scalei;
        z = (PB[i][j].z - II->minz)*II->scalei;

        PRB[i][j].x = rational(coord2int(x,cbuf),cbuf);
        PRB[i][j].y = rational(coord2int(y,cbuf),cbuf);
        PRB[i][j].z = rational(coord2int(z,cbuf),cbuf);
      }
    }
   // domain coordinates
    PRB_UV[0][0].x = rational(coord2int(B->u1+(T)1,cbuf),cbuf);
    PRB_UV[0][0].y = rational(coord2int(B->v1+(T)1,cbuf),cbuf);

    PRB_UV[0][1].x = rational(coord2int(B->u2+(T)1,cbuf),cbuf);
    PRB_UV[0][1].y = PRB_UV[0][0].y;

    PRB_UV[1][0].x = PRB_UV[0][0].x;
    PRB_UV[1][0].y = rational(coord2int(B->v2+(T)1,cbuf),cbuf);

    PRB_UV[1][1].x = PRB_UV[0][1].x;
    PRB_UV[1][1].y = PRB_UV[1][0].y;

    // form triangles for A
    TA[0].P1 = PRA[0][0];
    TA[0].P2 = PRA[0][1];
    TA[0].P3 = PRA[1][0];
    TAuv[0].P1 = PRA_UV[0][0];
    TAuv[0].P2 = PRA_UV[0][1];
    TAuv[0].P3 = PRA_UV[1][0];

    TA[1].P1 = PRA[1][1];
    TA[1].P2 = PRA[1][0];
    TA[1].P3 = PRA[0][1];
    TAuv[1].P1 = PRA_UV[1][1];
    TAuv[1].P2 = PRA_UV[1][0];
    TAuv[1].P3 = PRA_UV[0][1];

    // form triangles for B
    TB[0].P1 = PRB[0][0];
    TB[0].P2 = PRB[0][1];
    TB[0].P3 = PRB[1][0];
    TBuv[0].P1 = PRB_UV[0][0];
    TBuv[0].P2 = PRB_UV[0][1];
    TBuv[0].P3 = PRB_UV[1][0];

    TB[1].P1 = PRB[1][1];
    TB[1].P2 = PRB[1][0];
    TB[1].P3 = PRB[0][1];
    TBuv[1].P1 = PRB_UV[1][1];
    TBuv[1].P2 = PRB_UV[1][0];
    TBuv[1].P3 = PRB_UV[0][1];

    result = IntersectTrianglesUV(TA[0],TAuv[0],TB[0],TBuv[0],
                                  S1uv,S1flag,S2uv,S2flag,S);
    if( result )
    {
      // cout << "Intersection line segment: ((" << S[0] << "), (" 
      //     << S[1] << "))\n"; 
      StoreIntersectionSegment(S,II);
      StoreIntersectionSegmentUV(S1uv,S1flag,II->IA);
      StoreIntersectionSegmentUV(S2uv,S2flag,II->IB);
    }

    result = IntersectTrianglesUV(TA[1],TAuv[1],TB[0],TBuv[0],
                                  S1uv,S1flag,S2uv,S2flag,S);
    if( result )
    {
      StoreIntersectionSegment(S,II);
      StoreIntersectionSegmentUV(S1uv,S1flag,II->IA);
      StoreIntersectionSegmentUV(S2uv,S2flag,II->IB);
    }

    result = IntersectTrianglesUV(TA[0],TAuv[0],TB[1],TBuv[1],
                                  S1uv,S1flag,S2uv,S2flag,S);
    if( result )
    {
      StoreIntersectionSegment(S,II);
      StoreIntersectionSegmentUV(S1uv,S1flag,II->IA);
      StoreIntersectionSegmentUV(S2uv,S2flag,II->IB);
    }

    result = IntersectTrianglesUV(TA[1],TAuv[1],TB[1],TBuv[1],
                                  S1uv,S1flag,S2uv,S2flag,S);
    if( result )
    {
      StoreIntersectionSegment(S,II);
      StoreIntersectionSegmentUV(S1uv,S1flag,II->IA);
      StoreIntersectionSegmentUV(S2uv,S2flag,II->IB);
    }
  }
}
// template instantiation
template void DoIntersectSurfaces( 
     TessInfo<float, hpoint<float> > *A, TessInfo<float,hpoint<float> > *B,
     float tol,
     IntersectInfo<float> *II );


/*---------------------------------------------------------------------
  this structure is used to find intersection segments which are
  connected at their endpoints with the help of a kd-tree
----------------------------------------------------------------------*/
template< class T >
// class is_uv_point : public gul::kdpoint<T,point2<T> >
class is_uv_point : public gul::point2<T>
{
public:
  IntersectionLineInfo<T> *ii;
  point2<T>               *ip;

  is_uv_point *next;
  is_uv_point *prev;

  is_uv_point( point2<T> *inP, IntersectionLineInfo<T> *inii )
    : point2<T>( *inP )
  {
    ip = inP;
    ii = inii;
  }
};

/*---------------------------------------------------------------------
  intersects the linear approximation of two nurbs surfaces
----------------------------------------------------------------------*/
template< class T, class HP >
GULAPI void IntersectSurfaces( 
            int nu1, int pu1, const Ptr< T >& U1,
            int nv1, int pv1, const Ptr< T >& V1,
            const Ptr< Ptr < HP > >& Pw1,
            int nu2, int pu2, const Ptr< T >& U2,
            int nv2, int pv2, const Ptr< T >& V2,
            const Ptr< Ptr < HP > >& Pw2,
            T tol,
            List<ListNode<IntersectionLineInfo<T> > >& S1,
            List<ListNode<IntersectionLineInfo<T> > >& S2,
            List<ListNode<line<T> > >& S )
{
  TessInfo<T,HP> A,B;
  point<T> maxP,minP;
  int j;
  IntersectInfo<T> II;
  T maxx,maxy,maxz,scale,d;
  rkdtree<point2<T>,T,random_des> KDT(2);
  is_uv_point<T> *Puv,*Peq,*Puv_next;
  List< is_uv_point<T> > LPuv,LPeq,LPtmp;
  ListNode<kdrec<point2<T>,T> > *nn;
  ListNode<IntersectionLineInfo<T> > *iin;
  Ptr<kdnnrec<point2<T>,T> > N; 
  List< gul::ListNode<kdrec<point2<T>,T> > > outrecs;
  point2<T> M;

  // initialize tesselation info structure for A
  A.org->nu = nu1;
  A.org->pu = pu1;  
  A.org->U =  U1;
  A.org->nv = nv1;
  A.org->pv = pv1;  
  A.org->V =  V1;
  A.org->Pw = Pw1;

  guge::CalcBoundingBoxVerts<HP,point<T> >( 
                  A.org->nu+1, A.org->Pw[0], minP, maxP );
  for( j = 1; j < A.org->nv+1; j++ )
    guge::UpdateBoundingBoxVerts<HP,point<T> >( 
                  A.org->nu+1, A.org->Pw[j], minP, maxP );
  A.x1 = minP.x; A.x2 = maxP.x; 
  A.y1 = minP.y; A.y2 = maxP.y;
  A.z1 = minP.z; A.z2 = maxP.z;

  A.org->P00 = euclid( A.org->Pw[0][0] );
  A.org->P01 = euclid( A.org->Pw[0][A.org->nu] );
  A.org->P10 = euclid( A.org->Pw[A.org->nv][0] );
  A.org->P11 = euclid( A.org->Pw[A.org->nv][A.org->nu] );

  A.u1 = 0.0;
  A.u2 = 1.0; 
  A.v1 = 0.0;
  A.v2 = 1.0; 

  // initialize tesselation info structure for B
  B.org->nu = nu2;
  B.org->pu = pu2;  
  B.org->U =  U2;
  B.org->nv = nv2;
  B.org->pv = pv2;  
  B.org->V =  V2;
  B.org->Pw = Pw2;

  guge::CalcBoundingBoxVerts<HP,point<T> >( 
                  B.org->nu+1, B.org->Pw[0], minP, maxP );
  for( j = 1; j < B.org->nv+1; j++ )
    guge::UpdateBoundingBoxVerts<HP,point<T> >( 
                  B.org->nu+1, B.org->Pw[j], minP, maxP );
  B.x1 = minP.x; B.x2 = maxP.x; 
  B.y1 = minP.y; B.y2 = maxP.y;
  B.z1 = minP.z; B.z2 = maxP.z;

  B.org->P00 = euclid( B.org->Pw[0][0] );
  B.org->P01 = euclid( B.org->Pw[0][B.org->nu] );
  B.org->P10 = euclid( B.org->Pw[B.org->nv][0] );
  B.org->P11 = euclid( B.org->Pw[B.org->nv][B.org->nu] );

  B.u1 = 0.0;
  B.u2 = 1.0; 
  B.v1 = 0.0;
  B.v2 = 1.0;   

  // initialize float to rational conversion parameters
  II.minx = Min(A.x1,B.x1);
  maxx = Max(A.x2,B.x2);
  II.miny = Min(A.y1,B.y1);
  maxy = Max(A.y2,B.y2);
  II.minz = Min(A.z1,B.z1);
  maxz = Max(A.z2,B.z2);
  
  scale = maxx - II.minx;
  d = maxy - II.miny;
  if( d > scale ) scale = d;
  d = maxz - II.minz;
  if( d > scale ) scale = d;
  if( !test(scale) ) return;
  II.minx -= scale;
  II.miny -= scale;
  II.minz -= scale;
  II.scalei = (T)1/scale;

  DoIntersectSurfaces(&A,&B,tol,&II);

  // cout << "number of intersection segments = " 
  //      << II.IA.nElems << "\n";

  // insert endpoints of intersection segments into kd-tree
  for( iin = II.IA.First(); iin; iin = iin->next )
  {
    Puv = new is_uv_point<T>(&iin->el.Suv[0],&iin->el);
    LPuv.Append(Puv);
    KDT.insert(Puv);

    Puv = new is_uv_point<T>(&iin->el.Suv[1],&iin->el);
    LPuv.Append(Puv);
    KDT.insert(Puv);
  }

  // compress points which are approximately equal into a single point
  for(;;)
  {
    Puv = LPuv.First();
    if( !Puv ) break;

    LPuv.Remove(Puv);
    KDT.remove(Puv,outrecs);
    // cout << "removed " << outrecs.nElems << " record(s)\n";
    LPtmp.Append(Puv);

    for(;;)
    {
      Puv = LPtmp.First();
      if( !Puv ) break;
      LPtmp.Remove(Puv);
      LPeq.Append(Puv);

      N.reserve_pool(1);
      if( KDT.nearest_neighbors(  *Puv, 1, N ) )
      {
        nn = N[0].m_recs.First();

        if( rel_equal(*Puv,*nn->el.m_base,rtr<T>::zero_tol()) )
        {
          for(  ; nn; nn = nn->next ) // all records have the same distance
          {
            // cast back to derived class
            Peq = (is_uv_point<T> *)nn->el.m_base;
            LPuv.Remove(Peq);
            KDT.remove(Peq,outrecs);
            // cout << "removed " << outrecs.nElems << " record(s)\n";
            LPtmp.Append(Peq);
          }
        }
      }
    }

    if( LPeq.nElems > 1 )
    {
      // calculate center of "equal" points
      Puv = LPeq.First();
      M = *Puv;
      j = 1;
      Puv = Puv->next;
      for( ; Puv; Puv = Puv->next )
      {
        M += *Puv;
        j++;
      }
      M = ((T)1/(T)j)*M;
 
      // change points in affected intersection segments 
      Puv = LPeq.First();
      LPeq.ReleaseElems();
      for( ; Puv; Puv = Puv_next )
      {
        *Puv->ip = M;
        Puv_next = Puv->next;
        delete Puv;
      }
    }
    else
      LPeq.DeleteElems();
  }

  S1 = II.IA;
  S2 = II.IB;
  S = II.IS;
}
// template instantiation
template GULAPI void IntersectSurfaces( 
            int nu1, int pu1, const Ptr< float >& U1,
            int nv1, int pv1, const Ptr< float >& V1,
            const Ptr< Ptr < hpoint<float> > >& Pw1,
            int nu2, int pu2, const Ptr< float >& U2,
            int nv2, int pv2, const Ptr< float >& V2,
            const Ptr< Ptr < hpoint<float> > >& Pw2,
            float tol,
            List<ListNode<IntersectionLineInfo<float> > >& S1,
            List<ListNode<IntersectionLineInfo<float> > >& S2,
            List<ListNode<line<float> > >& S );


#if 0
/*---------------------------------------------------------------------
  Intersect linear approximation of two NURBS-Surfaces
----------------------------------------------------------------------*/
/*
template< class T >
void IntersectSurfaces( 
                 IntersectInfo<T> *A, IntersectInfo<T> *B,
                 T tol, const rtl<T>& t,
                 void    outfunc( IntersectionLineInfo<T> *, void * ),                                   
                 void   *outfunc_data )
{
  int result, nA, nB, i, j;
  IntersectInfo<T> **pA, **pB;  
  IntersectionLineInfo<T> linfo;
  T u,v;
  
  // --- Ueberschneiden sich die Bounding Boxen ? -----------

  if( (A->x1 > B->x2) || (A->x2 < B->x1) ||
      (A->y1 > B->y2) || (A->y2 < B->y1) ||
      (A->z1 > B->z2) || (A->z2 < B->z1) )
  {
    return;
  }   
    
  // --- A vierteln falls A nicht durch Dreiecke approximierbar:    

  if( !A->f.linearized && !A->f.divided )
  {      
    LinearizeOrDivide<T>( A, tol, t );
  }
  if( !B->f.linearized && !B->f.divided )
  {      
    LinearizeOrDivide<T>( B, tol, t );
  }

  // --- falls noch nicht beide Flaechen durch Dreiecke darstellbar Rekursion:

  if( !A->f.linearized || !B->f.linearized )
  {
    if( !A->f.linearized )
    {
      nA = 4;
      pA = A->sub;
    }
    else
    {
      nA = 1;
      pA = &A;        
    }

    if( !B->f.linearized )
    {
      nB = 4;
      pB = B->sub;
    }
    else
    {
      nB = 1;
      pB = &B;        
    }

    for( i = 0; i < nA; i++ )          // Rekursion
    {
      for( j = 0; j < nB; j++ )
      {
        IntersectSurfaces( pA[i], pB[j], tol, t, outfunc, outfunc_data );
      }      
    }
  }
  else       // ------- Dreiecke schneiden ------------------
  {
    point<T> P00,P01,P10,P11;
    triangle<T> Atri[2], Btri[2];
    
    P00 = A->org.P00;
    P01 = A->org.P01;
    P10 = A->org.P10;
    P11 = A->org.P11;

    Atri[0].P1 = P00;
    Atri[0].P2 = P10;
    Atri[0].P3 = P01;

    Atri[1].P1 = P10;
    Atri[1].P2 = P11;
    Atri[1].P3 = P01;

    P00 = B->org.P00;
    P01 = B->org.P01;
    P10 = B->org.P10;
    P11 = B->org.P11;

    Btri[0].P1 = P00;
    Btri[0].P2 = P10;
    Btri[0].P3 = P01;

    Btri[1].P1 = P10;
    Btri[1].P2 = P11;
    Btri[1].P3 = P01;


    result = guge::IntersectTriangles<T>( Atri[0], Btri[0],
                                    &linfo.P1, &linfo.P2 );
    if( result != 0 )
    {
      // --- Parameterraumdarstellung bzgl. beider Flaechen bestimmen: ---

      BarycentricUV<T>( Atri[0].P1, Atri[0].P3, Atri[0].P2,
                        linfo.P1, &u, &v ); 
      linfo.F1.P1.x = A->u1 + u * (A->u2 - A->u1);
      linfo.F1.P1.y = A->v1 + v * (A->v2 - A->v1);

      BarycentricUV( Atri[0].P1, Atri[0].P3, Atri[0].P2,
                     linfo.P2, &u, &v ); 
      linfo.F1.P2.x = A->u1 + u * (A->u2 - A->u1);
      linfo.F1.P2.y = A->v1 + v * (A->v2 - A->v1);

      // --- 2. Flaeche: ---

      BarycentricUV<T>( Btri[0].P1, Btri[0].P3, Btri[0].P2,
                        linfo.P1, &u, &v ); 
      linfo.F2.P1.x = B->u1 + u * (B->u2 - B->u1);
      linfo.F2.P1.y = B->v1 + v * (B->v2 - B->v1);

      BarycentricUV<T>( Btri[0].P1, Btri[0].P3, Btri[0].P2,
                        linfo.P2, &u, &v ); 
      linfo.F2.P2.x = B->u1 + u * (B->u2 - B->u1);
      linfo.F2.P2.y = B->v1 + v * (B->v2 - B->v1);

      // --- abspeichern: ---

      outfunc( &linfo, outfunc_data );
    }
    
// ------------------------------------------------------------------------

    result = guge::IntersectTriangles<T>( Atri[0], Btri[1],
                                    &linfo.P1, &linfo.P2 );
    if( result != 0 )
    {
      // --- Parameterraumdarstellung bzgl. beider Flaechen bestimmen: ---

      BarycentricUV<T>( Atri[0].P1, Atri[0].P3, Atri[0].P2,
                        linfo.P1, &u, &v ); 
      linfo.F1.P1.x = A->u1 + u * (A->u2 - A->u1);
      linfo.F1.P1.y = A->v1 + v * (A->v2 - A->v1);

      BarycentricUV<T>( Atri[0].P1, Atri[0].P3, Atri[0].P2,
                        linfo.P2, &u, &v ); 
      linfo.F1.P2.x = A->u1 + u * (A->u2 - A->u1);
      linfo.F1.P2.y = A->v1 + v * (A->v2 - A->v1);

      // --- 2. Flaeche: ---

      BarycentricUV<T>( Btri[1].P2, Btri[1].P1, Btri[1].P3,
                        linfo.P1, &u, &v ); 
      linfo.F2.P1.x = B->u2 - u * (B->u2 - B->u1);
      linfo.F2.P1.y = B->v2 - v * (B->v2 - B->v1);

      BarycentricUV<T>( Btri[1].P2, Btri[1].P1, Btri[1].P3,
                        linfo.P2, &u, &v ); 
      linfo.F2.P2.x = B->u2 - u * (B->u2 - B->u1);
      linfo.F2.P2.y = B->v2 - v * (B->v2 - B->v1);

      // --- abspeichern: ---

      outfunc( &linfo, outfunc_data ); 
    }

/* ------------------------------------------------------------------------

    result = guge::IntersectTriangles<T>( Atri[1], Btri[0],
                                    &linfo.P1, &linfo.P2 );
    if( result != 0 )
    {
      // --- Parameterraumdarstellung bzgl. beider Flaechen bestimmen: ---

      BarycentricUV<T>( Atri[1].P2, Atri[1].P1, Atri[1].P3,
                        linfo.P1, &u, &v ); 
      linfo.F1.P1.x = A->u2 - u * (A->u2 - A->u1);
      linfo.F1.P1.y = A->v2 - v * (A->v2 - A->v1);

      BarycentricUV<T>( Atri[1].P2, Atri[1].P1, Atri[1].P3,
                        linfo.P2, &u, &v ); 
      linfo.F1.P2.x = A->u2 - u * (A->u2 - A->u1);
      linfo.F1.P2.y = A->v2 - v * (A->v2 - A->v1);

      // --- 2. Flaeche: ---

      BarycentricUV<T>( Btri[0].P1, Btri[0].P3, Btri[0].P2,
                        linfo.P1, &u, &v ); 
      linfo.F2.P1.x = B->u1 + u * (B->u2 - B->u1);
      linfo.F2.P1.y = B->v1 + v * (B->v2 - B->v1);

      BarycentricUV<T>( Btri[0].P1, Btri[0].P3, Btri[0].P2,
                        linfo.P2, &u, &v ); 
      linfo.F2.P2.x = B->u1 + u * (B->u2 - B->u1);
      linfo.F2.P2.y = B->v1 + v * (B->v2 - B->v1);

      // --- abspeichern: ---

      outfunc( &linfo, outfunc_data ); 
    }

// ----------------------------------------------------------------------

    result = guge::IntersectTriangles<T>( Atri[1], Btri[1],
                                    &linfo.P1, &linfo.P2 );
    if( result != 0 )
    {
      // --- Parameterraumdarstellung bzgl. beider Flaechen bestimmen: ---

      BarycentricUV<T>( Atri[1].P2, Atri[1].P1, Atri[1].P3,
                        linfo.P1, &u, &v ); 
      linfo.F1.P1.x = A->u2 - u * (A->u2 - A->u1);
      linfo.F1.P1.y = A->v2 - v * (A->v2 - A->v1);

      BarycentricUV<T>( Atri[1].P2, Atri[1].P1, Atri[1].P3,
                        linfo.P2, &u, &v ); 
      linfo.F1.P2.x = A->u2 - u * (A->u2 - A->u1);
      linfo.F1.P2.y = A->v2 - v * (A->v2 - A->v1);

      // --- 2. Flaeche: ---

      BarycentricUV<T>( Btri[1].P2, Btri[1].P1, Btri[1].P3,
                        linfo.P1, &u, &v ); 
      linfo.F2.P1.x = B->u2 - u * (B->u2 - B->u1);
      linfo.F2.P1.y = B->v2 - v * (B->v2 - B->v1);

      BarycentricUV<T>( Btri[1].P2, Btri[1].P1, Btri[1].P3,
                        linfo.P2, &u, &v ); 
      linfo.F2.P2.x = B->u2 - u * (B->u2 - B->u1);
      linfo.F2.P2.y = B->v2 - v * (B->v2 - B->v1);

      // --- abspeichern: ---

      outfunc( &linfo, outfunc_data ); 
    }
  }
}
*/

/*---------------------------------------------------------------------
  Intersect linear approximation of two NURBS-Surfaces
----------------------------------------------------------------------*/
/*
template
void IntersectSurfaces( 
                 IntersectInfo<float> *A, IntersectInfo<float> *B,
                 float tol, const rtl<float>& t,
                 void    outfunc( IntersectionLineInfo<float> *, void * ),                                   
                 void   *outfunc_data );
template
void IntersectSurfaces( 
                 IntersectInfo<double> *A, IntersectInfo<double> *B,
                 double tol, const rtl<double>& t,
                 void    outfunc( IntersectionLineInfo<double> *, void * ),                                   
                 void   *outfunc_data );
*/
#endif

}

---