Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

CG: optimize projection of point on a triangle #221

Merged
merged 3 commits into from
Nov 12, 2021
Merged

CG: optimize projection of point on a triangle #221

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
7dd89b6
CG: optimize projection of points on a triangle
andrea-iob Aug 28, 2021
88f8a6d
CG: remove unused function
andrea-iob Aug 28, 2021
b49eac7
CG: remove dependency from LAPACKe library
andrea-iob Aug 28, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -463,7 +463,7 @@ set(IO_EXTERNAL_DEPS "LibXml2;RapidJSON")
set(COMMUNICATIONS_EXTERNAL_DEPS "MPI")
set(LA_EXTERNAL_DEPS "PETSc")
set(SA_EXTERNAL_DEPS "")
set(CG_EXTERNAL_DEPS "LAPACKE")
set(CG_EXTERNAL_DEPS "")
set(PABLO_EXTERNAL_DEPS "")
set(PATCHKERNEL_EXTERNAL_DEPS "")
set(SURFUNSTRUCTURED_EXTERNAL_DEPS "")
Expand Down
94 changes: 17 additions & 77 deletions src/CG/CG_elem.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,8 +29,6 @@

# include <assert.h>

# include "bitpit_private_lapacke.hpp"

# include "bitpit_operators.hpp"

# include "CG.hpp"
Expand All @@ -50,6 +48,8 @@ namespace CGElem{
/*!
* \private
* Projects points on a triangle.
* Transcribed from Christer Ericson's Real-Time Collision Detection book,
* it's effectively Cramer's rule for solving a linear system.
* Projection points are the closest points to the original points within the triangle.
*
* \param[in] nPoints number of points
Expand All @@ -63,90 +63,30 @@ namespace CGElem{
*/
void _projectPointsTriangle( int nPoints, array3D const *points, array3D const &Q0, array3D const &Q1, array3D const &Q2, array3D *proj, double *lambda )
{

assert( validTriangle(Q0,Q1,Q2) );

array3D s0 = Q1-Q0;
array3D s1 = Q2-Q0;

std::array<double, 4> A = {{ dotProduct(s0,s0), 0, dotProduct(s0,s1), dotProduct(s1,s1) }};

const int MAX_STACK_POINTS = 8;
BITPIT_CREATE_WORKSPACE(B, double, 2 * nPoints, 2 * MAX_STACK_POINTS);

for( int i=0; i<nPoints; ++i){
array3D rP = *points -Q0;
B[2*i] = dotProduct(s0,rP);
B[2*i+1] = dotProduct(s1,rP);
++points;
}

int info = LAPACKE_dposv_work( LAPACK_COL_MAJOR, 'U', 2, nPoints, A.data(), 2, B, 2 );
assert( info == 0 );
BITPIT_UNUSED( info );

for( int i=0; i<nPoints; ++i){

double *b = &B[2*i];

lambda[0] = 1. -b[0] -b[1];
lambda[1] = b[0];
lambda[2] = b[1];

*proj = restrictPointTriangle( Q0, Q1, Q2, lambda);
lambda +=3;
proj += 1;
}
}

/*!
* \private
* Project points on a plane described by a triangle
*
* \param[in] nPoints number of points
* \param[in] points pointer to points' coordinates
* \param[in] Q0 first triangle vertex
* \param[in] Q1 second triangle vertex
* \param[in] Q2 third triangle vertex
* \param[out] proj pointer to the projection point;
* \param[out] lambda pointer to barycentric coordinates of projection points
* \return distances
*/
void _projectPointsPlane( int nPoints, array3D const *points, array3D const &Q0, array3D const &Q1, array3D const &Q2, array3D *proj, double *lambda )
{

assert( validTriangle(Q0,Q1,Q2) );

array3D s0 = Q1-Q0;
array3D s1 = Q2-Q0;
array3D v0 = Q1 - Q0;
array3D v1 = Q2 - Q0;

std::array<double, 4> A = {{ dotProduct(s0,s0), 0, dotProduct(s0,s1), dotProduct(s1,s1) }};

const int MAX_STACK_POINTS = 8;
BITPIT_CREATE_WORKSPACE(B, double, 2 * nPoints, 2 * MAX_STACK_POINTS);

for( int i=0; i<nPoints; ++i){
array3D rP = *points -Q0;
B[2*i] = dotProduct(s0,rP);
B[2*i+1] = dotProduct(s1,rP);
++points;
}
double d00 = dotProduct(v0, v0);
double d01 = dotProduct(v0, v1);
double d11 = dotProduct(v1, v1);

int info = LAPACKE_dposv_work( LAPACK_COL_MAJOR, 'U', 2, nPoints, A.data(), 2, B, 2 );
assert( info == 0 );
BITPIT_UNUSED( info );
double denom = d00 * d11 - d01 * d01;

for( int i=0; i<nPoints; ++i){
array3D v2 = points[i] - Q0;

double *b = &B[2*i];
double d20 = dotProduct(v2, v0);
double d21 = dotProduct(v2, v1);

lambda[0] = 1. -b[0] -b[1];
lambda[1] = b[0];
lambda[2] = b[1];
double *pointLambda = lambda + 3 * i;
pointLambda[1] = (d11 * d20 - d01 * d21) / denom;
pointLambda[2] = (d00 * d21 - d01 * d20) / denom;
pointLambda[0] = 1. - pointLambda[1] - pointLambda[2];

*proj = reconstructPointFromBarycentricTriangle( Q0, Q1, Q2, lambda);
lambda +=3;
proj += 1;
array3D *pointProjections = proj + i;
*pointProjections = restrictPointTriangle( Q0, Q1, Q2, pointLambda);
}
}

Expand Down
1 change: 0 additions & 1 deletion src/CG/CG_private.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,6 @@ namespace CGElem{
typedef std::array<double,3> array3D ;

void _projectPointsTriangle( int, array3D const *, array3D const &, array3D const &, array3D const &, array3D *, double *);
void _projectPointsPlane( int, array3D const *, array3D const &, array3D const &, array3D const &, array3D *, double *);
bool _intersectSegmentBox( array3D const &, array3D const &, array3D const &, array3D const &, bool, bool, std::vector<array3D> *, std::vector<int> *, int dim = 3, const double tolerance = DEFAULT_DISTANCE_TOLERANCE ) ;
bool _intersectPlaneBox( array3D const &, array3D const &, array3D const &, array3D const &, std::vector<array3D> *, int dim=3, const double tolerance = DEFAULT_DISTANCE_TOLERANCE );
bool _intersectBoxTriangle( array3D const &, array3D const &, array3D const &, array3D const &, array3D const &, bool, bool, bool, std::vector<array3D> *, std::vector<int> *, int dim=3, const double tolerance = DEFAULT_DISTANCE_TOLERANCE ) ;
Expand Down