mirror of
https://github.com/cosmo-sims/MUSIC.git
synced 2024-09-19 17:03:46 +02:00
added support for convex hull based high-res regions
This commit is contained in:
parent
0c5cb001e9
commit
d3e48b3ec9
3 changed files with 440 additions and 1 deletions
293
plugins/convex_hull.hh
Normal file
293
plugins/convex_hull.hh
Normal file
|
@ -0,0 +1,293 @@
|
|||
#ifndef CONVEX_HULL_HH
|
||||
#define CONVEX_HULL_HH
|
||||
|
||||
#include <vector>
|
||||
#include <set>
|
||||
#include <cmath>
|
||||
|
||||
#include <omp.h>
|
||||
|
||||
#include "log.hh"
|
||||
|
||||
/***** Slow but short convex hull Implementation ******/
|
||||
/*
|
||||
* Finds the convex hull of a set of data points.
|
||||
* Very simple implementation using the O(nh) gift-wrapping algorithm
|
||||
* Does not properly take care of degeneracies or round-off problems,
|
||||
* in that sense, only 'approximate' convex hull.
|
||||
*
|
||||
* adapted and expanded from the 'ch3quad' code by Timothy Chan
|
||||
* (https://cs.uwaterloo.ca/~tmchan)
|
||||
*/
|
||||
|
||||
template< typename real_t >
|
||||
struct convex_hull{
|
||||
typedef const real_t *cpr_;
|
||||
|
||||
size_t npoints_;
|
||||
std::vector<int> faceidx_L_, faceidx_U_;
|
||||
std::vector<real_t> normals_L_, normals_U_;
|
||||
std::vector<real_t> x0_L_, x0_U_;
|
||||
real_t centroid_[3], volume_;
|
||||
real_t left_[3], right_[3];
|
||||
|
||||
inline double turn( cpr_ p, cpr_ q, cpr_ r ) const
|
||||
{ return (q[0]-p[0])*(r[1]-p[1]) - (r[0]-p[0])*(q[1]-p[1]); }
|
||||
|
||||
template< bool islower >
|
||||
inline double orient( cpr_ p, cpr_ q, cpr_ r, cpr_ s ) const
|
||||
{
|
||||
if( islower )
|
||||
return (q[2]-p[2])*turn(p,r,s) - (r[2]-p[2])*turn(p,q,s) + (s[2]-p[2])*turn(p,q,r);
|
||||
|
||||
return (p[2]-q[2])*turn(p,r,s) - (p[2]-r[2])*turn(p,q,s) + (p[2]-s[2])*turn(p,q,r);
|
||||
}
|
||||
|
||||
inline real_t dot( cpr_ x, cpr_ y ) const
|
||||
{
|
||||
return x[0]*y[0]+x[1]*y[1]+x[2]*y[2];
|
||||
}
|
||||
|
||||
inline real_t det3x3( cpr_ a ) const
|
||||
{
|
||||
return (a[0]*(a[4]*a[8]-a[7]*a[5])
|
||||
+ a[1]*(a[5]*a[6]-a[8]*a[3])
|
||||
+ a[2]*(a[3]*a[7]-a[6]*a[4]));
|
||||
}
|
||||
|
||||
void compute_face_normals( cpr_ points )
|
||||
{
|
||||
normals_L_.assign( faceidx_L_.size(), 0.0 );
|
||||
normals_U_.assign( faceidx_U_.size(), 0.0 );
|
||||
x0_L_.assign( faceidx_L_.size(), 0.0 );
|
||||
x0_U_.assign( faceidx_U_.size(), 0.0 );
|
||||
|
||||
#pragma omp parallel for
|
||||
for( int i=0; i<(int)faceidx_L_.size()/3; ++i )
|
||||
{
|
||||
real_t d1[3], d2[3];
|
||||
for( int j=0; j<3; ++j )
|
||||
{
|
||||
x0_L_[3*i+j] = points[ 3*faceidx_L_[3*i+0] + j ];
|
||||
d1[j] = points[ 3*faceidx_L_[3*i+1] + j ] - points[ 3*faceidx_L_[3*i+0] + j ];
|
||||
d2[j] = points[ 3*faceidx_L_[3*i+2] + j ] - points[ 3*faceidx_L_[3*i+0] + j ];
|
||||
}
|
||||
|
||||
normals_L_[3*i+0] = d1[1]*d2[2] - d1[2]*d2[1];
|
||||
normals_L_[3*i+1] = d1[2]*d2[0] - d1[0]*d2[2];
|
||||
normals_L_[3*i+2] = d1[0]*d2[1] - d1[1]*d2[0];
|
||||
|
||||
// negative sign for lower hull
|
||||
double norm_n = -sqrt(normals_L_[3*i+0]*normals_L_[3*i+0]+
|
||||
normals_L_[3*i+1]*normals_L_[3*i+1]+
|
||||
normals_L_[3*i+2]*normals_L_[3*i+2]);
|
||||
|
||||
normals_L_[3*i+0] /= norm_n;
|
||||
normals_L_[3*i+1] /= norm_n;
|
||||
normals_L_[3*i+2] /= norm_n;
|
||||
|
||||
|
||||
}
|
||||
|
||||
#pragma omp parallel for
|
||||
for( int i=0; i<(int)faceidx_U_.size()/3; ++i )
|
||||
{
|
||||
real_t d1[3], d2[3];
|
||||
for( int j=0; j<3; ++j )
|
||||
{
|
||||
x0_U_[3*i+j] = points[ 3*faceidx_U_[3*i+0] + j ];
|
||||
d1[j] = points[ 3*faceidx_U_[3*i+1] + j ] - points[ 3*faceidx_U_[3*i+0] + j ];
|
||||
d2[j] = points[ 3*faceidx_U_[3*i+2] + j ] - points[ 3*faceidx_U_[3*i+0] + j ];
|
||||
}
|
||||
|
||||
normals_U_[3*i+0] = d1[1]*d2[2] - d1[2]*d2[1];
|
||||
normals_U_[3*i+1] = d1[2]*d2[0] - d1[0]*d2[2];
|
||||
normals_U_[3*i+2] = d1[0]*d2[1] - d1[1]*d2[0];
|
||||
|
||||
double norm_n = sqrt(normals_U_[3*i+0]*normals_U_[3*i+0]+
|
||||
normals_U_[3*i+1]*normals_U_[3*i+1]+
|
||||
normals_U_[3*i+2]*normals_U_[3*i+2]);
|
||||
|
||||
normals_U_[3*i+0] /= norm_n;
|
||||
normals_U_[3*i+1] /= norm_n;
|
||||
normals_U_[3*i+2] /= norm_n;
|
||||
}
|
||||
}
|
||||
|
||||
void compute_center( cpr_ points )
|
||||
{
|
||||
real_t xc[3] = {0.0,0.0,0.0};
|
||||
real_t xcp[3] = {0.0,0.0,0.0};
|
||||
|
||||
real_t totvol = 0.0;
|
||||
|
||||
for( size_t i=0; i<3*npoints_; ++i )
|
||||
xc[i%3] += points[i];
|
||||
|
||||
xc[0] /= npoints_;
|
||||
xc[1] /= npoints_;
|
||||
xc[2] /= npoints_;
|
||||
|
||||
|
||||
for( size_t i=0; i<faceidx_L_.size()/3; ++i )
|
||||
{
|
||||
real_t a[9];
|
||||
real_t xct[3] = {xc[0],xc[1],xc[2]};
|
||||
|
||||
for( size_t j=0; j<3; ++j )
|
||||
{
|
||||
for( size_t k=0; k<3; ++k )
|
||||
{
|
||||
xct[k] += points[3*faceidx_L_[3*i+j]+k];
|
||||
a[3*j+k] = points[3*faceidx_L_[3*i+j]+k]-xc[k];
|
||||
}
|
||||
}
|
||||
|
||||
xct[0] *= 0.25;
|
||||
xct[1] *= 0.25;
|
||||
xct[2] *= 0.25;
|
||||
|
||||
real_t vol = fabs(det3x3(a))/6.0;
|
||||
|
||||
totvol += vol;
|
||||
xcp[0] += vol * xct[0];
|
||||
xcp[1] += vol * xct[1];
|
||||
xcp[2] += vol * xct[2];
|
||||
}
|
||||
|
||||
for( size_t i=0; i<faceidx_U_.size()/3; ++i )
|
||||
{
|
||||
real_t a[9];
|
||||
real_t xct[3] = {xc[0],xc[1],xc[2]};
|
||||
|
||||
for( size_t j=0; j<3; ++j )
|
||||
{
|
||||
for( size_t k=0; k<3; ++k )
|
||||
{
|
||||
xct[k] += points[3*faceidx_U_[3*i+j]+k];
|
||||
a[3*j+k] = points[3*faceidx_U_[3*i+j]+k]-xc[k];
|
||||
}
|
||||
}
|
||||
|
||||
xct[0] *= 0.25;
|
||||
xct[1] *= 0.25;
|
||||
xct[2] *= 0.25;
|
||||
|
||||
real_t vol = fabs(det3x3(a))/6.0;
|
||||
|
||||
totvol += vol;
|
||||
xcp[0] += vol * xct[0];
|
||||
xcp[1] += vol * xct[1];
|
||||
xcp[2] += vol * xct[2];
|
||||
}
|
||||
|
||||
|
||||
volume_ = totvol;
|
||||
centroid_[0] = xcp[0] / totvol;
|
||||
centroid_[1] = xcp[1] / totvol;
|
||||
centroid_[2] = xcp[2] / totvol;
|
||||
}
|
||||
|
||||
template< bool islower >
|
||||
void wrap( cpr_ points, int i, int j, std::vector<int>& idx )
|
||||
{
|
||||
int k,l,m;
|
||||
int h = (int)idx.size()/3;
|
||||
|
||||
for( m=0; m<h; ++m )
|
||||
if( ( idx[3*m+0]==i && idx[3*m+1]==j ) ||
|
||||
( idx[3*m+1]==i && idx[3*m+2]==j ) ||
|
||||
( idx[3*m+2]==i && idx[3*m+0]==j ) )
|
||||
return;
|
||||
|
||||
for( k=i, l=0; l < (int)npoints_; ++l )
|
||||
if( turn(&points[3*i],&points[3*j],&points[3*l]) < 0 &&
|
||||
orient<islower>(&points[3*i],&points[3*j],&points[3*k],&points[3*l]) >=0 )
|
||||
k = l;
|
||||
|
||||
if( turn(&points[3*i],&points[3*j],&points[3*k]) >= 0 ) return;
|
||||
|
||||
idx.push_back( i );
|
||||
idx.push_back( j );
|
||||
idx.push_back( k );
|
||||
|
||||
wrap<islower>( points, k, j, idx );
|
||||
wrap<islower>( points, i, k, idx );
|
||||
}
|
||||
|
||||
template< typename T >
|
||||
bool check_point( const T* x ) const
|
||||
{
|
||||
for( size_t i=0; i<normals_L_.size()/3; ++i )
|
||||
{
|
||||
double xp[3] = {x[0]-x0_L_[3*i+0],x[1]-x0_L_[3*i+1],x[2]-x0_L_[3*i+2]};
|
||||
if( dot( xp, &normals_L_[3*i])<0.0 ) return false;
|
||||
}
|
||||
|
||||
for( size_t i=0; i<normals_U_.size()/3; ++i )
|
||||
{
|
||||
double xp[3] = {x[0]-x0_U_[3*i+0],x[1]-x0_U_[3*i+1],x[2]-x0_U_[3*i+2]};
|
||||
if( dot( xp, &normals_U_[3*i])<0.0 ) return false;
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
void get_defining_indices( std::set<int>& unique ) const
|
||||
{
|
||||
unique.clear();
|
||||
|
||||
for( size_t i=0; i<faceidx_L_.size(); ++i )
|
||||
unique.insert( faceidx_L_[i] );
|
||||
for( size_t i=0; i<faceidx_U_.size(); ++i )
|
||||
unique.insert( faceidx_U_[i] );
|
||||
}
|
||||
|
||||
convex_hull( cpr_ points, size_t npoints )
|
||||
: npoints_( npoints )
|
||||
{
|
||||
faceidx_L_.reserve(npoints_*3);
|
||||
faceidx_U_.reserve(npoints_*3);
|
||||
|
||||
size_t i,j,l;
|
||||
for( i=0, l=1; l<npoints_; ++l )
|
||||
if( points[3*i] > points[3*l] ) i=l;
|
||||
for( j=i, l=0; l<npoints_; ++l )
|
||||
if( i!=l && turn(&points[3*i],&points[3*j],&points[3*l]) >= 0 ) j=l;
|
||||
|
||||
int nt = omp_get_num_threads();
|
||||
omp_set_num_threads( std::min(2,omp_get_max_threads()) );
|
||||
|
||||
#pragma omp parallel for
|
||||
for( int thread=0; thread<2; ++thread )
|
||||
{
|
||||
if( thread==0 )
|
||||
wrap<true>( points, i, j, faceidx_L_ );
|
||||
if( thread==1 )
|
||||
wrap<false>( points, i, j, faceidx_U_ );
|
||||
}
|
||||
|
||||
omp_set_num_threads(nt);
|
||||
|
||||
compute_face_normals( points );
|
||||
compute_center( points );
|
||||
|
||||
// finally compute AABB
|
||||
left_[0] = left_[1] = left_[2] = 1e30;
|
||||
right_[0] = right_[1] = right_[2] = -1e30;
|
||||
|
||||
for( size_t q=0; q<npoints_; ++q )
|
||||
for( size_t p=0; p<3; ++p )
|
||||
{
|
||||
if( points[3*q+p] > right_[p] )
|
||||
right_[p] = points[3*q+p];
|
||||
if( points[3*q+p] < left_[p] )
|
||||
left_[p] = points[3*q+p];
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
|
||||
#endif // CONVEX_HULL_HH
|
|
@ -19,7 +19,8 @@ struct point_reader{
|
|||
return iss.eof() && !iss.fail();
|
||||
}
|
||||
|
||||
void read_points_from_file( std::string fname, float vfac_, std::vector<float>& p )
|
||||
template< typename real_t >
|
||||
void read_points_from_file( std::string fname, float vfac_, std::vector<real_t>& p )
|
||||
{
|
||||
std::ifstream ifs(fname.c_str());
|
||||
if( !ifs )
|
||||
|
|
145
plugins/region_convex_hull.cc
Normal file
145
plugins/region_convex_hull.cc
Normal file
|
@ -0,0 +1,145 @@
|
|||
#include <vector>
|
||||
#include <iostream>
|
||||
#include <cmath>
|
||||
#include <cassert>
|
||||
#include <fstream>
|
||||
#include <sstream>
|
||||
#include <cctype>
|
||||
#include <algorithm>
|
||||
|
||||
#include "region_generator.hh"
|
||||
#include "convex_hull.hh"
|
||||
#include "point_file_reader.hh"
|
||||
|
||||
|
||||
//! Convex hull region plugin
|
||||
class region_convex_hull_plugin : public region_generator_plugin{
|
||||
private:
|
||||
|
||||
convex_hull<double> *phull_;
|
||||
int shift[3], shift_level, padding_;
|
||||
double vfac_;
|
||||
bool do_extra_padding_;
|
||||
|
||||
void apply_shift( size_t Np, double *p, int *shift, int levelmin )
|
||||
{
|
||||
double dx = 1.0/(1<<levelmin);
|
||||
LOGINFO("unapplying shift of previous zoom region to region particles :\n" \
|
||||
"\t [%d,%d,%d] = (%f,%f,%f)",shift[0],shift[1],shift[2],shift[0]*dx,shift[1]*dx,shift[2]*dx);
|
||||
|
||||
for( size_t i=0,i3=0; i<Np; i++,i3+=3 )
|
||||
for( size_t j=0; j<3; ++j )
|
||||
p[i3+j] = p[i3+j]-shift[j]*dx;
|
||||
}
|
||||
|
||||
public:
|
||||
explicit region_convex_hull_plugin( config_file& cf )
|
||||
: region_generator_plugin( cf )
|
||||
{
|
||||
std::vector<double> pp;
|
||||
|
||||
vfac_ = cf.getValue<double>("cosmology","vfact");
|
||||
padding_ = cf.getValue<int>("setup","padding");
|
||||
|
||||
|
||||
std::string point_file = cf.getValue<std::string>("setup","region_point_file");
|
||||
|
||||
point_reader pfr;
|
||||
pfr.read_points_from_file( point_file, vfac_, pp );
|
||||
|
||||
if( cf.containsKey("setup","region_point_shift") )
|
||||
{
|
||||
std::string point_shift = cf.getValue<std::string>("setup","region_point_shift");
|
||||
sscanf( point_shift.c_str(), "%d,%d,%d", &shift[0],&shift[1],&shift[2] );
|
||||
unsigned point_levelmin = cf.getValue<unsigned>("setup","region_point_levelmin");
|
||||
|
||||
apply_shift( pp.size()/3, &pp[0], shift, point_levelmin );
|
||||
shift_level = point_levelmin;
|
||||
}
|
||||
|
||||
// compute the convex hull
|
||||
phull_ = new convex_hull<double>( &pp[0], pp.size()/3 );
|
||||
|
||||
//expand the ellipsoid by one grid cell
|
||||
|
||||
unsigned levelmax = cf.getValue<unsigned>("setup","levelmax");
|
||||
double dx = 1.0/(1ul<<levelmax);
|
||||
//phull_->expand( dx );
|
||||
|
||||
// output the center
|
||||
float c[3] = { phull_->centroid_[0], phull_->centroid_[1], phull_->centroid_[2] };
|
||||
LOGINFO("Region center from convex hull centroid determined at\n\t (%f,%f,%f)",c[0],c[1],c[2]);
|
||||
|
||||
//-----------------------------------------------------------------
|
||||
// when querying the bounding box, do we need extra padding?
|
||||
do_extra_padding_ = false;
|
||||
|
||||
// conditions should be added here
|
||||
{
|
||||
std::string output_plugin = cf.getValue<std::string>("output","format");
|
||||
if( output_plugin == std::string("grafic2") )
|
||||
do_extra_padding_ = true;
|
||||
}
|
||||
}
|
||||
|
||||
~region_convex_hull_plugin()
|
||||
{
|
||||
delete phull_;
|
||||
}
|
||||
|
||||
void get_AABB( double *left, double *right, unsigned level )
|
||||
{
|
||||
for( int i=0; i<3; ++i )
|
||||
{
|
||||
left[i] = phull_->left_[i];
|
||||
right[i] = phull_->right_[i];
|
||||
}
|
||||
double dx = 1.0/(1ul<<level);
|
||||
double pad = (double)(padding_+1) * dx;
|
||||
|
||||
if( ! do_extra_padding_ ) pad = 0.0;
|
||||
|
||||
double ext = sqrt(3)*dx + pad;
|
||||
|
||||
for( int i=0;i<3;++i )
|
||||
{
|
||||
left[i] -= ext;
|
||||
right[i] += ext;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
void update_AABB( double *left, double *right )
|
||||
{
|
||||
// we ignore this, the grid generator must have generated a grid that contains the ellipsoid
|
||||
// it might have enlarged it, but who cares...
|
||||
}
|
||||
|
||||
bool query_point( double *x )
|
||||
{ return phull_->check_point( x ); }
|
||||
|
||||
bool is_grid_dim_forced( size_t* ndims )
|
||||
{ return false; }
|
||||
|
||||
void get_center( double *xc )
|
||||
{
|
||||
xc[0] = phull_->centroid_[0];
|
||||
xc[1] = phull_->centroid_[1];
|
||||
xc[2] = phull_->centroid_[2];
|
||||
}
|
||||
|
||||
void get_center_unshifted( double *xc )
|
||||
{
|
||||
double dx = 1.0/(1<<shift_level);
|
||||
float c[3] = { phull_->centroid_[0], phull_->centroid_[1], phull_->centroid_[2] };
|
||||
xc[0] = c[0]+shift[0]*dx;
|
||||
xc[1] = c[1]+shift[1]*dx;
|
||||
xc[2] = c[2]+shift[2]*dx;
|
||||
|
||||
}
|
||||
};
|
||||
|
||||
namespace{
|
||||
region_generator_plugin_creator_concrete< region_convex_hull_plugin > creator("convex_hull");
|
||||
}
|
||||
|
Loading…
Reference in a new issue