24 #include <boost/geometry/index/rtree.hpp> 25 #include <boost/mpl/range_c.hpp> 26 #include <boost/mpl/for_each.hpp> 31 namespace sps_clustering {
33 namespace bmpl = boost::mpl;
34 namespace bgi = boost::geometry::index;
40 template<
size_t I,
typename Po
intType>
43 inline static void calc_distance(
size_t dim,
const PointType& point1,
const PointType& point2,
double& d2)
45 d2 += pow( boost::geometry::get<I>(point1) - boost::geometry::get<I>(point2), 2);
46 DCalc<I-1, PointType>::calc_distance(dim-1, point1, point2, d2);
50 template<
typename Po
intType>
51 struct DCalc<0, PointType>
53 inline static void calc_distance(
size_t,
const PointType& point1,
const PointType& point2,
double& d2)
55 d2 += pow( boost::geometry::get<0>(point1) - boost::geometry::get<0>(point2), 2);
60 template<
typename Po
intType>
61 void d2_calc(PointType
const& point_1, PointType
const& point_2,
double& result)
63 static constexpr std::size_t Rank = boost::geometry::dimension<PointType>::value;
64 DCalc<Rank-1, PointType>::calc_distance(Rank-1, point_1, point_2, result);
69 template<
size_t I,
typename Po
intType>
72 inline static void add_scalar(
size_t dim, PointType& point,
double const& c)
74 double new_coord = boost::geometry::get<I>(point) + c;
75 boost::geometry::set<I>(point, new_coord);
76 AddScalar<I-1, PointType>::add_scalar(dim-1, point,c);
80 template<
typename Po
intType>
81 struct AddScalar<0, PointType>
83 inline static void add_scalar(
size_t, PointType& point,
double const& c)
85 double new_coord = boost::geometry::get<0>(point) + c;
86 boost::geometry::set<0>(point, new_coord);
91 template<
typename Po
intType>
92 void add_scalar_to_point(PointType& point,
double const& value)
94 static constexpr std::size_t Rank = boost::geometry::dimension<PointType>::value;
95 AddScalar<Rank-1, PointType>::add_scalar(Rank-1, point, value);
99 template<
typename NumRepType>
103 std::vector<std::vector<size_t>> groups;
105 typedef std::pair<PointType, size_t> ValueType;
106 typedef bgi::rtree<ValueType, bgi::linear<16>> Tree;
108 std::vector< std::pair<PointType, size_t> > points;
109 points.reserve(cands.
size());
111 typedef typename data::SpCcl<NumRepType>::Dm Dm;
112 std::pair<Dm, Dm> dm_range = cands.
dm_range();
113 Dm dm_step=(dm_range.second - dm_range.first)/static_cast<typename Dm::value_type>(std::numeric_limits<std::size_t>::max());
114 if (dm_step == 0 * pss::astrotypes::units::parsecs_per_cube_cm)
115 dm_step = (1.0 * pss::astrotypes::units::parsecs_per_cube_cm );
117 if(cands.
tf_blocks().size() == 0 || cands.
size() == 0 )
return groups;
119 auto tsamp =
static_cast<boost::units::quantity<data::MilliSeconds,double>
>(((*(cands.
tf_blocks().begin()))->sample_interval())).value();
121 for(
size_t i = 0 ; i < cands.
size() ; ++i)
124 point.template set<0>(
static_cast<double>((cands[i].tstart()/tsamp)/(_config.time_tolerance()/tsamp)));
127 point.template set<1>(
static_cast<double>((cands[i].dm()).value())/((_config.dm_tolerance())).value());
129 point.template set<2>(
static_cast<double>(std::log2(cands[i].width().value()/tsamp)/std::log2(_config.pulse_width_tolerance().value()/tsamp)));
132 points.emplace_back(point, i);
135 Tree tree(points.begin(), points.end());
137 while( !tree.empty())
139 std::vector<ValueType> to_add;
141 to_add.push_back( *tree.qbegin( bgi::satisfies([](ValueType
const &){ return true; })) );
142 tree.remove( to_add.begin(), to_add.end() );
143 std::vector<ValueType> added;
145 for( std::size_t to_add_i = 0L; to_add_i < to_add.size(); ++to_add_i )
147 auto it = to_add.begin() + to_add_i;
150 PointType lower = it->first;
151 add_scalar_to_point(lower, -_linking_length);
152 PointType upper = it->first;
153 add_scalar_to_point(upper, _linking_length);
155 boost::geometry::model::box<PointType> box( lower, upper );
157 auto within_ball = [
this, &it](ValueType
const &v)
160 d2_calc(it->first, v.first, d2);
161 return d2 < _linking_length_2;
165 tree.query( bgi::within(box) && bgi::satisfies(within_ball), std::back_inserter(added) );
167 to_add.reserve(to_add.size() + added.size());
169 std::copy(added.begin(), added.end(), std::back_inserter(to_add));
172 tree.remove( added.begin(), added.end() );
181 std::vector< size_t > group;
182 for(
auto p : to_add )
184 group.push_back(p.second);
186 groups.push_back(group);
std::pair< Dm, Dm > const & dm_range() const
the max and minimum Dispersion Measure values of the candidates
BlocksType const & tf_blocks() const
Access the TimeFrequency blocks associated with the single pulse candidtes.
std::size_t size() const
Retrieve the size of the underlying vector.
Some limits and constants for FLDO.
std::vector< std::vector< size_t > > operator()(data::SpCcl< NumRepType > const &cands)
Group the candidates using the fof algorithm.