VoxelComplex.ih

/**
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU Lesser General Public License as
 *  published by the Free Software Foundation, either version 3 of the
 *  License, or  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 **/
 
/**
 * @file VoxelComplex.ih
 * @author Pablo Hernandez-Cerdan (\c pablo.hernandez.cerdan@outlook.com)
 * Insitute of Fundamental Sciences, Massey University, New Zealand.
 *
 * @date 2018/01/01
 *
 * Implementation of inline methods defined in VoxelComplex.h
 *
 * This file is part of the DGtal library.
 */
 
//////////////////////////////////////////////////////////////////////////////
#include <DGtal/graph/ObjectBoostGraphInterface.h>
#include <DGtal/topology/CubicalComplexFunctions.h>
#include <DGtal/topology/NeighborhoodConfigurations.h>
#include <DGtal/helpers/StdDefs.h>
#include <boost/graph/connected_components.hpp>
#include <boost/graph/filtered_graph.hpp>
#include <boost/property_map/property_map.hpp>
#include <iostream>
#ifdef WITH_OPENMP
// #include <experimental/algorithm>
#include <omp.h>
#endif
//////////////////////////////////////////////////////////////////////////////
// Default constructor:
template <typename TKSpace, typename TCellContainer>
inline DGtal::VoxelComplex<TKSpace, TCellContainer>::VoxelComplex()
    : Parent(),
      myTablePtr(nullptr), myPointToMaskPtr(nullptr),
      myIsTableLoaded(false) {}
 
// Copy constructor:
template <typename TKSpace, typename TCellContainer>
inline DGtal::VoxelComplex<TKSpace, TCellContainer>::VoxelComplex(
    const VoxelComplex &other)
    : Parent(other),
      myTablePtr(other.myTablePtr),
      myPointToMaskPtr(other.myPointToMaskPtr),
      myIsTableLoaded(other.myIsTableLoaded) {}
 
///////////////////////////////////////////////////////////////////////////////
// IMPLEMENTATION of inline methods.
///////////////////////////////////////////////////////////////////////////////
//-----------------------------------------------------------------------------
template <typename TKSpace, typename TCellContainer>
inline DGtal::VoxelComplex<TKSpace, TCellContainer> &
DGtal::VoxelComplex<TKSpace, TCellContainer>::
operator=(const Self &other)
{
    if (this != &other) {
        this->myKSpace = other.myKSpace;
        this->myCells = other.myCells;
        myTablePtr = other.myTablePtr;
        myPointToMaskPtr = other.myPointToMaskPtr;
        myIsTableLoaded = other.myIsTableLoaded;
    }
    return *this;
}
//---------------------------------------------------------------------------
///////////////////////////////////////////////////////////////////////////////
// Voxel methods :
///////////////////////////////////////////////////////////////////////////////
 
///////////////////////////////////////////////////////////////////////////////
// Interface - Voxel :
//---------------------------------------------------------------------------
// template <typename TKSpace, typename TCellContainer>
// const typename DGtal::VoxelComplex<TKSpace,
//                                    TCellContainer>::Object::Point &
// DGtal::VoxelComplex<TKSpace, TCellContainer>::objPointFromVoxel(
//     const Cell &voxel) const {
//     ASSERT(isSpel(voxel) == true);
//     ASSERT(this->belongs(voxel));
//     const auto &ks = this->space();
//     return *(myObject.pointSet().find(ks.uCoords(voxel)));
// }
 
 
//-----------------------------------------------------------------------------
 
template <typename TKSpace, typename TCellContainer>
template <typename TDigitalSet>
inline void DGtal::VoxelComplex<TKSpace, TCellContainer>::construct(
    const TDigitalSet &input_set,
    const Alias<ConfigMap> input_table)
{
    Parent::construct(input_set);
    setSimplicityTable(input_table);
}
 
template <typename TKSpace, typename TCellContainer>
void DGtal::VoxelComplex<TKSpace, TCellContainer>::setSimplicityTable(
    const Alias<ConfigMap> input_table)
{
    this->myTablePtr = input_table;
    this->myPointToMaskPtr =
        functions::mapZeroPointNeighborhoodToConfigurationMask<Point>();
    this->myIsTableLoaded = true;
}
 
template <typename TKSpace, typename TCellContainer>
void DGtal::VoxelComplex<TKSpace, TCellContainer>::copySimplicityTable(
    const Self & other)
{
    myTablePtr = other.myTablePtr;
    myPointToMaskPtr = other.myPointToMaskPtr;
    myIsTableLoaded = other.myIsTableLoaded;
}
 
template <typename TKSpace, typename TCellContainer>
const typename DGtal::VoxelComplex<TKSpace,
                                   TCellContainer>::ConfigMap &
DGtal::VoxelComplex<TKSpace, TCellContainer>::table() const
{
    return *myTablePtr;
}
 
template <typename TKSpace, typename TCellContainer>
const bool &
DGtal::VoxelComplex<TKSpace, TCellContainer>::isTableLoaded() const
{
    return myIsTableLoaded;
}
 
template <typename TKSpace, typename TCellContainer>
const typename DGtal::VoxelComplex<TKSpace,
                                   TCellContainer>::PointToMaskMap &
DGtal::VoxelComplex<TKSpace, TCellContainer>::pointToMask() const
{
    return *myPointToMaskPtr;
}
//---------------------------------------------------------------------------
template <typename TKSpace, typename TCellContainer>
inline void DGtal::VoxelComplex<TKSpace, TCellContainer>::voxelClose(
    const Cell &kcell)
{
    const auto &ks = this->space();
    ASSERT(ks.uDim(kcell) == 3);
    Dimension l = 2;
    auto direct_faces = ks.uLowerIncident(kcell);
    for (typename Cells::const_iterator cells_it = direct_faces.begin(),
                                        cells_it_end = direct_faces.end();
         cells_it != cells_it_end; ++cells_it) {
        this->insertCell(l, *cells_it);
    }
    cellsClose(l, direct_faces);
}
//---------------------------------------------------------------------------
template <typename TKSpace, typename TCellContainer>
inline void DGtal::VoxelComplex<TKSpace, TCellContainer>::cellsClose(
    Dimension k, const Cells &cells)
{
    if (k <= 0)
        return;
    if (cells.size() == 0)
        return;
    const auto &ks = this->space();
    Dimension l = k - 1;
    for (auto const &kcell : cells) {
        auto direct_faces = ks.uLowerIncident(kcell);
        for (typename Cells::const_iterator cells_it = direct_faces.begin(),
                                            cells_it_end = direct_faces.end();
             cells_it != cells_it_end; ++cells_it) {
            this->insertCell(l, *cells_it);
        }
        cellsClose(l, direct_faces);
    }
}
//---------------------------------------------------------------------------
template <typename TKSpace, typename TCellContainer>
inline void
DGtal::VoxelComplex<TKSpace, TCellContainer>::insertVoxelCell(
    const Cell &kcell, const bool &close_it, const Data &data)
{
    const auto &ks = this->space();
    ASSERT(ks.uDim(kcell) == 3);
    this->insertCell(3, kcell, data);
    if (close_it)
        voxelClose(kcell);
}
 
//---------------------------------------------------------------------------
template <typename TKSpace, typename TCellContainer>
inline void
DGtal::VoxelComplex<TKSpace, TCellContainer>::insertVoxelPoint(
    const Point &dig_point, const bool &close_it, const Data &data)
{
    const auto &ks = this->space();
    insertVoxelCell(ks.uSpel(dig_point), close_it, data);
}
 
//---------------------------------------------------------------------------
template <typename TKSpace, typename TCellContainer>
template <typename TDigitalSet>
inline void
DGtal::VoxelComplex<TKSpace, TCellContainer>::dumpVoxels(
    TDigitalSet & in_out_set) const
{
    const auto &ks = this->space();
    for (auto it = this->begin(3), itE = this->end(3) ; it != itE ; ++it ){
        in_out_set.insertNew(ks.uCoords(it->first));
    }
}
//-----------------------------------------------------------------------------
 
template <typename TKSpace, typename TCellContainer>
void DGtal::VoxelComplex<TKSpace, TCellContainer>::pointelsFromCell(
    std::set<Cell> &pointels_out, const Cell &input_cell) const
{
    const auto input_dim = this->space().uDim(input_cell);
    if (input_dim == 0) {
        pointels_out.emplace(input_cell);
        return;
    } else {
        auto ufaces = this->space().uFaces(input_cell);
        for (auto &&f : ufaces)
            this->pointelsFromCell(pointels_out, f);
    }
}
 
//---------------------------------------------------------------------------
template <typename TKSpace, typename TCellContainer>
void DGtal::VoxelComplex<TKSpace, TCellContainer>::spelsFromCell(
    std::set<Cell> &spels_out, const Cell &input_cell) const
{
    const auto input_dim = this->space().uDim(input_cell);
    if (input_dim == this->dimension) {
        if (this->belongs(input_cell))
            spels_out.emplace(input_cell);
        return;
    }
    auto co_faces = this->space().uCoFaces(input_cell);
    for (auto &&f : co_faces) {
        auto f_dim = this->space().uDim(f);
        if (f_dim >= input_dim)
            this->spelsFromCell(spels_out, f);
    }
}
 
//---------------------------------------------------------------------------
template <typename TKSpace, typename TCellContainer>
typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Clique
DGtal::VoxelComplex<TKSpace, TCellContainer>::Kneighborhood(
    const Cell &input_cell) const
{
    auto spels_out = neighborhoodVoxels(input_cell);
 
    const auto &ks = this->space();
    Clique clique(ks);
    for (const auto &v : spels_out)
        clique.insertCell(v);
 
    return clique;
}
 
template <typename TKSpace, typename TCellContainer>
std::set<typename TKSpace::Cell>
DGtal::VoxelComplex<TKSpace, TCellContainer>::neighborhoodVoxels(
    const Cell &input_spel) const
{
    std::set<Cell> pointels_out;
    std::set<Cell> spels_out;
    pointelsFromCell(pointels_out, input_spel);
    for (const auto &p : pointels_out)
        spelsFromCell(spels_out, p);
    return spels_out;
}
 
template <typename TKSpace, typename TCellContainer>
std::set<typename TKSpace::Cell>
DGtal::VoxelComplex<TKSpace, TCellContainer>::properNeighborhoodVoxels(
    const Cell &input_spel) const
{
 
    auto spels_out = neighborhoodVoxels(input_spel);
    auto search = spels_out.find(input_spel);
    if (search != spels_out.end()) {
        spels_out.erase(search);
    }
    return spels_out;
}
 
//---------------------------------------------------------------------------
template <typename TKSpace, typename TCellContainer>
bool DGtal::VoxelComplex<TKSpace, TCellContainer>::isSpel(
    const Cell &b) const
{
    return (this->space().uDim(b) == this->space().DIM);
}
 
//---------------------------------------------------------------------------
template <typename TKSpace, typename TCellContainer>
typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Cell
DGtal::VoxelComplex<TKSpace,
    TCellContainer>::surfelBetweenAdjacentSpels(const Cell &A, const Cell &B)
    const
{
    ASSERT(isSpel(A) == true);
    ASSERT(isSpel(B) == true);
    const auto &ks = this->space();
    // Digital coordinates
    auto &&orientation_BA = ks.uCoords(B) - ks.uCoords(A);
    ASSERT(orientation_BA.norm1() == 1);
    // Khalimsky Coordinates
    return ks.uCell(A.preCell().coordinates + orientation_BA);
}
//---------------------------------------------------------------------------
///////////////////////////////////////////////////////////////////////////////
// Cliques
template <typename TKSpace, typename TCellContainer>
std::pair<bool, typename DGtal::VoxelComplex<TKSpace,
    TCellContainer>::Clique>
DGtal::VoxelComplex<TKSpace, TCellContainer>::K_2(
        const typename KSpace::Point &A,
        const typename KSpace::Point &B,
        bool verbose) const
{
    const auto &ks = this->space();
    using KPreSpace = typename TKSpace::PreCellularGridSpace;
    auto orientation_vector = B - A;
    ASSERT(orientation_vector.norm1() == 1);
    int direction{-1};
    for (auto i = 0; i != ks.DIM; ++i)
        if (orientation_vector[i] == 1 || orientation_vector[i] == -1) {
            direction = i;
            break;
        }
 
    Point right, up;
    if (direction == 0) {
        right = {0, 0, 1};
        up = {0, 1, 0};
    } else if (direction == 1) {
        right = {1, 0, 0};
        up = {0, 0, 1};
    } else {
        right = {0, 1, 0};
        up = {1, 0, 0};
    }
 
    const PreCell x0 = KPreSpace::uSpel(A + right);
    const PreCell x4 = KPreSpace::uSpel(A - right);
    const PreCell x2 = KPreSpace::uSpel(A + up);
    const PreCell x6 = KPreSpace::uSpel(A - up);
 
    const PreCell x1 = KPreSpace::uSpel(A + up + right);
    const PreCell x3 = KPreSpace::uSpel(A + up - right);
    const PreCell x7 = KPreSpace::uSpel(A - up + right);
    const PreCell x5 = KPreSpace::uSpel(A - up - right);
 
    const PreCell y0 = KPreSpace::uSpel(B + right);
    const PreCell y4 = KPreSpace::uSpel(B - right);
    const PreCell y2 = KPreSpace::uSpel(B + up);
    const PreCell y6 = KPreSpace::uSpel(B - up);
 
    const PreCell y1 = KPreSpace::uSpel(B + up + right);
    const PreCell y3 = KPreSpace::uSpel(B + up - right);
    const PreCell y7 = KPreSpace::uSpel(B - up + right);
    const PreCell y5 = KPreSpace::uSpel(B - up - right);
 
    const auto bx0 = this->belongs(KSpace::DIM, x0);
    const auto bx1 = this->belongs(KSpace::DIM, x1);
    const auto bx2 = this->belongs(KSpace::DIM, x2);
    const auto bx3 = this->belongs(KSpace::DIM, x3);
    const auto bx4 = this->belongs(KSpace::DIM, x4);
    const auto bx5 = this->belongs(KSpace::DIM, x5);
    const auto bx6 = this->belongs(KSpace::DIM, x6);
    const auto bx7 = this->belongs(KSpace::DIM, x7);
 
    const auto by0 = this->belongs(KSpace::DIM, y0);
    const auto by1 = this->belongs(KSpace::DIM, y1);
    const auto by2 = this->belongs(KSpace::DIM, y2);
    const auto by3 = this->belongs(KSpace::DIM, y3);
    const auto by4 = this->belongs(KSpace::DIM, y4);
    const auto by5 = this->belongs(KSpace::DIM, y5);
    const auto by6 = this->belongs(KSpace::DIM, y6);
    const auto by7 = this->belongs(KSpace::DIM, y7);
 
    Clique k2_crit(ks);
    if (bx0)
        k2_crit.insertCell(ks.uCell(x0));
    if (bx1)
        k2_crit.insertCell(ks.uCell(x1));
    if (bx2)
        k2_crit.insertCell(ks.uCell(x2));
    if (bx3)
        k2_crit.insertCell(ks.uCell(x3));
    if (bx4)
        k2_crit.insertCell(ks.uCell(x4));
    if (bx5)
        k2_crit.insertCell(ks.uCell(x5));
    if (bx6)
        k2_crit.insertCell(ks.uCell(x6));
    if (bx7)
        k2_crit.insertCell(ks.uCell(x7));
 
    if (by0)
        k2_crit.insertCell(ks.uCell(y0));
    if (by1)
        k2_crit.insertCell(ks.uCell(y1));
    if (by2)
        k2_crit.insertCell(ks.uCell(y2));
    if (by3)
        k2_crit.insertCell(ks.uCell(y3));
    if (by4)
        k2_crit.insertCell(ks.uCell(y4));
    if (by5)
        k2_crit.insertCell(ks.uCell(y5));
    if (by6)
        k2_crit.insertCell(ks.uCell(y6));
    if (by7)
        k2_crit.insertCell(ks.uCell(y7));
    // Note that input spels A,B are ommited.
 
    /////////////////////////////////
    // Critical Clique Conditions:
    using namespace DGtal::functions;
    // Intersection of k2-neighborhood with the object:
    // (i) k2_clique must be empty or NOT 0-connected
    bool is_empty{k2_crit.nbCells(KSpace::DIM) == 0};
 
    // Check connectedness using object if not empty
    bool is_disconnected{false};
    if (!is_empty) {
        using DigitalTopology = DGtal::Z3i::DT26_6;
        using DigitalSet = DGtal::DigitalSetByAssociativeContainer<
            DGtal::Z3i::Domain,
            std::unordered_set<typename DGtal::Z3i::Domain::Point>>;
        using NewObject = DGtal::Object<DigitalTopology, DigitalSet>;
        auto new_obj = objectFromSpels<NewObject, KSpace, CellContainer>(k2_crit);
        auto con = new_obj->computeConnectedness();
        is_disconnected = (con == DISCONNECTED);
    }
 
    bool conditionI = is_empty || is_disconnected;
 
    // (ii) Xi or Yi belongs to this for i={0,2,4,6}
    std::vector<bool> bb(4);
    bb[0] = bx0 || by0;
    bb[1] = bx2 || by2;
    bb[2] = bx4 || by4;
    bb[3] = bx6 || by6;
 
    bool conditionII = bb[0] && bb[1] && bb[2] && bb[3];
    // is_critical if any condition is true.
    bool is_critical = conditionI || conditionII;
 
    if (verbose) {
        trace.beginBlock(" K2 critical conditions ");
        trace.info() << " conditionI = " << conditionI
                     << " : is_empty || is_disconnected =  " << is_empty
                     << " && " << is_disconnected << std::endl;
        trace.info() << " conditionII = " << conditionII << " : " << bb[0]
                     << " && " << bb[1] << " && " << bb[2] << " && " << bb[3]
                     << std::endl;
        trace.info() << " is_critical = " << is_critical
                     << " conditionI || conditionII : " << conditionI << " || "
                     << conditionII << std::endl;
        trace.endBlock();
    }
 
    // Return the clique (A,B), not the mask k2_crit
    Clique k2_clique(ks);
    Cell Ac = ks.uSpel(A);
    Cell Bc = ks.uSpel(B);
    k2_clique.insertCell(Ac);
    k2_clique.insertCell(Bc);
    return std::make_pair(is_critical, k2_clique);
}
 
//---------------------------------------------------------------------------
template <typename TKSpace, typename TCellContainer>
std::pair<bool, typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Clique>
DGtal::VoxelComplex<TKSpace, TCellContainer>::K_2(const Cell &A,
        const Cell &B,
        bool verbose) const
{
    // Precondition:
    // A and B are contiguous spels.
    ASSERT(isSpel(A) == true);
    ASSERT(isSpel(B) == true);
    const auto &ks = this->space();
    auto B_coord = ks.uCoords(B);
    auto A_coord = ks.uCoords(A);
    return this->K_2(A_coord, B_coord, verbose);
}
//---------------------------------------------------------------------------
 
template <typename TKSpace, typename TCellContainer>
std::pair<bool, typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Clique>
DGtal::VoxelComplex<TKSpace, TCellContainer>::K_2(const Cell &face2,
        bool verbose) const
{
    const auto &ks = this->space();
    ASSERT(ks.uIsSurfel(face2));
    using KPreSpace = typename TKSpace::PreCellularGridSpace;
    const auto co_faces = KPreSpace::uCoFaces(face2);
    const auto nco_faces = co_faces.size();
    ASSERT(nco_faces == 2);
    const auto &cf0 = co_faces[0];
    const auto &cf1 = co_faces[1];
    // spels must belong to complex.
    if (this->belongs(cf0) && this->belongs(cf1))
        return this->K_2(ks.uCell(cf0), ks.uCell(cf1), verbose);
    else
        return std::make_pair(false, Clique(ks));
}
//---------------------------------------------------------------------------
 
template <typename TKSpace, typename TCellContainer>
std::pair<bool, typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Clique>
DGtal::VoxelComplex<TKSpace, TCellContainer>::K_1(const Cell &face1,
        bool verbose) const
{
    const auto &ks = this->space();
    ASSERT(ks.uDim(face1) == 1);
    using KPreSpace = typename TKSpace::PreCellularGridSpace;
    // Get 2 orth dirs in orient_orth
    std::vector<Point> dirs_orth;
    for (auto q = KPreSpace::uOrthDirs(face1); q != 0; ++q) {
        const Dimension dir = *q;
        Point positive_orth{0, 0, 0};
        for (Dimension i = 0; i != ks.DIM; ++i)
            if (i == dir)
                positive_orth[i] = 1;
 
        dirs_orth.push_back(positive_orth);
    }
 
    auto &kface = face1.preCell().coordinates;
    const Point a{kface + dirs_orth[0] + dirs_orth[1]};
    const Point b{kface + dirs_orth[0] - dirs_orth[1]};
    const Point c{kface - dirs_orth[0] + dirs_orth[1]};
    const Point d{kface - dirs_orth[0] - dirs_orth[1]};
 
    const PreCell A = KPreSpace::uCell(a);
    const PreCell B = KPreSpace::uCell(b);
    const PreCell C = KPreSpace::uCell(c);
    const PreCell D = KPreSpace::uCell(d);
 
    // Now we need the other spels forming the mask
    // Get the direction (positive) linel spans.
    Point dir_parallel{0, 0, 0};
    for (auto q = KPreSpace::uDirs(face1); q != 0; ++q) {
        const Dimension dir = *q;
        for (Dimension i = 0; i != ks.DIM; ++i)
            if (i == dir)
                dir_parallel[i] = 1;
    }
    // Note that C, B are interchangeable. Same in A,D. Same between X and Y
    // sets Changed notation from paper: XA=X0, XB=X1, XC=X2, XD=X3
    // X
    const Point xa{a + 2 * dir_parallel};
    const Point xb{b + 2 * dir_parallel};
    const Point xc{c + 2 * dir_parallel};
    const Point xd{d + 2 * dir_parallel};
    // Y
    const Point ya{a - 2 * dir_parallel};
    const Point yb{b - 2 * dir_parallel};
    const Point yc{c - 2 * dir_parallel};
    const Point yd{d - 2 * dir_parallel};
 
    // Cell of the mask from KCoords
    const PreCell XA = KPreSpace::uCell(xa);
    const PreCell XB = KPreSpace::uCell(xb);
    const PreCell XC = KPreSpace::uCell(xc);
    const PreCell XD = KPreSpace::uCell(xd);
 
    const PreCell YA = KPreSpace::uCell(ya);
    const PreCell YB = KPreSpace::uCell(yb);
    const PreCell YC = KPreSpace::uCell(yc);
    const PreCell YD = KPreSpace::uCell(yd);
 
    /////////////////////////////////
    // Critical Clique Conditions:
 
    /** is_critical = ConditionI && ConditionII
     * (i) ConditionI:
     * At least one the sets {A,D},{B,C} is subset of this complex
     */
    /** (ii) ConditionII = B1 OR B2
     *  B1) (U & *this != empty) AND (V & *this != empty)
     *  B2) (U & *this == empty) AND (V & *this == empty)
     */
 
    const bool A1{this->belongs(KSpace::DIM, A) &&
                  this->belongs(KSpace::DIM, D)};
    const bool A2{this->belongs(KSpace::DIM, B) &&
                  this->belongs(KSpace::DIM, C)};
    const bool conditionI{A1 || A2};
 
    const bool u_not_empty{
        this->belongs(KSpace::DIM, XA) || this->belongs(KSpace::DIM, XB) ||
        this->belongs(KSpace::DIM, XC) || this->belongs(KSpace::DIM, XD)};
 
    const bool v_not_empty{
        this->belongs(KSpace::DIM, YA) || this->belongs(KSpace::DIM, YB) ||
        this->belongs(KSpace::DIM, YC) || this->belongs(KSpace::DIM, YD)};
 
    const bool B1{u_not_empty && v_not_empty};
    const bool B2{!u_not_empty && !v_not_empty};
    const bool conditionII{B1 || B2};
 
    const bool is_critical{conditionI && conditionII};
 
    if (verbose) {
        trace.beginBlock(" K1 critical conditions ");
        trace.info() << "input linel: " << face1 << std::endl;
        trace.info() << "is_critical = " << is_critical
                     << " conditionI || condition II " << conditionI << " || "
                     << conditionII << std::endl;
        trace.info() << " conditionI = " << conditionI << " = A1 || A2 : " << A1
                     << " || " << A2 << std::endl;
        trace.info() << " conditionII = " << conditionII
                     << " = B1 || B2 : " << B1 << " || " << B2 << std::endl;
        trace.endBlock();
    }
 
    // out clique is the intersection between mask and object
    Clique k1(ks);
    if (this->belongs(KSpace::DIM, A))
        k1.insert(ks.uCell(A));
    if (this->belongs(KSpace::DIM, B))
        k1.insert(ks.uCell(B));
    if (this->belongs(KSpace::DIM, C))
        k1.insert(ks.uCell(C));
    if (this->belongs(KSpace::DIM, D))
        k1.insert(ks.uCell(D));
    return std::make_pair(is_critical, k1);
}
//---------------------------------------------------------------------------
 
template <typename TKSpace, typename TCellContainer>
std::pair<bool, typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Clique>
DGtal::VoxelComplex<TKSpace, TCellContainer>::K_0(const Cell &face0,
        bool verbose) const
{
    const auto &ks = this->space();
    ASSERT(ks.uDim(face0) == 0);
    using KPreSpace = typename TKSpace::PreCellularGridSpace;
    auto &kface = face0.preCell().coordinates;
    const Point z{0, 0, 1};
    const Point y{0, 1, 0};
    const Point x{1, 0, 0};
 
    const Point a{kface + x - y - z};
    const Point b{kface - x - y - z};
    const Point c{kface + x - y + z};
    const Point d{kface - x - y + z};
 
    const Point e{kface + x + y - z};
    const Point f{kface - x + y - z};
    const Point g{kface + x + y + z};
    const Point h{kface - x + y + z};
 
    const PreCell A{KPreSpace::uCell(a)};
    const PreCell B{KPreSpace::uCell(b)};
    const PreCell C{KPreSpace::uCell(c)};
    const PreCell D{KPreSpace::uCell(d)};
 
    const PreCell E{KPreSpace::uCell(e)};
    const PreCell F{KPreSpace::uCell(f)};
    const PreCell G{KPreSpace::uCell(g)};
    const PreCell H{KPreSpace::uCell(h)};
 
    /////////////////////////////////
    // Critical Clique Conditions:
    /** is_critical = B1 || B2 || B3 || B4
     * where:
     * B1 = isSubset{A, H}
     * B2 = isSubset{B, G}
     * B3 = isSubset{C, F}
     * B4 = isSubset{D, E}
     * @note that the subsets define the 4 longest diagonals between the 8
     * pixels.
     */
    const bool bA = this->belongs(KSpace::DIM, A);
    const bool bB = this->belongs(KSpace::DIM, B);
    const bool bC = this->belongs(KSpace::DIM, C);
    const bool bD = this->belongs(KSpace::DIM, D);
    const bool bE = this->belongs(KSpace::DIM, E);
    const bool bF = this->belongs(KSpace::DIM, F);
    const bool bG = this->belongs(KSpace::DIM, G);
    const bool bH = this->belongs(KSpace::DIM, H);
 
    const bool B1{bA && bH};
    const bool B2{bB && bG};
    const bool B3{bC && bF};
    const bool B4{bD && bE};
    const bool is_critical{B1 || B2 || B3 || B4};
 
    if (verbose) {
        trace.beginBlock(" K0 critical conditions ");
        trace.info() << "input pointel: " << face0 << std::endl;
        trace.info() << "is_critical = B1 || B2 || B3 || B4 " << std::endl;
        trace.info() << is_critical << " = " << B1 << " || " << B2 << " || "
                     << B3 << " || " << B4 << std::endl;
        trace.endBlock();
    }
    // out clique is the intersection between mask and object
    Clique k0_out(ks);
    if (bA)
        k0_out.insert(ks.uCell(A));
    if (bB)
        k0_out.insert(ks.uCell(B));
    if (bC)
        k0_out.insert(ks.uCell(C));
    if (bD)
        k0_out.insert(ks.uCell(D));
    if (bE)
        k0_out.insert(ks.uCell(E));
    if (bF)
        k0_out.insert(ks.uCell(F));
    if (bG)
        k0_out.insert(ks.uCell(G));
    if (bH)
        k0_out.insert(ks.uCell(H));
 
    return std::make_pair(is_critical, k0_out);
}
//---------------------------------------------------------------------------
 
template <typename TKSpace, typename TCellContainer>
std::pair<bool, typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Clique>
DGtal::VoxelComplex<TKSpace, TCellContainer>::K_3(const Cell &voxel,
        bool verbose) const
{
    const auto &ks = this->space();
    ASSERT(ks.uDim(voxel) == 3);
    const bool is_critical = !isSimple(voxel);
 
    if (verbose) {
        trace.beginBlock(" K3 critical conditions ");
        trace.info() << "input voxel: " << voxel << std::endl;
        trace.info() << "is_critical = " << is_critical << std::endl;
        trace.endBlock();
    }
 
    Clique clique(ks);
    clique.insertCell(voxel);
    return std::make_pair(is_critical, clique);
}
//---------------------------------------------------------------------------
 
/* BUG workaround: MSVC compiler error C2244.
 * It doesn't see the definition of these declarations (Moved to header)
template <typename TKSpace, typename TCellContainer>
std::array<
  typename DGtal::VoxelComplex<TKSpace,
TCellContainer>::CliqueContainer, DGtal::VoxelComplex<TKSpace,
TCellContainer>::dimension + 1
>
DGtal::VoxelComplex<TKSpace, TCellContainer>::criticalCliques(
    const Parent & cubical,
    bool verbose
    ) const
 
template <typename TKSpace, typename TCellContainer>
std::array<
  typename DGtal::VoxelComplex<TKSpace,
TCellContainer>::CliqueContainer, DGtal::VoxelComplex<TKSpace,
TCellContainer>::dimension + 1
>
DGtal::VoxelComplex<TKSpace, TCellContainer>::criticalCliques(
    bool verbose
    ) const
*/
//---------------------------------------------------------------------------
 
template <typename TKSpace, typename TCellContainer>
std::pair<bool, typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Clique>
DGtal::VoxelComplex<TKSpace, TCellContainer>::criticalCliquePair(
        const Dimension d, const CellMapConstIterator &cellMapIterator) const
{
    const auto &it = cellMapIterator;
    const auto &cell = it->first;
    // auto &cell_data = it->second;
    auto clique_p = std::make_pair(false, Clique(this->space()));
    if (d == 0)
        clique_p = K_0(cell);
    else if (d == 1)
        clique_p = K_1(cell);
    else if (d == 2)
        clique_p = K_2(cell);
    else if (d == 3)
        clique_p = K_3(cell);
    else
        throw std::runtime_error("Wrong dimension: " + std::to_string(d));
 
    return clique_p;
}
//---------------------------------------------------------------------------
 
template <typename TKSpace, typename TCellContainer>
typename DGtal::VoxelComplex<TKSpace, TCellContainer>::CliqueContainer
DGtal::VoxelComplex<TKSpace, TCellContainer>::criticalCliquesForD(
    const Dimension d, const Parent &cubical, bool verbose) const
{
#ifdef WITH_OPENMP
 
    ASSERT(dimension >= 0 && dimension <= 3);
    CliqueContainer critical;
 
    const auto nthreads = omp_get_num_procs();
    omp_set_num_threads(nthreads);
    std::vector<CliqueContainer> p_critical;
    p_critical.resize(nthreads);
#pragma omp parallel
    {
#pragma omp single nowait
        {for (auto it = cubical.begin(d), itE = cubical.end(d); it != itE; ++it)
#pragma omp task firstprivate(it)
             {auto clique_p = criticalCliquePair(d, it);
    const auto &is_critical = clique_p.first;
    const auto &clique = clique_p.second;
    // Push
    auto th = omp_get_thread_num();
    if (is_critical)
        p_critical[th].push_back(clique);
}
} // cell loop
#pragma omp taskwait
}
// Merge
std::size_t total_size = 0;
for (const auto &sub : p_critical)
    total_size += sub.size();
 
critical.reserve(total_size);
for (const auto &sub : p_critical)
    critical.insert(critical.end(), sub.begin(), sub.end());
 
if (verbose)
    trace.info() << " d:" << d << " ncrit: " << critical.size();
return critical;
 
#else
 
    ASSERT(dimension >= 0 && dimension <= 3);
    CliqueContainer critical;
    for (auto it = cubical.begin(d), itE = cubical.end(d); it != itE; ++it) {
        const auto clique_p = criticalCliquePair(d, it);
        auto &is_critical = clique_p.first;
        auto &clique = clique_p.second;
        if (is_critical)
            critical.push_back(clique);
    } // cell loop
    if (verbose)
        trace.info() << " d:" << d << " ncrit: " << critical.size();
    return critical;
 
#endif
}
//---------------------------------------------------------------------------
///////////////////////////////////////////////////////////////////////////////
template <typename TKSpace, typename TCellContainer>
bool DGtal::VoxelComplex<TKSpace, TCellContainer>::isSimpleByThinning(
        const Cell &input_spel) const
{
    // x = input_spel ; X = VoxelComplex ~ occupancy of khalimsky space
    // a) Get the neighborhood (voxels) of input_spel intersected
    // with the voxel complex. -- N^{*}(x) intersection X --
    ASSERT(this->space().uDim(input_spel) == 3);
    const auto spels_out = this->properNeighborhoodVoxels(input_spel);
    const auto &ks = this->space();
    Clique clique(ks);
    for (const auto &v : spels_out)
        clique.insertCell(v);
    clique.close();
    // b) Apply a thinning on the result of a)
    typename Parent::DefaultCellMapIteratorPriority default_priority;
    bool clique_is_closed = true;
    functions::collapse( clique, spels_out.begin(), spels_out.end(), default_priority, false /* spels_out is not closed */, clique_is_closed, false /*verbose*/);
    // c) If the result is a single pointel, it is reducible
    return clique.size() == 1;
}
 
// Object wrappers :
template <typename TKSpace, typename TCellContainer>
bool DGtal::VoxelComplex<TKSpace, TCellContainer>::isSimple(
    const Cell &input_cell) const
{
    ASSERT(isSpel(input_cell) == true);
 
    if (myIsTableLoaded) {
        auto conf = functions::getSpelNeighborhoodConfigurationOccupancy<Self>(
            *this, this->space().uCoords(input_cell), this->pointToMask());
        return (*myTablePtr)[conf];
    } else
        return isSimpleByThinning(input_cell);
}
//---------------------------------------------------------------------------
///////////////////////////////////////////////////////////////////////////////
// Interface - public :
 
/**
 * Writes/Displays the object on an output stream.
 * @param out the output stream where the object is written.
 */
template <typename TKSpace, typename TCellContainer>
inline void DGtal::VoxelComplex<TKSpace, TCellContainer>::selfDisplay(
    std::ostream &out) const
{
    out << "[VoxelComplex dim=" << this->dim() << " chi=" << this->euler();
    out << " isTableLoaded? " << ((isTableLoaded()) ? "True" : "False");
}
//---------------------------------------------------------------------------
 
/**
 * Checks the validity/consistency of the object.
 * @return 'true' if the object is valid, 'false' otherwise.
 */
template <typename TKSpace, typename TCellContainer>
inline bool
DGtal::VoxelComplex<TKSpace, TCellContainer>::isValid() const
{
    return true;
}
 
//-----------------------------------------------------------------------------
template <typename TKSpace, typename TCellContainer>
inline std::string
DGtal::VoxelComplex<TKSpace, TCellContainer>::className() const
{
    return "VoxelComplex";
}
 
///////////////////////////////////////////////////////////////////////////////
// Implementation of inline functions                                        //
 
template <typename TKSpace, typename TCellContainer>
inline std::ostream &DGtal::
operator<<(std::ostream &out,
           const VoxelComplex<TKSpace, TCellContainer> &object)
{
    object.selfDisplay(out);
    return out;
}
 
//                                                                           //
///////////////////////////////////////////////////////////////////////////////