\section{Introduction}\label{introduction}

Given a collection \(\mathcal{S}\) of geometric
objects\footnote{Examples include, but are not limited to: line segments, quads, triangles, polygons, meshes, pixels, voxels, volume images, B-reps, \emph{etc.}
In mathematical terms, a geometric object is a topological space embedded in some $\mathbb{E}^d$ \cite{Edelsbrunner:95}.},
the subject of this paper is computing the topology of their space
arrangement \(\mathcal{A}(\mathcal{S})\) as a \emph{chain complex},
\emph{i.e.},~as a short exact sequence of linear spaces \(C_i\) of (co)chains and linear
boundary/coboundary maps \(\partial_p\) and
\(\delta_p=\partial_{p+1}^\top\) between them:
\[
C_\bullet = (C_p, \partial_p) :=
C_3 \
\substack{
\delta_2 \\
\longleftarrow \\[-1mm]
\longrightarrow \\
\partial_3
}
\ C_2 \
\substack{
\delta_1 \\
\longleftarrow \\[-1mm]
\longrightarrow \\
\partial_2
}
\ C_1 \
\substack{
\delta_0 \\
\longleftarrow \\[-1mm]
\longrightarrow \\
\partial_1
}
\ C_0 .
\]

The $C_\bullet$ chain complex fully characterizes the topology of the
cellular arrangement induced within the Euclidean space by a collection
of geometric objects embedded in it. The data structures needed for such
computational program are sparse
multi-arrays, and their standard algebraic operations.
In this way, we introduce a novel approach based on piecewise-linear algebraic topology~\cite{hatcher:2002,Rourke:Sanderson:1982}, that allows us to treat rather  general cellular complexes, with cells homeomorphic to polyhedra, \emph{i.e.},~to triangulable spaces~\cite{Rourke:Sanderson:1982}, and hence possibly non convex and multiply connected.

We believe that basic geometric algebraic topology
provides the right set of mathematical concepts and tools to compute
and explore the cells of the space partition induced by a set of
geometric objects, as well as the related incidence/neighborhood relations.
The current escalation of quantity/quality of data, and their drift
through pipelines of micro/macro-services that need simple interfaces, along with the fast diffusion of hybrid architectures for advanced applications, contribute to
motivate this paper.
The notions we deal with include
geometric complexes, linear spaces of chains and cochains, the \emph{chain
complex} of linear maps between pairs of spaces, and their
compositions. The discussion is restricted to piecewise-linear topology
and to space dimensions less or equal to three. The paper formalizes the
algorithms to generate the matrices of (co)boundary operators on a
cellular arrangement in $\mathbb{E}^d$ space, with $d\in\{2,3\}$.

To our knowledge, algorithms for chain complex computing  only include  the recent paper~\cite{Alayrangues2015}, where combinatorial generalized maps~\cite{Lienhardt:2014}, a quite intricate data structure, are used for cycles/boundary calculation of the homology of a cellular complex.  Is is conversely well known~\cite{Munkres:84,Kannan:79,Paoluzzi:1993:DMS:169728.169719,Edelsbrunner:95} that for simplicial  complexes, \emph{i.e.},~triangulations, boundary operators are defined as linear extensions of basic boundary operators which act on simplices.

Construction of arrangements of lines, segments, planes and other geometrical objects is discussed in~\cite{fhktww-a-07}, with a description of CGAL software~\cite{Fabri:2000:DCC:358668.358687}, implementing 2D/3D arrangements with Nef polyhedra~\cite{Hachenberger:2007:BOS:1247750.1248141,bieri:95}. A review of papers and algorithms concerning construction and counting of cells may be found in the chapter on Arrangements in the `Handbook of Discrete and Computational Geometry'~\cite{Goodman:2017:HDC:285869}.
Arrangements of polytopes, hyperplanes and $d$-circles are discussed in~\cite{Ziegler:92}.

The standard way to look at combinatorial data structures is the IG
(Incidence Graph) data structure  described in the book
`Algorithms in Combinatorial Geometry' \cite{Edelsbrunner:1987:ACG:28905}.
IG is an implementation of the Hasse diagram~\cite{Birkhoff:1948} of the cells of a
\(d\)-complex \cite{ieee-tase}. Our (co)chain complex provides an
algebraic representation of the Hasse diagram with sparse matrices,
associating each two adjacent levels with a (co)boundary
map to traverse up and down the hierarchy.

Some early papers are concerned with the efficient representation of 3D
cellular decompositions~\cite{Baumgart:1972:WEP:891970,Muller:78}. In particular,
\cite{Dobkin:1987:PMT:41958.41967} defines the polygon-edge data
structure to represent orientable and contractible 3D decompositions
and their duals. Several other systems have been developed about three
decades ago, to handle the merging of complexes in the context of solid
modeling and manufacturing automation.  The merging of intersecting shells
of polyhedral solids is computable by finding all intersecting boxes of faces in $O(n\, log^2(n)+ k)$, where $n$ is the number of boxes and $k$ is the number of box-pair intersections~\cite{Hoffmann:1987:RSO:866286}.
The Selective Geometric Complex (SGC)
was introduced~\cite{Rossignac:89}  in 1989, allowing for lower dimensional cells contained inside the interior of cells. Later, scientists in UK proposed the Djin API
\cite{Woodwark:99,Middleditch:92}, with the theoretical
foundations for merging operation of cellular complexes, but (non-manifold) boundary representations remained the standard in the field.

Recently, a systematic procedure has been proposed in~\cite{Zhou:2016:MAS:2897824.2925901} for constructing a family of exact constructive solid geometry
operations, starting from a collection of boundary triangular meshes.

Most of earlier algorithms and procedures~\cite{4055948,
Ala:1992:PAB:616022.617736,
Baumgart:1972:WEP:891970,
bowyer1995introducing,
bowyer1995svlis,
Braid:1975:SSB:360715.360727,
Brisson:1989:RGS:73833.73858,
cadanda,
Dobkin:1987:PMT:41958.41967,
Gomes:1999:MMB:304012.304039,
Guibas:1985:PMG:282918.282923,
HoffmannK01,
ieee-tase,
Kalay:1989:HET:63718.63719,
Lee:2001:PES:376957.376976,
Lienhardt:1991:TMB:115604.115610,
Mantyla:1988:ISM:60949,
Paoluzzi:1989:BAO:70248.70249,
Paoluzzi:1993:DMS:169728.169719,
Paoluzzi:1995:GPP:212332.212349,
Pascucci:1995:DCB:218013.218055,
Pratt94ashape,
Raghothama:1999:CUD:304012.304019,
Rap97,
Requicha:1980:RRS:356827.356833,
RequichaVoelcker:77,
Rossignac:1991:CNG:115604.115606,
Rossignac:SGC:90,
Shapiro:1991:RSS:124951,
Shapiro:1995:PFS:218013.218029,
Silva:81,
Weiler:86,
Weiler:88,
Woo:85,
wozny1990geometric,
Yamaguchi:85,
yamaguchi1995ntb,
Zhou:2016:MAS:2897824.2925901,
bieri:95,
Rossignac:89,
Hoffmann:91,
Hoffmann:1989:GSM:74803,
Hoffmann:1987:RSO:866286}
work with data
structures optimized for selected classes of geometric
objects. Contrary-wise, our formulation, cast in
terms of (co)chain complexes of (co)boundary maps, may be
applied to very different geometric objects, ranging from solid models
to engineering meshes, geographical systems, biomedical images.

Numerical methods aiming to integrate domain modeling, differential
topology and mathematical modeling with physical simulations are also based
on chains and cochains~\cite{PALMER1995733,Palmer1993}. In particular, Discrete Exterior Calculus
(DEC) with simplicial complexes was introduced by
\cite{Hirani:2003:DEC:959640} and made popular by
\cite{Desbrun:2006:DDF:1185657.1185665,Elcott:2006:BYO:1185657.1185666}.
FEEC is a recent advance~\cite{arnold_falk_winther_2006,Arnold:2010,Arnold:2018} in the mathematics of finite element methods that employs differential complexes to construct stable numerical schemes.
The Cell Method~(CM)
is a purely algebraic computational method
for modeling and simulation~\cite{Tonti:1975,Tonti:2013,Ferretti:2014} based on boundary/coboundary maps and
 a direct discrete formulation of field laws.
Our own research in geometrical and physical modeling  with chain and cochain complexes was introduced in \cite{DiCarlo:2009:DPU:1629255.1629273,ieee-tase,Dicarlo:2014:TNL:2543138.2543294}.

All merge algorithms, by nature, follow the same steps. The
advantages of formulating them in terms of (co)chain complexes and operations on sparse matrices are that (1)  the common and general algebraic topological
nature of this operation is revealed; (2) implementation specific low-level
details and algorithms are hidden; (3) explicit connection to SpMV kernels (sparse matrix-vector  multiplication) and to sparse
numerical linear algebra systems~\cite{Bell:SpMV:NVIDIA:2008,gemmexp,Davis:2006:DMS:1196434,BEKS14}
on GPU and HPC platforms is provided; (4) systematic and rigorous development of the
algorithms, that are correct by construction, is supported.

\begin{figure}[h]
   \centering
   \includegraphics[height=0.43\linewidth,width=0.5\linewidth]{figs/arxiv1.pdf}%
   \includegraphics[height=0.43\linewidth,width=0.5\linewidth]{figs/arxiv2.pdf}%

   {\footnotesize\hspace{.23\textwidth}(a)\hfill(b)\hspace{.23\textwidth}}
   \vspace{-2mm}
   \caption{ Examples:
   (a) Regularized 2D arrangement $X_2$ of the plane generated by a set of random line segments. The Euler characteristic is $\chi = \chi_0 - \chi_1 + \chi_2 = 11361 - 20813 + 9454 = 2$; (b) merge of two 3-complexes with $2\times 10^3$ 3-cells. The Euler characteristic (of the non-exploded resulting 3-complex) is $\chi = \chi_0 - \chi_1 + \chi_2 - \chi_3 = 8787 - 26732 + 26600 - 8655 = 0$. This  count includes the outer (unbounded) 2-cell or 3-cell, respectively, that are also computed by theTopological Gift Wrapping  algorithm TGW (see Section~\ref{tgw-algorithm}).
   The Euler characteristic of the $d$-sphere is $\chi = 1 + (-1)^d =$ 2 or 0 for either even or odd space dimension $d$. }
   \label{fig:arxiv}
\end{figure}

The paper is organized as follows. A brief introduction to the proposed
computational pipeline is given in Section \ref{sec:pipeline}, going
from the arrangement of \(\mathbb{E}^2\) induced by a set of line segments, to
the arrangement of \(\mathbb{E}^3\) induced by a collection of open/closed
piecewise-linear surfaces and/or meshes.
In the short subsection~\ref{sec:prologue}, we show that representing chains as sparse arrays is compact and flexible.
In Section \ref{sec:algorithms},
the algorithms for computing (co)boundary operators are presented in pseudocode format. In Section \ref{other-works}, our approach is compared with the relevant literature, remarking differences and merits. Past development and current prospects of this project are outlined in Section \ref{sec:projects}. The closing section presents a summary of contents, and outlines
 possible applications of ideas.
The Appendix  gives the theoretical minimum on chains and examples of topology computation.

\section{Computational pipeline}
\label{sec:pipeline}

Let us start with an input collection $\mathcal{S}$ of piecewise-linear geometric objects of  $(d-1)$ dimension, embedded in $\mathbb{E}^d$ space, with $d \in \{2,3\}$.  Examples include soups of lines or  polygons, triangled surfaces,  quads from cubical meshes, 1-, 2-, or 3-cells from 2D or 3D image elements (pixels or voxels, respectively), 2-skeletons/boundaries of triangulated polyhedra, non manifold B-reps or decompositive reps of solid models. These objects are \emph{geometric complexes}, \emph{i.e.}, pairs $(X,\mu)$, where $X$ is a cellular complex\footnote{See Appendix~\ref{chain-complexes-1-page-definitions} for this and related definition(s).} specifying the topology and $\mu: X_0 \to \mathbb{E}^d$ is the embedding function of 0-cells, sufficient for a piecewise-linear geometry.
The data may contains both $(d-1)$- and $d$-complexes: the combinatorial union of their $(d-1)$-skeletons is selected as the actual input to the pipeline.
An admissible input collection $\mathcal{S}$ of geometric complexes will mutually intersect and partition $\mathbb{E}^d$ into a cellular complex $X = \bigcup X_p$ $(0\leq p\leq d)$, called the \emph{arrangement} $\mathcal{A}(\mathcal{S})$ induced by $\mathcal{S}$.

The object of this paper is the computation of the chain complex $C_\bullet(X) = (C_p,\partial_p)$, starting from some representation\footnote{Our prototype implementation in \\ \href{https://github.com/cvdlab/LinearAlgebraicRepresentation.jl/tree/julia-1.0}{https://github.com/cvdlab/LinearAlgebraicRepresentation.jl/tree/julia-1.0}, makes use of LAR representation~\cite{Dicarlo:2014:TNL:2543138.2543294}, on which this approach strongly relies.} of $\mathcal{S}$. In particular, we compute the matrices of the linear maps $\partial_p$ (and their duals $\delta_{p-1}$)  between chain spaces $C_p$. Definitions and examples are given in Appendix~\ref{sec:appendix}.
Since the matrix of a linear map $C_p \to C_{p-1}$ between linear spaces contains by columns the target space representation of the domain space basis elements, the paper also provides constructive algorithms to generate a sparse matrix representation of basis elements $u_p\in C_p$, which are one-to-one with $p$-cells $\sigma_p$ in $X_p$ skeletons. Note that cells in $X$, specifically in $X_d$, are not known in advance.

\begin{figure}[htbp]
   \centering
   \includegraphics[width=\linewidth]{figs/cartoon.png}%

   \caption{Cartoon display of the computational pipeline: (a) two solids in $\mathcal{S}$; (b) the exploded input collection $\mathcal{S}_2$ in $\mathbb{E}^3$; (c) 2-cell $\sigma$ (red) and the set $\Sigma(\sigma)$ (blue) of possible intersection; (d)  $\sigma\cup\Sigma$ affinely mapped on $z=0$;
   (e)~reduction to a set of 1D segments in $\mathbb{E}^2$ via intersection with $z=0$; (f) pairwise intersections; (g) exploded $U_2$ basis of $C_2$ generated as columns of $\partial_2: C_2 \to C_1$; (h) exploded $U_3$ basis of $C_3$ generated as columns of operator's $\partial_3: C_3 \to C_2$ sparse matrix, via the TGW algorithm in 3D.    }
   \label{fig:process}
\end{figure}

The computation is correct and robust  because the boundaries of adjacent decomposed 2-cells are compatible as cellular complexes by construction. This fact is induced by continuity of topological spaces, mathematically modeled here by chains of cells of a complex. A requirement of the standard definition of a cellular complex~\cite{Munkres:84,hatcher:2002} demands boundary compatibility to hold. This fact is guaranteed here, since abutting subsets of 1-cells have non-empty intersection, so they generate congruent 0-, 1-cells.
A summary of the computational steps for $d=3$ follows.

\subsection{2D arrangements generated by 2-cells (Merge)}\label{facet-arrangements}

Let $\mathcal{S}_2 \subseteq \mathcal{S}$ be the set of 2-cells of input geometric complexes\footnote{Cell complex embedded in Euclidean space via association of position vectors to 0-cells.}, embedded in $\mathbb{E}^3$.
Note that $\mathcal{S}_2$ is not required to be a cellular complex, since cells may intersect outside of their boundaries.
It is only required that each cell is connected and manifold. Each $\sigma\in S_2$ is mapped to subspace $z=0$ by an affine transformation $\T{Q}_\sigma$, together with the set $\mathcal{I}(\sigma)\subset S_2$ of cells potentially intersecting it. The set $\Sigma = X_1(\sigma \cup \mathcal{I}(\sigma))$ is intersected with $x=0$, producing a set $\mathcal{S}_1(\sigma)$ of line segments in $\mathbb{E}^2$.

First, these are mutually intersected, producing the chain complex $C_\bullet(\sigma) = (C_{2,1,0},\,\partial_{2,1})$ generated by $\mathcal{A}(\mathcal{S}_1)$ (see Figure~\ref{fig:2D}d).
Care must be taken to identify the 1-cycles around holes within partitioned 2-cells, in order to remove their outer boundary cycles, by removing their columns, as well the outer cycle, from the operator matrix. Identification is easy: each hole produces two opposite columns summing to 0.
Finally, the geometric 2-complex $X_\sigma$ is transformed back in $\mathbb{E}^3$ by $\T{Q}^{-1}_\sigma$.
In summary, Algorithm~\ref{alg:two} executes the one-to-one map $\sigma\mapsto X_\sigma$, by computing the maps $\sigma\mapsto C_\bullet(\sigma)$ independently from each other. The output is a  set $\varmathbb{C} := \{ C_\bullet(\sigma), \sigma\in \mathcal{S}_2\}$.

\begin{figure}[htbp]
   \centering
   \includegraphics[width=0.65\linewidth]{figs/ch2-planararrangement.png}
   \caption{Computation of the regularized arrangement of a set of lines in $\mathbb{E}^2$: (a) the input, \emph{i.e.},~the 2-cell $\sigma$ (signed blue) and the line segment intersections of $\Sigma(\sigma)$ with $z=0$; (b) all pairwise intersections of 1-cells; (c)~removal of the 1-subcomplex external to $\partial\sigma$; (d) the 2-chain generated by $\sigma\cup\Sigma$ via TGW in 2D.}
   \label{fig:2D}
\end{figure}

\subsection{Quotient set computations (Congruence)}\label{quotient-set}

The idea allowing us to compute independent fragmentations of 2-cells comes from a similitude between homology and congruence. Two $(d-1)$-spaces (curves, surfaces, \emph{etc.}) embedded in $\mathbb{E}^d$  are topologically \emph{homologous} when their boundaries can be glued, enclosing a portion of the ambient space, and subdivide $\mathbb{E}^d$ in two parts, inner and outer.
Two geometric figures are geometrically \emph{congruent} iff one can be transformed into the other by an isometry~\cite{Coxeter:1967}. The congruences $R_p$ between $p$-cells of geometric complexes in $\varmathbb{C}$ are equivalence relations, so we may compute the chain complex of  quotient chain spaces,
\[
C_2(U_2/R_2) \overset{\partial_2}{\longrightarrow} C_1(U_1/R_1) \overset{\partial_1}{\longrightarrow} C_0(U_0/R_0),
\]
over which subsequently build the yet unknown basis of $C_3$. Note that  $U_p = \bigcup\,U_p^\sigma$, with $\sigma\in\mathcal{S}_2$, is the union of bases of fragmented $p$-chains in $\varmathbb{C}$, module the congruence relations $R_p$. \
$C_p(U_p/R_p)$ stands for the chain space generated by $X_p = U_p/R_p$.
 In this stage of the computational pipeline, we compute, for each $\sigma\in\mathcal{S}_2$, the quotient sets and the maps $\partial_p$ in-between, for $p=0,1,2$.

As usual, we proceed by  an inductive procedure: each stage consists in glueing cells of given dimension to the result of the previous stage~\cite{Baladze:EoM}. The construction of the sparse matrix of the signed operator $\partial_1: C_1\to C_0$ is straightforward: for each $u_1^h = u_0^{k_2} - u_0^{k_1}$, just write $\partial_1[k_2,h]=1$ and $\partial_1[k_1,h]=-1$, by convention for $k_2>k_1$. For details of quotient operations, see Section~\ref{quotients}.

\begin{example}[Merge with incompatible boundaries]
\label{ex:square}

Here we discuss the merge result when the input collection is defined by a pair of adjacent complexes with incompatible sub-complexes along their interface. Such a 2D example is displayed in Figure~\ref{fig:abutting}.
Another simple example may include two tetrahedralized unit cubes that are incident on a planar face triangulated into two triangles but along different diagonals. Similarly to the 2D case, each boundary triangle would be fragmented against all the incident ones, producing  fragmented output faces, so that the result on the common affine support (say, the vertical plane) would be exactly four triangles, because each input triangle is fragmented by the other diagonal. A larger  example is shown in Figure~\ref{fig:arxiv}a, where the 2D arrangement generated by a number of random line segments is displayed.
\begin{figure}[htbp]
   \centering
   \includegraphics[width=\linewidth]{figs/square-3}
   \\
   {\footnotesize \hspace{0.06\linewidth}(a)\hspace{0.38\linewidth}(b)\hspace{0.38\linewidth}(c)\hspace{0.04\linewidth}}\vspace{-3mm}
   \caption{Merge of two 2D complexes with incompatible boundaries: (a) Input data $\mathcal{S}_2 = \Sigma^1 \cup \Sigma^2$; (b)~independent fragmentation of 2-cells $\sigma\in\mathcal{S}_2$ induced by $\mathcal{I}(\sigma)$; (c)  \emph{local} arrangement $X_2 = \mathcal{A}(\mathcal{S}_1(\mathcal{S}_2))$ generated by \emph{Merge},  2D \emph{TGW}, and \emph{Congruence} algorithm pipeline.
   }
   \label{fig:abutting}
\end{figure}
\end{example}
\vspace{-3mm}
Note in Figure~\ref{fig:abutting} that: (i) each 2-cell $\sigma$ is processed \emph{independently} as $\Sigma = \{\sigma\}\cup\mathcal{I}(\sigma)$; (ii) the bigger diagonal is fragmented by the normal 1-cell; (iii) the unit 2-chains are reconstructed in 2D (before reduction by congruence) as 1-cycles by the TGW algorithm; (iv) 0-cells and 1-cells from fragmented diagonal are finally identified \emph{mod congruence}.

\subsection{Topological gift wrapping (TGW)}\label{tgw-algorithm}

The  algorithm discussed here is used to compute topologically in 2D/3D, respectively, the sparse matrices of signed operators $\partial_2$ and $\partial_3$, starting from $\partial_1$ and $\partial_2$ input, respectively.  The matrix $[\partial_d]$ of the linear map $C_d \to C_{d-1}$ between linear spaces contains by columns the target space representation of domain space basis elements, as closed $(d-1)$-chains, \emph{i.e.},~$(d-1)$-cycles.
Within the computational pipeline discussed in this paper, TGW is used \emph{locally} for each 2-cell to be decomposed, and \emph{globally} to generate the 3-cells of the arrangement of the ambient  space $\mathbb{E}^3$. The algorithm pseudocode is given and discussed in Section~\ref{partial_3-computation}. The needed extensions for handling holes and non-connected components are detailed in Section~\ref{delta_2-computation}.

\begin{figure*}[htbp]
   \centering
\includegraphics[width=0.2\textwidth]{figs/1.png}\includegraphics[width=0.2\textwidth]{figs/2.png}\includegraphics[width=0.2\textwidth]{figs/3.png}\includegraphics[width=0.2\textwidth]{figs/4.png}\includegraphics[width=0.2\textwidth]{figs/5.png}

{\footnotesize\hspace{.09\textwidth}(a)\hfill(b)\hfill(c)\hfill(d)\hfill(e)\hspace{.09\textwidth}}
\vspace{-2mm}
\caption{Extraction of a minimal 1-cycle from $\mathcal{A}(X_1)$: (a) the initial value for $c\in C_1$ and the signs of its oriented boundary; (b) cyclic subgroups on $\delta\partial c$; (c) new (coherently oriented) value of $c$ and $\partial c$; (d) cyclic subgroups on $\delta\partial c$; (e) final value of $c$, with $\partial c = \emptyset$. The step-by-step computation is discussed in Example~\ref{example1}.}
\label{fig:step-by-step}
\end{figure*}

\begin{figure*}[htbp]
\hfill\includegraphics[height=0.15\textwidth,width=0.1428\textwidth]{figs/3D1a.png}%
\includegraphics[height=0.15\textwidth,width=0.1428\textwidth]{figs/3D2.png}%
\includegraphics[height=0.15\textwidth,width=0.1428\textwidth]{figs/3D3b.png}%
\includegraphics[height=0.15\textwidth,width=0.1428\textwidth]{figs/3D4a.png}%
\includegraphics[height=0.15\textwidth,width=0.1428\textwidth]{figs/3D5a.png}%
\includegraphics[height=0.15\textwidth,width=0.1428\textwidth]{figs/3D6b.png}%
\includegraphics[height=0.15\textwidth,width=0.1428\textwidth]{figs/3D7a.png}

{\footnotesize\hspace{.06\textwidth}(a)\hfill(b)\hfill(c)\hfill(d)\hfill(e)\hfill(f)\hfill(g)\hspace{.06\textwidth}}

\caption{Extraction of a minimal 2-cycle from $\mathcal{A}(X_2)$: (a) initial (0-th) value for $c\in C_2$; (b) cyclic subgroups on  $\delta\partial c$; (c)~1-st value of $c$; (d) cyclic subgroups on $\delta\partial c$; (e) 2-nd value of $c$; (f) cyclic subgroups on $\delta\partial c$; (g) 3-rd value of $c$, such that $\partial c=0$, hence stop.}
   \label{fig:3D}
\end{figure*}

The topological method introduced here, reminiscent of the ``gift-wrapping'' algorithm~\cite{Cormen:2009:IAT:1614191,Jarvis:1973:ICH} for computing convex hulls of 2D and 3D discrete sets of points, is detailed and formalized in Section~\ref{partial_3-computation}. The TGW algorithm takes a sparse matrix $[\partial_{d-1}]$ as input and produces in output a sparse matrix  $[\partial_{d}^+]$, augmented with the outer cell. A geometric embedding function $\mu: X_0\to\mathbb{E}^d$ is used to compute the angular ordering, around some $(d-2)$-cells, of $(d-1)$-basis elements in the boundary's coboundary, while wrapping up a $(d-1)$-cycle, as illustrated in Figures~\ref{fig:step-by-step} and \ref{fig:3D}. The built (minimal) cycles are set as columns of $[\partial_d]$, in the construction of the $C_d$ basis.
Note also that, once the ordered sets of basis elements is fixed, columns contain the \emph{coordinate representation} of $(d-1)$-cycles, built from group coefficients $(\{-1,0,1\},\,+) \simeq \mathbb{Z}/{3}\mathbb{Z} = \mathbb{Z}_3$. Analogously,  boundaries and coboundaries of chains are calculated by multiplication of operator matrices times proper coordinate vectors of such coefficients.

\subsection{Non-connected components (Holes)}

The outer cell of the space arrangement $X = \mathcal{A}(\mathcal{S})$ might have a non-connected boundary, compraising more than one $(d-1)$ cycle\footnote{Called a \emph{shell} in the literature of solid modeling.}, like a  3-ball minus a smaller concentric 3-ball. Analogously, $X$  might contain both non-connected and possibly nested components. The TGW algorithm actually computes all of the boundary cycles, that must be properly handled in order to produce a single $[\partial_d]$ matrix:  first decompose the input $[\partial_{d-1}]$ into connected components; then assemble/remove the empty cycles.

\subsubsection{Decomposition of 2-skeleton}\label{decomposition}

Consider a bipartite graph {$G = (N, A)$, with $N\simeq\Lambda_2 \cup \Lambda_0$, and $A \subseteq \Lambda_2 \times \Lambda_0$,} associated with the sparse characteristic matrix encoding the \emph{incidence} relation. $G$ has one node for each 2-cell, one node for each 0-cell, and one arc for each incident pair. Therefore, the arcs in $G$ are one-to-one with the nonzero elements of the $A$ matrix. By computing the {maximal} point-connected components of $G$, we subdivide the $X_2$ skeleton into $h$ connected components: $\mathcal{X}_2 = \{ {X}_2^p \}$, such that {$1\leq p\leq h$.}
For each component ${X}_2^p$, repeat the following actions. First, assemble the $[\partial_2]^p$ sparse matrix, and compute the corresponding $[\partial_3^+]^p$ generated by Algorithm~\ref{alg:one}. Then, subdivide $[\partial_3^+]^p$ into the boundary operator $\partial_3^p: C_3^p\to C_2$ and  the column matrix $c^p = \partial_3^+[\sigma^p] \in C_2$ of the outer cell $\sigma^p\in X_3$.
The set $S = \{ c^p \}$ of $h$ disjoint 2-cycles, is the initialization of the set of $X_d$ \emph{shells}.
Other (empty) shells of $X_d$ can be discovered later from mutual containment of $S$ elements.

\begin{figure}[htbp]
   \centering
   \includegraphics[width=\linewidth]{figs/nested_cycles}
   \caption{Non intersecting cycles within a 2D cellular complex with three connected components and only three cells, denoted by the image colors: (a) cellular complex; (b) graph of the \emph{reduced} containment relation $R$ between shells, with dashed arcs of even depth index. Removing the dashed arcs produces a \emph{forest} of small \emph{trees}; (c) matrix of transitively reduced $R$. Note that the ones {equate} the number of edges in the graph. }
   \label{fig:shells}
\end{figure}

\subsubsection{Shell discovering and assembling/removing}\label{shells}

Then add to the set $S$ the empty cycles (shells) already included within some non-contractible cell. We  call them \emph{input holes}. We discover the (unmodified) input holes by direct inspection of matrices $[\partial_3^p]$, since holes are represented there as \emph{pairs of columns with zero sum}. In fact, each input hole, if non intersected by other data, returns unmodified, and produces \emph{two opposite  columns} in the component matrix $[\partial_3^p]$ it belongs to. Each corresponding pair of columns has all non-zero elements on the same rows, but with opposite signs (orientations).
The inspection is done by a sort of scan-line algorithm working on the rows of each $[\partial_3^p]$ matrix, that recognizes the emerging pairs of opposite columns, and stores each pair (candidate hole) until its contents eventually differ. The  algorithm proceeds by moving from a row to the  next one, until the bottom row is reached and the set, possibly empty, of holes is returned. One column from each pair (with proper sign---see~\ref{poset}) is removed from the component matrix $[\partial_3^p]$, and its opposite 2-cycle is added to the set $S$ of shells.

\subsubsection{Transitive reduction of shell poset}\label{poset}

The antisymmetric \emph{containment relation} between 2-cycles in $S$ is computed for $i,j\in [1,\ldots,h], i<j$, by the containment test $\mathit{PointInPolygon}(u_0^i,c^j)$ between a single unit chain $u_0^i \in c^i$ and the cycle $c^j$, with $(c^i,c^j)\in S$. Then, the transitive reduction $R$ is extracted.
If the edge-set of the $R$ tree is empty, no disjoint component of $X_3$ is contained inside another one, and both $X_3$ and $\partial_3$ may by assembled by disjoint union of 3-cells of $X_3^p$ and columns of $[\partial_3]^p$, $(1\leq p\leq h)$, respectively.
If, conversely,  the above is not true, then reduce $R$ to a forest of small trees by cancelling the arcs at even distance from root (see Figure~\ref{fig:shells}) and, for each arc $(i,j)$, discover which \emph{cell} of the container component $X^i$ actually contains the contained component $X^j$, \emph{i.e.}, its shell $c^j$ and its possibly non-empty interior. More complex intersecting situations are impossible by construction, since the components are a priori disjoint. Therefore, in case of containment, one component is necessarily contained in some \emph{empty cell} of the other.

\subsection{Computational interlude}\label{sec:prologue}

For reader's convenience, we give in Appendix~\ref{sec:appendix} definitions and facts about computing with chains and cochains, and provide some examples of elementary topological computation with chains, cochains and their operators.
Be it noted that our \emph{topology representations}, \emph{i.e.},
\begin{itemize}
\item 1-,2-,3-cells as either 0-chains or p-cycles ($p=1,2$), as \emph{sparse 1-arrays} of signed numbers,
\item 1-,2-,3-complexes as bases of linear spaces of $p$-chains and graded linear transformations, represented as \emph{sparse 2-arrays} of signed numbers,
\end{itemize}
have smaller space complexity than common data structures of well-known efficient  representations~\cite{Woo:85,Dicarlo:2014:TNL:2543138.2543294} of Solid Modeling. For example, with regard to the B-rep of a closed 2-manifold $A$, we have $\mathit{Space(A)} = \mathit{Space([\partial_2])}+\mathit{Space([\partial_1])} = 2\#\texttt{E}+2\#\texttt{E}$, where $[\partial_2], [\partial_1]$ are sparse matrices of boundary operators, and $\#\texttt{E}$ is the number of unit 1-chains (edges). Hence the storage required by $\mathit{Space(A)}$ is equal to $2/3$ of \emph{half-edge}~\cite{Muller:78}, largely used in Computational Geometry, and $1/2$ of \emph{winged-edge}~\cite{Baumgart:1972:WEP:891970}, often used as a reference representation for manifold Solid Modeling~\cite{Ala:1992:PAB:616022.617736}.

The cardinality of \emph{all} the incidence/adjacency relations between p-cells and q-cells $(1\leq p,q \leq 3)$ is also minimal, according to~\cite{Woo:85}.  For example, $\#\texttt{FV} = O(Space([\partial_2] [\delta_1]) = 2\#\texttt{E}$), where \texttt{V} are the vertices of a complex and \texttt{EV}, \texttt{FV} are binary incidences of edges and faces with vertices. Therefore, every set of local queries about the $3\times 3$ incidences/adjacencies between $p$-cells can be answered by multiplication, via software kernels for \emph{sparse matrix} product and transposition, just by collecting the coordinate vectors of unit chains, ``subject''  of elementary queries, as \emph{columns} of a sparse $Q$ matrix, and by left-multiplying $Q$ times one/two operator matrices $[\partial_1]$ and/or $[\partial_2]$, suitably ordered and/or transposed~\cite{Dicarlo:2014:TNL:2543138.2543294}, to get the algebraic equivalent of multiple database queries at once.

\section{Chain-based arrangement algorithms}
\label{sec:algorithms}

In this section, we provide a slightly simplified pseudocode\footnote{The actual implementation is available as open-sourced Julia package in Github: \href{https://github.com/cvdlab/LinearAlgebraicRepresentation.jl}{ https://github.com/cvdlab/LinearAlgebraicRepresentation.jl}.} of the main algorithms introduced in the previous section,
and discuss their worst-case complexity.  The pseudocode style {is a blend of} Python and Julia styles.
A few words about  notations:  greek letters are used for the \emph{cells} of a space partition, and roman letters for  \emph{chains} of cells, all coded as either signed integers or sparse arrays of signed integers.  $[\partial_d]$ or $[c]$ stand for general matrices or column matrices, {respectively}, whereas $\partial_d[h,k]$ or $c[\sigma]$ stand for their indexed elements. Also, $|c|$ stands for \emph{unsigned} (nonzero) indices of the (sparse) array $[c]$.
The \emph{accumulated assignment} statement $A\ $+=$\ B$ stands for $A = A + B$, where the meaning of ``$+$'' symbol depends on the {context}, \emph{e.g.} may stand either  for sum (of chains), or for union (of sets), or for concatenation (of matrix columns). Analogously, $A\ $-=$\ B$ stands for $A = A - B$.

\subsection{Arrangements of 2-cells (Merge algorithm)}\label{arrangement-of-2-cells}

The sequence of computations performed on each 2-cell $\sigma\in\mathcal{S}_2\subseteq \mathcal{S}$ is discussed in the following. A visualization of the decomposition process, discussed in Example~\ref{ex:square}, is shown in Figure~\ref{fig:abutting}.

\subsubsection{Fragmentation of a 2-cell}\label{fragmentation}

In a first stage, the subset $\mathcal{I}(\sigma)$ of 2-cells of \emph{potential intersection} with $\sigma$ is computed (see Figure~\ref{fig:process}c), by intersecting the results of three queries upon 1D \emph{interval trees}, generated at the beginning of the pipeline.
Each 1D interval-tree is based on one dimension of the 3D containment boxes of input 2-cells. Then, the set $\Sigma = \{\sigma\}\cup\mathcal{I}(\sigma)$ is transformed so that $\sigma$ is mapped into the $z=0$ subspace (Figure~\ref{fig:process}d).
In $\mathbb{E}^3$ the mapped 2-cells in $\Sigma\cap\{(x,y,z)|z=0\}$ are used to compute a \emph{set of line segments} in $\mathbb{E}^2$, generated by intersection of edges of 2-cells in $\Sigma$ with the 2D plane, and by join (convex combination) of alternate pairs\footnote{The 2-cell being intersected with $z=0$ may be non-convex, and its 1-cells may produce and even number $k>2$ of intersection points along the same line.} of intersection points of boundary edges, along the intersection line of each 2-cell in $\Sigma$ (Figures~\ref{fig:process}e and~\ref{fig:2D}a).

\begin{algorithm}[h]
\SetAlgoNoLine
\DontPrintSemicolon
\KwIn{$\mathcal{S}_2\subseteq \mathcal{S}_{d-1}$ $\qquad$ { \emph{collection} of all 2-cells from $\mathcal{S}_{d-1}$ input in $\mathbb{E}^{d}$}}
\KwOut{$[\partial_2]$ $\qquad$ {signed CSC (Compressed Sparse Column)  matrix}}
$\widetilde{\mathcal{S}_2} = \emptyset$ $\qquad$ {{ initialisation of collection of local fragments }}\;
\For($\qquad$ {{ for each 2-cell $\sigma$ in the input set}}){$\sigma \in \mathcal{S}_2$}{
	$M = \mathit{SubManifoldMap}(\sigma)$ $\qquad$ {{ affine transform \emph{s.t.} $\sigma\mapsto x_3=0$ subspace}}\;
    $\Sigma = M\,(\mathcal{I}(\sigma)\cup\{\sigma\})$ $\qquad$ {{ apply the map $M$ to (possible) incidencies to $\sigma$}}\;
	$\mathcal{S}_1(\sigma) = \emptyset$ $\qquad${{ collection  of line segments in $x_3=0$}}\;
    \For($\qquad$ {{ for each 2-cell $\tau$ in $\Sigma$}}){$\tau \in \Sigma$}{
        $\mathcal{P}(\tau),\mathcal{L}(\tau) = \emptyset,\emptyset$ $\qquad$ {{ intersection points and segment(s) with $x_3=0$}}\;
        \For($\qquad$ {{ for each 1-cell $\lambda$ in $X_1(\tau)$}}){$\lambda \in X_1(\tau)$}{
        	\lIf($\quad$ {{ save intersection of $\lambda$ and $x_3=0$}}){$\lambda\not\subset\{\p{q}\,|\,x_3(\p{q})=0\}$}{
				$\mathcal{P}(\tau)\ $+=$\ \{\p{p}\}$}
		}
		$\mathcal{L}(\tau) = \mathit{Points2Segments}(\mathcal{P}(\tau))$ $\quad$ {{ set of collinear intersections}}\;
		$\mathcal{S}_1(\sigma)\ $+=$\ \mathcal{L}(\tau)$ $\quad$ {{ accumulate intersection segments generated by $\tau$}}\;
	}
	$X_2(\sigma) = \mathcal{A}(\mathcal{S}_1(\sigma))$ $\quad$ {{ arrangement of $\sigma$ induced by a set of 1-complexes}} \;
	$\widetilde{\mathcal{S}_2}\ $+=$\ M^{-1}\,X_2$  $\qquad$ {{ accumulate local fragments, back transformed in $\mathbb{E}^d$ }}\;
}
$[\partial_1] = \mathit{QuotientBases}(\widetilde{\mathcal{S}_2})$  $\qquad$ {{ identification of 0- and 1-cells using $kd$-trees  }}\;
$[\partial_2] = \mathit{TGW}([\partial_1])$ $\qquad$ {{ output computation via TGW algorithm in 2D}}\;
\Return{$[\partial_2]$}
\caption{ \emph{{Subdivision} of 2-cells }}
\label{alg:two}
\end{algorithm}

The {planar processing of the 2-cell $\sigma$} continues by  pairwise executing the line segment intersection, producing a linear graph, as shown in Figure~\ref{fig:2D}b. The dangling~\footnote{In a $d$-complex, \emph{dangling edges} are $p$-cells, p<d, that are not contained in any boundary cycle of a $d$-cell. They are the interior structures of SGC cells. In Solid Modeling terminology, they are called non-regular subsets, whence the term \emph{regularized} Boolean operation.} edges are removed, by computing the maximal 2-vertex-connected subgraphs\footnote{A connected graph G is \emph{2-vertex-connected} if it has at least three vertices and no articulation points.
A vertex is an \emph{articulation point} if its removal increases the number of connected components of G.}~\cite{vialar2016handbook}, with the Hopcroft's and Tarjan's algorithm~\cite{Hopcroft:1973:AEA:362248.362272}.
Only the non-external biconnected components enter the following computations, since the other graph parts are either external to $\sigma$ or certainly dangling (1-connected subgraphs), and will  contribute separately to the space arrangement (Figure~\ref{fig:2D}b).
{Finally}, the oriented 2-chain of the partition $\mathcal{A}(\Sigma)$ is computed  as shown in Figures~\ref{fig:process}g and ~\ref{fig:2D}d, using the TGW in 2D, so generating the $\partial_2(\sigma)$ matrix from $X(\sigma)$.
The fragmentation process is repeated for each $\sigma\in \mathcal{S}_2$, with each geometric map $\mu(\sigma): X_0(\sigma)\to\mathbb{E}^2$ composed with its inverse transformation back to $\mathbb{E}^3$.

\subsubsection{Complexity of 2-cells {subdivision}}

The time complexity of Algorithm~\ref{alg:two} is bounded by the number $n$ of 2-cells in the input collection $\mathcal{S}_{d-1}$ times the worst-case cost required by the {subdivision} of one of them. In turn, this cost depends on the size of the actual input, \emph{i.e.}, on the number of potentially intersecting 2-cells in $ \mathcal{I}(\sigma)$. In all the regular cases we usually meet in computer graphics, CAD meshes and engineering applications, the number of 2-cells incident in (even on the boundaries of) a given cell is bounded by a constant number $k_1$. If the maximum number of 1-cells on the boundary of a 2-cell is $k_2$, then the whole computation of  Algorithm~\ref{alg:two} requires time $O(k_1 k_2 n + A)$, where  $A$ is the time needed to compute the quotient sets, \emph{i.e.},~to glue all  $X_2(\sigma)$ in $\mathbb{E}^d$ space.  When $d=3$, the affine transformations $\T{Q}_\sigma$ of each set $\Sigma$ (see~\ref{facet-arrangements}) are computable in $O(1)$ time; building a static $kd$-tree generated by $m$ points requires $O(m \log^2 m)$; and each query for finding the nearest neighbor in a balanced $kd$-tree requires $O(\log m)$ time on average. The number of occurrences of the same vertex on incident 2-cells is certainly bounded by a small constant $k_3$, approximately equal to $m/v$, where $v=\#X_0$ is the number of 0-cells after the identification processing. The transformation of LAR output to canonical form (sorted 1-array of integers)  is done in $O(1)$ for each edge, so giving $A = O(m \log^2 m) + O(m\log m) + O(1) = O(m \log m)$. In conclusion, {the worst-case running time of Algorithm~\ref{alg:two} is}  $O(k_1 k_2 n + m \log m)$.

\subsection{Quotient sets computation (Congruence algorithm)}\label{quotients}

Small sparse matrices of signed operators $\partial_2(\sigma): C_2(\sigma)\to C_1(\sigma)$ have  already been assembled \emph{independently} in 2D for each fragmented $\sigma$, \emph{i.e.}, for each $X_2^\sigma$, as detailed in the previous Section~\ref{fragmentation}. The output of that pipeline stage is a collection  $\varmathbb{C} := \{ C_\bullet(\sigma), \sigma\in \mathcal{S}_2\}$ of small chain complexes, one for each input 2-cell, embedded in $\mathbb{E}^3$.
They were buil by repeatedly applying in 2D the TGW algorithm (see Section~\ref{partial_3-computation}) and mapping back the results in 3D.
The quotients of chain spaces modulo the $p$-congruence relations are calculated at this point, starting from $p=0$. Therefore, the unit 0-chains are identified numerically via their geometric maps and snap rounded by numerical identification of nearby-coincident points using a $kd$-tree.  The congruent unit 1-chains and 2-chain are identified symbolically, making use of their {unique} canonical indexed representation. The canonical representation of a unit $d$-chain is the array of sorted indices of the unit elements of its $(d-1)$-cycle.  The 2-cells of the output 2-complex $X_2(\mathcal{S}_2)$ embedded in $\mathbb{E}^3$, written as 1-cycles, \emph{i.e.},~as linear combinations of signed 1-cells, are stored by column in the matrix of the operator $\partial_2: C_2\to C_1$. A 1-cell $\tau$, is written\footnote{As a 0-chain of signed 0-cells in the matrix representation of $\partial_1: C_1\to C_0$.} by convention as $1u_0^{i_k} -1u_0^{i_h}$ {when $k>h$},  and is oriented from $u_0^{i_h}$ to $u_0^{i_k}$.  The conventional rules used in this paper about sign and orientation of cells are summarized at the end of Section~\ref{sec:thechains}.

\begin{algorithm}[h]
\SetAlgoNoLine
\DontPrintSemicolon
\tcc{{\rm\textbf{Pre-condition:} $d$ equals the space $\mathbb{E}^d$ dimension, such that~$(d-1)$-cells are shared by \emph{two} $d$-cells}}
\tcc{{\rm ~}}
\KwIn{$[\partial_{d-1}]\qquad~$ { signed CSC matrix} $(a_{ij})$, where $a_{ij}\in\{ -1,0,1 \}$}
\KwOut{$[\partial_d^+]\qquad~$ { signed CSC matrix of $(d-1)$-cycles}}
$[\partial_d^+] = []\ $; $m,n = [\partial_{d-1}].shape\ $; $marks = Zeros(n)\ $ $\qquad~$ {{ initializations}}\;
\While{\emph{Sum}$(marks) < 2n$}{
	$\sigma = Choose(marks)$ $\qquad$ {{ select the $(d-1)$-cell seed of column extraction\;}}
	\lIf{$marks[\sigma] == 0$}{$[c_{d-1}] = [\sigma] $}
	\lElseIf{$marks[\sigma] == 1$}{$ [c_{d-1}] = [-\sigma]$}
	$[c_{d-2}] = [\partial_{d-1}]\,[c_{d-1}]$ $\qquad$ {{ compute boundary $c_{d-2}$ of seed cell\;}}
	\While($\qquad$ {{ loop until boundary becomes empty}}) {$[c_{d-2}] \not= []$}{
		$corolla = []$\;
		\For($\qquad$ {{ for each ``hinge'' $\tau$ cell}}){$\tau \in c_{d-2}$}{
			$[b_{d-1}] = [\tau]^t [\partial_{d-1}]$  $\qquad$ {{ compute the $\tau$ coboundary}}\;
			$pivot = \{|b_{d-1}|\} \cap \{|c_{d-1}|\}$  $\qquad$ {{ compute the $\tau$ support}}\;
			\lIf($\qquad$ {{ compute the new adj cell}}){$\tau > 0$}{$adj = Next(pivot, Ord(b_{d-1}))$}
			\lElseIf{$\tau < 0$}{$adj = Prev(pivot, Ord(b_{d-1}))$}
			\lIf($\qquad$ {{ orient adj}}){$\partial_{d-1}[\tau,adj] \not= \partial_{d-1}[\tau,pivot]$}{$corolla[adj] = c_{d-1}[pivot]$}
			\lElse{$corolla[adj] = -(c_{d-1}[pivot])$}
		}
		$[c_{d-1}]\, $+$= corolla$  $\qquad$ {{ insert $corolla$ cells in current $c_{d-1}$\;}}
		$[c_{d-2}] = [\partial_{d-1}]\,[c_{d-1}]$  $\qquad$ {{ compute again the boundary of $c_{d-1}$\;}}
	}
	\lFor($\qquad$ {{ update the counters of used cells}}) {$\sigma\in c_{d-1}$}{$marks[\sigma]\ $+=$\ 1$}
	[$\partial_d^+]\ $+$= [c_{d-1}]$  $\qquad$ { append a new column to $[\partial_d^+]$ }
}
\Return{$[\partial_d^+]$}
\caption{ \emph{Computation of signed $[\partial_d^+]$ matrix} }
\label{alg:one}
\end{algorithm}

\subsection{\texorpdfstring{Computation of \(\partial_2\) and \(\partial_3\)  (TGW Algorithm)}{\textbackslash{}delta\_1 computation (2 pages)}}\label{partial_3-computation}

\subsubsection{Topological Gift Wrapping}

The algorithm was introduced in Section~\ref{tgw-algorithm}. Here we provide a readable pseudocode, with the only \emph{caveat} that it actually computes a redundant set of generators for $C_3$ (resp.~$C_2$),  as minimal connected $2$-cycles (resp.~$1$-cycles) from a $\partial_2$ (resp.~$\partial_1$)  matrix.
A step-by-step formalized example of computation of a unit 2-chain as 1-cycle, using the TGW algoritm, is discussed in Example~\ref{example1}.
The given pseudocode makes use of math symbols and high-level math operations;  its actual implementation in scientific languages like Python or Julia uses sparse arrays and discrete coordinates in $\{-1,0,1\}$, to achieve an efficient execution in terms of storage space and computation time.

Note the precondition of Algorithm~\ref{alg:one}, warning that the method used will compute the $\partial_d$ matrix only for a cell decomposition of $d$-space. In fact, only in this case the $(d-1)$-cells are shared by exactly \emph{two} $d$-cells, including the outer cell.
This condition implies that the input cellular complex it applies to should be a (possibly non-connected) CW-complex, with all cells homeomorphic to spheres. Note also that the matrix of the boundary operator for the $d$-chain space of a cellular complex with holes as well as inner and outer components will be built starting from output of Algorithm~\ref{alg:one} in a later stage. The termination predicate of Algorithm~\ref{alg:one} is a consequence of the above  property: the algorithm terminates when all incidence numbers in the $marks$ array are 2, so that their sum is exactly $2n$, where $n$ is the number of $(d-1)$-cells, equal to the number of columns in the input matrix $[\partial_{d-1}]$.

\subsubsection{Valid input and unique output}\label{sec:limitations}

{The algorithm works properly with legitimate input. In particular, input $(d-1)$-skeletons must be regular, \emph{i.e.}, without dangling parts, so that every $(d-1)$-cell belongs at most to \emph{two} $(d-1)$-cycles. In 2D, this fact is guaranteed by applying the algorithm separately to each 2-maximally connected component of the 1-skeleton, considered as a graph, and then by merging the results (clearly disconnected). Analogously, in 3D, the adjacency graph of 2-cells should not contain dangling subgraphs. }

The validity set of the input may contain 2-skeletons of 3-complexes, boundaries of solid models, sets of manifold boundary components of non-manifold solid models. Of course, to apply the algorithm to data which do not determine a partition of the embedding space does not make sense and produces an empty result.  We call them \emph{illegal data}.
When applied to valid input, as described above, the TGW algorithm produces valid output, because always produces a set of generators for $C_d$ that satisfies the Eq.~\ref{eq:valid} below:
\begin{equation}
[\partial_d] = (a_{ij})\quad \mbox{where}\quad
\sum_{i=1}^{\#X_{d-1}} \sum_{j=1}^{\#X_{d}} |a_{ij}|= 2\,(\#X_{d-1}),
\label{eq:valid}
\end{equation}
The results are also \emph{unique}, modulo reordering, since otherwise two different bases for the linear space $C_d$ would produce two boundary operators that, applied to the total 3-chain (vector of all ones) would return the same boundary cycle, which is impossible.
There are no ambiguities in the algorithm, since in a $d$-complex every two $d$-cells share at most two $(d-1)$-cells, or exactly two if the outer $d$-cell is considered. Also, it halts when this last condition is exactly reached.  Note that the suitable choice of the next ``petals'' from  ``corolla'' (see the pseudocode in Algorithm~\ref{alg:one}) implies that a 2-cell cannot be used more that twice.

\subsubsection{Complexity of 3-cell extraction}
\label{sec:complexity1}

In three dimensions, Algorithm~\ref{alg:one} constructs iteratively (outer \textbf{while}) one unit 3-chain (represented as a 2-cycle, \emph{i.e.}, as a closed 2-chain), building the corresponding column of the matrix $[\partial_3^+]$, and so adding one outer boundary column for each connected component of the input complex, as detailed in \ref{extension}.

The space complexity of a 3-cell is measured by a set of triples (row, column, value) for each non-zero values\footnote{{The coordinate (\texttt{COO}) representation of sparse matrices~\cite{gemmexp} is an array of triples $(i,j,value)$.}} in its column, \emph{i.e.},~with its representation as cycle of unit 2-chains. Hence, the total number of triples, \emph{i.e.},~the space complexity of the \texttt{COO} representation of $[\partial_3^+]$, is exactly $2n$, where $n$ is the number of 2-cells in the $X_2$ skeleton.

The construction of a single 3-cell requires the search of the adjacent $\mathit{adj}$ 2-cell for each $\mathit{pivot}$ unit 2-chain in the boundary shell. The search for $\mathit{next}$ or $\mathit{prev}$ 2-cell as $\mathit{adj}$ for each $\mathit{pivot}$ requires the circular sorting of this permutation subgroup of 2-cells incident to each 1-cell on each boundary of an incomplete 2-cycle. Consequently, we have several sorts of small sets, where each set is normally bounded by a very small integer, hence each sort is $O(1)$ timewise, and their total number is upper bounded by the number of $(d-1)$-cells on the $d$-cell boundary (equal to 6 for cubical 3-complexes, and to 4 for simplicial 3-complexes, and a small integer in general).

The subsets to be sorted are encoded in the columns of the incidence matrix from 2-cells to 1-cells, \emph{i.e.}, by the $i,j$ indices of  non-zero elements of $[\partial_2]$. The computation of the (unsigned) $[\partial_2]$ may be performed through SpGEMM\footnote{SpGEMM is a subroutine for matrix multiplication between two general sparse matrices~\cite{gemmexp}, \emph{i.e.}, no banded, nor Hermitian, \emph{etc.} Name derived from BLAS rules~\cite{Ballard:2015}.} multiplication of two sparse matrices (see~\cite{Dicarlo:2014:TNL:2543138.2543294}), hence in time linear with the size of the output, \emph{i.e.}, with the number of non-zero elements of the   $[\partial_2]$ matrix.  Summing up, if $n$  is the number of $d$-cells and $m$ is the number of $(d-1)$-cells, the time complexity of this algorithm is $O(n m\log m)$ in the worst case of unbounded complexity of $d$-cells, and roughly $O(n k\log k)$ if their $(d-1)$-cycle complexity is bounded by $k$.

\subsection{Isolated shells (Holes algorithm)}\label{extension}

In the general case, both the outer cell and the inner cells may contain holes and/or isolated components. Talking of isolated holes is improper, since holes  are not empty within an arrangement, \emph{i.e.}, a partition of the ambient space, but contain an isolated component within their boundary represented as a ($d-1$)-cycle. The aim of this section is to discuss the handling of isolated components and their boundaries, to be considered holes within their container cells.

The TGW Algorithm~\ref{alg:one} produces CW-complexes, despite the fact that the subdivision of Merge Algorithm~\ref{alg:two} may generate non contractible $d$-cells, \emph{i.e.}, cells with holes {(but without isolated points). These spaces are handled by combining standard CW-complexes, \emph{i.e.}, with cells homeomorphic to spheres, and by adding $d$-cells to the interior of $d$-cells. }
In other words, the boundary of holes in a cell coming from disconnected sub-complexes is merged in the container cell. The orientation is handled depending on the parity of relative containment relation. The management of isolated boundaries concerns essentially the adjoining/removal of columns in the final boundary matrix.

\subsubsection{Synthesis of the whole pipeline}\label{preview-tgw}

We need to consider two main issues: (a) the computation of {maximal} connected components of $X_{d-1}$ may produce $h>1$ {disconnected} \emph{$d$-components} of the output complex $X_d$; (b) the inclusion of components within single cells of the output $d$-complex: see Figure~\ref{fig:shells}.  In the following we list the main stages of  the \emph{Holes} algorithm to take care of these issues. Our goal is the computation of both the $X_d$ skeleton, and the $\partial_d$ operator for spaces with multiple components nested into holes.  Note that exactly the same points apply (scaled in dimension) before and after TGW execution, for both \emph{local} arrangements in 2D and the \emph{global} arrangement in 3D, respectively. This fact gives a cue for a possible multidimensional extension.

\subsubsection{\texorpdfstring{Non-intersecting shells}{\textbackslash{}delta\_2 computation (2 pages)}}\label{delta_2-computation}

If the shell-set $S \not= \emptyset$, then the $h$ \emph{isolated} boundary components $(0\leq p\leq h)$ in $S$ must be compared with each other, to determine their relative containment, if any, and consequently their orientation.
The $h\times h$ binary and antisymmetric matrix $M=(m_{ij})$ of the relation is built, by computing each element $m_{ij}$ $(i<j)$, with a single point-cycle ray firing, because the two corresponding cycles (columns $i$ and $j$) are guaranteed not to intersect.
The attribute of $c_j$ as outer/inner, and hence its relative orientation is given by the parity of $c^j$ in $R$.
When the cancellation \ref{preview-tgw}.4 of empty cells has been performed for all ``solid'' arcs of $R$ (see Figure~\ref{fig:shells}), the updated matrices $[\partial_d]^p$ can be assembled into the final $\partial_d$ operator matrix.

\begin{algorithm}[htbp]
\SetAlgoNoLine
\DontPrintSemicolon
\KwIn{\texttt{LAR}${}_{d-1}$, $[\partial_{d-1}]\qquad\qquad~$ { for $d=3$: \texttt{FV}, $\partial_2$}}
\KwOut{\texttt{LAR}${}_{d}$, $[\partial_d]\qquad\qquad~$ { for $d=3$: \texttt{CV}, $\partial_3$}}
$N = \Lambda_2 \cup \Lambda_0$; $\quad A \subseteq \Lambda_2\times\Lambda_0$; $\quad G = (N,A)$ $\qquad\qquad~$ {{ initializations}}\;
$\mathcal{G} = \{ G^p\ |\ 1\leq p\leq h \} \leftarrow \mathit{ConnectedComponents}(G)$  $\qquad\qquad~$ {{ partition of $G$ into $h$ connected components}}\;
$\mathcal{X}_{d-1} = \{ (X_{d-1}^p, \partial_{d-1}^p)\ |\ 1\leq p\leq h \} \leftarrow \mathit{Rearrange}(\mathcal{G})$  {$\qquad\qquad~$ {{ partition of $X_{d-1}$ into $h$ connected components}}}\;
$S = []$ $\qquad\qquad~$ {{ initialize the sparse array of \emph{shells}}}\;
\For({$\qquad${ for each connected component of $(d-1)$-skeleton}}){$p \in \{ 1,\ldots, h \}$}{
	$[\partial_d^+]^p = \mathit{Algorithm\_1}([\partial_{d-1}]^p)$ $\qquad$ {{ compute the minimal $d$-cycles of a component of complex}}\;
	$(c^p, \partial_d^p) = \mathit{Split}([\partial_{d}^+]^p)$ $\quad$ {{ split the component into the exterior $(d-1)$-cycle and the boundary $\partial_d^p$}}\;
	$S\ $+=$\ [c]^p$ $\qquad\qquad~$ {{ append the boundary shell to the shell array}}\;
}
\For({$\qquad${ for each shell pair $(c^i,c^j)\in R$,\ }}){$i,j \in \{ 1,\ldots, h \},i<j$}{
	$(R[i,j],R[j,i])::\mathit{Bool}\times\mathit{Bool} \leftarrow \mathit{PointSet}(u_0^i\in c^i,c^j)$ {$\qquad${ containment test of $u_0^i$ in $c^j$ }}\;
}
$R = \{ (i,j) \} \leftarrow \mathit{Tree(TransitiveReduction}(R\,))$ {$\qquad${ set of arcs of reduced containment tree of shells }}\;
\If({$\qquad${ if the containment tree of shells is not empty}}){$R \not=\emptyset$}{
	\For({$\qquad${ for each shell pair $(c^i,c^j)$ such that $dist(c^j)\%2\ $!=$\ 0$}}){$(i,j) \in R$}{
		$\rho = \mathit{FindContainerCell}(u_0^i,c^j,\texttt{LAR}_{d-1})$ {$\qquad${ look for a $d$-cell $\rho$ such that $u_0^i\in  |c^i| \subseteq |\rho|\subseteq |c^j|$ }}\;
		$[\partial_d]^j\ $-=$\ \partial_d^j[\rho]$ {$\qquad${ remove $\rho$ from $\partial_d^j$ }}\;
	}
	}
$\partial_d = [\partial_d^1 \cdots \partial_d^p \cdots \partial_d^h]$ {$\qquad${ return the aggregate $\partial_d$ operator }} \;
$\texttt{LAR}_{d} = [\cup_k \texttt{LAR}_{d-1}(c^k = \partial_d[\cdot,k]),\ \mathit{for\ } k\in\mathit{Range}(Cols(\partial_d))]\qquad$ {{ for $d=3$: \texttt{LAR}${}_{d} = \texttt{CV}$}}\;
\Return \texttt{LAR}${}_{d}$, $[\partial_d]$

\caption{ \emph{Non-intersecting shells}}
\label{alg:four}
\end{algorithm}

\subsubsection{Complexity of shell management}
\label{sec:inclusion}

The computation of the connected components of a graph $G$ can be performed in linear time~\cite{Hopcroft:1973:AEA:362248.362272}. The recognition of the $h$ shells requires the computation of $[\partial_{d}^+]^p$ ($1\leq p\leq h$) and the extraction of the boundary of each connected component $X_{d}^p$ . To compute the reduced relation $R$ we execute $O(h^2)$ point-cycle containment tests, linear in the size of a cycle, so spending a time $O(h^2 nr)$, with $h$ number of shells, and $n$ average size a cycle. Actually, the point-cycle containment test can be easily computed in parallel, with a minimal transmission overhead of the arguments.
The restructuring of boundary submatrices has the same cost of the read/rewrite of columns of a sparse matrix, depending on the number of non-zeros of $[\partial_3]$, and hence is $O(n \#X_d)$, \emph{i.e.}, linear with the product of the number of $d$-cells and their average size $n$ as chains of $(d-1)$-cells, with $n$ size of the average isolated cycle.

\subsection{The whole picture}\label{whole-picture}
A short synthesis of sequential steps of the whole computational pipeline, from input collection to chain complex output, follows in the more general setting, with both isolated components (within the outer cell), and possibly nested isolated components (within holes in inner cells).
\begin{description}
    \item[Input]  Facet selection, \emph{i.e.}, construction of the collection $\mathcal{S}_{d-1}$ from $\mathcal{S}_d$, using \texttt{LAR}.
    \item[Indexing]	Spatial index made by intersection of  $d$ interval-trees on bounding boxes of $\sigma\in\mathcal{S}_{d-1}$.
    \item[Decomposition]  Pairwise $z=0$ intersection of line segments in $\sigma\cup\mathcal{I}(\sigma)$, for each  $\sigma\in\mathcal{S}_{d-1}$.
    \item[Congruence]	Graded bases of equivalence classes $C_k(U_k)$, with $U_k=X_k/R_k$ for $0\leq k\leq 2$.
    \item[Connection]	Extraction of $(X_{d-1}^p,\partial_{d-1}^p)$,  maximal connected components  of $X_{d-1}$ ($0\leq p\leq h$).
    \item[Bases]	Computation of redundant cycle basis $[\partial_{d}^+]^p$ for each $p$-component, via TGW.
    \item[Boundaries]	Accumulation into $H +{}\!\!\!= [o]^p$ (hole-set) of outer boundary cycle from each $[\partial_d^+]^p$.
    \item[Containment]	Computation of antisymmetric containment relation $S$  between $[o]^p$ holes in $H$.
    \item[Reduction]	Transitive $R$ reduction of $S$ and generation of forest of flat trees $\langle[o_d]^p, [\partial_{d}]^p\rangle$.
    \item[Adjoining]	of  roots $[o_{d}]^r$ to (unique) outer cell, and non-roots $[\partial_{d}^+]^q$ to container cells.
    \item[Assembling]	 Quasi-block-diagonal assembly of matrices relatives to isolated components $[\partial_d]^p$.
    \item[Output] Global boundary map $[\partial_d]$ of $\mathcal{A}(\mathcal{S}_{d-1})$, and reconstruction of 0-chains of $d$-cells in $X_d$.
\end{description}

\section{Relevant literature}\label{other-works}

In this section we mention some relevant connections of the present approach with recent papers concerning similar topics, and discuss some remarks in relation to our own work.

In~\cite{Alayrangues2015} a topological approach to homology is introduced for subclasses of subdivided spaces constructed by combinatorial and generalized maps.  Generalized map  (\emph{Gmap}) is a combinatorial model which allows for representing and handling subdivided objects~\cite{Lienhardt:2014} via  ``connecting darts'' between cell pairs. Gmaps are used to describe the topology of manifold-like cellular objects where $p$-cells are homeomorphic to $p$-spheres.
Alayrangues, Damiand, Lienhardt, and Peltier give an algorithm to build signed boundary maps for $0\leq i\leq 3$, focusing on the equivalence between computing homology via Gmaps and via simplicial complexes.

The time complexity of boundary maps in~\cite{Alayrangues2015} is linear in the number of incidence numbers. In~\cite{Dicarlo:2014:TNL:2543138.2543294}, we already obtained the same result, linear in the sparse output size, for computation via SpGEMM multiplication (see, \emph{e.g.},~the Example~\ref{ex:example0}) when the input complex is known. The   complexity of TGW for computing the \emph{unknown} $X_d = \mathcal{A}(\mathcal{S}_{d-1})$ is obviously higher (see Section~\ref{sec:complexity1}), and equates the standard in Solid Modeling. The main difference with our approach is that Alayrangues and colleagues start from a \emph{given} Gmap cellular model, whose construction is quite complex and requires interactive operations with a graphical user interface or a symbolic logic systems (such as INRIA's Coq~\cite{DEHLINGER2014869}) with a formal specification language. If simplicity metric matters, our linear algebraic representation of chains with sparse arrays compares well with chains of Gmaps.

With Selective Geometric Complex~\cite{Rossignac:1991:CNG:115604.115606}, Rossignac and O'Connor proposed a significant extension of topological CW-complexes, allowing for $p$-cells with  structures of dimension $0\leq k<p$ internal to cells, \emph{i.e.}~not necessarily embedded in a  cell boundary. The  selection bit associated to each cell allows selective choice of substructures. The association of incidence and local ordering among incident cells is maintained via hierarchical links, analogous to the Hasse diagram between vertices, edges, and faces, as well with geometric extents of non-linear surface patches.
Two attributes for \emph{c.boundary} and \emph{c.star} of a cell return the cells on its boundary and those it is on the boundary of. The search for boundary of more complex substructures is algorithmic.

In the present paper, we conversely use graded and combinable linear operators for  boundary and coboundary to traverse, both locally and globally, the incidence hierarchy. Hence we obtain, via multiplication of sparse matrices and vectors, a complete linear characterization of the space topology. Even local updates to topology, via Euler operators, can be done algebraically~\cite{ieee-tase}.  Another significant difference concerns the large amount of information and pointers associated by SGC to each cell, including
extent, dimension, boundary, activity bit, and extendable attributes. Conversely, in the present paper an oriented unit $p$-chain is characterized only by a signed integer index to $U_p$ basis, and by its signed $(p-1)$ boundary cycle, stored as a sparse column in $[\partial_{p}]$. Of course, all topological queries, both local and global, are allowed by suitable SpGEMM multiplication. We  allow for nesting inner cycles (holes) and substructures to the cells, but not for explicit containment of internal edges and points.  The most part of topological algorithms is algebraic in nature.

\cite{Zhou:2016:MAS:2897824.2925901}, by Zhou, Grinspun, Zorin, and Jacobson,  computes mesh arrangements for solid geometry, takes as input any number of triangle meshes, resolves triangle intersections in 3D, and assigns a winding number vector to subdivided cells, to evaluate variadic Boolean expressions.
Their data are represented by (small) BSPs enriched with explicit convex surface patches on nodes, and adjacency structure between nodes, together with a large amount of additional information.
``The crux of our method is construction of the mesh arrangement data structure, consisting of cells annotated with winding numbers, patches and their adjacency graph, that allows us to extract results of a variety of operations from the arrangement'' (Zhou et al., page 3).
Just consider their example where 10$K$ geometric models are intersected. It seems reasonable to assume that each of them weights for at least 1$K$ triangles, producing winding number vectors, a forest of BSP trees, incidence structures, \emph{one-to-one} with cells of the resulting arrangement. Telling how the models are distributed in space would give the reader a clue about  the amount of subdivided 3-cells, and its growth function.

Then, compare the approach~\cite{Zhou:2016:MAS:2897824.2925901}  with our topological method, where each unit $d$-chain (cell) is described just by the sparse array of signed indices of cells on its boundary cycle, and each input 2-cell produces its \emph{(local) fragmented} chain complex of maps, for a total number of 2D complexes equal to that of input 2-cells before 2D fragmentation.
The approach~\cite{Zhou:2016:MAS:2897824.2925901} to variadic Boolean expressions, while very bright, actually requires a scanning of the full stack of data to extract the result of every Boolean expression. Zhou, Grinspun, Zorin, and Jacobson claim that ``the method is variadic, operating on any number of input meshes.'' Actually, their approach can be split into two stages: first, \emph{adding iteratively meshes} to an arrangement;  second, executing all classic Boolean operations. Contrariwise, we do not add each input to the previous result but, in decomposition stage, operate \emph{independently} on each input 2-cell (see~\ref{quotient-set}), according to an embarrassingly parallel data-driven approach. It is also remarkable that the present approach works with more general meshes: sets of 2-manifolds  with- and/or without-boundary, sets of non-manifolds, sets of 3-manifolds, \emph{etc.}, versus just sets of triangle meshes. We do not discuss it here, but extending  our approach to Boolean operations and to Boolean functions is straightforward.

An extremely fast mesh repairing algorithm with guaranteed topology is described in
\cite{Attene:2014:DRS:2953208.2953514}, mostly based on floating point arithmetics, and requiring exact arithmetics only in relatively few situations. At variance with Attene~\cite{Attene:2014:DRS:2953208.2953514}, we not discriminate between manifold and non-manifold case, and do not use any special data structure in any steps of the pipeline, except 1D interval-trees and $kd$-trees for acceleration. By identifying the conditions that make floating-point arithmetics not reliable, J.R.~Shewchuk, the author of the Triangle library\footnote{\href{https://www.cs.cmu.edu/~quake/triangle.html}{https://www.cs.cmu.edu/~quake/triangle.html}}~\cite{shewchuk96b,SHEWCHUK200221} used in our method, identifies the key for fast robust geometric predicates in adaptive precision floating-point arithmetic~\cite{Shewchuk:1996:RAF:237218.237337}. He wrote that: ``the techniques Priest and I have developed are simple enough to be coded directly in numerical algorithms, avoiding function call overhead and conversion costs.''  In fact, the numerical results we obtained on triangulations with large arrangements of 2D line segments are very fast. We ported the CDT (Constrained Delaunay Triangulation) functions from his C library to Julia language\footnote{\href{https://github.com/cvdlab/Triangle.jl}{https://github.com/cvdlab/Triangle.jl}}, and used it to triangulate on-the-fly each non-convex 2-cell, in order to correctly compute the ordering of ``corollas'' 2-cells around ``pivot'' 1-cells in the 3D TGW.

In \cite{Campen:2010}, Campen and Kobbelt present a technique to implement operators that modify the topology of polygonal meshes at intersections and self-intersections, by combining an adaptive octree with nested binary space partitions. An analogous decompositive technique was introduced in the geometric language PLaSM~\cite{Paoluzzi:1995:GPP:212332.212349,Paoluzzi2003a}\footnote{\href{https://github.com/plasm-language/pyplasm}{https://github.com/plasm-language/pyplasm}} since 2004 by Scorzelli, Paoluzzi and Pascucci, in the contest of progressive geometry detailing allowing parallel modeling with BSP trees~\cite{Paoluzzi:2004:PDB:1217875.1217907,ScorzelliPP-PSM2008}. The technique is now being substituted by methods given in the present paper, since it does not guarantee sufficient robustness and speed.
In \cite{Guibas:98}, Guibas and Marimont describe a dynamic algorithm to compute the arrangement of a set of line segments in the \emph{digital plane}, and to snap the intersection points at the pixel centers. At variance with them, we snap small clusters of very close (numerically ``quasi-congruent'') points, rounding at the center of their $\epsilon$-neighboroughs in 2D, with average diameter of $10^{-16}$, close to the resolution of IEEE-754 binary floating point.

Barki, Guennebaud, and Foufou
claim that~\cite{BARKI20151235} presents  an exact, robust, and efficient method to execute regularized Booleans on general 3D meshes. They use a triangulation of all faces, and reduce the intersection of two surfaces to the 3D intersection of two triangles. Their simple decomposition process for intersecting faces is very similar to the old paper~\cite{Paoluzzi:1989:BAO:70248.70249} by Paoluzzi and his students.  Note that both~\cite{BARKI20151235} and~\cite{Paoluzzi:1989:BAO:70248.70249} contain a procedure for computation of regularized Boolean operations including isolated shells. The novel point about this matter here is that  in the present paper outer and inner oriented shells, as well as isolated components, are handled through signed sum of closed chains and implemented as sum or difference of their vector coordinates in $\mathbb{Z}$ or $\mathbb{Z}/3\mathbb{Z}$.

Finally, we recall that
\emph{Half-edge}, the smallest known efficient representation~\cite{Muller:78} of topology of  planar graphs and closed 2-manifolds by  Muller and Preparata, largely used in Computational Geometry for triangulations and Voronoi diagrams, as well as in meshes for games, requires $6\#\texttt{E}$ space. With the $[\partial_1]$ \emph{and} $[\delta_1]$ operators given in the present paper, we obtain for this class the optimal size $\Omega(n) = {4\#\texttt{E}}$, equal to the input size: two vertices and two faces per edge. It is well known~\cite{Woo:85} that both \texttt{EV} and \texttt{EF} relations weight for $2\#\texttt{E}$, i.e., equal to the space occupied by $[\partial_1]$ {or} $[\partial_2]$ maps, which also allow for the algebraic equivalent of multiple database queries at once---via product of the operator matrix times a matrix of unit binary columns, corresponding to single elementary queries..

\section{Past development and prospects of this project}\label{sec:projects}

This last section sketches birth and development of  contents discussed in this paper, provides references to past and present uses, outlines possible extensions, and introduces the prospects of research that are opened by our linear algebraic approach.

The first three authors started this project about computing with chains, cycles, cochains, and (co)boundary in a seminar series on novel algebraic methods for physical simulation and optimization of geometric design, during the sabbatical of V.S. in Rome [year 2000]. This project was awarded  the IBM SUR award in 2003. LAR sparse arrays, big geometric data and geometric services were discussed in many meetings in Rome, Paris, Madison, Berkeley, and Berlin. Our conversations produced some papers~\cite{DiCarlo:2009:DPU:1629255.1629273,ieee-tase,Dicarlo:2014:TNL:2543138.2543294} that started a lasting sequence of web and face-to-face discussions, algebraic experiments and software tests, producing in recent years three open-sourced partial implementations in Python and Julia. Partial implementations were used for software-based experiments of user-tracking and interior geo-mapping in LAR-based  Building Information Modeling (BIM), meta-design of a general hospital, and delivery of web services aiming at deconstruction and reuse of buildings~%
\cite{stag.20151290, visigrapp17:cvdlab, SpiniMDCP-WEB3D2016.bib, paoluzziMS:2014}.

Currently, some of the authors work to materialize a Julia package~\cite{DBLP:journals/corr/abs-1710-07819} for topological and geometric design, already including a preliminary implementation of the algorithms in this paper.
A~full implementation, covering the handling of isolated components and holes, is under development. A~second version will include vectorization on the GPU using native Julia~\cite{Besard:2017,besard:2017a}. Next, our plan is to port the chain-maps pipeline on Nvidia's GPX-1 (A.P.~has recently obtained the hardware in-house).
We are already using the modeling approach introduced here relying on the Julian open-source library, for rapid development of building models from analysis of Italian Cadastre  documents and 3D models from pictures or 3D images scanned by flying of drones.

We hope that the basic structures and algorithms discussed in this paper may also find some  appropriate use when combined with representations for convolutional neural
networks, based as well on tensors and linear algebra, in order to properly combine image understanding and geometric modeling.
In particular, our approach to compute the
chain complex\footnote{Of very general type, with basis cells possibly non convex and multiply connected.} of an unknown space arrangement
should match well with deep NNs~\cite{Goodfellow:2016:DL:3086952,Boxel:2016:DLT:3019358}.
Also our first experiments with topological methods in medical imaging~\cite{Paoluzzi-DCFJ2016,ClementiSSPP-CAD16} look promising.

\section*{Acknowledgements}
We would like to thank the anonymous reviewers who carefully read our manuscript and provided us with many useful comments. Alberto Paoluzzi and Antonio DiCarlo gratefully acknowledge the support of Antonio Bottaro CTO of R\&D at Sogei, now CEO at Geoweb, who believed in our work, and the EU project \href{https://www.medtrain3dmodsim.eu}{medtrain3dmodsim}. {Vadim Shapiro is supported in part by National Science Foundation grant CMMI-1344205 and National Institute of Standards and Technology.

\section{Summary of results and conclusion}\label{conclusion}

We have introduced a novel view on topology computation of space arrangements, that may find good use in disparate subdomains of geometric and visual computing, discussed an original computational architecture based on linear topological algebra, and claim that our approach is in tune with current trends towards hybrid hardware and its more advanced software applications. In particular, in this paper we provide a pseudocode implementation of the full computational pipeline: from a collection of virtual geometric objects to the chain complex $(C_p,\partial_p)$ of their partition of space, giving a full characterization of the topology induced by the input. This result is obtained going beyond simplicial complexes, and working with general piecewise-linear topology with non-contractible cells. A full porting of this approach to Julia, ``the fresh approach to technical computing''~\cite{BEKS14}, is ongoing (a preliminary open-source implementation is available); some parts are already parallelized and others are open to be. Among the strong points we cite: the compact representation; the combinable nature of maps, allowing for \emph{multiple} queries about the $3\times 3$ local topology relations, via fast sparse kernels for multiplication and transposition; the independent fragmentation of input cells through cell congruence; and the topological gift wrapping algorithm. Last but not least, the whole approach seems to be extendable to higher dimension.
We believe that a full real-time implementation of our algorithms on GPUs will generate new techniques for image understanding, in particular when inputs come from next-generation 3D cameras, already on the market and going to be installed on self-driving mobile vehicles. Part of this work was developed within the framework of the IEEE standardization of model extraction from medical images~\cite{10.1109/MC.2014.103}.