diff --git a/cgogn/modeling/skeleton_sampling.h b/cgogn/modeling/skeleton_sampling.h index 6ea31254..bd001101 100644 --- a/cgogn/modeling/skeleton_sampling.h +++ b/cgogn/modeling/skeleton_sampling.h @@ -25,7 +25,7 @@ #define CGOGN_RENDERING_SKELETON_SAMPLER_H_ #include -// #include +#include namespace cgogn { @@ -199,7 +199,6 @@ class SkeletonSampler inline void inter_skeleton(const VEC3& P, const VEC3& Du, std::vector& I, SCALAR ds_max) { - // std::cout << "INTER! P=" << P.transpose() << std::endl; SCALAR d = eval_skeleton(P); SCALAR ds = SCALAR(0); SCALAR d_prec = SCALAR(0);; @@ -208,7 +207,6 @@ class SkeletonSampler VEC3 Q; while (ds < ds_max) { - // std::cout << "ds:" << ds << " < s_max:" << ds_max << std::endl; while ((d >= SCALAR(0)) && (ds < ds_max)) { ds_prec = ds; @@ -216,103 +214,166 @@ class SkeletonSampler Q = P + ds * Du; d_prec = d; d = eval_skeleton(Q); - // std::cout <<"d="<< d<< " -> ds: " << ds << " -> Q: " << Q.transpose() << std::endl; } if (ds >= ds_max) break; - // std::cout << "ENTER OBJ" << d<< std::endl; SCALAR a = (ds * d_prec - ds_prec * d)/(d_prec-d); I.push_back(P + Du * a); - // std::cout << "d: " << d_prec << " => " << d << "ds: " << ds_prec << " => " << ds << "=> a = " << a << std::endl; - // std::cout << "PUSH FRONT " << I.back().transpose() << std::endl; -// ds += Epsilon; while ((d < SCALAR(0)) && (ds < ds_max)) { - // std::cout << "IN OBJ" << d<< std::endl; ds_prec = ds; ds += std::max(Epsilon,-d); Q = P + ds * Du; d_prec = d; d = eval_skeleton(Q); - // std::cout << "-> ds: " << ds << " -> Q: " << Q.transpose() << " d: "< " << d << "ds: " << ds_prec << " => " << ds << "=> a = " << a << std::endl; - // std::cout << "PUSH BACK " << I.back().transpose() << std::endl; } } - // inline void sample(SCALAR step) - // { - // std::cout << "BB: [" << bb_min_.transpose() << " - " << bb_max_.transpose() << "]" << std::endl; - // Epsilon = step / SCALAR(10); - // samples_.clear(); - // samples_.reserve(8192); - // // Z - // for (SCALAR y = bb_min_[1] + step / 2; y < bb_max_[1] - step / 2; y += step) - // for (SCALAR x = bb_min_[0] + step / 2; x < bb_max_[0] - step / 2; x += step) - // inter_skeleton(VEC3{x, y, bb_min_[2]}, VEC3{0, 0, 1}, samples_, bb_max_[2] - bb_min_[2]); - - // // Y - // for (SCALAR x = bb_min_[0] + step / 2; x < bb_max_[0] - step / 2; x += step) - // for (SCALAR z = bb_min_[2] + step / 2; z < bb_max_[2] - step / 2; z += step) - // inter_skeleton(VEC3{x, bb_min_[1], z}, VEC3{0, 1, 0}, samples_, bb_max_[1] - bb_min_[1]); - - // // X - // for (SCALAR y = bb_min_[1] + step / 2; y < bb_max_[1] - step / 2; y += step) - // for (SCALAR z = bb_min_[2] + step / 2; z < bb_max_[2] - step / 2; z += step) - // inter_skeleton(VEC3{bb_min_[0], y, z}, VEC3{1, 0, 0}, samples_, bb_max_[0] - bb_min_[0]); - - // std::cout << samples_.size() << "points genrated" << std::endl; - // } - - void sample(SCALAR step, SCALAR epsi=SCALAR(0)) + //inline void sample(SCALAR step) + //{ + // std::cout << "BB: [" << bb_min_.transpose() << " - " << bb_max_.transpose() << "]" << std::endl; + // Epsilon = step / SCALAR(10); + // samples_.clear(); + // samples_.reserve(8192); + // // Z + // for (SCALAR y = bb_min_[1] + step / 2; y < bb_max_[1] - step / 2; y += step) + // for (SCALAR x = bb_min_[0] + step / 2; x < bb_max_[0] - step / 2; x += step) + // inter_skeleton(VEC3{x, y, bb_min_[2]}, VEC3{0, 0, 1}, samples_, bb_max_[2] - bb_min_[2]); + + // // Y + // for (SCALAR x = bb_min_[0] + step / 2; x < bb_max_[0] - step / 2; x += step) + // for (SCALAR z = bb_min_[2] + step / 2; z < bb_max_[2] - step / 2; z += step) + // inter_skeleton(VEC3{x, bb_min_[1], z}, VEC3{0, 1, 0}, samples_, bb_max_[1] - bb_min_[1]); + + // // X + // for (SCALAR y = bb_min_[1] + step / 2; y < bb_max_[1] - step / 2; y += step) + // for (SCALAR z = bb_min_[2] + step / 2; z < bb_max_[2] - step / 2; z += step) + // inter_skeleton(VEC3{bb_min_[0], y, z}, VEC3{1, 0, 0}, samples_, bb_max_[0] - bb_min_[0]); + + // std::cout << samples_.size() << "points genrated" << std::endl; + //} + + //void sample(SCALAR step, SCALAR epsi=SCALAR(0)) + //{ + // std::cout << "BB: [" << bb_min_.transpose() << " - " << bb_max_.transpose() << "]" << std::endl; + // Epsilon = epsi==SCALAR(0) ? step / SCALAR(20): epsi; + // std::vector samp1; + // samp1.reserve(8192); + // std::vector samp2; + // samp2.reserve(8192); + // std::vector samp3; + // samp3.reserve(8192); + + // std::random_device rd; // Will be used to obtain a seed for the random number engine + // std::mt19937 gen(rd()); + // std::uniform_real_distribution dis(SCALAR(0), step / SCALAR(2)); + + // // Z + // std::thread th1([&]() + // { + // for (SCALAR y = bb_min_[1] + step / 2; y < bb_max_[1] - step / 2; y += step) + // for (SCALAR x = bb_min_[0] + step / 2; x < bb_max_[0] - step / 2; x += step) + // inter_skeleton(VEC3{x + dis(gen), y + dis(gen), bb_min_[2]}, VEC3{0, 0, 1}, samp1, + // bb_max_[2] - bb_min_[2]); + // }); + + // // Y + // std::thread th2([&]() + // { + // for (SCALAR x = bb_min_[0] + step / 2; x < bb_max_[0] - step / 2; x += step) + // for (SCALAR z = bb_min_[2] + step / 2; z < bb_max_[2] - step / 2; z += step) + // inter_skeleton(VEC3{x + dis(gen), bb_min_[1], z + dis(gen)}, VEC3{0, 1, 0}, samp2, + // bb_max_[1] - bb_min_[1]); + // }); + + // // X + // std::thread th3([&]() + // { + // for (SCALAR y = bb_min_[1] + step / 2; y < bb_max_[1] - step / 2; y += step) + // for (SCALAR z = bb_min_[2] + step / 2; z < bb_max_[2] - step / 2; z += step) + // inter_skeleton(VEC3{bb_min_[0], y + dis(gen), z + dis(gen)}, VEC3{1, 0, 0}, samp3, + // bb_max_[0] - bb_min_[0]); + // }); + // th1.join(); + // th2.join(); + // th3.join(); + + // samples_.clear(); + // samples_.insert(samples_.end(),samp1.begin(),samp1.end()); + // samples_.insert(samples_.end(),samp2.begin(),samp2.end()); + // samples_.insert(samples_.end(),samp3.begin(),samp3.end()); + + // std::cout << samples_.size() << "points genrated" << std::endl; + //} + + + void sample(SCALAR step, SCALAR epsi = SCALAR(0)) { std::cout << "BB: [" << bb_min_.transpose() << " - " << bb_max_.transpose() << "]" << std::endl; - Epsilon = epsi==SCALAR(0) ? step / SCALAR(20): epsi; - std::vector samp1; - samp1.reserve(8192); - std::vector samp2; - samp2.reserve(8192); - std::vector samp3; - samp3.reserve(8192); - - // Z - std::thread th1([&]() - { - for (SCALAR y = bb_min_[1] + step / 2; y < bb_max_[1] - step / 2; y += step) - for (SCALAR x = bb_min_[0] + step / 2; x < bb_max_[0] - step / 2; x += step) - inter_skeleton(VEC3{x, y, bb_min_[2]}, VEC3{0, 0, 1}, samp1, bb_max_[2] - bb_min_[2]); - }); + Epsilon = epsi == SCALAR(0) ? step / SCALAR(10) : epsi; - // Y - std::thread th2([&]() - { - for (SCALAR x = bb_min_[0] + step / 2; x < bb_max_[0] - step / 2; x += step) - for (SCALAR z = bb_min_[2] + step / 2; z < bb_max_[2] - step / 2; z += step) - inter_skeleton(VEC3{x, bb_min_[1], z}, VEC3{0, 1, 0}, samp2, bb_max_[1] - bb_min_[1]); - }); + SCALAR s2 = step / 2; + VEC3 min_bb = bb_min_ + VEC3(s2, s2, s2); + VEC3 max_bb = bb_max_ - VEC3(s2, s2, s2); - // X - std::thread th3([&]() + int nbthr = std::thread::hardware_concurrency() / 3; + SCALAR nstep = step * nbthr; + + std::vector th(3 * nbthr, nullptr); + std::vector> samp(3 * nbthr); + int k = 0; + + for (int i = 0; i < nbthr; ++i) { - for (SCALAR y = bb_min_[1] + step / 2; y < bb_max_[1] - step / 2; y += step) - for (SCALAR z = bb_min_[2] + step / 2; z < bb_max_[2] - step / 2; z += step) - inter_skeleton(VEC3{bb_min_[0], y, z}, VEC3{1, 0, 0}, samp3, bb_max_[0] - bb_min_[0]); - }); - th1.join(); - th2.join(); - th3.join(); + samp[k].reserve(8192); + th[k] = new std::thread([this,i, k, step, nstep,min_bb,max_bb,&samp]() { + thread_local std::random_device rd; + thread_local std::mt19937 gen(rd()); + thread_local std::uniform_real_distribution dis(-step / SCALAR(2), step / SCALAR(2)); + for (SCALAR y = min_bb[1] + step * i; y < max_bb[1]; y += nstep) + for (SCALAR x = min_bb[0] + step / 2; x < max_bb[0]; x += step) + inter_skeleton(VEC3{x + dis(gen), y + dis(gen), min_bb[2]}, VEC3{0, 0, 1}, samp[k], + bb_max_[2] - bb_min_[2]); + }); + k++; + samp[k].reserve(8192); + th[k] = new std::thread([this,i, k, step, nstep, min_bb, max_bb, &samp]() { + thread_local std::random_device rd; + thread_local std::mt19937 gen(rd()); + thread_local std::uniform_real_distribution dis(-step / SCALAR(2), step / SCALAR(2)); + for (SCALAR x = min_bb[0] + step * i; x < max_bb[0]; x += nstep) + for (SCALAR z = min_bb[2] + step / 2; z < max_bb[2]; z += step) + inter_skeleton(VEC3{x + dis(gen), min_bb[1], z + dis(gen)}, VEC3{0, 1, 0}, samp[k], + bb_max_[1] - bb_min_[1]); + }); + k++; + samp[k].reserve(8192); + th[k] = new std::thread([this,i, k, step, nstep, min_bb, max_bb, &samp]() { + thread_local std::random_device rd; + thread_local std::mt19937 gen(rd()); + thread_local std::uniform_real_distribution dis(-step / SCALAR(2), step / SCALAR(2)); + for (SCALAR y = min_bb[1] + step * i; y < max_bb[1]; y += nstep) + for (SCALAR z = min_bb[2] + step / 2; z < max_bb[2]; z += step) + inter_skeleton(VEC3{min_bb[0], y + dis(gen), z + dis(gen)}, VEC3{1, 0, 0}, samp[k], + bb_max_[0] - bb_min_[0]); + }); + k++; + } + std::cout << k << " threads launched" << std::endl; samples_.clear(); - samples_.insert(samples_.end(),samp1.begin(),samp1.end()); - samples_.insert(samples_.end(),samp2.begin(),samp2.end()); - samples_.insert(samples_.end(),samp3.begin(),samp3.end()); + for (int i = 0; i < k; ++i) + { + th[i]->join(); + delete th[i]; + samples_.insert(samples_.end(), samp[i].begin(), samp[i].end()); + } std::cout << samples_.size() << "points genrated" << std::endl; } diff --git a/cgogn/rendering/apps/test_skel_shapes.cpp b/cgogn/rendering/apps/test_skel_shapes.cpp index ab08ee96..1d0dbf8c 100644 --- a/cgogn/rendering/apps/test_skel_shapes.cpp +++ b/cgogn/rendering/apps/test_skel_shapes.cpp @@ -59,6 +59,8 @@ class MySkelRender : public ui::ViewModule GLVec4 C1 = {-0.5f, 0.0f, 0.0f, 0.4f}; GLVec4 C2 = {0.4f, 1.0f, 0.2f, 0.7f}; GLVec4 C3 = {1.0f, -1.0f, -0.5f, 0.2f}; + + skel_drawer_.add_vertex(C00.topRows<3>(), C00[3]); // just to test version with radius parameter skel_drawer_.add_vertex(GLVec3{-2.0f, 0.0f, 0.4f}, 0.2f); skel_drawer_.add_vertex(C1); @@ -70,35 +72,8 @@ class MySkelRender : public ui::ViewModule skel_drawer_.add_edge(C2, C3); skel_drawer_.add_edge(C3, C1); skel_drawer_.add_triangle(C1, C2, C3); - - skel_drawer_.update(); - -// skel_sampler_.add_vertex(GLVec4(5, 5, 5, 2)); - //skel_sampler_.add_edge(GLVec3(0, 0, 0), 2, GLVec3(9, 0, 0), 4 ); - //for (float f = 0.0f; f<=11.0f; f+=0.5) - //std::cout << "D = " << skel_sampler_.eval_skeketon(GLVec3(f, 10, 0)) << std::endl; - // skel_sampler_.add_triangle(GLVec4(0, 0, 0, 2), GLVec4(4, 0, 0, 2.5), GLVec4(0, 4, 0, 3)); - - //skel_sampler_.add_triangle(GLVec4(0, 0, 0, 2), GLVec4(4, 0, 0, 2.65), GLVec4(0, 4, 0,3.1)); - //skel_sampler_.add_vertex(GLVec4(0, 0, 0, 2)); - //skel_sampler_.add_vertex(GLVec4(4, 0, 0, 2.65)); - //skel_sampler_.add_vertex(GLVec4(0, 4, 0, 3.1)); - - // for (int i=0; i<200; i++) - // { - // float alpha = 0.314f*i; - // float h = -2.0f+0.04f*i; - // float r = (i%2==1)?0.4f:0.3f; - // float rr = (i%2==0)?0.4f:0.3f; - // //skel_sampler_.add_vertex(GLVec4(7.0f*std::cos(alpha),7.0f*std::sin(alpha),h,r)); - // skel_sampler_.add_edge(GLVec4(1.0f*std::cos(alpha),1.0f*std::sin(alpha),h,r),GLVec4(2.5f*std::cos(alpha),2.5f*std::sin(alpha),h+2,rr)); - // } - - //skel_sampler_.add_vertex(GLVec4(0,0,0,1)); - // skel_sampler_.add_edge(GLVec4(0,0,0,1),GLVec4(5,3,2,2)); - // skel_sampler_.add_vertex(GLVec4(5,3,2,2)); skel_sampler_.add_vertex(C00.topRows<3>(), C00[3]); // just to test version with radius parameter skel_sampler_.add_vertex(GLVec3{-2.0f, 0.0f, 0.4f}, 0.2f); @@ -112,12 +87,33 @@ class MySkelRender : public ui::ViewModule skel_sampler_.add_edge(C3, C1); skel_sampler_.add_triangle(C1, C2, C3); + + // for (int i=0; i<100; i++) + //{ + // float alpha = 0.34f*i; + // float h = -2.0f+0.04f*i; + // float r = (i%2==1)?0.3f:0.2f; + // float rr = (i%2==0)?0.3f:0.2f; + // skel_sampler_.add_edge(GLVec4(1.0f * std::cos(alpha), 1.0f * std::sin(alpha), h, r), + // GLVec4(3.5f * std::cos(alpha), 3.5f * std::sin(alpha), h + 0.5, rr)); + // skel_sampler_.add_vertex(GLVec4(1.0f * std::cos(alpha), 1.0f * std::sin(alpha), h, r)); + // skel_sampler_.add_vertex(GLVec4(3.5f * std::cos(alpha), 3.5f * std::sin(alpha), h + 0.5, rr)); + + // skel_drawer_.add_edge(GLVec4(1.0f * std::cos(alpha), 1.0f * std::sin(alpha), h, r), + // GLVec4(3.5f * std::cos(alpha), 3.5f * std::sin(alpha), h + 0.5, rr)); + // skel_drawer_.add_vertex(GLVec4(1.0f * std::cos(alpha), 1.0f * std::sin(alpha), h, r)); + // skel_drawer_.add_vertex(GLVec4(3.5f * std::cos(alpha), 3.5f * std::sin(alpha), h + 0.5, rr)); + //} + //skel_drawer_.update(); + + + + // skel_sampler_.add_vertex(GLVec4(0, 0, 0, 1)); -// skel_sampler_.add_vertex(GLVec4(0, 0, 0, 1)); GLVec3 bbw = skel_sampler_.BBwidth(); - float step = std::min(std::min(bbw.x(),bbw.y()),bbw.z())/20; - skel_sampler_.sample(step,step/10); + float step = std::min(std::min(bbw.x(),bbw.y()),bbw.z())/200; + skel_sampler_.sample(step); // for(auto p:skel_sampler_.samples()) // std::cout << p.norm()<< std::endl;