37 #ifndef GMM_DOMAIN_DECOMP_H__ 38 #define GMM_DOMAIN_DECOMP_H__ 50 template <
typename Matrix,
typename Po
int>
54 std::vector<Matrix> &vB) {
55 typedef typename linalg_traits<Matrix>::value_type value_type;
56 typedef abstract_null_type void_type;
57 typedef std::map<size_type, void_type> map_type;
59 size_type nbpts = pts.size();
60 if (!nbpts || pts[0].size() == 0) { vB.resize(0);
return; }
61 int dim = int(pts[0].size());
64 Point pmin = pts[0], pmax = pts[0];
65 for (size_type i = 1; i < nbpts; ++i)
66 for (
int k = 0; k < dim; ++k) {
67 pmin[k] = std::min(pmin[k], pts[i][k]);
68 pmax[k] = std::max(pmax[k], pts[i][k]);
71 std::vector<size_type> nbsub(dim), mult(dim);
72 std::vector<int> pts1(dim), pts2(dim);
73 size_type nbtotsub = 1;
74 for (
int k = 0; k < dim; ++k) {
75 nbsub[k] =
size_type((pmax[k] - pmin[k]) / msize)+1;
76 mult[k] = nbtotsub; nbtotsub *= nbsub[k];
79 std::vector<map_type> subs(nbtotsub);
81 std::vector<size_type> ns(dim), na(dim), nu(dim);
82 for (size_type i = 0; i < nbpts; ++i) {
83 for (
int k = 0; k < dim; ++k) {
84 register double a = (pts[i][k] - pmin[k]) / msize;
86 pts1[k] = int(a + overlap); pts2[k] = int(ceil(a-1.0-overlap));
91 for (
int k = 0; k < dim; ++k)
92 if ((ns[k] >= nbsub[k]) || (pts1[k] <
int(ns[k]))
93 || (pts2[k] >
int(ns[k]))) { ok =
false;
break; }
95 size_type ind = ns[0];
96 for (
int k=1; k < dim; ++k) ind += ns[k]*mult[k];
97 subs[ind][i] = void_type();
99 for (
int k = 0; k < dim; ++k) {
100 if (na[k] < 2) { na[k]++; ns[k]++; ++sum;
break; }
101 na[k] = 0; ns[k] -= 2; sum -= 2;
106 size_type nbmaxinsub = 0;
107 for (size_type i = 0; i < nbtotsub; ++i)
108 nbmaxinsub = std::max(nbmaxinsub, subs[i].size());
110 std::fill(ns.begin(), ns.end(),
size_type(0));
111 for (size_type i = 0; i < nbtotsub; ++i) {
112 if (subs[i].size() > 0 && subs[i].size() < nbmaxinsub / 10) {
114 for (
int k = 0; k < dim; ++k) nu[k] = ns[k];
115 size_type nbmax = 0, imax = 0;
117 for (
int l = 0; l < dim; ++l) {
119 for (
int m = 0; m < 2; ++m, nu[l]+=2) {
121 for (
int k = 0; k < dim && ok; ++k)
122 if (nu[k] >= nbsub[k]) ok =
false;
124 size_type ind = ns[0];
125 for (
int k=1; k < dim; ++k) ind += ns[k]*mult[k];
126 if (subs[ind].size() > nbmax)
127 { nbmax = subs[ind].size(); imax = ind; }
133 if (nbmax > subs[i].size()) {
134 for (map_type::iterator it=subs[i].begin(); it!=subs[i].end(); ++it)
135 subs[imax][it->first] = void_type();
139 for (
int k = 0; k < dim; ++k)
140 { ns[k]++;
if (ns[k] < nbsub[k])
break; ns[k] = 0; }
145 for (size_type i = 0; i < nbtotsub; ++i) {
146 if (subs[i].size() > 0)
147 {
if (i != effnb) std::swap(subs[i], subs[effnb]); ++effnb; }
153 for (size_type i = 0; i < effnb; ++i) {
154 clear(vB[i]);
resize(vB[i], nbpts, subs[i].size());
156 for (map_type::iterator it=subs[i].begin(); it!=subs[i].end(); ++it, ++j)
157 vB[i](it->first, j) = value_type(1);
size_t size_type
used as the common size type in the library
void rudimentary_regular_decomposition(std::vector< Point > pts, double msize, double overlap, std::vector< Matrix > &vB)
This function separates into small boxes of size msize with a ratio of overlap (in [0...
void resize(V &v, size_type n)
*/
Include the base gmm files.
void clear(L &l)
clear (fill with zeros) a vector or matrix.