extfieldsofstudy:
Computer Science
Engineering

text:
Constructing High-Order Signed Distance Maps from Computed Tomography Data with Application to Bone Morphometry

An algorithm is presented for constructing high-order signed distance fields for two phase materials imaged with computed tomography. The signed distance field is high-order in that it is free of the quantization artifact associated with the distance transform of sampled signals. The narrowband is solved using a closest point algorithm extended for implicit embeddings that are not a signed distance field. The high-order fast sweeping algorithm is used to extend the narrowband to the remainder of the domain. The order of accuracy of the narrowband and extension methods are verified on ideal implicit surfaces. The method is applied to ten excised cubes of bovine trabecular bone. Localization of the surface, estimation of phase densities, and local morphometry is validated with these subjects. Since the embedding is high-order, gradients and thus curvatures can be accurately estimated locally in the image data.

I. INTRODUCTION
This paper develops a method of constructing signed distance fields from computed tomography datasets. Primarily, the method is concerned with two-phase materials such as trabecular bone, synthetic biomaterials, and soil samples. Being biphasic defines the interface between the two materials as an orientable surface, advantageous for modelling changes over time and measuring metrics of the surface such as volume, area, curvatures, and thickness.
Classically, a signed distance field would be constructed in two steps. First, an appropriate segmentation method would generate a binary image defining the object in space. Second, the exact signed distance transform would produce the distance map from the binary image [1]. However, the exact signed distance transform computes a quantized version of the true signed distance field [2]. This produces errors for geometric flow and curvature computation algorithms because the algorithm is only initialized first-order accurate, resulting in large errors in measurements and curvature computation. Furthermore, constructing high-order signed distance maps from binary sampled signals is a non-unique problem [3], permitting multiple embeddings that satisfy the same binary image.
The task presented here is to construct the signed distance field with a high-order of accuracy from data measured using computed tomography. The field will be constructed directly from grayscale data where gradients are well-defined. Once constructed, the densities of the respective phases can be estimated. Finally, it is demonstrated that morphometry can be performed accurately on these embeddings, a task that is too noisy when computed from the exact signed distance transform.

II. METHOD A. Problem Definition
The problem is to construct a signed distance field, φ, from a density-calibrated computed tomography image, ρ. Both will be represented on a sampled, rectilinear n-dimensional grid Ω ⊂ Z n . Let the domain Ω be split between the two phases, Ω 1 and Ω 2 such that Ω = Ω 1 ∪ Ω 2 . Associate to each phase a density, ρ 1 and ρ 2 . Γ is the interface between the two materials. We represent the interface implicitly as the zero level set of an embedding function u.
If the embedding is known, the density field can be reconstructed knowing the density of the two phases.
where θ is the Heaviside function and the convention inside-isnegative is used. This construction originates in the simulation of incompressible two-phase flow [4]. In many cases, we want the embedding function to satisfy an additional constraint, namely the Eikonal equation: where ∇ is the gradient operator, | · | is the vector norm, and F is a spatially varying speed function which yields the distance function when F = 1. Here, we used the symbol φ to distinguish the signed distance field from any general embedding u.
A signed distance field is highly desirable for many reasons. First, it can be used directly for many morphometric measures such as measuring areas, volumes, curvatures, and thicknesses [5], [6], [7]. Second, the signed distance field can represent arbitrary topology surfaces and evolve them without splitting and merging rules as would be needed for parametric meshes [8]. Third, during surface evolution problems, the Eikonal equation provides a way of rapidly redistancing the embedding during evolution [4], [9]. Constructing φ from ρ is the central challenge of this paper.

B. Solver
We present a method of constructing φ from ρ in computed tomography datasets given the assumptions of Section II-A. An overview of the method is given in Figure 1.
1) Constructing an Implicit Surface: An implicit representation of the surface is constructed by shifting the image by a threshold. This causes the zero level set to be at the level of the threshold. Standard preprocessing can also be performed on the image such as smoothing to remove noise. A simple expression is used in this work: where G σ is a Gaussian kernel with standard deviation σ, * is convolution, and T is a threshold. The symbol ψ is used to denote that this is not a signed distance transform. Importantly, this is an implicit representation of the surface Γ that is subvoxel. This construction is inspired by a sub-pixel distance map method [10]. Of note, the threshold function can be made spatially varying, T = T (x), important to high resolution bone scans where cortical and spongy bone compartments are often segmented with different thresholds [11]. Similarly, alternative preprocessing filters can be used that match the specific application.
2) Solving the Narrowband: A modification of the closest point method [12] is used to construct the signed distance function in the narrowband. The concept is to take each point on the rectilinear grid, x, that is in the narrowband and find the closest point, y, on the implicit surface: where d(·, ·) is the Euclidean distance. The distance between these two points gives the embedding function: where the sign is determined from the sign of ψ(x). This is demonstrated graphically in Figure 2. The closest point method uses gradients of ψ to iteratively move a point in the narrowband onto the zero level set.
However, the closest point method originally used a signed distance field, which is not available in this case. The algorithm is modified for the embedding in Equation 4.
Consider a point in the narrowband. This point is moved in the negative direction of the gradient until it is sufficiently close to the zero level set. However, the embedding gradient points "uphill" (away from the medial axis) in both the positive and negative sections of the curve. As such, the direction is modified by multiplying the sign of the embedding so both interior and exterior points move towards the zero level set. The update equation is then: where λ is a step size, sgn(·) is the sign function, and n(y) is the normal at the point y. As is standard in image processing and level set methods, the normal is a function of the gradient: The process is repeated till the embedding at the point y is sufficiently close to zero. Equation 7 can be rewritten as a continuous partial differential equation: Visual representation of the problem. To construct the signed distance field φ, the distance to the implicit curve Γ is found. This is done first in the narrowband (gray) and then marched outwards. The surface Γ is shown as the gray level set. Vector length approaches zero as they get closer to the curve due to the regularized sign function.
where τ is a virtual time parameter. The closest point method is then a forward Euler approximation for the initial value integration problem where the point is integrated through a stationary vector field: This vector field is plotted in Figure 3. To solve Equation 9, an interpolator is needed such that the field ψ can be assessed between sample points. We construct the interpolator using basis splines: where β n is the n th order normalized basis spline, c(q) are the coefficients of the interpolant, andψ(x) is the interpolating function for ψ. Conveniently, once the coefficients have been constructed, derivatives can be evaluated quickly: This allows the computation of normals and magnitude gradients without having to compute finite differences in the image. B-splines have been well studied [13], [14], [15] and are particularly appropriate for this task. First, they have a finite support, making them rapid to evaluate. Second, the smoothness of the interpolating function is the order of the normalized basis spline, allowing a simple selection of the order of accuracy. The drawback is that the construction of the basis splines can take some time and requires additional memory to store the coefficients. However, the time cost of construction is made up by the speed of evaluation.
Combining the interpolator and vector field integrator gives the closest point algorithm presented in Algorithm 1. The integration time step is selected as λ = h, making the time step on the order of a voxel. Third order basis splines are used in the interpolant for accurate computation of normals. The integration stops when the point is sufficiently close to the implicit surface, given by h 3 . Since a signed distance function Algorithm 1 Closest Point Algorithm, CP . The vector field v is defined in Equation 10 and computed using the regularized sign equation and basis spline interpolator. Input: ψ, x, h Output: y Initialization 1: y = x 2:ψ = BuildBSpline(ψ) Process 3: while |ψ(y)| > h 3 |∇ψ(y)| do 4: y ← y − h sgn(ψ(y)) n(y) 5: end while 6: return y is not used, the stopping criterion is modified as described in Appendix A.
A numerical approximation is made for the sign function: where δ is a small constant preventing the denominator from being zero, selected to be the sample spacing, h. The regularization constant δ is weighted by the gradient magnitude for the same reasons presented in Appendix A. Such a regularization has been used previously [9]. Regularization of the sign function has the additional benefit of causing the rate of Equation 7 to approach zero as the point y n approaches the implicit curve, preventing oscillations around the zero crossing. For images with anisotropic spacing, h is selected as the minimum of the voxel edge lengths. An issue with Algorithm 1 is that it does not enforce collinearity with the surface [16] where the vector from x to y must be collinear with the surface normal at y. Even with high-order Runge-Kutta methods, integration through the vector field does not satisfy the collinearity property. In two dimensions, the candidate point y can be corrected by moving along the tangent [16]. However, in high dimensions, there is not a single tangent but a tangent plane. We extend the collinearity method to arbitrary dimensions here.
Let x be the candidate point and define p to be the current solution from CP . Define the vector q as pointing from p to x. Define T p M to be the tangent plane at the point p on the manifold. The candidate point is moved by a vector being the projection of q onto T p M .
This is shown graphically in Figure 4. The new candidate point is now the original plus this projection.
The projected vector z will be zero when it is collinear with the normal of the tangent plane, giving a natural stopping condition. At areas of large curvature, Equation 16 can cause oscillations in the candidate point. To avoid oscillations, the tangent plane vector is weighted: where β ∈ [0, 1] is the step size, taken to be 0.5 in this work. The magnitude of z can be understood as the size of the open ball at p n and enforcing the radius to be less than a power of the spacing enforces an order of accuracy. This produces slower convergence in flat areas where β = 1 would be exact, but greatly improves the convergence in high curvature areas.
The surfaces of interest have very few flat locations, making this trade-off advantageous. The order of accuracy is selected the same of CP , h 3 . This is termed the Collinear Closest Point Algorithm [16], CP ⊥ , and is described in Algorithm 2.
In rare cases, CP or CP ⊥ will not converge. This can occur when the initial point is at the medial axis of the image, the gradient of the edge is too flat causing noise in the normal estimation, or the noise in the density image is too large. A maximum iteration is set (100 in this work for both algorithms) and if the maximum iterations are reached, the point is solved with the method of Section II-B3. This causes a degradation of accuracy locally, but keeps the implementation stable.
z ← proj TpM ( yx) 5: y ← CP (ψ, y + 0.5z, h) 6: end while 7: return y 3) Extending off the Narrowband: After solving the narrowband, the signed distance field must be extended to the rest of the domain. This is done using the fast sweeping method of Zhao [17], [18]. The fast marching method [19], [20], [21] and fast iterative method [22] are also suitable algorithms for extending the narrowband, but the authors prefer the simplicity of the fast sweeping method.
As signed distance fields are non-smooth surfaces (technically, Lipschitz continuous), they cannot be better in general than first-order accurate at discontinuities [18], such as the medial axis where characteristics flow into each other. However, away from the medial axis, arbitrary order of accuracy can be achieved.
The sweeping method proceeds in three steps. First, the absolute value of the narrowband voxels are taken to remove the sign condition. This solves for the unsigned distance field from the surface. Next, an upwind quadratic solver is used to solve the quadratic equations at each voxel. Importantly, as the method is upwind, the solution is propagated from the narrowband to the rest of the domain. The solver is fast because alternating sweeping directions cascades the solution from the surface to the domain, Finally, the sign condition is reintroduced by multiplying the solved distance field by the sign of ψ. There are techniques to include the sign condition in the solver [3], but they are considerably slower due to the branching requirements and are considerably more difficult to implement.

C. Component Estimation
Now that a non-parametric representation of the interface is available, we would like to measure the parameters ρ 1 and ρ 2 . We follow here the method of Chan and Vese [6]: where ρ 1 is the density of the phase inside the embedding and ρ 2 is the density of the phase outside the embedding. In the context of trabecular bone, ρ 1 is the tissue mineral density and ρ 2 the marrow mineral density. A nice feature of this method of estimating the phase densities is that the reconstructed density image (Equation 2) will have an average density near the original image, allowing for the computation of average density from the embedded image and the phase densities. Of note, components are estimated from the smoothed density image, G σ * ρ, not the implicit surface, ψ, or the original density image, ρ. This keeps the phase density estimation and surface localization consistent with any filtering. Numerical approximations are needed for the integration and Heaviside function. Numerical integration is performed with Simpson's rule, particularly important for surfaces with few samples across their thickness. A regularized version of the Heaviside function is used: where 2 is the finite support of the numerical Heaviside θ . A finite support is advantageous for structures with few samples across their thickness as it avoids extending into the medial axis.

D. Narrowband Size
The voxels in the narrowband must be determined. Since the initial embedding is not a signed distance transform, the standard absolute value threshold cannot be used: where |·| is the absolute value function and γ is the narrowband size. Instead, the set of voxels in the narrowband is determined by thresholding the embedding below zero and performing morphological operators: where ⊕ is morphological dilation, is morphological erosion, and S is a suitable structuring element. This is considerably slower than the thresholding approach but selects a technically correct narrowband. The cross structuring element is used in this work and the radius is set to half the stencil size minus one.
This makes the narrowband large enough that a finite difference stencil placed outside the narrowband does not extend across the narrowband. The image J then has the value one for all voxels in the narrowband and zero for all elements outside the narrowband.

A. Correctness on Ideal Surfaces
A series of experiments are performed on ideal surfaces to validate the correctness of the algorithms. The signed distance field of a sphere and torus are known analytically: where r is the radius of the sphere, a is the radius of the torus, and c is the distance to the center of the torus. This allows us to construct a ground truth φ, create a density image ρ, and assess the ability of the algorithm to recover φ from ρ. Additionally, two intersecting spheres are constructed, producing a continuous but non-smooth implicit surface. The surface parameters are r = 5 for the sphere, a = 2, c = 3 for the torus, and r = 5 for the double spheres with separation distance 3 √ 3. Images were instantiated on a 20 × 20 × 20 grid (nominal physical units).
We are particularly concerned with the order of accuracy of the method. Decreasing the grid size will always increase the accuracy of the method, but with appropriate methods the accuracy improves to some power of the sampled spacing. This is important for measuring morphometric quantities such as curvature that require a certain amount of smoothness in the embedding.
In general, the order of accuracy, m, is estimated by: where φ is the ground truth, φ h is the solution at spacing h, φ h/2 is the solution at spacing h/2, and · p is an p norm.
Since the distance transform is only Lipschitz continuous at the medial axis, different norms are used to approximate measures of the global (everywhere) and local (away from medial axis) order of accuracy. Because the signed distance map is Lipschitz continuous, the method can only be firstorder accurate when measured with the ∞ norm but can achieve arbitrary order of accuracy when measured with the 1 norm [18]. 1) Order of Accuracy in the Narrowband: The order of accuracy for solving the narrowband is measured. The analytic signed distance transform of the sphere, torus, and doublesphere are instantiated on a grid with a prescribed spacing. Using Equation 2, phase densities ρ 1 = 100, ρ 2 = 0 (in arbitrary density units), and Heaviside smoothness constant = 2, the density image is instantiated. The narrowband is solved and the 1 and ∞ norm measured on the narrowband voxels. The grid size is divided by two and repeated, allowing an estimation of the order of accuracy of the narrowband initialization. No smoothing is performed such to preserve the zero crossing location and the threshold is set to the average of the phase densities: One can see this threshold is correct from Equations 2 and 4 noting that θ(0) = 0.5.
2) Order of Accuracy in Sweeping: The order of accuracy of the signed fast sweeping method is measured. The analytic signed distance transform of the sphere, torus, and doublesphere are instantiated at the narrowband with a prescribed spacing. The remainder of the domain is solved from the narrowband and the 1 and ∞ norms are measured excluding narrowband voxels. The grid size is divided by two and the process repeated, allowing an estimation of the order of accuracy of the extension method.
3) Visualizing Result: The ideal signed distance transform of a double sphere is sampled at h = 0.25 and mapped into density as in Section III-A1. From the density image, the signed distance transform is solved using the narrowband and subsequent fast sweeping method. Different level sets of the ideal and solved distance transform are visualized using the Marching Cubes algorithms [23] to demonstrate the visual smoothness and correctness of the method.

B. Trabecular Bone
10 previously imaged bovine trabecular bone cubes were used to verify the correctness of the algorithm [24]. Bone cubes were sawed to 10 mm in edge length and imaged at a nominal resolution of 20 µm. The bones were embedded with Equation 4 using a standard deviation σ = 20 µm and threshold T = 400 mg HA/cc. The signed distance transform was computed with the proposed method.
To increase the execution speed, the first-order fast sweeping method was used for extending the narrowband. A high accuracy is only needed at the narrowband for visualization, phase density estimation, and morphometry. The models here had an average of 143 million voxels, requiring roughly 20 minutes to compute the narrowband and roughly 40 minutes per iteration in the fast sweeping method on a consumer grade MacBook Pro (32 GB random access memory, 2.3 GHz processor with 8 cores). The fast sweeping method can be accelerated by only solving to a certain distance, but a highorder cannot be parallelized due to the lack of a monotonicity condition [18]. It should be noted that the narrowband method is a so-called "embarrassingly" parallelizable algorithm, where every narrowband voxel could be solved on a separate thread.
1) Embedding and Surface Localization: A visual comparison is made between the density field, the implicit surface, and the solved signed distance transform to verify the location of the level set coincides between images. The exact signed distance transform is also computed from the implicit surface to demonstrate the quantization artifact [25]. Meshes are extracted using the Marching Cubes algorithms [23] at the appropriate level set for the given field.
2) Phase Densities: Next, density estimation is validated. The phase components are estimated using Equations 18 and 19. Once estimated, the density image is reconstructed using Equation 2 and compared visually. The bone mineral density (vBMD [mg HA/cc]) between the original density image and the constructed density image are compared using regression and Bland-Altman analysis [26]. The regularization parameter was set to = 2h.
χ =K 2π (35) where V is the object volume, A is the object surface area, H is the mean curvature, K is the Gaussian curvature,H,K are the total curvatures, H , K are the surface average curvature, χ is the Euler-Poincaré characteristic, θ is the Heaviside function, and δ is the Dirac delta function. The curvatures are estimated using fourth-order central differences [7] and integration is performed using Simpson's rule from Section II-C. Furthermore, the numerical estimation of θ in Equation 20 is used. Taking the derivative of Equation 20 gives the numerical approximation to the Dirac delta function: From these definitions, traditional bone morphometrics [27] can be derived: can be computed simply as the ratio of BV to TV. The link between mean curvature and SMI was first described by Jinnai et al. [31] and extended to TBPf by Stauber and Müller [32]. The proposed method is compared to the corresponding standard methods for each measure. The original density image is blurred with the same Gaussian filter and thresheld to binarize the volume. It is standard in bone morphometry to perform a largest connected component filter to force one bone component and one marrow component in the image. In fact, connected component filtering is required to meet the assumptions on Betti numbers in connectivity density estimation [30]. As a result, the standard methods are applied to the binary image after largest connected component filtering. However, connected component filter cannot be performed in the proposed method in a straight forward manner. As a result, Euler-Poincaré characteristic and not connectivity density is used for validation. Additionally, since SMI and TBPf are both related to mean curvature, only SMI is used for comparison. BV/TV, BS, SMI, and χ are compared between the proposed and standard methods using regression and Bland-Altman plots. Standard analysis was performed with Image Processing Language (IPL v5.42, SCANCO Medical AG, Brüttisellen, Switzerland). Finally, curvature is visualized across the surface. The Marching Cubes algorithm [23] is run

A. Correctness on Ideal Surfaces
1) Order of Accuracy in the Narrowband: Measured order of accuracy is reported in Table I. For the sphere and torus, the narrowband has an order of accuracy of 4 in both the 1 and ∞ norms. This is an order higher than expected from analysis, a result seen in similar studies [21]. Since the narrowband is smooth, the ∞ and 1 norms have the same order of accuracy.
For sphere-sphere intersection, the 1 order of accuracy is 2 while the ∞ order of accuracy is one. The degenerate 1 order of accuracy comes from the implicit surface itself not being smooth, so the order of accuracy of the method cannot be achieved. Similarly, a degenerate 1 order of accuracy is seen, reflecting the presence of a non-smooth solution in the narrowband. The ∞ order of accuracy is worse than one at higher resolutions. This is caused by points in the narrowband that are non-smooth being mapped to a wrong but collinear point on the surface. At the sphere-sphere intersection, nonsmooth locations are mapped onto one of the two spheres and not the intersection, while still satisfying the collinear requirement.
2) Order of Accuracy in Sweeping: Measured order of accuracy in sweeping is reported in Table II. For the sphere and torus, an order of accuracy of 3 is seen in the 1 norm with an order of accuracy of 1 in the ∞ norm. The reduced ∞ corresponds to a shock at their center. Across all grid sizes, the number of iterations was less than 40 The non-smooth sphere-sphere intersection demonstrated an order of accuracy near 2 for the 1 norm and an order of accuracy of 1 in the ∞ norm. Errors were seen in areas of rarefaction, corresponding to the lens shape of the intersecting spheres, degrading the 1 norm order of accuracy. Outside the area of rarefaction the result was as expected. This result is well known in the sweeping method [18]. The number of iterations required for convergence was larger for the nonsmooth surface.
3) Visualizing Result: The ideal and solved double sphere surface is rendered at four level sets in Figure 6. The zero level set surfaces correspond very well with no visual indication of the quantization artifact. In the positive domain, the level set   in rarefaction. Finally, the dumbbell separates into two teardrop shapes, some errors can be seen in the roundness of the solution.
B. Trabecular Bone 1) Embedding and Surface Localization: The surface contour of ρ, ψ, and φ for a single slice are overlaid in Figure 7. Gaussian blurring reduces noise and smooths the surface while slightly modifying the level set location. However, no further modification to the level set is seen in the solved signed distance transform. The exact signed distance transform produces large errors as expected. The meshed level sets are displayed in Figure 8. Noise has been removed by blurring and the solved distance transform agrees with the filtered image. Importantly, the exact signed distance transform applied to the binarized image exhibits artifacts due to quantization of the metric.
2) Phase Densities: Measured phase densities are reported in Table III. The average bone phase density is 646.6 ± 34.0 mg HA/cc (3734.8 ± 192.0 HU) and the average marrow phase density is 3.8 ± 26.5 mg HA/cc (100.9 ± 149.8 HU), reported mean ± standard deviation. Average densities are separated on either side of the 400 mg HA/cc threshold. Negative effective densities are seen in the marrow space because the calibration model does not account for densities below water. However, the marrow space Hounsfield units are in appropriate ranges between fat (−100 HU) and general soft tissue (100 HU). The Hounsfield units of the bone phase are very large, but are reasonable given the high resolution of the imager. When weighting the bone and marrow components by their volume fraction (Section IV-B3), the image Hounsfield units range from 535.5 HU to 1575.9 HU, consistent with clinical computed tomography where partial voluming would average the components together. Consistency in the phase densities is expected since these were the same species at the same location and there were no known diseases affecting mineralization or marrow composition.
The filtered image and constructed image are shown in Figure 9. Noise has been removed due to replacing the phases with their average values. Similarly, because a regularized Heaviside function is used and the high-order signed distance transform has gradients, the constructed image has gradients at the edge between the phases. However, the constructed image has a thicker appearance than the original image with gradients not matching local hypo-intensities in the bone phase. This suggests the regularization parameter could be included in optimization. It should be noted that the zero-crossings of ψ, φ, and the constructed image correspond (Figure 7), so the object is not actually thicker. The object just appears thicker because of discrepancies in threshold, local intensities, and the regularization parameter.
A regression and Bland-Altman plot between the volumetric bone minearal density of the original (ψ) and constructed image are given in Figure 10. A small proportional bias is seen in the Bland-Altman analysis (R 2 = 66.1%) caused by selecting the same density threshold independent of the object densities. Both the bias and proportional bias is tiny (1.41%−0.69%) in comparison to the volumetric bone mineral density. The component estimation method of Equations 18 and 19 has the property that the constructed volume will always agree with the original image on average density. As a result, vBMD errors originate in the ground truth data not matching the assumptions of a biphasic material, either because of noise or that the phases do not exhibit a single density.
3) Morphometry: Regression and Bland-Altman plots are displayed in Figure 11. Agreement in bone volume fraction (BV/TV) is very high suggesting very few voxels were removed by component labelling. A slight proportional bias is seen, consistent with the proportional bias in vBMD ( Figure 9). Agreement in bone surface area (BS) is high (R 2 = 100.0%) with a small bias of 2.03 [−29.35, 33.42] (95% LoA). The proposed method underestimated surface area, consistent with a sub-voxel localization of the surface that is smaller than the binarized volume. Excellent agreement (R 2 = 100.0%) and small bias (0.02 [−0.06, 0.10]) is seen in structure model index (SMI). The standard method estimates SMI using a finite sized dilation of the surface on the order of the sample spacing. This appears to result in an overestimation of SMI relative to the infinitesimal definition of the proposed method. Finally, Euler-Poincaré characteristic (χ) shows reduced agreement (R 2 = 96.3%) with a modest bias of −1857.78 [−3241.42, −474.1265]. The standard method requires component labelling to estimated χ, removing many small components which would add to the total Euler-Poincaré characteristic. As a result, the proposed method overestimated χ relative to the standard method.
Mean and Gaussian curvature are displayed across the surface of a single sample in Figure 12. Rendered mean curvature is smooth across the surface, showing negative and positive curvature across rods and plates. The mean curvature histogram shows a spike around zero where the image border clips the bone flat. Rendered Gaussian curvature is less smooth, owing to the distribution being fat tailed across the surface. The fat tailed distribution is confirmed in the histogram. The smoothness of both the mean and Gaussian curvature agrees with the smoothness of ψ in Figure 8b.

V. DISCUSSION
A method is presented for constructing high-order signed distance transforms of two phase materials imaged with computed tomography. The narrowband is solved by integrating voxel positions through a stationary velocity field. The narrowband is extended to the remainder of the domain using a fast sweeping method. The phase densities are estimated from the original density image and computed signed distance field.  (13c) At shocks, a point in the domain has two or more closest points, consistent with Blum's grassfire burn [33]. All other points have a one-to-one and onto mapping to the implicit curve for a given level set. Rarefaction and shocks should be understood as regions in the related problem of simulating hyperbolic conservative laws.
The utility of the method is demonstrated by measuring bone morphometrics directly from the embedding. The novelty of the method is to solve the narrowband from an embedding that is not the signed distance transform. This is done by integrating narrowband voxels through a stationary velocity field until they are both sufficiently close to the zero level set and sufficiently collinear with the normal. The modification of the convergence criterion for an image without a gradient magnitude of one was sufficient to obtain the prescribed order of accuracy. The computation of a distance map from a grayscale image was inspiried by, but varies considerably from, the sub-pixel distance map method of Kimmel et al. [10] where a redistancing algorithm is used to construct the high-order signed distance map (called sub-pixel in that work). Originally, the narrowband method was inspired by closest point methods [12], [16] but extended for arbitrary dimensions and embeddings that are not already a signed distance field. The same nomenclature for the CP and CP ⊥ algorithms are used. The narrowband algorithm is most closely related to the redistancing algorithm of Chopp [21] with a different interpolator and modified convergence criterion.
The foundational assumption of the method was the construction of an embedding function that had the proper zero crossing. A smoothing filter and threshold standard to bone analysis was modified to construct the implicit surface, relying on density calibration. In general, any implicit representation will suffice as long as the zero crossing corresponds to the surface, the embedding is monotone in the narrowband, and the embedding has well-defined gradients. However, many image analysis tasks cannot make this assumption. Integration with alternative definitions of edges [34], [35], [36], [37] would increase the utility of the technique. Such methods can even be directly integrated into active contour segmentation algorithms which have a penalty term forcing the embedding to be a signed distance field [38].
The limitation of the method is that the surface must be smooth, coming from smoothness assumptions in the interpolator and velocity field. This is a good assumption for the exact task of embedding biphasic materials imaged with computed tomography where filtering has smoothed the surface. However, the accuracy degraded at continuous but non-smooth corners of the surface. As such, we do not recommend our method for redistancing [4], [9], [21].
Establishing the embedding for non-smooth surfaces is very challenging. An attempt to alleviate the issue could be to force the divergence of the vector field to zero. This is the case in cortical thickness mapping where Laplace's equation provides a one-to-one mapping between two surfaces [39]. As a result, the integration paths from one surface to another cannot cross. However, the current problem does not permit a one-to-one mapping between level sets of the embedding because there can be multiple points on the surface that are closest (in the shock regions of advection) or multiple points on the same level set that map to a single point on the surface ( Figure 13).
An improvement to the method would be to estimate the regularization parameter, , from the image data. One approach would be to estimate in conjunction with ρ 1 and ρ 2 based on a non-linear least squares. It is expected that the magnitude of the regularization parameter corresponds to the size of the point spread function of the imager. Such a result would improve the accuracy of the gradients in the constructed image. Additionally, the method could be fully automated by including a threshold estimation step. Otsu's method [40] is particularly appropriate for biphasic materials, and would give a distribution instead of point estimate for phase densities.

VI. CONCLUSION
A method is presented for computing high-order signed distance transforms of biphasic materials imaged with computed tomography. Phase densities and general morphometry can be estimated from these embeddings. The algorithm is highorder, free of the quantization artifact associated with distance transforms of sampled binary images.

APPENDIX APPROXIMATING DISTANCE FROM THE SURFACE
It is desirable to be able to select the convergence criterion based on a distance from the surface. This allows measuring the convergence of CP in physical units, permitting one to prescribe an order of accuracy.
The geometric picture is given in Figure 14. A Taylor series approximation to the implicit contour is made at the zero crossing: where δx is a directional vector. The direction is selected as the normal vector of a prescribed length: where γ is the length. Second, we know ψ(x) = 0 by definition of the zero crossing. As such, we get the following approximation for the distance from the contour: We modify the convergence criterion to be: This exact expression avoids division by zero in flat areas of ψ. Furthermore, it prevents early stopping away from the interface since |ψ| is greater than zero but |∇ψ| approaches zero near the medial axis. Additionally, the expression |∇ψ| is already available from the computation of the normal in Algorithm 1, avoiding any additional computation. Finally, the error in the approximation decreases as the point approaches the surface at a rate O(γ 2 ). If ψ is a signed distance function, Equation 46 becomes the same convergence criterion as before since the gradient magnitude is unity. In general, |∇ψ| can be seen as a normalizing factor bringing the gradient magnitude close to one such that the distance travelled in the embedding is the same as the distance travelled in the domain.


Explaining Time Series Predictions with Dynamic Masks

How can we explain the predictions of a machine learning model? When the data is structured as a multivariate time series, this question induces additional difficulties such as the necessity for the explanation to embody the time dependency and the large number of inputs. To address these challenges, we propose dynamic masks (Dynamask). This method produces instance-wise importance scores for each feature at each time step by fitting a perturbation mask to the input sequence. In order to incorporate the time dependency of the data, Dynamask studies the effects of dynamic perturbation operators. In order to tackle the large number of inputs, we propose a scheme to make the feature selection parsimonious (to select no more feature than necessary) and legible (a notion that we detail by making a parallel with information theory). With synthetic and real-world data, we demonstrate that the dynamic underpinning of Dynamask, together with its parsimony, offer a neat improvement in the identification of feature importance over time. The modularity of Dynamask makes it ideal as a plug-in to increase the transparency of a wide range of machine learning models in areas such as medicine and finance, where time series are abundant.

Introduction and context
What do we need to trust a machine learning model? If accuracy is necessary, it might not always be sufficient. With the application of machine learning to critical areas such as medicine, finance and the criminal justice system, the black-box nature of modern machine learning models has appeared as a major hindrance to their large scale deployment (Caruana et al., 2015;Lipton, 2016;Ching et al., 2018). With the necessity to address this problem, the field Proceedings of the 38 th International Conference on Machine Learning, PMLR 139, 2021. Copyright 2021 by the author(s). of explainable artificial intelligence (XAI) thrived (Barredo Arrieta et al., 2020;Das & Rad, 2020;Tjoa & Guan, 2020).
Saliency Methods Among the many possibilities to increase the transparency of a machine learning model, we focus here on saliency methods. The purpose of these methods is to highlight the features in an input that are relevant for a model to issue a prediction. We can distinguish them according to the way they interact with a model to produce importance scores.
Gradient-based: These methods use the gradient of the model's prediction with respect to the features to produce importance scores. The premise is the following: if a feature is salient, we expect it to have a big impact on the model's prediction when varied locally. This translates into a big gradient of the model's output with respect to this feature. Popular gradient methods include Integrated Gradient (Sundararajan et al., 2017), DeepLIFT (Shrikumar et al., 2017) and GradSHAP (Lundberg et al., 2018).
Perturbation-based: These methods use the effect of a perturbation in the input on the model's prediction to produce importance scores. The premise is similar to the one used in gradient based method. The key difference lies in the way in which the features are varied. If gradient based methods use local variation of the features (according to the gradient), perturbation-based method use the data itself to produce a variation. A first example is Feature Occlusion (Suresh et al., 2017) that replaces a group of features with a baseline. Another example is Feature Permutation that performs individual permutation of features within a batch.
Attention-based: For some models, the architecture allows to perform simple explainability tasks, such as determining feature saliency. A popular example, building on the success attention mechanisms (Vaswani et al., 2017), is the usage of attention layers to produce importance scores (Choi et al., 2016;Song et al., 2017;Xu et al., 2018;Kwon et al., 2018).
Other: There are some methods that don't clearly fall in one of the above categories. A first example is SHAP (Lundberg & Lee, 2017), which attributes importance scores based on Shapley values. Another popular example is LIME in which the importance scores correspond to weights in a local linear model (Ribeiro et al., 2016a). Finally, some methods such as INVASE (Yoon et al., 2019a) or ASAC (Yoon et al., 2019b) train a selector network to highlight important features.
Time Series Saliency Saliency methods were originally introduced in the context of image classification (Simonyan et al., 2013). Since then, most methods have focused on images and tabular data. Very little attention has been given to time series (Barredo Arrieta et al., 2020). A possible explanation for this is the increasing interest in model agnostic methods (Ribeiro et al., 2016b).
Model agnostic methods are designed to be used with a very wide range of models and data structures. In particular, nothing prevents us from using these methods for Recurrent Neural Networks (RNN) trained to handle multivariate time series. For instance, one could compute the Shapley values induced by this RNN for each input x t,i describing a feature i at time t. In this configuration, all of these inputs are considered as individual features and the time ordering is forgotten. This approach creates a conceptual problem illustrated in Figure 1. Context matters. This graph shows a fictional security price over time. A simplistic model recommends to sell just after a local maximum, to buy just after a local minimum and to wait otherwise. When recommending to sell or to buy, 3 consecutive time steps are required as a context.
In this example, the time dependency is crucial for the model. At each time step, the model provides a decision based on the two previous time steps. This oversimplified example illustrates that the time dependency induces a context, which might be crucial to understand the prediction of some models. In fact, recurrent models are inherently endowed with this notion of context as they rely on memory cells. This might explain the poor performances of model agnostic methods to identify salient features for RNNs (Ismail et al., 2020). These methods are static as the time dependency, and hence the context, is forgotten when treating all the time steps as separate features.
If we are ready to relax the model agnosticism, we can take a look at attention based methods. In methods such as RE-TAIN (Choi et al., 2016), attention weights are interpreted as importance scores. In particular, a high attention weight for a given input indicates that this input is salient for the model to issue a prediction. In this case, the architecture of the model allows us to provide importance scores without altering the dynamic nature of the data. However, it has been shown recently that the attention weights of a model can significantly be changed without inducing any effect on the predictions (Jain & Wallace, 2019). This renders the parallel between attention weights and feature importance rather unclear. This discrepancy also appears in our experiments.
Another challenge induced by time series is the large number of inputs. The number of features is multiplied by the number of time steps that the model uses to issue a prediction. A saliency map can therefore become quickly overwhelming in this setting. To address these challenges, it is necessary to incorporate some treatment of parsimony and legibility in a time series saliency method. By parsimony, we mean that a saliency method should not select more inputs than necessary to explain a given prediction. For instance, a method that identifies 90% of the inputs as salient is not making the interpretation task easier. By legibility, we mean that the analysis of the importance scores by the user should be as simple as possible. Clearly, for long time sequences, analysing feature maps such as (c) and (d) on Figure 3 can quickly become daunting.
To address all of these challenges, we introduce Dynamask. It is a perturbation-based saliency method building on the concept of masks to produce post-hoc explanations for any time series model. Masks have been introduced in the context of image classification (Fong & Vedaldi, 2017;Fong et al., 2019). In this framework, a mask is fitted to each image. The mask highlights the regions of the image that are salient for the black-box classifier to issue its prediction. These masks are obtained by perturbing the pixels of the original image according to the surrounding pixels and study the impact of such perturbations on the black-box prediction. It has been suggested to extend the usage of masks beyond image classification (Phillips et al., 2018;Ho et al., 2019). However, to our knowledge, no available implementation and quantitative comparison with benchmarks exists in the context of multivariate time series. Moreover, few works in the literature explicitly mention explainability in a multivariate time series setting (Siddiqui et al., 2019;Tonekaboni et al., 2020). Both of these considerations motivate our proposition.
Contributions By building on the challenges that have been described, our work is, to our knowledge, the first saliency method to rigorously address the following questions in a time series setting.
(1) How to incorporate the context? In our framework, this is naturally achieved by studying the effect of dynamic perturbations. Concretely, a perturbation is built for each feature at each time by using the value of this feature at adjacent times. This allows us to build meaningful perturbations that carry contexts such as the one illustrated in Figure 1.
(2) How to be parsimonious? A great advantage of using masks is that the notion of parsimony naturally translates into the extremal property. Conceptually, an extremal mask selects the minimal number of inputs allowing to reconstruct the black-box prediction with a given precision.
(3) How to be legible? To make our masks as simple as possible, we encourage them to be almost binary. Concretely, this means that we enforce a polarization between low and high saliency scores. Moreover, to make the notion of legibility quantitative, we propose a parallel with information theory 1 . This allows us to introduce two metrics : the mask information and the mask entropy. As illustrated in Figure 3, the entropy can be used to assess the legibility of a given mask. Moreover, these metrics can also be computed for other saliency methods, hence allowing insightful comparisons.
The paper is structured as follows. In Section 2, we outline the mathematical formalism related to our method. Then, we evaluate our method by comparing it with several benchmarks in Section 3 and conclude in Section 4.

Mathematical formulation
In this section, we formalize the problem of feature importance over time as well as the proposed solution. For the following, it is helpful to have the big picture in mind. Hence, we present the blueprint of Dynamask in Figure 2. The rest of this section makes this construction rigorous.

Preliminaries
Let X ⊂ R d X be an input (or feature) space and Y ⊂ R d Y be an output (or label) space, where d X and d Y are respectively the dimension of the input and the output space. For the sake of concision, we denote by [n 1 : n 2 ] the set of natural numbers between the natural numbers n 1 and n 2 with n 1 < n 2 . We assume that the data is given in terms of time series (x t ) t∈[1:T ] , where the inputs x t ∈ X are indexed by a time parameter t ∈ [1 : T ] with T ∈ N * . We consider the problem of predicting a sequence (y t ) t∈[ty:T ] of outputs y t ∈ Y indexed by the same time parameter, but starting at a time t y ≥ 1, hence allowing to cover sequence to vector predictions (t y = T ) as well as sequence to sequence predictions (t y < T ). For the following, it is convenient to introduce a matrix notation for these time series. In this way, X = (x t,i ) (t,i)∈[1:T ]×[1:d X ] denotes the matrix in R T ×d X whose rows correspond to time steps and whose columns correspond to features. Similarly 2 , we denote Y = (y t,i ) (t,i)∈[ty:T ]×[1:d Y ] for the matrix in R (T +1−ty)×d Y . 1 Many connections exist between explainability and information theory, see for instance (Chen et al., 2018). 2 In the following, we do not make a distinction between R T ×d X and X T as these are isomorphic vector spaces. The same goes for R (T +1−ty )×d Y and Y T +1−ty .
Our task is to explain the prediction Y = f (X) of a black box f that has been pre-trained for the aforementioned prediction task. In other words, our method aims to explain individual predictions of a given black box. More specifically, our purpose is to identify the parts of the input X that are the most relevant for f to produce the prediction Y. To do this, we use the concept of masks introduced in (Fong & Vedaldi, 2017;Fong et al., 2019), but in a different context. Here, we adapt the notion of mask to a dynamic setting. Definition 2.1 (Mask). A mask associated to an input sequence X ∈ R T ×d X and a black box f : X T → Y T +1−ty is a matrix M = (m t,i ) ∈ [0, 1] T ×d X of the same dimension as the input sequence. The element m t,i of this matrix represents the importance of feature i at time t for f to produce the prediction Y = f (X). Remark 2.1. For a given coefficient in the mask matrix, a value close to 1 indicates that the feature is salient, while a value close to 0 indicates the opposite.
Now how can we obtain such a mask given a black box and an input sequence? To answer this question, we should let the mask act on the inputs and measure the effect it has on a black box prediction. As the name suggests, we would like this mask to hide irrelevant inputs contained in X. To make this rigorous, it is useful to spend some time to think about perturbation operators.

Perturbation operators
The method that we are proposing in this paper is perturbation-based. Concretely, this means that a mask is used to build a perturbation operator. Examples of such operators are given in (Fong et al., 2019) in the context of image classification. Here, we propose a general definition and we explain how to take advantage of the dynamic nature of the data to build meaningful perturbations. By recalling that the mask coefficients indicate the saliency, we expect the perturbation to vanish for features x t,i whose mask coefficient m t,i is close to one. This motivates the following definition. Definition 2.2 (Perturbation operator). A perturbation operator associated to a mask M ∈ [0, 1] T ×d X is a linear operator acting on the input sequence space Π M : R T ×d X → R T ×d X . It needs to fulfil the two following assumptions for any given (t, i) ∈ [1 : T ] × [1 : d X ]: 1. The perturbation for x t,i is dictated by m t,i : where π is differentiable for m ∈ (0, 1) and continuous at m = 0, 1.
2. The action of the perturbation operator is trivial when the mask coefficient is set to one :

Loss
Black Box Figure 2. Diagram for Dynamask. An input matrix X, extracted from a multivariate time series, is fed to a black-box to produce a prediction Y. The objective is to give a saliency score for each component of X. In Dynamask, these saliency scores are stored in a mask M of the same shape as the input X. To detect the salient information in the input X, the mask produces a perturbed version of X via a perturbation operator Π. This perturbed X is fed to the black-box to produce a perturbed prediction YM. The perturbed prediction is compared to the original prediction and the error is backpropagated to adapt the saliency scores contained in the mask.
Moreover, we say that the perturbation is dynamic if for any given feature at any given time, the perturbation is constructed with the values that this feature takes in neighbouring times. More formally, the perturbation is dynamic if there exist a couple (W 1 , Remark 2.2. The first assumption ensures that the perturbations are applied independently for all inputs and that our method is suitable for gradient-based optimization. The second assumption 3 translates the fact that the perturbation should have less effect on salient inputs. Remark 2.3. The parameters W 1 and W 2 control how the perturbation depends on neighbouring times. A non-zero value for W 1 indicates that the perturbation depends on past time steps. A non-zero value for W 2 indicates that the perturbation depends on future time steps. For the perturbation to be dynamic, at least one of these parameters has to be different from zero.
This definition gives the freedom to design perturbation operators adapted to particular contexts. The method that we develop can be used for any such perturbation operator. In our case, we would like to build our masks by taking advantage the dynamic nature of the data into account. This is crucial if we want to capture local time variations of the features, such as a quick increase of the blood-pressure or 3 Together with the continuity of π with respect to the mask coefficients.
an extremal security price from Figure 1. This is achieved by using dynamic perturbation operators. To illustrate, we provide three examples: where π g is a temporal Gaussian blur 4 with Similarly, π m can be interpreted as a fade to moving average perturbation with where W ∈ N is the size of the moving window. Finally, π p is similar to π m with one difference: the former only uses past values of the features to compute the perturbation 5 : Note that the Hadamard product M X used in (Ho et al., 2019) is a particular case of π m with µ t,i = 0 = so that this perturbation is static. In the following, to stress that a mask is obtained by inspecting the effect of a dynamic perturbation operator, we shall refer to it as a dynamic mask (or Dynamask in short).
In terms of complexity, the computation of π g requires O(d X · T 2 ) operations 6 while π m requires O(d X · T · W ) operations 7 . When T is big 8 , it might therefore be more interesting to use π m with W T or a windowed version of π g . With this analysis of perturbation operators, everything is ready to explain how dynamic masks are obtained.

Mask optimization
To design an objective function for our mask, it is helpful to keep in mind what an ideal mask does. From the previous subsection, it is clear that we should compare the black-box predictions for both the unperturbed and the perturbed input. Ideally, the mask will identify a subset of salient features contained in the input that explains the black-box prediction. Since this subset of features is salient, the mask should indicate to the perturbation operator to preserve it. More concretely, a first part of our objective function should keep the shift in the black-box prediction to be small. We call it the error part of our objective function. In practice, the expression for the error part depends on the task done by the black-box. In a regression context, we minimize the squared error between the unperturbed and the perturbed prediction: Similarly, in a classification task, we minimize the crossentropy between the predictions: Now we have to make sure that the mask actually selects salient features and discards the others. By remembering that m t,i = 0 indicates that the feature x t,i is irrelevant for the black-box prediction, this selection translates into imposing sparsity in the mask matrix M. A first approach for this, used in (Fong & Vedaldi, 2017;Ho et al., 2019), is to add a l 1 regularisation term on the coefficients of M. However, it was noted in (Fong et al., 2019) that this produces mask that vary with the regularization coefficient λ in a way that renders comparisons between different λ difficult. To solve this issue, they introduce a new regularization term 6 One must compute a perturbation for each dX · T component of the input and the sums in each perturbation have T terms 7 One must compute a perturbation for each dX · T component of the input and each moving average involves 2W + 1 terms 8 In our experiments, however, we never deal with T > 100 so that both approaches are reasonable. to impose sparsity: where . denotes the vector 2-norm, vecsort is a function that vectorizes M and then sort the elements of the resulting vector in ascending order. The vector r a contains (1 − a) · d X · T zeros followed by a · d X · T ones, where a ∈ [0, 1]. In short, this regularization term encourages the mask to highlight a fraction a of the inputs. For this reason, a can also be referred to as the area of the mask. A first advantage of this approach is that the hyperparameter a can be modulated by the user to highlight a desired fraction of the input. For instance, one can start with a small value of a and slide it to higher values in order to see features gradually appearing in order of importance. Another advantage of this regularization term is that it encourages the mask to be binary, which makes the mask more legible as we will detail in the next section.
Finally, we might want to avoid quick variations in the saliency over time. This could either be a prior belief or a preference of the user with respect to the saliency map. If this is relevant, we can enforce the salient regions to be connected in time with the following loss that penalizes jumps of the saliency over time: For a fixed fraction a, the mask optimization problem can therefore be written as Note that, in this optimization problem, the fraction a is fixed. In some contexts, one might want to find the smallest fraction of input features that allows us to reproduce the black-box prediction with a given precision. Finding this minimal fraction corresponds to the following optimization problem: where ε is the threshold that sets the acceptable precision. The resulting mask M * a * is then called extremal mask. The idea of an extremal mask is extremely interesting by itself: it explains the black-box prediction in terms of a minimal number of salient features. This is precisely the parsimony that we were referring to in the introduction.
Finally, it is worth mentioning that we have here presented a scheme where the mask preserves the features that minimize the error. There exists a variant of this scheme where the mask preserves the features that maximize the error. This other scheme, together with the detailed optimization algorithm used in our implementation, can be found in Appendix B.

Masks and information theory
Once the mask is obtained, it highlights the features that contain crucial information for the black-box to issue a prediction. This motivates a parallel between masks and information theory. To make this rigorous, we notice that the mask admits a natural interpretation in terms of information content. As aforementioned, a value close to 1 for m t,i indicates that the input feature x t,i carries information that is important for the black box f to predict an outcome. It is therefore natural to interpret a mask coefficient as a probability that the associated feature is salient for the black box to issue its prediction. It is tempting to use this analogy to build the counterpart of Shannon information content in order to measure the quantity of information contained in subsequences extracted from the time series. However, this analogy requires a closer analysis. We recall that the Shannon information content of an outcome decreases when the probability of this outcome increases (Shannon, 1948;MacKay, 2003;Cover & Thomas, 2005). In our framework, we would like to adapt this notion so that the information content increases with the mask coefficients (if m t,i gets closer to one, this indicates that x t,i carries more useful information for the black-box to issue its prediction). To solve this discrepancy, we have to use 1−m t,i as the pseudoprobability appearing in our adapted notion of Shannon information content.
Definition 2.3 (Mask information). The mask information associated to a mask M and a subsequence ( Remark 2.4. Conceptually, the information content of a subsequence measures the quantity of useful information it contains for the black-box to issue a prediction. It allows us to associate a saliency score to a group of inputs according to the mask. Remark 2.5. Note that, in theory, the mask information diverges when a mask coefficient is set to one. In practice, this can be avoided by imposing M ∈ (0, 1) T ×d X .
As in traditional information theory, the information content is not entirely informative on its own. For instance, consider two subsequences indexed by, respectively, A, B with |A| = |B| = 10. We assume that the submask extracted from M with A contains 3 coefficients m = 0.9 and 7 coefficients m = 0 so that the information content of A is I M (A) ≈ 6.9. Now we consider that all the coefficient extracted from M with B are equal to 0.5 so that I M (B) ≈ 6.9 and hence I M (A) ≈ I M (B). In this example, A clearly identifies 3 important features while B gives a mixed score for the 10 features. Intuitively, it is pretty clear that the information provided by A is sharper. Unfortunately, the mask information by itself does not allow to distinguish these two subsequences. Hopefully, a natural distinction is given by the counterpart of Shannon entropy.
Definition 2.4 (Mask entropy). The mask entropy associated to a mask M and a subsequence ( Remark 2.6. We stress that the mask entropy is not the Shannon entropy for the subsequence (x t,i ) (t,i)∈A . Indeed, the pseudo-probabilities that we consider are not the probabilities p(x t,i ) for each feature x t,i to occur. In particular, there is no reason to expect that the probability of each input decouple from the others so that it would be wrong to sum the individual contribution of each input separately 9 like we are doing here.
Clearly, it is desirable for our mask to provide explanations with low entropy. This stems from the fact that the entropy is maximized when mask coefficients m t,i are close to 0.5. In this case, given our probabilistic interpretation, the mask coefficient is ambiguous as it does not really indicate whether the feature is salient. This is consistent with our previous Since masks coefficients take various values in practice, masks with higher entropy appear less legible, as illustrated in Figure 3. Therefore, we use the entropy as a measure of the mask's sharpness and legibility in our experiments. In particular, the entropy is minimized for perfectly binary masks M ∈ {0, 1} T ×d X , which are easy to read and contain no ambiguity. Consequently, the regularization term L a (M) that we used in Section 2.3 has the effect of reducing the entropy. Since our adapted notions of information and entropy rely on individual contributions from each feature, they come with natural properties.  Masks with high entropy appear less legible, especially for long sequences.
Proof. See Appendix A.
Together, these properties guarantee that I M and S M define measures for the discrete σ-algebra

Experiments
In this section 10 , we evaluate the quality of our dynamic masks. There are two big difficulties to keep in mind when evaluating the quality of a saliency method. The first one is that the true importance of features is usually unknown with real-world data. The second one is that the the performance of a saliency method in identifying relevant features depends on the black-box and its performances. In order to illustrate these challenges, we propose three different experiments in ascending order of difficulty. In the first experiment, a white box with known feature importance is used so that both difficulties are avoided. In the second experiment, a blackbox trained on a dataset with known feature importance is used and hence only the second difficulty is encountered. In the third experiment, a black-box trained on a real-world clinical dataset is used and hence both difficulties are encountered. For each experiment, more details are given in Appendix C.

Feature importance for a white-box regressor
Experiment In this experiment, we work with a trivial white-box regressor whose predictions only rely on a known 10 Our implementation can be found at https://github. com/JonathanCrabbe/Dynamask.
where A T and A X respectively give the salient times and features: Note that, in this experiment, d Y = 1 so that we can omit the second index for f (X). In our experiment, we consider two scenarios that are known to be challenging for state of the art saliency methods (Ismail et al., 2020). In the first one, the white-box depends on a small portion of salient features |A X | d X . In the second one, the white-box depends on a small portion of salient times |A T | T . We fit a mask to this black-box by using the squared error loss together with the regularization. In our experiment: T = d X = 50. The salient time and features are selected randomly, the input features are generated with an ARMA process. We repeat the experiment 10 times.
Metrics Since the true salient features are known unambiguously in this setup, we have a large set of metrics at hand to evaluate the performance of a saliency method. To measure the proportion of identified features that are indeed salient, we use the area under the precision curve (AUP, higher is better). To measure the portion of salient features that have indeed been identified, we use the area under the recall curve (AUR, higher is better). To measure how much information the saliency method predicts for the salient region, we use the mask information content (I M (A), higher is better). To measure the sharpness of the explanation in this region, we use the mask entropy (S M (A), lower is better).
Benchmarks The fact that we are dealing with a whitebox regressor already disqualifies some methods such as FIT (Tonekaboni et al., 2020) or DeepLIFT. We compare our method with Feature Occlusion (FO), Feature Permutation (FP), Integrated Gradient (IG) and Shapley Value Sampling (SVS) (Castro et al., 2009). For fair comparison, we use these baselines to evaluate the importance of each feature at each given time. As explained in Section 2.4, a mask can be associated to each saliency method.

Feature importance for a black-box classifier
Experiment We reproduce the state experiment from (Tonekaboni et al., 2020) but in a more challenging setting. In this case, the data is generated according to a 2-state hidden Markov model (HMM) whose state at time t ∈ [1 : T ] (T = 200) is denoted s t ∈ {0, 1}. At each time, the input feature vector has three components (d X = 3) and is generated according to the current state via x t ∼ N µ st , Σ st . To each of these input vectors is associated a binary label y t ∈ {0, 1}. This binary label is conditioned by one of the three component of the feature vector, based on the state: The label is then emitted via a Bernoulli distribution y t ∼ Ber(p t ). The experiment proposed in (Tonekaboni et al., 2020) focuses on identifying the salient feature at each time where a state transition occurs. However, it is clear from the above discussion that there is exactly one salient feature at any given time. More precisely, the set of salient indexes is A = { (t, 2 + s t ) | t ∈ [1 : T ] }. We generate 1000 such time series, 800 of them are used to train a RNN black-box classifier f with one hidden layer made of 200 hidden GRU cells 11 . We then fit an extremal mask to the black-box by minimizing the cross entropy error for each test time series. We repeat the experiment 5 times.
Metrics We use the same metrics as before.

Feature importance on clinical data
Experiment We reproduce the MIMIC mortality experiment from (Tonekaboni et al., 2020). We fit a RNN blackbox with 200 hidden GRU cells to predict the mortality of a patient based on 48 hours (T = 48) of patient features (d X = 31). This is a binary classification problem for which the RNN estimates the probability of each class. For each patient, we fit a mask with area a to identify the most important observations x t,i . Since the ground-truth feature saliency is unknown, we must find an alternative way to assess the quality of this selection. If these observations are indeed salient, we expect them to have a big impact on the black-box's prediction if they are replaced. For this reason, we replace the most important observations by the time average value for each associated feature x t,i →x t,i = 1 T T t=1 x t,i . We then compare the prediction for the original input f (X) with the prediction for the input where the most important observations have been replaced f X . A big shift in the prediction indicates that these observations are salient. We repeat this experiment 3 times for various values of a.
Dataset We use the MIMIC-III dataset , that contains the health record of 40, 000 ICU deidentified patients at the Beth Israel Deaconess Medical Center. The selected data and its preprocessing is the same as the one done by (Tonekaboni et al., 2020).
Metrics To estimate the importance of the shift induced by replacing the most important observations, we use the crossentropy between f (X) and f X (CE, higher is better). To evaluate the number of patient whose prediction has been flipped, we compute the accuracy of f X with respect to the initial predictions f (X) (ACC, lower is better).
Benchmarks We use the same benchmarks as before.

Discussion
The results are shown in Figure 4 & 5. The observations selected by Dynamask have the most significant impact when replaced. The high accuracy suggests that the perturbation we use (replacing the most important observations by their time average) rarely produces a counterfactual input.

Conclusion
In this paper, we introduced Dynamask, a saliency method specifically designed for multivariate time series. These masks are endowed with an insightful information theoretic interpretation and offer a neat improvement in terms of performance. Dynamask has immediate applications in medicine and finance, where black-box predictions require more transparency.
For future works, it would be interesting to design more sophisticated consistency tests for saliency methods in a dynamic setting, like the ones that exist in image classification (Adebayo et al., 2018). This could be used to study the advantages or disadvantages of our method in more details.
Another interesting avenue would be to investigate what the dynamic setting can offer to provide richer explanations with some treatment of causality (Moraffah et al., 2020).

A. More details on the mathematical formulation A.1. Proofs
In this subsection section, we prove the propositions form the main paper. Proof. Let us proof all the properties one by one.
Positivity By definition, all coefficients from the mask are normalized: Positivity follows trivially from the properties of the logarithm function The same goes for the entropy Additivity We first note that the proposition follows trivially if the sets are distinct A ∩ B = ∅: (2) We shall now prove the additivity property in general by using these two ingredients. First we note that the set A ∪ B can be written as the disjoint union A [B \ (A ∩ B)]. It follows that where we have used the additivity property for disjoint sets in the first equality and the fact that (A ∩ B) ⊂ B in the second equality. The same reasoning holds for the entropy.
Monotonicity To prove the monotonicity property, it is useful to note that if A ⊂ B, we can use (2) to write where we have used the information positivity to produce the inequality. The same reasoning holds for the entropy.

A.2. Normalized information and entropy
In this subsection, we introduce the normalized counterparts of our information theoretic metrics. It is important to keep in mind that all the available information for issuing a black-box prediction is in [1 : T ] × [1 : d X ]. Therefore, the monotonicity property allows to introduce normalized counterparts of the information and the entropy.
Definition A.1 (Normalized metrics). The normalized mask information associated to a mask M and a subsequence .
The same goes for the related normalized mask entropy .
Remark A.1. By the monotonicity and the positivity properties, it is clear that This gives a natural interpretation of these quantities as being, respectively, the fraction of the total information and entropy contained in the subsequence A according to the mask M.
The normalized version of the metrics allow to measure what percentage of the total mask information/entropy is contained in a given subsequence.

A.3. Definition of a mask for other saliency methods
In this section, we explain how to associate a mask to any saliency method. Suppose that a given method produces a score matrix R ∈ R T ×d X that assigns an importance score r t,i for each element x t,i of the input matrix X. Then, if we normalize the coefficients of the score matrix, we obtain an associated mask: where (1) T ×d X denotes a T × d X matrix with all elements set to 1. This mask can subsequently be used to compute the mask information content and entropy. In our experiments, we use this correspondence to compare our method with popular saliency methods.
B. More details on the implementation B.1. Algorithm

Algorithm 1 Dynamask
Input: input sequence X ∈ R T ×d X , black-box f , perturbation operator Π, mask area a ∈ [0, 1], learning rate η ∈ R + , momentum α ∈ R + , initial size regulator The mask optimization algorithm is presented in Algorithm 1. In the algorithm, we used the notation (0.5) T ×d X for a T × d X matrix with all elements set 13 to 0.5. Similarly, (0) T ·d X ·(1−a) denotes a vector with T · d X · (1 − a) components set to 0 and (1) T ·d X ·a denotes a vector with T · d X · a components set to 1. The symbol ⊕ denotes the direct sum between two vector spaces, which is equivalent to the concatenation in Algorithm 1. The error part of the loss L e depends on the task (regression or classification), as explained in Section 3 of the paper. The momentum and the learning rate are typically set to 1, the number of epoch is typically 1000. We also use the clamp function, which is defined component by component as Finally, we note that the mask size regularization coefficient λ a grows exponentially during the optimization to reach a maximum value of δ · λ 0 at the end of the optimization. In practice, it is initialized to a small value (typically λ 0 = 0.1) and dilated by several order of magnitude during the optimization (typically δ = 1000). In this way, the optimization procedure works in two times. At the beginning, the loss is dominated by the error L e (M) so that the mask increases the mask coefficients of salient features. As the regulation coefficient λ a increases, the regulation term becomes more and more important so that the mask coefficients are attracted to 0 and 1. At the end of the optimization, the mask is almost binary.

B.2. Deletion variant
We notice that Algorithm 1 produces a mask that highlights the features that allow to reproduce the black-box prediction by keeping the error part of the loss L e (M) to be small. However, it is possible to highlight important features in another way. For instance, we could try to find the features that maximizes the prediction shift when perturbed. In this alternative formulation, the mask is obtained by solving the following optimization problem: Note that, in this case, the sign of the error part is flipped in order to maximizes the shift in the prediction. Moreover, the error is now evaluated for 1 − M rather than M. This is because important features are maximally perturbed in this case. In this way, a salient feature x t,i can keep a mask coefficient m t,i close to 1 while being maximally perturbed. The regulator stays the same in this deletion variant, as the mask area still corresponds to the number of mask coefficients set to 1. We use the deletion variant to obtain the masks in the experiment with clinical data in the main paper.

C. More details on the experiments C.1. Metrics
We give the precise definition of each metric that appears in the experiments. Let us start with the metrics that are defined when the true importance is known. Definition C.1 (AUP,AUR). Let Q = (q t,i ) (t,i)∈[1:T ]×[1:d X ] be a matrix in {0, 1} T ×d X whose elements indicate the true saliency of the inputs contained in X ∈ R T ×d X . By definition, q t,i = 1 if the feature x t,i is salient and 0 otherwise. Let M = (m t,i ) (t,i)∈[1:T ]×[1:d X ] be a mask in [0, 1] T ×d X obtained with a saliency method. Let τ ∈ (0, 1) be the detection threshold for m t,i to indicate that the feature x t,i is salient. This allows to convert the mask into an estimatorQ(τ ) = (q t,i (τ )) (t,i)∈[1:T ]×[1: Consider the sets of truly salient indexes and the set of indexes selected by the saliency method We define the precision and recall curves that map each threshold to a precision and recall score: The AUP and AUR scores are the area under these curves Remark C.1. Roughly speaking, we consider the identification of salient features as a binary classification task. Each saliency method can thus be seen as a binary classifier for which we compute the AUP and the AUR. Remark C.2. Integrating over several detection thresholds allows to evaluate a saliency method with several levels of tolerance on what is considered as a salient feature.
In our experiment with MIMIC-III, since the ground true feature importance is unknown, we use the following metrics defined for a binary classification problem 14 . Definition C.2 (CE, ACC). Consider a classifier f that maps the input X to a probability f (X) ∈ [0, 1]. LetX be a perturbed input produced by a saliency method 15 . We define the function that converts a probability into a class class(p) = 0 if p < 0.5 1 else.
14 In our experiment, each input corresponds to a patient. Class 0 indicates that the patient survives and class 1 indicates that the patient dies. 15 In our experiment, we replace the most important features by the time average of the corresponding feature.
To measure the shift in the classifier's prediction caused by the perturbation of the input for several test examples {X k | k ∈ [1 : K]}, we use the binary cross-entropy (or logloss) To measure the number of prediction flipped by the perturbation, we use the accuracy We reproduce our experiment several times to get an average and a standard deviation for all of these metrics.

C.2. Computing infrastructure
All our experiments have been performed on a machine with Intel(R) Core(TM) i5-8600K CPU @ 3.60GHz [6 cores] and Nvidia GeForce RTX 2080 Ti GPU.

C.3. Details on the rare experiment
Data generation Since this experiment relies on a whitebox, we only have to generate the input sequences. As we explain in the main paper, each feature sequence is generated with an ARMA process: with ϕ 1 = 0.25, ϕ 2 = 0.1, ϕ 3 = 0.05 and t ∼ N (0, 1).
We generate one such sequence with t ∈ [1 : 50] for each feature i ∈ [1 : 50] by using the Python statsmodels library.
In the rare feature experiment, 5 features are selected as salient. Their indices are contained in A X and drawn uniformly without replacement from Mask fitting For each time series, we fit a mask by using the temporal Gaussian blur π g as a perturbation operator with σ max = 1 and by using the squared error loss. A mask is fitted for each value of a ∈ (n + 1) · 10 −3 | n ∈ [0 : 49] .
In our experiments, we don't consider a > 0.05. This is because we found experimentally that the error L e (M * a )  Runtime For rare time, finding the best mask takes on average 15.7s. For rare feature, finding the best mask takes on average 20.7s.
Illustrations To illustrate the results of our experiments, we show the saliency masks produced by various methods for the rare feature experiment in Figure 8 & 9 and for the rare time experiment in Figure 10 & 11. For all of these examples, we notice that Dynamask identifies a bigger portion of the truly salient inputs, which illustrates the bigger AUR reported in the main paper.

C.4. Details on the state experiment
Data generation The data generation is governed by a Hidden Markov Model (HMM). The initial distribution vector for this HMM is given by π = (0.5, 0.5) and its transition matrix is C = 0.1 0.9 0.1 0.9 .
At each time, the input feature vector has three components (d X = 3) and is generated according to the current state via x t ∼ N µ st , Σ st with mean vectors depending on To each of these input vectors is associated a binary label y t ∈ {0, 1}. This binary label is conditioned by one of the three component of the feature vector, based on the state: The length of each time series is fixed to 200 (T = 200).
We generate 1000 such time series, 800 are used for model training and 200 for testing.

Model training
We train a RNN with one layer made of 200 GRU cells trained using the Adam optimizer for 80 epochs (lr = 0.001, β 1 = 0.9, β 2 = 0.999 and no weight decay).
Mask fitting For each test time series, we fit a mask by using the temporal Gaussian blur π g as a perturbation operator with σ max = 1. A mask is optimized for each value of a ∈ 0.15 + 2n · 10 −2 | n ∈ [0 : 10] . We keep the extremal mask for a threshold set to ε = 0.9 · L e M = (1) T ×d X .
Runtime Finding the extremal mask for a given input takes 49.8s on average.
Illustrations To illustrate the results of our experiments, we show the saliency masks produced by various methods on Figure 12 & 13. By inspecting these figures, we notice that only Dynamask and Integrated Gradients seem to produce saliency maps where the imprint of the true saliency can be distinguished. One advantage of Dynamask is the contrast put between these salient inputs and the rest. In the case of Integrated Gradients, we see that many irrelevant inputs are assigned an important saliency score, although smaller than the truly salient inputs. This is because gradients are computed individually, without the goal of achieving parsimonious feature selection.
In addition, we have reported the mask entropy for each of the methods. As claimed in the main paper, Dynamask produces mask that have significantly lower entropy. Among the methods that produce masks with high entropy, we notice two trends. Some methods, such as RETAIN, produce masks where a significant portion of the inputs are assigned a mask entropy close the 0.5. As discussed in the main paper, these significance of the saliency scores is limited in this situation, since no clear contrast can be drawn between the saliency of different inputs. On the other hand, some methods like FIT produce masks with many different masks coefficients, which renders the saliency map somewhat fuzzy. In both cases, the high entropy detects these obstructions for legibility.

C.5. Details on the mimic experiment
Data preprocessing The data preprocessing used here is precisely the same as the one described in (Tonekaboni et al., 2020), we summarize it here for completeness. We use the adult ICU admission data from the MIMIC-III dataset   Mask fitting For each test patient, we simply fit a mask with a = 0.1 by maximizing the cross-entropy loss in the deletion variant formulation of Dynamask. We use the fadeto-moving average perturbation π m with W = 48. The hyperparameters for this optimization procedure are η = 1, α = 1, λ 0 = 0.1, λ c = 0, δ = 1000, N = 1000.
Runtime Fitting a mask for a given patient takes 3.58 s on average.
Illustrations To illustrate the results of our experiments, we show the 10% most important features for patients that are predicted to die on Figure 14 & 15 and for patients that are predicted to survive on Figure 16 & 17. In each case, we indicate the cross-entropy between the unperturbed and the perturbed prediction, as defined in Definition C.2. We note that Dynamask identifies the features that create the biggest shift. Qualitatively, we note that Dynamask seems to focus much more on the input that appear at latter time. This is consistent with the observations in (Ismail et al., 2019): these inputs are the most important for the black-box, since it is trained to predict the mortality after the 48 hours and RNNs have short term memory.

C.6. Influence of the perturbation operator
To study the effect of the perturbation operator choice, we have performed the following experiment: in the setup of the state experiment, we optimize 100 masks on distinct examples by using a Gaussian blur perturbation (π g , σ max = 1) and a fade-to-moving average perturbation (π m , W = 3). We do the same with a fade-to-past average perturbation that only uses past values of the features: (π p , W = 6). For each pair of perturbation operators, we compute the average accuracy between the associated masks (i.e. the fraction of inputs where both masks agree). For each method, we report the AUROC for the identification of true salient inputs. The results are reported in Table 4. We observe that π g , π m , π p generally agree and offer similar performances. The black-box predicts that this patient will survive. In each case, the cross entropy (CE) between the unperturbed and the perturbed prediction is reported.