synaps/topology/TopSbdCurve3d.h

00001 /***************************************************************************
00002  *   Copyright (C) 2005 by Chen Liang                                      *
00003  *   cliang@cs.hku.hk                                                      *
00004  *                                                                         *
00005  *   This program is free software; you can redistribute it and/or modify  *
00006  *   it under the terms of the GNU General Public License as published by  *
00007  *   the Free Software Foundation; either version 2 of the License, or     *
00008  *   (at your option) any later version.                                   *
00009  *                                                                         *
00010  *   This program is distributed in the hope that it will be useful,       *
00011  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
00012  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
00013  *   GNU General Public License for more details.                          *
00014  *                                                                         *
00015  *   You should have received a copy of the GNU General Public License     *
00016  *   along with this program; if not, write to the                         *
00017  *   Free Software Foundation, Inc.,                                       *
00018  *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
00019  ***************************************************************************/
00020 
00021 #ifndef synaps_topology_TopSbdCurve3d_h
00022 #define synaps_topology_TopSbdCurve3d_h
00023 
00024 //--------------------------------------------------------------------
00025 #include <vector>
00026 #include <list>
00027 #include <math.h>
00028 #include <synaps/upol/bezier/BEZIER.m>
00029 #include <synaps/mpol/MPol.h> 
00030 #include <synaps/mpol/MPolDse.h> 
00031 #include <synaps/msolve/sbdslv.h> 
00032 #include <synaps/topology/Cell3dN.h>
00033 #include <synaps/topology/Octree.h>
00034 //--------------------------------------------------------------------
00035 __BEGIN_NAMESPACE_SYNAPS
00036 //--------------------------------------------------------------------
00037 template<class C, 
00038          class CELL= Cell3dN<5,C>, 
00039          class MSLV= SbdSolver<typename CELL::coeff_t, SBDSLV_RDL> >
00040 struct TopSbdCurve3d
00041 {
00042 
00043   //type definitions
00044   typedef TopSbdCurve3d<C, CELL,MSLV>         self_t;
00045 
00046   typedef typename CELL::coeff_t         coeff_t;
00047 
00048   typedef typename CELL::MPOL_T          MPOL_T;
00049   typedef MPolDse<C, mpoldse::bernstein<C> > mpolb_t;
00050 
00051   typedef topology::box3d<coeff_t>       box_t;
00052   typedef topology::point<C>             point_t;
00053   typedef topology::point_graph<C>       pointgraph_t;
00054 
00055   
00056   typedef CELL                           cell_t;
00057   typedef OCTREE::Node<cell_t*>           node_t;
00058   typedef OCTREE::Octree<cell_t*>         octree_t;
00059   
00060   typedef intersection3d<C>       intersection_t;
00061   typedef shared_object<intersection_t*>  intersection_ptr;
00062   //  typedef intersection3d<coeff_t>         intersection_t;
00063   typedef Seq<intersection3d<coeff_t>*>   intersectionlist_t;
00064   //  typedef shared_object<intersection_t*>  intersection_ptr;
00065   typedef Seq<intersection_ptr>           intersection_ptr_list;
00066   typedef topology::point_graph<coeff_t,  intersection_ptr_list> pointgraph_int_t;
00067 
00068 protected:
00069 // data members
00070 
00071     C           m_dEpsilon; // maximum tolerance
00072     int         m_nMinSbd;
00073     octree_t    m_cOcTree;  // the octree structure
00074 
00075     int         m_nCount; // traveral count
00076 
00077 public:
00078 // constructors/distructor
00079 
00080     TopSbdCurve3d():m_dEpsilon(1e-3), m_nMinSbd(0),m_nCount(0) {}
00081     TopSbdCurve3d(coeff_t e):m_dEpsilon(e), m_nMinSbd(0),m_nCount(0) {}
00082     TopSbdCurve3d(coeff_t e, int d):m_dEpsilon(e), m_nMinSbd(d),m_nCount(0) {}
00083 
00084     ~TopSbdCurve3d();
00085 
00086     //int regularity(cell_t & c);
00087 
00088     // compute intersection on the dedicated face for the given polynomial
00089     void     computeinter_face(intersectionlist_t& inters, cell_t& c, int nAxis, int nFace);
00090  
00091     void     computeintersections(intersectionlist_t& inters, cell_t& c);
00092 
00093     template<class POINTGRAPH>
00094     void run(POINTGRAPH& g, CELL* rc);
00095 
00096 protected:
00097 // helpers
00098 
00099     // this function subdivide a given node
00100     void subdivide(node_t* n, std::vector<node_t*>& aNodeList);
00101     // this function subdivide a cell
00102     void subdivide(CELL* cell, std::vector<CELL*>& subCells);
00103 
00104     // housekeeping: octree
00105     void deleteNode(node_t* n, std::vector<intersection_t*>& inters);
00106 };
00107 
00109 // template instantiation       //
00111 
00112 //inline template class TopSbdCurve3d< double, double, bezier::bzslv<double> >;
00113 
00115 // constructor/distructor       //
00117 template< class X, class CELL, class MSLV>
00118 TopSbdCurve3d<X, CELL, MSLV>::~TopSbdCurve3d()
00119 {
00120   // release the memory occupied by the octree
00121   m_nCount++;
00122   std::vector<intersection_t*> inters;
00123   deleteNode(m_cOcTree.GetRoot(), inters);
00124   for (int i = 0; i < inters.size(); i++)
00125     delete inters[i];
00126 }
00127 
00128 //--------------------------------------------------------------------
00130 /***************************************/
00131 /*               Z^                    */
00132 /*                |_____               */
00133 /*               /|    /|              */
00134 /*     N  ^     /_|_B_/ |              */
00135 /*        |    |  |___|_|___> Y        */
00136 /*        |    | /F   | /              */
00137 /*             |/_____|/               */
00138 /*            |/_                      */
00139 /*           X                         */
00140 /*                                     */
00141 /*    Fig.1  The coordinate system     */
00142 /***************************************/
00144 
00145 template< class C, class CELL, class MSLV>
00146 void TopSbdCurve3d<C, CELL, MSLV>::deleteNode(node_t* n, std::vector<intersection_t*>& inters)
00147 {
00148     if (n->GetF_NEchild())
00149         deleteNode(n->F_NEchild, inters);
00150     if (n->GetF_NWchild())
00151         deleteNode(n->F_NWchild, inters);
00152     if (n->GetF_SEchild())
00153         deleteNode(n->F_SEchild, inters);
00154     if (n->GetF_SWchild())
00155         deleteNode(n->F_SWchild, inters);
00156     if (n->GetB_NEchild())
00157         deleteNode(n->B_NEchild, inters);
00158     if (n->GetB_NWchild())
00159         deleteNode(n->B_NWchild, inters);
00160     if (n->GetB_SEchild())
00161         deleteNode(n->B_SEchild, inters);
00162     if (n->GetB_SWchild())
00163         deleteNode(n->B_SWchild, inters);
00164 
00165     for (int i = 0; i < n->GetCell()->m_aInterList.size(); i++)
00166     {
00167         if (n->GetCell()->m_aInterList[i]->nCount != m_nCount)
00168         {
00169             inters.push_back(n->GetCell()->m_aInterList[i]);
00170             n->GetCell()->m_aInterList[i]->nCount=m_nCount;
00171         }
00172     }
00173 
00174     delete n->GetCell();
00175     delete n;
00176 }
00177 //--------------------------------------------------------------------
00178 template<class C, class CELL,  class MSLV>
00179 void TopSbdCurve3d<C, CELL, MSLV>::subdivide(node_t* n, std::vector<node_t*>& aNodeList)
00180 {
00181   using namespace OCTREE;
00182   cell_t* c = n->GetCell();
00183   std::vector<cell_t*> aSubCells;
00184   
00185   // subdivide a cell into 8 subcells
00186   subdivide(c,aSubCells);
00187   
00188   // encapsulate the cells with the node of the octree
00189   aNodeList.push_back(new node_t(B_SW, n, aSubCells[0]));
00190   aNodeList.push_back(new node_t(F_SW, n, aSubCells[1]));
00191   aNodeList.push_back(new node_t(B_SE, n, aSubCells[2]));
00192   aNodeList.push_back(new node_t(F_SE, n, aSubCells[3]));
00193   aNodeList.push_back(new node_t(B_NW, n, aSubCells[4]));
00194   aNodeList.push_back(new node_t(F_NW, n, aSubCells[5]));
00195   aNodeList.push_back(new node_t(B_NE, n, aSubCells[6]));
00196   aNodeList.push_back(new node_t(F_NE, n, aSubCells[7]));
00197 }
00198 //--------------------------------------------------------------------
00199 
00201 /***************************************/
00202 /*          ^ Z                        */
00203 /*          |                          */
00204 /*          |                          */
00205 /*          |                          */
00206 /*          |4----6                    */
00207 /*         /||___/|_______> Y          */
00208 /*        5/-0--7 2                    */
00209 /*        /|/   |/                     */
00210 /*       / 1----3                      */
00211 /*     \/_                             */
00212 /*     X                               */
00213 /*                                     */
00214 /*    Fig.2  The layout of the         */
00215 /*        subdivided cells             */
00216 /***************************************/
00218 
00220 /***************************************/
00221 /*          ^ Z                        */
00222 /*          |   5     0                */
00223 /*          |  |     /                 */
00224 /*          |  v   |/_                 */
00225 /*          |------                    */
00226 /*         /||___/|_______> Y          */
00227 /*  2-->  -/ ---- | <--3               */
00228 /*        /|/ 1 |/                     */
00229 /*       / ------                      */
00230 /*     \/_    ^                        */
00231 /*     X      |                        */
00232 /*             4                       */
00233 /*    Fig.3  Faces index of a cell     */
00234 /***************************************/
00236 // this function spawns 8 evenly subdivided child cell of this call
00237 template < class C, class CELL, class MSLV >
00238 void TopSbdCurve3d<C, CELL, MSLV>::subdivide(CELL* cell,
00239                                             std::vector<CELL*>& subCells)
00240 {
00241   typedef typename CELL::intersection_t intersection_t;
00242   typedef typename CELL::intersectionlist_t intersectionlist_t;
00243   typedef typename CELL::MPOL_T MPOL_T;
00244     // define the subdivision points
00245 #ifdef SYNAPS_SBD3D_RNDFACTORWHENSUBDIVISION
00246     // a perturbation ranged from 0.005-0.1
00247     C t[] = {(C(rand())/C(RAND_MAX)-C(0.5))*C(0.095)+C(0.5),
00248              (C(rand())/C(RAND_MAX)-C(0.5))*C(0.095)+C(0.5),
00249              (C(rand())/C(RAND_MAX)-C(0.5))*C(0.095)+C(0.5)};
00250     t[0]+= (t[0]>0.5?0.005:-0.005);
00251     t[1]+= (t[1]>0.5?0.005:-0.005);
00252     t[2]+= (t[2]>0.5?0.005:-0.005);
00253     C mid[] = {cell->m_cBox.x0*(C(1)-t[0])+cell->m_cBox.x1*t[0], cell->m_cBox.y0*(C(1)-t[1])+cell->m_cBox.y1*t[1], cell->m_cBox.z0*(C(1)-t[2])+cell->m_cBox.z1*t[2]};
00254 #else
00255     C mid[] = {(cell->m_cBox.x0+cell->m_cBox.x1)/C(2), (cell->m_cBox.y0+cell->m_cBox.y1)/C(2), (cell->m_cBox.z0+cell->m_cBox.z1)/C(2)};
00256 #endif
00257     C box[] = {cell->m_cBox.x0, cell->m_cBox.x1, cell->m_cBox.y0, cell->m_cBox.y1, cell->m_cBox.z0, cell->m_cBox.z1};
00258 
00259     // create new cells
00260     subCells.push_back(new CELL(box[0], mid[0], box[2], mid[1], box[4], mid[2]));
00261     subCells.push_back(new CELL(mid[0], box[1], box[2], mid[1], box[4], mid[2]));
00262     subCells.push_back(new CELL(box[0], mid[0], mid[1], box[3], box[4], mid[2]));
00263     subCells.push_back(new CELL(mid[0], box[1], mid[1], box[3], box[4], mid[2]));
00264     subCells.push_back(new CELL(box[0], mid[0], box[2], mid[1], mid[2], box[5]));
00265     subCells.push_back(new CELL(mid[0], box[1], box[2], mid[1], mid[2], box[5]));
00266     subCells.push_back(new CELL(box[0], mid[0], mid[1], box[3], mid[2], box[5]));
00267     subCells.push_back(new CELL(mid[0], box[1], mid[1], box[3], mid[2], box[5]));
00268 
00269     // update the depth count for each subdivided cells
00270     for (int i = 0; i < 8; i++)
00271     {
00272         subCells[i]->m_nDepth = cell->m_nDepth+1;
00273     }
00274 
00275     // the subcells inherits the previous computed intersections on its parent cell
00276     for (int i = 0; i < cell->m_aInterList.size(); i++)
00277     {
00278         intersection_t* inter = cell->m_aInterList[i];
00279 
00280         // assign the intersections on the parent cell to subdivided cells, according to the plane
00281         // they are on and the planar coordinates
00282         int nidx = inter->nPlane; // testing the plane YZ (nidx=0) or XZ (nidx=1) or XZ (nidx=2)
00283         int *cidx=subcellidx+(nidx<<3); // the right set of subcell indices to be assigned each intersection
00284         int axisidx0 = plane_axisidx[nidx<<1], axisidx1 = plane_axisidx[(nidx<<1)+1]; // the first and the second axis indices, {1, 2} for YZ plane, {0, 2} for XZ plane and {0, 1} for XY plane
00285 
00286         assert(_is_verysmall(inter->pnt[nidx]-box[nidx<<1]) || _is_verysmall(box[(nidx<<1)+1]-inter->pnt[nidx]));
00287         assert(inter->pnt[axisidx0] >= box[axisidx0<<1] && inter->pnt[axisidx0] <= box[(axisidx0<<1)+1] &&
00288                inter->pnt[axisidx1] >= box[axisidx1<<1] && inter->pnt[axisidx1] <= box[(axisidx1<<1)+1]);
00289 
00290         if (inter->pnt[axisidx0] < mid[axisidx0]) // if smaller than the first mid-point
00291         {
00292             if (inter->pnt[axisidx1] < mid[axisidx1]) // if smaller than the second mid-point
00293             {
00294                 if (inter->pnt[nidx]-box[nidx<<1] < box[(nidx<<1)+1]-inter->pnt[nidx]) // if the plane corresponds to x (or y, or z) = a, where x (or y, or z) is in [a, b]
00295                     subCells[cidx[0]]->m_aInterList.push_back(inter);
00296                 else // if the plane corresponds to x (or y, or z) = b, where x (or y, or z) is in [a, b]
00297                     subCells[cidx[1]]->m_aInterList.push_back(inter);
00298             }
00299             else // if larger than the first mid-point
00300             {
00301                 if (inter->pnt[nidx]-box[nidx<<1] < box[(nidx<<1)+1]-inter->pnt[nidx]) // if the plane corresponds to x (or y, or z) = a, where x (or y, or z) is in [a, b]
00302                     subCells[cidx[2]]->m_aInterList.push_back(inter);
00303                 else
00304                     subCells[cidx[3]]->m_aInterList.push_back(inter);
00305             }
00306         }
00307         else // if lager than the first mid-point
00308         {
00309             if (inter->pnt[axisidx1] < mid[axisidx1]) // if smaller than the second  mid-point
00310             {
00311                 if (inter->pnt[nidx]-box[nidx<<1] < box[(nidx<<1)+1]-inter->pnt[nidx]) // if the plane corresponds to x (or y, or z) = a, where x (or y, or z) is in [a, b]
00312                     subCells[cidx[4]]->m_aInterList.push_back(inter);
00313                 else // if the plane corresponds to x (or y, or z) = b, where x (or y, or z) is in [a, b]
00314                     subCells[cidx[5]]->m_aInterList.push_back(inter);
00315             }
00316             else // if larger than the second mid-point
00317             {
00318                 if (inter->pnt[nidx]-box[nidx<<1] < box[(nidx<<1)+1]-inter->pnt[nidx])
00319                     subCells[cidx[6]]->m_aInterList.push_back(inter);
00320                 else // if the plane corresponds to x (or y, or z) = b, where x (or y, or z) is in [a, b]
00321                     subCells[cidx[7]]->m_aInterList.push_back(inter);
00322             }
00323         }
00324     }
00325 
00326     // update the 'computed faces' flag for all subcells
00327     for (int i = 0; i < 8*MAX_FACEIDX; i++)
00328     {
00329         int nCellIdx=_cellidx[i], nInheritIdx = _inheritidx[i];
00330         // assert that a face never inherit from more than one parent cell's face
00331         assert((subCells[nCellIdx]->m_aFaceComputed[nInheritIdx]) == false);
00332         /*if ((subCells[nCellIdx]->m_aFaceComputed[nInheritIdx]))
00333             std::cout << "assertion failed at depth " << cell->m_nDepth << std::endl;*/
00334         if (cell->m_aFaceComputed[nInheritIdx]) // if the parent cell's face this cell inherited from has 'computed' flag
00335             subCells[nCellIdx]->m_aFaceComputed[nInheritIdx] = true;
00336     }
00337 
00338 
00339     // split the representations for each surface
00340     for (int i = 0; i < 5; i++)
00341     {
00342         MPOL_T splited[8];
00343         int vr_idx[3] = {VARIDX(cell->m_F[i], 0), VARIDX(cell->m_F[i], 1), VARIDX(cell->m_F[i], 2)};
00344 
00345         subCells[7]->m_F[i] = cell->m_F[i];
00346 #ifdef SYNAPS_SBD3D_RNDFACTORWHENSUBDIVISION
00347         // subdivide along z direction with a perturbation
00348         mpoldse::split(subCells[3]->m_F[i].rep(), subCells[7]->m_F[i].rep(), vr_idx[2], t[2]);
00349 
00350         // subdivide along y direction with a perturbation
00351         mpoldse::split(subCells[1]->m_F[i].rep(), subCells[3]->m_F[i].rep(), vr_idx[1], t[1]);
00352         mpoldse::split(subCells[5]->m_F[i].rep(), subCells[7]->m_F[i].rep(), vr_idx[1], t[1]);
00353 
00354         // subdivide along x direction with a perturbation
00355         mpoldse::split(subCells[0]->m_F[i].rep(), subCells[1]->m_F[i].rep(), vr_idx[0], t[0]);
00356         mpoldse::split(subCells[2]->m_F[i].rep(), subCells[3]->m_F[i].rep(), vr_idx[0], t[0]);
00357         mpoldse::split(subCells[4]->m_F[i].rep(), subCells[5]->m_F[i].rep(), vr_idx[0], t[0]);
00358         mpoldse::split(subCells[6]->m_F[i].rep(), subCells[7]->m_F[i].rep(), vr_idx[0], t[0]);
00359 #else
00360         // subdivide along z direction
00361         mpoldse::split(subCells[3]->m_F[i].rep(), subCells[7]->m_F[i].rep(), vr_idx[2]);
00362 
00363         // subdivide along y direction
00364         mpoldse::split(subCells[1]->m_F[i].rep(), subCells[3]->m_F[i].rep(), vr_idx[1]);
00365         mpoldse::split(subCells[5]->m_F[i].rep(), subCells[7]->m_F[i].rep(), vr_idx[1]);
00366 
00367         // subdivide along x direction
00368         mpoldse::split(subCells[0]->m_F[i].rep(), subCells[1]->m_F[i].rep(), vr_idx[0]);
00369         mpoldse::split(subCells[2]->m_F[i].rep(), subCells[3]->m_F[i].rep(), vr_idx[0]);
00370         mpoldse::split(subCells[4]->m_F[i].rep(), subCells[5]->m_F[i].rep(), vr_idx[0]);
00371         mpoldse::split(subCells[6]->m_F[i].rep(), subCells[7]->m_F[i].rep(), vr_idx[0]);
00372 #endif
00373     }
00374 
00375     // compute the intersections on the cell borders (Dynamic Programming applied here
00376     // to avoid repetitive computation, see arrangements of _cellidx/_faceidx/_shareidx)
00377     for (int i = 0; i < 8*MAX_FACEIDX; i++)
00378     {
00379         intersectionlist_t aInters;
00380         int nCellIdx=_cellidx[i], nFaceIdx = _faceidx[i], nShareIdx = _shareidx[i];
00381 
00382         if (nFaceIdx >= 0)
00383         {
00384             computeinter_face(aInters, *subCells[nCellIdx], nFaceIdx>>1, nFaceIdx%2);
00385             // the cell shares a face with this cell (i.e. subcells[_cellidx[i]]) should also
00386             // share the some intersections (on a particular face) with this cell
00387             cell->mergelist(subCells[nCellIdx]->m_aInterList, aInters);
00388             cell->mergelist(subCells[nShareIdx]->m_aInterList, aInters);
00389             // mark the this face (i.e. face nFaceIdx of cell nCellIdx) as 'computed'
00390             subCells[nCellIdx]->m_aFaceComputed[nFaceIdx] = true;
00391             // for two subcells that shares a face, the shared face should of coz be parallel,
00392             // so their face indices should have same rounded value when devided by 3 (since
00393             // the face indices are (YZ_Plane_0, YZ_Plane_1, XZ_Plane0, XZ_Plane1, XY_Plane0,
00394             // XY_Plane1)). And of coz, they have different mod against 2.
00395             // by a simple deduction, f(x)=x+(x%2==1?-1:1) gives the correct face index for face x;
00396             // now we mark the shared face also as 'computed'
00397             assert(subCells[nShareIdx]->m_aFaceComputed[nFaceIdx+(nFaceIdx%2?-1:1)] == false);
00398             //if (subCells[nShareIdx]->m_aFaceComputed[nFaceIdx+(nFaceIdx%2?-1:1)])
00399             //    std::cout << "assertion failed at depth " << cell->m_nDepth << " nCell " << nCellIdx << " nFace " << nFaceIdx << " nShare " << nShareIdx << " i"  << i << ", cell " << nShareIdx << " face " << (nFaceIdx+(nFaceIdx%2?-1:1)) << " has already been computed." << std::endl;
00400             subCells[nShareIdx]->m_aFaceComputed[nFaceIdx+(nFaceIdx%2?-1:1)] = true;
00401         }
00402     }
00403 }
00404 //====================================================================
00405 template <class C, class CELL, class MSLV >
00406 void TopSbdCurve3d<C, CELL, MSLV>::computeintersections(intersectionlist_t& inters, CELL& c)
00407 {
00408   // for each facet of the cell, compute its intersections with the curve
00409   for (int v = 0; v < 3; v++) // for each axis
00410     {
00411       for (int f = 0; f < 2; f++) // for each facet perpendicular to that axis
00412         {
00413           // compute the intersections on this facet if it has not yet been computed
00414           if (!c.m_aFaceComputed[v*2+f])
00415             {
00416               computeinter_face(inters, c, v, f);
00417               c.m_aFaceComputed[v*2+f] = true;
00418             }
00419         }
00420     }
00421 }
00422 //====================================================================
00423 template < class C, class CELL, class MSLV >
00424 void TopSbdCurve3d<C, CELL, MSLV>::computeinter_face(intersectionlist_t& inters, 
00425                                                     CELL& c, 
00426                                                     int nAxis, int nFace)
00427 {
00428     point_t l(3u);
00429 
00430     // obtain the polynomail representing the face
00431     MPOL_T facepol[2]; // = {m_F[0], m_F[1]};
00432     face(facepol[0].rep(), c.m_F[0].rep(), nAxis, nFace);
00433     face(facepol[1].rep(), c.m_F[1].rep(), nAxis, nFace);
00434 
00435     // build a system of equations
00436     std::list<MPOL_T> eqnList;
00437     eqnList.push_back(facepol[0]);
00438     eqnList.push_back(facepol[1]);
00439 
00440     // solve the equation using a 2-variate solver
00441     MSLV solver(facepol, facepol+2, 1e-4);
00442     std::vector<C> results;
00443     int nSol, nVar;
00444 
00445     solver.solve(results, nSol, nVar, eqnList.begin(), eqnList.end());
00446     assert(nVar==2);
00447 
00448     // rescale the solution to the coordinates of the cell
00449     // also update the multiplicity accordingly
00450     for (int s = 0; s < nSol; s++)
00451     {
00452         C x[3] = {(C)nFace,(C)nFace,(C)nFace}; // the 3D coord of the point
00453 
00454         switch(nAxis)
00455         {
00456         case PLANEYZ: // yz-planes
00457             x[1] = (C)(results[nVar*s]);
00458             x[2] = (C)(results[nVar*s+1]);
00459             break;
00460         case PLANEXZ: // xz-planes
00461             x[0] = (C)(results[nVar*s]);
00462             x[2] = (C)(results[nVar*s+1]);
00463             break;
00464         case PLANEXY: // xy-planes
00465             x[0] = (C)(results[nVar*s]);
00466             x[1] = (C)(results[nVar*s+1]);
00467             break;
00468         }
00469 
00470         // determine the multiplicity
00471         //int nMaxDeg = c.m_T[nAxis].rep().env.sz(VARIDX(c.m_T[nAxis], nAxis)); 
00472         // the degree of nAxis_th variable of nAxis_th component of t = {t1, t2, t3} = gradF cross gradG
00473         int nMulti = 1;
00474 
00475         // currently, only account for multiplicity 1 or 2...
00476         //for (int v = 0; v< nMaxDeg; v++)
00477         //{
00478             C val[3] = {-1,-1,-1};
00479             eval(val[0], c.m_T[0], x); // evaluate the t_nAxis(x0, y0, z0) at the intersection (x0, y0, z0)
00480             eval(val[1], c.m_T[1], x); // evaluate the t_nAxis(x0, y0, z0) at the intersection (x0, y0, z0)
00481             eval(val[2], c.m_T[2], x); // evaluate the t_nAxis(x0, y0, z0) at the intersection (x0, y0, z0)
00482 
00483             /*if (abs(val[nAxis]) < 0.1)
00484             {
00485                 std::cout << "sin val: " << val[nAxis];
00486                 val[nAxis] /= sqrt(val[0]*val[0]+val[1]*val[1]+val[2]*val[2]); // normalize the vector
00487                 std::cout << " nor: " << val[nAxis] << std::endl;
00488             }
00489             else
00490             {
00491                 val[nAxis] /= sqrt(val[0]*val[0]+val[1]*val[1]+val[2]*val[2]); // normalize the vector
00492             }*/
00493 
00494             val[nAxis] /= sqrt(val[0]*val[0]+val[1]*val[1]+val[2]*val[2]); // normalize the vector
00495 
00496             if (_is_verysmall(std::abs(val[nAxis])))
00497             {
00498                 nMulti = 2;
00499                 //std::cout << "multiple point: val =" << val << ",  f(" << nAxis << ")=" << c.m_T[nAxis] << std::endl;
00500             }
00501         //}
00502 
00503 
00504         // scale to the cell's coordinate
00505         l[0] = c.m_cBox.x0*((C)1-x[0]) + c.m_cBox.x1*x[0];
00506         l[1] = c.m_cBox.y0*((C)1-x[1]) + c.m_cBox.y1*x[1];
00507         l[2] = c.m_cBox.z0*((C)1-x[2]) + c.m_cBox.z1*x[2];
00508 
00509         intersection_t *ni = new intersection_t(l, nAxis, nMulti);
00510         inters.push_back(ni);
00511         ni->nDepth = c.m_nDepth;
00512     }
00513 }
00514 
00515 //====================================================================
00516 template< class C, class CELL, class MSLV>
00517 template< class POINTGRAPH>
00518 void TopSbdCurve3d<C, CELL, MSLV>::run(POINTGRAPH& g, CELL* rc)
00519 {
00520     // initialization
00521     node_t* root = m_cOcTree.GetRoot();
00522 
00523     // step0:
00524    
00525     computeintersections(rc->m_aInterList,*rc);
00526     std::cout << "no. of roots at root: " << rc->m_aInterList.size() << std::endl;
00527     root->SetCell(rc);
00528 
00529     // step3:
00530     // subdivision
00531     std::list<node_t*> aNodeQueue;
00532     std::vector<node_t*> aSingular;
00533     pointgraph_int_t g_tmp;
00534     int nCellProc = 0;
00535 
00536     aNodeQueue.push_front(root);
00537 
00538     while(aNodeQueue.size())
00539     {
00540         // process the first node in the queue
00541         node_t* n = *aNodeQueue.begin();
00542         aNodeQueue.pop_front();
00543 
00544         bool bSubdivide = false;
00545         bool bBuildGraph = false;
00546 
00547         // if the cell is regular, i.e. the topology of the surface within the cell is certified
00548         if (regularity(*n->GetCell()))
00549         {
00550           if (is_variety_in(*n->GetCell())) // it may contains the polar curve
00551             {
00552                 if (n->GetCell()->m_nDepth < m_nMinSbd) // if we did less than the minimum number of subdivision
00553                 {
00554                     // subdivide further
00555                     bSubdivide = true;
00556                 }
00557                 else // build the graph
00558                 {
00559                     bBuildGraph = true;
00560                 }
00561             }
00562             else // if regular yet does not contain the surface, of no interest for further subdivision
00563             {
00564                 // its safe to say that not any facet of this cell should contain a intersection point
00565                 memset(n->GetCell()->m_aFaceComputed, true, sizeof(bool)*6);
00566             }
00567         }
00568         else // if the cell is not regular
00569         {
00570             if (n->GetCell()->testsize(m_dEpsilon)) // if the cell size is smaller than epsilon
00571             {
00572                 // keep a count of singular cells
00573               if (is_variety_in(*n->GetCell()))
00574                 {
00575                     bBuildGraph = true;
00576                     aSingular.push_back(n);
00577                 }
00578             }
00579             else // if the cell size is still larger than the set threshold
00580             {
00581                 // further subdivision if it may contains the intersection
00582               if (is_variety_in(*n->GetCell()))
00583                     bSubdivide = true;
00584             }
00585         }
00586 
00587         // in fact, bSubdivide and bBuildGraph cannot be both true (but can be both false)
00588         // assert(!(bSubdivide && bBuildGraph))
00589         if (bSubdivide)
00590         {
00591             // append the subdivided node at the back of the queue
00592             // in doing this, we can guarantee that larger nodes are processed first
00593             std::vector<node_t*> aNodeList;
00594             subdivide(n, aNodeList);
00595             for (int i = 0; i < (int)aNodeList.size(); i++)
00596                 aNodeQueue.push_front(aNodeList[i]);
00597         }
00598         else if (bBuildGraph)
00599         {
00600           build_graph(g_tmp,*n->GetCell());
00601         }
00602 
00603         // a counter
00604         nCellProc++;
00605     }
00606 
00607     // convert g_tmp to g
00608     // calcuate the size
00609     int nVert = (int)rep(rep(g_tmp).Vertices).size();
00610     int nEdges = 0;
00611 
00612     for(int i=0; i<(int)rep(rep(g_tmp).Edges).size(); ++i)
00613         nEdges += (int)rep(rep(g_tmp).Edges)[i].size();
00614 
00615     if (nVert)
00616     {
00617         int nVertIdx = 0;
00618         int nBase = (int)rep(rep(g).Vertices).size();; // base index of g when g is not empty
00619         m_nCount++;
00620         for(int i=0; i<nVert; ++i)
00621         {
00622             if (rep(rep(rep(g_tmp).Vertices)[i])->nCount != m_nCount)
00623             {
00624                 rep(rep(rep(g_tmp).Vertices)[i])->nCount = m_nCount;
00625                 rep(rep(rep(g_tmp).Vertices)[i])->nFlag = nVertIdx++;
00626                 g.push_back(rep(rep(rep(g_tmp).Vertices)[i])->pnt);
00627             }
00628         }
00629 
00630         for(int i=0; i<(int)rep(rep(g_tmp).Edges).size(); ++i)
00631         {
00632             for(int j=0; j<(int)rep(rep(g_tmp).Edges)[i].size(); ++j)
00633             {
00634                 g.insert(rep(rep(rep(g_tmp).Vertices)[i])->nFlag+nBase,
00635                          rep(rep(rep(g_tmp).Vertices)[rep(rep(g_tmp).Edges)[i][j]])->nFlag+nBase);
00636             }
00637         }
00638 
00639         std::cout << "No. of Vertices: " << (int)rep(rep(g).Vertices).size() - nBase << std::endl;
00640         std::cout << "No. of Edges: " << nEdges << std::endl;
00641     }
00642 
00643     std::cout << "Total Cell Processed: " << nCellProc << std::endl;
00644     std::cout << "Sing Cells: " << aSingular.size() << std::endl;
00645 }
00646 
00647 
00648 //--------------------------------------------------------------------
00649 namespace topology {
00650 
00651   template<class C, class U, class CELL, class SLV>
00652   void assign(point_graph<C>& g, 
00653               const MPol<U>& p1, 
00654               const MPol<U>& p2, 
00655               TopSbdCurve3d<C,CELL,SLV> slv,
00656               double x0, double x1,
00657               double y0, double y1,
00658               double z0, double z1)
00659   {
00660     typedef TopSbdCurve3d<C,CELL,SLV> MTH;
00661     typedef typename MTH::cell_t  cell_t;
00662     typedef typename MTH::mpolb_t mpolb_t;
00663     
00664     // step1:
00665     // computing grad F1 cross grad F2
00666     mpolb_t p1dse, p2dse;
00667 
00668     toBernstein(p1, p1dse, x0,x1,y0,y1,z0,z1);
00669     toBernstein(p2, p2dse, x0,x1,y0,y1,z0,z1);
00670 
00671     // debug only
00672 //     std::cout << "F("
00673 //               << degree(p1,0) << ", "
00674 //               << degree(p1,1) << ", "
00675 //               << degree(p1,2) << "): "
00676 //               << p1 << std::endl;
00677 
00678 //     std::cout << "G("
00679 //               << degree(p2,0) << ", "
00680 //               << degree(p2,1) << ", "
00681 //               << degree(p2,2) << "): "
00682 //               << p2 << std::endl;
00683 
00684 //     std::cout << "Fb("
00685 //               << p1dse.rep().env.sz(VARIDX(p1dse, 0)) << ", "
00686 //               << p1dse.rep().env.sz(VARIDX(pd1dse, 1)) << ", "
00687 //               << p1dse.rep().env.sz(VARIDX(pd1dse, 2)) << "): "
00688 //               << p1dse << std::endl;
00689 
00690 //     std::cout << "Gb("
00691 //               << p2dse.rep().env.sz(VARIDX(pd2dse, 0)) << ", "
00692 //               << p2dse.rep().env.sz(VARIDX(pd2dse, 1)) << ", "
00693 //               << p2dse.rep().env.sz(VARIDX(pd2dse, 2)) << "): "
00694 //               << p2dse << std::endl;
00695 
00696     mpolb_t p1_dx, p1_dy, p1_dz;
00697     mpolb_t p2_dx, p2_dy, p2_dz;
00698 
00699     diff(p1_dx, p1dse, VARIDX(pd1dse, 0));
00700     diff(p1_dy, p1dse, VARIDX(pd1dse, 1));
00701     diff(p1_dz, p1dse, VARIDX(pd1dse, 2));
00702 
00703     diff(p2_dx, p2dse, VARIDX(pd1dse, 0));
00704     diff(p2_dy, p2dse, VARIDX(pd2dse, 1));
00705     diff(p2_dz, p2dse, VARIDX(pd2dse, 2));
00706 
00707 //     std::cout << "grad Fb = {" << p1_dx << ", " <<  std::endl
00708 //               << "           " << p1_dy << ", " <<  std::endl
00709 //               << "           " << p1_dz << "} " <<  std::endl;
00710 
00711 //     std::cout << "grad Gb = {" << p2_dx << ", " <<  std::endl
00712 //               << "           " << p2_dy << ", " <<  std::endl
00713 //              << "           " << p2_dz << "} " <<  std::endl;
00714 
00715     mpolb_t c1, c2, c3; // grad F1 cross grad F2
00716     using arithm::operator*;
00717     c1 = p1_dy * p2_dz - p1_dz * p2_dy;
00718     c2 = p2_dx * p1_dz - p2_dz * p1_dx;
00719     c3 = p1_dx * p2_dy - p1_dy * p2_dx;
00720 
00721     // debug only
00722 //     std::cout << "grad(F) X grad(G) = {" << std::endl
00723 //               << "                     " << c1 << "," << std::endl
00724 //               << "                     " << c2 << "," << std::endl
00725 //               << "                     " << c3 << "}" << std::endl;
00726 
00727     // step2:
00728     // assign to the octree
00729     cell_t* rc = new cell_t(x0,x1,y0,y1,z0,z1);
00730     rc->m_F[0] = p1dse;
00731     rc->m_F[1] = p2dse;
00732     rc->m_T[0] = c1;
00733     rc->m_T[1] = c2;
00734     rc->m_T[2] = c3;
00735 
00736     slv.run(g,rc);
00737   }
00738 }
00739 
00740 //====================================================================
00741 template < class C >
00742 bool is_variety_in(Cell3dN<5,C>& c)
00743 {
00744   return c.containsVariety(0) && c.containsVariety(1);
00745 }
00746 
00747 //--------------------------------------------------------------------
00748 template < class POINTGRAPH, class CELL>
00749 void build_graph(POINTGRAPH& g, CELL& c)
00750 {
00751 
00752   typedef typename CELL::coeff_t         coeff_t;
00753   typedef typename CELL::MPOL_T          MPOL_T;
00754   typedef intersection3d<coeff_t>        intersection_t;
00755   typedef topology::point<coeff_t>       point_t;
00756   typedef shared_object<intersection_t*> intersection_ptr;
00757   // determine the regularity if it is not known
00758   if (c.m_nRegularIdx == -2)
00759     regularity(c);
00760   
00761   if (c.m_nRegularIdx == -1) // a singular cell
00762     {
00763       if (c.m_aInterList.size()) // if the cell has other point connecting to this cell
00764         {
00765         buildsingular:
00766           ;//     build_box(g,c.m_cBox);
00767         }
00768     }
00769     else // a regular cell
00770     {
00771         Seq<intersection_t*> inters = c.m_aInterList;
00772 
00773         // compute the intersections of the curve and the cell facets
00774         // computeintersections(inters); // no longer necessary
00775 
00776         // sort according to the regular direction
00777         switch(c.m_nRegularIdx)
00778         {
00779         case 0:
00780             sort(inters.begin(), inters.end(), XMAJOR());
00781             break;
00782         case 1:
00783             sort(inters.begin(), inters.end(), YMAJOR());
00784             break;
00785         case 2:
00786             sort(inters.begin(), inters.end(), ZMAJOR());
00787             break;
00788         }
00789 
00790         // construct the graph
00791         bool bConnectWithLast = false;
00792         int nMulti = 1;
00793         int nLastIdx = -1;
00794         for (int v = 0; v < (int)inters.size(); v++)
00795         {
00796             int nBIdx = g.push_back(intersection_ptr(inters[v]));
00797 
00798             if (bConnectWithLast)
00799                 g.insert(nLastIdx, nBIdx);
00800 
00801             if (!bConnectWithLast && nMulti < inters[v]->nMulti)
00802                 bConnectWithLast = false;
00803             else
00804                 bConnectWithLast = !bConnectWithLast;
00805 
00806             if (nMulti < inters[v]->nMulti)
00807             {
00808                 nMulti++;
00809                 v--;
00810             }
00811             else
00812             {
00813                 nMulti = 1;
00814             }
00815 
00816             nLastIdx = nBIdx;
00817         }
00818 
00819         if (inters.size() > 2)
00820         {
00821           //            std::cout << "a cell with " << inters.size() << " intersection points." << std::endl;
00822 //             for (int i = 0; i < inters.size(); i++)
00823 //             {
00824 //                 std::cout << "(" << inters[i]->pnt[0] << ", " << inters[i]->pnt[1]
00825 //                           << ", " << inters[i]->pnt[2] << ") d=" << inters[i]->nDepth << std::endl;
00826 //             }
00827             goto buildsingular;
00828         }
00829     }
00830 }
00831 //--------------------------------------------------------------------
00832 template < class C >
00833 bool regularity(Cell3dN<5,C>& c)
00834 {
00835     // check if the curve tangent, i.e. t = (t0, t1, t2) = grad F cross grad G
00836     // has no sign change for at least one of its component (or say, along one axis)
00837     for (int i = 0; i < 3; i++)
00838     {
00839         C* coeff = c.m_T[i].begin();
00840         C lastC = coeff[0];
00841         bool bSignChange = false;
00842 
00843         for (int j = 1; j < c.m_T[i].rep().env.sz(); j++)
00844         {
00845             C cirC = coeff[j];
00846             if (_is_singular(cirC)) // if cirC is a zero value
00847             {
00848                 continue;
00849             }
00850             else
00851             {
00852                 if (lastC*cirC < 0) // a sign change is detected
00853                 {
00854                     bSignChange = true;
00855                     lastC = cirC;
00856                     break;
00857                 }
00858 
00859                 lastC = cirC; // always remember the last non-singular coefficient
00860             }
00861         }
00862 
00863         if (_is_singular(lastC))  //all coefficients are singular (or 0)
00864             continue;
00865         else if (!bSignChange) // no sign changes
00866         {
00867             c.m_nRegularIdx = i;
00868             return TRUE;
00869         }
00870     }
00871     c.m_nRegularIdx = -1;
00872     return FALSE;
00873 }
00874 //--------------------------------------------------------------------
00875 template<typename C, typename MPOL, typename MPOLB>
00876 void toBernstein(const MPOL& p, MPOLB& bp, 
00877                  const C& a0, const C& b0, 
00878                  const C& a1, const C& b1, 
00879                  const C& a2, const C& b2)
00880 {
00881     MPOLB tmp(p);
00882     //    PRINT_DEBUG(p);
00883     mpoldse::restrict(tmp.rep(),VARIDX(tmp,0),a0,b0); // RESCALE(p,v,a,b) ==> p = p((1-t)*a+t*b)
00884     mpoldse::restrict(tmp.rep(),VARIDX(tmp,1),a1,b1); // RESCALE(p,v,a,b) ==> p = p((1-t)*a+t*b)
00885     mpoldse::restrict(tmp.rep(),VARIDX(tmp,2),a2,b2); // RESCALE(p,v,a,b) ==> p = p((1-t)*a+t*b)
00886 
00887     bp = tmp;
00888 }
00889 //--------------------------------------------------------------------
00890 template < class POINTGRAPH, class CELL>
00891 void build_box(POINTGRAPH& g, CELL* c)
00892 {
00893   typedef typename CELL::coeff_t         coeff_t;
00894   typedef topology::point<coeff_t>       point_t;
00895   int nBIdx;
00896   point_t l(3u);
00897   
00898   // form the graph representing the shape of the cell
00899   for (int i = 0; i < 8; i++)
00900     {
00901       l[0] = ((i&0x1) ? c->m_cBox.x1:c->m_cBox.x0);
00902       l[1] = ((i&0x2) ? c->m_cBox.y1:c->m_cBox.y0);
00903       l[2] = ((i&0x4) ? c->m_cBox.z1:c->m_cBox.z0);
00904       //nBIdx = g.push_back(shared_object(l));
00905       //nBIdx = g.push_back(intersection_ptr(new intersection_t(l,0,1)));
00906       nBIdx = g.push_back(l);
00907     }
00908   
00909   nBIdx -= 7;
00910   
00911   g.insert(nBIdx,nBIdx+1);
00912   g.insert(nBIdx,nBIdx+2);
00913   g.insert(nBIdx+1,nBIdx+3);
00914   g.insert(nBIdx+2,nBIdx+3);
00915   
00916   g.insert(nBIdx+4,nBIdx+5);
00917   g.insert(nBIdx+4,nBIdx+6);
00918   g.insert(nBIdx+5,nBIdx+7);
00919   g.insert(nBIdx+6,nBIdx+7);
00920   
00921   g.insert(nBIdx,nBIdx+4);
00922   g.insert(nBIdx+1,nBIdx+5);
00923   g.insert(nBIdx+2,nBIdx+6);
00924   g.insert(nBIdx+3,nBIdx+7);
00925 }
00926 //--------------------------------------------------------------------
00927 __END_NAMESPACE_SYNAPS
00928 /********************************************************************/
00929 #endif // synaps_topology_TopSbdCurve3d_h
00930 

SYNAPS DOCUMENTATION
logo