InShaDe: Invariant Shape Descriptors for visual 2D and 3D cellular and nuclear shape analysis and classiﬁcation

,


Introduction
The last decades have witnessed the rapid improvement and proliferation of high-throughput digital acquisition technology. As a result, high-quality digital representations of real-world * Corresponding authors: magus@hbku.edu.qa, jeschneider@hbku.edu.qa. scenes and objects have become commonplace in many application domains. In biology and medicine in particular, the rise of whole-slide scanners and the digitization of traditional, confocal, and electron microscopy has led to both fully digital analysis and the creation of large image databases [1]. Early uses of this technology mostly included tele-pathology, solicitation of second opinions, and education in research and clinical practice, as well as visual ultra-structural analysis in neu- roscience. In most cases, digital workflows were designed to closely mimic the traditional investigation process. Only recent years saw a shift towards exploiting the large amount of information in the acquired images and collections by means of novel data-driven analysis methods [2]. In this context, a wide array of basic tools are employed, ranging from handcrafted feature descriptors over fully data-driven approaches to combinations of matching of various approaches [3]. Machine learning, and, especially, deep learning approaches have become popular in the context of digital pathology and biology. Their significant success stems from their ability to provide automatic tools for tasks such as segmentation and labeling of cellular entities from large individual microscope images [4,5] or for supporting connectomics investigations by reconstructing the neural connections in large portions of brain tissue samples [6,7]. One drawback of the current purely data-driven frameworks, however, is that they often disregard specific domain knowledge and taxonomies [8]. As a result, tools for domain-specific proofreading of segmented images, classification according to taxonomies, filtering, visual exploration, and subsequent computation are still lacking. For this reason, many applications require the use of descriptors that, by design, preserve some domain-specific characteristics. To better exploit the capabilities of novel learning frameworks we thus advocate not to dismiss designing features a priori. Instead, we believe that it is necessary to research the development of methods integrating such designed features into powerful descriptive models. One benefit of taking a hybrid design+data-driven approach is that domain knowledge can be integrated into the design process. This potentially leads to models that are easier to explain, and, thus, result in increased discrimination performance for human analysis. In addition, since features are designed, the training efforts in terms of the amount of data and computational power required can be eased.
Inspired by these considerations, we propose a novel visual analysis pipeline based on discrete differential geometry concepts, whose recent findings provide very powerful theoretical formulations for describing 2D and 3D shapes [9,10]. The framework is based on a 2D/3D shape descriptor that can be used in various application domains to complement or enhance generic deep learning networks, such as U-Net [11]. The descriptor, dubbed InShaDe, is based on the concept of discrete curvature along closed, resampled contours (see Fig. 1 for an overview of the proposed pipeline). In 2D, discrete curvature is computed using vertex and edge osculating circles. Interleaving the resulting edge and vertex curvatures produces a highresolution curvature vector. These curvature vectors are naturally invariant under rigid body transformations (translations and rotations). Invariance under parametric shifts is ensured by using energy-based elliptic Fourier descriptors. Cellular shapes in 2D tissue samples are sliced, implying that their apparent size may be smaller than the real radius [12]. We therefore also propose to achieve optional invariance under uniform scaling by replacing components in the curvature vectors by standard (z-)scores. This manuscript is an extended version of the conference contribution [13] recently presented at the Eurographics Workshop on Visual Computing for Biology and Medicine (EG VCBM 2020) held in Tübingen. The conference paper presented: i) a robust geometry processing pipeline for computing 2D invariant shape descriptors exploiting shifted linear interpolation and discrete differential geometry schemes, and ii) visual mapping schemes from 2D cellular contours to shape descriptor embeddings based on modern dimension reduction schemes such as UMAP [14].
This manuscript extends the original pipeline by: iii) extending the shape descriptors to 3D shapes represented by triangle meshes, and iv) augmenting the shape descriptors by additional features based on sparse coding to improve analysis and classification performance.
In addition, we study the use of the proposed pipeline for proofreading and visual, unsupervised classification of various histology images. We also present results for various 2D and 3D data sets stemming from histology and neuroscience.

Related Work
Our work is concerned with shape feature extraction from closed contours and surfaces and with the analysis of histopathology images and electron microscopy stacks. These are very broad topics and a full coverage of the state-of-theart is beyond the scope of this paper. For a comprehensive overview of all related fields, we refer the reader to various surveys on 2D shape analysis [15,3], digital histopathology analysis [16,17,18], and recent deep learning methods for cellular analysis [19].
In the following, we discuss the methods that are most closely related to our approach.

Shape feature descriptors
During the last two decades, significant research efforts have been carried out on both the theoretical and the practical aspects of the shape-based image retrieval problem [3]. For an overview of the seminal methods for shape-based invariant feature extraction for object recognition, we refer also to Yang et al. [20].
In general, there are two main modeling strategies for representing shapes: region-based methods and boundary-based ones. Region-based techniques use moment descriptors to describe shapes, like geometrical moments [21], Zernike moments [22,23], Legendre moments [24], and Tchebichef moments [25]. Although region-based approaches are global in nature and can be applied to generic shapes, boundary-based techniques appear to be more efficient for handling objects that can be described by their object contours. In this latter category, a number of boundary-based techniques have been proposed, including Fourier descriptors [26], curvature scale space [27], and wavelet descriptors [28]. Our descriptor combines the features of curvature analysis and Fourier analysis, similarly to the technique proposed by El Ghazal et al. [29,30]. Differently from them, our method is based on recent findings in discrete differential geometry [31], thus resulting in a more robust formulation with respect to the sampling strategy and better classification results. Moreover, Osjanikov et al. applied the concepts of invariant features in the 3D world to the problem of non-rigid shape search and retrieval in large databases [32].
Complementing advances for the general classification problem, machine learning strategies exploiting the existence of large amounts of data have led to significant advances [33,34]. Many current efforts attempt to work directly on raw data images [35,36], by designing deep neural networks in which the modeling is hidden in the network design and training strategy and the feature computation and filtering of information is automatically performed by the network. At the same time and in an attempt to simplify classification and automatic shape generation, techniques to reduce the depth of networks by introducing meaningful parameterizations or embeddings of input shapes are gaining interest, since such parameterizations can simplify the automatic classification or shape generation (model-based or "shallow" learning) and reduce the number of training examples [37]. Our work goes towards that direction, since we propose a simplified contour description that can be used either for supporting machine learning frameworks or for supervised visual analysis. In this work, we focus on the latter aspect.

Histology analysis
Digital pathology and microscopy-image analysis is widely used in the biomedical domain for comprehensive studies of cell morphology or tissue structure. In most cases, analysis is carried out through manual assessment, which is labor-intensive and prone to inter-observer variations. Computer-aided systems have recently attracted significant interest since they can dramatically reduce the manual efforts and increase reproducibility [38,39,17].
Among the various parts composing a computer-aided diagnostic system, nucleus or cell detection and segmentation play a key role to describe the molecular and morphological information underlying the investigated samples [17,40]. In the past few decades, many efforts have been devoted to automated nucleus/cell detection and segmentation, and an independent field named computational pathology emerged simultaneously to the rapid proliferation of deep learning (DL) models for quantitative analysis of spatial patterns in digitized whole-slide images (WSIs) of cancerous tissue [41]. To this end, various techniques for the detection, extraction, recognition of pathological patterns at various scales have been recently established [42,43,44]. Various medical studies have since demonstrated the potential of DL models in detecting neoplastic tissue and recognizing diagnostically relevant structures [45,44].
One of the most successful and widely used architecture is U-Net, introduced in 2015 by Ronneberger et al. [11]. U-Nets operate on the entire image and jointly segment and provide per-pixel labels, leading to an improvement in spatial segment and label coherence. The same authors also demonstrate that U-Nets improve accuracy on several bio-image segmentation tasks, even when the data set is relatively small [11]. In the context of nuclear segmentation of histopathology images, Chidester et al. [4] enhance U-Nets by enforcing rotationequivariance to groups, similar in style to group-equivariant CNNs (GCNNs) [46]. Moreover, in order to attract efforts to particular tasks in medical imaging, various challenge contests and public data sets have been published [47,41]. However, as of this writing, such methods are still far from being accepted in fully automated clinical workflows [48]. Proofreading efforts from domain scientists are, thus, still required to double-check labeling consistency and segmentation accuracy [8].
Consequently, the work presented herein provides a visual analysis framework that supports digital histologists to efficiently carry out investigations on labeling and segmentation quality. Our input data is the automatic segmentation obtained from networks of the U-Net family [4]. The proposed framework then allows for visual analysis in a reduced parameter space obtained by performing dimension reduction on our Fourier-based contour shape descriptor.
To properly capture the visual variance of nuclear shapes under dimension reduction, autoencoders [49] provide a convenient way to effectively uncover latent feature spaces. It is thus no surprise that their use is increasing in popularity [50,51]. For instance, Xu et al. propose Stacked Sparse Autoencoders [50] to learn high-level features from pixel intensities. They are then applied to high resolution breast cancer histopathology images. Hou et al. [51] modify the general autoencoder scheme by applying adaptive convolutional filters to match the size of the nuclei to be represented. In this work, we use k-sparse autoencoders [52] to produce feature vectors that describe the inner visual features of nuclei. We then use these feature vectors to augment our proposed shape descriptor, resulting in a description of the exterior (shape descriptor) and interior (auto-encoder) of each nucleus.

Shape analysis in neuroscience
Recent advances in imaging technology have led to the availability of 3D sparse and dense reconstructions of brain cells at high resolution. This, in turn, has fueled the development of various methods for shape analysis in the context of automatic classification to aid studying the variability associated with different structures and conditions [53,54]. Likewise, the availability of high-resolution imaging data has also triggered shape analysis studies of brain structures at the nanometer scale [55,56,54]. For instance, Queisser et al. [57] propose a method to reconstruct the 3D view of cell nuclear envelopes from laser scanning confocal microscopy data. Wittmann et al. [58] later use this method to show how synaptic activity induces significant modifications in the geometry of the cell nucleus. To study heterogeneities in nuclear shapes obtained through optical projection tomographic microscopy, Nandukumar et al. [59] use conformal mapping to extract rotationinvariant shape descriptors. Finally, Agus et al. [60,61] perform classification of nuclear brain cells through implicit and explicit shape representations of cell nuclei obtained from electronicimaging data. They demonstrate an improvement in terms of classification accuracy over previous approaches based on simple spherical or ellipsoidal fittings.
In this paper, we improve existing work [61] by introducing feature vectors based on the spherical harmonics spectrum of mean curvatures. This reduces the amount of data serving as the descriptor from three complex vectors to a single real-valued vector, resulting in faster implementation and processing times. Fig. 1 schematically summarizes the InShaDe framework. As can be seen, we use two separate yet similar pipelines for 2D and 3D closed nuclear envelopes that use different microscope imaging techniques as data source. The two methods have important similarities:

Methodology Overview
• they use the curvature signal (planar curvature in the 2D case, mean curvature in the 3D case), • they involve parametrization (circular parametrization in the 2D case, and spherical parametrization in the 3D case), • they involve Fourier analysis (Elliptic Fourier Analysis in the 2D case, and Spherical Harmonics decomposition in the 3D case), It is worth noticing that the Spherical Harmonics framework is a 3D generalization of the Elliptic Fourier Analysis [62]., • they use the same strategy for computing energy descriptors according to the harmonic frequencies [63].
Overview of InShaDe 2D. The input to the InShaDe 2D pipeline are segmented nuclear envelopes of cells obtained by applying an U-Net [4] on microscopic histopathology images (see Fig. 1 top). We then extract closed contours (Sec. 4.1) from each segmentation mask and perform the following processing steps (see also Contour smoothing serves to reduce pixelation noise, whereas geodesically uniform resampling removes sampling biases and is a pre-requisite to computing discrete curvatures using discrete differential geometry formulations using osculating circles. Embedding the resulting descriptors in constant dimensions helps in removing noise and spurious frequencies during the EFA stage, but is also necessary to allow for easy comparison between shapes using, e.g., cosine or Euclidean metrics. The Fourier analysis is used to remove shift (i.e., choice of origin) from the parameterization of the closed curve. So far the resulting descriptor is invariant under translation and rotation (3) and invariant under parameterization shift (6). In addition, the optional feature scaling step (4) ensures invariance under uniform scaling. In the result section, we furthermore show how the final descriptor can be used in combination with dimension reduction schemes for visualizing clusters of nuclear shapes with similar geometric characteristics.
Overview of InShaDe 3D. The input of the InShaDe 3D pipeline are closed triangular meshes extracted from image stacks obtained through Serial Section Electron Microscopy acquisition of samples from rodent brains [54]. These shapes represent the envelopes of brain cell nuclei and are obtained from images through a processing pipeline involving automatic segmentation tools as well as manual proof-reading tools [64]. We then process the 3D meshes by performing the following operations (also see Fig. 1, bottom): (1) discrete mean curvature computation (Sec. 5.1), (2) spherical parameterization using Willmore flow (Sec. 5.2), We use discrete mean curvature (1) as the basis of our embedding and to represent the features of 3D shapes. In contrast to previous formulations of shape decomposition [60,61], the proposed embedding is based on a single scalar-and real-valued

InShaDe 2D
In this section, we provide details for the various processing steps for computing the descriptor for closed shapes extracted from 2D images.

Contour Extraction & Chordal Parameterization
Given a segmentation mask, we extract a closed contour enveloping each nucleus using iso-contouring (specifically Marching Squares, which is a special case of the Marching Cubes algorithm [65]). We reject open contours (i.e., the nucleus intersects the image boundary) and contours falling into the lowest 5% with respect to their number of samples. Let C := {p i } N 1 , a closed curve with N vertices p i . We let ∆ i := p i+1 − p i , the i th edge, consistent with Bobenko [31], and abbreviate l i := ∆ i 2 (edge length). We then obtain an initial chordal parameterization t (C) with t 1 := t (p 1 ) = 0 and

Contour Smoothing
The discrete nature of binary segmentation masks may lead to pixelation artifacts in the extracted contour. To prevent the resulting high spikes in curvature, we pre-smooth contours iteratively, using a superscript (k) to denote quantities at iteration k. The process is shown in Fig. 2. Specifically, we replace each vertex with a length-weighted average of the bisector of adjacent edges, (1) As shown by Gottschalk [66], this sum of of length-weighted edge bisectors computes the barycenter of the points on the piecewise linear curve segment , it is numerically stable and robust. Similar to virtually all smoothing operators, this does not yet preserve area. We therefore compute the area a (0) enclosed by the curve prior to smoothing and the area a (k) after each iteration. We then scale the curve by

Geodesically Uniform Resampling
In order to remove sampling bias and to employ discrete differential geometry formulations for vertex and edge curvature, we perform geodesically uniform resampling. We do so by placing equidistant samples p on the piece-wise linear curve C, thereby yielding a new piece-wise linear curve C that is Arclength parameterized with respect to a unit scale u. Starting at a point p 1 = p 1 and u = 1, we intersect the edges of C with a unit circle around p 1 . This yields between zero and two intersection. If we find two intersections, we select one intersection as p 2 and keep track of the last edge, ∆ 1 = p 2 − p 1 . We then continue intersecting linear segments with unit spheres, but when deciding on p i , we chose the intersection that maximizes This enforces progress along the curve and prevents jumping back and forth on the curve. For our data, we did not encounter the case of finding less than two circle-curve intersections. A total of zero intersections would correspond to extremely small contours that cover less than a few pixels after processing; and we remove the bottom 5% shortest curves. One crossing would arise if part of the contour degenerates into a double line segment; Marching Squares does not extract such pathological curves.
Once the best intersection p N "laps" past p 1 , we use p N = p 1 instead to close the loop. This means that the last edge ∆ N−1 may be shorter than unit length. In order to resolve this issue, we now calculate the length L of the curve. Knowing that ∆ i 2 = u for all but the last edge, we have L = N − 2 u + ∆ N 2 . To obtain an u for which u −1 L is approximately integral, we round u −1 L to the nearest integer L and update u ← L −1 L.
We then revert to placing samples along the original curve C with the updated spacing u. We repeat this process until the rounding error ρ = |u −1 L − L | (using the old u and the updated L ) becomes negligibly small. While we do not have a proof of convergence of this heuristic at the time, we note that the rest of our method is orthogonal to this Arc-length parameterization. This means that, in the future, as more robust methods become available, this step can be exchanged. In all of our experiments, three to five iterations reduced ρ to less than 10 −4 . Given any number x ∈ R + , rounding to the nearest integer changes x by 0.25 on average. We therefore expect that |1 − u| ≈ 0.25L −1 , which we see confirmed in our experiments with typical contour lengths of more than 100 pixel widths (for reference, L = 100 → |1 − u| ≈ 2.5 × 10 −3 ). The result of this step is a new piece-wise linear curve C that is Arc-length parameterized with respect to a close-to-unit scale u.

Discrete Curvatures
For a discrete Arc-length parameterized curve, there are two definitions of discrete curvature based on osculating circles [31] (Sec. 2.3 therein). By defining the turning angle at vertex p i as and by embedding the planar curve in the z = 0 plane (see also Fig. 3), we obtain, assuming for now an Arc-length parameterization with ∆ i = 1 for all i, the (unsigned) vertex curvature: For the edge curvature we use the standard equation [31]: The choice to use unsigned vertex curvature was made to be consistent with the unsigned edge curvature. Using such a discrete differential geometry approach results in much more robust and stable curvature estimates than by using an intermediate interpolating spline. A reason may be that splines tend to over-and undershoot near vertices, and are thus not representative of the curvature in these points. Since one of our goals for the final shape descriptor is optional scale-invariance, we still have to scale curvatures back from our arbitrary unit length u to u = 1 in case scale-invariance is not desired. This is achieved by dividing each κ v and κ e by u 2 . Finally, we interleave vertex and edge curvatures to obtain a high-resolution, coherent descriptor. After this step, we also abandon the notion of curvature "living" on vertices and edges and transition to the notion that the shape descriptor computed so far is a vector in a highdimensional vector space. We also adopt the notion that this vector represents a 1D periodic signal on a uniform grid on the 2D circle. This interpretation is crucially supported by the fact that all edges have the same length prior to computing curvature. The descriptor computed so far is invariant under translation and rotation, but neither parametric shift nor scale. We now establish the optional scale-invariance followed by shiftinvariance.

Feature Scaling
Given a sequence of curvatures, {κ i } 2 N i=1 , we compute standard scores (also called z-scores) by mapping Fig. 4: Resampling to constant dimension: to reduce noise and spurious frequencies during the Fourier analysis and to enforce constant dimensionality of our descriptor, we apply uniform resampling through shifted linear interpolation [67]. In this example, we show resampling with 64, 32 and 16 points respectively. where are the empiric mean and variance, respectively. Such a scaling is commonly employed in statistics as well as in training convolutional neural networks. However, normally standard scores are computed using global moments derived from the entire data set. This, in turn, does not provide full scale-invariance, since vectors with pre-dominantly small components will stay small. In contrast, by computing individual standard scores we enforce the optional scale-invariance of our descriptor. Assuming that the curvature components of each vector are normaldistributed results in the expectation that all but 0.2% of the data is represented by z-scores in the range [−3, 3].

Constant Dimensionality
Resampling the contour to a constant dimensionality as depicted in Fig. 4 allows us to control the number of elliptic harmonics in our Elliptic Fourier Analysis in a way to agree with the Nyquist sampling constraint. It is also a pre-requisite for easy comparison of shape descriptors using, e.g., cosine and Euclidean metrics. As an added side-benefit, it also allows us to eliminate remaining traces of noise on the curve. In this paper, we perform this resampling step based on shifted-linear interpolation [67] for the following reasons: (i) shifted linear interpolation achieves performances that compare favourably to cubic interpolation at a much lower computational cost, (ii) shiftedlinear interpolation is still convex, albeit with respect to shifted samples. It is thus free of oscillations and the amount of foreign frequencies introduced by resampling can be computed easily.
The basic idea of shifted-linear interpolation is to sample the original signal at positions other than the original underlying sampling grid, followed by standard linear interpolation. Blu et al. prove, somewhat surprisingly, that there is a data-independent and thus constant shift τ ≈ 0.21 that results in L 2 −optimal reconstruction of the unknown original signal given only the known samples [67]. Samples κ i at shifted positions t i are obtained using the infinite impulse response scheme described by Blu et al., It should be noted, however, that linear interpolation on κ is literally shifted "to the right" by τ, meaning that a sample κ (t) corresponds to κ(t − τ). The resulting interpolation thus becomes a shifted discrete convolution of the hat kernel with the shifted discrete signal κ :

Elliptic Fourier Analysis (EFA)
To achieve shift-invariance (i.e., invariance under choice of parametric origin), we consider the Fourier spectrum of each given curve. In particular, we compute elliptic Fourier descriptors [62], similarly to what was proposed by Khazhdan et al. [63] and what has been successfully used in various applications [60,61,12].
For a piecewise linear, periodic function κ(t)t ∈ [0, 2π] representing the curvature of a contour, its Fourier elliptic expansion is obtained through linear combination of elliptic harmonics functions which provide a complete orthonormal basis for the decomposition In order to compute the coefficients for the curvature function κ (t) representing closed contours, we normalize the parameterization t to the interval [0, 2π]. As we are concerned with closed contours, the assumption of periodicity, t = 0 ≡ 2π is naturally supported. We then consider the classic method proposed by Kuhl and Giardina [68] for piecewise linear contours. This method essentially equates the discrete time derivative of Eqn. (11), at locations p i , with a Fourier expansion of the time derivative of the curvature, Noting that in Eqn. (13), the coefficients α n and β n can be computed as Kuhl and Giardina derive the following for the n th harmonic, by equating the two different derivative expressions in Eqn. (13) and (12): We would like to remind here that, according to the Nyquist theorem, the number N s of contour samples after smoothing and resampling limits the number N h of harmonics necessary to reconstruct the contour curvature without adding noise N h ≤ N s 2 . Finally, in order to obtain shift-invariance, we compute harmonic energies through the Euclidean norm of the harmonic coefficients [63], resulting in the following Curvature Fourier Descriptor K with which provides a vector of shape features that can be used for various machine learning applications. Like the more commonly employed traditional Fourier transform, the elliptic Fourier transform results in a space-agnostic spectrum, thereby making our descriptor invariant under parameter shift (translation of the underlying domain). In this paper, we chose the elliptic Fourier transform over the traditional Fourier transform since its additional expressiveness resulted in better results. Fig. 5 demonstrates both rotation-and shift-invariance.
Sorted curvatures. We also consider another, much simpler scheme for obtaining shift-invariance, namely, to sort the individual (unsigned) curvatures from highest to lowest (see also   6). The feature vector obtained in this way can be further augmented with energy-based coefficients to obtain a composite feature vector. We would like to note that, in our experience, sorted descriptors are outperformed by spectral descriptors if used in stand-alone fashion. Sparse-coding based image descriptors. Finally, we consider a compact descriptor of images representing cellular structures obtained through sparse coding. To provide additional shape cues for this case, we place the segmented nuclei on a black background, as shown in Fig. 7, left. Sparse coding methods are typically composed of two steps. Firstly, an offline learning process for finding a dictionary W that sparsely represents the image data {I i } N i=0 , and, secondly, an encoding step that maps a given input image I to a compressed feature vectorx using W , normally through a pursuit algorithm for minimizing the constrained least squares problem For obtaining the codebook (see also Fig. 7, center) and creating the approximated sparse representation of nuclei images (Fig. 7, right), we use the K-sparse autoencoder proposed by Makhzani and Frey [52]. The technique uses linear activation functions and tied weights. In contrast to other autoencoders, only the k largest codes are used while the others are set to zero. The resulting code k-vector gives us additional cues that, albeit not rotation invariant, can be combined with the InShaDe descriptor (see also Fig. 6, right). Section 6 evaluates various descriptors obtained by composing the three different feature vectors: sorted curvatures, energy coefficient and sparse coding weights. Our 3D pipeline is a natural adaptation of the 2D case. From a mathematical perspective, the Laplace-Beltrami operator on either 1-or 2-manifold induces a Fourier space. That is, the eigenfunctions of the Laplace-Beltrami operator constitute the "classic" Fourier space in the 2D case and Spherical Harmonics for the 3D spherical case. Albeit we use the elliptic Fourier transform for its superior performance in the 2D case [12], the two pipelines share the same mathematical foundations. We then compute discrete differential geometric attributes of the manifold and express them in this Fourier space. It is imaginable to generalize this even further by utilizing manifold harmonics [69], at the likely expense of higher computational complexity. The bottom half of Fig. 1 depicts the pipeline schematically for easy comparison with the 2D pipeline (Fig. 1, top). The details of the processing steps are provided in this section.

Mean curvature
According to the differential geometry theory of surfaces, for every twice-differentiable surface we can find the tangent plane for a point on the surface. We can then proceed to define a quadratic form, that is a polynomial containing only terms of degree two, using the two tangent directions x, y [70, chapter 19 therein]. This quadratic form, sometimes called the shape tensor describes extrinsic invariants of the surface, such as principal curvatures, at the point where manifold and tangent plane touch. This form is called the second fundamental form II, The second fundamental form approximates the surface z = f (x, y) with z = 0 the plane tangent to the surface (informally, z is the "height over tangent plane") in a neighborhood around the touching point. Therefore, the idea of the second fundamental form is to measure, in R 3 , how a surface curves away from its tangent plane at a given point. The eigenvectors of the 2 × 2 matrix II are called principal directions, and the eigenvalues are called principal curvatures, denoted κ 1 , κ 2 [9]. Given the principal curvatures, the mean curvature H = κ 1 +κ 2 2 provides a meaningful and natural description of 3D surfaces, and it can be computed on triangular mesh using a discretization of the Laplace-Beltrami operator [9]. In this work, we compute the mean curvature at each vertex of a closed triangle mesh. Given a parameterization (s,t) of the triangle mesh, we thus obtain a discrete, scalar function H(s,t) (see also Fig. 8).

Spherical parametrization
As hinted at in the last section, our 3D pipeline relies on a surface parameterization. We use a spherical parameterization, which, despite its apparent simplicity, is much harder to obtain than circular parameterizations of planar shapes. From a geometrical point of view, the second fundamental form can be used for classifying surface points according to the signs and values of principal curvatures κ 1 and κ 2 . Of particular interest are so-called umbilic ("locally spherical") points for which κ 1 = κ 2 . A measure for "local sphericity" can thus be defined based on κ 1 , κ 2 , such as the Willmore energy of a surface S, The geometric flow associated with this energy, will evolve any genus-0 surface S to a sphere, providing a way to obtain a spherical parameterization. In this work, we use the discrete Willmore flow formulation proposed by Crane et al. [71] that was also previously used by Agus et al. [61]. This spherical parameterization maps each vertex of the 3D input shape to a corresponding point (θ , φ ) on the unit sphere S 2 , thereby also providing a parameterization for the discrete scalar function H (θ , φ ).

Spherical Harmonics decomposition
The Spherical harmonic basis provides a Fourier basis for functions defined over a sphere. We can thus approximate a generic function defined over a closed surface as a finite linear combination of spherical harmonics Y m l (θ , φ ) up to a given maximum frequency L: where the weights w m l can be found through least-square error minimization with respect to the samples computed on the original 3D shape. To this end, we used a method similar to [61], with the main difference that, since the mean curvature signal is scalar, we only consider the real part of the spherical harmonics. As a result, the weight coefficients are also real. Specifically, given the real part of spherical harmonics R m l (θ , φ ) = ℜ(Y m l (θ , φ )), a spherical parameterization of the surface S, Σ S = (θ i , φ i ) = (θ (v i ), φ (v i )) ∈ S 2 , ∀v i ∈ S , and the mean curvature values computed across the surface H S = {H i := H(p i ), ∀p i ∈ S}, the spherical harmonic decomposition is obtained by computing the coefficients w = {w m l , 0 ≤ l ≤ L, −l ≤ m ≤ l} that minimize the square error: leading to the linear system We solve this system using LDL T factorization, a robust, symmetric pivoting variant of the Cholesky decomposition [72, chapter 4.1.7 therein], in combination with Tikhonov regularization [61]. This regularization add a diagonal matrix T depending on the Spherical Harmonics order and weighted by a small numeric value ν (see [61] for details). Hence, the final linear system has the following form.
We set the Tikhonov regularization weight to ν = 10 −5 in all experiments reported in this paper.

Energy coefficients
Similarly to the 2D formulation we use the harmonic energies, to gain rotation-invariance. Energies are defined as the Euclidean norm of the Spherical Harmonic coefficients w for each harmonic frequency separately. Specifically, the 3D curvature Fourier descriptor Φ is defined by These coefficients provide a compact descriptor of genus-0 shapes, and can be used for the analysis of nuclear envelopes extracted from Serial Section Electron Microscopy stacks.

Results
We implemented our general shape descriptor pipelines and tested them on several challenging use cases. In this section, we first provide details on our implementation (Sec. 6.1) and then provide an evaluation on general shape analysis, on the analysis of histopathological images, and nuclear shapes extracted from Serial Section Electron Microscopy (SSEM) stacks. We separate the evaluation of the two pipelines: for the 2D framework, we report on consistency evaluation performed on classic shape collections commonly used in literature for testing shape retrieval methods (Sec. 6.2). We also provide results obtained with our pipeline on various histology samples for medical diagnostics and neuroscience investigations (Sec. 6.3). For the 3D framework, we report on the usage of the pipeline for the classification of neural cells reconstructed from a rodent brain sample (Sec. 6.4). In both evaluations, we involve expert domain scientists, for providing a qualitative evaluation of the framework, and for getting suggestions for designing a full visual analytics framework for histology images.

Implementation notes
The code used to generate the results presented in this paper for the InShaDe 2D pipeline is available in GitHub 1 (Python scripts & Jupyter notebooks). After further testing and cleaning, we plan to also release the C++ code for the InShaDe 3D pipeline.
InShaDe 2D. We implemented the 2D geometry processing pipeline in Python using the following building blocks & modules: G-U-Net [4] for automatic segmentation, sklearn, skimage for contour processing and dimension reduction, interactive matplotlib for visualization. For testing the pipeline, we developed simple interactive widgets in which users can compare the clustering visualization in the parameter space to the reconstructed cellular shapes in the histology images. In order to attenuate the amplitudes of high frequencies, we used a frequency equalization scheme weighting of the InShaDe coefficients according to the square root of their order w(k) = √ k . Our geometry processing pipeline can be used in combination with different dimension reduction schemes and clustering methods. In this work we use the recent Uniform Manifold Approximation and Projection (UMAP) method, which is based on Riemannian geometry and algebraic topology [14]. For clustering, we used HDBSCAN [73] or k-Means [74,75,76] (depending on the case).
InShaDe 3D. We implemented the InShaDe 3D pipeline in C++, using conformal curvature flow based on spin transformations 2 . We furthermore used the Eigen library [77] for the LDL T solver arising in the SPH least squares optimization. For the SPH functions we used boost and libigl for geometry processing. We then fed the parameters derived from the SPH decomposition to standard machine learning methods (e.g., support vector machines, SVMs) using Python's sklearn to classify the nuclei.

Consistency validation
We first performed a consistency validation of the InShaDe descriptor. For this, we used MPEG-7 and Animals [78], which are among the most popular data sets for evaluating and comparing the accuracy of shape retrieval methods [78]. The MPEG-7 shape collection is composed of 1,400 binary images containing objects of 70 different classes [79] (see Fig. 9 left), while the Animals shape collection (see Fig. 9 right) is an even more challenging data set containing 2,000 binary images grouped in 20 classes of 100 animals each one. To test 1 https://github.com/HBKUVisCommunity/inshade/ 2 https://github.com/nitronoid/flo the InShaDe 2D descriptor we considered three different assessment criteria: (1) Retrieval accuracy of a basic SVM scheme, trained on an augmented data set. We triple the size of the input data set by adding randomly rotated and shifted copies of original images. Moreover, we use a hyperparameter optimization scheme to find the best SVM linear parameters with respect to cross-correlation accuracy, and we test the obtained model over the original collection. We also show the accuracy of shape retrieval in the form of a confusion matrix to highlight the accuracy differences between classes.
(2) Bull's Eye accuracy, which is commonly used to score shape retrieval tasks when the number of objects is limited: First, a similarity distance between objects represented by feature vectors is defined. For this, we use the L 1 norm (3) Qualitative visual assessment of the reduced parameter space obtained by projecting the feature vectors on a 2D plane through dimension reduction techniques and observing how objects cluster together. In our experiments we used UMAP [14].
For the composition of the feature vector we consider three different contributions: the sorted local curvature signal Γ, the elliptic Fourier analysis energy coefficients Φ, and the weight values of k-sparse autoencoding [52] (see the examples in Fig. 10). The purpose of using sparse coding features is not to precisely reconstruct the original image, but rather to extract the important information using very few parameters. The reconstructed images, using only 256 coefficients, clearly represent the main contents of the original images.
We investigate all possible composition permutations: sorted local curvatures alone (Γ), Fourier energy coefficients alone (Φ), sparse coding weights alone (Ω), local curvatures and energy coefficients (ΓΦ ), local curvatures and sparse coding  weights (ΓΩ), energy coefficients and sparse coding weights (ΦΩ), local curvatures together with energy coefficients and sparse coding weights (ΓΦΩ). For merging heterogeneous feature vectors, we perform pre-normalization of the various feature vectors. Fig. 11 shows the Bull's Eye score obtained for both shape collections on top of four retrievals for various feature vector compositions at varying number of coefficients. The highest accuracy was obtained for 240 coefficients and with the descriptor ΓΦΩ (composing sorted curvatures, energy coefficients, and sparse coefficients). Considering the feature coefficients alone, sparse coding descriptors outperform sorted curvatures and energy descriptors: we suspect that this is due to the fact that sparse coding is able to provide a valid description of both the boundary shape and the inner part of the objects. Nonetheless, incorporating the energy coefficients proved to be beneficial since the composed descriptors can take into account transformations like rotations and shifts. The obtained values are in line with current state of the art methods (see tables in [78]): for example the obtained Bull's Eye score on MPEG-7 data set is 0. 87 [80]). Better performances can be obtained through post-processing retrieval schemes that are orthogonal to our method and can be incorporated successively: for example the Online to Offline O2O scheme [78] applied on top of different descriptors can achieve Bull's Eye score up to 0.99 for the MPEG-7 data set and 0.66 for the Animals data set. In our experiments, we also noticed a slight degradation of accuracy for feature dimensions higher than 500, indicating that the curse of dimensionality can affect the proposed descriptor. Fig. 12: Processing time for descriptors: we compared the CPU processing time of InShaDe 2D energy descriptor with the sparse coding scheme computed through k-sparse autoencoders [52] with respect to the number of coefficients. We report the processing time per shape for MPEG-7 (on the left), and Animals (on the right).
We also compared the CPU computation times per shape between the various descriptors on the MPEG-7 dataset (the average number of vertices per shape is 1917) and Animals (the average number of vertices per shape is 961). For timing measurements, we used a workstation equipped with an Intel i9-9900 CPU (8 cores, 3.1 Ghz), 64 GB of RAM, and an Nvidia RTX 2080 GPU with 8GB RAM. In Fig. 12, we report the processing times per shape as function of number of coefficients, and for both datasets (left: MPEG-7, right: Animals). For the k-sparse autoencoder, we report the training time for 500 epochs using the total images of the dataset for training but averaging the times reported by the number of images. The input image resolution is 256 × 256 for MPEG-7 and 640 × 432 for Animals. It is worth noticing that, according to the number of coefficients, the memory resources needed for using k-sparse autoencoder are proportional to the size of input images and the size of output sparse descriptor and can easily reach the limits of available RAM in many systems. Moreover, the InShaDe pipeline can be further accelerated through GPU-friendly implementations that would be able to manage batches of images in parallel. Fig. 13 and 14 show the UMAP projection of the composed descriptor for the various shapes for both MPEG-7 and Animals data sets as a visual reference (separated in groups of at most 12 labels to reduce clutter).
It appears evident that the proposed descriptor is, in most cases, able to discriminate the shapes of MPEG-7 data set. The clusters for the Animals data set appear more confused, thus confirming the retrieval rates in this work and prior literature. Fig. 15 depicts a typical failure case of our scheme on the Animals data set. A leopard is considered very similar to a cow, a cat and another cow. The middle row shows the local curvature Fig. 13: MPEG-7 UMAP clustering: we test UMAP dimension reduction on the composed descriptor ΓΦΩ for MPEG-7 collections (256 feature elements per descriptor for a total of 768 features). The Bull's Eye score obtained for this collection was 0.898. In order to reduce visual cluttering, the various shapes are separated in 6 groups of maximum 12 labels.
Fig. 14: Animals UMAP clustering: we test UMAP dimension reduction on the composed descriptor ΓΦΩ for Animals collection. The Bull's Eye score obtained for this collection was 0.543. In order to reduce visual cluttering, the various shapes are separated in 3 groups of maximum 7 labels. signals computed over the shape contour, while the bottom row shows the composed feature vector containing the sorted curvatures (left), the energy coefficients (middle), and the sparse coefficients (right). Fig. 16 shows the confusion matrix for linear SVM classification obtained for our composite descriptor ΓΦΩ. In terms of accuracy, we obtain results aligned with state of the art methods (88.4% for our descriptor on the MPEG-7 data set versus 66% for Curvature based Fourier descriptor [29], 78% for blurred shape models [82], 78% Morphological Pattern Spectrum [83] and 90% for Zernicke moments with geometric features [79]). It is important here to note that the composition of shape features and image features significantly improved retrieval accuracy in agreement with prior work [79] (88.4% for the composed descriptor versus 87.5% for the sparse coding descriptor and 78% of the simple energy shape descriptor proposed in the conference paper [13]). Given that the results show the proposed descriptor to be consistent for classifying shapes of natural objects, we will now proceed to analyze its performance for the analysis of biomedical images.

Histopathology analysis
For the analysis of histopathology images, we use public domain data from the MoNuSeg contest [47], and the very recent PanNuke data set [41]. The former contains 30 images from seven organs with unclassified annotations of more than 20k individual nuclei. The latter contains more than 220K labeled nuclei from 19 different tissue types and, as of writing, is the largest open pan-cancer histology data set for nuclei instance segmentation and classification. Finally, we apply the pipeline to a whole slide image representing a paediatric appendix specimen. Fig. 15: Bull's Eye testing: Here, we show a typical failure case using the L1 Bull's Eye score for the descriptor ΓΦω on the Animals data set. Top row, left to right: a leopard is considered similar to a cow, a cat, and another cow. We also show the curvature signal of each curve (middle row) and composite descriptor (bottom row). Fig. 17 shows examples of images from the MoNuSeg data set [47] classified by our framework: we reduce the dimensions of the feature descriptors using UMAP, followed by k-Means for clustering. We then color-code contours by cluster. We notice that cells recognized as having similar shape features according to UMAP do not necessarily relate to intuitive discrimination through visible attributes such as length or thickness or smoothness. However, they do not only form feature clusters (same color) but also tend to form spatial clusters. The latter fact can provide additional visual information to digital pathologists for diagnosis through spatial aggregation of such clusters. While further investigation is needed to understand and evaluate the clinical value and to find explainable taxonomies, initial feedback from pathologists confirmed that in many cases nuclear features and clusters can provide decisive information for recognizing specific conditions. Histopathology shallow classification. We also tested whether our descriptor could be used for shallow classification of nuclear cells for diagnostic purposes. To this end, we tested various composed descriptors (Φ, Ω, ΦΩ, and ΓΦΩ) with varying number of features (64,128,256), and we trained a linear SVM classifier on PanNuke data set [41] for discriminating between three classes of nuclei: neoplastic cells, inflammatory cells, and others. Fig. 18, left, shows the accuracy performance of the various descriptors considered in this work. The highest accuracy (0.601) is obtained for the composite descriptor ΓΦΩ with 256 features per component. However, the improvement of the composed descriptor with respect to the sparse coding descriptor is almost imperceptible (0.578 for the descriptor Ω). In this case, sparse coding captures not just the shape but also texture information, which might be helpful for classification. Thus, the improvement brought by InShaDe is less pronounced compared to Fig. 11, where the input images only contain shape and not texture. In this case, the proposed shape descriptor cannot adequately discriminate the various cell classes accord-ing to the proposed taxonomy. It is still far from being sufficiently accurate for reliable classification of individual nuclei. However, this preliminary accuracy performance was obtained with a simple SVM classifier, and it can be improved by considering more sophisticated classifiers, like ANNs. In general, the dimension reduction plots show that cells of same type do not cluster together when using the InShaDe descriptor (see in Fig. 18 right some examples). Nonetheless, the presence of outliers in the parameter space can provide pathologists visual hints for proofreading the labeling of nuclei or evaluating the accuracy of contours (see an example in Fig. 18 right, in which a group of images is processed together to obtain a parametric scatter plot to be used for proof-reading patches). Fig. 19 shows an example for the visual analysis of whole slide images (WSIs). We trained an SVM model on InShaDe feature vectors derived from PanNuke data. Then, we used the SVM to classify nuclei in a large-scale, 80, 986 × 99, 328-pixel WSI of a paediatric appendix specimen. All nuclei are classified as either neoplastic (red), inflammatory (blue), or other (green). Inflamed nuclei cluster together, providing a clear indication of specific affected areas. Neoplastic nuclei are very rare and do not form structured clusters. They are therefore considered classification errors by the domain scientists. Histopathologists can use the processing pipeline for preliminary analysis targeted at the individuation of inflammatory areas. We believe (a hypothesis supported by the domain scientists in our team) that spatial aggregation of classes, i.e., density estimations of the nuclei distribution in space (also see Fig. 19), could become a valuable diagnostic tool in discriminating and individuating different tissue regions. Since a full study into the usefulness of such spatial density estimates is beyond the scope of this work, it is left as a future research direction.
Qualitative evaluation. We tested our InShaDe processing pipeline also on histology images of rodent brain samples stemming from neuroscience. Two expert neuroscientists aided this study by providing a qualitative evaluation of the framework as applied to images obtained with different staining techniques. As a general outcome, the domain scientists particularly appreciated the fact that they could try to map specific features in the shape features space to specific patterns in the histology images.
Specifically, Fig. 20, left, shows the outcomes of Nissl staining of mice brain sections. clustering was obtained with k-Means. The Nissl staining is not specific for particular cell types and is commonly used for cell counting, since it provides an excellent contrast between the cellular and extracellular space. On the other hand, it does not provide a very good contrast between the cytoplasm and the cell nucleus. In the example reported, the contrast allowed the automated algorithm to efficiently segment cell profiles, but only few nuclei were segmented correctly (mostly in light blue, some of them highlighted with blue arrows). In this case, the usage of the parameter space for highlighting the contour shapes in the image space provides visual hints for recognizing particular features, like blurred segmentations of soma mixed with dendrites, appearing as irregular and elongated shapes (see red arrows in Fig. 20 top right). So far, neuroscientists consider the framework potentially useful for proofreading the quality of the staining, and Fig. 16: Shape retrieval experiments: a simple Support Vector Machine classifier using our descriptor is able to obtain classification accuracy on par with standard geometry-based classification methods (88.4% for the MPEG-7 and 55% for Animals over the complete shape collection). We also show the full confusion matrix obtained on the testing data (left: MPEG-7, right: Animals collection). Fig. 17: Examples of combining the InShaDe pipeline with dimension reduction and clustering for visual classification in histology: Color-coding of shape clusters in the MoNuSeg data set [47] result in recognizable spatial patterns. Similar shape features according to UMAP do not necessarily relate to intuitive discrimination through visible attributes such as length or thickness or smoothness. Fig. 18: Accuracy on PanNuke data set: we trained a linear SVM model on our descriptor and we compared the accuracy with respect to different feature vectors obtained by composing energy coefficients and/or sparse coding coefficients (left). The maximum obtained accuracy over 3-classes is 0.601. PanNuke represents the largest open pan-cancer histology data set for nuclei instance segmentation and classification [41]. InShaDe can be used for proof-reading the quality and the accuracy of labelling (right).  : Visual analysis of mouse brain sections: our visual analysis pipeline is used for a neuroscience investigation. Top: a brain section fixed with paraformaldehyde is stained with Cresyl Violet, which highlights Nissl substance in the cytoplasm of neurons. Only few nuclei are segmented correctly (in light blue and highlighted with blue arrows on the right) and in various cases soma are mixed with dendrites (in red and highlighted with red arrows). Bottom: toluidine blue is used in an attempt to discriminate pyramidal neurons nuclei (in pink and highlighted by pink arrows) from blood vessels (in red and highlighted with red arrows) and artifacts. filtering some information even in case of wrong staining.
Finally, Fig. 20, right, shows a portion of somatosensory cortex from an ultrastructural work on ageing [84,85,86]. Nuclei were stained with toluidine blue on semithin sections prepared for electron microscopy in order to count cells. The extracted contours were clustered through k-Means. In this case, the shape feature space enabled scientists to distinguish immediately between blood vessels (in red, with some of them highlighted by arrows), and nuclei from different kind of neurons (pyramidal neurons mostly in pink, highlighted by arrows and with different distribution according to the layer).

Evaluation of 3D pipeline
For the evaluation of the InShaDe 3D pipeline, we used two collections of 3D reconstructions of brain cells nuclei, extracted from reconstructions of nanometric scale electron microscopy stacks, obtained after imaging a volume of brain parenchyma from layer II/III (see Fig. 21 left )and layer VI (see Fig. 21 right) somatosensory cortex of a P14 rat [54]. The nuclear shapes were manually assigned to known cell types, namely neurons, astrocytes, microglia, pericytes, unknown cells (most likely oligodendrocytes), and endotelium cells for both collections. We used InShaDe 3D feature vectors as the input for Fig. 21: Data sets: we tested the InShaDe 3D pipeline on two collections of brain cells nuclei extracted from layer II/III (left) and layer VI (right) of somatosensory cortex of a P14 rat. a kernel SVM with radial basis functions. To assess the classification performance, we considered four cases for SPH decomposition with order L max = 8, 16, 24, 32, corresponding to the number of rotation-invariant energy descriptors (see also Sec. 5.1). For each case, we performed a grid-search to configure the two hyperparameters in the SVM model: the constant γ of the Gaussian radial basis function, and the weight C for the soft margin regularization function. We chose a grid logarith- Fig. 22: Accuracy for 3D nuclei classification: we trained SVMs with InShaDe 3D features on nuclei shape collections from somatosensory cortex of a juvenile rat in layer III and layer VI. Left: cross accuracy of the SVM model for layer III collection, layer VI collection and full collection. Right: we compare the accuracy performance of InShaDe with respect to WISH [61] on the same data. Despite the simplified formulation, the accuracy is similar. mic in C (ranging from 10 −2 to 10 10 ) and γ (ranging from 10 −9 to 10 3 ). We then trained the model on 95 nuclear shapes for the layer VI shape collection, and 82 shapes for the layer II/III shape collection. We performed a 5-fold cross-validation using sklearn's StratifiedShuffleSplit function. This partitions the input data into five image sets while maintaining the relative ratio of classes in each set. Four sets were used for training and validation (using an 80/20 split) and the remaining set was used for a blind test. Fig. 22, left, summarizes the best cross-accuracy among the five folds for the SVM models trained on different shape collections and varying feature dimensions.
We also compared the performances of this new formulation with respect to our previous framework WISH [61] (also see Fig. 22 right). The InShaDe accuracy is similar to WISH, with a best score of 83% versus 84% for the layer VI data set. It is worth noting that this is despite the new formulation proposed here containing fewer coefficients (one real-valued vector here as opposed to three complex vectors [61]) and is easier to compute numerically. In Fig. 23 we report the average processing times of the Spherical Harmonics coefficients for varying orders of components. It is evident that the simplified formulation results in a dramatic reduction of computation times. These tim- Fig. 23: Processing time for SH computation: we compare the processing time for coefficients computation between InShaDe 3D and WISH [61]. The simplified formulation results in a dramatic reduction of processing times. ings were measured on a Razor Stealth laptop equipped with an Intel i7-8565U CPU (4 cores, 1.8GHz) , 16GB RAM and connected through USB-C to an e-GPU NVIDIA Titan RTX with 24GB RAM. Moreover, the spherical parametrization step is identical between the two pipelines InShaDe and WISH, and the mean curvature signal is obtained as free by-product from the usage of Willmore Flow [87]. Given the availability of shape collections extracted from different layers, we tested whether models trained on one collection could be generalized for inference on the shapes of another collection. The resulting performance was poor, in particular for neurons (below 60%). This confirms the hypothesis from domain scientists [54] that nu-clear envelopes exhibit different shape features according to the layer from which they are extracted. To confirm this point, Fig. 24 shows the full shape collection under different dimension reduction schemes. It appears that it is not possible to cluster together cells extracted from different layers (II & III vs VI). From these preliminary results, it appears evident that: • it is difficult to find models using InShaDe 3D descriptors that can generalize the classification of cell types regardless of the layers from which they are extracted; • neurons of different layers appear to form separate clusters, suggesting a shape variability depending on the area from which they are extracted [88]. This is an interesting hypothesis worthy of further investigations which we plan to carry out in the future.

Discussion
We summarize the main outcomes of this study as follows.
• Relationships between shape parameter space and image space: in various cases we notice that spatial clusters of cells exhibit closer shape features in the reduced parameter space. Further investigation is needed to understand whether and in which cases spatial patterns or clusters in the image space correspond to patterns or clusters in the parameter space, and to associate shape clusters to specific taxonomies. In this context, we would like to remind that performing clustering on parameter space obtained after dimension reduction is still considered a complex task prone to producing unreliable results [89]. Therefore, we plan to explore different automatic and manual dimension reduction techniques to support domain scientists during their analysis.
• Coupling with image descriptors: we integrated the InShaDe 2D with sparse coding for decomposing the inner part of nuclei as function of specific texture patterns with different physical and molecular characteristics. Sparse coding features outperform the proposed descriptor (we suspect that the reason is related to the fact that they can also describe shape boundary features), while composed descriptors provide slightly better discrimination capabilities for standard evaluation datasets because they recover invariance with respect to transformations. However, the composed descriptors are not reliable yet for fine-grain diagnosis on histopathology images. Spatial aggregation (i.e., class density estimation) could be used to alleviate this problem, but future research is needed. We also suspect that we have hit a performance wall for engineered and model-based descriptors. Therefore, we plan to integrate model-based descriptors into more general deep learning architectures.
• Caveats due to staining techniques: depending on the structure to be identified within a cell, or the type of tissue, a large plethora of immunohistochemical staining techniques are available. The proposed analysis framework can provide effective proof-reading tools for checking the quality of staining methods and semi-automatically individuating the structures of interest.
• Taxonomy-based visual analytics system: a real challenge in the analysis of histology images is the difficulty to individuate correct taxonomies of nuclei in order to simplify understanding and diagnosis. A visual analytics framework incorporating contour analysis, image analysis, and expert domain knowledge would help digital pathologists in labeling and proof-reading, and would provide fast ways for creating labeled data for more sophisticated artificial intelligence frameworks.To this end, our processing pipeline provides encouraging results and can be easily integrated in such systems.
• 2D Arc-length parameterization: while we have yet to observe our Arc-length parameterization algorithm to diverge, we do not have a formal proof of convergence at the time of writing. We believe it works so well since changes in u happen very gradually and the original curve remains untouched. Each reparameterization attempt therefore slides vertices around the input curve. While formal analysis is hindered by the fact that our method is discontinuous at original vertices, we believe a full treatise to be an interesting direction for future work.
• Benefits of scale-invariance: for the MPEG-7 and animals data sets, utilizing the optional scale-invariance boosts performance by up to 10% in many of our experiments. This should come as no surprise, since, e.g., the outline of a butterfly stays the outline of a butterfly under magnification and minification. Consequently, deep learning-based pipelines have made rescaling a main step of their data augmentation stage, in an attempt the achieve de-facto rather than by-design scale-invariance. What is surprising, however, is that scale-invariance added only marginal and in many cases statistically insignificant improvements for cell nuclei classification. We believe that this is due to the "apparent size" problem, in which cell nuclei always appear smaller than the original size due to slicing. It seems that having plenty of slices under different angles at the classifiers disposal is more important than to remove scale-variance. A full analysis of this problem is beyond the scope of this paper and offers an interesting future research direction.
• Limitations of 3D pipeline: the encouraging results obtained with our 3D formulation are counterbalanced by two important limiting bottlenecks. Firstly, the process for producing nuclear surfaces from electron microscopy image stacks is still time consuming and requiring highly specialized human efforts. Even though important progresses in automatic segmentation of EM stacks has recently been made [19], custom models for automatic extraction of nuclei are not available to our knowledge. We plan to focus future efforts towards this direction. Secondly, the spherical parameterization task is complex and can be unstable. One of its limitations is that it cannot be applied to arbitrary closed shapes but only genus-0 and (if the flow is appropriately regularized) genus-1 surface (that is, surfaces either homeomorphic to spheres or torii with at most 1 hole). To overcome these limitations, we plan to investigate more general invariant formulations based on manifold harmonics [90].

Conclusion
We have presented a general shape processing framework rooted in a novel differential-geometry-based descriptor of closed contours and surfaces. Our descriptor provides an embedding into a fixed-dimensional feature space that can be utilized for various applications, which range from serving as input feature for deep and shallow learning techniques to supporting dimension reduction schemes for providing a visual reference for clustering collection of shapes. While our methods are of general use, our work is motivated by the study of cellular nuclear envelopes extracted from histopathological images and serial section electron microscopy stacks. In this context, we have shown the capabilities of the proposed framework for visual analysis and unsupervised classification. Our results are very encouraging and we identify several major areas of future work in the previous section. In particular, we plan to develop, on top of our pipeline, a taxonomy-based visual analytics system to simplify study and diagnosis.