Simplified the usage of traits' functors in NN_coordinates_3

There is no need to define local versions of functors that are already available
through the traits / kernel.

Also do not default construct traits.
This commit is contained in:
Mael Rouxel-Labbé 2017-04-18 12:43:46 +02:00
parent b638cb4d66
commit b5233c93d7
1 changed files with 159 additions and 171 deletions

View File

@ -41,32 +41,18 @@ namespace CGAL {
template <class Gt>
typename Gt::FT
compute_squared_distance(const typename Gt::Point_3 &p, const typename Gt::Point_3 &q);
template <class Gt>
typename Gt::FT
compute_signed_volume(const typename Gt::Point_3 &p, const typename Gt::Point_3 &q,
const typename Gt::Point_3 &r, const typename Gt::Point_3 &s);
// positive when s is on the positive side of the plane defined
// by p, q, and r
template <class Gt>
typename Gt::FT
compute_squared_area(const typename Gt::Point_3 &p, const typename Gt::Point_3 &q,
const typename Gt::Point_3 &r);
template <class Gt>
typename Gt::FT
signed_area(const typename Gt::Point_3 &p, const typename Gt::Point_3 &q,
const typename Gt::Point_3 &r, const typename Gt::Point_3 &point_of_vue);
//signed area of the triangle determined by p q r
signed_area(const typename Gt::Point_3& p, const typename Gt::Point_3& q,
const typename Gt::Point_3& r, const typename Gt::Point_3& point_of_view,
const Gt& gt = Gt());
// ====================== Delaunay Triangulation utilities ==========================
// === Declarations
template < class DT>
typename DT::Geom_traits::Point_3
construct_circumcenter(const typename DT::Facet &f,const typename DT::Geom_traits::Point_3 &Q);
construct_circumcenter(const typename DT::Facet& f,
const typename DT::Geom_traits::Point_3& Q,
const typename DT::Geom_traits& gt = typename DT::Geom_traits());
// ====================== Natural Neighbors Querries ==========================
// === Definitions
@ -125,51 +111,56 @@ laplace_natural_neighbor_coordinates_3(const Dt& dt,
std::map<Vertex_handle,Coord_type> coordinate;
typename std::map<Vertex_handle,Coord_type>::iterator coor_it;
for (bound_it=bound_facets.begin();
bound_it!=bound_facets.end(); ++bound_it)
{//for each facet on the boundary
Facet f1=*bound_it;
Cell_handle cc1=f1.first;
if (dt.is_infinite(cc1))
return make_triple(nn_out,norm_coeff=Coord_type(1), false);//point outside the convex-hull
CGAL_triangulation_assertion_code(Cell_handle cc2=cc1->neighbor(f1.second);)
CGAL_triangulation_assertion(std::find(cells.begin(),cells.end(),cc1)!=cells.end());//TODO : Delete
CGAL_triangulation_assertion(std::find(cells.begin(),cells.end(),cc2)==cells.end());//TODO : Delete
Point C_1 = construct_circumcenter<Dt>(f1,Q);
for(int j=1;j<4;j++)
{//for each vertex P of the boundary facet
Vertex_handle vP=cc1->vertex((f1.second+j)&3);
Vertex_handle vR=cc1->vertex(dt.next_around_edge(f1.second,(f1.second+j)&3));
// turn around the oriented edge vR vP
Cell_handle cc3=cc1;
int num_next=dt.next_around_edge((f1.second+j)&3,f1.second);
Cell_handle next=cc3->neighbor(num_next);
while (std::find(cells.begin(),cells.end(),next)!=cells.end())
{
CGAL_triangulation_assertion( next != cc1 );
cc3=next;
num_next=dt.next_around_edge(cc3->index(vR),cc3->index(vP));
next=cc3->neighbor(num_next);
}
Point C_3=construct_circumcenter<Dt>(Facet(cc3,num_next),Q);
Point midPQ = midpoint(vP->point(),Q);
Coord_type coor_add = signed_area<Gt>(C_3,C_1,midPQ, vP->point());
((coor_it=coordinate.find(vP))==coordinate.end())?
coordinate[vP]=coor_add : coor_it->second+=coor_add;// Replace by a function call
}
}//end : for each facet on the boundary
for (bound_it = bound_facets.begin(); bound_it != bound_facets.end(); ++bound_it)
{
//for each facet on the boundary
Facet f1 = *bound_it;
Cell_handle cc1 = f1.first;
if (dt.is_infinite(cc1))
return make_triple(nn_out,norm_coeff=Coord_type(1), false);//point outside the convex-hull
norm_coeff=0;
for (coor_it = coordinate.begin();
coor_it != coordinate.end();
++coor_it)
CGAL_triangulation_assertion_code(Cell_handle cc2 = cc1->neighbor(f1.second);)
CGAL_triangulation_assertion(std::find(cells.begin(),cells.end(),cc1) != cells.end());//TODO : Delete
CGAL_triangulation_assertion(std::find(cells.begin(),cells.end(),cc2) == cells.end());//TODO : Delete
Point C_1 = construct_circumcenter<Dt>(f1, Q, dt.geom_traits());
for(int j=1; j<4; j++)
{
Coord_type co = coor_it->second/
(CGAL_NTS sqrt(compute_squared_distance<Gt>(coor_it->first->point(),Q)));
*nn_out++= std::make_pair(coor_it->first,co);
norm_coeff+=co;
//for each vertex P of the boundary facet
Vertex_handle vP = cc1->vertex((f1.second+j)&3);
Vertex_handle vR = cc1->vertex(dt.next_around_edge(f1.second,(f1.second+j)&3));
// turn around the oriented edge vR vP
Cell_handle cc3 = cc1;
int num_next = dt.next_around_edge((f1.second+j)&3,f1.second);
Cell_handle next = cc3->neighbor(num_next);
while (std::find(cells.begin(),cells.end(),next) != cells.end())
{
CGAL_triangulation_assertion( next != cc1 );
cc3 = next;
num_next = dt.next_around_edge(cc3->index(vR),cc3->index(vP));
next = cc3->neighbor(num_next);
}
Point C_3 = construct_circumcenter<Dt>(Facet(cc3,num_next), Q, dt.geom_traits());
Point midPQ = midpoint(vP->point(),Q);
Coord_type coor_add = signed_area<Gt>(C_3,C_1,midPQ, vP->point(), dt.geom_traits());
((coor_it = coordinate.find(vP)) == coordinate.end())?
coordinate[vP] = coor_add : coor_it->second += coor_add; // Replace by a function call
}
return make_triple(nn_out,norm_coeff,true);
} //end : for each facet on the boundary
norm_coeff = 0;
for (coor_it=coordinate.begin(); coor_it!=coordinate.end(); ++coor_it)
{
Coord_type co = coor_it->second /
(CGAL_NTS sqrt(dt.geom_traits().compute_squared_distance_3_object()(
coor_it->first->point(),Q)));
*nn_out++ = std::make_pair(coor_it->first,co);
norm_coeff += co;
}
return make_triple(nn_out, norm_coeff, true);
}
template <class Dt, class OutputIterator>
@ -217,73 +208,92 @@ sibson_natural_neighbor_coordinates_3(const Dt& dt,
typename std::map<Vertex_handle,Coord_type>::iterator coor_it;
for (cit = cells.begin(); cit != cells.end(); ++cit)
{// for each cell cc1 in conflict
Cell_handle cc1=*cit;
CGAL_triangulation_assertion(std::find(cells.begin(),cells.end(),cc1)!=cells.end());//TODO : Delete
if (dt.is_infinite(cc1))
return make_triple(nn_out,norm_coeff=Coord_type(1), false);//point outside the convex-hull
Point C1 = dt.dual(cc1);
for(int i=0;i<4;i++)
{//for each neighboring cell cc2 of cc1
Cell_handle cc2=cc1->neighbor(i);
if(std::find(cells.begin(),cells.end(),cc2)==cells.end())
{// cc2 outside the conflict cavity
Point C_1 = construct_circumcenter<Dt>(Facet(cc1,i),Q);
for(int j=1;j<4;j++)
{//for each vertex P of the boundary facet
Vertex_handle vP=cc1->vertex((i+j)&3);//&3 in place of %4
Vertex_handle vR=cc1->vertex(dt.next_around_edge(i,(i+j)&3));
// turn around the oriented edge vR vP
Cell_handle cc3=cc1;
int num_next=dt.next_around_edge((i+j)&3,i);
Cell_handle next=cc3->neighbor(num_next);
while (std::find(cells.begin(),cells.end(),next)!=cells.end())
{ //next is in conflict
CGAL_triangulation_assertion( next != cc1 );
cc3=next;
num_next=dt.next_around_edge(cc3->index(vR),cc3->index(vP));
next=cc3->neighbor(num_next);
}
if (dt.is_infinite(cc3))
return make_triple(nn_out,norm_coeff=Coord_type(1), false);//point outside the convex-hull
Point C3=dt.dual(cc3);
Point C_3=construct_circumcenter<Dt>(Facet(cc3,num_next),Q);
Point midPQ = midpoint(vP->point(),Q);
Point midPR = midpoint(vP->point(),vR->point());
Coord_type coor_add = compute_signed_volume<Gt>(C_1,C1,midPR,midPQ);
coor_add -= compute_signed_volume<Gt>(C_1,C_3,midPR,midPQ);
coor_add += compute_signed_volume<Gt>(C3,C_3,midPR,midPQ);
((coor_it=coordinate.find(vP))==coordinate.end())?
coordinate[vP]=coor_add : coor_it->second+=coor_add;// Replace by a function call
}
}
else // cc2 in the conflict cavity
{
CGAL_triangulation_assertion(std::find(cells.begin(),cells.end(),cc2)!=cells.end());//TODO : Delete
if (dt.is_infinite(cc2))
return make_triple(nn_out,norm_coeff=Coord_type(1), false);//point outside the convex-hull
Point C2=dt.dual(cc2);
for(int j=1;j<4;j++)
{//for each vertex P of the internal facet
Vertex_handle vP=cc1->vertex((i+j)&3);
Vertex_handle vR=cc1->vertex(dt.next_around_edge(i,(i+j)&3));
Point midPQ = midpoint(vP->point(),Q);
Point midPR = midpoint(vP->point(),vR->point());
Coord_type coor_add = compute_signed_volume<Gt>(C2,C1,midPR,midPQ);
((coor_it=coordinate.find(vP))==coordinate.end())?
coordinate[vP]=coor_add : coor_it->second+=coor_add;// Replace by a function call
}
}
}
}
norm_coeff=0;
for (coor_it = coordinate.begin();
coor_it != coordinate.end();
++coor_it)
{
// for each cell cc1 in conflict
Cell_handle cc1 = *cit;
CGAL_triangulation_assertion(std::find(cells.begin(),cells.end(),cc1)!=cells.end());//TODO : Delete
if (dt.is_infinite(cc1))
return make_triple(nn_out,norm_coeff=Coord_type(1), false);//point outside the convex-hull
typename Dt::Geom_traits::Compute_volume_3 vol =
dt.geom_traits().compute_volume_3_object();
Point C1 = dt.dual(cc1);
for(int i=0; i<4; i++)
{
*nn_out++= std::make_pair(coor_it->first,coor_it->second);
norm_coeff+=coor_it->second;
//for each neighboring cell cc2 of cc1
Cell_handle cc2 = cc1->neighbor(i);
if(std::find(cells.begin(),cells.end(),cc2) == cells.end())
{
// cc2 outside the conflict cavity
Point C_1 = construct_circumcenter<Dt>(Facet(cc1,i), Q, dt.geom_traits());
for(int j=1; j<4; j++)
{
//for each vertex P of the boundary facet
Vertex_handle vP = cc1->vertex((i+j)&3);//&3 in place of %4
Vertex_handle vR = cc1->vertex(dt.next_around_edge(i,(i+j)&3));
// turn around the oriented edge vR vP
Cell_handle cc3 = cc1;
int num_next = dt.next_around_edge((i+j)&3,i);
Cell_handle next = cc3->neighbor(num_next);
while (std::find(cells.begin(),cells.end(),next) != cells.end())
{ //next is in conflict
CGAL_triangulation_assertion( next != cc1 );
cc3 = next;
num_next = dt.next_around_edge(cc3->index(vR),cc3->index(vP));
next = cc3->neighbor(num_next);
}
if (dt.is_infinite(cc3))
{
//point outside the convex-hull
return make_triple(nn_out,norm_coeff = Coord_type(1), false);
}
Point C3 = dt.dual(cc3);
Point C_3 = construct_circumcenter<Dt>(Facet(cc3,num_next), Q, dt.geom_traits());
Point midPQ = midpoint(vP->point(),Q);
Point midPR = midpoint(vP->point(),vR->point());
Coord_type coor_add = vol(C_1,C1,midPR,midPQ);
coor_add -= vol(C_1,C_3,midPR,midPQ);
coor_add += vol(C3,C_3,midPR,midPQ);
((coor_it = coordinate.find(vP)) == coordinate.end())?
coordinate[vP] = coor_add : coor_it->second += coor_add;// Replace by a function call
}
}
else // cc2 in the conflict cavity
{
CGAL_triangulation_assertion(std::find(cells.begin(),cells.end(),cc2)!=cells.end());//TODO : Delete
if (dt.is_infinite(cc2))
{
//point outside the convex-hull
return make_triple(nn_out,norm_coeff = Coord_type(1), false);
}
Point C2 = dt.dual(cc2);
for(int j=1;j<4;j++)
{
//for each vertex P of the internal facet
Vertex_handle vP=cc1->vertex((i+j)&3);
Vertex_handle vR=cc1->vertex(dt.next_around_edge(i,(i+j)&3));
Point midPQ = midpoint(vP->point(),Q);
Point midPR = midpoint(vP->point(),vR->point());
Coord_type coor_add = vol(C2,C1,midPR,midPQ);
((coor_it=coordinate.find(vP))==coordinate.end())?
coordinate[vP]=coor_add : coor_it->second+=coor_add;// Replace by a function call
}
}
}
}
norm_coeff=0;
for (coor_it=coordinate.begin(); coor_it!=coordinate.end(); ++coor_it)
{
*nn_out++ = std::make_pair(coor_it->first,coor_it->second);
norm_coeff += coor_it->second;
}
return make_triple(nn_out,norm_coeff,true);
}
@ -319,39 +329,14 @@ bool is_correct_natural_neighborhood(const Dt& /*dt*/,
template <class Gt>
typename Gt::FT
compute_squared_distance(const typename Gt::Point_3 &p, const typename Gt::Point_3 &q)
signed_area(const typename Gt::Point_3& p, const typename Gt::Point_3& q,
const typename Gt::Point_3& r,
const typename Gt::Point_3& point_of_view,
const Gt& gt /* = Gt() */)
{
return Gt().compute_squared_distance_3_object()(p,q);
}
template <class Gt>
typename Gt::FT
compute_signed_volume(const typename Gt::Point_3 &p, const typename Gt::Point_3 &q,
const typename Gt::Point_3 &r, const typename Gt::Point_3 &s)
{
return Gt().compute_volume_3_object()
(Gt().construct_tetrahedron_3_object()(p,q,r,s));
}// positive when s is on the positive side of the plane defined
// by p, q, and r
template <class Gt>
typename Gt::FT
compute_squared_area(const typename Gt::Point_3 &p, const typename Gt::Point_3 &q,
const typename Gt::Point_3 &r)
{
return Gt().compute_squared_area_3_object()
(Gt().construct_triangle_3_object()
(p,q,r));
}
template <class Gt>
typename Gt::FT
signed_area(const typename Gt::Point_3 &p, const typename Gt::Point_3 &q,
const typename Gt::Point_3 &r, const typename Gt::Point_3 &point_of_vue)
//signed area of the triangle determined by p q r
{
return sqrt(compute_squared_area<Gt>(p,q,r))
*(orientation(p, q, r, point_of_vue) == COUNTERCLOCKWISE?+1:-1);
// signed area of the triangle p q r
return gt.compute_area_3_object()(p,q,r)
* (gt.orientation_3_object()(p, q, r, point_of_view) == COUNTERCLOCKWISE?+1:-1);
}
// ====================== Delaunay Triangulation utilities ==========================
@ -359,21 +344,24 @@ signed_area(const typename Gt::Point_3 &p, const typename Gt::Point_3 &q,
template < class DT>
typename DT::Geom_traits::Point_3
construct_circumcenter(const typename DT::Facet &f,const typename DT::Geom_traits::Point_3 &Q)
{
construct_circumcenter(const typename DT::Facet& f,
const typename DT::Geom_traits::Point_3& Q,
const typename DT::Geom_traits& gt /* = typename DT::Geom_traits() */ )
{
CGAL_triangulation_precondition(//&3 in place of %4
!coplanar(f.first->vertex((f.second+1)&3)->point(),
f.first->vertex((f.second+2)&3)->point(),
f.first->vertex((f.second+3)&3)->point(),
Q));
!gt.coplanar_3_object()(
f.first->vertex((f.second+1)&3)->point(),
f.first->vertex((f.second+2)&3)->point(),
f.first->vertex((f.second+3)&3)->point(),
Q));
// else the facet is not on the enveloppe of the conflict cavity associated to P
return typename DT::Geom_traits().construct_circumcenter_3_object()
(f.first->vertex((f.second+1)&3)->point(),
f.first->vertex((f.second+2)&3)->point(),
f.first->vertex((f.second+3)&3)->point(),
Q);
return gt.construct_circumcenter_3_object()(
f.first->vertex((f.second+1)&3)->point(),
f.first->vertex((f.second+2)&3)->point(),
f.first->vertex((f.second+3)&3)->point(),
Q);
}
} //namespace CGAL
#endif // CGAL_NATURAL_NEIGHBORS_3_H