include/shape/ssi/ssi_dsearch.hpp

Go to the documentation of this file.

#ifndef shape_ssi_dsearch_hpp
#define shape_ssi_dsearch_hpp
#include <shape/ssi/ssi_def.hpp>
#include <shape/ssi/ssi_qsegment.hpp>
#include <list>
#include <queue>
#include <shape/ssi/ssi_ordered_pair.hpp>
#include <shape/ssi/ssi_tri_tri.hpp>
#include <shape/ssi/ssi_math.hpp>

#define TMPL template<class C,class V>

#define SHIFT 3
#define __up__   0
#define __down__ 1
/* return the up triangle of quad q */
#define __up__triangle__(q)         q[0],q[3],q[2]
/* return the down triangle of quad q */
#define __down__triangle__(q)       q[0],q[1],q[2]
/* return the parametric plane of up triangle of q */
#define __up__param_triangle__(q)   q[3],q[2],q[0]
/* return the parametric plane of down triangle of q */
#define __down__param_triangle__(q) q[1],q[0],q[2]
#define __convert_order__(i,point) { point[(i+1)%2] = 1.0-point[(i+1)%2]; }
/* check for the intersection beetween {UP,DOWN}/{UP,DOWN} triangles */
#define __triangle_triangle_case__(trig0,trig1)                           \
{                                                                     \
  if ( geom::intersectp_triangles3_isegment                           \
       ( coplanar, seg[0], seg[1],                                    \
         trig0##triangle__(a),                                        \
         trig1##triangle__(b), point3::value_type(1e-12) ))           \
    {                                                                 \
      if ( !coplanar )                                                \
        {                                                             \
          space2prm(/* points in parameter space */                    \
                 seg0[0], seg0[1],                                    \
                 /* corresponding points in R3 */                     \
                  seg[0],  seg[1],                                    \
                 /* parametric plane = base + 2 vectors  */           \
                 trig0##param_triangle__(a)                           \
                 );                                                   \
          space2prm( seg1[0], seg1[1],                                 \
                    seg[0], seg[1],                                   \
                    trig1##param_triangle__(b) );                     \
          /* retrieve the correct offset coordinates in quad */       \
        __convert_order__(trig0,seg0[0]);                             \
        __convert_order__(trig0,seg0[1]);                             \
        __convert_order__(trig1,seg1[0]);                             \
        __convert_order__(trig1,seg1[1]);                             \
        }                                                             \
    else                                                              \
      {\
        seg0[0][0] = 1.0/3.0; seg0[0][1] = 1.0/3.0;                   \
        seg0[1][0] = 1.0/3.0; seg0[1][0] = 1.0/3.0;                   \
        seg1[0][0] = 1.0/3.0; seg1[0][1] = 1.0/3.0;                   \
        seg1[1][0] = 1.0/3.0; seg1[1][1] = 1.0/3.0;                   \
        __convert_order__(trig0,seg0[0]);                             \
        __convert_order__(trig0,seg0[1]);                             \
        __convert_order__(trig1,seg1[0]);                             \
        __convert_order__(trig1,seg1[1]);                             \
      };                                                              \
        {                                                             \
          f->convert(seg0,this,2);                                    \
          s->convert(seg1,this,2);                                    \
          push( seg0[0], seg0[1], seg1[0], seg1[1], seg[0],seg[1] );  \
        }                                                             \
    };                                                                \
}                                                                     \




namespace mmx {
namespace ssi {

void space2prm( vector2 & pa,
                vector2 & pb,
                const vector3& sa,
                const vector3& sb,
                const vector3& base,
                const vector3& pu,
                const vector3& pv );

typedef std::list< point2 >  list2_t;
typedef std::list< point3 >  list3_t;
typedef list2_t::iterator    p2ref_t;

struct dcurve
{
    unsigned size;
    dcurve  *   owner;
    list2_t     left,right;
    list3_t     inter;
    dcurve( const point2& l, const point2& r, const point3& i ) ;
    dcurve( const point2& l0, const point2& l1, const point2& r0, const point2& r1, const point3& i0, const point3& i1 );
};

template<class DPoint>
struct dim_cmp
{
    unsigned dim;
    dim_cmp(unsigned d = 0): dim(d) {};
    dim_cmp next() const { return dim_cmp((dim+1)%4); };
    bool operator()(  DPoint*const& a, DPoint *const& b ) const { return (*a)[dim] < (*b)[dim]; };
    bool operator()( const DPoint& a, const DPoint&  b ) const { return    a[dim] < b[dim];    };
};

template<class DP>
inline double distance( DP * pa, DP * pb )
{
    double acc = 0.0;
    for ( unsigned i = 0; i < 4; i++ )
        acc += ((*pa)[i]-(*pb)[i])*((*pa)[i]-(*pb)[i]);
    return acc;
};

struct sdpoint_t
{
    typedef double value_type;
    dcurve  *  curve;
    p2ref_t idl;
    p2ref_t idr;

    inline const double& operator[]( unsigned i ) const
    {
        if ( i < 2 ) return (*idl)[i];
        return (*idr)[i-2];
    };

    inline double& operator[]( unsigned i )
    {
        if ( i < 2 ) return (*idl)[i];
        return (*idr)[i-2];
    };
};

struct pdpoint_t
{
    typedef double value_type;
    dcurve * curve;
    p2ref_t ref;
    pdpoint_t * a;

    inline value_type& operator[]( unsigned i )
    {
        if ( i < 2 ) return (*ref)[i];
        return (*(a->ref))[i-2];
    };
    inline const value_type& operator[]( unsigned i ) const
    {
        if ( i < 2 ) return (*ref)[i];
        return (*(a->ref))[i-2];
    };
};

void append( dcurve *  a, dcurve   * b ) ;
void prepend( dcurve * a, dcurve * b ) ;
void reverse( dcurve * c ) ;


struct sdknode
{
    sdpoint_t   * data;
    sdknode * l, * r;
    sdknode( sdpoint_t * d, sdknode * _l, sdknode * _r ) : data(d), l(_l), r(_r){};
    ~sdknode() { if ( l ) delete l; if ( r ) delete r; };
};

struct pdknode
{
    pdpoint_t   * data;
    pdknode * l, * r;
    pdknode( pdpoint_t * d, pdknode * _l, pdknode * _r ) : data(d), l(_l), r(_r){};
    ~pdknode() { if ( l ) delete l; if ( r ) delete r; };
};



inline  bool   empty( dcurve * c ) { return c->owner != c; };
inline  pdpoint_t * sibble( pdpoint_t * p )  { return p->a; };

void extremums( sdpoint_t* dst, dcurve* src );
void extremums( pdpoint_t * dst, dcurve * c );
dcurve * curve( pdpoint_t * p );
dcurve *  curve( sdpoint_t * p );

inline  bool isfirst ( sdpoint_t* p ) { return curve(p)->left.begin() == p->idl;     };
inline  bool islast  ( sdpoint_t* p ) { return (--(curve(p)->left.end())) == p->idl; };

inline bool isfirst ( pdpoint_t* p ) { return curve(p)->left.begin() == p->ref || curve(p)->right.begin() == p->ref; };
inline bool islast ( pdpoint_t* p ) { return (--(curve(p)->left.end())) == p->ref ||(--(curve(p)->right.end())) == p->ref;};

inline bool isextrem( pdpoint_t* p  ) { return  isfirst(p) || islast(p); };
inline bool isextrem( sdpoint_t*  p ) { return  isfirst(p) || islast(p); };

void link( sdpoint_t* a, sdpoint_t* b);
void link( pdpoint_t* a, pdpoint_t* b);

inline void swap( dcurve * c ) { c->left.swap( c->right ); };

inline bool isright( pdpoint_t* p )
{ return p->ref == curve(p)->right.begin() || p->ref == --(curve(p)->right.end());};


template<class DPoint>
inline void _link(  DPoint* a, DPoint* b )
{
    if ( islast(a)   )
    {
        if( islast(b) )
            reverse(curve(b));
        // last a + first b
        append(curve(a),curve(b));
    }
    else
    {
        if ( isfirst(b) )
            reverse(curve(b));
        // last b + first a
        prepend(curve(a),curve(b));
    };
};


sdknode * make( sdpoint_t * first, sdpoint_t * last, const dim_cmp<sdpoint_t>& f ) ;
pdknode * make( pdpoint_t ** first, pdpoint_t ** last, const dim_cmp<pdpoint_t>& f ) ;



template<class DPoint>
struct assoc_t
{
    std::pair<DPoint*,DPoint*>         pp;
    double                              d;
    // comparing association on pointers
    struct pair_cmp
    { bool operator()( const assoc_t& a1, const assoc_t& a2 ) const { return a1.pp < a2.pp;  }; };
    // comparing association on distances
    struct dist_cmp
    { bool operator()( const assoc_t& a1, const assoc_t& a2 ) const { return a1.d > a2.d;    }; };
    assoc_t(){};
    assoc_t( DPoint * a, DPoint * b, double _d ):
        pp(a,b), d(_d){};
};

typedef std::priority_queue< assoc_t<pdpoint_t>,
std::vector< assoc_t<pdpoint_t> >,
assoc_t< pdpoint_t>::dist_cmp > pheap_t;
typedef std::priority_queue< assoc_t<sdpoint_t>,
std::vector< assoc_t<sdpoint_t> >,
assoc_t< sdpoint_t>::dist_cmp > sheap_t;

typedef std::set< assoc_t< pdpoint_t > ,  assoc_t<pdpoint_t>::pair_cmp > pset_t;

typedef assoc_t<sdpoint_t> sdassc_t;
typedef assoc_t<pdpoint_t> pdassc_t;

void solve_sheap( sheap_t& h );
void solve_pheap( pheap_t& h );

typedef std::pair<pdpoint_t*,pdpoint_t*> ppair_t;
typedef std::list< ppair_t > links_t;
typedef std::map< dcurve*, links_t > curves_links_t;
typedef ordered_pair<dcurve*> cpair_t;

template <class A, class B>
inline void container_add( A& h, const B& a ) { h.push(a); };
inline void container_add( pset_t& s, const pdassc_t& a ) { s.insert(a); };
inline cpair_t curves( const pdassc_t& ass ) { return cpair_t(curve(ass.pp.first),curve(ass.pp.second)); };
inline pdassc_t sibbles( const pdassc_t& ass )
{
    pdpoint_t * a = sibble(ass.pp.first);
    pdpoint_t * b = sibble(ass.pp.second);
    if ( a < b )
        return pdassc_t(a,b,0);
    else
        return pdassc_t(b,a,0);
};

void satisfy_links(  curves_links_t::iterator it );
void solve_pset( curves_links_t& links, pset_t& s );
std::list<dcurve*> * reduce2( std::vector< dcurve * > * conflicts );
typedef  pset_t::iterator pset_iterator;

void build_sheap( sheap_t& h, sdpoint_t * first, sdpoint_t * last, double prec );

TMPL
struct dsearch : qsegment<C,V>
{
    typedef typename qsegment<C,V>::ParametricSurface ParametricSurface;
    typedef qnode<C,V>* qtree;
    bool m_cnf;

    //------------------------------------------------------------------------------
    std::vector< dcurve * >  * curr_config;           // represent the current region/region conflict resultings segments
    std::list  < dcurve * >  ** _branches;            // represent the reduced set of curves for each region/region conflict

    void cnfpush( qtree f, qtree s );
    void push( qtree f, qtree s );
    inline  std::list  < dcurve * >  *& branches( unsigned i, unsigned j )
    { return _branches[i*this->regions.size()+j]; };
    // add a pair of segment to the segment soup
    void push (  const point2& l0,  const point2& l1,
                 const point2& r0,  const point2& r1,
                 const point3& i0,  const point3& i1 )
    {
        curr_config->push_back( new dcurve( l0,l1,r0, r1, i0, i1 ) );
    };

    // use of the qsegment_t bounding box hierarchy to retrieve unit quad that intersects
#define push_s push
#define push_f push
    // search specialisation s is known to be a leaf
    void search_f(  qtree f, qtree const s );
    // same thing here f is a leaf
    void search_s(   qtree const f, qtree s );
    // search generic case: both f and s are root node
    void search(qtree f, qtree s );
    
    // launch the search for region i vs region j
    void search( rid_t i, rid_t j );

    // try to reduce the number of double curves
    std::list<dcurve*> * reduce( std::vector< dcurve * > * conflicts );
    inline void find_branches( rid_t id0, rid_t id1 )
    {
        branches( id0, id1 ) = reduce( curr_config );
        delete curr_config;
    };

    // stats structure
    struct
    {
        double search_time;
        double recons_time;
    } stats;

    std::list< dcurve * >  *result;

    point2 ** lefts;
    point2 ** rights;

    unsigned * sizes;
    unsigned nbcurves;
    unsigned nbpoints;
    unsigned conflict_count;

    void GnuPlotOutput( char * name );

    std::vector< point3 > cnfdata3;
    std::vector< point2 > cnfdata2;


    double _st;
    dsearch( const ParametricSurface * s , unsigned w, unsigned h, bool cnf = false  );
    ~dsearch();
};

TMPL
dsearch<C,V>::dsearch( const ParametricSurface * s , unsigned w, unsigned h, bool cnf  ) : qsegment<C,V>(s, w, h )
{
    m_cnf = cnf;
    //    _st = time();
    unsigned i,j;
    conflict_count = 0;
    nbpoints = 0;
    _branches = new std::list< dcurve * > *[this->regions.size()*this->regions.size()];
    memset(_branches,0,sizeof( std::list< dcurve * > * )*
           this->regions.size()*this->regions.size() );

    for ( i = 0; i < this->regions.size() ; i++ )
        for ( j = i+1; j < this->regions.size() ; j++   )
        {
#define VOISINS_AUSSI
#ifndef VOISINS_AUSSI
            if ( !this->grp->neighbors(i,j) )
#endif
            {
                //cout << i << " " << j << endl;
                search(i,j);
                find_branches(i,j);
            };
        };

    //      cout << "search: ok\n";
    std::vector< dcurve * > * tmp = new std::vector<dcurve*>();

    for ( i = 0; i < this->regions.size(); i++ )
        for ( j = i+1; j < this->regions.size(); j++ )
        {
            if ( branches(i,j) )
            {
                for (  std::list<dcurve*>::iterator it = branches(i,j)->begin(); it != branches(i,j)->end(); it ++ )
                {
                    tmp->push_back( *it );
                };
                delete branches(i,j);
            }
        };

    result = reduce2( tmp );
    delete tmp;
    nbcurves = 0;
    lefts = rights = 0;
    sizes = 0;
    if ( result )
    {
        nbpoints = 0;
        std::list<dcurve *>::iterator it;
        nbcurves = result->size();

        sizes  = new unsigned[nbcurves];

        for ( i = 0, it = result->begin(); it != result->end(); it++ , i++ )
        {
            sizes[i] = (*it)->left.size();
            nbpoints+= sizes[i];
        };

        lefts  = new point2*[nbcurves];
        rights = new point2*[nbcurves];
        unsigned c_n = 0;
        unsigned k;
        for (  it = result->begin(); it != result->end(); it ++, c_n++  )
        {
            lefts[c_n]  = new point2[sizes[c_n]];
            rights[c_n] = new point2[sizes[c_n]];
            list2_t::iterator itcr, itcl;
            for ( k = 0, itcl = (*it)->left.begin(), itcr = (*it)->right.begin();
                  itcl != (*it)->left.end();
                  itcl++ , itcr++, k++ )
            {

                lefts[c_n][k] = *itcl;
                rights[c_n][k] = *itcr;
            };
        };

        for (  std::list<dcurve*>::iterator it = result->begin(); it!= result->end(); it++ )
        {
            delete *it;
        };
        delete result;
    };
    delete[] _branches;

    //std::cout << "d = " << (time()-_st) << std::endl;
    //    std::cout << nbcurves << std::endl;

    //    GnuPlotOutput("dsearch.gmv");
};

TMPL
dsearch<C,V>::~dsearch()
{
    if ( sizes ) delete[] sizes;
    if ( lefts )
    {
        for ( unsigned i = 0; i < nbcurves; i++ )
            delete[] lefts[i];
        delete[] lefts;
    };
    if ( rights )
    {
        for ( unsigned i = 0; i < nbcurves; i++ )
            delete[] rights[i];
        delete[] rights;
    };

};

TMPL
void dsearch<C,V>::GnuPlotOutput( char * name )
{
    /*
    gmv::stream out(name);

    std::vector< vector3 > v;
    for ( unsigned i = 0; i < nbcurves; i++ )
      {
        for ( unsigned k = 0; k < sizes[i]; k++ )
          {
            v.push_back(vector3());
            v.back()[0] = lefts[i][k][0];
            v.back()[1] = lefts[i][k][1];
            v.back()[2] = 0;
            v.push_back(vector3());
            v.back()[0] = rights[i][k][0];
            v.back()[1] = rights[i][k][1];
            v.back()[2] = 0;
          };
      };
    out.points(&(v[0][0]),v.size());
    */
}

TMPL
void dsearch<C,V>::cnfpush( qtree f, qtree s )
{
    point3 a[4],b[4]; //      point3 a[4]; point3 b[4];
    f->fill(a,this);
    s->fill(b,this);
    cnfdata3.push_back( ((double)(1.0/4.0))*(a[0]+a[1]+a[2]+a[3]) );
};

// push an quad / quad conflict
// intersection is computed if any and added to the current conflict set of segments
TMPL
void dsearch<C,V>::push( qtree f, qtree s )
{
    if ( std::abs(f->umin-s->umin)+std::abs(f->vmin-s->vmin) < SHIFT ) return;
    conflict_count ++;
    point3 a[4],b[4];
    qsegment<C,V>::scale_conflict(a,b,f,s);
    {
        point3  seg[2];
        point2 seg0[2], seg1[2];
        bool   coplanar;
        __triangle_triangle_case__(__up__,__up__);
        __triangle_triangle_case__(__up__,__down__);
        __triangle_triangle_case__(__down__,__down__);
        __triangle_triangle_case__(__down__,__up__);
    }
};
// search specialisation s is known to be a leaf
TMPL
void dsearch<C,V>::search_f(  qtree f, qtree const s )
{

    if (intersectp(f->box,s->box))
    {
        f->split(this);
        if ( leaf(f) ){ push_f(f,s); return; };
        search_f(f->l,s);
        search_f(f->r,s);
    };
};

// same thing here f is a leaf
TMPL
void dsearch<C,V>::search_s(   qtree const f, qtree s )
{

    if ( intersectp( f->box, s->box ) )
    {
        s->split(this);
        if ( leaf(s) ) { push_s(f,s); return; };
        search_s(f,s->l);
        search_s(f,s->r);
    };
};

// search generic case: both f and s are root node
TMPL
void dsearch<C,V>::search(  qtree f, qtree s )
{


    // if there bounding boxes intersect
    if ( intersectp(f->box,s->box) )
    {
        // compute if needed the next level in the hierarchy
        f->split(this);
        s->split(this);
        // if f and s are root node
        if ( !(f->leaf()) && !(s->leaf()) )
        {
            // check all the possible intersection between  sons
            search(f->l,s->l), search(f->l,s->r);
            search(f->r,s->l), search(f->r,s->r);
            return;
        };

        // s is leaf (not f) switch to specialized case
        if ( !leaf(f) )
        {
            search_f(f->l,s);
            search_f(f->r,s);
            return;
        };
        // same thing with f
        if ( !leaf(s) )
        {
            search_s(f,s->l);
            search_s(f,s->r);
            return;
        };
        // two unit quad had been reached push the conflict
        // the function push( qtree, qtree ) compute the intersection on the fly if any
        push(f,s);
    };
};

// launch the search for region i vs region j
TMPL
void dsearch<C,V>::search( rid_t i, rid_t j )
{
    curr_config = new std::vector< dcurve* >();
    search( this->trees[i], this->trees[j] );
};


TMPL
std::list<dcurve*> * dsearch<C,V>::reduce( std::vector< dcurve * > * conflicts )
{
    unsigned i;
    if ( conflicts->size() == 0 ) return 0;

    unsigned   endsize  = 2*conflicts->size();
    sdpoint_t * ends = new sdpoint_t[endsize];

    for ( i = 0; i < endsize; i+= 2 )
        extremums( ends + i, (*conflicts)[i/2] );

    sheap_t h;
    build_sheap(h,ends, ends + endsize, 0.04 );
    solve_sheap(h);
    delete[] ends;
    std::list< dcurve* > * result = new std::list< dcurve * >();
    for ( i = 0; i < conflicts->size() ; i++ )
    {
        if ( empty((*conflicts)[i]) ) delete (*conflicts)[i];
        else result->push_front( (*conflicts)[i] );
    };
    return result;
};


} // namespace shape_ssi

} //namespace mmx
# undef ParametricSurface
# undef TMPL
# endif

---