arbor/fvm_layout.cpp
// Construct cv_geometry for cell from locset describing CV boundary points.
Construct cv_geometry for cell from locset describing CV boundary points.
Construct cv_geometry for cell from locset describing CV boundary points .
cv_geometry cv_geometry_from_ends(const cable_cell& cell, const locset& lset) {
    auto pop = [](auto& vec) { auto h = vec.back(); return vec.pop_back(), h; };
    cv_geometry geom;
    const auto& mp = cell.provider();
    const auto& m = mp.morphology();
    if (m.empty()) {
        geom.cell_cv_divs = {0, 0};
        return geom;
    }
    mlocation_list locs = thingify(lset, mp);
    auto neither_root_nor_terminal = [&m](mlocation x) {
        return !(x.pos==0 && x.branch==(m.branch_children(mnpos).size()>1u? mnpos: 0)) 
            && !(x.pos==1 && m.branch_children(x.branch).empty()); 
    };
    locs.erase(std::partition(locs.begin(), locs.end(), neither_root_nor_terminal), locs.end());
    util::sort(locs);
    util::unique_in_place(locs);
    constexpr fvm_index_type no_parent = -1;
    std::vector<std::pair<mlocation, fvm_index_type>> next_cv_head; 
    next_cv_head.emplace_back(mlocation{mnpos, 0}, no_parent);
    mcable_list cables;
    std::vector<msize_t> branches;
    geom.cv_cables_divs.push_back(0);
    fvm_index_type cv_index = 0;
    while (!next_cv_head.empty()) {
        auto next = pop(next_cv_head);
        mlocation h = next.first;
        cables.clear();
        branches.clear();
        branches.push_back(h.branch);
        geom.cv_parent.push_back(next.second);
        while (!branches.empty()) {
            msize_t b = pop(branches);
            auto it = locs.end();
            if (b!=mnpos && b==h.branch) {
                it = std::upper_bound(locs.begin(), locs.end(), h);
            }
            else if (b!=mnpos) {
                it = std::lower_bound(locs.begin(), locs.end(), mlocation{b, 0});
            }
            if (it!=locs.end() && it->branch==b) {
                cables.push_back({b, b==h.branch? h.pos: 0, it->pos});
                next_cv_head.emplace_back(*it, cv_index);
            }
            else {
                if (b!=mnpos) {
                    cables.push_back({b, b==h.branch? h.pos: 0, 1});
                }
                for (auto& c: m.branch_children(b)) {
                    branches.push_back(c);
                }
            }
        }
        sort(cables);
        append(geom.cv_cables, std::move(cables));
        geom.cv_cables_divs.push_back(geom.cv_cables.size());
        ++cv_index;
    }
    auto n_cv = cv_index;
    arb_assert(n_cv>0);
    arb_assert(geom.cv_parent.front()==-1);
    arb_assert(util::all_of(util::subrange_view(geom.cv_parent, 1, n_cv),
            [](auto v) { return v!=no_parent; }));
    assign(geom.cv_children, util::make_span(1, n_cv));
    stable_sort_by(geom.cv_children, [&geom](auto cv) { return geom.cv_parent[cv]; });
    geom.cv_children_divs.reserve(n_cv+1);
    geom.cv_children_divs.push_back(0);
    auto b = geom.cv_children.begin();
    auto e = geom.cv_children.end();
    auto from = b;
    for (fvm_index_type cv = 0; cv<n_cv; ++cv) {
        from = std::partition_point(from, e,
            [cv, &geom](auto i) { return geom.cv_parent[i]<=cv; });
        geom.cv_children_divs.push_back(from-b);
    }
    geom.cv_to_cell.assign(cv_index, 0);
    geom.cell_cv_divs = {0, (fvm_index_type)cv_index};
    geom.branch_cv_map.resize(1);
    std::vector<pw_elements<fvm_size_type>>& bmap = geom.branch_cv_map.back();
    for (auto cv: util::make_span(n_cv)) {
        for (auto cable: geom.cables(cv)) {
            if (cable.branch>=bmap.size()) {
                bmap.resize(cable.branch+1);
            }
            bmap[cable.branch].push_back(cable.prox_pos, cable.dist_pos, cv);
        }
    }
    return geom;
}
cv_geometry cv_geometry_from_ends ( const cable_cell & cell , const locset & lset ) { auto pop = [ ] ( auto & vec ) { auto h = vec . back ( ) ; return vec . pop_back ( ) , h ; } ; cv_geometry geom ; const auto & mp = cell . provider ( ) ; const auto & m = mp . morphology ( ) ; if ( m . empty ( ) ) { geom . cell_cv_divs = { 0 , 0 } ; return geom ; } mlocation_list locs = thingify ( lset , mp ) ; auto neither_root_nor_terminal = [ & m ] ( mlocation x ) { return ! ( x . pos == 0 && x . branch == ( m . branch_children ( mnpos ) . size ( ) > 1u ? mnpos : 0 ) ) && ! ( x . pos == 1 && m . branch_children ( x . branch ) . empty ( ) ) ; } ; locs . erase ( std :: partition ( locs . begin ( ) , locs . end ( ) , neither_root_nor_terminal ) , locs . end ( ) ) ; util :: sort ( locs ) ; util :: unique_in_place ( locs ) ; constexpr fvm_index_type no_parent = -1 ; std :: vector < std :: pair < mlocation , fvm_index_type > > next_cv_head ; next_cv_head . emplace_back ( mlocation { mnpos , 0 } , no_parent ) ; mcable_list cables ; std :: vector < msize_t > branches ; geom . cv_cables_divs . push_back ( 0 ) ; fvm_index_type cv_index = 0 ; while ( ! next_cv_head . empty ( ) ) { auto next = pop ( next_cv_head ) ; mlocation h = next . first ; cables . clear ( ) ; branches . clear ( ) ; branches . push_back ( h . branch ) ; geom . cv_parent . push_back ( next . second ) ; while ( ! branches . empty ( ) ) { msize_t b = pop ( branches ) ; auto it = locs . end ( ) ; if ( b != mnpos && b == h . branch ) { it = std :: upper_bound ( locs . begin ( ) , locs . end ( ) , h ) ; } else if ( b != mnpos ) { it = std :: lower_bound ( locs . begin ( ) , locs . end ( ) , mlocation { b , 0 } ) ; } if ( it != locs . end ( ) && it -> branch == b ) { cables . push_back ( { b , b == h . branch ? h . pos : 0 , it -> pos } ) ; next_cv_head . emplace_back ( * it , cv_index ) ; } else { if ( b != mnpos ) { cables . push_back ( { b , b == h . branch ? h . pos : 0 , 1 } ) ; } for ( auto & c : m . branch_children ( b ) ) { branches . push_back ( c ) ; } } } sort ( cables ) ; append ( geom . cv_cables , std :: move ( cables ) ) ; geom . cv_cables_divs . push_back ( geom . cv_cables . size ( ) ) ; ++ cv_index ; } auto n_cv = cv_index ; arb_assert ( n_cv > 0 ) ; arb_assert ( geom . cv_parent . front ( ) == -1 ) ; arb_assert ( util :: all_of ( util :: subrange_view ( geom . cv_parent , 1 , n_cv ) , [ ] ( auto v ) { return v != no_parent ; } ) ) ; assign ( geom . cv_children , util :: make_span ( 1 , n_cv ) ) ; stable_sort_by ( geom . cv_children , [ & geom ] ( auto cv ) { return geom . cv_parent [ cv ] ; } ) ; geom . cv_children_divs . reserve ( n_cv + 1 ) ; geom . cv_children_divs . push_back ( 0 ) ; auto b = geom . cv_children . begin ( ) ; auto e = geom . cv_children . end ( ) ; auto from = b ; for ( fvm_index_type cv = 0 ; cv < n_cv ; ++ cv ) { from = std :: partition_point ( from , e , [ cv , & geom ] ( auto i ) { return geom . cv_parent [ i ] <= cv ; } ) ; geom . cv_children_divs . push_back ( from - b ) ; } geom . cv_to_cell . assign ( cv_index , 0 ) ; geom . cell_cv_divs = { 0 , ( fvm_index_type ) cv_index } ; geom . branch_cv_map . resize ( 1 ) ; std :: vector < pw_elements < fvm_size_type > > & bmap = geom . branch_cv_map . back ( ) ; for ( auto cv : util :: make_span ( n_cv ) ) { for ( auto cable : geom . cables ( cv ) ) { if ( cable . branch >= bmap . size ( ) ) { bmap . resize ( cable . branch + 1 ) ; } bmap [ cable . branch ] . push_back ( cable . prox_pos , cable . dist_pos , cv ) ; } } return geom ; }
Construct cv_geometry for cell from locset describing CV boundary points.
Construct cv_geometry for cell from locset describing CV boundary points .
// Filter out root, terminal locations and repeated locations so as to // avoid trivial CVs outside of fork points. (This is not necessary for // correctness, but is for the convenience of specification by lset.) // root? // terminal? // Collect cables constituting each CV, maintaining a stack of CV // proximal 'head' points, and recursing down branches in the morphology // within each CV. // head loc, parent cv index // Find most proximal point in locs on this branch, strictly more distal than h. // If found, use as an end point, and stop descent. // Otherwise, recurse over child branches. // Construct CV children mapping by sorting CV indices by parent. // Fill cv/cell mapping for single cell (index 0). // Build location query map. // Ordering of CV ensures CV cables on any given branch are found sequentially.